Re: [R] how to check linearity in Cox regression

From: Steven McKinney <smckinney_at_bccrc.ca>
Date: Fri, 09 May 2008 18:33:21 -0700

> -----Original Message-----
> From: r-help-bounces_at_r-project.org on behalf of array chip
> Sent: Fri 5/9/2008 2:54 PM
> To: r-help_at_stat.math.ethz.ch
> 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_at_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_at_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. Received on Sat 10 May 2008 - 01:39:22 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Sat 10 May 2008 - 03:30:37 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.

list of date sections of archive