Biases Introduced by Choosing Controls to Match Risk Factors of Cases in Biomarker Research

Data for bivariate binormal scenarios

Let X denote the baseline risk factor variables, Y the biomarker and D the binary outcome with D=1 for a case and D=0 for a control. For a population of N=1,000,000 subjects we generated D with probability D=1 being 10%. For controls, (X,Y) were generated as bivariate normal with means (0,0), standard deviations (1,1) and correlation ρ = 0.1, 0.5 or 0.9. For cases, (X,Y) were bivariate normal with the same standard deviations and correlation as for controls but with means (μ_X, μ_Y), where μ_X>0 and μ_Y>0. In this scenario, the magnitude of μ_X determines the performance of X alone and μ_Y determines the performance of Y alone. For example, the area under the ROC curve for X is ${\mathrm{AUC}}_X=\Phi \left({\mu }_X{∕}^{\sqrt[]{2}}\right)$ where Φ is the standard normal cumulative distribution function, and the true positive rate corresponding to the threshold that allows a false positive rate f, denoted by ROC_X(f), is ROC_X(f)= Φ(μ_X+ Φ⁻¹(f)). We fixed ROC_X(0.2) (and ROC_Y(0.2) ) and used the corresponding values of μ_X (and μ_Y) to generate data. For example, in Figure 1 we wanted a setting where ROC(0.2)=0.64 for X alone and for Y alone so we used values μ_X = μ_Y=1.19 because the associated value of ROC(0.2)= Φ(μ+ Φ⁻¹(0.2)) is 0.64.We selected 20,000 cases at random. For the ‘unmatched controls’ analyses we selected 20,000 controls at random from the population. For the ‘matched controls’ analyses we generated a matched control for each of the 20,000 cases. Specifically, with X_i the clinical risk factor predictor for the i^th case, we selected a random control with the same value of X= X_i as the case. Technically we implemented this by generating the matched control’s biomarker, Y_i, from the conditional distribution of Y given X=X_i and D=0, i.e. as a normal variable with mean = ρ × X_i and standard deviation = √(1-ρ²) where ρ is the correlation between X and Y.

To demonstrate that the magnitude of bias can affect conclusions drawn in studies with moderate sample sizes, we generated datasets with 100 cases and 100 controls rather than 20,000 cases and 20,000 controls.

We implemented two sets of simulations to investigate if matching can distort comparisons between markers. In the first scenario, where we compared markers Y_A and Y_B used on their own, (X,Y_A,Y_B) is multivariate normal in controls with means(0,0,0) and standard deviations (1,1,1) and correlation(X,Y_A)=0.7, correlation(X,Y_B)=0.1 and correlation (Y_A,Y_B)=0. In cases the means are (1.19, 1.19, 0.8) and the standard deviations and correlations are the same as in controls. In the second scenario, where we compared the incremental values of markers Y_A and Y_B over use of X alone, (X,Y_A,Y_B) is standard multivariate normal in controls. In cases (X,Y_A,Y_B) is multivariate normal with means(1,1,1), standard deviations (1,1,1) and correlation(X,Y_A)=0.1, correlation(X,Y_B)=0.7 and correlation(Y_A,Y_B)=0.

Data for predicting critical illness

We investigated the impact of matching controls to cases in the context of a study to identify patients in need of critical care services who undergo care in the prehospital setting (e.g. ambulance). We used a validation dataset of 57,647 patients from a population-based cohort study of an emergency medical services (EMS) system in King County, Washington. Methods and results of this study have been published [11]. In brief, the study developed a prediction model for critical illness in non-trauma non-cardiac-arrest adult patients transported to a hospital by EMS from 2002 through 2006. Critical illness was defined as severe sepsis, delivery of mechanical ventilation, or death during hospitalization. Predictors of risk in the model included older age, male gender, lower systolic blood pressure, abnormal respiratory rate, lower Glasgow Coma Scale, lower pulse oximetry and residence in a nursing home. Studies investigating biomarkers for use in this clinical context are ongoing but data are not yet available. We simulated a predictive biomarker informed by characteristics of whole blood lactate, which is a strong candidate biomarker [12-14]. Specifically we generated the marker as a log-normal variable with (mean, standard deviation) equal to (4.0, 2.6) and uncorrelated with risk factors in the 3121 patients with critical illness and equal to (2.5, 2.0) and uncorrelated with risk factors in the 54,526 controls without critical illness. Because of their common association with outcome, the biomarker and risk factors are correlated in the population that combines cases and controls. The log transformed marker was used in data analysis [15].

