[Rd] acepack
Frank E Harrell Jr
fharrell at virginia.edu
Sun Jun 1 21:23:16 MEST 2003
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
More information about the R-devel
mailing list