1 Introduction↩︎

Consider the setting of \(2n\) observations each with \(p\) covariate measurements denoted as the row vectors \(\boldsymbol{x}_1, \ldots, \boldsymbol{x}_{2n} \in \mathbb{R}^p\) and stacked rowwise into the matrix \(\boldsymbol{X}\in \mathbb{R}^{n \times p}\). Each observation is assigned (or observed to have) a treatment denoted as \(w_i\). Herein we will be considering the case of a binary treatment \(w_1, \ldots, w_{2n} \in \left\{0,1\right\}\) stacked into a column vector \(\boldsymbol{w}\) where \(w_i = 0\) will be synonymously termed “control” and \(w_i = 1\) will be synonymously termed “treatment”. Further, the measured response for each observation is binary, i.e., \(y_1, \ldots, y_{2n} \in \left\{0,1\right\}\) stacked into a column vector \(\boldsymbol{y}\) where \(y_i = 1\) will be synonymously termed a “positive response”. We will denote \(Y_1, \ldots, Y_{2n}\) stacked into a column vector \(\boldsymbol{Y}\) as the Bernoulli random variables responsible for realizing the responses. We then assume a model for the log-odds of the positive response to be linear in the conditioned covariate measurements and the treatment effect, \[\begin{align} p_i := \mathbb{P}\left(Y_i~|~\beta_w, \beta_0, \boldsymbol{\beta}, w_i, \boldsymbol{x}_i\right) = (1 + e^{-(\beta_0 + \boldsymbol{x}_i \boldsymbol{\beta}+ \beta_w w_i)})^{-1}. \end{align}\] Assuming all \(2n\) responses are independent, the likelihood function can be written as \[\begin{align} \label{eq:logistic95regression95likelihood} \mathbb{P}\left(\boldsymbol{Y}~|~\beta_w, \beta_0, \boldsymbol{\beta}, \boldsymbol{w}, \boldsymbol{X}\right) = \prod_{i=1}^{2n} p_i^{y_i} \left(1 - p_i\right)^{1 - y_i} \end{align}\tag{1}\] where \(\beta_0 \in \mathbb{R}, \boldsymbol{\beta}\in \mathbb{R}^p\) denote the \(p+1\) nuisance parameters and \(\beta_w \in \mathbb{R}\) represents the parameter of interest: the additive effect on log-odds of positive response of the treatment. Estimation of all \(p+2\) parameters is via the iterative weighted least squares procedure in standard logistic regression [1].

In both observational and experimental research, data often emerge in inherent pairings where the two units are so closely linked that they share a vast majority of their underlying characteristics. This can occur when data are measured on human siblings (especially monozygotic twins), before/after (pre-test / post-test) designs where each individual’s binary response is measured at two different times, contra-lateral studies where response measurements are made on an individual’s left/right eyes (or left/right limbs, etc), social science studies where measurements are made on both the husband and wife or on each roommate in a two-person living arrangement, etc.

Let our data now be observed as \(n\) pairs. Without loss of generality, let the pairs (sometimes referred to by the general term “strata”) be organized as observation indices \(<1,2>, <3,4>, \ldots, <2n-1,2n>\) where the first observation in the pair receives \(w=1\) (i.e., the odd numbered indices) and the second observation in the pair receives \(w=0\) (i.e., the even numbered indices). In this setting, the likelihood function in 1 is incorrect as it omits the many more covariate measurements shared among the intrapair subjects (beyond the \(p\) observed covariates \(\boldsymbol{x}_i\) measured on each subject). These unobserved characteristics induce an intrapair dependence.

