gof: *g*oodness *o*f *f*it test for a 'coxph' object In survMisc: Miscellaneous Functions for Survival Data

Description

goodness of fit test for a `coxph` object

Usage

 ```1 2 3 4``` ```gof(x, ...) ## S3 method for class 'coxph' gof(x, ..., G = NULL) ```

Arguments

 `x` An object of class `coxph` `...` Additional arguments (not implemented) `G` Number of groups into which to divide risk score. If `G=NULL` (the default), uses closest integer to G = max(2, min(10, ne/40)) where ne is the number of events overall.

Details

In order to verify the overall goodness of fit, the risk score r[i] for each observation i is given by

r[i] = B.X[i]

where B is the vector of fitted coefficients and X[i] is the vector of predictor variables for observation i.
This risk score is then sorted and 'lumped' into a grouping variable with G groups, (containing approximately equal numbers of observations).
The number of observed (e) and expected (exp) events in each group are used to generate a Z statistic for each group, which is assumed to follow a normal distribution with Z \sim N(0,1).
The indicator variable `indicG` is added to the original model and the two models are compared to determine the improvement in fit via the likelihood ratio test.

Value

A `list` with elements:

 `groups` A `data.table` with one row per group G. The columns are nNumber of observations eNumber of events expNumber of events expected. This is exp = ∑ e_i - M_i where e_i are the events and M_i are the martingale residuals for each observation i zZ score, calculated as Z = (e - exp) / exp^0.5 pp-value for Z, which is p = 2 * pnorm(-|z|) where `pnorm` is the normal distribution function with mean 0 and standard deviation 1 and |z| is the absolute value. `lrTest` Likelihood-ratio test. Tests the improvement in log-likelihood with addition of an indicator variable with G-1 groups. This is done with `survival:::anova.coxph`. The test is distributed as chi-square with G-1 degrees of freedom

Note

The choice of G is somewhat arbitrary but rarely should be > 10.
As illustrated in the example, a larger value for G makes the Z test for each group more likely to be significant. This does not affect the significance of adding the indicator variable `indicG` to the original model.

The Z score is chosen for simplicity, as for large sample sizes the Poisson distribution approaches the normal. Strictly speaking, the Poisson would be more appropriate for e and exp as per Counting Theory.
The Z score may be somewhat conservative as the expected events are calculated using the martingale residuals from the overall model, rather than by group. This is likely to bring the expected events closer to the observed events.

This test is similar to the Hosmer-Lemeshow test for logistic regression.

Source

Method and example are from:
May S, Hosmer DW 1998. A simplified method of calculating an overall goodness-of-fit test for the Cox proportional hazards model. Lifetime Data Analysis 4(2):109–20. Springer (paywall)

References

Default value for G as per:
May S, Hosmer DW 2004. A cautionary note on the use of the Gronnesby and Borgan goodness-of-fit test for the Cox proportional hazards model. Lifetime Data Analysis 10(3):283–91. Springer (paywall)

Changes to the `pbc` dataset in the example are as detailed in:
Fleming T, Harrington D 2005. Counting Processes and Survival Analysis. New Jersey: Wiley and Sons. Chapter 4, section 4.6, pp 188. Wiley (paywall)

Examples

 ``` 1 2 3 4 5 6 7 8 9 10 11``` ```data("pbc", package="survival") pbc <- pbc[!is.na(pbc\$trt), ] ## make corrections as per Fleming pbc[pbc\$id==253, "age"] <- 54.4 pbc[pbc\$id==107, "protime"] <- 10.7 ### misspecified; should be log(bili) and log(protime) instead c1 <- coxph(Surv(time, status==2) ~ age + log(albumin) + bili + edema + protime, data=pbc) gof(c1, G=10) gof(c1) ```

Example output

```Loading required package: survival
\$groups
n  e       exp          z          p
1: 31  4  2.835712  0.6914001 0.48931411
2: 31  4  4.897323 -0.4054799 0.68512478
3: 31  1  7.092041 -2.2875847 0.02216172
4: 31  5  6.979415 -0.7492509 0.45370597
5: 32 11  8.455328  0.8751180 0.38150969
6: 31  8 13.104004 -1.4099673 0.15854935
7: 31 18 13.660470  1.1741129 0.24034982
8: 31 19 11.709924  2.1303697 0.03314110
9: 31 24 19.235577  1.0863199 0.27733748
10: 32 31 37.030206 -0.9909554 0.32170735

\$lrTest
Analysis of Deviance Table
Cox model: response is  Surv(time, status == 2)
Model 1: ~ age + log(albumin) + bili + edema + protime
Model 2: ~ age + log(albumin) + bili + edema + protime + indicG
loglik  Chisq Df P(>|Chi|)
1 -553.84
2 -543.14 21.402  9   0.01098 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

\$groups
n  e      exp          z         p
1: 104 11 16.69725 -1.3942573 0.1632399
2: 104 35 37.18701 -0.3586375 0.7198663
3: 104 79 71.11574  0.9349284 0.3498252

\$lrTest
Analysis of Deviance Table
Cox model: response is  Surv(time, status == 2)
Model 1: ~ age + log(albumin) + bili + edema + protime
Model 2: ~ age + log(albumin) + bili + edema + protime + indicG
loglik  Chisq Df P(>|Chi|)
1 -553.84
2 -549.23 9.2235  2  0.009935 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

survMisc documentation built on May 29, 2017, 3:50 p.m.