Pooling controls from nested case–control studies with the proportional risks model

Abstract

1 Introduction

Cohort studies often include multiple failure types, and it is common to consider one failure type as the main event of interest and treat others as censoring, where such censoring is considered as competing risks. For example, when the main interest is to understand risk factors associated with colorectal cancer incidence, diagnoses of other cancers are considered as competing risks, since the occurrence and treatment of such cancers may impact the risk for subsequent cancers. The analysis of time to first cancer is natural in the competing risks framework, as adopted by studies on the association between vitamin D and cancer incidence in the Prostate, Lung, Colorectal, and Ovarian Cancer Screening Trial (PLCO) (Mondul et al. 2012; Weinstein et al. 2015).

The impact of covariates on a cause or type of failure can be estimated by a Cox (1972) proportional hazards regression model on the cause-specific hazard function (Kalbfleisch and Prentice 2002). The usual approach fits regression models separately to each failure type. We focus on an alternative proportional risks model for cause-specific hazards (Holt 1978; Lunn and McNeil 1995), which assumes that the baseline hazard functions are proportional to each other. This encompassing model permits a formal comparison of covariate effects across causes, e.g. in comparing the effects of risk factors on different subtypes of a disease (rgp120 HIV Vaccine Study Group 2005; Song et al. 2016). Moreover, there is potentially an efficiency gain by simultaneously fitting a single proportional hazards model for all causes. With prospective cohort data, the proportional risks model can be fitted easily using standard software by double coding the data and using the cause as a model covariate (Lunn and McNeil 1995).

If information on all covariates are available in a cohort, regression coefficients of the proportional risks model can be estimated by maximizing a partial likelihood. In contrast to the standard approach that conditions on a failure from a particular cause, the Lunn and McNeil approach conditions on a failure from any cause, which leads to a single partial likelihood for all failure types. However, resources are often limited such that covariates are available only for a subset of cohort participants and partial likelihood analyses are not applicable. For example, it may be too expensive to genotype, measure plasma biomarkers, or extensively interview all cohort participants. The nested case–control design (NCC), one of the most popular cohort sampling designs, samples a small number of controls for each case from subjects who are uncensored and at risk at the case’s failure time (Langholz and Thomas 1990). Covariates are only collected on the cases and the selected controls, and a matched set of a case and its control(s) is called a sampled risk set. For a standard analysis, the regression coefficients can be estimated with NCC data by maximizing the product of conditional probabilities that condition on the cause of failure in each sampled risk set (Lubin 1985).

In this paper, we investigate a conditional likelihood based on the proportional risks model, which has not been studied in the literature. This conditional likelihood can pool controls from multiple NCC studies and potentially give more efficient estimates. Current methods for pooling or reusing NCC data are primarily based on inverse probability weighting (IPW) or maximum likelihood estimation (MLE) (Saarela et al. 2008; Salim et al. 2012; Støer and Samuelsen 2013), which may be more efficient than conditional likelihoods because of the use of all-cause controls or by virtue of distributional assumptions. However, they may be more cumbersome to implement, not as robust to batch effects, or require larger cohort sizes (see Section 2.3). Thus, the proposed new conditional likelihood provides a novel and robust alternative to reusing existing NCC data.

In Section 2, we formally derive the new conditional likelihood under the proportional risks model with NCC data. Scenarios are considered where NCC data for different failure types are either obtained simultaneously in a single study, or pooled from different NCC studies. Unlike the standard analyses, matching variables and different follow-up lengths may influence the conditional likelihood and bias may result if ignored. We propose simple adjustments which give unbiased estimates. As with full cohort data, double coding may be used with standard software. We also discuss settings where the new conditional likelihood may be preferred over other methods for pooling NCC data. Section 3 reports simulation studies showing that limited efficiency gains are observed with prospective cohort data, but large gains are evident with NCC data. In Section 4, a comprehensive analysis of colorectal and prostate cancer incidence from the PLCO trial demonstrates the practical utility of the Lunn and McNeil approach to combining multiple NCC studies. A discussion concludes the paper in Section 5.