To correct the likelihood function, there are many strategies. One could add different correlations for each pair \(\rho_1, \ldots, \rho_n\) but then the number of parameters expands at the same rate as the number of observations. This is called the “incidental parameter problem” and results in biased and inconsistent estimation. To shrink the parameter space to not grow with \(n\), one commonly assumes the simplified setting of one shared intrapair correlation \(\rho\) for all \(n\) pairs. Estimation then follows from the Generalized Estimating Equations (GEE) of [2], [3]. Equivalently, the effect of the unobserved covariates can be distilled into a separate intercept for each pair \(\alpha_1, \ldots, \alpha_n\). Again, the number of parameters expands at the same rate as the number of observations. To shrink the parameter space, one can assume that the intercepts are drawn randomly, e.g. from a normal distribution centered at 0 with unknown variance and this variance becomes the one nuisance parameter. Estimation then uses the generalized linear mixed model (GLMM) of [4], [5]. A fourth strategy, and our focus herein, is the Conditional Logistic Regression (CLR) model of [6] explained below. Consider pairs of observations \(<j,k>\) where the response of one observation is positive and the response of the other observation is negative, i.e. \(y_j +y_k = 1\). The two possible likelihoods for such pairs can be shown through elementary algebra [7] to be \[\begin{align} \label{eq:discordant95pairs95clr} \mathbb{P}\left(Y_j = 1, Y_k = 0~|~Y_j + Y_k = 1, \beta_w, \boldsymbol{\beta}\right) &=& \frac{1}{1 + e^{(\boldsymbol{x}_k - \boldsymbol{x}_j) \boldsymbol{\beta}+ \beta_w}}, \nonumber \\ \mathbb{P}\left(Y_j = 0, Y_k = 1~|~Y_j + Y_k = 1, \beta_w, \boldsymbol{\beta}\right) &=& \frac{1}{1 + e^{(\boldsymbol{x}_j - \boldsymbol{x}_k) \boldsymbol{\beta}+ \beta_w}} ~~~~~~ \end{align}\tag{2}\] and are aggregated together by multiplying across all such pairs as pairs are independently realized. Such pairs above are called “discordant” as their two responses are not equal (literally, their responses do not accord with one another). The likelihoods of the “concordant” pairs (i.e., those where \(y_j + y_k \in \left\{0,2\right\}\)) conditional on their sum \(Y_j + Y_k\) can be shown to be 1. Thus, observations in the concordant pairs do not factor into the likelihood function and, thus, traditionally are considered ancillary adding no information to the inference of \(\beta_w\).

The conceptual origins of CLR can be traced to [8] who used a conditional likelihood approach to isolate the effects of parental age on Down syndrome risk. By conditioning the analysis on the sibship, there was effective control for birth order and other familial confounders. This idea seemed to remain dormant until the 1970s, when it was independently reinvented and formalized. In the field of economics, [9] developed the conditional logit model to address discrete choice behavior, where the “stratum” represented a decision-maker and the observations were the various alternatives available to them. In epidemiology, researchers match cases to controls on a one-to-one or one-to-many basis, they create numerous strata. To solve the incidental parameter problem we described above, [6] adapted the partial likelihood principles of [10]’s proportional hazards model to retrospective data. By conditioning the likelihood on the total number of cases within each matched set, they found the stratum-specific intercepts \(\alpha_i\) cancel, yielding the elegant expressions found in 2 . The practical standardization of this method was achieved by [11], whose definitive work established CLR as standard methodology. While primarily associated with case-control studies, the utility of CLR extends naturally into experimental design, particularly in within-subject experiments and randomized matched-pair trials (our setting).

Our contribution is conceptually simple but powerful as it is costless. Instead of discarding the observations within the concordant pairs, we salvage the information in this data via a two stage estimation: (a) pre-modeling with the concordant data to estimate the nuisance parameters \(\boldsymbol{\beta}\) and (b) use these estimates to specify an empirical Bayes prior for a Bayesian conditional logistic regression (BCLR). This prior information on \(\boldsymbol{\beta}\) sharpens the estimation for the parameter of interest \(\beta_w\). Section 2 details our methodological contribution, Section 3 demonstrates our methods’ advantages over models common in standard practice and Section 4 concludes, offering advice to practitioners and discusses further directions.

2 Methodology↩︎

2.1 Setup↩︎

Let \(\mathcal{C}\) denote the set of concordant pairs of observations and let \(\mathcal{D}\) denote the set of disconcordant pairs of observations so that \(\mathcal{C} \cup \mathcal{D} = \left\{<1,2>, \ldots, <2n-1,2n>\right\}\) and \(|\mathcal{C}| + |\mathcal{D}| = n\). Denote \(\boldsymbol{B}_{\mathcal{C}}\) and \(\boldsymbol{b}_{\mathcal{C}}\) as an estimator and estimate of \(\boldsymbol{\beta}\) respectively. In a premodeling step, we estimate \(\boldsymbol{b}_{\mathcal{C}}\) and \(\boldsymbol{\Sigma}_{\mathcal{C}}:= \mathbb{V}\text{ar}[\boldsymbol{B}_{\mathcal{C}}]\) using data \(\boldsymbol{y}_{\mathcal{C}}\), \(\boldsymbol{X}_{\mathcal{C}}\). Each pair in \(\mathcal{C}\) is treated as a stratum, \(\boldsymbol{w}\) is dropped (as its parameter is zero) and the nuisance intercept parameter is included.

