[Rd] acepack

[Frank E Harrell Jr]

On 01 Jun 2003 12:53:05 +0200
Peter Dalgaard BSA <p.dalgaard at biostat.ku.dk> wrote:

> Notice that one of your (FEH's) fixes is only easy if the statement
> "whereas I believe that the ace algorithm does allow for categorical
> response variables" is taken at face value. Checking that it is
> actually correct could take a bit longer...

The following verifies that my local version of ace (see below) with the fixes works when y is categorical.

# Verify that ace works for categorical response variable, giving
# a y-transformation that is a linear translation of Fisher's optimum scores
# (y-specific mean of x) when there is one predictor that is forced to be linear
# For now using a local version of ace
library(acepack)
set.seed(1)
y <- rep(1:3,100)
x <- -3*(y==1) -7*(y==2) + 30*(y==3) + runif(300) - .5
xbar <- tapply(x, y, mean)
xbar
        1         2         3 
-3.010843 -7.021050 30.002227 

z <- ace(x, y, cat=0, lin=1)
table(y, z$ty)
   
y   -0.824639636875907 -0.582656111909047 1.40729574878496
  1   0                100                  0             
  2 100                  0                  0             
  3   0                  0                100             

plot(xbar[y], z$ty)  # perfect linearity
cor(xbar[y], z$ty)

[1] 1

Here is ace with fixes.  Changes are marked with 'FEH'.

ace <- function (x, y, wt = rep(1, nrow(x)), cat = NULL, mon = NULL, 
    lin = NULL, circ = NULL, delrsq = 0.01) 
{
    x <- as.matrix(x)
    if (delrsq <= 0) {
        cat("delrsq must be positive")
        return()
    }
    iy <- ncol(x) + 1
    l <- matrix(1, ncol = iy)
    if (!is.null(circ)) {
        for (i in 1:length(circ)) {
            if (circ[i] < 0 || circ[i] > ncol(x)) {  # FEH nrow -> ncol
                cat("bad circ= specification")
                return()
            }
            if (circ[i] == 0) {
                cat("response spec can only be lin or ordered (default)")
                return()
            }
            else {
                nncol <- circ[i]
                if (l[nncol] != 2 & l[nncol] != 1) {
                  cat("conflicting transformation specifications")
                  return()
                }
                l[nncol] <- 2
            }
        }
    }
    if (length(mon)) {
        for (i in 1:length(mon)) {
            if (mon[i] < 0 || mon[i] > ncol(x)) {  # FEH nrow -> ncol
                cat("bad mon= specification")
                return()
            }
            if (mon[i] == 0) {
                cat("response spec can only be lin or ordered (default)")
                return()
            }
            else {
                nncol <- mon[i]
                if (l[nncol] != 3 && l[nncol] != 1) {
                  cat("conflicting transformation specifications")
                  return()
                }
                l[nncol] <- 3
            }
        }
    }
    if (length(lin)) {
        for (i in 1:length(lin)) {
            if (lin[i] < 0 || lin[i] > ncol(x)) {  # FEH nrow -> ncol
                cat("bad lin= specification")
                return()
            }
            if (lin[i] == 0) {
                nncol <- iy
            }
            else {
                nncol <- lin[i]
            }
            if (l[nncol] != 4 && l[nncol] != 1) {
                cat("conflicting transformation specifications")
                return()
            }
            l[nncol] <- 4
        }
    }
    if (length(cat)) {
        for (i in 1:length(cat)) {
            if (cat[i] < 0 || cat[i] > ncol(x)) {  # FEH nrow -> ncol
                cat("bad cat= specification")
                return()
            }
            if (cat[i] == 0) {
#  Next 2 lines commented out FEH
#                cat("response spec can only be lin or ordered (default)")
#                return()
            }
            else {
                nncol <- cat[i]
                if (l[nncol] != 4 && l[nncol] != 1) {
                  cat("conflicting transformation specifications")
                  return()
                }
                l[nncol] <- 4
            }
        }
    }
    tx <- x
    ty <- y
    m <- matrix(0, nrow = nrow(x), ncol = iy)
    z <- matrix(0, nrow = nrow(x), ncol = 12)
    z <- as.matrix(z)
    ns <- 1
    mode(x) <- "double"
    mode(y) <- "double"
    mode(tx) <- "double"
    mode(ty) <- "double"
    mode(wt) <- "double"
    mode(delrsq) <- "double"
    mode(z) <- "double"
    junk <- .Fortran("mace", p = as.integer(ncol(x)), n = as.integer(nrow(x)), 
        x = t(x), y = y, w = as.double(wt), l = as.integer(l), 
        delrsq = delrsq, ns = as.integer(ns), tx = tx, ty = ty, 
        rsq = double(1), ierr = integer(1), m = as.integer(m), 
        z = z, PACKAGE = "acepack")
    return(junk)
}


---
Frank E Harrell Jr              Prof. of Biostatistics & Statistics
Div. of Biostatistics & Epidem. Dept. of Health Evaluation Sciences
U. Virginia School of Medicine  http://hesweb1.med.virginia.edu/biostat