Abstract
Reuse of controls from nested case-control designs can increase efficiency in many situations, for instance with competing risks or in other multiple endpoints situations. The matching between cases and controls must be broken when controls are to be used for other endpoints. A weighted analysis can then be performed to take care of the biased sampling from the cohort. We present the R package multipleNCC for reuse of controls in nested case-control studies by inverse probability weighting of the partial likelihood. The package handles right-censored, left-truncated and additionally matched data, and varying numbers of sampled controls and the whole analysis is carried out using one simple command. Four weight estimators are presented and variance estimation is explained. The package is illustrated by analyzing health survey data from three counties in Norway for two causes of death: cardiovascular disease and death from alcohol abuse, liver disease, and accidents and violence. The data set is included in the package.The nested case-control (NCC) design (Thomas 1977) (incidence density sampling/risk set sampling) is popular within epidemiology due to its cost-effectiveness and efficiency. At each event time, \(m\) controls are sampled from the subjects at risk. The controls must be event-free at the time their case experienced the event. They might also be matched to the cases on additional factors. Covariates are then obtained for cases and sampled controls.
The analysis of nested case-control data has traditionally been carried out using stratified Cox-regression, where the stratification is with respect to matched case-control sets. That method, however, does not allow for breaking the matching, i.e., analyze the data without directly considering the matched sets. Hence, the sampled controls cannot be used for other cases than they originally were sampled for. In many situations one may want to reuse the controls, for instance when analyzing a subset of the original cases or when there exist more than one endpoint of interest. A few examples of such studies are Hultman et al. (1999; Parsonnet et al. 1991; Floderus et al. 1993; Øyen et al. 1997; Tynes and Haldorsen 1997; Hankinson et al. 1998; Grimsrud et al. 2002; Levine et al. 2004; Clendenen et al. 2011; Meyer et al. 2013).
We assume a competing risks situation in this paper for ease of
presentation. Thus there are different types of endpoints, and the
subjects can experience at most one of them. A special instance of this
is a setting with only one type of endpoint, the single event situation
is therefore covered by the competing risks situation. The R package multipleNCC
(Nathalie C. Støer and Samuelsen 2016) is
built with competing risks in mind, however it is not limited to such
situations. Even though more complex event history settings would often
call for more advanced multi-state modelling, reuse of controls may be
handled by using multipleNCC together with coxph()
from package survival
(T. M. Therneau 2016; T. M. Therneau and Grambsch
2000), see Section 5.
A method for breaking the matching in NCC designs was introduced by Samuelsen (1997) and further studied by Chen (2001; Samuelsen, Ånestad, and Skrondal 2007; Saarela et al. 2008; Salim et al. 2009; Cai and Zheng 2012; Salim, Yang, and Reilly 2012; N. C. Støer and Samuelsen 2012, 2013; N. C. Støer, Meyer, and Samuelsen 2014). This method bases the estimation on a weighted partial likelihood, thus weighted Cox-regressions are carried out. The weights are inverse sampling probabilities, which must be estimated from the data, and different estimators have been suggested.
Even though it is fairly easy to estimate the weights, it is an extra
step in the analysis. Having a more automatic estimation procedure,
i.e., a one-line call in R with similar syntax as coxph()
will make this way of analyzing NCC data more generally available.
We present the R package multipleNCC in this paper. The
function wpl() estimates weights and carries out weighted
Cox-regressions. The users can choose between four options for weight
estimation. The function handles both right-censored and left-truncated
data and has some possibilities for variance estimation, apart from
robust variances. Additional matching is incorporated for three of the
four weight estimators, and varying number of controls for all four.
This is the first statistical software that performs inverse probability
weighting (IPW) aimed at NCC data.
The outline of the paper is as follows; we introduce the general framework of inverse probability weighting for nested case-control data in Section 2. Then Sections 3 and 4 follow. The package is described in Section 5 and illustrated in Section 6. A comparison between the traditional estimator and the IPW estimators is given in Section 7, followed by Section 8.
We have a cohort consisting of \(n\) subjects, where the \(i\)-th individual is followed from the left-truncation time \(l_i\), which may be zero, to time of event \(\tilde{t}_i\) or time of censoring \(c_i\). Thus the cohort members are followed from \(l_i\) to \(t_i = \min\left(\tilde{t_i},c_i\right)\). There are \(K\) competing endpoints, and each subject can at most experience one of them.
At each event time, \(m\) controls
are sampled for the case experiencing the event and \(m = m\left(t\right)\) may depend on time.
Finally, let us introduce what we call sampling-status indicator \(S_i\), which is required input to
wpl(). This indicator takes values in \(\{0,1,\ldots,K,K+1\}\), zero indicates a
non-sampled subject in the cohort, 1 indicates sampled controls (and not
cases) for any of the endpoints in question, while “\(2,\ldots,K+1\)” indicate cases of one of
\(K\) types.
The NCC design is tightly connected to the Cox proportional hazards model and our model for endpoint \(k\) can be written as \[\begin{aligned} \label{EQ:IPW_mod} h_{ki}\left(t|x_i,z_i\right) = h_{k0}\left(t\right)\exp\left(\beta_k' x_i + \gamma_k' z_i\right). \end{aligned} (\#eq:EQIPW-mod)\] Here \(h_{k0}\left(t\right)\) is the baseline hazard for endpoint \(k\), \(x_i\) are covariates and confounders while \(z_i\) are additional matching variables (additional to matching on time), if additional matching has been carried out, otherwise \(z_i\) will be zero. The \(\beta_k\) and \(\gamma_k\) are the log-hazard ratios connected to \(x\) and \(z\) for the \(k\)-th endpoint.
Since the matching is broken with IPW, it will generally be important to adjust for the matching variables (N. C. Støer and Samuelsen 2013). They are included as linear functions in Equation @ref(eq:EQIPW-mod) for simplicity, however more general models with other types of functions are possible.
The controls in a NCC design are matched to the cases on at risk status and possibly additional factors. Due to this matching, it is not straightforward to reuse the controls for other endpoints/cases since the matching must be broken. Samuelsen (1997) suggested a weighted partial likelihood which resembles the standard Cox-likelihood \[\begin{aligned} \label{EQ:wpl} L_k\left(\beta,\gamma\right) = \prod_j\frac{\exp\left(\beta_k' x_j + \gamma_k' z_j\right)}{\sum_{i \in \mathcal{R}_j}\exp\left(\beta_k' x_i + \gamma_k' z_i\right)w_i}. \end{aligned} (\#eq:EQwpl)\] This likelihood enables the controls to be used for other endpoints and Saarela et al. (2008) and Salim et al. (2009) were the first to discuss this likelihood in connection to competing risks. The product is over all cases of type \(k\), while the sum is over a set \(\mathcal{R}_j\), defined as all cases (of all types) and all controls at risk at \(t_j\). The weight, \(w_i = 1/p_i\), is the inverse probability that individual \(i\) is ever being sampled. This probability will be 1 for cases since all of them are sampled by design, and it must be estimated from the data for the controls. We assume time invariant covariates, although time-dependent covariates are in theory possible as long as they are known at all event times at which the subject is at risk. This has, however, not yet been implemented in package multipleNCC.
The fundamental idea behind inverse probability weighting is to adjust for the biased sample from the cohort. The sample is biased, first and foremost, with respect to the proportion of cases and controls, but with additional matching it can also be biased with respect to matching variables. The idea is then to let each control represent a number of subjects in the cohort by giving them weights larger than 1. The less probable it was for a given subject to be sampled, the more subjects in the cohort it should represent since that “type” of controls likely are under-represented in the NCC sample. By using inverse sampling probabilities as weights, this is accomplished. The analysis is then carried out “as if the data were from a cohort study”, thus by a weighted Cox-regression. This idea was first proposed by Hansen and Hurwitz (1943) with a survey sampling perspective for sampling with replacement, and later generalised to sampling without replacement by Horvitz and Thompson (1952). Inverse probability weighting is also commonly used in the context of missing data (Robins, Rotnitzky, and Zhao 1994).
To increase efficiency and adjust for confounding, the controls in a nested case-control design are often matched on additional factors than at risk status. This can for instance be year of birth, sex or years since first employment. We divide such matching into two groups: category matching and caliper matching (Cochran and Rubin 1973).
With category matching, the controls are required to have the same value on the matching variable as the case, and the matching variable is often a categorical covariate. As an example, the controls can be matched to the cases on sex and a male case is then required to have a male control. With caliper matching, the matching variable is typically continuous and the control’s value of the matching variable must lie within a specified interval around the case’s value. For instance the controls could be matched on year of birth plus/minus 2 years, i.e., the birth year of a control must be within two years of the case’s birth year.
The weights in Equation @ref(eq:EQwpl) must be estimated from the data at hand, and three types of estimators have been considered: Kaplan-Meier (KM) type of weights (Samuelsen 1997; Salim et al. 2009; Cai and Zheng 2012), more model based logistic regression type (Mark and Katki 2006; Samuelsen, Ånestad, and Skrondal 2007; Saarela et al. 2008; N. C. Støer and Samuelsen 2013) and local averaging (Chen 2001), referred to as Chen-weights.
The Kaplan-Meier type of estimator without additional matching can be formulated as \[\begin{aligned} \label{EQ:KMuM} p_i = 1-\prod_{l_i <t_j < t_i}\left\{ 1 - \frac{m}{n\left(t_j\right)-1}\right\}. \end{aligned} (\#eq:EQKMuM)\] Here \(n\left(t_j\right)\) is the number at risk in the cohort at time \(t_j\) and \(m\) the number of sampled controls per case. The estimator in Equation @ref(eq:EQKMuM) resembles the Kaplan-Meier estimator. The KM weights can be generalized to situations with additional matching by taking the product only over subjects that meet the matching criteria, and letting the denominator only consist of subjects at risk that meet the matching criteria. Let \(n_j\left(t_j\right)\) count the subjects at risk at time \(t_j\) who meet the matching criteria of case \(j\). Then the formula with additional matching can be expressed as \[\begin{aligned} \label{EQ:KMM} p_i = 1-\prod_{j} \left\{1-\frac{m}{n_j\left(t_j\right)-1}I\left(\text{Control } i \text{ could be sampled for case } j\right)\right\}. \end{aligned} (\#eq:EQKMM)\]
By replacing \(m\) with \(m\left(t_j\right)\), where \(m\left(t_j\right)\) is the number of sampled controls for the case at \(t_j\), the situation with varying number of controls is covered.
A more model based approach is to use logistic regression models, either traditional logistic regression or the more flexible generalized additive model [GAM; T. J. Hastie and Tibshirani (2009) Chap. 9] with logit-link. The sampling indicator, \(O_i\), is used as outcome and the left-truncation time and censoring time as covariates with \(\xi\) being an intercept term: \[\begin{aligned} \label{glm-weights} p_i = \mathbb{E}[O_i|t_i,l_i] = \frac{\exp\left(\xi+f\left(t_i,l_i \right)\right)}{1+\exp\left(\xi+f\left(t_i,l_i \right)\right)}. \end{aligned} (\#eq:glm-weights)\] It is important to note that this regression is carried out on the cohort excluding cases of all types. The reason for this is that all cases are sampled with a known probability of 1, thus including them in the regression would interfere with the estimation of the sampling probabilities for the controls.
With \(f\left(t_i,l_i\right) = f_1\left(t_i\right)+f_2\left(l_i\right)\) where \(f.\left(.\right)\) are linear functions, the estimator in Equation @ref(eq:glm-weights) is the traditional logistic regression model, and the inverse of those probabilities are referred to as GLM-weights. When \(f.\left(.\right)\) are smooth functions, the result is a GAM-model. In situations with additional matching, the matching variables should also be included in the regression model. Category matched variables are included as categorical factors while caliper matched variables are included as continuous covariates with GLM-weights and as smoothed functions with GAM-weights.
The number of controls for each case is not an explicit part of the estimation procedure for logistic regression weights and therefore no extra care must be taken in situations with time- and endpoint-dependent number of controls.
Before we can introduce the Chen-weights method (Chen 2001), which is also known as local
averaging, we need some additional notation. Let \(0 = l^0 < l^1 < \ldots < l^A\) be
a partition of the range of left-truncation times and \(0 = t^0 < t^1 < \ldots < t^B\) a
partition of the range of follow-up times where \(t^A\) and \(t^B\) are the upper limit of the
left-truncation times and censoring times respectively. We also define
\(\mathcal{J}_a =
\left(l^{a-1},l^a\right]\) and \(\mathcal{I}_b = \left(t^{b-1},t^b\right]\).
The intervals in each direction are taken to be of the same length in
wpl(), however the interval length can in principle vary.
The local averaging weights can then be expressed as \[\begin{aligned}
w_{ab} = \frac{\sum_{i=1}^n I\left(l_i \in \mathcal{J}_a,t_i \in
\mathcal{I}_b,i \text{ is not a case}\right)}{\sum_{i=1}^n I\left(l_i
\in \mathcal{J}_a,t_i \in \mathcal{I}_b,i \text{ is a sampled control
and not a case}\right)}.
\end{aligned}\] All controls included in the study in \(\mathcal{J}_a\) with a censoring time in
\(\mathcal{I}_b\) are given weight
\(w_{ab}\). Hence, all subjects sampled
within the same combination of intervals will be given the same weight.
Samuelsen, Ånestad, and Skrondal (2007)
noted that this amounts to post-stratifying on grouped left-truncation
and censoring times.
Generalizing these weights to handle additional matching would require partitioning of the matching variables in addition to the range of left-truncation and right-censoring times, and this will introduce a large number of combination of intervals. Due to this, the weights will likely be unstable since a small number of subjects would belong to each group. The Chen-weights are therefore only implemented for situations without additional matching.
As with logistic regression weights, the number of controls per case is not an explicit part of the estimation procedure for local averaging. Situations with time- and endpoint-dependent number of controls are therefore handled with no modification of the estimator.
Since the subjects enter the weighted partial likelihood
(Equation @ref(eq:EQwpl)) whenever they are at risk, the likelihood
contributions are not independent and the variance estimation cannot be
based on the inverse of the information matrix only. A simple, yet
sometimes conservative solution (Cai and Zheng
2012), is to use robust variances (Lin and
Wei 1989; Barlow 1994). This is the default option in
wpl(). A variance estimator for the KM-weights can be found
in Samuelsen (1997) when there is no
additional matching. A variance estimator for Chen-weights is given in
Chen (2001), but since Samuelsen, Ånestad, and Skrondal (2007) argues
that this estimator can be considered as a post-stratified case-cohort
estimator, one may apply the variance estimator of Borgan et al. (2000). We have implemented the
variance estimator of Samuelsen (1997),
which exactly accounts for the procedure used to estimate the weights,
and extended it to allow for additional matching, see details below and
web-appendix to Cai and Zheng (2012). For
the Chen-weights, we have implemented the variance estimator based on
post-stratification, but only without additional matching since, as
discussed in Section 3.3, the
approach will be difficult to extend to additional matching.
Samuelsen (1997) showed that the covariance matrix with KM-weights (without additional matching) can be estimated by \[\begin{aligned} \label{EQ:Gammahat} \hat\Gamma={\hat I}^{-1} + {\hat I}^{-1}\hat\Delta {\hat I}^{-1}. \end{aligned} (\#eq:EQGammahat)\] Here \({\hat I}^{-1}\) is the inverse of the information matrix returned from the weighted Cox-regression and \[\begin{aligned} \label{EQ:Deltahat} \hat\Delta=\sum_i U_i U_i^\top \frac{1-p_i}{p_i^2} + \sum_{i \neq j} U_i U_j^\top \frac{\widehat{COV}_{ij}}{p_{ij} p_i p_j}. \end{aligned} (\#eq:EQDeltahat)\] \(U_i\) is the score contribution for individual \(i\) and \(U_i^\top\) its transpose, \(p_{ij}\) the estimated probability that both individual \(i\) and \(j\) were sampled and \(\widehat{COV}_{ij}\) the estimated covariance between sampling indicators for these two individuals. Note that the score contributions \(U_i\) and the \(W_i\)’s considered by Samuelsen (1997) are asymptotically equivalent. The indices \(i\) and \(j\) in the sums run over all non-cases in the study. Let also \(U\) be a matrix consisting of all \(U_i\) for non-cases. It was argued by Samuelsen (1997) that the term \(p_{ij}\) can be replaced by \(p_i p_j\). A formula for \(\widehat{COV}_{ij}\) is given in Samuelsen (1997). We have implemented this variance formula both with and without additional matching, although in the latter case the \(\widehat{COV}_{ij}\) needs modification. If two individuals \(i\) and \(j\) cannot both be sampled for any same case, then these two individuals are sampled independently and the covariance of their sampling indicators equals zero. If they can both be sampled for one or more of the same cases then \(\widehat{COV}_{ij}\) is obtained using Equation (3.1) in Samuelsen (1997), replacing the number at risk at time \(t_j\) by the number at risk \(n_j\left(t_j\right)\) at \(t_j\) satisfying the matching criteria, and taking the product only over the cases where both \(i\) and \(j\) can be sampled as controls.
The covariance matrix estimator in Equation @ref(eq:EQGammahat) has been considered numerically hard to calculate. However, it can be simplified by noting that Equation @ref(eq:EQDeltahat) can be expressed as a matrix product. Let \(R\) be a matrix with the terms \(\left(1-p_i\right)/p_i^2\) for all non-cases along the diagonal and terms \(\widehat{COV}_{ij}/\left(p_i p_j\right)^2\) otherwise. Let furthermore \(U_{\left(l\right)}\) be a vector of the \(l\)-th component of \(U\). Then the component \(\hat\Delta_{kl}\) of \(\hat\Delta\) can be expressed as \(U_{\left(l\right)}^\top R U_{\left(k\right)}\) and taken together this means that with \(U\) being the matrix consisting of all \(U_i\) for the non-cases we get \(\hat\Delta = U^\top R U\). Thus the calculation of the second sum in Equation @ref(eq:EQDeltahat) as a double for-loop is avoided. However, with additional matching and many sampled non-cases \(\widehat{COV}_{ij}\) and the matrix \(R\) may still be computationally demanding.
The multipleNCC package provides one main function:
wpl() which carries out the weighted Cox-regressions. It
calls one of four internal functions, pKM(),
pGAM(), pGLM() and pChen(), for
weight estimation and may call two additional functions,
ModelbasedVar() and PoststratVar(), for
variance estimation. The functions that estimate the sampling
probabilities can be accessed through the wrapper functions
KMprob(), GAMprob(), GLMprob()
and Chenprob(). The function wpl() is used for
estimation of hazard ratios. It has a number of mandatory arguments and
some which are set to default values if not specified by the user, see
Table 1. An important part of the
wpl() function is the weighted Cox-regression which is
carried out by coxph() from the survival
package.
| Argument | Description | Default value |
|---|---|---|
Surv object |
A survival object. | Mandatory |
formula |
A formula object, with the response on the left of a
~ operator, and the terms on the right. The response must
be a survival object. The status variable going into
Surv should have one for cases and zero for controls and
non-sampled subjects. Cohort dimension. |
Mandatory |
data |
A data.frame used for interpreting the
variables named in the formula. Cohort dimension. |
Mandatory |
samplestat |
A vector containing sampling and status information: 0 represents non-sampled subjects in the cohort; 1 represents sampled controls; 2,3,… indicate different events. Cohort dimension. | Mandatory |
m |
Number of sampled controls. A scalar if equal number of
controls for all cases. If unequal number of controls per case: a vector
of length equal to the number of cases. The vector must be in the same
order as the cases in the samplestat vector. |
1 |
weight.method |
One of four weights: "KM",
"gam", "glm", "Chen". |
"KM" |
no.intervals |
Number of intervals for censoring times for Chen-weights with only right-censoring. | 10 |
variance |
"robust", or "Modelbased" for
KM-weights and "Poststrat" for Chen-weights. |
"robust" |
left.time |
Entry time if the survival times are left-truncated. Cohort dimension. | Zero |
no.intervals.left |
Number of intervals for Chen-weights with left-truncation. A vector of the form \[number of intervals for left-truncated time, number of intervals for survival time\]. | \[3, 4\] |
match.var |
If the controls are matched to the cases (on other
variables than time), match.var is the vector or matrix of
matching variables. Cohort dimension. |
Zero |
match.int |
A vector of length \(2
\times\) number of matching variables. For caliper matching
(matched on value plus/minus epsilon) match.int should
consist of c(-epsilon, epsilon). For exact
matching match.int should consist of
c(0,0). |
Zero |
It is important to note that most arguments should have cohort
dimension (Table 1). By that we mean that the
arguments should have the same length as the number of subjects \(n\) in the cohort. Thus, partially known
covariate information must be imputed. These imputed values are not
included in the estimation, hence any values can be chosen, however
NA should be avoided.
The GLM-weights are estimated with the glm() function in
R with family = binomial. For GAM-weights, the
gam() function in the mgcv package
(S. N. Wood 2006; S. Wood 2015) is used,
also with family = binomial. The KM- and Chen-weights are
estimated as explained in Sections 3.1 and
3.3.
For GLM-weights, the matching variables are included in the logistic
regression as continuous covariates with caliper matching and as
categorical covariates for category matching. For GAM-weights the
matching variables are included as smooth functions with caliper
matching and as categorical covariates for category matching. If a
matching variable has many levels, a large number of parameters must be
estimated and some sort of grouping of the levels of the category
matching variable(s) could be sensible. Such grouping must be applied to
the matching variable before it enters wpl().
If the method of variance estimation is not specified, the robust
option is chosen by default. The "Modelbased" option can be
chosen for KM-weights. Using the "Modelbased" option for
other weights will result in an error message. Similarly if the option
"Poststrat" is chosen together with other weights than
"Chen", an error message will be displayed.
The code for fitting a wpl model is
R> wpl(formula, data, samplestat)The formula is a formula object and has the same syntax as the
formula in the coxph() function. The minimal code for
fitting a wpl model is therefore
R> wpl(Surv(survival_time, status) ~ X, data, samplestat)This is a model with one control per case, no additional matching or left-truncation and KM-weights. With left-truncation and GLM-weights it would be
R> wpl(Surv(left_time, survival_time, status) ~ X, data, samplestat,
+ weight.method = "glm")and with additional matching
R> wpl(Surv(left_time, survival_time, status) ~ X + M, data, samplestat,
+ weight.method = "glm", match.var = M, match.int = c(-2, 2))In this model the controls are matched to the cases on only one
additional factor M, which for instance could be year of
birth, plus/minus two years, represented by
match.int = c(-2, 2). Generally, it would be important to
adjust for the additional matching variables, since the matching has
been broken, thus M is also included as an ordinary
covariate. If there are more than one matching variable, for instance
\(M1\) and \(M2\), both of them should be included as
covariates and the match.var should be a matrix consisting
of the \(M1\) and \(M2\). The match.int argument
is still a vector with the intervals in the same order as in
match.var. Thus the syntax for a model with one caliper
matching criterion and one category criterion could be
R> wpl(Surv(left_time, survival_time, status) ~ X + M1 + M2, data, samplestat,
+ weight.method = "gam", match.var = cbind(M1, M2), match.int = c(-2, 2, 0, 0))Note that the interval should be c(0,0) for category
matching. As an example, if the controls were matched on year of birth
plus/minus 2 years, county of residence and years since first employment
plus/minus 6 months, with year of birth and year since first employment
measured in years, match.int = c(-2,2,0,0,-0.5,0.5).
In a traditional analysis of NCC data, subjects appear in the data set as many times as they are sampled. For instance if a subject is first sampled as a control and later itself becomes a case, it appears in the data set first as a control and then a second time as a case. With IPW, the subjects should only appear in the data set once, and it is important to let all subjects who at some point became a case, be a case in the data set.
The event indicator included in the Surv() function
could have a one for cases and a zero for non-cases. However, this event
indicator is not actually used by the program. With only right-censored
data an event indicator is not required by Surv and it can
thus be omitted. When the data is left-truncated the event indicator
must be included and we therefore suggest to always include it even
though it is not used by the program. All information regarding events,
possibly of different types, controls and non-sampled subjects in the
cohort is included through the sampling-status indicator
samplestat. Non-sampled subjects should be given value
zero, sampled controls (who are not cases of any type) should have 1,
while cases of the first type should be given 2, cases of the second
type 3, etc. By default all subjects with non-zero values (except for
the cases in question) are used as controls. If this is not desirable,
the sampling-status indicator can be modified, see Section 5.1.
The estimated weights can be extracted directly from the wpl object
by $weights, it is, however important to note that the data
is sorted by inclusion time inside wpl and the ordering of the weights
is therefore not the same as the ordering of the original data. The
sampling probabilities can also be estimated directly with
KMprob, GAMprob, GLMprob and
Chenprob which also return them in the same order as the
input data. These four functions have similar syntax, although with some
differences in required arguments.
R> KMprob(survtime, samplestat, m, left.time = 0, match.var = 0, match.int = 0)
R> GLMprob(survtime, samplestat, left.time = 0, match.var = 0, match.int = 0)
R> GAMprob(survtime, samplestat, left.time = 0, match.var = 0, match.int = 0)
R> Chenprob(survtime, samplestat, no.intervals = 10, left.time = 0,
+ no.intervals.left = c(3, 4))Some arguments to these functions are mandatory, while some arguments
should only be supplied in given situations. Those arguments are
left.time and no.intervals.left which should
only be included with left-truncated data, and match.var
and match.int only with additionally matched data. All
arguments to the functions above can be found in Table 1, except for survtime which is the
follow-up time. It is important to note that survtime,
samplestat, left.time and
match.var should have length or number of rows equal to the
cohort size \(n\).
When the controls are matched on additional factors it may happen
that there are none, or too few, eligible controls for some cases.
Normal practice is then to widen the matching criteria somewhat for
those particular cases. But by doing this there will be fewer subjects
at risk that meet the original matching criterium than there are sampled
controls. wpl() will therefore print out a warning telling
the user for which case this happens and that the controls for that
particular case are given weight = 1. This is reasonable since those
controls are sampled with probability close or equal to 1.
multipleNCC does not carry out a traditional analysis using a stratified Cox-regression. The main reason for this is that each subject should only appear once in the data set provided to multipleNCC. Without explicit information about the original case-control sets it is impossible to reconstruct the original nested case-control data where each subject appears as many times as they are sampled. Additionally, it is so simple to carry out a traditional nested case-control analysis with a stratified Cox-regression that there is no reason to include it in multipleNCC.
A main endpoint can often be divided into sub-endpoints. For instance
a cancer endpoint can be divided into metastatic and non-metastatic
cancer. Analysis of sub-endpoints using all sampled controls can be
carried out with wpl() by making a new sampling-status
indicator with unique values for each sub-type. So instead of having a
samplestat indicator from 0–2 (\(0=\) non-sampled subjects, \(1=\) sampled controls, \(2=\) cancer cases), it will have values
0–3, where 2 corresponds to metastatic cancer and 3 to non-metastatic
cancer.
Many multiple endpoint applications do not fit into a competing risks
framework. However, reuse of controls can still be of interest. The
solution can then be to estimate the sampling probabilities with
KMprob, GAMprob, GLMprob or
Chenprob and use the inverse of those estimated
probabilities as weights in the ordinary coxph() function
with robust variances. An example of a multiple endpoint situation which
does not fit into the competing risks framework is the case of
subsequent events (N. C. Støer, Meyer, and
Samuelsen 2014). This is a situation with two types of events and
the second type may only occur after the first, e.g. incidence of
prostate cancer and death from prostate cancer. Using all sampled
controls when analyzing the subsequent endpoint may increase the
efficiency substantially.
Inverse probability weighting is a standard approach for handling missing data and case-control data can be considered as data missing by design. Thus the functions for estimating inclusion probabilities described above can be used more generally for addressing other models than proportional hazards. For instance Samuelsen (1997) discussed the possibility of fitting parametric survival models, Suissa, Edwardes, and Boivin (1998) considered estimation of standardized mortality ratios and Cai and Zheng (2012) discussed estimation from estimating equations in general with IPW from NCC studies.
Often software allows for weighting and then it is straightforward to obtain the weighted estimators. Examples are the Nelson-Aalen and Breslow estimators for cumulative hazards or additive hazards models. With respect to variance estimation of weighted estimators we believe that robust variances often will give results that closely match the model based, although theoretically they are conservative (Samuelsen 1997; Cai and Zheng 2012). In general, however, theory has not been carefully developed for such estimators and implementation of robust procedures has not generally been validated, thus care should be taken when interpreting the variance estimates.
Furthermore, we believe that the IPW weights for NCC are in particular useful when considering time to event data with the original cohort follow-up scheme. We are less convinced that the IPW approach is the best choice when for instance estimating the population means of the exposures obtained in the NCC or considering regression models for how such exposures depend on other cohort information.
There exists a number of packages for weighted analysis for different
purposes. Some are aimed at causal inferences such as ipw (van der Wal and Geskus 2011) and MatchIt
(Ho et al. 2011). Others have a missing
data or survey sampling perspective. Two such packages are NestedCohort
(Mark and Katki 2006; Katki and Mark 2008)
and survey
(T. Lumley 2004; T. Lumley 2014). The
survey package was developed for analyses of survey data, but
include functions for two-phase designs. A nested case-control design
can be seen as a two-phase design where the cohort is collected at Phase
1 and the Phase 2 data is the cases and sampled controls with additional
covariate information. The Phase 2 sampling scheme is of course somewhat
more complex than what is usually seen in survey sampling. However, when
Chen-weights are used for the nested case-control design, a stratified
Phase 2 sampling is implicitly assumed, which is a well-known survey
sampling design. Hence, a weighted analysis using Chen-weights can be
carried out by specifying a stratified two-phase design with the
function twophase and carry out the Cox-regression using
this design with svycoxph.
The NestedCohort package is a general package for cohort
analyses with a missing data/two-phase design perspective. The theory
behind it is based to some extent on Robins,
Rotnitzky, and Zhao (1994) and the variance estimators also build
on this paper. Weighted Cox-regressions are carried out with
nested.coxph where a working model is assumed for the
sampling probabilities in Phase 2. This working model is a logistic
regression model for the sampling indicator with all variables
contributing to the missingness as covariates. This is identical to
specifying weight.method = "glm" in wpl.
However, the variance estimator in nested.coxph may be
somewhat better in special situations where the robust variances are
conservative. The NestedCohort package can however only be used
for glm-weights so that the more commonly used KM-weights are not
applicable with this package. The authors of the package also state that
fine matching is problematic with NestedCohort. In the data
example below nested.coxph gave identical estimates and
standard errors as using glm-weights with wpl.
We use a collection of cardiovascular health screenings as our cohort, also used in Aalen, Borgan, and Gjessing (2008). All men and women aged 35–49 from three Norwegian counties: Oppland, Sogn og Fjordane and Finnmark were invited to participate in health screenings from 1974–1978. The screening consisted of a health examination including measurement of height and weight. The participants were also asked to respond to a questionnaire which among other things contained questions regarding smoking status.
More than 90% of the invited subjects chose to participate which
resulted in a cohort of about 50 000 subjects. The cohort was linked to
the Causes of Death Registry kept by Statistics Norway and followed up
for deaths until the end of year 2000. Since there were few subjects at
risk younger than 40 years, the survival times were left-truncated at
age 40 or age at health screening. The left-truncation time is named
agestart in the data. The survival time is named
agestop and is either the age at death or censoring.
The data set included in multipleNCC and used as
illustration here, is a random sample of 3933 subjects from the cohort.
It is a cohort study, but we have carried out a synthetic nested
case-control study within this smaller cohort. Having full cohort
information enables comparison between wpl() and
corresponding analyses carried out on the full cohort. The synthetic
data generation was performed in two steps. In the first step a nested
case-control sampling was carried out by sampling one control per case
matched sex and BMI \(\pm\) 2. This was
done in a for-loop over event times sampling one subject at random from
those still alive at that time in addition to fitting the matching
criterium. When the controls are not additionally matched or matched
only on a categorical variable the ccwc-function from the
Epi
package (Carstensen et al. 2016) offers a
simple way to create a nested case-control design from a cohort.
However, if some of the matching variables are continuous it still has
to be done as described above. The second step is necessary only for the
weighted Cox-regression and it involves removing duplicate subjects, due
to that some controls are sampled multiple times and some controls later
become cases, from the sampled data. For those controls who later become
cases it is important to keep the “case-entry”, while for those controls
sampled multiple times, the entries will be identical and one of them is
kept at random.
Four types of deaths are recorded in the data: (1) cancer, (2)
cardiovascular disease, (3) alcohol abuse, liver disease, accidents and
violence, and (4) other medical causes. For simplicity we only use
cardiovascular disease (\(n=236\)) and
alcohol abuse, liver disease, accidents and violence, henceforth
referred to as ALAV (\(n=60\)). The
data set is included in multipleNCC and can be loaded by
data("CVD_Accidents").
When we analyze this data set with IPW, we will use all sampled controls for both endpoints, hence supplement the controls for ALAV deaths with the cardiovascular disease controls and cases. And vice versa, supplement controls for cardiovascular disease cases with cases and controls for ALAV cases. For ALAV deaths, the number of controls increases substantially when including the controls for cardiovascular deaths.
The matching variables are BMI (\(M_1\) = bmi) and sex (\(M_2\) = sex). The matching
variables are adjusted for in the model to remove confounding due to
breaking the matching. We included BMI as continuous variable and sex is
a categorical variable. Smoking (\(X\)
= smoking3gr) is our explanatory variable and is
categorized into never smoked, former smoker and current smoker.
We consider the two models \[\begin{aligned} h_{1i}\left(t_i|X,M_1,M_2\right) &= h_{10}\left(t_i\right)\exp\left(\beta_1' X_i + \gamma_{11} M_{1i} + \gamma_{12} M_{2i}\right)\\ h_{2i}\left(t_i|X,M_1,M_2\right) &= h_{20}\left(t_i\right)\exp\left(\beta_2' X_i + \gamma_{21} M_{1i} + \gamma_{22} M_{2i}\right). \end{aligned}\] Here \(h_1\left(.\right)\) and \(h_2\left(.\right)\) are the hazard for death from CVD and ALAV, respectively. To fit the weighted partial likelihoods corresponding to those models, the command below can be used
R> fit <- wpl(Surv(agestart, agestop, dead24) ~ factor(X) + bmi + sex,
+ data = CVD_Accidents, m = 1, samplestat, weight.method = "glm",
+ match.var = cbind(CVD_Accidents$BMI, CVD_Accidents$sex),
+ match.int = c(-2, 2, 0, 0))The status argument to the Surv function
(in this example dead24) can in practice contain anything
since the information regarding who are cases and controls, and also
which type of endpoint the cases experienced are provided through
samplestat. The match.int argument will in
this situation be a vector of two elements,
match.int = c(-2, 2, 0, 0), since the controls are matched
to the cases on BMI \(\pm\)2,
and sex.
The summary and print commands can be used
to display the results of the analysis. The output resembles output from
traditional Cox-regression (see next page), and the results for the
different endpoints are printed below each other. Both the naive and
robust/estimated standard errors are printed, although only the
robust/estimated should be reported.
The $-operator is usually used when specific parts of an
object (i.e., fit$coefficients) are to be extracted. When
there is more than one endpoint this will only give you the
corresponding element for the first endpoint, i.e.,
fit$coefficients will give you (0.47107, 1.34245, 0.08051,
\(-\)1.22307), thus only the
estimates for the cardiovascular disease endpoint. To extract any
element from the object, [[]] can be used, i.e.,
fit[[23]] will give you the estimates from the second
analysis of death from ALAV. For a full list of elements in the wpl
object, the names() function can be used. Each element will
occur the same number of times as there are different endpoints, and the
first occurrence corresponds to the first endpoint, the second
occurrence to the second endpoint, etc.
R> summary(fit)
Endpoint 1 :
Call:
wpl.formula(formula = Surv(agestart, agestop, dead24) ~ factor(smoking3gr) +
bmi + factor(sex), data = CVD_Accidents,
samplestat = CVD_Accidents$samplestat,
match.var = cbind(CVD_Accidents$bmi, CVD_Accidents$sex),
match.int = c(-2, 2, 0, 0), weight.method = "glm")
n= 566, number of events= 236
coef exp(coef) se(coef) robust se z Pr(>|z|)
factor(smoking3gr)2 0.47107 1.60171 0.22099 0.26057 1.808 0.07062 .
factor(smoking3gr)3 1.34245 3.82842 0.19400 0.23424 5.731 9.98e-09 ***
bmi 0.08051 1.08384 0.01799 0.02562 3.143 0.00167 **
factor(sex)2 -1.22307 0.29433 0.15980 0.22475 -5.442 5.27e-08 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
factor(smoking3gr)2 1.6017 0.6243 0.9611 2.6692
factor(smoking3gr)3 3.8284 0.2612 2.4190 6.0591
bmi 1.0838 0.9226 1.0308 1.1396
factor(sex)2 0.2943 3.3976 0.1895 0.4572
Endpoint 2 :
Call:
coxph(formula = Surv(left.time.ncc, survtime.ncc, status.ncc ==
i) ~ x + cluster(ind.no.ncc), weights = 1/p)
n= 566, number of events= 60
coef exp(coef) se(coef) robust se z Pr(>|z|)
factor(smoking3gr)2 -0.61629 0.53994 0.45484 0.48347 -1.275 0.202413
factor(smoking3gr)3 0.92343 2.51792 0.32202 0.34402 2.684 0.007270 **
bmi 0.08383 1.08744 0.03813 0.04587 1.828 0.067619 .
factor(sex)2 -1.42549 0.24039 0.32702 0.36888 -3.864 0.000111 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
factor(smoking3gr)2 0.5399 1.8520 0.2093 1.3928
factor(smoking3gr)3 2.5179 0.3972 1.2829 4.9417
bmi 1.0874 0.9196 0.9939 1.1897
factor(sex)2 0.2404 4.1599 0.1167 0.4954The sampling probabilities can be directly estimated using one of the
functions: KMprob, GAMprob,
GLMprob and Chenprob, for instance
R> glmp <- GLMprob(CVD_Accidents$agestop, CVD_Accidents$samplestat,
+ left.time = CVD_Accidents$agestart, match.var = cbind(CVD_Accidents$sex,
+ CVD_Accidents$bmi), match.int = c(0, 0, -2, 2))
R> kmp <- KMprob(CVD_Accidents$agestop, CVD_Accidents$samplestat, 1,
+ left.time = CVD_Accidents$agestart, match.var = cbind(CVD_Accidents$sex,
+ CVD_Accidents$bmi), match.int = c(0, 0, -2, 2))It is of interest to examine the weights for the sampled controls only, as the cases have weight 1 and the non-sampled subjects will not affect estimates. We can for instance inspect summary statistics
R> summary(1/glmp[CVD_Accidents$samplestat == 1])
Min. 1st Qu. Median Mean 3rd Qu. Max.
4.575 6.545 9.201 13.240 15.080 73.180
R> summary(1/kmp[CVD_Accidents$samplestat == 1])
Min. 1st Qu. Median Mean 3rd Qu. Max.
2.481 6.839 8.772 13.150 15.020 95.810The two types of weights correspond well, although KM-weights have a somewhat heavier tail. It is worth noting that the subjects with the largest weights are those with the shortest follow-up time and therefore those subjects do not have a too large impact on the analysis even though their weight is large, since they are included in few risk sets. GAM-weights were similar to GLM-weights, but with a somewhat heavier tail (not shown) and Chen-weights are not applicable here since they are not implemented for additional matching.
coxph and
wplTable 2 displays a comparison between the full
cohort analysis, the traditional estimator for nested case-control data
and the IPW-estimator using wpl() with GLM- and KM-weights.
Being a smoker significantly increases the risk of death from
cardiovascular disease. Being a former smoker also increases the risk of
death from cardiovascular disease, although only significantly in the
cohort analysis. Being a former smoker has a non-significant protective
effect on death from ALAV, while being a smoker significantly increases
the risk of dying from ALAV.
The hazard ratios estimated with the traditional estimator and with
wpl() are fairly similar to the cohort estimates taking
into account the size of the standard errors. The most pronounced
difference is between the cohort hazard rate and the hazard rate for the
traditional estimator for being a smoker on death from ALAV, 2.05
vs. 2.90. However, considering the size of the standard errors this is
not a large difference.
For the cardiovascular disease endpoint, the standard errors of the IPW analyses are somewhat smaller than the standard errors of the traditional estimator, resulting in a little bit higher efficiency for the IPW-estimators. For the ALAV endpoint, the standard errors are substantial lower and a large efficiency gain is obtained with the IPW-estimators. The reason for this is that the IPW-estimators make use of a number of extra controls as all cases and controls from the cardiovascular disease endpoint are used as additional controls. On average the number of controls per case increase from 1 to more than 8. We have chosen to include cardiovascular disease cases as additional controls for the cases who died from ALAV, and also the cases who died from ALAV as additional controls for the cases who died from cardiovascular disease. However, the cases who are used as additional controls will contribute little to the analysis since they are non-cases (for the particular endpoint) with weight equal to 1. We have reported robust standard errors for KM-weights in Table 2, however the estimated standard errors are very similar, for CVD endpoint 0.27 and 0.24 and for ALAV 0.48 and 0.35 for former and current smoker respectively.
| CVD | |||||||||
| Former smoker | Current smoker | ||||||||
| 2-4 | HR (95 % CI) | SE(\(\beta\)) | Eff. | HR (95 % CI) | SE(\(\beta\)) | Eff. | |||
| 2-4 Cohort | 1.72 (1.11–2.66) | 0.22 | 1.00 | 3.25 (2.21–4.79) | 0.20 | 1.00 | |||
| Trad. | 1.66 (0.94–2.93) | 0.29 | 0.59 | 3.52 (2.09–5.94) | 0.27 | 0.55 | |||
| IPW-GLM | 1.60 (0.96–2.67) | 0.26 | 0.73 | 3.83 (2.42–6.06) | 0.23 | 0.72 | |||
| IPW-KM | 1.65 (0.98–2.78) | 0.26 | 0.66 | 3.97 (2.48–6.35) | 0.24 | 0.69 | |||
| ALAV | |||||||||
| Former smoker | Current smoker | ||||||||
| 2-4 | HR (95 % CI) | SE(\(\beta\)) | Eff. | HR (95 % CI) | SE(\(\beta\)) | Eff. | |||
| 2-4 Cohort | 0.56 (0.23–1.38) | 0.46 | 1.00 | 2.05 (1.07–3.90) | 0.33 | 1.00 | |||
| Trad. | 0.48 (0.13–1.69) | 0.65 | 0.50 | 2.90 (1.11–7.61) | 0.49 | 0.45 | |||
| IPW-GLM | 0.54 (0.21–1.39) | 0.48 | 0.90 | 2.52 (1.28–4.94) | 0.34 | 0.92 | |||
| IPW-KM | 0.55 (0.21–1.40) | 0.49 | 0.90 | 2.56 (1.28–5.02) | 0.35 | 0.89 |
We have demonstrated the multipleNCC package in R which
allows for breaking the matching and reusing controls in NCC designs.
The main function wpl() estimates sampling probabilities
and perform weighted Cox-regressions. It handles right-censored,
left-truncated and additionally matched data, and varying number of
sampled controls. We have also explained how the variance can be
estimated without additional matching (KM- and Chen-weights) and with
additional matching (KM-weights).
The package is particularly useful in situations with multiple
outcomes. It has a competing risks perspective, in the sense that with
more than one type of endpoint, each endpoint is estimated separately
and all controls and cases of other types are used as additional
controls. In many situations the competing risks framework may not be
suitable, although reusing controls can still be of interest. The
solution can then be to estimate the sampling probabilities for all
cohort members, using one of the four functions KMprob,
GAMprob, GLMprob or Chenprob, and
carry out weighted analysis that fit the situation at hand.
The gam() function from the mgcv package is
used for estimation of the GAM-weights. We could alternatively have used
the gam() function from the gam package
(T. Hastie 2015) or even a different form
of smoothing. It has however become evident that the exact value of the
weights are not too important as long as they are fairly reasonable,
thus the choice of smoothing does probably not affect final hazard
ratios and standard errors. For the same reason there are usually only
minor differences with regards to final hazard ratios and standard
errors between the four weight estimators discussed in Section 3 (N. C. Støer
and Samuelsen 2012, 2013).
Sometimes the cases and controls are matched closer together than is
strictly necessary. Very close matching can lead to small KM-weights
since the probability of being sampled will be large for the controls
that were sampled (and very small for most of the non-sampled subjects)
and this could lead to biased estimates (N. C.
Støer and Samuelsen 2013). A solution could be to increase the
length of match.int. For example in the data example above,
we could replace match.int = c(-2, 2) with
match.int = c(-4, 4). Widening the matching interval could
introduce bias, thus it is important to carry this out with caution. The
equivalence of this for category matching is to reduce the number of
levels of the matching variable by some sort of grouping.