2 Materials and Methods

2.1 Proportional risks model for cause-specific hazards in full cohort studies

The estimates $β^$ from $L1(β)$ and $L2(β)$ are consistent for $β$⁠, and their variance can be estimated by the inverse of the observed information matrices.

Simultaneous estimation of $βk$ ’s using $L2(β)$ under model (2.2) may be more efficient than using $L1(β)$ under model (2.1). However, if the baseline hazard functions do not satisfy the proportionality assumption, estimates from $L2(β)$ may be biased. This assumption can be evaluated by checking if the estimated cumulative baseline hazards $Λ^k0(t)$ under model (2.1) are proportional. $Λ^k0(t)=∑τj≤tdΛ^k0(τj)=∑τj≤tI(ϵj=k)/∑l∈R(τj) exp (ZlTβ^k)$⁠, where $I(·)$ is an indicator function, and $β^k$ is obtained from fitting model (2.1) (Kalbfleisch and Prentice 2002). $Λ^k0(t)$ is the Breslow estimator for cause-k cumulative baseline hazard where failures from other causes are treated as censoring (Breslow 1972). Alternatively, $λk0(t)$’s can be estimated by smoothing $dΛ^k0(τj)$’s, the discrete increments of cumulative hazards (Gray 1990). $Λ^k0(t)$ or smoothed $dΛ^k0(τj)$ may be depicted graphically for a visual model check.

When the baseline hazards are not proportional, Lunn and McNeil (1995) recommend using model (2.1), or breaking down the follow-up time into intervals in which the proportionality assumption is tenable, and fit model (2.2) in a piecewise fashion. Alternatively, the proportionality factors can be time-dependent, e.g. piecewise constants or polynomials, to approximate the true ratios of baseline hazards. The plotted $Λ^k0(t)$ and the smoothed $dΛ^k0(τj)$ should inform reasonable functions of time for the proportionality factors. Such time-dependent proportionality factors can reduce bias and improve the efficiency of model (2.2) estimates, which will be shown in our simulation study. A formal test on the proportionality assumption can be achieved (Belot et al. 2010), for instance, by testing whether all piecewise constants are equal, or whether the linear and higher order terms in the polynomials are all equal to zero.

To assess whether $γk$’s are specified appropriately in model (2.2), one can check graphically if the estimated baseline hazards based on model (2.1) are close to those based on model (2.2), which is $Λ^k0(t)=∑τj≤t exp {I(ϵj=k)γ^k(τj)}dΛ^0(τj)$ where $dΛ^0(τj)=1/[∑l∈R(τj)∑k=1K exp {γ^k(τj)+ZlTβ^k}]$ is the increment of the joint cumulative baseline at $τj$ (Tai et al. 2001), and $β^k$’s and $γ^k(τj)$’s are from model (2.2), with $γ^1(t)=0$ for identifiability.

Model (2.1) may be fitted using standard software for the proportional hazards model, treating failure from other causes as independent censoring. Model (2.2) can be fitted using standard software following the non-standard data augmentation method developed by Lunn and McNeil (1995), outlined below. Each row of data is replicated K times, and a different failure indicator variable is created to take value 1 at the k-th row if the subject fails from cause k, and 0 otherwise. Models for all causes are simultaneously fitted using $L2(β)$ based on the augmented dataset. As with model (2.2), standard output based on data augmentation may be used for inference.

2.2 Proportional risks model for cause-specific hazards in NCC studies

2.2.1 NCC design for multiple outcomes

