What is known?

• • Quasi-experimental methods, including instrumental variable (IV) approaches, are underutilized in IPD meta analyses of observational data despite their increased utilization in single studies.

• • One barrier to their implementation in pooled analyses is that the classical F-statistics are unavailable in mixed effects models, leaving analysts without a practical way to assess instrument strength in these settings.

What is new?

• • We adapt and evaluate a frequentist mixedeffects two-stage residual inclusion (2SRI) framework for survival IPD meta analyses, extending IV methods to handle study-level and temporal clustering with Cox proportional hazards models for time-to-event outcomes.

• • We propose the Wald X2 statistic as a practical instrument-strength diagnostic and show that under weak confounding naive models dominate, while under moderate-to-strong confounding with realized Wald X2 exceeding 150–200, mixed-effects 2SRI substantially reduces bias and achieves nearnominal coverage.

• • We provide empirical guideposts linking realized first-stage instrument strength to expected estimator performance, giving analysts decision rules for when mixed-effects 2SRI will outperform conventional approaches.

Potential impact for RSM readers

• • These guideposts offer a concrete diagnostic that may lower the practical barrier to adopting IV-based 2SRI methods in IPD meta analyses of observational data.

2.1 Instrumental variable models

IV methods address unmeasured confounding by exploiting quasi-random variation in treatment assignment or exposure. A valid IV must satisfy three assumptions: (i) relevance (the IV predicts treatment), (ii) independence (the IV shares no unmeasured common causes with the outcome), and (iii) exclusion restriction (the IV affects the outcome only through treatment).Reference Angrist, Imbens and Rubin²¹

Figure 1 illustrates these assumptions as a directed acyclic graph (DAG). Of the three, only the relevance assumption (i) is directly testable: the strength of the instrument–exposure association can be evaluated using first-stage diagnostics, most commonly the partial F-statistic in linear models.Reference Staiger and Stock²² The dashed red arrows depict the two pathways that must be assumed absent—a direct path from the instrument to unmeasured confounders (assumption ii, independence) and a direct path from the instrument to the outcome, bypassing treatment (assumption iii, exclusion restriction). Neither assumption (ii) nor assumption (iii) can be empirically verified from the observed data. The plausibility of both must therefore be defended through contextual reasoning and domain knowledge, including careful consideration of the processes through which the instrument is thought to operate, and careful consideration of alternative pathways through which the instrument could plausibly influence the outcome or correlate with unmeasured characteristics of individuals or study sites.

Figure 1

Directed acyclic graph (DAG) illustrating the three assumptions required for a valid instrumental variable (IV): (i) relevance, (ii) independence, and (iii) exclusion restriction. Dashed red arrows indicate the two pathways that must be assumed absent but cannot be empirically verified from observed data.

2.2 Mixed effects two-stage residual inclusion

In medicine, when the outcome of interest is survival, researchers conducting IPDMAs may adopt a mixed-effects Cox proportional hazards model in a one-step IPDMA.Reference de Jong, Moons and Riley²³

Formally, the naive model (without accounting for endogeneity) is

(2.1) $$ \begin{align} \log h_{ijk}(t) = \log h_0(t) + \beta_T T_{ijk} + \beta_{\text{X}}^{\mathsf{T}} X_{ijk} + u_{j}^{(S)} + v_{k}^{(S)}, \end{align} $$

where $h_{0}(t)$ is the unspecified baseline hazard, $T_{ijk}$ is the treatment for individual i in study j entering at time k (e.g., month), $X_{ijk}$ is a vector of observed covariates, $\beta _T$ is the treatment effect, $\beta _{\text {X}}$ is a vector of covariate effects, and $(u_{j}^{(S)}, v_{k}^{(S)})$ are random intercepts for study and time. The random intercepts for study and time capture cluster-level heterogeneity in baseline risk; their interpretation in the context of survival IPDMAs is detailed in Section 3.2.

Among IV approaches, 2SRI is appropriate when the second-stage model is nonlinear, as in the Cox proportional hazards model for survival data, where alternative approaches, such as two-stage least squares (2SLS) or two-stage predictor substitution (2SPS), can be biased.Reference Terza, Basu and Rathouz²⁰ In 2SRI, the first-stage models treatment as a function of the instrument and observed covariates, then residuals are calculated, and both treatment and residuals enter the second-stage survival model. Formally, the mixed-effects 2SRI model is estimated in two stages.

First, a mixed-effects logistic regression models treatment assignment:

(2.2) $$ \begin{align} \operatorname{logit}\Pr(T_{ijk}=1) = \alpha_{0} + \alpha_{1} Z_{ijk} + \alpha_{2}^{\mathsf T} X_{ijk} + u_{j}^{(T)} + v_{k}^{(T)}, \end{align} $$

where $T_{ijk}$ is the treatment, $Z_{ijk}$ is the instrument, $X_{ijk}$ is the vector of observed covariates, and $(u_{j}^{(T)}, v_{k}^{(T)})$ are random effects for study and time.

The predicted probability is then estimated:

(2.3) $$ \begin{align} \widehat p_{ijk} = \operatorname{logit}^{-1}\!\Big( \hat\alpha_{0} + \hat\alpha_{1} Z_{ijk} + \hat\alpha_{2}^{\mathsf T} X_{ijk} + \hat u_{j}^{(T)} + \hat v_{k}^{(T)} \Big). \end{align} $$

The residuals are calculated:

(2.4) $$ \begin{align} \hat r_{ijk} = T_{ijk} - \widehat p_{ijk}. \end{align} $$

The second stage is a mixed-effects Cox model with residuals:

(2.5) $$ \begin{align} \log h_{ijk}(t) = \log h_0(t) + \beta_T T_{ijk} + \delta \hat{r}_{ijk} + \beta_{\text{X}}^{\mathsf{T}} X_{ijk} + u_{j}^{(S)} + v_{k}^{(S)}. \end{align} $$

Here, the same covariates (X) from the first stage (2.2) enter the second-stage model. The random intercepts $u_j^{(S)}$ and $v_k^{(S)}$ mirror the structure of the first stage but are estimated separately; the superscripts (T) and (S) distinguish the two stages. The coefficient $\beta _T$ represents the causal effect of treatment on the log-hazard of the outcome, while $\delta $ captures the adjustment for unmeasured confounding through the first-stage residual.

2.3 Diagnosing first-stage strength in nonlinear mixed-effects models

In linear first-stage models, weak-instrument diagnostics are based on the partial F-statistic testing whether the instrument predicts treatment $H_{0}:\alpha _{1}=0$ . With a single instrument, the partial F equals the square of the t-statistic on the instrument coefficient, and the conventional “strong IV” rule of thumb is $F \gtrsim 10$ .Reference Staiger and Stock²² Software packages like ivreg,Reference Fox, Kleiber, Zeileis and Suschnigg²⁴ fixest,Reference Berge, Krantz, McDermott, Lenth and Butts²⁵ and plmReference Croissant and Millo²⁶ provide built-in tests for linear models.

When the treatment variable is binary, requiring a logistic first-stage model, no direct analog to the F-statistic exists. Recommended alternatives include the Wald statistic for the instrument coefficient or likelihood ratio tests.Reference Wooldridge²⁷ However, no standard R packages currently provide weak-instrument diagnostics for logistic mixed-effects first-stage models. Of the available alternatives, the Wald statistic is the most direct analog to the F-statistic.

We can test the same null $H_{0}:\alpha _{1}=0$ with a Wald statistic. Let $\hat \alpha _{1}$ be the estimated IV coefficient and $\mathrm {SE}(\hat \alpha _{1})$ its standard error; define

$$\begin{align*}z_{\mathrm{IV}}=\frac{\hat\alpha_{1}}{\mathrm{SE}(\hat\alpha_{1})}, \qquad W=z_{\mathrm{IV}}^{\,2}. \end{align*}$$

Under $H_{0}$ , W follows a $\chi ^{2}_{1}$ distribution, from which the p-value can be obtained as

$$\begin{align*}p_W = 1 - F_{\chi^{2}_{1}}(W), \end{align*}$$

where $F_{\chi ^{2}_{1}}$ is the cumulative distribution function of the chi-squared distribution with one degree of freedom. We report W (df $=1$ ) as the first-stage relevance diagnostic for the logistic mixed model. However, the operating characteristics of this W diagnostic in IPDMA settings are not well-documented. Our simulations therefore systematically vary design knobs that affect instrument strength (instrument–exposure association strength and cluster-level instrument–exposure allocation patterns), and then stratify 2SRI performance metrics by the realized first-stage Wald $\chi ^2$ . This links a quantity analysts can compute in practice to expected bias, coverage, and precision in mixed-effects 2SRI.

2.4 Hypotheses

Based on the theoretical framework and existing literature, we propose the following hypotheses:

1. 1. Under weak confounding, mixed-effects 2SRI will not produce improved causal estimation over mixed-effects Naïve Cox.

2. 2. Under moderate to strong confounding, mixed-effects 2SRI will improve coverage of the true effect estimate.