To do this premodeling step, we consider fitting with three modeling choices: LR, GEE and GLMM. The first two of these choices have systemic problems. The LR ignores the dependence structure of \(\mathcal{C}\). The GEE model is slightly incoherent as its population-level interpretation of \(\boldsymbol{\beta}\) will clash with the subject-level interpretation of \(\boldsymbol{\beta}\) in CLR. We further discuss our recommendation choice in the concluding Section 4 after first seeing the simulations of Section 3.

Regardless of the modeling choice employed, we then use the premodel’s estimates \(\boldsymbol{b}_{\mathcal{C}}\) and \(\boldsymbol{\Sigma}_{\mathcal{C}}\) as hyperparameters in one of four informative prior distributions (see following subsections). BCLR is then defined as performing Bayesian estimation of \(\boldsymbol{\beta}\), \(\beta_w\) using one of those priors, denoted \(\pi\). We will partition the posterior’s kernel as \[\begin{align} f(\beta_w, \boldsymbol{\beta}\,|\,\boldsymbol{y},\,\boldsymbol{X}) \propto p(\boldsymbol{y}_{\mathcal{D}}\,|\,\beta_w, \boldsymbol{\beta}, \boldsymbol{X}_{\mathcal{D}}) \pi(\beta_w, \boldsymbol{\beta}\,|\,\boldsymbol{y}_{\mathcal{C}},\boldsymbol{X}) \end{align}\] where the discrete likelihood above is given in 2 . We then sample \(f(\beta_w, \boldsymbol{\beta}\,|\,\boldsymbol{y},\,\boldsymbol{X})\) via a Markov Chain Monte Carlo (MCMC). By considering only the MCMC samples of \(\beta_w\), we have inference for our parameter of interest averaged over the uncertainty in the nuisance parameters \(\boldsymbol{\beta}\) [12]. For confidence set construction for \(\beta_w\), we will consider an equal-tailed \(1 - \alpha = 95\%\) credible region denoted CR\(_{\beta_w, 1 - \alpha}\), an idea first found in [13], which was noted to be a transformation-invariant decision rule by [14]. For two-sided hypothesis testing of \(H_0: \beta_w = \theta_0\), we consider a rejection defined as \(\theta_0 \notin\) CR\(_{\beta_w, 1 - \alpha}\).

The rest of this section is devoted to four different specifications of the prior \(\pi\), which describe four iterative improvements to CLR.

2.2 The Naive Prior↩︎

We employ a multivariate normal (MVN) as recommended by [15] as they considered the MVN a good starting point for a weakly informative prior, as it provides some regularization without being “stuck” in heavy tails. We use a flat prior for the treatment effect parameter and assume independence between it and the covariate effects, \[\begin{align} \label{eq:Normal95prior95without95intercept95flat95treatment} \pi\left(\beta_w, \boldsymbol{\beta}\,|\,\boldsymbol{y}_{\mathcal{C}},\boldsymbol{X}\right) = \mathcal{N}_{p+1}\left( \left[\begin{array}{c} \beta_w \\ \boldsymbol{\beta} \end{array}\right];\,\left[\begin{array}{c} 0 \\ \boldsymbol{b}_{\mathcal{C}} \end{array}\right] ,\,\begin{bmatrix} \tau^2& 0\\ 0 & \boldsymbol{\Sigma}_{\mathcal{C}} \end{bmatrix}\right) \end{align}\tag{3}\] where \(\tau^2\) is chosen to be a large positive number during implementation.

2.3 The Mixture of g’s prior↩︎