For a single outcome, the NCC design samples controls for a case from its risk set. Complete covariates for the intended analyses are only collected for the cases and their controls. To study multiple competing outcomes, one may prospectively design a single study that samples cases and their matched controls (Borgan 2002), where a case is defined as a failure from any cause. We present the conditional likelihoods for model (2.1) and (2.2) using such prospective NCC samples from the same cohort without matching in Section 2.2.2. When combining multiple NCC studies for different failure types from the same cohort, the studies may have different matching variables, start and end of follow-up, and inclusion-exclusion criteria. We show in Sections 2.2.3 to 2.2.6 the modifications needed to address these issues.

2.2.2 Models and likelihoods

(Appendix 1 of the Supplementary Materials), which resembles $L2(β)$ with $R(τj)$ replaced by $R˜(τj)$⁠. Note that $L˜1(β)$ and $L˜2(β)$ differ in what is conditioned on, similar to their full cohort counterparts. When model (2) holds, the denominator of $L˜2(β)$ implies that all-cause controls are used to estimate all $βk$ ’s, thus $L˜2(β)$ is expected to yield more efficient estimates than $L˜1(β)$⁠, where controls for cause k are only used to estimate $βk$⁠. In contrast, both $L1(β)$ and $L2(β)$ use all at-risk subjects as controls for estimation of $β$⁠, so the efficiency gain of the latter over the former is expected to be less, but the latter still has the advantage of allowing comparison of covariate effects across causes.

Inference on $β$ using NCC data may be based on estimates and variance from $L˜1(β)$ and $L˜2(β)$⁠, similar to using standard output from $L1(β)$ and $L2(β)$ for inference in a full cohort analysis. Note that contrasts of covariate effects across causes are not interpretable with $L1(β)$ and $L˜1(β)$⁠, because the baseline hazard functions are free to take any shape under model (2.1) and not necessarily proportional as in model (2.2).

2.2.3 Additional matching variables

In practice, there could be additional matching variables other than at-risk status such as age or gender. Regardless of what variables are matched on in each NCC study, they cancel out in $L˜1(β)$ because of conditioning. However, if these variables are associated with risks of one or more failure types, their association should be accounted for in $L˜2(β)$ because they do not cancel out.

This conditional likelihood corresponds to including the matching variables $V$ in the model of the second failure type, thus it can be maximized with standard software. The individual effects of the matching variables on the risks of the two failure types, $b22$ and $b12$⁠, are not identifiable. However, their difference, $b22−b12$⁠, is identifiable and $V$ should be included in the model, as $b22−$$b12$ only vanishes from (2.3) when $b12$ and $b22$ are equal, which is generally unknown a priori.

In this setup, the two studies are assumed to match on variables with the same forms. If the two studies match on a variable with different forms, we suggest choosing the form with wider groups, or a continuous form if available. For example, if one study matches on 5-yr age groups and the other matches on 10-yr age groups, formula (2.3) may be applied by including 10-yr age groups or continuous age in the model. In reality, different NCC studies are often matched on different variables even if they arise from the same cohort. For example, the PLCO trial conducted a series of NCC studies on the association between vitamin D and risk of different types of cancers, all of which were matched on sex and race, but the lung cancer study was the only study that matched on smoking history (Muller et al. 2018), and all except the colorectal adenoma study matched on age or year of birth (Peters et al. 2004). In this scenario, the same principle of including matching variables in model (2.2) still applies, but the correct likelihood is no longer (2.3).

Unlike (2.3), (2.4) does not take the form of a usual Cox partial likelihood. Therefore, standard software for Cox models does not readily maximize (2.4). To maximize (2.4), one may still double code the data, and find the vector $(a^,b^11,b^12,b^13,b^21,b^21,b^23)T$ that maximizes (2.4) through algorithms such as the Newton-Raphson method.

Construction of (2.3) and (2.4) requires all matching variables to be observed in all contributing NCC studies, which is usually the case for NCC studies from a single cohort, like the PLCO vitamin D NCC studies. Although (2.3) and (2.4) are based on two competing failure types, extension to three or more failure types is straightforward.

2.2.4 Different follow-up intervals

2.2.5 Different inclusion–exclusion criteria