All 3121 cases that developed critical illness were included in data analyses. We generated a dataset of unmatched controls in which 3121 non-critical illness encounters were selected at random from the 54,526 available controls. For the ‘matched controls’ dataset, we selected controls who, on the basis of the clinical predictors alone, would be considered at similar risk of critical illness as the cases. Specifically, we adopted the risk score categories defined previously for these data [11], collapsing categories ≥5 into a single category. For each case, we selected a control at random from the same risk category as the case. Sampling was done without replacement in all categories but the highest where the number of controls available (n=247) was less than the number of cases (n=261). Sampling with replacement was used for this category.

Standard data analyses

We fit standard unconditional logistic regression models to derive linear combinations of risk factor predictors (denoted by X) and the biomarker (Y) [26]. Unconditional logistic regression is appropriate for the settings considered in this paper where data are frequency matched on explicit variables that can be included in the logistic model [27, chapter 6]. ROC curves were estimated empirically from the case-control data for X alone, for Y alone, and for the linear combination. These analyses were carried out separately for the cases and unmatched controls dataset and for the cases and matched controls dataset. We calculated the empirical AUC statistic for each ROC curve in addition to the ROC(0.2) statistic, which is the true positive rate corresponding to the threshold that allows a false positive rate of 0.2.

These analysis procedures were carried out for the bivariate binormal and critical illness datasets. In addition, confidence intervals were constructed and statistical tests were performed using bootstrap resampling methods for the smaller bivariate binormal dataset that included 100 cases and controls using the Stata software package [16].

Bootstrapping involved resampling separately from case and control subgroups and refitting the logistic regression models to derive the combination of X and Y in each resampled dataset. P-values were calculated with the Wald test using the bootstrap standard error estimate. We note recent concerns about the validity of P-values when testing for performance improvement [17,18] with changes in AUC or other metrics. By refitting the models we greatly improve performance of the P-values but further improvements are still warranted [19].

Analyzing cohort and matched case-control data together

In the Supplementary Data we describe in detail methods for combining cohort data on the matching factors with the matched case-control data to arrive at bias corrected estimates of ROC curves. Briefly, since ROC curves are derived from the case and control distributions of Y, the key task is to estimate the distribution of Y in the population of controls, which is distorted in the matched set by selecting controls matched to cases on X. This is accomplished by using the distribution of X for controls in the cohort to reweight the controls in the matched set. The ROC curve for the combination of (X,Y) requires the additional step of adjusting the combination derived from the matched set for associations between X and D in the population that are nullified by matching the controls to the cases on X.

Bivariate binormal scenarios

For data simulated from a bivariate binormal model, performances of the marker and the clinical risk factor variables are displayed in Figure 1. The true performances of the marker in the population are shown in the upper panels. These are also the values that are estimated from an unmatched study. We see that use of the marker on its own provides good prediction performance (upper left panel). For example, the AUC is 0.80 and when the threshold is chosen so that 20% of controls are allowed to be biomarker positive, 64% of cases are detected (i.e. ROC(0.2)=0.64). In the matched study however, performance appears to be much less (lower left panel), with AUC=0.66 and only 40% of cases detected, ROC(0.2)=0.40. See the first column of Table 1 for a range of other scenarios. Biomarker performance is severely underestimated in the matched study when standard simple estimation methods are used. We see from the first two columns of Table 1 that the bias is most severe when the marker is highly correlated with the matching factor (ρ=0.9) and less so when the correlation is small (P=0.1). Note that, by definition, if the marker is uncorrelated with the matching factor, matching has no effect on the marker distribution of selected controls and no bias results from matching.

