Dear list members,
I'm unable to fit a parametric survival regression using survreg() in the
survival package with data in "counting-process" ("long") form.
To illustrate using a scaled-down problem with 10 subjects (with data placed on
the web):
--------------- snip ------------
> library(survival)
> RW <- read.table("http://socserv.mcmaster.ca/jfox/.Pickup/RW.txt";)
> RL <- read.table("http://socserv.mcmaster.ca/jfox/.Pickup/RL.txt";)
> RW # "wide" data
week arrest age
1 20 1 27
2 17 1 18
3 25 1 19
4 52 0 23
5 52 0 19
6 52 0 24
7 23 1 25
8 52 0 21
9 52 0 22
10 52 0 20
> head(RL, 20) # "long" data, counting-process form
start stop arrest.time age
1.1 0 1 0 27
1.2 1 2 0 27
1.3 2 3 0 27
1.4 3 4 0 27
1.5 4 5 0 27
1.6 5 6 0 27
1.7 6 7 0 27
1.8 7 8 0 27
1.9 8 9 0 27
1.10 9 10 0 27
1.11 10 11 0 27
1.12 11 12 0 27
1.13 12 13 0 27
1.14 13 14 0 27
1.15 14 15 0 27
1.16 15 16 0 27
1.17 16 17 0 27
1.18 17 18 0 27
1.19 18 19 0 27
1.20 19 20 1 27
--------------- snip ------------
I have no trouble fitting a Cox model to both the wide and long forms of the
data, obtaining (as should be the case) identical results:
--------------- snip ------------
> coxph(Surv(week, arrest) ~ age, data=RW) # works
Call:
coxph(formula = Surv(week, arrest) ~ age, data = RW)
coef exp(coef) se(coef) z p
age 0.0963 1.1011 0.2073 0.46 0.64
Likelihood ratio test=0.21 on 1 df, p=0.643
n= 10, number of events= 4
> coxph(Surv(start, stop, arrest.time) ~ age, data=RL) # works, same
Call:
coxph(formula = Surv(start, stop, arrest.time) ~ age, data = RL)
coef exp(coef) se(coef) z p
age 0.0963 1.1011 0.2073 0.46 0.64
Likelihood ratio test=0.21 on 1 df, p=0.643
n= 397, number of events= 4
--------------- snip ------------
But when I try to fit a parametric survival regression with survreg(), I get an
error with the long form of the data:
--------------- snip ------------
> survreg(Surv(week, arrest) ~ age, data=RW) # works
Call:
survreg(formula = Surv(week, arrest) ~ age, data = RW)
Coefficients:
(Intercept) age
6.35386771 -0.08982624
Scale= 0.7363196
Loglik(model)= -22.1 Loglik(intercept only)= -22.2
Chisq= 0.3 on 1 degrees of freedom, p= 0.58
n= 10
> survreg(Surv(start, stop, arrest.time) ~ age, data=RL) # fails
Error in survreg(Surv(start, stop, arrest.time) ~ age, data = RL) :
Invalid survival type
--------------- snip ------------
I expect that there's something about survreg() that I'm missing. I first noted
this problem quite some time ago but didn't look into it carefully because I
didn't really need to use survreg().
Any help would be appreciated.
Thanks,
John
-----------------------------
John Fox, Professor
McMaster University
Hamilton, Ontario
Canada L8S 4M4
Web: socserv.mcmaster.ca/jfox
______________________________________________
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.
