[R] Help w/ termplot & predict.coxph/ns

[John Fieberg]

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