When studies have different inclusion-exclusion criteria, the combined analysis using model (2.2) can be confined to subjects who meet all the criteria, but this may lead to a significant drop in sample size and difficulty in generalizing the results to a broader population. For example, when combining the colorectal and the prostate cancer studies, women from the colorectal cancer study are not at risk of prostate cancer, but confining the analysis to men means loss of efficiency and generalizability. One way to retain subjects in the combined analysis is to condition on which sets of inclusion criteria each sampled risk set satisfies, rather than to require all sampled risk sets to meet the superset of all inclusion criteria.

The augmented dataset for $L˜IC(β)$ and $L˜FU,IC(β)$ can be set up similar to Section 2.2.4 by replicating data K times, creating cause-specific failure indicators, removing sampled risk sets with different eligibility to the K studies, and removing rows where $∏l∈R˜(τj)r˜lk=0$ or $I(τj∈[lk,rk])∏l∈R˜(τj)l˜lk=0$⁠. Standard output based on the augmented dataset can be used for inference on $β$⁠.

2.2.6 The proportionality assumption

Similar to a full cohort analysis, the adequacy of the time-dependent proportionality factors can be evaluated graphically, using the cause-k cumulative baseline hazard estimates$Λ^k0(t)=∑τj≤t exp {I(ϵj=k)γ^k(τj)}dΛ^0(τj)$ where $dΛ^0(τj)=1/∑l∈R˜(τj)∑k=1K exp {γ^k(τj)+ZlTβ^k}wk(τj)$⁠, $β^k$’s and $γ^k(τj)$’s are from model (2.2), and $γ^1(t)=0$⁠.

2.3 Comparison to current methods for pooling NCC data

As seen in Section 2.2.2, both $L˜1(β)$ and $L˜2(β)$ can be used to pool NCC samples for different outcomes into a competing risks analysis. $L˜2(β)$ and its variations are more efficient because they pool subjects within each sampled risk set to estimate the entire $β$ vector. We review MLE and IPW methods for pooling NCC samples into competing risks analyses, and compare them to $L˜1(β)$ and $L˜2(β)$ in terms of efficiency, robustness, and computation.

Current MLE and IPW methods assume model (2.1) and fit regression models to each failure type separately. The MLE approach requires specifications of baselines and the distribution of covariates only observed in the NCC sample given covariates observed in all cohort members. When correctly specified, MLE is more efficient than IPW, and the efficiency advantage is mainly attributed to utilization of full cohort information (Saarela et al. 2008). However, MLE is computationally intensive, sensitive to starting values and misspecification of the conditional distribution of the partially observed covariates (Støer and Samuelsen 2012). This approach is not commonly used in practice, and no off-the-shelf software is available.

The IPW approach weights a subject’s log-likelihood or log-partial likelihood contribution by the inverse of his probability to be included in at least one of the NCC studies (Saarela et al. 2008; Salim et al. 2012). Weighting of log-partial likelihood is also known as weighted partial likelihood (WPL), and it does not require specification of the baselines. WPL is more efficient than $L˜1(β)$⁠, as sampled risk sets from all studies are pooled together for estimation of each $βk$⁠. WPL is reported to be more robust to misspecification of the Cox models than MLE, but the variance estimation for WPL is less straightforward or lacks theoretical justification depending on the type of weights used (Støer and Samuelsen 2012). In the presence of additional matching, variance estimation may be computationally intensive and estimates may be biased if matching variables are omitted from the Cox models and in the estimation of inclusion probabilities. WPL may also break down under close matching or strong batch effects (Støer and Samuelsen 2013). The R package multipleNCC (Støer and Samuelsen 2016) may be used to conduct a WPL analysis.

