The impact of unmeasured within- and between-cluster confounding on the bias of effect estimatorsof a continuous exposure

Abstract

Unmeasured confounding almost always exists in observational studies and can bias estimates of exposure effects. Instrumental variable methods are popular choices in combating unmeasured confounding to obtain less biased effect estimates. However, we demonstrate that alternative methods may give less biased estimates depending on the nature of unmeasured confounding. Treatment preferences of clusters (e.g. physician practices) are the most frequently used instruments in instrumental variable analyses. These preference-based instrumental variable analyses are usually conducted on data clustered by region, hospital/facility, or physician, where unmeasured confounding often occurs within or between clusters. We aim to quantify the impact of unmeasured confounding on the bias of effect estimators in instrumental variable analysis, as well as several common alternative methods including ordinary least squares regression, linear mixed models, and fixed-effect models to study the effect of a continuous exposure (e.g. treatment dose) on a continuous outcome. We derive closed-form expressions of asymptotic bias of estimators from these four methods in the presence of unmeasured within- and/or between-cluster confounders. Simulations demonstrate that the asymptotic bias formulae well approximate bias in finite samples for all methods. The bias formulae show that instrumental variable analyses can provide consistent estimates when unmeasured within-cluster confounding exists, but not when between-cluster confounding exists. On the other hand, fixed-effect models and linear mixed models can provide consistent estimates when unmeasured between-cluster confounding exits, but not for within-cluster confounding. Whether instrumental variable analyses are advantageous in reducing bias over fixed-effect models and linear mixed models depends on the extent of unmeasured within-cluster confounding relative to between-cluster confounding. Furthermore, the impact of unmeasured between-cluster confounding on instrumental variable analysis estimates is larger than the impact of unmeasured within-cluster confounding on fixed-effect model and linear mixed model estimates. We illustrate the use of these methods in estimating the effect of erythropoiesis stimulating agents on hemoglobin levels. Our findings provide guidance for choosing appropriate methods to combat the dominant types of unmeasured confounders and help interpret statistical results in the context of unmeasured confounding.

Keywords

Bias formula causal inference instrumental variables linear mixed model observational study unmeasured confounders

1 Introduction

Unmeasured confounding almost always exists in observational studies and can often bias the estimates of exposure effects. For example, patients’ comorbid conditions are common confounders for treatment effects. We may know of patients’ comorbid conditions but we usually do not have detailed information on the levels of disease severity, which can lead to unmeasured confounding and subsequent bias. Observational studies are often conducted with individuals clustered by region, hospital, facility, or physician. In these studies, confounding frequently occurs within or between clusters.^1–4 In our motivating example, we study the effect of erythropoietin-stimulating agent (ESA) administration on raising hemoglobin (Hgb) levels among patients receiving hemodialysis for end-stage kidney disease from multiple dialysis facilities. Facility indicators of quality of clinical care, such as the percentage of facility patients receiving dialysis via a central venous catheter, are likely between-cluster confounders, while patients’ responsiveness to ESA (i.e. change in Hgb level) is likely a within-cluster confounder. When within- or between-cluster confounding is unmeasured or unadjusted for, most statistical methods are invalid and can give biased effect estimates. However, it is less known that some statistical methods are more robust to unmeasured confounding and give less biased effect estimates than others. In this study, we aim to quantify the bias of effect estimates obtained from the instrumental variable analysis (IVA) and several alternative methods in the presence of between- or within-cluster unmeasured confounding.

IVAs are popular choices for obtaining effect estimates robust to unmeasured confounding. Treatment preferences of clusters (e.g. physician practices) are the most frequently used instruments in the literature (i.e. preference-based instruments).^5–8 Out of the 187 comparative effectiveness research studies that used IVA between years 1990 and 2011, about one-half of the instruments were preference-based instruments.⁹ Different practices prefer different treatment dose (percentages) even after sufficient adjustments for patient/practice heterogeneity. Part of the variation in treatment dose preferences may arise from differential group policies, insurance coverage, patient/physician knowledge/preference, and may be random such that it is independent of unmeasured confounders (e.g. patients’ disease severity) after adjustment for measured confounders. The preference-based IVA aims to use the “random” component of variation to obtain treatment effect estimates. For an IV to be valid, it must be associated with the exposure of interest, independent of unmeasured confounders and must have no direct effect on the outcome, conditional on measured covariates (Figure 1). Previously, we derived a bias formula when assumptions for valid preference-based IVs are met.¹ In this study, we generalize and derive bias formulae of preference-based IV estimators in the more common and realistic scenarios when IV assumptions may or may not be satisfied in the presence of unmeasured between-cluster and/or within-cluster confounders.

It is known that IVA is advantageous in handling unmeasured confounding; however, it is less known that other popular analytic methods may give less biased effect estimates than IVA depending on the nature of unmeasured confounding. For comparisons with IVA, we choose three commonly used models: ordinary least squares (OLS) regression, linear mixed models (LMM), and fixed-effect models (FE) to estimate the effect of a continuous exposure (e.g. treatment dose) on a continuous outcome. The previous work has focused on the importance of assumptions being satisfied in order for these statistical methods to generate valid effect estimates. However, as we demonstrate in this paper, when between- and within-cluster confounding are not fully controlled, most of these methods are likely to be invalid. And the robustness of these methods toward bias differs depending on the extent of unmeasured within-cluster confounding relative to unmeasured between-cluster confounding. Our prior work derived bias formula for IVA when unmeasured within-cluster confounding exists¹; in this study, we derive bias formulae for all four methods when between-cluster and/or within-cluster confounding exits. Our study of the bias patterns and factors that impact the magnitudes of the bias for these methods will (a) provide evidence in selecting better methods to combat the dominant types of unmeasured confounders, (b) help appropriately interpret statistical results in the context of unmeasured confounding, and (c) assist in detecting the presence of unmeasured confounders.

As an illustrating example, we use these methods to estimate the effect of ESA on Hgb levels among patients receiving hemodialysis for end-stage kidney disease using data from the Dialysis Outcomes and Practice Patterns Study (DOPPS), an international prospective cohort study.¹⁰

2 Between- and within-cluster confounders

We aim to estimate the effect of an exposure (T) on outcome (Y) free of confounding. We categorize all unmeasured confounders of the T – Y relationship into between-cluster confounders (denoted by B) and within-cluster confounders (denoted by W).^1,2 Between-cluster confounders B_i are cluster-specific and identical for any patient j in each cluster i; but B_is likely differ across different clusters. Within-cluster confounders (W_ij) likely have different values for different subjects js within each cluster i; but their cluster-specific means are identical for all i. For example, when the cluster is the dialysis facility, between-cluster confounders can be for-profit status, or nurse/patient ratio of dialysis facilities. Many confounders tend to be both between-cluster and within-cluster confounders (e.g. patients’ age). For these confounders, we decompose them into between-cluster components (e.g. facility mean age) and within-cluster components (e.g. a patient’s age – facility mean age). The decomposition of unmeasured confounders helps us examine how each type affects the T – Y relationship and the assumptions required for valid statistical methods. We let C_ij denote all adjusted (measured) confounders. Here W_ij, B_i, and C_ij are K_w, K_b, and K_c dimensional vectors, respectively. We assume W_ij and B_i represent residual within-cluster and between-cluster components of the unadjusted confounders after controlling for C_ij (i.e. $(W_{i j}, B_{i}) ⊥ C_{i j}$ ).

3 Methods

Our objective is to examine the effect of a continuous exposure (i.e. treatment dose) on a continuous outcome that is free of confounding. First, we describe the preference-based IVA method and three other methods to compare their estimators and assumptions for valid effect estimates. We then assess possible violations of these assumptions when between-cluster or within-cluster unmeasured confounding exists.

3.1 Preference-based instrumental variable approach

The models for the preference-based IVA are expressed below as two simultaneous equations