The prior in 3 does not allow for calibrating the strength of \(\boldsymbol{\beta}\) in the concordant data on the posterior. To do so, we consider the mixture-of-\(g\)’s-prior [16], [17], \[\begin{align} \label{eq:g95prior} \pi(\beta_w,\boldsymbol{\beta},\,g\,|\,\boldsymbol{y}_{\mathcal{C}},\boldsymbol{X}) &=& \pi(\beta_w, \boldsymbol{\beta}\,|\,g,\,\,\boldsymbol{y}_{\mathcal{C}},\boldsymbol{X}) \pi(g\,|\,\boldsymbol{y}_{\mathcal{C}}, \boldsymbol{X}) \nonumber\\ &=& \mathcal{N}_{p+1}\left( \left[\begin{array}{c} \beta_w \\ \boldsymbol{\beta} \end{array}\right];\,\left[\begin{array}{c} 0 \\ \boldsymbol{b}_{\mathcal{C}} \end{array}\right] ,\,\begin{bmatrix} \tau^2& 0\\ 0 & g\boldsymbol{\Sigma}_{\mathcal{C}} \end{bmatrix}\right) \times \nonumber\\ && \text{Inv}\mathcal{G}\left(g;\,\frac{1}{2},\,\frac{|\mathcal{D}|}{2}\right). \end{align}\tag{4}\] The choice of the inverse gamma hyperprior implies a marginal multivariate Cauchy distribution for \(\boldsymbol{\beta}\) (as originally found in [18]). Fitting \(g\) using both the concordant data and the size of the discordant data can provide a more flat prior than the MVN prior of 3 .

2.4 The Probability-Matching Prior (PMP)↩︎

As properly-sized hypothesis testing is essential in our upgraded CLR and retainment of our hypothesis test is defined as 0 \(\in\) CR\(_{\beta_w, 1-\alpha}\), we must ensure that this credible region has \(1-\alpha\) frequentist coverage. At low sample sizes the credible region is prone to being dominated by the prior. To mitigate this issue, we consider the PMP [19] for the parameter of interest which was proven by [20] to match frequentist coverage to within \(O(|{\mathcal{D}|}^{-1})\).

To employ the PMP in our setting, we do a preliminary step for computational simplicity: we first orthogonalize the treatment vector with respect to the observed covariates, i.e., \(\tilde{\boldsymbol{w}}_\mathcal{D} := (\boldsymbol{I}_{|\mathcal{D}|} - \boldsymbol{X}_{\mathcal{D}}(\boldsymbol{X}_{\mathcal{D}}^\top \boldsymbol{X}_{\mathcal{D}})^{-1} \boldsymbol{X}_{\mathcal{D}}^\top)\boldsymbol{w}_\mathcal{D}\). This step ensures that the treatment vector is uncorrelated with the other covariates [21]. Doing so also allows us to split the prior into (a) the uninformative PMP for the parameter of interest and (b) the informative prior for the nuisance parameters without any overlapping information [22]. The noninformative PMP treatment effect has a Jeffreys-style prior proportional to the square root of its entry on the main diagonal of the Fisher information matrix for the discordant pairs. In the CLR setting, this value is computed as \[\begin{align} I_{ww}(\beta_w, \boldsymbol{\beta}, \boldsymbol{X}_{\mathcal{D}}) = \sum_{i=1}^{|\mathcal{D}|}(\tilde{w}_{{\mathcal{D}}_i})^2 p_i(1-p_i) \end{align}\] where \(p_i\) is the likelihood of each pair respectively in 2 and thus a function of both \(\beta_w, \boldsymbol{\beta}\) and the two relevant rows of \(\boldsymbol{X}_{\mathcal{D}}\). For the informative part (the prior or the nuisance parameters \(\boldsymbol{\beta}\)), we will consider the naive MVN prior of 3 . Putting these two pieces together yields a prior proportional to \(\sqrt{I_{ww}(\beta_w, \boldsymbol{\beta}, \boldsymbol{X}_{\mathcal{D}})} \times [\text{Equation}~ \ref{eq:Normal95prior95without95intercept95flat95treatment}]\).

2.5 The Hybrid Prior: A Mixture of g’s + PMP↩︎

The final prior we consider combines the advantages of the \(g\)-prior of and the PMP and thus is a prior proportional to \(\sqrt{I_{ww}(\beta_w, \boldsymbol{\beta}, \boldsymbol{X}_{\mathcal{D}})} \times [\text{Equation}~\ref{eq:g95prior}]\).

2.6 Our Software↩︎

Our BCLR flavors require a choice of model for the concordant pairs and a choice of the prior. Thus we have 12 flavors given by the combinations {LR, GLMM, GEE} \(\times\) . Posterior inference for each flavor are implemented in our new R [23] package bclogit on CRAN. Our package makes use of use R’s Rcpp package [24] for its ability to optimize speed and stan [25] via rstan [26] for its posterior sampling engine that implements Hamiltonian MCMC [27] with the no U-turn improvement [28]. Our package bclogit also has an option for the highest posterior density (HPD) inclusion test of [29] for both contiguous and disjoint intervals using functionality from the CODA [30] and ggdist [31] packages respectively.

