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
[EMAIL PROTECTED]
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
