I am fitting a cox PH model w/ 2 predictors, x1 = 0/1 treatment variable
and x2=continuous variable. I am using natural splines (ns) to model
the effect of x2.
I would like to examine the estimated effect of x2 on the hazard. I
have tried various approaches (below; let model.fit= fitted model using
coxph in survival library):
1. The simplest method appears to be using termplot(model.fit). This
appears to work fine as long as I include the treatment variable in the
model. However, without the treatment effect in the model termplot
returns a "?" prompt without plotting the effect of x2. This happens
whenever I consider any model containing only one predictor variable
modeled using ns(). Any suggestions? Also, I am presuming that these
plots indicate the effect of x2 averaged over the effect of x1 (is this
correct?). Ultimately, I would like to be able to produce similar plots
for both x1=0 and x1=1.
2. Using predict(model.fit, newdata, type="lp", safe=T)...as far as I
can tell, this does not appear to give results that are consistent w/
termplot. For my model, the effect of x1 is VERY small (not
statistically significant) and when I overlay the results w/ those
produced by termplot (using "lines") they do not line up at all:
# Predict linear predictor vs. x2 for x1 =1
newdata<-data.frame(x1=rep(1,60), x2=seq(0,30, length=60))
newpred<-predict(model.fit, newdata, type="lp", safe=T)
termplot(model.fit)
lines(newdata[,2], newpred)
Interestingly enough, predict(model.fit) does give back the correct
values for the actual data set used in the fitting:
max(predict(model.fit)-model.fit$linear.predictors)=0. Am I missing
something here?
3. Using the fitted coeficients:
# Coefficients
fitc<-coef(model.fit)
# Predictors
basis <- ns(x2, df = 3) ; # df= 3 were used to fit the model
newx2<- seq(0,30,length=60)
# new data in the coords of the basis and x1=1 for all obs
newdata2<-cbind(rep(1,60),predict(basis, newx2))
newpred<-newdata2%*%fitc
termplot(model.fit, ylim=c(0,3))
lines(newx2, newpred)
This method gives predictions that appear to be proportional to the
results in termplot - but all predicted values are higher. I am missing
something here?
Ideally, I would like to use this method as it appears the most
flexible -allowing one to obtain the linear predictors for various
combinations of x1 and x2.
4. Using Frank Harrell's Design library and the cph function and
plot.Design. This appears to work - giving results that agree w/
termplot (although the results are slightly different because of the
different basis functions used). However, I am more comfortable
(currently) using the coxph function than cph and have had problems w/
cph when trying to fit more complicated (conditional) models w/ multiple
obs per subject.
Any help would be greatly appreciated. I am using R version 1.7.1 on
Windows 2000. I would be glad to share the data set and code directly if
it would be helpful.
John
John Fieberg, Ph.D.
Wildlife Biometrician, MN DNR
5463-C W. Broadway
Forest Lake, MN 55434
Phone: (651) 296-2704
