[R] survest() for cph() in Design package

[array chip]

Hi, I am trying to run a conditional logistic model on a nested case-control 
study using cph() and then estimate survival based on the model. The data came 
from Prof Bryan Langholz website where he also has the SAS code to this, so I am 
trying to replicate the SAS results.

The data attached. Basically, the variables are:

rstime: risk set age
rsentry: fake entry time, just before rstime
setno: risk set identifier
cc: lung cancer death (0/1)
cr500: independent variable 1 (0/1), cumulative radon > 500WLM
smoke25: independent variable 2 (0/1), 1/2+ pack/dat at age 25
logw: weight to be used as an offset in the model
rstime2: indicator variable to indicate whether age is below 50, 50-69 or >=70, 
used as strata in the model

The SAS code is (if this helps to address my question/problem)

proc phreg data=uminers;
  model rstime*cc(0)=cr500 smoke25/ entry=rsentry offset=logw rl;
baseline out=ch covariates=covs cumhaz=cumhaz stdcumhaz=secumhaz lowercumhaz=lci 
uppercumhaz=uci/ nomean;
  strata rstime2;
run;


My R code is:

library(survival)
library(Design)

> uminers<-read.table("uminers.txt",sep='\t',header=T,row.names=NULL)
> fit <- cph(Surv(rsentry,rstime,cc)~cr500+smoke25+strat(rstime2) +offset(logw), 
>uminers, x=T, y=T)
> fit
Cox Proportional Hazards Model

cph(formula = Surv(rsentry, rstime, cc) ~ cr500 + smoke25 + strat(rstime2) + 
offset(logw), data = uminers, x = T, y = T)

       Obs     Events Model L.R.       d.f.          P 
       774        258      89.98          2          0 
     Score    Score P         R2 
     82.61          0      0.211 

           Status
Stratum     No Event Event
  rstime2=1      100    50
  rstime2=2      348   174
  rstime2=3       68    34

         coef se(coef)    z        p
cr500   1.491    0.184 8.08 6.66e-16
smoke25 0.461    0.182 2.53 1.16e-02

These results are the same as what SAS reported!

Next, I want to estimate the survival for stratum age 50-69 (rstime2=2), so I 
used the following, this is where I got the error message:

> survest(fit, expand.grid(cr500=c(0,1),smoke25=c(0,1),rstime2=2, logw=0), 
>times=unique(uminers$rstime[uminers$rstime2=='2']), conf.int=.95)

Error in factor(rep(1:nstrat, scount), labels = names(fit$strata)) : 
  invalid labels; length 3 should be 1 or 2

Anyone has any suggestions what went wrong?

Thanks!

John


      
-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: uminers.txt
URL: <https://stat.ethz.ch/pipermail/r-help/attachments/20110307/2f5d748c/attachment.txt>