3 Numerical Studies↩︎

3.1 Synthetic Data↩︎

We investigate performance of our BCLR methods and compare to four canonical methods appropriate in our setting: LR, CLR (both implemented in our package), GLMM (with a random intercept for each pair) implemented in R’s package glmmTMB [32] and GEE (with exchangeable correlation matrices for each pair) implemented in R’s package geepack [33]. LR should do poorly as it ignores the dependence structure. For our method BCLR, we consider all 12 flavors (see Section 2.6). We use the rstan defaults of random initialization, 1,000 burn-in iterations, 2,000 post-burn iterations and no thinning. We measure “performance” as power when \(\alpha = 5\%\). Power was estimated by the proportion of rejections of \(H_0: \beta_w = 0\) where a rejection is defined as 0 \(\notin\) CR\(_{1 - \alpha}\) or 0 \(\notin\) HPD\(_{1 - \alpha}\). We examine \(2n = 100, 250, 500\) and \(p = 6\). Our \(\boldsymbol{X}\) is generated via sampling a \(n \times p\) matrix from \(\mathrm{U}\left(-1,\,1\right)\) and then duplicating each row with a small \(\epsilon\) drawn from \(\mathcal{N}\left(0,\,0.05\right)\). We do this to simulate natural pairing environments, such as twin or pre-post studies, in which most covariates are nearly identical within each pair. The original rows and the duplicate ones are then zipped together such that row \(\boldsymbol{x}_{i+1} = \boldsymbol{x}_{i} + \mathcal{E}\). We then shuffle the first column so that it is uncorrelated within the pair, resulting in a non-trivial inference for \(\beta_1\). This final \(\boldsymbol{X}\) remains fixed throughout all simulation settings.

We then simulate the log-odds of positive response under a linear and nonlinear model. The linear model, \(p_i = (1+\mathrm{exp}\left(-(\beta_0 + \boldsymbol{\beta}\boldsymbol{x}_{i,\cdot} + \beta_w w_i)\right))^{-1}\), serves to assess performance without model misspecification nor covariate estimator bias. The parameter values we chose are \(\beta_0 = -0.5\), \(\boldsymbol{\beta}= \left[1.25\,1.25\,\ldots\right]\). The non-linear logit-linked [34] model, \[\begin{align} p_i = \left(1+\mathrm{exp}\left(-(\text{sin}\left(\pi x_{i,1} x_{i,2}\right) + x_{i,3}^3 + x_{i,4}^2+x_{i,5}^2 + \beta_w w_i)\right)\right)^{-1}, \end{align}\] is designed to represent a more complex, realistic scenario. Here, fitting the linear model is merely a best linear approximation of the function above. In both probability response models, we chose \(\beta_w = 0.5\), which was picked to reveal power differences among the methods for our range of considered sample sizes. (We additionally consider \(\beta_w = 0\) to assess the size of the tests in each setting).

We further consider two covariate settings where (1) only the first covariate is observed or (2) the first two covariates are observed. The setting of (1) simulates the case in which most of the predictive power lies in latent covariates that are virtually identical within pairs, while considering only a variable independent of the pair structure that has a significant effect on the response. The setting of (2) allows the model to peer into the underlying pair structure.

To reiterate, we have 3 (sample sizes)\(\times\) 12+4 (our model flavors + canonical models) \(\times\) 2 (power vs size) \(\times\) 2 (covariate settings) \(\times\) 2 (response models) = 384 settings. For each setting, we run \(N_{sim} = 10,000\) simulation iterations. In each iteration, we draw \(\boldsymbol{w}\) by uniformly setting each pair to either \((0,1)\) or \((1,0)\). Then we draw the \(y_i\)’s from independent Bernoullis with parameter \(p_i\) given by the linear / nonlinear response models. Then we record the estimate of \(\beta_w\) and whether the null hypothesis is rejected at a significance level of 5%. Over the \(N_{sim}\) simulations, we approximate power (when \(\beta_w = 0.5\)) and size (when \(\beta_w = 0\)) as the proportion of rejections.

