Sebastian Pölsterl
2015-Aug-06 13:44 UTC
[R] Problem with survfit when model has interaction with strata
Dear all, I encountered a problem when running survfit using a coxph model that contains an interaction with a strata variable (see the attached example). I hope this is the right place to report problems. > source("~/survfit-example.R") Error in scale.default(x2, center = xcenter, scale = FALSE) : length of 'center' must equal the number of columns of 'x' > traceback() 9: stop("length of 'center' must equal the number of columns of 'x'") 8: scale.default(x2, center = xcenter, scale = FALSE) 7: scale(x2, center = xcenter, scale = FALSE) 6: survfit.coxph(fit, newdata = newd) 5: survfit(fit, newdata = newd) at survfit-example.R#10 4: eval(expr, envir, enclos) 3: eval(ei, envir) 2: withVisible(eval(ei, envir)) 1: source("~/survfit-example.R") > sessionInfo() R version 3.2.1 (2015-06-18) Platform: x86_64-unknown-linux-gnu (64-bit) Running under: Fedora 22 (Twenty Two) locale: [1] C attached base packages: [1] stats graphics grDevices utils datasets methods base other attached packages: [1] survival_2.38-3 foreign_0.8-65 loaded via a namespace (and not attached): [1] tools_3.2.1 splines_3.2.1 I believe the problem is occurring on line 263 of survfit.coxph.R, where a model.matrix from newdata is created. The values passed to scale look like: > x2 RXplacebo RXtreatment:strata(SEX)male RXplacebo:strata(SEX)male 1 0 1 0 2 0 0 0 > xcenter RXplacebo strata(SEX)male:RXplacebo 0.5000000 0.2380952 Clearly, model.matrix is adding an additional (redundant) column that was omitted when fitting the original model. Since interactions with strata variables are a special cases in coxph, I believe similar care should be taken when creating the model matrix for newdata. Best regards, Sebastian P?lsterl -------------- next part -------------- library(foreign) library(survival) data <- data.frame(read.spss("http://web1.sph.emory.edu/dkleinb/allDatasets/surv2datasets/anderson.sav")) S <- with(data, Surv(SURVT, I(STATUS == "relapse"))) fit <- coxph(S ~ strata(SEX)*RX, data=data) newd <- data.frame(RX=c("treatment", "treatment"), SEX=c("male", "female")) sfit <- survfit(fit, newdata=newd)