Therneau, Terry M., Ph.D.
2015-Jun-11 13:44 UTC
[R] Different behavior of model.matrix between R 3.2 and R3.1.1
Frank,
I'm not sure what is going on. The following test function works for me
in both 3.1.1
and 3.2, i.e, the second model matrix has fewer columns. As I indicated to you
earlier,
the coxph code removes the strata() columns after creating X because I found it
easier to
correctly create the assign attribute.
Can you create a worked example?
require(survival)
testfun <- function(formula, data) {
tform <- terms(formula, specials="strata")
mf <- model.frame(tform, data)
terms2 <- terms(mf)
strat <- untangle.specials(terms2, "strata")
if (length(strat$terms)) terms2 <- terms2[-strat$terms]
X <- model.matrix(terms2, mf)
X
}
tdata <- data.frame(y= 1:10, zed = 1:10, grp =
factor(c(1,1,1,2,2,2,1,1,3,3)))
testfun(y ~ zed*grp, tdata)
testfun(y ~ strata(grp)*zed, tdata)
Terry T.
----- original message --
For building design matrices for Cox proportional hazards models in the
cph function in the rms package I have always used this construct:
Terms <- terms(formula, specials=c("strat", "cluster",
"strata"), data=data)
specials <- attr(Terms, 'specials')
stra <- specials$strat
Terms.ns <- Terms
if(length(stra)) {
temp <- untangle.specials(Terms.ns, "strat", 1)
Terms.ns <- Terms.ns[- temp$terms] #uses [.terms function
}
X <- model.matrix(Terms.ns, X)[, -1, drop=FALSE]
The Terms.ns logic removes stratification factor "main effects" so
that
if a stratification factor interacts with a non-stratification factor,
only the interaction terms are included, not the strat. factor main
effects. [In a Cox PH model stratification goes into the nonparametric
survival curve part of the model].
Lately this logic quit working; model.matrix keeps the unneeded main
effects in the design matrix. Does anyone know what changed in R that
could have caused this, and possibly a workaround?
-------