T_{i j} = γ_{i} + C_{i j}^{' *} α_{I c}^{*} + e_{i j}^{t}

(1)

Y_{i j} = β_{I} T_{i j} + C_{i j}^{'} β_{I c} + v_{i} + e_{i j}^{y}

(2)

where

j = 1, 2, \dots, n_{i}

for each i and

i = 1, 2, \dots, m

, β_I is the parameter of interest denoting the effect of T on Y; γ_i represents the random exposure level at cluster i;

α_{I c}^{*} = {(α_{2 I c}, \dots, α_{K_{c} I c})}^{'}

represents the effects of

C_{i j}^{*} = {(C_{2 i j}, \dots, C_{K_{c} i j})}^{'}

on T_ij; and

β_{I c} = {(β_{1 I c}, \dots, β_{K_{c} I c})}^{'}

corresponds to the effects of

C_{i j} = {(C_{1 i j}, \dots, C_{K_{c} i j})}^{'}

on Y_ij. Note that

C_{1 i j}

denotes the intercept

(C_{1 i j} = 1)

and

C_{2 i j}, \dots, C_{K_{c} i j}

for measured confounders. Also note that different from C_ij,

C_{i j}^{*}

does not include intercept. In addition, different from LMM, γ_i can be arbitrarily correlated with

C_{i j}^{' *}

. Without loss of generality, we assume cluster size n_i = n for any i for notational simplicity in the following.

We assume the between-cluster errors (or random effects) v_i and within-cluster errors $e_{i j}^{t}$ and $e_{i j}^{y}$ are identically, independently distributed (i.i.d.) and follow normal distributions such that $v_{i} \sim N (0, σ_{v}^{2}), e_{i j}^{t} \sim N (0, σ_{e t}^{2})$ , and $e_{i j}^{y} \sim N (0, σ_{e y}^{2})$ , where $v_{i} ⊥ (e_{i j}^{t}, e_{i j}^{y})$ . Here v_i captures the intra-cluster correlation for the outcome. Let $ξ_{i} = v_{i} J_{n} + e_{i}^{y}$ where J_n is an $n \times 1$ vector of ones, and $e_{i}^{y} = {(e_{i 1}^{y}, \dots, e_{i n}^{y})}^{'}$ . Then we have $ξ_{i} \sim N (0, Ω)$ with $Ω = σ_{v}^{2} I_{n} + σ_{e y}^{2} J_{n} J_{n}^{'}$ where I_n is an identity matrix with rank n. Other assumptions for valid IVA model equations include $(v_{i}, e_{i j}^{t}, e_{i j}^{y}) ⊥ C_{i j}$ and $e_{i j}^{y} ⊥ T_{i j}$ . Note that the assumption of $e_{i j}^{t} ⊥ e_{i j}^{y}$ is not required for valid IVA.

Let $η_{I} = {(β_{I}, β_{1 I c}, \dots, β_{K_{c} I c})}^{'}$ . With preference-based IVA, the two-stage generalized least squares estimator of η_I is given by

{\hat{η}}_{I} = {(\sum_{i = 1}^{m} {\hat{O}}_{i}^{'} {\hat{Ω}}^{- 1} {\hat{O}}_{i})}^{- 1} (\sum_{i = 1}^{m} {\hat{O}}_{i}^{'} {\hat{Ω}}^{- 1} Y_{i})

(3)

where

\hat{Ω}

is an estimate of Ω and

{\hat{O}}_{i} = {({\hat{O}}_{i 1}, \dots, {\hat{O}}_{i n})}^{'}

with

{\hat{O}}_{i j} = {({\hat{T}}_{i j}, C_{1 i j}, \dots, C_{K_{c} i j})}^{'}

{\hat{T}}_{i j}

is the predicted T obtained from equation (

1

) using OLS estimation by regressing T_ij on

C i j *'

and

Z_{i}^{'}

, where Z_i is an

m \times 1

indicator vector with its elements being I

(ℓ = i) = 1

(if

ℓ = i

) or 0 (if

ℓ \neq i

) where

ℓ = 1, \dots, m

representing any potential cluster memberships. The estimator for the coefficients of Z_i is the least squares dummy variable estimator in economics^11,12,14; it is also the same estimator of the coefficients in the fixed-effect regression model in biostatistics and economics.^11–14 The estimator of β_I is given by

{\hat{β}}_{I} = (1, 0, \dots, 0) \times {\hat{η}}_{I}

. The variance of

{\hat{η}}_{I}

is estimated by

Var ({\hat{η}}_{I}) = {(\sum_{i = 1}^{m} {\hat{O}}_{i}^{'} {\hat{Ω}}^{- 1} {\hat{O}}_{i})}^{- 1} .

Hence,

Var ({\hat{β}}_{I}) = (1, 0, \dots, 0) Var ({\hat{η}}_{I}) {(1, 0, \dots, 0)}^{'}

. We obtain

\hat{Ω}

using the procedure provided in our prior work.¹

The preference-based IVA aims at taking advantage of the random component of the treatment dose assignment to obtain valid effect estimates. Conditional on measured confounders, this random component is independent of unmeasured confounders, may arise from differential cluster policies or preference, and is represented by γ_i in equation (1). The IVA assumptions require that γ_i does not have a direct effect on Y and that γ_i must be independent of the unmeasured confounders (conditional on measured covariates). This implies that γ_i needs to be independent of $v_{i}, e_{i j}^{t}$ , and $e_{i j}^{y}$ . Additional IVA assumption includes that γ is positively associated with the treatment dose received.¹

When unmeasured within-cluster confounders W_ij are present, W_ij are absorbed by $e_{i j}^{t}$ and $e_{i j}^{y}$ such that $e_{i j}^{t} = W_{i j}^{'} α_{w} + ε_{i j}^{t}$ and $e_{i j}^{y} = W_{i j}^{'} β_{w} + ε_{i j}^{y}$ with α_w and β_w representing the effects of W on T and Y, respectively. This results in Cov $(e_{i j}^{t}, e_{i j}^{y}) \neq 0$ and subsequently Cov $(T_{i j}, e_{i j}^{y}) \neq 0$ , even if $ε_{i j}^{t} ⊥ ε_{i j}^{y}$ . This non-zero correlation between T_ij and $e_{i j}^{y}$ can result in biased estimates of β_I when we only fit the single model equation (2). However, in IVA, when W_ij are present, because of $γ_{i} ⊥ W_{i j}, γ_{i} ⊥ (e_{i j}^{t}, e_{i j}^{y}, v_{i})$ remains true after adjusting for C_ij, and hence IVA assumptions are not violated. When unmeasured between-cluster confounder B_i is present, B_i is absorbed by γ_i and v_i such that $γ_{i} = B_{i}^{'} α_{b} + r_{0 i}$ and $v_{i} = B_{i}^{'} β_{b} + u_{0 i}$ with α_b and β_b representing the effect of B on T and Y, respectively. Subsequently, this leads to Cov $(γ_{i}, v_{i}) \neq 0$ , even if $r_{0 i} ⊥ u_{0 i}$ . This violates the $γ_{i} ⊥ v_{i}$ assumption required for a valid IVA and can lead to invalid IVA estimates. Here, B is also termed as IV-outcome confounders in literature.⁹

The causal interpretations of IVA estimators are described in further details by us and others.^1,15–17 Briefly, to point identify the treatment effect, we need to make the additional assumptions of either monotonicity or homogeneous effects. With the monotonicity assumption, we can interpret the IVA estimate as the causal treatment effect for the compliers, the subgroup of patients who would adopt the treatment dose suggested by the instrument. With the homogeneous effect assumption, we can interpret the IVA estimate as the average causal effect for the whole population. In this article, we assume homogeneous treatment dose effect for simulations, but make no such an assumption for the data analyses since the truth is unknown.

3.2 Alternative methods

3.2.1 Ordinary least squares regression

The model for OLS regression can be expressed as

Y_{i j} = β_{O} T_{i j} + C_{i j}^{'} β_{O c} + ς_{i j}^{y}

(4)

where β_O is the parameter of interest representing the effect of T on Y,

ς_{i j}^{y}

is the random error, and

β_{O c} = {(β_{1 O c}, \dots, β_{K_{c} O c})}^{'}

for the effects of C_ij on Y_ij. Let

η_{O} = {(β_{O}, β_{1 O c}, \dots, β_{K_{c} O c})}^{'}

O_{i j} = {(T_{i j}, C_{1 i j}, \dots, C_{K_{c} i j})}^{'}, Y_{i} = {(Y_{i 1}, Y_{i 2}, \dots, Y_{i n})}^{'}

and

T_{i} = {(T_{i 1}, T_{i 2}, \dots, T_{i n})}^{'}

. We assume that

ς_{i j}^{y} \sim N (0, σ_{ς}^{2})

. The OLS estimator (

{\hat{η}}_{O}

) of η_O is

{\hat{η}}_{O} = {(\sum_{i = 1}^{m} O_{i}^{'} O_{i})}^{- 1} (\sum_{i = 1}^{m} O_{i}^{'} Y_{i})

(5)

The OLS estimator of β_O is then given by ${\hat{β}}_{O} = (1, 0, \dots, 0) \times {\hat{η}}_{O}$ . The variance of ${\hat{η}}_{O}$ is estimated by $Var ({\hat{η}}_{O}) = {(\sum_{i = 1}^{m} O_{i}^{'} O_{i})}^{- 1} σ_{ς}^{2}$ . Hence, $Var ({\hat{β}}_{O}) = (1, 0, \dots, 0) Var ({\hat{η}}_{O}) {(1, 0, \dots, 0)}^{'}$ .

The assumptions required for the OLS regression to obtain a valid estimate of β_O are $ς_{i j}^{y}$ s are i.i.d. and $ς_{i j}^{y} ⊥ (T_{i j}, C_{kij})$ for any k. In clustered data, $ς_{i j}^{y}$ s are correlated, and hence the i.i.d. assumption is violated. Furthermore, when unmeasured between-cluster or within-cluster confounders exist, the assumption $ς_{i j}^{y} ⊥ T_{i j}$ is no longer valid.

3.2.2 Fixed-effect regression model

The fixed-effect regression model is expressed as

Y_{i j} = β_{F} T_{i j} + C_{i j}^{'} β_{F c} + μ_{2} + \dots + μ_{m} + ε_{i j}^{y}

(6)

where β_F is the parameter of interest for the effect of T on Y;

μ_{2}, \dots, μ_{m}

are fixed effects, representing cluster-specific effects for cluster 2 through m; and

ε_{i j}^{y}

is the random error with

ε_{i j}^{y} \sim N (0, σ_{ε}^{2})

. Let

{\bar{Y}}_{i} = \frac{1}{n_{i}} \sum_{j = 1}^{n_{i}} Y_{i j}, {\bar{T}}_{i} = \frac{1}{n_{i}} \sum_{j = 1}^{n_{i}} T_{i j}, {\bar{C}}_{k i} = \frac{1}{n_{i}} \sum_{j = 1}^{n_{i}} C_{kij}

, and

{\bar{ε}}_{i} = \frac{1}{n_{i}} \sum_{j = 1}^{n_{i}} ε_{i j}

. To obtain the fixed-effect estimator, we first subtract

{\bar{Y}}_{i}

from both sides of the model (6) and obtain the following model:

Y_{i j}^{†} = β_{F} T_{i j}^{†} + C_{i j}^{' †} β_{F c} + ε_{i j}^{y †}

where

Y_{i j}^{†} = Y_{i j} - {\bar{Y}}_{i}, T_{i j}^{†} = T_{i j} - {\bar{T}}_{i}, C_{kij}^{†} = C_{kij} - {\bar{C}}_{k i}

, and

ε_{i j}^{y †} = ε_{i j}^{y} - {\bar{ε}}_{i}

. Note that here

μ_{2}, \dots, μ_{m}

, the intercept

C_{1 i j}

, and any measured between-cluster confounders are eliminated by this transformation. Hence,

C_{i j}^{†}

excludes the intercept and measured between-cluster confounders. Let

O_{i j}^{†} = (T_{i j}^{†}, \dots, C_{kij}^{†}, \dots)

and

{\hat{δ}}_{F} = (β_{F}, β_{1 F c}, \dots, β_{K_{c} F c})

. The fixed-effect estimator is given by

{\hat{δ}}_{F} = {(\sum_{i = 1}^{m} O_{i}^{' †} O_{i}^{' †})}^{- 1} (\sum_{i = 1}^{m} O_{i}^{' †} Y_{i}^{' †})

The fixed-effect estimator is also named as within-cluster estimator and is numerically identical to least squares dummy variable estimator. The fixed-effect estimator of β_F is given by ${\hat{β}}_{F} = (1, 0, \dots, 0) \times {\hat{δ}}_{F}$ . The variance of ${\hat{β}}_{F}$ is estimated by $Var ({\hat{β}}_{F}) = {(\sum_{i = 1}^{m} O_{i}^{' †} O_{i}^{' †})}^{- 1} σ_{ε}^{2}$ . Hence, $Var ({\hat{β}}_{F}) = (1, 0, \dots, 0) Var ({\hat{δ}}_{F}) {(1, 0, \dots, 0)}^{'}$ .

The assumptions required for the fixed-effect regression to obtain a valid estimate of β_F are $ε_{i j}^{y}$ are i.i.d. and $ε_{i j}^{y} ⊥ (T_{i j}, C_{kij})$ for any k. Since the statistical inference for the fixed-effect regression model is built upon conditioning being in the same cluster, subjects in the same cluster are still i.i.d. even if there is a correlation among subjects who belong to the same cluster. Subsequently, unmeasured between-cluster confounders B do not have any impact on the model assumptions and are absorbed by the fixed effects $μ_{2}, \dots, μ_{m}$ . However, when unmeasured within-cluster confounders W exist, the assumption $ε_{i j}^{y} ⊥ T_{i j}$ no longer holds.

3.2.3 Linear mixed model

The LMM is written as follows

Y_{i j} = β_{L} T_{i j} + C_{i j}^{'} β_{L c} + d_{0 i} + χ_{i j}^{y}

(7)

where β_L and

β_{L c} = {(β_{1 L c}, \dots, β_{K_{c} L c})}^{'}

correspond to the effects of T_ij and

C_{i j}^{'}

, respectively. Note that

d_{0 i}

refers to the random-effect or between-cluster error, and

χ_{i}^{y} = {(χ_{i 1}^{y}, \dots, χ_{i n}^{y})}^{'}

refers to the within-cluster error.

χ_{i j}^{y}

and

d_{0 i}

represent the totality of within-cluster covariates and between-cluster covariates, respectively, omitted from the model that are orthogonal to covariates already in the model.¹⁸ The random-effect

d_{0 i}

accommodates the intra-cluster correlation. We assume that

d_{0 i}

and

χ_{i j}^{y}

are i.i.d., and

d_{0 i} \sim N (0, σ_{d}^{2})

and

χ_{i j}^{y} \sim N (0, σ_{χ}^{2})

. Let

ζ_{i} = d_{0 i} J_{n} + χ_{i}^{y}

. We have

ζ_{i} \sim N (0, Φ)

with

Φ = σ_{d}^{2} I_{n} + σ_{χ}^{2} J_{n} J_{n}^{'}

Let $η_{L} = {(β_{L}, β_{1 L c}, \dots, β_{K_{c} L c})}^{'}$ . The maximum likelihood estimator of η_L is given by

{\hat{η}}_{L} = {(\sum_{i = 1}^{m} O_{i}^{'} {\hat{Φ}}^{- 1} O_{i})}^{- 1} (\sum_{i = 1}^{m} O_{i}^{'} {\hat{Φ}}^{- 1} Y_{i})

(8)

where

\hat{Φ}

is an estimate of

Φ

. The estimator of β_L is given by

{\hat{η}}_{L} = (1, 0, \dots, 0) \times {\hat{η}}_{L}

. The variance estimate is

Var ({\hat{η}}_{L}) = {(\sum_{i = 1}^{m} O_{i}^{'} {\hat{Φ}}^{- 1} O_{i})}^{- 1}

. Hence,

Var ({\hat{β}}_{L}) = (1, 0, \dots, 0) Var ({\hat{η}}_{L}) {(1, 0, \dots, 0)}^{'}

The assumptions required for the LMM to obtain unbiased estimates of β_L include $(T_{i j}, C_{kij}) ⊥ χ_{i j}^{y}$ (namely, the level-1 independence¹⁹), $(T_{i j}, C_{kij}) ⊥ d_{0 i}$ (namely, the level-2 independence¹⁹) for any k, and $d_{0 i} ⊥ χ_{i j}^{y}$ . When unmeasured within-cluster confounding exists, it is absorbed by $χ_{i j}^{y}$ . This induces correlation between T_ij and $χ_{i j}^{y}$ and hence violates the $T_{i j} ⊥ χ_{i j}^{y}$ assumption. When unmeasured between-cluster confounding exits, it is absorbed by $d_{0 i}$ . This induces correlation between T_ij and $d_{0 i}$ and hence violates the $T_{i j} ⊥ d_{0 i}$ assumption. When both unmeasured between-cluster and within-cluster confounders exist, both level-1 and level-2 independence assumptions for valid LMM are no longer valid.

4 Bias in the presence of unmeasured confounders

In this section, we first specify true models for T and Y and then derive the expression of asymptotic bias of the four effect estimators described above when the exposure and outcome are continuous and unmeasured within-cluster and/or between-clustering confounders exist. Note that the asymptotic bias is derived assuming the number of independent units (i.e. clusters) goes to infinity ( $m \to \infty$ ). However, we will examine how well our asymptotic bias formulae approximate bias in finite samples when the number of clusters is finite through simulations in the next section.

4.1 True models

LMMs are commonly used in clustered data settings to estimate the effect of a continuous exposure (i.e. treatment dose) on a continuous outcome. We assume that T and Y are generated by LMMs with measured confounders C_ij and unmeasured confounders W_ij and B_i as follows

T_{i j} = a_{0 i} + C_{i j}^{'} α_{c} + W_{i j}^{'} α_{w} + B_{i}^{'} α_{b} + ϵ_{i j}^{t}

(9)

Y_{i j} = b_{0 i} + β T_{i j} + C_{i j}^{'} β_{c} + W_{i j}^{'} β_{w} + B_{i}^{'} β_{b} + ϵ_{i j}^{y}

(10)

where β is the parameter of interest for the effect of T on Y; α_c, α_w, α_b, β, β_c, β_w, and β_b are fixed effects;

a_{0 i}

and

b_{0 i}

are between-cluster random errors (or random effects);

ϵ_{i j}^{t}

and

ϵ_{i j}^{y}

are within-cluster random errors. We assume

a_{0 i}, b_{0 i}, ϵ_{i j}^{t}

, and

ϵ_{i j}^{y}

are i.i.d. and

a_{0 i} \sim N (0, σ_{a}^{2}), b_{0 i} \sim N (0, σ_{b}^{2}), ϵ_{i j}^{t} \sim N (0, σ_{ϵ t}^{2})

and

ϵ_{i j}^{y} \sim N (0, σ_{ϵ y}^{2})

. We also make the standard assumptions required in valid LMMs:

a_{0 i} ⊥ (C_{kij}, ϵ_{i j}^{t}, ϵ_{i j}^{y}, W_{kij}, B_{i}), b_{0 i} ⊥ (a_{0 i}, C_{kij}, T_{i j}, ϵ_{i j}^{t}, ϵ_{i j}^{y}, W_{kij}, B_{i}), ϵ_{i j}^{t} ⊥ (C_{kij}, W_{kij}, B_{i})

, and

ϵ_{i j}^{y} ⊥ (ϵ_{i j}^{t}, C_{kij}, T_{i j}, W_{kij}, B_{i})

for any k.

Let $T_{i} = {(T_{i 1}, T_{i 2}, \dots, T_{i n})}^{'}, Y_{i} = {(Y_{i 1}, Y_{i 2}, \dots, Y_{i n})}^{'}$ , $C_{i} = {(C_{i 1}, C_{i 2}, \dots, C_{i n})}^{'}, W_{i} = {(W_{i 1}, W_{i 2}, \dots, W_{i n})}^{'}, ϵ_{i}^{t} = {(ϵ_{i 1}^{t}, ϵ_{i 2}^{t}, \dots, ϵ_{i n}^{t})}^{'}$ , and $ϵ_{i}^{y} = {(ϵ_{i 1}^{y}, ϵ_{i 2}^{y}, \dots, ϵ_{i n}^{y})}^{'}$ . The true models can be expressed in a matrix form as follows

T = A_{0} + C α_{c} + W α_{w} + B α_{b} \otimes J_{n} + ϵ^{t}

(11)

Y = B_{0} + β T + C β_{c} + W β_{w} + B β_{b} \otimes J_{n} + ϵ^{y}

(12)

where Y, T, C, W,

ϵ^{t}

ϵ^{y}

, A₀, and B₀ consist of stacked elements of Y_i, T_i, C_i, W_i,

ϵ_{i}^{t}, ϵ_{i}^{y}, a_{0 i} J_{n}

, and

b_{0 i} J_{n}

, respectively. Note that Y, T,

ϵ^{t}

ϵ^{y}

, A₀ and B₀ are

m n \times 1

vectors, W and C are

m n \times K_{w}

and

m n \times K_{c}

matrices, respectively, and

B = {(B_{1}, B_{2}, \dots, B_{m})}^{'}

is an

m \times K_{b}

matrix.

4.2 Asymptotic bias

With the true models described above, we derive and compare close-form expressions of asymptotic bias of the effect estimates obtained from the four methods when their model assumptions may or may not hold in the presence of unmeasured between-cluster or within-cluster confounders. Table 1 summarizes the asymptotic bias formulae for these four methods under two scenarios: (1) when the number of clusters $m \to \infty$ and the cluster size is fixed at n; (2) when both $m \to \infty$ and the cluster size $n \to \infty$ .

Table 1.

Asymptotic bias formulae for preference-based IVA, OLS, FE, and LMM in the presence of unmeasured within- and/or between-cluster confounders.

Unmeasured	Asymptotic bias formulae
Confounder	IV	OLS	FE	LMM
	$m \to \infty$

W	$\frac{α_{w}^{'} V_{w} β_{w} / n}{σ_{a}^{2} + (α_{w}^{'} V_{w} α_{w} + σ_{ϵ t}^{2}) / n}$	$\frac{α_{w}^{'} V_{w} β_{w}}{σ_{a}^{2} + α_{w}^{'} V_{w} α_{w} + σ_{ϵ t}^{2}}$	$\frac{α_{w}^{'} V_{w} β_{w}}{α_{w}^{'} V_{w} α_{w} + σ_{ϵ t}^{2}}$	$\frac{α_{w}^{'} V_{w} β_{w}}{σ_{a}^{2} \frac{σ_{χ e}^{2}}{σ_{χ e}^{2} + (n - 1) σ_{d e}^{2}} + (α_{w}^{'} V_{w} α_{w} + σ_{ϵ t}^{2})}$
B	$\frac{α_{b}^{'} V_{b} β_{b}}{(σ_{a}^{2} + α_{b}^{'} V_{b} α_{b}) + σ_{ϵ t}^{2} / n}$	$\frac{α_{b}^{'} V_{b} β_{b}}{σ_{a}^{2} + α_{b}^{'} V_{b} α_{b} + σ_{ϵ t}^{2}}$	0	$\frac{α_{b}^{'} V_{b} β_{b}}{(σ_{a}^{2} + α_{b}^{'} V_{b} α_{b}) + σ_{ϵ t}^{2} \frac{σ_{χ e}^{2} + (n - 1) σ_{d e}^{2}}{σ_{χ e}^{2}}}$
W and B	$\frac{α_{b}^{'} V_{b} β_{b} + α_{w}^{'} V_{w} β_{w} / n}{(σ_{a}^{2} + α_{b}^{'} V_{b} α_{b}) + (α_{w}^{'} V_{w} α_{w} + σ_{ϵ t}^{2}) / n}$	$\frac{α_{b}^{'} V_{b} β_{b} + α_{w}^{'} V_{w} β_{w}}{σ_{a}^{2} + α_{b}^{'} V_{b} α_{b} + α_{w}^{'} V_{w} α_{w} + σ_{ϵ t}^{2}}$	$\frac{α_{w}^{'} V_{w} β_{w}}{α_{w}^{'} V_{w} α_{w} + σ_{ϵ t}^{2}}$	$\frac{α_{b}^{'} V_{b} β_{b} \frac{σ_{χ e}^{2}}{σ_{χ e}^{2} + (n - 1) σ_{d e}^{2}} + α_{w}^{'} V_{w} β_{w}}{(σ_{a}^{2} + α_{b}^{'} V_{b} α_{b}) \frac{σ_{χ e}^{2}}{σ_{χ e}^{2} + (n - 1) σ_{d e}^{2}} + (α_{w}^{'} V_{w} α_{w} + σ_{ϵ t}^{2})}$

	$m \to \infty$ and $n \to \infty$

W	0	$\frac{α_{w}^{'} V_{w} β_{w}}{σ_{a}^{2} + α_{w}^{'} V_{w} α_{w} + σ_{ϵ t}^{2}}$	$\frac{α_{w}^{'} V_{w} β_{w}}{α_{w}^{'} V_{w} α_{w} + σ_{ϵ t}^{2}}$	$\frac{α_{w}^{'} V_{w} β_{w}}{α_{w}^{'} V_{w} α_{w} + σ_{ϵ t}^{2}}$
B	$\frac{α_{b}^{'} V_{b} β_{b}}{σ_{a}^{2} + α_{b}^{'} V_{b} α_{b}}$	$\frac{α_{b}^{'} V_{b} β_{b}}{σ_{a}^{2} + α_{b}^{'} V_{b} α_{b} + σ_{ϵ t}^{2}}$	0	0
W and B	$\frac{α_{b}^{'} V_{b} β_{b}}{σ_{a}^{2} + α_{b}^{'} V_{b} α_{b}}$	$\frac{α_{b}^{'} V_{b} β_{b} + α_{w}^{'} V_{w} β_{w}}{σ_{a}^{2} + α_{b}^{'} V_{b} α_{b} + α_{w}^{'} V_{w} α_{w} + σ_{ϵ t}^{2}}$	$\frac{α_{w}^{'} V_{w} β_{w}}{α_{w}^{'} V_{w} α_{w} + σ_{ϵ t}^{2}}$	$\frac{α_{w}^{'} V_{w} β_{w}}{α_{w}^{'} V_{w} α_{w} + σ_{ϵ t}^{2}}$

Note: m and n refer to the number of clusters and the cluster size. The true model is $T_{i j} = a_{0 i} + C_{i j}^{'} α_{c} + W_{i j}^{'} α_{w} + B_{i}^{'} α_{b} + ϵ_{i j}^{t}$ and $Y_{i j} = b_{0 i} + β T_{i j} + C_{i j}^{'} β_{c} + W_{i j}^{'} β_{w} + B_{i}^{'} β_{b} + ϵ_{i j}^{y}$ . Note that T is for exposure, Y is for outcome, C is for measured confounders, W and B are for unmeasured within- and between-cluster confounders, respectively, $σ_{a}^{2}, σ_{b}^{2}, σ_{ϵ t}^{2}, σ_{ϵ y}^{2}$ , V_b, V_w denote the variances of $a_{0 i}, b_{0 i}, ϵ_{i j}^{t}, ϵ_{i j}^{y}$ , B_i, and W_ij for any subject j in cluster i, $σ_{χ e}^{2}$ and $σ_{d e}^{2}$ refer to asymptotic variances of the random intercept and within-cluster error, respectively, for linear mixed models.

FE: fixed-effect models; IVA: instrumental variable analysis; LMM: linear mixed models; OLS: ordinary least squares.

To obtain these bias formulae, we make two assumptions to simplify the derivation process without lessening the generality of the formulae. First, we assume that the means of W, B, Y, and T are all zeros because the means can only influence the estimates of the intercept but not the estimates for the parameter of interest, β. Second, we assume there are no measured confounders C (including the intercept). To examine the impact of C on the bias derivation, we first transform Y and T by pre-multiplying $M_{C} = I_{m n} - C {(C^{'} C)}^{- 1} C^{'}$ to equations (11) and (12) and obtain

\begin{array}{l} T^{‡} = A_{0} + W α_{w} + B α_{b} \otimes J_{n} + ϵ^{t}, \\ Y^{‡} = B_{0} + β T^{‡} + W β_{w} + B β_{b} \otimes J_{n} + ϵ^{y} \end{array}

where

T^{‡} = M_{c} T

and

Y^{‡} = M_{c} Y

are projection errors of Y and T on the space spanned by C. Note that

M_{c} ϵ^{y} = ϵ^{y}, M_{c} ϵ^{t} = ϵ^{t}

, M_cW = W, and M_cB = B because

C_{i j} ⊥ (ϵ_{i j}^{t}, ϵ_{i j}^{y}, W_{i j}, B_{i}, a_{0 i}, b_{0 i})

. The derived bias formulae only consist of the second moments of W, B,

ϵ^{t}

, and

ϵ^{y}

and will not change by the process of projection. Hence, the transformation will not change any element in the bias formulae, and the bias formulae should be the same with or without adjusting C. This technique was successfully implemented in our prior paper.¹ Through simulations, we will further confirm that the assumptions of mean zeros and no C will not influence the bias formulae.

4.2.1 Preference-based instrumental variable approach

In the absence of C, the two-stage generalized least squares estimator of β in (3) can then be simplified to

{\hat{β}}_{I} = {(\sum_{i = 1}^{m} {\hat{T}}_{i}^{'} {\hat{Ω}}^{- 1} {\hat{T}}_{i})}^{- 1} (\sum_{i = 1}^{m} {\hat{T}}_{i}^{'} {\hat{Ω}}^{- 1} Y_{i})

(13)

where

{\hat{T}}_{i} = J_{n} J_{n}^{'} T_{i} / n

and

{\hat{Ω}}^{- 1}

is an estimator for

Ω^{- 1}

. When both unmeasured between- and within-cluster confounders exist,

Ω = (β_{b}^{'} V_{b} β_{b} + σ_{b}^{2}) I_{n} + (β_{w}^{'} V_{w} β_{w} + σ_{ϵ y}^{2}) J_{n} J_{n}^{'}

and

Ω^{- 1} = H_{I} I_{n} - H_{I} H_{J} J_{n} J_{n}^{'}

, where

H_{I} = \frac{1}{β_{w}^{'} V_{w} β_{w} + σ_{ϵ y}^{2}}

H_{J} = \frac{β_{b}^{'} V_{b} β_{b} + σ_{b}^{2}}{(β_{w}^{'} V_{w} β_{w} + σ_{ϵ y}^{2}) + n (β_{b}^{'} V_{b} β_{b} + σ_{b}^{2})}

, with

σ_{b}^{2}, σ_{ϵ y}^{2}

, V_b, and V_w representing the variance of

b_{0 i}, ϵ_{i j}^{y}

, B_i, and W_ij for any i, j, respectively.

When unmeasured between- and within-cluster confounders exist, as $m \to \infty$ , as proved in Appendix 1, the bias of ${\hat{β}}_{I}$ can be approximated as

{\hat{β}}_{I} - β \to_{p} \frac{α_{b}^{'} V_{b} β_{b} + α_{w}^{'} V_{w} β_{w} / n}{(σ_{a}^{2} + α_{b}^{'} V_{b} α_{b}) + (α_{w}^{'} V_{w} α_{w} + σ_{ϵ t}^{2}) / n}

When both $m \to \infty$ and $n \to \infty, {\hat{β}}_{I} - β \to_{p} \frac{α_{b}^{'} V_{b} β_{b}}{σ_{a}^{2} + α_{b}^{'} V_{b} α_{b}}$ , a function of unmeasured between-cluster confounders but not within-cluster confounders.

When unmeasured within-cluster confounders but not between-cluster confounders exist, as $m \to \infty$ , the bias of ${\hat{β}}_{I}$ can be simplified as

{\hat{β}}_{I} - β \to_{p} \frac{α_{w}^{'} V_{w} β_{w} / n}{σ_{a}^{2} + (α_{w}^{'} V_{w} α_{w} + σ_{ϵ t}^{2}) / n}

where

σ_{a}^{2}

and

σ_{ϵ t}^{2}

denote the variance of

a_{0 i}

and

ϵ_{i j}^{t}

for any

i, j

, respectively. This bias formula was first derived in our prior paper.¹ As stated previously, when only unmeasured within-cluster confounders exist, the assumptions for the IVA are not violated and the IVA is valid. Nonetheless, finite sample bias still exits.¹ However, when both

m \to \infty

and

n \to \infty, {\hat{β}}_{I} - β \to_{p} 0

. Hence, unmeasured within-cluster confounders lead to finite-sample bias but not asymptotic bias. When the number of clusters and cluster sizes are large, the bias of instrumental variable estimates becomes negligible.

When unmeasured between-cluster but not within-cluster confounding exists, as $m \to \infty$ , the bias of ${\hat{β}}_{I}$ can be simplified as

{\hat{β}}_{I} - β \to_{p} \frac{α_{b}^{'} V_{b} β_{b}}{(σ_{a}^{2} + α_{b}^{'} V_{b} α_{b}) + σ_{ϵ t}^{2} / n}

When both $m \to \infty$ and $n \to \infty, {\hat{β}}_{I} - β \to_{p} \frac{α_{b}^{'} V_{b} β_{b}}{σ_{a}^{2} + α_{b}^{'} V_{b} α_{b}}$ , which is the same as the asymptotic bias when both unmeasured between- and within-cluster confounders exist. Hence, unmeasured between-cluster confounders result in both finite and asymptotic bias.

4.2.2 Ordinary least squares regression

In the absence of C, the OLS estimator in ( $5$ ) can be simplified to

{\hat{β}}_{O} = {(\sum_{i = 1}^{m} T_{i}^{'} T_{i})}^{- 1} (\sum_{i = 1}^{m} T_{i}^{'} Y_{i})

When unmeasured between- and within-cluster confounders exist, as shown in Appendix 1, with $m \to \infty$ , the asymptotic bias of ${\hat{β}}_{O}$ is given as

{\hat{β}}_{O} - β \to_{p} \frac{α_{b}^{'} V_{b} β_{b} + α_{w}^{'} V_{w} β_{w}}{σ_{a}^{2} + α_{b}^{'} V_{b} α_{b} + α_{w}^{'} V_{w} α_{w} + σ_{ϵ t}^{2}}

Note that this bias formula does not depend on the cluster size n. Both unmeasured between- and within-cluster confounders contribute to the bias.

When unmeasured within-cluster confounders but not between-cluster confounders exist, as $m \to \infty$ , the bias formula is simplified to

{\hat{β}}_{O} - β \to_{p} \frac{α_{w}^{'} V_{w} β_{w}}{σ_{a}^{2} + α_{w}^{'} V_{w} α_{w} + σ_{ϵ t}^{2}}

When unmeasured between-cluster but not within-cluster confounding exist, as $m \to \infty$ , the bias formula is simplified to

{\hat{β}}_{O} - β \to_{p} \frac{α_{b}^{'} V_{b} β_{b}}{σ_{a}^{2} + α_{b}^{'} V_{b} α_{b} + σ_{ϵ t}^{2}}

Note that, when only unmeasured between-cluster confounders exist, the absolute value of the asymptotic bias of ${\hat{β}}_{O}$ is smaller than that of ${\hat{β}}_{I}$ for any $n \geq 2$ and the difference between asymptotic bias of ${\hat{β}}_{O}$ and ${\hat{β}}_{I}$ becomes larger as the cluster size n grows.

4.2.3 Fixed-effect model

In the absence of C, the FE estimator in (6) can be simplified to

{\hat{β}}_{F} = {(\sum_{i = 1}^{m} T_{i}^{'} (I_{n} - Q_{n}) T_{i})}^{- 1} (\sum_{i = 1}^{M} T_{i}^{'} (I_{n} - Q_{n}) Y_{i})

where

Q_{n} = J_{n} {(J_{n}^{'} J_{n})}^{- 1} J_{n}^{'} = J_{n} J_{n}^{'} / n

as the orthogonal projection matrix of J_n.

When unmeasured between- and within-cluster confounders exist, as shown in Appendix 1, as $m \to \infty$ , the bias of ${\hat{β}}_{F}$ can be approximated as

{\hat{β}}_{F} - β \to_{p} \frac{α_{w}^{'} V_{w} β_{w}}{α_{w}^{'} V_{w} α_{w} + σ_{ϵ t}^{2}} .

Like the OLS case, this expression does not depend on the size of n and only depends on unmeasured within-cluster confounders.

When unmeasured within-cluster but not between-cluster confounders exist, the bias formula remains unchanged.

When unmeasured between-cluster but not within-cluster confounders exist, ${\hat{β}}_{F} - β \to_{p} 0$ .

This implies that only unmeasured within-cluster confounders will result in asymptotic bias in the FE estimator (i.e. as $m \to \infty$ ). It is quite intuitive since μ_i and the intercept term in equation (6) completely control for unobserved between-cluster confounders and exhaust the degree of freedom at the cluster level.

4.2.4 Linear mixed model

In the absence of C, the LMM estimator in equation (8) can be simplified as

{\hat{β}}_{L} = {(\sum_{i = 1}^{M} T_{i}^{'} {\hat{Φ}}^{- 1} T_{i})}^{- 1} (\sum_{i = 1}^{M} T_{i}^{'} {\hat{Φ}}^{- 1} Y_{i})

where

Φ

is the same as Ω in equation (13). When unmeasured between- and within-cluster confounders exist, the cluster-level error term in (7)

d_{0 i} = B_{i}^{'} β_{b} + b_{0 i}

with its variance

σ_{d}^{2} = β_{b}^{'} V_{b} β_{b} + σ_{b}^{2}

, and the individual-level error term

χ_{i j}^{y} = W_{i j}^{'} β_{w} + ϵ_{i j}^{y}

with its variance

σ_{χ}^{2} = β_{w}^{'} V_{w} β_{w} + σ_{ϵ y}^{2}

. Note that

Φ = σ_{d}^{2} I_{n} + σ_{χ}^{2} J_{n} J_{n}^{'}

. Now

d_{0 i}

is correlated with T_ij such that

Cov (d_{0 i}, T_{i j}) = β_{b}^{'} V_{b} α_{b}

, and

χ_{i j}^{y}

is also correlated with T_ij such that

Cov (χ_{i j}^{y}, T_{i j}) = α_{w}^{'} V_{w} β_{w}

. These correlations lead to the violation of LMM assumptions and subsequently we cannot obtain consistent estimates of

σ_{d}^{2}

σ_{χ}^{2}

Φ

Hence, when unmeasured between- and within-cluster confounders exist, we assume positive and bounded constants $σ_{d e}^{2}$ and $σ_{χ e}^{2}$ such that ${\hat{σ}}_{d}^{2} \to_{p} σ_{d e}^{2}$ and ${\hat{σ}}_{χ}^{2} \to_{p} σ_{χ e}^{2}$ as $m \to \infty$ for some estimators ${\hat{σ}}_{d}^{2}$ and ${\hat{σ}}_{χ}^{2}$ . Then, when $m \to \infty$ , as shown in Appendix 1, the bias of ${\hat{β}}_{L}$ can be approximated as

{\hat{β}}_{L} - β \to_{p} \frac{α_{b}^{'} V_{b} β_{b} \frac{σ_{χ e}^{2}}{σ_{χ e}^{2} + (n - 1) σ_{d e}^{2}} + α_{w}^{'} V_{w} β_{w}}{(σ_{a}^{2} + α_{b}^{'} V_{b} α_{b}) \frac{σ_{χ e}^{2}}{σ_{χ e}^{2} + (n - 1) σ_{d e}^{2}} + (α_{w}^{'} V_{w} α_{w} + σ_{ϵ t}^{2})}

When both $m \to \infty$ and $n \to \infty, {\hat{β}}_{L} - β \to_{p} \frac{α_{w}^{'} V_{w} β_{w}}{α_{w}^{'} V_{w} α_{w} + σ_{ϵ t}^{2}}$ . Note that the asymptotic bias here is the same as that of the FE estimator. This is consistent with the fact that the LMM estimator becomes the FE estimator when the cluster size increases. This is because the LMM estimator is a weighted combination of within- and between-cluster estimators, with weights depending on variance components and cluster size.²⁰

When unmeasured within-cluster confounders but not between-cluster confounders exist, the bias of ${\hat{β}}_{L}$ can be simplified as

{\hat{β}}_{L} - β \to_{p} \frac{α_{w}^{'} V_{w} β_{w}}{σ_{a}^{2} \frac{σ_{χ e}^{2}}{σ_{χ e}^{2} + (n - 1) σ_{d e}^{2}} + (α_{w}^{'} V_{w} α_{w} + σ_{ϵ t}^{2})}

When both $m \to \infty$ and $n \to \infty, {\hat{β}}_{L} - β \to_{p} \frac{α_{w}^{'} V_{w} β_{w}}{α_{w}^{'} V_{w} α_{w} + σ_{ϵ t}^{2}}$ , the same as above.

When unmeasured between-cluster but not within-cluster confounding exists, the bias of ${\hat{β}}_{L}$ can be simplified as

{\hat{β}}_{L} - β \to_{p} \frac{α_{b}^{'} V_{b} β_{b}}{(σ_{a}^{2} + α_{b}^{'} V_{b} α_{b}) + σ_{ϵ t}^{2} \frac{σ_{χ e}^{2} + (n - 1) σ_{d e}^{2}}{σ_{χ e}^{2}}}

When both $m \to \infty$ and $n \to \infty, {\hat{β}}_{L} - β \to_{p} 0$ .

5 Simulations

We conduct simulations to examine (1) how well the asymptotic bias formulae we derived approximate bias in finite samples for all four methods and (2) the bias patterns in the effect estimates from four methods when unmeasured between- and/or within-cluster confounders (B and/or W) exist. We simulate T_ij and Y_ij using the true models (9) and (10) specified below

\begin{array}{l} T_{i j} = α_{1 c} + α_{2 c} C_{2 i j} + α_{3 c} C_{3 i} + α_{1 b} B_{1 i} + α_{1 w} W_{1 i j} + a_{0 i} + ϵ_{i j}^{t}, \\ Y_{i j} = β_{1 c} + β T_{i j} + β_{2 c} C_{2 i j} + β_{3 c} C_{3 i} + β_{1 b} B_{1 i j} + β_{1 w} W_{1 i j} + b_{0 i} + ϵ_{i j}^{y} \end{array}

Previously, in order to simplify the bias derivation, we assumed mean zeros for W, B, T, and Y and no presence of C. We argued that the bias formulae should be the same with or without these assumptions. In simulations, we do not make these assumptions. The default parameter specifications are as follows: $m = 200, n = 20, α_{1 c} = 18, α_{2 c} = - 1, α_{3 c} = - 1, α_{1 b} = 0.6, α_{1 w} = 0.6$ , $β_{1 c} = 3, β = 0.7, β_{2 c} = 1, β_{3 c} = 1, β_{1 b} = 0.6, β_{1 w} = 0.6$ . We let $C_{2 i j} \sim N (0, 1), C_{3 i} \sim N (11, 1), B_{1 i} \sim N (1, 1), W_{1 i j} \sim N (1, 1), ϵ_{i j}^{t}, ϵ_{i j}^{y} \sim N (0, 1), a_{0 i} \sim N (0, {0.3}^{2})$ , and $b_{0 i} \sim N (0, 1)$ . Various simulations with a wide range of parameter values were performed; all demonstrated similar patterns. Here we present the results when we vary one parameter while holding other parameters constant. For each set of parameter specifications, we simulate 5000 and 1000 data sets for m = 10 and m = 200, respectively, since it requires more simulations to reach stable estimates when the number of clusters is small. We estimate the treatment dose effect β using the OLS, LMM, FE, and IVA methods for each data set and then report the empirical bias of the estimates averaged over simulations. We also present analytical bias calculated directly from the asymptotic bias formulae we derived for the four estimators.

5.1 When unmeasured within-cluster confounders (W) exist

When only unmeasured within-cluster confounders exist, we present the simulation results in the top panels of both Figures 2 and 3 to examine the bias of the estimates from the four methods. Figure 2 top panel shows the formulae approximate the finite-sample bias well for all methods, even when the number of clusters is very small (i.e. m = 10). When m is larger (i.e. m = 200), the approximation is even better. Overall, when only unmeasured within-cluster confounders exist, IVA estimates are much less biased than OLS, FE, and LMM estimates. As the cluster size (n) approaches 1, the bias of the IVA estimates approaches the OLS, FE, and LMM estimates; as n increases to around 400, the bias of IVA estimates becomes negligible. In contrast, n has little influence on the OLS, FE, and LMM estimates. Figure 3 top panel shows the bias patterns of the four estimators based on the bias formulae for n = 200 when either the effect of unmeasured within-cluster confounding on treatment dose ( $α_{1 w}$ ) or on outcome ( $β_{1 w}$ ) varies. When $α_{1 w}$ departs further away from 0, the magnitude of the bias for all estimates increases before decreasing. When the magnitude of $β_{1 w}$ increases, the magnitude of bias of all estimates increases monotonically. Overall, the OLS, FE, and LMM estimates are similar, and the IVA estimates are much less biased. These simulation results are consistent with what the derived asymptotic formulae indicate.

5.2 When unmeasured between-cluster confounders (B) exist

When only unmeasured between-cluster confounders B exist, Figure 2 middle panel demonstrates that the asymptotic bias formulae approximate the bias in finite samples very well for all four methods. The cluster size n has the biggest impact on IVA estimates, and some impact on LMM estimates, but very little on OLS or FE estimates. As n approaches 1, the IVA and LMM estimates approach OLS estimates. In Figure 3 bottom panel, as the effect of B on treatment $α_{1 b}$ departs further away from 0, the biases of IVA and OLS estimates first increase and then decrease. As the effect of B on outcome $β_{1 b}$ becomes more different from 0, the magnitudes of bias in both IVA and OLS estimates increase monotonically. Both Figures 2 and 3 demonstrate that when only unmeasured between-cluster confounders exist, IVA estimates have the largest bias in magnitudes, OLS estimates have some bias, and FE estimates have negligible bias. When n is relatively small (n < 20), LMM estimates have small finite-sample bias; as n increases, their bias becomes negligible. Overall, the performance of FE and LMM estimators is similar, except when n is small. These observed bias patterns are consistent with the derived bias formulae.

5.3 When both unmeasured within-cluster and between-cluster confounders (W and B) exist

When both unmeasured within-cluster and between-cluster confounders exist, we summarize the simulation results in the bottom panel of Figure 2 and in Figure 4. Figure 2 bottom panel shows that the asymptotic bias formulae approximate the finite-sample bias very well for all estimators, particularly when m = 200. Here, we assume the respective effects of W and B on T and Y are the same and set to 0.6. We find that all estimates are biased with the IVA estimates having the largest biases, followed by OLS estimates and then LMM and FE estimates. The cluster size n has the biggest impact on IVA estimates, some impact on LMM estimates, and very little impact on OLS and FE estimates. As n approaches 1, the bias of both IVA and LMM estimates approaches the bias of OLS estimates.

Figure 4 shows the bias patterns as the effects of W and B on T and Y increase based on the bias formulae. The magnitudes of the bias have complex relationships with the effects of W and B on T ( $α_{1 w}$ and $α_{1 b}$ ) for all four methods. On the other hand, the biases of all four estimators have linear relationships with both effects of W and B on Y ( $β_{1 w}$ and $β_{1 b}$ ). As $β_{1 w}$ departs further away from 0, the rate of increase in the magnitude of bias is the fastest among LMM and FE estimates, moderate among OLS estimates, and negligible among IVA estimates. As $β_{1 b}$ becomes more different from 0, the rate of increase in the magnitude of bias is the fastest among IVA estimates, moderate among OLS estimates, and negligible among FE and LMM estimates. The IVA estimates are more immune to the impact of unmeasured within-cluster confounding while the FE and LMM estimates are more robust to the impact of unmeasured between-cluster confounding. When both W and B are present, which method is the least biased depends on the interplay between W and B. When both W and B have the same effects on T and Y, IVA estimates have largest bias because IVA estimates are much more sensitive to B than what OLS, FE, and LMM estimates are to W. IVA estimates can have smallest bias among these four estimators (results not shown), only when the effect of W on Y is overwhelmingly larger than the effect of B on Y.

6 A case study

We compare these four methods in handling potential unmeasured between- and within-cluster confounders using DOPPS phase 3 data (2005–2008)^10,21 to estimate the effect of ESA on Hgb levels. The outcome of interest is the Hgb at the 14th month, and the exposure of interest is the dose of ESA measured at 1 month prior to the outcome measure. These are chosen because it typically takes about 4 weeks for ESA to be fully effective in raising Hgb levels,²² and DOPPS started to collect monthly Hgb and ESA data in study phase 3. Our analytical sample includes 1434 dialysis patients in 67 facilities. We consider the following covariates: patients’ characteristics (age, sex, race, years on dialysis, history of coronary artery disease, congestive heart failure, cancer, cerebrovascular disease, diabetes, gastrointestinal bleeding, peripheral vascular disease, hypertension, intravenous iron, psychiatric disorder, intravenous iron use) as well as facility quality indicators (i.e. percentage of patients using a central venous catheter for dialysis, serum albumin level < 3.5 g/dL, serum phosphorus level >5.5 mg/dL, and single pool kt/V (a measure of dialysis dose) < 1.2). We also consider patients’ ESA responsiveness, an indicator of ESA dose required to raise or sustain Hgb level, and effectively captured by the combination of ESA dose at 2 months prior to the outcome Hgb and Hgb level at 1 month prior to the outcome Hgb (as demonstrated in our prior work¹). Prior work by others has often omitted these adjustments.²³ We conduct three sets of analyses. The first set of analyses includes all variables, the second set excludes patients’ ESA responsiveness variables, and the third excludes ESA responsiveness variables and facility quality indicators. The results of three sets of analyses are summarized in Table 2. For the IVA, we choose the ESA dose preferences of dialysis facilities as the instrument. An instrument is called strong when the association between this instrument and the exposure of interest is strong.^24,25 The stronger the IV is, the smaller the finite sample bias is. The partial F statistic in the regression of the exposure on the instrument is usually used as a measure of the strength of an instrument.²⁵ An instrument is considered weak if the F statistic is <10.²⁵ The F statistics for IVA in these three sets of analyses are 5, 6, and 6, respectively, indicating that a stronger instrument would be preferred.^24–26

Table 2.

Three sets of analyses to estimate the effect of ESA on Hgb.

	All Vars		–ESA Resp		–ESA Resp/Fac Vars
Methods	Est	95% CI	Est	95% CI	Est	95% CI
IVA	0.039	(0.015, 0.063)	0.047	(0.016, 0.078)	0.10	(0.072, 0.127)
OLS	0.018	(0.012, 0.024)	–0.003	(–0.010, 0.004)	0.011	(0.004, 0.018)
FE	0.015	(0.009, 0.022)	–0.011	(–0.018, –0.005)	–0.011	(–0.018, –0.005)
LMM	0.017	(0.010, 0.023)	–0.010	(–0.016, –0.003)	–0.008	(–0.015, –0.002)

Note that ESA is known to increase Hgb levels and the effect of ESA on Hgb should be positive. All Vars: adjust for all covariates listed in text; ESA Resp: adjust for all but patient ESA responsiveness variables; ESA Resp/Fac Vars: adjust for all but patient ESA responsiveness variables and facility quality indicators.

ESA: erythropoietin-stimulating agent; FE: fixed-effect models; Hgb: hemoglobin; IVA: preference-based instrumental variable analysis; LMM: linear mixed models; OLS: ordinary least squares.

In the first set of analyses adjusting for all variables, although the magnitudes of the estimated ESA effects differ across the four methods, all effects are positive and qualitatively consistent with the well-known fact that ESA increases Hgb levels.²² In the second set of analyses excluding ESA responsiveness variables, IVA continues to produce positive ESA effects, but the other three methods give falsely negative estimates. In the third set of analyses excluding both ESA responsiveness and facility quality indicators, the largest change in the estimated effect of ESA on Hgb, compared with the second set, is observed in the IVA method, while there is some change in the OLS estimate, a very small change in the LMM estimate, and no change in the FE estimate.

As demonstrated in our prior work,¹ the ESA responsiveness variables are largely within-cluster confounders of the ESA–Hgb relationship. Hence, without adjusting for these variables, it is not surprising that the IVA estimate of ESA effect on Hgb does not change much and remains positive, given the robustness of the IVA method to unmeasured within-cluster confounders. It is also not surprising that all other three methods are noticeably impacted and give falsely negative effect estimates, given their sensitivity toward unmeasured within-cluster confounders. When we further exclude facility quality indicators, which are between-cluster confounders, the biggest change occurs in the IVA estimate, with little or no change in the LMM and FE estimates. Again, this is consistent with our finding that LMM and FE are robust to unmeasured between-cluster confounding, but IVA is not.

7 Discussion

IVA is known for its advantages in combating unmeasured confounders. However, it is less known that certain commonly used alternative methods may be more robust to unmeasured confounding than IVA depending on the nature of the confounding. In this study, we focused on preference-based IVA (the most popular IV methods) with clustered data (the most common data structure for observational studies). We derived asymptotic bias for IVA and three alternative methods including OLS, FE, and LMM, when unmeasured between-cluster and/or within-cluster confounders exist. Simulations demonstrated that all bias formulae perform well in finite samples. We also examined the validity of the assumptions required for each method and the degree of bias when either type of unmeasured confounding exists. While almost none of the methods are valid when either type of unmeasured confounding exists, we found that some methods are more robust than others depending on the type of unmeasured confounding. When unmeasured within-cluster confounding is more probable, IVAs are least biased compared with other methods; on the other hand, when unmeasured between-cluster confounding is more probable, FE methods or LMMs are less biased than other methods. When both types of unmeasured confounding exist, it depends on which unmeasured confounders are dominant. Our results provide guidance for researchers because the effect of each type of unmeasured confounding is specified directly with explicit and closed-form bias formulae.

Our analysis of bias gives insight into the robustness and assumption of each method in handling unmeasured confounders. In practice, it is more likely that both unmeasured between- and within-cluster confounders exist and that subsequently none of the methods have all their model assumptions met. When all methods give similar analysis results, it is likely that these results are robust toward unmeasured confounding. When effect estimates from these four methods differ substantially, as in our data analyses, it is likely that unmeasured confounders exist. We can examine which potential unmeasured within- or between-cluster confounders may exist with our substantive knowledge about the data and the research questions. If we believe that unmeasured within-cluster confounding strongly dominates over unmeasured between-cluster confounding, IVA estimates should be closer to the truth than the other methods. On the other hand, if we believe that unmeasured between-cluster confounding dominates within-cluster confounding, LMM or FE estimates should be closer to the truth. IVA methods are more sensitive to unmeasured between-cluster confounding than FE or LMMs toward unmeasured within-cluster confounding. Hence, our research results are most useful when we have substantive knowledge about whether between- or within-cluster confounding is of major concern. For example, in our data, ESA responsiveness is known to be a dominating within-cluster (patient-level) confounder of the ESA–Hgb relationship and is of a major concern if not adjusted for. In summary, our analysis of bias provides insights into the potential impact of unmeasured confounders and how we can interpret our analysis results accordingly. It also helps us become more aware of the strength and weakness of the four different methods. It gives us the opportunities to detect the presence of unmeasured confounding and subsequently adjust for more confounding. We can then translate our finding to broader audience better and more precisely.

Previous research has demonstrated that the effect of treatment dose itself can be partitioned into between- and within-cluster components.² In this article, we assume the between- and within-cluster treatment effects are the same. If the between-cluster level treatment effect is different from the within-cluster level treatment effect, it often implies there may be unmeasured confounders, selection bias, or measurement errors.²⁷ All of these can be characterized as omitted covariate problems either at within-cluster or between-cluster level. Without any omitted covariates, it often requires the between-cluster and within-cluster effects to be equal for the effects to have causal interpretations.²⁷ It is worthwhile to point out that there is a connection between the treatment effect estimates from the preference-based IVA and the between-cluster treatment effect estimates as well as a connection between the treatment effect estimates from the fixed-effect models and the within-cluster treatment effect estimates. These connections are worthy of further investigations.

Previous studies focused on the impact of between-cluster unmeasured confounding on mixed models^2,27–29 or the validity of assumptions and finite sample bias for IVA.^1,24,30 To our knowledge, this study is the first to derive bias formulae for these four methods when unmeasured within-cluster and/or between-cluster confounding may exist, a common scenario in medical studies using clustered data sets. Our data generating models are LMMs, which are most commonly used in clustered data structures. We also assumed that unmeasured confounders are normally distributed and independent of other confounders in deriving the bias. In practice, it is likely that multiple unmeasured between- and within-cluster confounders exist, and many of them are correlated with other confounders. We simulated unmeasured confounding to represent the overall residual confounding after conditioning on other measured confounders and used the normal distribution to approximate the residual confounders, as often done.³¹ Our results are general and form a foundation to perform sensitivity analyses. Our research focuses on continuous outcomes, continuous exposures, and linear association between them. Our finding that between- and within-cluster unmeasured confounders have different impacts on the bias of the effect estimators from these four types of models likely extends to other types of outcomes and exposures. It will be important to investigate the bias expression for other types of outcomes, exposures, or nonlinear associations between them, where the nonlinearity could incur additional bias.³²

Our research quantifies the bias in the presence of unmeasured confounders. The results are also applicable to measurement error, selection bias, or selective drop-out issues because they can be cast into the same general framework of omitted variables. We assume the treatment effect is homogeneous, as done in most regression analyses. However, this assumption may not hold. When the effect is heterogeneous, IVA measures the treatment effect among compilers, that is, complier average treatment effect,¹⁷ which may also explain why the IVA effect can be different from other model estimates. However, when there is no reason to believe the effect may be different for compliers, we can interpret treatment effect estimates in the context of unmeasured confounding.

In summary, we derived the bias formulae for IVA, LMM, FE, and OLS in the presence of unmeasured between- and within-cluster confounding. Our findings provide evidence to support future selection of methods to combat the dominant types of unmeasured confounders, and facilitate the interpretation of statistical analysis results in the context of unmeasured confounding, and help detect the presence of unmeasured confounders.

Supplemental Material

SMM883323 Supplemetal Material - Supplemental material for The impact of unmeasured within- and between-cluster confounding on the bias of effect estimatorsof a continuous exposure

Supplemental material, SMM883323 Supplemetal Material for The impact of unmeasured within- and between-cluster confounding on the bias of effect estimatorsof a continuous exposure by Yun Li, Yoonseok Lee, Friedrich K Port and Bruce M Robinson in Statistical Methods in Medical Research

Footnotes

Acknowledgements

The authors thank Drs Min Zhang, Lu Wang, and the reviewers for helpful suggestions, as well as Shauna Leighton for providing editorial assistance on this article.

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported in part by National Institutes of Health grants 5R01-DK070869 and UL1TR002240. The DOPPS is administered by Arbor Research Collaborative for Health and is supported by Amgen, Kyowa Hakko Kirin, AbbVie, Sanofi Renal, Baxter Healthcare, and Vifor Fresenius Medical Care Renal Pharma. Additional support for specific projects and countries is provided by Keryx Biopharmaceuticals, Merck Sharp & Dohme, Proteon Therapeutics, Relypsa, and F Hoffmann-LaRoche; in Canada by Amgen, BHC Medical, Janssen, Takeda, and Kidney Foundation of Canada (for logistics support); in Germany by Hexal, DGfN, Shire, and WiNe Institute; and for PDOPPS in Japan by the Japanese Society for Peritoneal Dialysis. All support is provided without restrictions on publications.

Appendix 1

References

Lee

Wolfe

, et al. On a preference-based instrumental variable approach in reducing unmeasured confounding-by-indication. Stat Med 2015; 34(7): 1150–1168.

Neuhaus

Kalbfleisch

JD.

Between- and within-cluster covariate effects in the analysis of clustered data.

Biometrics 1998; 54: 638–645.

Salas

Hofman

Stricker

BH.

Confounding by indication: An example of variation in the use of epidemiologic terminology.

Am J Epidemiol 1999; 149: 981–983.

Bosco

Silliman

Thwin

, et al. A most stubborn bias: No adjustment method fully resolves confounding by indication in observational studies. J Clin Epidemiol 2010; 63: 64–74.

Korn

Baumrind

Clinician preferences and the estimation of causal treatment differences.

Stat Sci 1998; 13(3): 209–235.

Chen

Briesacher

BA.

Use of instrumental variable in prescription drug research with observational data: A systematic review.

J Clin Epidemiol 2011; 64(6): 687–700.

Davies

Smith

Windmeijer

, et al. Issues in the reporting and conduct of instrumental variable studies: A systematic review. Epidemiology 2013; 24: 363–369.

Brookhart

Wang

Solomon

, et al. Evaluating short-term drug effects using a physician-specific prescribing preference as an instrumental variable. Epidemiology 2006; 17(3): 268–275.

Garabedian

Chu

Toh

, et al. Potential bias of instrumental variable analyses for observational comparative effectiveness research potential bias of instrumental variable analyses for observational research. Ann Intern Med 2004; 161: 131–138.

10.

Pisoni

Gillespie

Dickinson

, et al. The Dialysis Outcomes and Practice Patterns Study (DOPPS): Design, data elements, and methodology. Am J Kidney Dis 2004; 44: 7–15.

11.

Baltagi

BH.

Econometric analysis of panel data. 3rd ed. New York: John Wiley & Sons Inc., 2005.

12.

Hsiao

Analysis of panel data. Cambridge: Cambridge University Press, 2014.

13.

Allison

PD.

Fixed effects regression models. Thousand Oaks: Sage Publications, 2009.

14.

Wooldridge

JM.

Econometric analysis of cross section and panel data. Cambridge: MIT Press, 2002.

15.

Baiocchi

Cheng

Small

DS.

Instrumental variable methods for causal inference.

Stat Med 2014; 33: 2297–2340.

16.

Sargan

JD.

The estimation of economic relationships using instrumental variables.

Econometrica 1958; 26(3): 393–415.

17.

Imbens

Angrist

Identification and estimation of local average treatment effects.

Econometrica 1994; 62(2): 467–475.

18.

Neuhaus

Kalbfleisch

Hauck

WW.

A comparison of cluster-specific and population-averaged approaches for analyzing correlated binary data.

Int Stat Rev 1991; 59: 25–35.

19.

Ebbes

Bockenholt

Wedel

Regressor and random-effects dependencies in multilevel models. Statistica Neerlandica 2004; 58: 161–178. 10.1046/j.0039-0402.2003.00254.x.

20.

Fitzmaurice

Laird

Ware

JH.

Applied longitudinal analysis. 2nd ed. Hoboken: Wiley, 2011.

21.

Young EW, Goodkin DA, Mapes DL, et al. The Dialysis Outcomes and Practice Patterns Study (DOPPS): An international hemodialysis study. Kidney Int 2000; 57(suppl 74): S-74--S-81.

22.

Eschbach

Egrie

Downing

, et al.

Correction of the anemia of end-stage renal disease with recombinant human erythropoietin. Results of a combined phase I and II clinical trial.

New Engl J Med 1987; 316(2): 73–78.

23.

Madore

Lowrie

Brugnara

, et al. Anemia in hemodialysis patients: Variables affecting this outcome predictor. J Am Soc Nephrol 1997; 8(12): 1921–1929.

24.

Burgess

Thompson

and CRP CHD Genetics Collaboration.

Avoiding bias from weak instruments in Mendelian randomization studies.

Int J Epidemiol 2011; 40(3): 755–764.

25.

Stock

Wright

Yogo

A survey of weak instruments and weak identification in generalized method of moments.

J Am Stat Assoc 2002; 20(4): 518–529.

26.

Staiger

Stock

Instrumental variables regression with weak instruments.

Econometrica 1997; 65(3): 557–586.

27.

Palta

Seplaki

. Causes, problems and benefits of different between and within effects in the analysis of clustered data.

Health Serv Outcomes Res Methodol 2002; 3: 177.

28.

Chao

Palta

Young

Effect of omitted confounders on the analysis of correlated binary data.

Biometrics 1997; 53(2): 678–689.

29.

Palta

Yao

TJ.

Analysis of longitudinal data with unmeasured confounders.

Biometrics 1991; 47(4): 1355–1369.

30.

Martens

Pestman

de Boer

, et al. Instrumental variables: Application and limitations. Epidemiology 2006; 17(3): 260–267.

31.

Cessie

SL.

Bias formulas for estimating direct and indirect effects when unmeasured confounding is present. Epidemiology 2016; 27(1): 125–132. doi: 10.1097/EDE.0000000000000407.

32.

Lee

Phillips

PCB.

Model selection in the presence of incidental parameters.

J Econometrics 2015; 188(2): 474–489.

Supplementary Material

Please find the following supplemental material available below.

For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.

For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.

0.00 MB

0.48 MB