The performance of the marker in the population may be partly attributed to its association with risk factors that are themselves associated with the outcome. How much does the marker add to the information already contained in the clinical data? This question concerning incremental value of the marker is best addressed by comparing the ROC curve for the combination of the marker and clinical risk factor variables with that for the clinical variables alone [20]. See the upper right panel of Figure 1 for an example where by adding the marker we only increase the proportion of cases detected by 6%, from 64% with the clinical predictors alone to 70% when the marker is added to them. The incremental value of the marker is small in the scenario of Figure 1. In the lower right panel of Figure 1 we see the ROC curves expected if controls are selected to exactly match the cases in regards to the clinical variables. The clinical variables on their own are completely uninformative in the matched study, a simple consequence of the design that chooses controls with the same clinical risk factor values as the cases. We see that the difference between the ROC curve for the marker combined with risk factors and the ROC curve for the clinical risk factors alone is much larger in the matched study than in the population. For example, again setting thresholds to allow 20% of controls positive, the case detection rate with the combination is 44% and when compared with that for the clinical predictors alone (20%), we see the improvement appears to be 24%. This is much larger than the actual increment in performance in the population, namely 6%. The apparent change in the AUC with addition of the marker is 0.183 in the matched study compared with the true value of 0.033. Therefore the apparent increment in performance due to the marker is severely biased by matching, this time towards over-estimation of performance. Table 1 shows numerical results for a variety of other scenarios in the column labeled ‘Combination versus Risk Factors’. In all the scenarios we studied, increments calculated from matched studies were many times larger than the true increments, ranging from 2 to more than 10 times larger than those observed in unmatched studies. The bias was more severe when the marker was highly correlated with risk factors (Table 1, column 6).

We used very large sample sizes in the simulated studies reported in Table 1 and Figure 1 in order to quantify the expected magnitudes of bias. Real studies with moderate samples have the additional issue of variability that may or may not overwhelm the biases. Figure 2 shows results for a simulated study from the same scenario as Figure 1 but with only 100 cases and 100 random or matched controls included.

We found that for the biomarker alone (Figure 2 left panels) with random unmatched controls we estimated AUC = 0.76 (95%CI=(0.69,0.83)) and ROC(0.2)=0.53 (95%CI=(0.43,0.74)) while with matched controls AUC= 0.64 (95%CI=(0.56,0.71)) and ROC(0.2)=0.35 (95%CI=(0.23,0.49)). The downward bias in the matched data estimates resulted in confidence intervals that do not contain the true values (AUC=0.80, ROC(0.2)=0.64). This contrasts with the unmatched study estimates that do contain the true values. Moreover, the confidence intervals for the matched study barely overlap with those for the unmatched study.

Considering the increment in performance gained with the biomarker (Figure 2 right panels), with unmatched data we found that the change in the AUC is 0.019 and non-significant (P=0.12) while with matched data the estimated change in the AUC is much larger and statistically significant, AUC=0.17 (P =0.003). Similarly, we found that the change in the proportion of cases detected, ROC(0.2), is only 6% (P =0.36) with unmatched data but much inflated, 19% (P =0.08) with matched data. Substantially different conclusions are drawn about the magnitude and statistical significance of the marker in the matched and unmatched studies.

Predictors of critical illness

In the dataset concerning predictors of critical illness the performance of the biomarker alone is not biased by matching. The explanation is that the biomarker is uncorrelated with the matching factors among controls in this particular setting, so matching does not distort the distribution of the marker in study controls. However, the ROC curves for the risk factors and for the combination of marker with risk factors are biased by matching. In particular the true proportion of cases detected at threshold allowing 20% of controls positive is 67% without the marker and 75% with the marker, a gain of 8% in the detection rate. In the matched analysis, however, the apparent gain is much larger, 19% (Table 2, Figure 3). Similarly, the improvement in the AUC is 0.05 with use of the marker in the population, but the matched study indicates an improvement of 0.14, almost 3 times the true value. Note that the clinical predictors are somewhat predictive in this matched study, in contrast to results for the bivariate binormal setting. The reason is that we used a categorized score for selecting the matched controls here. Since we did not match controls exactly to cases on the values of the clinical risk factors, there is some residual predictive information in the risk factors even after matching.

We implemented our proposed methodology for calculating bias corrected ROC curves from matched studies when cohort data on the matching variable is available. Results are summarized in Table 2. We found that the curves were indistinguishable from those calculated with random unmatched controls (Supplemental Data Figure 1).

Comparisons between markers

Turning now to possible effects of matching on comparisons between markers we illustrate in the left panels of Figure 4 an example where in unmatched data (and in the population) the ROC curve for marker A is higher than for marker B but the reverse holds when controls are chosen to match cases on risk factors. In a second example (right panels of Figure 4) the combination of marker A and risk factor variables has a higher ROC curve than that of the combination of marker B and risk factor variables but again the reverse holds when controls are chosen to match cases.