[R] coxph fails to survfit

[Bond, Stephen]

Responding to the suggestion by D Winsemius :

> s1 <- survfit(mod1,newdata=inc[50050:50100,c("strt","stp","incpost", "amt","rate", "termfac")],
+               se.fit=F,individual=T,type="aa")
Error in Surv(time = strt, time2 = stp, event = (resp == 1)) : 
  object 'resp' not found

it appears it wants to fit a survival curve instead of predicting a survival curve for the subject using the explanatory vars. I think I probably need 
predict.coxph(...,type="expected")
or write code for matching time dependent risk to the proper time index of the survival curve. 

Thanks David.

Stephen 

-----Original Message-----
From: David Winsemius [mailto:dwinsemius at comcast.net] 
Sent: Thursday, February 03, 2011 3:10 PM
To: Bond, Stephen
Cc: r-help at r-project.org
Subject: Re: [R] coxph fails to survfit


On Feb 3, 2011, at 2:14 PM, Bond, Stephen wrote:

> I have a model with quant vars only and the error message does not  
> make sense:
>
> (mod1 <- coxph(Surv(time=strt,time2=stp,event=(resp==1))~ +incpost 
> +I(amt/1e5)+rate+strata(termfac),
>               subset=dt<"2010-08-30", data=inc,method="efron"))
> Call:
> coxph(formula = Surv(time = strt, time2 = stp, event = (resp ==
>    1)) ~ +incpost + I(amt/1e+05) + rate + strata(termfac), data = inc,
>    subset = dt < "2010-08-30", method = "efron")
>
>
>                coef exp(coef) se(coef)      z       p
> incpost       0.2563     1.292  0.02479  10.34 0.0e+00
> I(amt/1e+05) -0.0532     0.948  0.00487 -10.92 0.0e+00
> rate         -0.0507     0.951  0.00945  -5.36 8.2e-08
>
> Likelihood ratio test=295  on 3 df, p=0  n= 1192634
>
>> length(mod1$xlevels)
> [1] 0
>
> # now calling survfit with just a few rows from the data
>> s1 <-  
>> survfit 
>> (mod1,newdata=inc[50050:50100,],se.fit=F,individual=T,type="aa")
> Error in `contrasts<-`(`*tmp*`, value = "contr.treatment") :
>  contrasts can be applied only to factors with 2 or more levels

I obviously cannot test my theory, but would have thought the call  
would be:

  s1 <- survfit(mod1, newdata=inc[50050:50100,], c("incpost", "amt,  
"rate", "termfac") ],
                      se.fit=F,individual=T,type="aa")

I.e., a data.frame for newdata that matched the variables used in the  
RHS of the formula.

>
>> inc[50050:50100,]
>              dt      inc incpost strt stp resp rate    amt p1 matfac  
> termfac      vol
> 50050 2009-02-05  0.00000    0.00    0   1    0 5.35 266833  C       
> 5       3 14229571
> 50051 2009-02-06 -0.09575    0.00    1   2    0 5.35 266833  C       
> 5       3 14229571
> 50052 2009-02-09 -0.03875    0.00    4   5    0 5.35 266833  C       
> 5       3 14229571
> 50053 2009-02-10  0.00850    0.00    5   6    0 5.35 266833  C       
> 5       3 14229571
> 50054 2009-02-11  0.00250    0.00    6   7    0 5.35 266833  C       
> 5       3 14229571
> 50055 2009-02-12 -0.04725    0.00    7   8    0 5.35 266833  C       
> 5       3 14229571
> 50056 2009-02-13 -0.00075    0.00    8   9    0 5.35 266833  C       
> 5       3 14229571
> 50057 2009-02-16 -0.10425    0.00   11  12    0 5.35 266833  C       
> 5       3 14229571
> 50058 2009-02-17 -0.10425    0.00   12  13    0 5.35 266833  C       
> 5       3 14229571
> 50059 2009-02-18 -0.02925    0.00   13  14    0 5.35 266833  C       
> 5       3 14229571
> 50060 2009-02-19 -0.01300    0.00   14  15    0 5.35 266833  C       
> 5       3 14229571
> 50061 2009-02-20 -0.04075    0.00   15  16    0 5.35 266833  C       
> 5       3 14229571
> 50062 2009-02-23 -0.03100    0.00   18  19    0 5.35 266833  C       
> 5       3 14229571
> 50063 2009-02-24  0.01975    0.00   19  20    0 5.35 266833  C       
> 5       3 14229571
> 50064 2009-02-25  0.13050    0.00   20  21    0 5.35 266833  C       
> 5       3 14229571
> 50065 2009-02-26  0.12725    0.00   21  22    0 5.35 266833  C       
> 5       3 14229571
>
> ... further output truncated
>
> Please, comment if you see what's the issue.
> Thank you.
> Stephen
>
>
> 	[[alternative HTML version deleted]]
>
> ______________________________________________
> R-help at r-project.org mailing list
> 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.

David Winsemius, MD
West Hartford, CT