# 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

> 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.