Sample size and power determination for assessing overall SNP effects in joint modeling of longitudinal and time-to-event data

1. Problem Statement

The paper addresses a critical gap in the statistical design of genetic association studies involving joint modeling of longitudinal biomarkers (e.g., HbA1c levels) and time-to-event outcomes (e.g., onset of diabetic complications).

• Context: While methods for jointly modeling these data types exist, there is a lack of rigorous sample size and power calculation tools specifically tailored for Single Nucleotide Polymorphisms (SNPs).
• Specific Challenges:
  • Genetic Complexity: Unlike binary treatment groups, SNPs are categorical (0, 1, or 2 copies of a risk allele) and often have unbalanced genotype distributions (depending on Minor Allele Frequency, $p$), which affects statistical efficiency.
  • Mechanism of Action: Genetic variants influence event risk through two pathways: direct effects (on the hazard) and indirect effects (mediated through the longitudinal biomarker). Standard survival analysis formulas (e.g., Schoenfeld) do not account for this mediation.
  • Testing Direction: Genetic effects can be bidirectional (increasing or decreasing risk), necessitating two-sided hypothesis testing, unlike the often directional assumptions in treatment trials.
• Goal: To derive a closed-form sample size formula to test the overall SNP effect ($\gamma_g + \alpha\beta_g$) within a joint modeling framework and to validate its accuracy.

2. Methodology

A. Statistical Framework

The authors utilize a joint modeling framework linking two sub-models:

1. Longitudinal Sub-model:
   $$Y_i(t) = \eta_i(t) + \epsilon_i(t)$$
   $$\eta_i(t) = (\beta_0 + b_{i0}) + (\beta_1 + b_{i1})t + \dots + \beta_g \text{SNP}_i$$
   Here, $\beta_g$ represents the SNP's effect on the longitudinal trajectory, and $b_i$ are random effects.
2. Survival Sub-model (Proportional Hazards):
   $$\lambda_i(t) = \lambda_0(t) \exp\{ \gamma_g \text{SNP}_i + \alpha \eta_i(t) \}$$
   Here, $\gamma_g$ is the direct effect of the SNP on the hazard, and $\alpha$ links the true longitudinal trajectory to the event risk.

The Overall Effect: The total effect of the SNP on the event risk is defined as the sum of the direct and indirect effects:
$$\theta_g = \gamma_g + \alpha\beta_g$$
The hypothesis test is $H_0: \theta_g = 0$ vs. $H_a: \theta_g \neq 0$.

B. Derivation of the Sample Size Formula

Building on Schoenfeld (1983) and Chen et al. (2011), the authors derive the number of events ($D$) required to achieve power $1-\tilde{\beta}$ at significance level $\tilde{\alpha}$:
$$D = \frac{(Z_{1-\tilde{\beta}} + Z_{1-\tilde{\alpha}/2})^2}{2pq(\gamma_g + \alpha\beta_g)^2}$$

• Key Adaptation: The formula replaces the binary treatment variance with the SNP genotype variance ($2pq$, where $p$ is the minor allele frequency).
• Power Inversion: The formula is inverted to estimate power given a fixed sample size and event count.

C. Estimation Strategy (Two-Stage Approach)

Recognizing that full joint modeling via the Expectation-Maximization (EM) algorithm is computationally intensive for large-scale genetic studies, the authors propose a two-stage method for practical power estimation and validation:

1. Stage 1: Fit the full longitudinal model to estimate $\hat{\beta}_g$.
2. Adjustment: Compute adjusted longitudinal responses $y'_i(t) = y_i(t) - \hat{\beta}_g \text{SNP}_i$ to remove SNP-related variation.
3. Stage 2: Fit a reduced longitudinal model to estimate the adjusted trajectory $\hat{\eta}'_i(t)$, then fit a survival model: $\lambda_i(t) = \lambda_0(t) \exp\{ \theta_g \text{SNP}_i + \alpha' \hat{\eta}'_i(t) \}$.
4. Testing: Assess significance of $\theta_g$ under $H_0: \theta_g = 0$.

3. Key Contributions

1. Closed-Form Formula: Derivation of the first closed-form sample size formula specifically for testing the overall effect of a SNP (direct + indirect) in joint models.
2. Handling Genetic Specifics: Explicitly accounts for SNP dosage coding (0, 1, 2), unbalanced genotype frequencies ($2pq$), and the necessity of two-sided testing.
3. Computational Efficiency: Proposes a robust two-stage estimation approach that avoids the heavy computational burden of full joint likelihood maximization, making it feasible for Genome-Wide Association Studies (GWAS).
4. Practical Tool: Development of an interactive Shiny app (PowerSNP) to facilitate sample size and power calculations for researchers.

4. Results

Simulation Studies

• Accuracy: The proposed formula accurately predicts the required number of events ($D$) across various scenarios, including different effect sizes ($\alpha, \beta_g, \gamma_g$), random effect structures, and significance levels.
• Robustness: Power estimates remain stable across different random effect covariance structures.
• Limitations:
  • Measurement Frequency: The formula performs well with sufficient measurement frequency but underestimates power when the number of visits is very small (interval censoring issues).
  • Model Misspecification: If the longitudinal model is misspecified (e.g., fitting a linear model when the true trajectory is quadratic), empirical power deviates from calculated power, highlighting the need for correct trajectory modeling.
  • Bias: Ignoring the longitudinal trajectory in survival models leads to biased estimates of the direct SNP effect ($\gamma_g$), whereas the proposed two-stage and joint methods yield unbiased results.

Retrospective Analysis: DCCT Data

The authors applied the method to the Diabetes Control and Complications Trial (DCCT) to assess the power of detecting SNP effects on retinopathy.

• Findings: For both the conventional and intensive treatment arms, the DCCT sample size (approx. 526–527 subjects) is insufficient to achieve adequate power (80%) at conventional GWAS significance levels ($\tilde{\alpha} = 5 \times 10^{-8}$) for modest effect sizes.
• Sensitivity: The intensive arm (fewer events) showed lower power and higher sensitivity to minor allele frequency ($p$) and significance thresholds compared to the conventional arm.
• Implication: Detecting modest genetic effects mediated by longitudinal biomarkers in stratified treatment groups requires significantly larger sample sizes than those currently available in many existing cohorts.

5. Significance and Conclusion

• Study Design Optimization: The paper provides essential tools for researchers to determine feasible sample sizes for genetic studies involving longitudinal mediators, preventing underpowered studies that fail to detect true associations.
• Methodological Advancement: By extending joint modeling theory to include specific genetic effect calculations, the authors bridge the gap between complex statistical modeling and practical study design.
• Real-World Impact: The application to DCCT data serves as a cautionary example, demonstrating that even well-conducted historical trials may lack the power for modern GWAS-level discovery when analyzing complex genetic pathways.
• Future Directions: The authors note that while the formula is robust, it assumes non-informative censoring and correct specification of the longitudinal trajectory. Future work should address model misspecification and informative dropout.

In summary, this paper offers a rigorous, validated, and computationally efficient framework for planning genetic association studies where the mechanism of action involves both direct and longitudinal-mediated pathways.