Accounting for overdispersion and clustering in binomial data from N-of-1 trials

Technical Summary: Accounting for Overdispersion and Clustering in Binomial Data from N-of-1 Trials

Problem Statement
N-of-1 trials are patient-centered, randomized, double-blinded crossover experiments designed to evaluate interventions for individual patients with chronic conditions. While the primary objective is to determine efficacy for a specific patient, there is significant value in pooling results across a set of N-of-1 trials to obtain population-level estimates of efficacy and to assess the success of treatment individualization. However, standard statistical methods for pooling binomial data often fail to account for two critical features inherent in N-of-1 designs:

1. Overdispersion: The variance of the binary outcomes often exceeds the mean-variance relationship prescribed by the standard binomial distribution.
2. Hierarchical Clustering: Measurements are correlated within patients due to repeated treatment episodes (within-patient correlation).

While Bayesian methods have been proposed for this context, particularly for binary outcomes, there is a lack of frequentist approaches that rigorously address both overdispersion and clustering simultaneously. This paper addresses this gap by presenting and comparing two frequentist analytical strategies for pooling binomial data from N-of-1 trials.

Methodology
The authors propose two distinct frequentist approaches to estimate the pooled proportion ($\theta$) of treatment preference or success.

• Method 1: Quasi-Likelihood (QL) Approach (Overdispersion Only)
  This method models the data using a Beta-Binomial distribution to account for overdispersion ($\rho$) but treats patients as independent units regarding the clustering structure.

  • Model: $X_j \sim \text{BetaBinom}(n_j, \theta, \rho)$, where $X_j$ is the number of successes for patient $j$ out of $n_j$ episodes.
  • Estimation: The pooled proportion $\theta$ is estimated using a quasi-likelihood estimator: $\hat{\theta}_{QL} = \frac{\sum c_j X_j}{\sum c_j n_j}$, where weights $c_j$ depend on the overdispersion parameter $\rho$.
  • Variance Estimation: Standard errors are derived using the Jackknife method to ensure consistency.
  • Estimating $\rho$: The paper evaluates several estimators for the intra-cluster correlation (ICC) parameter $\rho$, including Analysis of Variance (AOV), Fleiss-Cuzick (FC), Pearson (PE), Extended Quasi-Likelihood (EQL), Pseudo-Likelihood (PL), and Method of Moments (MM). The authors note that closed-form variance formulas for some estimators (like AOV) can yield negative values or bias, leading them to rely on Jackknife-based variance estimation for comparison.

• Method 2: Random Effects Approach (Overdispersion and Clustering)
  This method, adapted from Molenberghs et al., utilizes a hierarchical random effects model to simultaneously account for overdispersion and the correlation among repeated measures within patients.

  • Model: $X_j \sim \text{Binom}(n_j, \theta_j)$, where $\theta_j = \pi_j \kappa_j$. Here, $\pi_j$ (modeled via a Beta distribution) accounts for remaining overdispersion, and $\kappa_j$ (a random effect derived from a normal distribution) accounts for the hierarchical clustering.
  • Estimation: Parameters are estimated via Maximum Likelihood (ML) using numerical integration. The final pooled estimate is derived as the product of the expected values of the random components: $\hat{\theta} = \hat{E}(\pi_j) \times \hat{E}(\kappa_j)$.
  • Approximation: The expectation of the sigmoid function of the normal random variable is approximated using a standard normal CDF with a scaling factor ($\lambda \approx 0.588$).

Simulation Study and Results
The authors conducted a simulation study with 1,000 iterations to compare the performance of these methods under various scenarios of true proportion ($\theta$) and overdispersion ($\rho$), with varying sample sizes ($J$ patients and $n_j$ episodes).

• Performance of $\rho$ Estimators: In Method 1, the Extended Quasi-Likelihood (EQL) estimator for $\rho$ exhibited high Mean Squared Error (MSE) and bias, particularly when $\theta = 0.50$. However, despite the poor performance of the EQL estimator for $\rho$, the resulting pooled proportion estimate ($\hat{\theta}$) performed comparably to, or sometimes better than, other estimators. The authors attribute this to a cancellation effect where the bias in $\hat{\rho}$ reduces the variance contribution of the $\rho$ term in the Delta method expansion of $\text{Var}(\hat{\theta})$.
• Performance of Method 1 vs. Method 2:
  • MSE: Method 1 estimators generally showed stable MSEs. Method 2 showed competitive MSEs in some scenarios but performed worse when $\theta = 0.75$.
  • Coverage Probability: Method 1 maintained coverage probabilities close to the nominal 95% level across all scenarios. Method 2 (Random Effects) showed significantly poor coverage (as low as ~78%) when $\theta = 0.75$. The authors attribute this to the skewed sampling distribution of variance estimates in Generalized Linear Mixed Models (GLMMs), which renders symmetric confidence intervals inappropriate in those specific scenarios.
• Sample Size Effects: Increasing the number of patients ($J$) reduced the MSE more effectively than increasing the number of episodes per patient ($n_j$).

Real Data Applications
The methods were applied to three real-world N-of-1 trial datasets:

1. Amitriptyline for Fibromyalgia: All methods produced similar pooled estimates of preference ($\hat{\theta} \approx 0.67$). No significant overdispersion was detected (confidence intervals for $\rho$ included zero).
2. NSAIDs vs. Paracetamol for Osteoarthritis: Results were consistent across methods ($\hat{\theta} \approx 0.42$ for total diary score). No significant overdispersion was found.
3. Theophylline for Irreversible CAL: Pooled estimates were consistent ($\hat{\theta} \approx 0.51$). Again, no significant evidence of overdispersion was detected.

In all real data examples, the point estimates for the pooled proportion were remarkably consistent across the different methods, despite variations in the estimated overdispersion parameters.

Significance and Claims
The paper claims to provide the first frequentist-based methods for pooling binary data from N-of-1 trials that explicitly address overdispersion and hierarchical clustering.

• Methodological Contribution: The authors introduce a quasi-likelihood framework for overdispersion and adapt a random effects framework for simultaneous overdispersion and clustering.
• Interpretational Distinction: The paper highlights a conceptual difference between the two methods. Method 1 yields "population-averaged" (PA) estimates, while Method 2 yields "subject-specific" (SS) estimates. Given that N-of-1 trials are inherently focused on individual differences, the authors suggest that the interpretation of Method 2 (SS) may be slightly more appropriate for the N-of-1 context, despite its computational complexity and occasional coverage issues in skewed scenarios.
• Practical Utility: The study demonstrates that while the estimation of the overdispersion parameter ($\rho$) can be unstable (particularly with EQL), the resulting pooled proportion estimates ($\hat{\theta}$) are robust across different estimation strategies. The authors conclude that for practical applications, the choice of method may depend on whether the inferential goal is population-level or individual-level, and that Method 1 offers a simpler, robust alternative when coverage properties of Method 2 are a concern.