3. 3. Under moderate to strong confounding, relative bias will be lower in mixed-effects 2SRI.

4. 4. Under weak to moderate confounding, mixed-effects 2SRI will have higher variance and worse MSE than naïve survival.

3.1 Aim

The primary aim of this study is to evaluate the performance of mixed-effects 2SRI model, relative to a mixed-effects naïve Cox proportional hazards model, and whether it yields valid estimates of the treatment log-hazard effect in a one-step IPDMA with two random effects under varying levels of unmeasured confounding.

The secondary aim is to characterize the relationship between first-stage Wald statistics and 2SRI performance across simulation scenarios, and explore whether approximate regions of Wald values correspond to acceptable levels of bias and inference. In order to do this, we systematically vary instrument–exposure association levels and instrument–exposure prevalence patterns across studies while holding the IPDMA structure fixed.

3.2 Data-generating mechanisms

We generate data from parametric models anchored to the sizes and structure of a pooled Ebola Virus Disease dataset held by IDDO.Reference Merson and Strudwick³¹ We fix the number of studies at 10, with study-specific sample sizes matching those observed in the Ebola dataset (100, 148, 111, 738, 1096, 185, 1230, 214, 90, 445).Reference Merson and Strudwick³¹ Sample sizes are fixed to ensure external relevance to infectious disease outbreak contexts, where studies are generally small, and computational feasibility. While sample sizes are fixed, the number of events arises from the competing risks process.

Treatment assignment is generated from a logistic regression model, including instrument assignment, unmeasured confounding, and covariates. Survival times are generated from a Weibull model with log-hazard determined by treatment, covariates, unmeasured confounding, and random effects for study and calendar month. In Ebola and similar severe infectious diseases, patients typically die or are discharged within 30 days of admission. Thus, our Weibull parameters reflect this short survival window typical of real-world outbreaks. The true treatment effect is set to log-hazard $= -0.4$ .

We vary three design factors in a fully factorial grid:

• • Instrument assignment pattern (three levels): balanced (50%), moderately unbalanced (70%–30%), and a distribution based on a real-world Ebola treatment distribution (76%, 97%, 98%, 5%, 77%, 5%, 52%, 10%, 29%, 15%).Reference Merson and Strudwick³¹

• • Instrument–treatment association strength (four levels, log-odds scale): $iv_{\text{strength}} \in \{0.3, 0.5, 0.8, 1.0\}$

• • Unmeasured confounding (four levels, affecting both treatment and outcome):

  • – Weak: $(\beta _U^{(t)}, \beta _U^{(d)}) = (0.5, 0.5).$

  • – Moderate: $(1.0, 1.0).$

  • – Strong: $(1.5, 1.5).$

  • – Very strong: $(2.0, 2.0).$

3.2.2 Instrument–exposure association (first-stage design strength)

We manipulate the design-time instrument–exposure association via the coefficient on the instrument in the first-stage logistic mixed model:

$$\begin{align*}\operatorname{logit}\Pr(T{=}1)=\alpha_0+\underbrace{\alpha_1 Z}_{\text{design strength}}+\text{covariates}+u_{\text{study}}+v_{\text{time}}. \end{align*}$$

Varying $\alpha _1$ directly controls the odds ratio $\exp (\alpha _1)$ for Z and thus the induced shift in treatment propensity (holding other terms fixed). This follows the standard approach used in simulation studies that tune instrument relevance through the first-stage slope on Z (e.g., $\alpha _1\in \{1,2,3\}$ Reference Koladjo, Escolano and Tubert-Bitter³²). We use a slightly finer grid to improve resolution:

$$\begin{align*}\alpha_1 \in \{0.3,\ 0.5,\ 0.8,\ 1.0\}. \end{align*}$$

Realized first-stage diagnostics (e.g., Wald $\chi ^2$ ) are reported descriptively but are not used to set design strength.

3.2.3 Instrument allocation patterns across studies

We generate a binary instrument $Z_{ij}\sim \mathrm {Bernoulli}(p_j)$ with study-specific probabilities $p_j$ . Patterns include balanced ( $p_j{=}0.50$ ), moderately unbalanced ( $0.70/0.30$ ), and a heterogeneous empirical pattern derived from treatment assignment patterns in the pooled Ebola dataset.Reference Merson and Strudwick³¹ This construction is consistent with preference-style instruments in the literature, which use either a single overall prevalence (e.g., $P(Z{=}1)\approx 0.7$ )Reference Koladjo, Escolano and Tubert-Bitter³² or cluster-specific propensities spanning a realistic range (e.g., hospital preferences between roughly 17% and 58%).Reference Ceyisakar, van, Steyerberg and Lingsma³³ In our design, $\alpha _1$ governs the instrument signal, while $\{p_j\}$ governs the overlap and information; we vary them orthogonally.

