1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
|
% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/gof.R
\name{gof}
\alias{gof}
\alias{gof.coxph}
\title{\bold{g}oodness \bold{o}f \bold{f}it test for a \code{coxph} object}
\source{
Method and example are from: \cr
May S, Hosmer DW 1998.
A simplified method of calculating an overall goodness-of-fit test
for the Cox proportional hazards model.
\emph{Lifetime Data Analysis} \bold{4}(2):109--20.
\doi{10.1023/A:1009612305785}
}
\usage{
gof(x, ...)
\method{gof}{coxph}(x, ..., G = NULL)
}
\arguments{
\item{x}{An object of class \code{coxph}}
\item{...}{Additional arguments (not implemented)}
\item{G}{Number of \bold{g}roups into which to divide risk score.
If \code{G=NULL} (the default), uses closest integer to
\deqn{G = \max(2, \quad \min(10, \quad \frac{ne}{40}))}{
G = max(2, min(10, ne/40))}
where \eqn{ne} is the number of events overall.}
}
\value{
A \code{list} with elements:
\item{groups}{A \code{data.table} with one row per group \eqn{G}.
The columns are \describe{
\item{n}{Number of observations}
\item{e}{Number of events}
\item{exp}{Number of events expected. This is
\deqn{exp = \sum e_i - M_i}
where \eqn{e_i} are the events and
\eqn{M_i} are the martingale residuals
for each observation \eqn{i}}
\item{z}{\eqn{Z} score, calculated as
\deqn{ Z = \frac{e - exp}{\sqrt{exp}}}{
Z = (e - exp) / exp^0.5}
}
\item{p}{\eqn{p}-value for \eqn{Z}, which is
\deqn{ p = 2. \code{pnorm}(-|z|)}{
p = 2 * pnorm(-|z|)}
where \code{pnorm} is the normal distribution function
with mean \eqn{\mu =0}{0} and standard deviation \eqn{\sigma =1}{1}
and \eqn{|z|} is the absolute value.}
}}
\item{lrTest}{Likelihood-ratio test.
Tests the improvement in log-likelihood with addition
of an indicator variable with \eqn{G-1} groups.
This is done with \code{survival:::anova.coxph}.
The test is distributed as chi-square with \eqn{G-1} degrees of freedom}
}
\description{
\bold{g}oodness \bold{o}f \bold{f}it test for a \code{coxph} object
}
\details{
In order to verify the overall goodness of fit,
the risk score \eqn{r_i}{r[i]} for each observation \eqn{i} is given by
\deqn{r_i = \hat{\beta} X_i}{r[i] = B.X[i]}
where \eqn{\hat{\beta}}{B} is the vector of fitted coefficients
and \eqn{X_i}{X[i]} is the vector of predictor variables for
observation \eqn{i}.
\cr
This risk score is then sorted and 'lumped' into
a grouping variable with \eqn{G} groups,
(containing approximately equal numbers of observations).
\cr
The number of observed (\eqn{e}) and expected (\eqn{exp}) events in
each group are used to generate a \eqn{Z} statistic for each group,
which is assumed to follow a normal distribution with
\eqn{Z \sim N(0,1)}.
\cr
The indicator variable \code{indicG} is added to the
original model and the two models are compared to determine the
improvement in fit via the likelihood ratio test.
}
\note{
The choice of \eqn{G} is somewhat arbitrary but rarely should
be \eqn{> 10}.
\cr
As illustrated in the example, a larger value for
\eqn{G} makes the \eqn{Z} test for each group more likely to be significant.
This does \emph{not} affect the significance of adding the
indicator variable \code{indicG} to the original model.
\cr \cr
The \eqn{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 \eqn{e} and \eqn{exp}{exp} as
per Counting Theory.
\cr
The \eqn{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.
\cr \cr
This test is similar to the Hosmer-Lemeshow test for logistic regression.
}
\examples{
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)
}
\references{
Default value for \eqn{G} as per: \cr
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.
\emph{Lifetime Data Analysis} \bold{10}(3):283--91.
\doi{10.1023/B:LIDA.0000036393.29224.1d}
Changes to the \code{pbc} dataset in the example are as detailed in: \cr
Fleming T, Harrington D 2005.
\emph{Counting Processes and Survival Analysis}.
New Jersey: Wiley and Sons. Chapter 4, section 4.6, pp 188.
\doi{10.1002/9781118150672}
}
|
