[R] predicting from coxph with pspline

[Hjellvik, Vidar]

Hello.

I get a bit confused by the output from the predict function when used
on an object from coxph in combination with p-spline, e.g. 

fit <- coxph(Surv(time1, time2, status)~pspline(x), Data)
predict(fit, newdata=data.frame(x=1:2))

It seems like the output is somewhat independent of the x-values to
predict at. For example x=1:2 gives the same result as x=21:22. Does the
result span the x-values in the original data set so that if I want to
predict a value at x=25, than x=c(min(Data$x), 25, max(Data$x)) would be
the right choice of x in newdata?

Below is some code and the output it produces. (An example with the
linear model "fit <- coxph(Surv(time1, time2, status)~x, Data)" which
gives more intuitive results is also included).

By the way, running the example in S-plus gives the same results in the
p-spline case, and slightly different results in the linear case.

require(survival)
Data <- data.frame(matrix(nc=4, byrow=T, c(
   37.6, 67.3, 0, 17.0,
   40.0, 69.7, 0, 22.7,
   47.5, 75.6, 0, 23.0,
   43.4, 62.9, 0, 23.1,
   43.1, 69.1, 1, 23.1,
   44.5, 73.1, 0, 23.1,
   41.7, 60.6, 0, 23.9,
   47.4, 75.5, 0, 23.9,
   42.9, 70.8, 0, 24.1,
   44.5, 63.8, 0, 24.1,
   44.8, 61.7, 0, 25.2,
   42.1, 55.8, 1, 25.9,
   41.0, 68.9, 1, 26.0,
   43.3, 61.2, 0, 26.2,
   38.0, 66.6, 0, 27.9,
   48.3, 77.9, 0, 28.3,
   42.5, 59.2, 0, 28.7,
   46.2, 74.3, 0, 29.3,
   39.9, 65.2, 1, 30.0,
   42.8, 73.2, 0, 44.0)))

names(Data) <- c("time1","time2","status","x")

print("Linear model")
fit <- coxph(Surv(time1, time2, status)~x, Data)
print(data.frame(x=21, pred=predict(fit, data.frame(x=21))))
print(data.frame(x=c(21,31), pred=predict(fit, data.frame(x=c(21,31)))))
print(data.frame(x=1:2, pred=predict(fit, data.frame(x=1:2))))
print(data.frame(x=1:3, pred=predict(fit, data.frame(x=1:3))))
print(data.frame(x=1:4, pred=predict(fit, data.frame(x=1:4))))
print(data.frame(x=c(1,2,4), pred=predict(fit, data.frame(x=c(1,2,4)))))
print(data.frame(x=c(1,3,4), pred=predict(fit, data.frame(x=c(1,3,4)))))

print("Pspline model")
fit <- coxph(Surv(time1, time2, status)~pspline(x), Data)
#print(data.frame(x=21, pred=predict(fit, data.frame(x=21)))) doesn't
work
print(data.frame(x=c(21,31), pred=predict(fit, data.frame(x=c(21,31)))))
print(data.frame(x=1:2, pred=predict(fit, data.frame(x=1:2))))
print(data.frame(x=1:3, pred=predict(fit, data.frame(x=1:3))))
print(data.frame(x=1:4, pred=predict(fit, data.frame(x=1:4))))
print(data.frame(x=c(1,2,4), pred=predict(fit, data.frame(x=c(1,2,4)))))
print(data.frame(x=c(1,3,4), pred=predict(fit, data.frame(x=c(1,3,4)))))



Output:
[1] "Linear model"
   x       pred
1 21 0.03244105
   x        pred
1 21  0.03244105
2 31 -0.03276709
  x      pred
1 1 0.1628573
2 2 0.1563365
  x      pred
1 1 0.1628573
2 2 0.1563365
3 3 0.1498157
  x      pred
1 1 0.1628573
2 2 0.1563365
3 3 0.1498157
4 4 0.1432949
  x      pred
1 1 0.1628573
2 2 0.1563365
3 4 0.1432949
  x      pred
1 1 0.1628573
2 3 0.1498157
3 4 0.1432949

[1] "Pspline model"
   x      pred
1 21 -4.023847
2 31 -4.012924
  x      pred
1 1 -4.023847
2 2 -4.012924
  x      pred
1 1 -4.023847
2 2  2.476954
3 3 -4.012924
  x      pred
1 1 -4.023847
2 2  1.537351
3 3  9.559708
4 4 -4.012924
  x      pred
1 1 -4.023847
2 2  1.537351
3 4 -4.012924
  x      pred
1 1 -4.023847
2 3  9.559708
3 4 -4.012924

Best regards,

Vidar Hjellvik
Norwegian Institute of Public Health
P.O.Box 4404 Nydalen
N-0403 Oslo
+47 21 07 82 66
vidar.hjellvik at fhi.no