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? -------