Since WPL is more utilized and robust among current methods, we focus on comparing WPL with $L˜1(β)$ and $L˜2(β)$ for data generated under model (2.2). Due to page limit, we only state the main findings here and leave the details in Appendix 3 of the Supplementary Materials. Our simulations show that WPL is more efficient than both conditional likelihoods in simple settings, e.g. when there is no additional matching or when the cohort size is large (Appendix 3 Study A). The better efficiency of WPL over $L˜2(β)$ is likely because WPL pools subjects across sampled risk sets but $L˜2(β)$ does not. With close matching or a small cohort size, weights may be too small to construct the pseudopopulation for inference, which leads to biased WPL estimates. In contrast, $L˜1(β)$ and $L˜2(β)$ do not depend on weights and are thus more robust insuch situations.

WPL estimates may also be biased when covariates are measured in batches and are more similar within than between batches (Appendix 3 Study B). Unlike $L˜1(β)$⁠, batch effects do not cancel out in WPL since WPL is unconditional. On the other hand, batch effects are absorbed by the proportionality factor in $L˜2(β)$ as long as subjects in the same sampled risk set are placed in the same batch. The proportionality factor estimate may be undercovered as a result, but estimates for $β$ are valid. Thus, $L˜2(β)$ and $L˜1(β)$ may be preferred over WPL if batch effectsare suspected.

Another scenario where WPL estimates may be biased is when matching variables affect risks of different outcomes in the same way, but cannot be adequately modeled as Cox model covariates. For example, consider cause-specific hazards $λk(t|U,V)=α(V)tα(V)−1 exp (ak+$$bk1TU)$⁠, where $α(V)$ is a positive function of matching variables $V$⁠, and $U$ is the exposure vector. Specifying $V$ as covariates may not be adequate to model their relationships with the outcomes, resulting in biased WPL and full cohort estimates. Conditional likelihoods may be more robust in such settings (Appendix 3 Study C).

3 Simulation studies

3.1 Data generation

We consider two competing causes and generate failure times using cause-specific hazards $λk(t|Z)=exp (γk+βk1Z1+βk2Z2)αktαk−1$⁠, in which $λk0(t)=exp (γk)αktαk−1,k=1,2$⁠. The covariates $Z1$ and $Z2$ are simulated from the standard normal distribution. We fix $α1=1$ and let $α2$ take different positive values so that $α2=1$ represents proportional baseline hazards of the two causes, while $α2$ away from 1 represents deviation from that assumption. The entry time for each subject follows Uniform$(0,2)$⁠. A random censoring time is independently generated from an exponential distribution such that the probability of dropping out or lost to follow-up is 0.1 at the administrative censoring time $(t=10)$⁠. The on-study event time is the minimum of the failure time and the censoring time (the minimum of $(10−$entry time) and the random censoring time).

3.2 Scenarios

We simulate data under five different scenarios and compare the bias and efficiency of analyses using models (2.1) and (2.2). When a scenario allows both full cohort and NCC analyses, we show how efficiency gain of model (2.2) over model (2.1) differ between the full cohort and the NCC analyses. Time-dependent proportionality factors, bias corrections, and inclusion of matching variables for model (2.2) are applied when appropriate.

For all scenarios, $β$ are chosen such that the hazard ratios associated with 1 unit increase in $Z1$ is 2 for cause 1 and 1.5 for cause 2, and the hazard ratios associated with 1 unit increase in $Z2$ is 2.5 for cause 1 and 3 for cause 2. We set $γ1$ to $log(−log(0.9)/10)$ and $γ2$ to $log(−log(1−p2)/10α2$⁠) so the survival probabilities for $z=(0,0)T$ at the administrative censoring time of 10 yrs is 0.9 for cause 1 and $1−p2$ for cause 2. Unless specified otherwise, $α2=1,p2=0.05$⁠, and 500 cohorts of 10000 subjects are simulated for each true parameter value in each scenario. A full cohort analysis includes all subjects in a simulated cohort, while an NCC analysis includes only the cases and their matched controls. Each case is matched to one control on at-risk status, and additionally on age in Scenario 3.3.3. In all simulations, relative efficiency (RE) is the ratio of mean square errors (MSEs) of the model (2.2) estimates to the model (2.1) estimates.