For all 12 flavors of BCLR, we display comprehensive power and size results in Tables [tbl:tab:power95sub] and [tbl:tab:size95sub] respectively for the CR method. Tables ¿tbl:tab:power95sub95HPD? and [tbl:tab:size95sub95HPD] in the supplementary information show power and size using the disjoint HPD testing method. The HPD testing method was found to have lower power. This is possibly due to our posteriors being skewed as HPD intervals overextend on the short tail. We also include posterior estimates of \(g\) (for the \(g\) and hybrid priors) in Tables [tbl:tab:g95values950465] and [tbl:tab:g95values950] in the supplementary information for the cases of \(\beta_w = 0.5\) and \(\beta_w = 0\) respectively. Additionally, we include results on two more performance metrics in the supplementary information: the approximate mean squared error (MSE) in Table 3 and approximate coverage in Table 4. We use the posterior expectation (estimated by the average of the MCMC samples) as the definition of the point estimate \(b_w\) and compute the squared error \((b_w - \beta_w)^2\). MSE is then approximated as the average of those squared errors over the \(N_{sim}\) iterations. In each iteration we record if \(\beta_w \in \text{CR}_{\beta_w, 95\%}\). The proportion of these set comparisons being true over the \(N_{sim}\) iterations approximates the coverage.

To simplify the results illustration in the manuscript, we will display only the combination choice we believe is “best” as determined by (1) being mostly properly sized across the various simulation settings and (2) being numerically stable. The winning flavor we believe to be LR+naive and henceforth termed BCLR in this section’s results henceforth. Although GEE+\(g\) had the highest power on average over our settings, both GEE and \(g\) we believe have separate problems which prevents us from recommending it to a practitioner. GEE suffered from occasional unacceptably large MSE. In some of our iterations, a covariate perfectly predicts the outcome. Here, the iterative reweighted least squares algorithm will fail to converge as the parameter estimates run towards infinity. Additionally, the sandwich estimator used for robust standard errors in GEE can produce singular matrices. In our software implementation, we have many types of robust fallbacks to LR. GLMM also experienced such instabilities, but less so. The \(g\) prior (and hence the hybrid as well) experiences inflated size. This is likely due to \(g\) being large, yielding flat priors. The PMP is observed to be improperly sized being both too small and too large. Although the naive prior cannot weight the concordant data explicitly like in a \(g\) prior (a weighting we conjectured at the outset to be critical), the sample size \(|\mathcal{C}|\) is built implicitly into \(\boldsymbol{\Sigma}_{\mathcal{C}}\), its variance estimate.

The power and size results for BCLR and the canonical models are found in Table 1 for each sample size \(\times\) covariate setting \(\times\) response model. BCLR’s power performance dominated CLR’s power performance in all settings. BCLR is also on par with the canonical methods GEE and GLMM in the linear logodds response setting. However, in the nonlinear log odds response setting, BCLR dominates GEE and GLMM, methods that rely heavily on the linear model assumption. BCLR uses matched pair covariate data differences as seen in 2 which have a long history of being robust in nonlinear response models. These differences reduce bias because the nonlinearities “cancel out” or are minimized [35], [36]. As expected LR underperforms as it is misspecified.

3.2 Real Data Application↩︎

The Framingham Heart Study [37], [38] is one of the most significant and longest-running longitudinal studies in medical history, dating to 1948 with its initial \(\approx\)​5,000 adult subjects. We make use of the study’s publicly available dataset as found in the R package riskCommunicator [39].

There are many endpoints in this study. We use the binary endpoint “prevalent coronary heart disease” where the positive level was defined as preexisting angina pectoris, myocardial infarction (hospitalized, silent or unrecognized), or coronary insufficiency (unstable angina). We considered eight covariates: total cholesterol (numeric), systolic blood pressure (numeric), diastolic blood pressure (numeric), heart rate (numeric), number of cigarettes smoked per day (count), body mass index (numeric), diabetes (binary), and use of anti-hypertensive medication (binary). These covariates varied longitudinally across patients and had low missingness. We defined a pair as the first measurement (\(w = 0\)) and the last measurement (\(w = 1\)) made on the same patient. Intermediate measurements were dropped. We then dropped missing data and retained only the patients who had these two measurements. For these patients, the time between first and last measurement was always 10-14 years elapsed, and 82% of the subjects had 12 years elapsed. Hence, the treatment effect of interest is the additive effect of approximately 12 years of aging on the log-odds of a prevalent coronary heart disease event. As aging increases the risk of heart disease, the treatment effect is real and positive.