Treatment assignment: Treatment assignment follows a logistic model

$$\begin{align*}\operatorname{logit} \Pr(T_{ijk} = 1) = -0.2 + \texttt{iv\_strength}\cdot Z_{ijk} + \beta^{(T)}_{U}\, U_i + \beta_{\text{me}}\,\varepsilon^{(\text{me})}_i + u^{(T)}_{j} + v^{(T)}_{k}, \end{align*}$$

with $Z_{ijk} \sim \mathrm {Bernoulli}(p_j)$ given the instrument assignment pattern across studies, $U_i \sim \mathcal {N}(0,0.5^2)$ , $u^{(T)}_{j} \sim \mathrm {Uniform}(-0.1,0.1)$ , $v^{(T)}_{k} \sim \mathcal {N}(0,0.15^2)$ , and $\beta _{\text {me}}=0.1$ . Here, $\varepsilon ^{(\mathrm {me})}_i$ is a mean-zero individual-level noise (measurement-error) term entering the first stage via $\beta _{\mathrm {me}}\,\varepsilon ^{(\mathrm {me})}_i$ .

Outcomes: In severe infectious disease outbreaks such as Ebola, survival times are short, with discharge or death typically occurring under 30 days, so the baseline hazard operates over a narrow window. We define the cause-specific hazards for death and discharge by $h^{\text {death}}_{ijk}(t)$ and $h^{\text {disch}}_{ijk}(t)$ , respectively, where i indexes individuals, j studies, and k time. The corresponding log-hazard models are

$$\begin{align*}\begin{aligned} \log h^{\text{death}}_{ijk}(t) &= -4.2 + u^{(S)}_{j} + v^{(S)}_{k} + \beta_T\, T_{ijk} + \beta_{\text{age}}\,\text{age}_i + \beta_{\text{sex}}\, 1\{\text{Male}_i\} + \beta^{(S)}_{U}\, U_i,\\[6pt] \log h^{\text{disch}}_{ijk}(t) &= -1.1 - u^{(S)}_{j} - v^{(S)}_{k} - \beta_T\, T_{ijk} - \beta_{\text{age}}\,\text{age}_i - \beta_{\text{sex}}\, 1\{\text{Male}_i\} - \beta^{(S)}_{U}\, U_i, \end{aligned} \end{align*}$$

where $u^{(S)}_{j}\sim \mathrm {Uniform}(-0.2,0.2)$ , $v^{(S)}_{k}\sim \mathcal {N}(0,0.15^2)$ , $\beta _T=-0.4$ , $\beta _{\text {age}}=0.005$ , $\beta _{\text {sex}}=0.1$ , and $(\beta ^{(T)}_{U},\beta ^{(S)}_{U})\in \{(0.5,0.5),(1.0,1.0),(1.5,1.5),(2.0,2.0)\}$ by confounding level.

Random effects structure: Study-level random effects were drawn from Uniform distributions to induce bounded heterogeneity, while calendar-time effects were drawn from Normal distributions to allow occasional larger deviations. These choices ensure simulation stability across the factorial grid. This structure is mirrored in the analysis models, which include random intercepts for study and admission month.

Times are drawn

$$\begin{align*}T^{\text{death}}_{i} \sim \mathrm{Weibull}(1.2,\; 0.45 / h^{\text{death}}_{ijk}(t)), \quad T^{\text{disch}}_{i} \sim \mathrm{Weibull}(6.0,\; 3.0 / h^{\text{disch}}_{ijk}(t)), \end{align*}$$

then $Y_i = \min \{T^{\text {death}}_{i}, T^{\text {disch}}_{i}\}$ and $\Delta _i = 1\{T^{\text {death}}_{i} \le T^{\text {disch}}_{i}\}$ .

3.3 Estimand

The primary estimand is the cluster-conditional treatment log-hazard effect in the structural hazard model, with true value $\theta = \beta _T = -0.4$ . We evaluate whether each estimator (naïve mixed-effects Cox and mixed-effects 2SRI) recovers $\theta $ under varying design factors.

3.4 Methods

