> -----Original Message-----
> From: [EMAIL PROTECTED] on behalf of array chip
> Sent: Fri 5/9/2008 2:54 PM
> To: [EMAIL PROTECTED]
> Subject: [R] how to check linearity in Cox regression
> Hi, I am just wondering if there is a test available for testing if a linear
> fit of an independent variable in a Cox regression is enough? Thanks for any
> suggestions.
> John Zhang
> require("survival")
> require("splines")
>
> #
> # You could fit a model with the covariate of interest,
> # then an additional fit with other functional forms.
> # I'll use the stanford2 data and age as the covariate.
> # The first fit will have age as a linear term.
> # The second fit will have a quadratic term as well.
> #
>
> stanford2df <- stanford2
> stanford2df$agectr <- stanford2$age - mean(stanford2$age)
> stanford2df$age2 <- stanford2df$agectr^2
> fit1 <- coxph(Surv(time, status) ~ agectr, data = stanford2df)
> fit2 <- coxph(Surv(time, status) ~ agectr + age2, data = stanford2df)
> summary(fit1)
Call:
coxph(formula = Surv(time, status) ~ agectr, data = stanford2df)
n= 184
coef exp(coef) se(coef) z p
agectr 0.0292 1.03 0.0106 2.74 0.0061
exp(coef) exp(-coef) lower .95 upper .95
agectr 1.03 0.971 1.01 1.05
Rsquare= 0.044 (max possible= 0.996 )
Likelihood ratio test= 8.27 on 1 df, p=0.00403
Wald test = 7.51 on 1 df, p=0.00613
Score (logrank) test = 7.6 on 1 df, p=0.00583
> summary(fit2)
Call:
coxph(formula = Surv(time, status) ~ agectr + age2, data = stanford2df)
n= 184
coef exp(coef) se(coef) z p
agectr 0.04100 1.04 0.010237 4.01 6.2e-05
age2 0.00197 1.00 0.000689 2.86 4.2e-03
exp(coef) exp(-coef) lower .95 upper .95
agectr 1.04 0.960 1.02 1.06
age2 1.00 0.998 1.00 1.00
Rsquare= 0.08 (max possible= 0.996 )
Likelihood ratio test= 15.4 on 2 df, p=0.00046
Wald test = 17.4 on 2 df, p=0.00017
Score (logrank) test = 17.3 on 2 df, p=0.000175
> anova(fit1, fit2, test = "Chisq")
Analysis of Deviance Table
Model 1: Surv(time, status) ~ agectr
Model 2: Surv(time, status) ~ agectr + age2
Resid. Df Resid. Dev Df Deviance P(>|Chi|)
1 183 1017.94
2 182 1010.84 1 7.10 0.01
> # the likelihood ratio test suggests that a linear term for
> # age alone may not be adequate. The quadratic form yields
> # a better fit.
>
> # You can also examine functional forms using spline fits.
> fit0 <- coxph(Surv(time, status) ~ age, data = stanford2df)
> fit3 <- coxph(Surv(time, status) ~ ns(age, df = 4), data = stanford2df)
> pred3 <- predict(fit3, type="terms", se=TRUE)
>
> hfit <- pred3$fit[,1]
> hse <- pred3$se[,1]
> hmat=cbind(hfit, hfit+2*hse,hfit-2*hse)
> o <- order(stanford2$age)
> matplot(stanford2$age[o], hmat[o, ], pch="*",col=c(2,4,4), xlab = "age",
+ ylab="Log Relative Risk",main="Cox Model: Survival",type="l")
> # The plot of the spline fit for age shows a non-linear form
>
> anova(fit0, fit3, test = "Chisq")
Analysis of Deviance Table
Model 1: Surv(time, status) ~ age
Model 2: Surv(time, status) ~ ns(age, df = 4)
Resid. Df Resid. Dev Df Deviance P(>|Chi|)
1 183 1017.94
2 180 1009.45 3 8.49 0.04
> # The likelihood ratio test suggests the linear model
> # can be improved upon.
Hope this helps
Steve McKinney
____________________________________________________________________________________
[[elided Yahoo spam]]
______________________________________________
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.
______________________________________________
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.