After this initial data curation, there were \(n=2,971\) patient pairs. In the approximate 12yr elapse, most patients’ endpoints remained the same; there were \(|\mathcal{C}|\) = 2,718 concordant pairs and only \(|\mathcal{D}|\) = 253 discordant pairs. Hence, this is an example where CLR’s performance is likely to suffer as 91% of the subject data is discarded. Conversely, this is an example where our BCLR will perform well as we use that 91% of the data to estimate accurately the approximate linear effect of the eight informative covariates on the log odds of prevalent coronary heart disease. This is exactly what is found in the observed estimates and significant tests (see Table 2) where CLR failed to find the real effect and BCLR comfortably found the effect. As an analogy, in the FHS data BCLR is using covariate adjustment to find a noisy effect while CLR is using a model without covariates at all.

4 Discussion↩︎

We developed an improvement to the classic CLR model in the setting of estimating a binary treatment effect along with control covariates. The source of our improvement is the information-salvaging of the concordant pairs which are traditionally viewed as ancillary to the likelihood. Our method BCLR is a two-step procedure: (1) we use the information about the control covariates found in the concordant pairs by creating a pre-model, and (2) we use this information as a prior in a Bayesian implementation of CLR, which we call BCLR. Thus we provide a principled way to recover efficiency that is lost in the standard workflow.

We investigated many different flavors of pre-models, but we found through simulation that the most reliable was the prosaic LR whose coefficient estimates and variance-covariance matrix of the coefficients’ estimator were rolled into a prosaic multivariate normal prior. By “most reliable” we mean it provided tests of the parameter of interest that were properly sized and protected from numeric instability, which is a risk due to the structure of the concordant pair data.

Our BCLR method was found to be superior to CLR as measured by the power of the test of no treatment effect, and it was more numerically stable. (Our method also translates to improvements for inference on the nuisance parameters, but this is unshown herein). The BCLR method was also competitive with other canonical models, such as GEE and GLMM, although GEE seems to perform better. In the case of a non-linear log odds of positive response model, BCLR was found superior. We hesitate to state numerically how much better BCLR’s power is over the classic CLR as improvements are setting-dependent. But we can say that in the simulated synthetic data settings of Section 3.1 where we examined both linear and nonlinear log-odds of positive response models, improvements ranged between 58-116% for \(n=100\), between 7.1-8.1% for \(n=250\) and between 2.0-2.9% for \(n=500\). Additionally, in the case of a real dataset (the FHS of Section 3.2), BCLR was able to decisively identify a known effect while CLR could not.

As for GEE performing better than BCLR in the linear case, we do not believe this to be an apples-to-apples comparison as GEE provides estimates of the population-level effects, thus answering the question, “What is the change in the log odds of positive response for the whole population”? where BCLR provides estimates of subject-specific effects, thus answering the question, “What is the change in the log odds of positive response for the average specific subject”? and hence the parameters are defined differently.

We have many further paths to explore. We believe our methods could be extended to strata with more than two subjects, continuously-measured treatments, and categorical treatments, and it is our plan to update the software to make such accommodations. We would also like to explore the flavors more, e.g., we can correct the large MSE due to numerical instability in GEE and GLMM by reverting to LR for those few cases. There are also other flavors worth exploring such as power priors [40]. Additionally, we believe we can prove the finite efficiency of BCLR vs CLR for the estimation of the additive treatment effect parameter. This would involve making assumptions on the distribution of covariates and the distribution of the parameters \(\boldsymbol{\beta}\). Then we would need to carefully compute the variance of the naive+LR estimator (and possibly other flavors), which is likely doable.

5 Conflicts of interest↩︎

The authors declare that they have no competing interests.

6 Funding↩︎

This research was supported by Grant No 2018112 from the United States-Israel Binational Science Foundation (BSF).

7 Data availability↩︎

The data and code to reproduce the results herein can be found at https://github.com/Tennenbaum-J/bclogit_package_and_paper_repo/.

8 Author contributions statement↩︎

J.T. and A.K. shared in conception and manuscript writing and review. J.T. did the majority of the numerical studies. J.T and A.K. shared in the development of the CRAN package bclogit.

9 Acknowledgments↩︎

The authors thank the support of Grant No 2018112 from the United States-Israel Binational Science Foundation (BSF). We thank the Gemini CLI for aiding with software development.

References↩︎