For each simulated dataset, we estimate the treatment effect using (i) a naïve mixed-effects Cox model (fit via coxme Reference Therneau³⁴) adjusting for age and sex with random effects for study and time and (ii) mixed-effects 2SRI, with a first-stage logistic mixed-effects model (fit via lme4 Reference Bates, Mächler, Bolker and Walker³⁵) adjusting for age and sex with random effects for study and time, and a second-stage Cox mixed-effects model (fit via coxme) with residual inclusion, adjusting for age and sex with random effects for study and time. Analyses are run in R (RStudio 2024.12.1+563, “Kousa Dogwood”),³⁶ with R code available on GitHub and archived on Zenodo: https://doi.org/10.5281/zenodo.20022540.

3.5 Performance measures and Monte Carlo precision

For the primary aim, we evaluate bias, coverage, empirical SE, and MSE; for the secondary aim, we examine how these same metrics vary across realized first-stage Wald $\chi ^2$ values.

Terminology: We distinguish (i) the model-based standard error $\widehat {\mathrm {SE}}_i$ returned by each fitted model in replication i (used to form Wald CIs and compute coverage), from (ii) the empirical variability across replications, $s_{\hat \theta } = \sqrt {\frac {1}{n-1}\sum _{i=1}^{n}(\hat \theta _i-\overline {\hat \theta })^2}$ , which we label “Empirical SE” in figures.

Relative bias (%):

(3.1) $$ \begin{align} \mathrm{RB} = 100 \times \frac{\left|\overline{\hat{\theta}} - \theta\right|}{|\theta|}. \end{align} $$

Coverage (95% CI):

(3.2) $$ \begin{align} \hat{p} \;=\; \frac{1}{n}\sum_{i=1}^{n} \mathbf{1}\!\left\{\theta \in \left[\hat{\theta}_i \pm 1.96\,\widehat{\mathrm{SE}}_i\right]\right\}, \end{align} $$

where $\widehat {\mathrm {SE}}_i$ is the model-based standard error from replication i.

Empirical SD of $\hat {\theta }$ (“Empirical SE” in figures):

(3.3) $$ \begin{align} \widehat{\mathrm{SD}}(\hat{\theta}) \;=\; s_{\hat{\theta}} = \sqrt{\frac{1}{n-1}\sum_{i=1}^{n}\,(\hat{\theta}_i - \overline{\hat{\theta}})^2 }. \end{align} $$

Mean squared error (MSE):

(3.4) $$ \begin{align} \mathrm{MSE} = \frac{1}{n}\sum_{i=1}^{n}\,(\hat{\theta}_i - \theta)^2. \end{align} $$

We also record the realized distribution of first-stage instrument strength (Wald statistic $W = z_{\text {IV}}^{\,2}$ ).

3.5.2 Monte Carlo precision and choice of 1,000 replications per scenario

We quantify MC uncertainty for each scenario-level summary using analytic MC standard errors (MC-SEs) that scale as $n^{-1/2}$ , where n is the number of replications.

Proportions (coverage): With $n=1,000$ and coverage near 95%, the binomial MC-SE is

$$\begin{align*}\sqrt{\tfrac{0.95\ \times\ 0.05}{1000}} \approx 0.0069, \end{align*}$$

yielding a 95% MC interval half-width of $1.96\times 0.0069 \approx 1.35$ percentage points. The worst case for proportions ( $p=0.5$ ) gives MC-SE $=0.0158$ (half-width $\approx 3.10$ percentage points).

Continuous summaries (bias): For the scenario-level bias, $\overline {\hat \theta }-\theta $ , the MC-SE is

$$\begin{align*}\mathrm{MC\text{-}SE}(\overline{\hat\theta}-\theta)=\frac{s_{\hat\theta}}{\sqrt{n}}, \end{align*}$$

where $s_{\hat \theta }$ is the across-replication standard deviation of $\hat \theta $ . In our runs $s_{\hat \theta }\in [0.15,\,0.25]$ , so at $n=1,000,$ the MC-SE is $\approx 0.005$ – $0.008$ (about 1%–2% of $|\theta |=0.4$ ).

Empirical SD and MSE: We use standard analytic formulas for the MC uncertainty of the across-replication SD and of the MSE; these appear as error bars on the forest plots.

Rationale for $n=1,000$ : MC uncertainty decreases as $1/\sqrt {n}$ . At 1,000 replications, the 95% MC interval for coverage is approximately $\pm 1.35$ percentage points; increasing to 2,000 and 4,000 would reduce this to $\sim \pm 0.95$ and $\sim \pm 0.68$ percentage points, respectively, at 2 $\times $ and 4 $\times $ the computational cost. We therefore adopt 1,000 replications per scenario as a balanced choice that yields decision-grade precision without disproportionate computational cost.