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>
