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
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)
______________________________________________
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