3.2.1 Scenario 1: proportional risks and different incidence rates

Data are simulated under $p2=0.02,0.1$⁠, and 0.2 to explore the performance of the two models under proportional risks and different baseline incidence rates. The average sizes of the pooled NCC samples are $29%,41%$⁠, and $53%$ that of the full cohort for $p2=0.02,0.1$⁠, and 0.2, respectively.

3.2.2 Scenario 2: non-proportional risks

Data are simulated under $α2=0.5$ and 5 so the baseline hazards of the two causes are not proportional. Supplementary Figure S1 of the Supplementary Materials shows the two theoretical baseline hazards, where the baseline hazard of cause 1 is constant due to $α1=1$⁠, and the baseline hazard of cause 2 is monotone decreasing or monotone increasing under $α2=0.5$ or 5, respectively. The pooled NCC samples are on average $35%$ and $31%$ the full cohort size for $α2=0.5$ and 5, respectively. The proportionality factor is specified in three ways: a constant (model (2.2)), two piecewise constants, or a linear function of time. The boundaries for the piecewise constants are $t=2.5$ yrs for $α2=0.5$ and $t=7.5$ yrs for $α2=5$⁠, so that the cumulative baseline hazards are roughly proportional before and after.

3.2.3 Scenario 3: matching on continuous age

We generate age from Uniform$(18,65)$⁠, and let the two cause-specific hazards additionally depend on normalized age $Z3$⁠. Specifically, $λk(t|Z)=exp (γk+βk1Z1+βk2Z2+βk3Z3)$⁠, where $β13=log 1.6$ and $β23=log 1.1$⁠. Each case is matched to one control at risk and with an age difference under 2 yrs. On average, the pooled NCC samples are $35%$ of the full cohort size. We focus on the NCC analysis in this scenario to show how omitting matching variables affects the estimates.

3.2.4 Scenario 4: different follow-up intervals

The NCC study for cause 1 is administratively censored at year 7, earlier than year 10 of the full cohort data and the NCC study for cause 2. On average, the pooled NCC sample includes $29%$ of the cohort subjects and $78%$ of failures from cause 1.

3.2.5 Scenario 5: different inclusion-exclusion criteria

For each cohort member, we generate two independent variables. $W1$ is drawn from the standard normal distribution, and subjects are eligible to the NCC study for cause 1 if $W1≤0.68$⁠. $W2$ is drawn from Bernoulli$(0.8)$⁠, and subjects with $W2=1$ are eligible to the NCC study for cause 2. For each cause, eligible cases are matched to eligible controls on at-risk status. Only an NCC analysis is performed. An average of $19%$ cohort subjects are included in the pooled NCC sample.

3.3 Result summary

3.3.1 Scenario 1

Under proportional risks and different relative incidence rates, model (2.2) yields unbiased estimates close to those from model (2.1) for both the full cohort (Table 1) and the NCC analyses (Table 2). In the full cohort analysis, the Monte Carlo standard error (SE) and MSE are similar for the two models so the RE of model (2.2) to model (2.1) is around 1. In the NCC analysis, model (2.2) estimates have much smaller SEs and MSEs than model (2.1) estimates, with RE ranging from 1.13 to 3.24. With model (2.2), it is observed that the coefficients of covariates associated with the rarer outcome gain more efficiency than those associated with the more common outcome. This is because $βk$ of the rarer failure type gains more controls from the more common failure type. In summary, the efficiency gain of model (2.2) over model (2.1) is minimal in the full cohort setting as both $L2(β)$ and $L1(β)$ use all cohort members. However, it is substantial in the NCC setting because $L˜2(β)$ uses all cases and controls for estimation of each $βk$⁠, whereas $L˜1(β)$ only uses cases and controls from study k for estimation of $βk$⁠.