Causal variance decompositions for institutional comparisons in healthcare

TL;DR

This paper introduces a causal variance decomposition method to analyze and compare healthcare institutions based on disease-specific quality indicators, focusing on overall variation and causal effects rather than pairwise hospital comparisons.

Contribution

It develops a novel three-way variance decomposition framework with causal interpretation, including model-based estimators for analyzing hospital performance variations.

Findings

Effective in decomposing variation into causal components

Model estimators perform well in simulations

Applied successfully to real healthcare data

Abstract

There is increasing interest in comparing institutions delivering healthcare in terms of disease-specific quality indicators (QIs) that capture processes or outcomes showing variations in the care provided. Such comparisons can be framed in terms of causal models, where adjusting for patient case-mix is analogous to controlling for confounding, and exposure is being treated in a given hospital, for instance. Our goal here is to help identifying good QIs rather than comparing hospitals in terms of an already chosen QI, and so we focus on the presence and magnitude of overall variation in care between the hospitals rather than the pairwise differences between any two hospitals. We consider how the observed variation in care received at patient level can be decomposed into that causally explained by the hospital performance adjusting for the case-mix, the case-mix itself, and residual variation. For this purpose, we derive a three-way variance decomposition, with particular attention to its causal interpretation in terms of potential outcome variables. We propose model-based estimators for the decomposition, accommodating different link functions and either fixed or random effect models. We evaluate their performance in a simulation study and demonstrate their use in a real data application.

Tables

Table 1. Table 1: A summary table of the effects of case-mix factors for the four different QIs. The numbers are log-odds ratios (z-scores), while m 𝑚 m and n 𝑛 n represent the number of hospitals and patients. dx stands for diagnosis and tx for treatment (surgery).

Process/outcome
Covariate	Partial	Partial-CKD	MIS	Readmission
	$n=4738$	$n=3028$	$n=4100$	$n=7362$
	$m=112$	$m=110$	$m=72$	$m=115$
	${\chi }^2=377$	${\chi }^2=216$	${\chi }^2=825$	${\chi }^2=16$
Male sex	0.21 (2.81)	0.22 (2.39)	0.05 (0.67)	-0.02 (-0.18)
Age	-0.01 (-4.11)	-0.01 (-2.67)	-0.01 (-2.06)	0.01 (2.74)
Income quintile	0.03 (1.14)	0.03 (1.07)	0.01 (0.47)	<|0.01|(-0.05)
Rural vs urban	0.18 (1.46)	0.20 (1.37)	-0.14 (-1.17)	-0.18 (-1.33)
Charlson score	-0.14 (-4.93)	-0.16 (-5.02)	-0.11 (-4.37)	0.12 (6.56)
ACG score	-0.01 (-1.20)	<|0.01|(-0.19)	0.01 (1.61)	0.02 (4.64)
log(dx to tx days+1)	0.02 (1.02)	0.03 (1.41)	0.02 (0.65)	0.09 (3.64)
Year of dx	0.19 (21.87)	0.18 (17.13)	0.20 (15.28)	-0.01 (-1.43)
Tumor size (cm)	-0.82 (-17.25)	-0.75 (-12.87)	-0.16 (-7.21)	0.02 (1.44)
CKD risk factors	-0.02 (-0.26)	-	-0.03 (-0.36)	0.15 (1.31)
T2 vs T1 stage	-	-	-0.31 (-2.16)	-0.03 (-0.21)
T3 or 4 vs T1 stage	-	-	-	0.19 (-1.60)

Table 2. Table 2: Three-way variance decomposition for the observed variance in the four different QIs. The numbers are absolute and proportional variances (95% CI = 95% credible interval).

Source of variation
Process/outcome		Patient case-mix	Between-hospital	Residual
Partial	Variance	$0.067$	$0.017$	$0.164$
	$95\%$ CI	$(0.062,0.072)$	$(0.014,0.020)$	$(0.158,0.169)$
	Proportion	$26.9\%$	$6.9\%$	$66.2\%$
	$95\%$ CI	$(24.9\%,29.0\%)$	$(5.5\%,8.1\%)$	$(64.0\%,68.5\%)$
Partial-CKD	Variance	$0.058$	$0.016$	$0.173$
	$95\%$ CI	$(0.052,0.065)$	$(0.012,0.020)$	$(0.165,0.179)$
	Proportion	$23.7\%$	$6.3\%$	$70.0\%$
	$95\%$ CI	$(21.2\%,26.3\%)$	$(4.7\%,8.1\%)$	$(67.1\%,72.7\%)$
MIS	Variance	$0.023$	$0.044$	$0.164$
	$95\%$ CI	$(0.019,0.027)$	$(0.039,0.050)$	$(0.158,0.170)$
	Proportion	$9.9\%$	$19.0\%$	$71.1\%$
	$95\%$ CI	$(8.1\%,11.5\%)$	$(16.8\%,21.6\%)$	$(68.4\%,73.7\%)$
Readmission	Variance	$0.0020$	$0.0001$	$0.0679$
	$95\%$ CI	$(0.0013,0.0027)$		$(0.0670,0.0685)$
	Proportion	$2.81\%$	$0.16\%$	$97.03\%$
	$95\%$ CI	$(1.92\%,3.92\%)$		$(95.74\%,97.94\%)$

Equations

$$E[Y_i\mid Z_i=z,X_i]=g^{-1}\left({\alpha }_0+{\alpha }_z+{\beta }^{\prime }{X}_i\right),$$

$$V[Y(Z)]=V_X[E(Y(Z)\mid X)]+E_X[V(Y(Z)\mid X)].$$

$$E(Y(Z)\mid X)=E_{Z\mid X}[E(Y(Z)\mid Z,X)],$$

$$V(Y(Z)\mid X)=V_{Z\mid X}[E(Y(Z)\mid Z,X)]+E_{Z\mid X}[V(Y(Z)\mid Z,X)].$$

$$\begin{array}{rl}V[Y(Z)]=& V_X\{E_{Z\mid X}[E(Y(Z)\mid Z,X)]\}+E_X\{V_{Z\mid X}[E(Y(Z)\mid Z,X)]\}\\ & +E_{Z,X}[V(Y(Z)\mid Z,X)].\end{array}$$

$$E(Y(z)\mid Z=z,X)=E(Y(z)\mid X)\text{ and }V(Y(z)\mid Z=z,X)=V(Y(z)\mid X).$$

$$V[Y]$$

$$+E_X\left\{\sum _z{\left[E(Y(z)\mid X)-\sum _{z^{\prime }}E(Y(z^{\prime })\mid X)P(Z=z^{\prime }\mid X)\right]}^{2}P(Z=z\mid X)\right\}$$

$$+E_X\left\{\sum _{z}V(Y(z)\mid X)P(Z=z\mid X)\right\}.$$

$$V[Y]=\text{total observed variance in care received}$$

$$=\text{variance explained by the patient case-mix}$$

$$+\text{average variance causally explained by the between-hospital differences in performance}$$

$$+\text{unexplained (residual) variance}.$$

$$V[Y]$$

$$V_X\left\{\sum _{z}E(Y(z)\mid X)P(Z=z\mid X)\right\}$$

$$\begin{array}{rl}& E_X\left\{\sum _{z=1}^2{\left[E(Y(z)\mid X)-\sum _{z^{\prime }=1}^{2}E(Y(z^{\prime })\mid X)P(Z=z^{\prime }\mid X)\right]}^{2}P(Z=z\mid X)\right\}\\ & =E_X\left\{\pi (X)(1-\pi (X)){\left[E(Y(1)\mid X)-E(Y(2)\mid X)\right]}^2\right\}\\ & =E_X\left\{V(Z\mid X){\left[E(Y(1)\mid X)-E(Y(2)\mid X)\right]}^2\right\}.\end{array}$$

$$E(Y(a)\mid A=a,X)=E(Y(a)\mid X)=E(Y(a)\mid Z=a,X),$$

$$V[Y(A)]$$

$$+E_X\left\{\sum _a{\left[E(Y(a)\mid X)-\sum _{a^{\prime }}E(Y_i(a^{\prime })\mid X)P(A=a^{\prime }\mid X)\right]}^{2}P(A=a\mid X)\right\}$$

$$+E_X\left\{\sum _{a}V(Y(a)\mid X)P(A=a\mid X)\right\}.$$

$$E(Y(a)\mid Z=a,X)=E(Y(a)\mid X)={\alpha }_0+{\alpha }_a+{\beta }^{\prime }X,$$

$$E_X\left\{\sum _a{\left[E(Y(a)\mid X)-\sum _{a^{\prime }}E(Y(a^{\prime })\mid X)P(A=a^{\prime }\mid X)\right]}^{2}P(A=a\mid X)\right\}$$

$$=E_X\left\{\frac{1}{m}\sum _a{\left[E(Y(a)\mid X)-\frac{1}{m}\sum _{a^{\prime }}E(Y(a^{\prime })\mid X)\right]}^2\right\}$$

$$=\frac{1}{m}\sum _a{\left[{\alpha }_a-\frac{1}{m}\sum _{a^{\prime }}{\alpha }_{a^{\prime }}\right]}^2.$$

$$V(Y(A)\mid X)$$

$$\to {\tau }^2+{\sigma }^2\text{when}m\to \infty ,$$

$$E[Y(z)\mid X;\theta ]=E[Y\mid Z=z,X;\theta ]=g^{-1}\left({\alpha }_0+{\alpha }_z+{\beta }^{\prime }X\right),$$

$${\omega }_1(\theta ,\eta )$$

$$\approx \frac{1}{n-1}\sum _{i=1}^n{\left\{\sum _z{\mu }_i(z;\hat{\theta })P(Z_i=z\mid x_i;\hat{\eta })-\frac{1}{n}\sum _{i^{\prime }=1}^n\sum _z{\mu }_{i^{\prime }}(z,\hat{\theta })P(Z_{i^{\prime }}=z\mid x_{i^{\prime }};\hat{\eta })\right\}}^2,$$

$$P(Z=z\mid X;\eta )=\{\begin{array}{ll}\frac{1}{1+{\sum }_{z^{\prime }=2}^m\exp ({\gamma }_{z^{\prime }}+{\varphi }_{z^{\prime }}^{\prime }X)}& z=1\\ \frac{\exp ({\gamma }_z+{\varphi }_z^{\prime }X)}{1+{\sum }_{z^{\prime }=2}^m\exp ({\gamma }_{z^{\prime }}+{\varphi }_{z^{\prime }}^{\prime }X)}& z\ne 1,\end{array}$$

$${\omega }_2(\theta ,\eta )$$

$$\displaystyle\approx\frac{1}{n}\sum_{i=1}^{n}\left\{\sum_{z}\mu_{i}(z;\hat{\theta})^{2}P(Z_{i}=z\mid x_{i};\hat{\eta})-\Bigg{[}\sum_{z}\mu_{i}(z;\hat{\theta})P(Z_{i}=z\mid x_{i};\hat{\eta})\Bigg{]}^{2}\right\}.$$

$${\omega }_3(\theta ,\eta )$$

$$\approx \frac{1}{n}\sum _{i=1}^n\left\{\sum _z{\mu }_i(z;\hat{\theta })[1-{\mu }_i(z;\hat{\theta })]P(Z=z\mid x_i;\hat{\eta })\right\}.$$

$$p(\theta ,\eta \mid \mathbf{X},\mathbf{Y},\mathbf{Z})$$

$$=\frac{p(\mathbf{Y}\mid \theta ,\mathbf{X},\mathbf{Z})p(\theta \mid \mathbf{X})}{p(\mathbf{Y}\mid \mathbf{X},\mathbf{Z})}\frac{p(\mathbf{Z}\mid \eta ,\mathbf{X})p(\eta \mid \mathbf{X})}{p(\mathbf{Z}\mid \mathbf{X})}$$

$$=p(\theta \mid \mathbf{X},\mathbf{Y},\mathbf{Z})p(\eta \mid \mathbf{X},\mathbf{Z}),$$

$$g(E[Y(z)\mid X;\theta ])={\alpha }_z+X_1+2X_2,$$

Full text

Keywords: Causal inference, hospital profiling, intraclass correlation, quality indicators, variance decomposition

1 Introduction

1.1 Background

Striving for accountability and transparency in healthcare has increased the interest in quantifying and comparing the performance of institutions in terms of the quality of care that the patients receive. This is commonly done in terms of disease-specific quality indicators (QIs) that measure structural, process or outcome elements related to the care of a particular condition [1]. For instance, in the context of surgical care for kidney cancer, examples would be (i) proportion of partial versus radical nephrectomies for early-stage patients with chronic kidney disease risk factors, (ii) proportion of minimally invasive versus open radical nephrectomies for early-stage patients, (iii) average length of stay after a radical nephrectomy or (iv) proportion of unplanned readmissions within 30 days of a radical nephrectomy [2]. Process measures may be preferred for measuring quality for the reason that they are directly related to clinical decisions and thus are more actionable [3]. However, even before considering improving quality, the chosen indicator has to demonstrate existing differences in the quality of care. In doing so, it already answers a causal question, namely “would a patient receive different care if treated in a different hospital?”. The goal of the present work is developing metrics and statistical methods that can help to identify processes or outcomes that can capture variation in the quality of care between hospitals. Following Varewyck et al. [4], this can be framed in a causal inference framework, where adjusting for patient case-mix factors is analogous to controlling for confounding. However, since our focus is in choosing QIs, rather than comparing or ranking hospitals in terms of an already chosen metric, we are interested in the presence and magnitude of overall variation between the hospitals, rather than contrasts between any two hospitals. While we restrict the discussion to between hospital comparisons, without loss of generality the methods would also apply to other comparisons in the healthcare system, such as administrative subregions or individual providers (such as surgeons).

Variance decompositions have been used in the hospital profiling context to measure reliability of rankings in terms of a given quality indicator (rankability)[5, 6, 7], though this decomposition is constructed at aggregate (hospital) level rather than at individual patient level. This approach broadly corresponds to the random effect meta-analysis model, where the case-mix adjusted hospital-specific quality indicators $\hat{\theta}_{z}$ are assumed to follow the model $\hat{\theta}_{z}=\theta_{0}+\theta_{z}+\xi_{z}$, where $z=1,\ldots,m$ is used to index the hospital and $m$ stands for the total number of hospitals. In this model, $\theta_{z}\sim N(0,\tau^{2})$ is a normally distributed hospital effect (‘heterogeneity’) around the average performance $\theta_{0}$, and $\xi_{z}\sim N(0,s_{z}^{2})$ reflects the sampling variation under the ‘null’ of no hospital effect ($\tau^{2}=0$). For instance, following van Houwelingen and Brand [5] and Racz and Sedransk [8], to obtain an indirectly standardized proportion type quality indicator, one could take $\hat{\theta}_{z}=p\frac{O_{z}}{E_{z}}$, where the overall proportion $p$ in the entire patient population is multiplied by the standardized ratio of observed to expected counts. The observed count for hospital $z$ is given by $O_{z}=\sum_{i=1}^{n}\mathbf{1}_{\{Z_{i}=z\}}Y_{i}$, where $Y_{i}$ is the process/outcome measure of interest, and $Z_{i}$ indicates the hospital where patient $i\in\{1,\ldots,n\}$ was treated, with $n$ being the size of the overall patient population serving as the standard population (e.g. combination of patient populations from hospitals in a given administrative region). The expected count is given by $E_{z}=\sum_{i=1}^{n}\mathbf{1}_{\{Z_{i}=z\}}\hat{\pi}_{i}$, where $\hat{\pi}_{i}=\textrm{expit}(\hat{\alpha}_{0}+\hat{\beta}^{\prime}X_{i})$ is a prediction from a logistic regression model adjusted for the case-mix factors $X_{i}$, but not including the hospital indicators. Supposing that the standard population is large enough to ignore estimation uncertainty in the expected counts, we can take $s_{z}^{2}=\frac{p^{2}}{E_{z}^{2}}\sum_{i=1}^{n}\mathbf{1}_{\{Z_{i}=z\}}\hat{\pi}_{i}(1-\hat{\pi}_{i})$, with $\tau^{2}$ reflecting the between-hospital variation after adjusting for patient case-mix. The latter in turn can be estimated using the DerSimonian & Laird method [9] or maximum likelihood, with the meta-analysis $I^{2}$ statistic indicating the proportion of variation due to between-hospital heterogeneity.

Closely related, but defined at individual level, is the concept of intra-class correlation (ICC), commonly estimated through mixed effect models [10, 11, 12, 13, 14, 15, 16]. Consider a generalized linear mixed effect model (GLMM)

[TABLE]

where $g()$ is the link function, and $Y_{i}$, $Z_{i}$, and $X_{i}$ are as before. The hospital level intercept terms are taken to follow $\alpha_{z}\sim N(0,\tau^{2})$, $z=1,\ldots,m$. For instance, with an identity link and residual variation distributed as $\varepsilon_{i}\sim N(0,\sigma^{2})$, we have that the conditional correlation $\textrm{cor}(Y_{i},Y_{i^{\prime}}\mid Z_{i}=Z_{i^{\prime}}=z,X_{i},X_{i^{\prime}})=\tau^{2}/(\tau^{2}+\sigma^{2})$. That is, the conditional correlation of the outcomes of two individuals treated in the same hospital is equal to the proportion of between-hospital variance of the total variance left after controlling for the patient case-mix, and can be estimated by fitting the random intercept model to obtain estimates for $\tau^{2}$ and $\sigma^{2}$. However, the connection between the ICC and the random intercept variance is more complicated with non-identity links, as the variance decomposition is different at the linear predictor and link function scale [17, 18]. An attempt for interpretation of the between cluster variance in terms of more familiar effect measures are the median odds ratio and median hazard ratio for binary outcomes and time-to-event outcomes, respectively [19]. We opt to take another route of considering the causal interpretation of the variance decomposition, as well as estimation methods that are invariant to the choice of the link function.

Furthermore, two-way variance decompositions generalize mathematically into arbitrary multi-way decompositions, though these depend on the chosen order of factorization [20]. Aligned with our goal of identifying QIs that can capture between-hospital variation in quality without being completely driven by the case-mix, we are interested in decomposing the observed variation in care received to that explained by patient case-mix, causal between-hospital variation given the case-mix and unexplained (residual) variation.

1.2 Objectives

Summarizing the objectives motivated above, the structure of the paper is as follows. In Section 2.1 we recapitulate the potential outcomes (Rubin’s) causal model [21], the related assumptions, and its extensions to between hospital comparisons following Varewyck et al. (2014) [4]. In Section 2.2, we introduce the causal variance decomposition to partition observed variation in care received, and discuss its interpretation in Section 2.3. In Section 2.4, we investigate the connection between our causal variance decomposition and the ICC. In Section 3, we propose model-based estimation methods for the variance decompositions, and investigate their performance in a simulation study in Section 4 and illustrate them in a real data application in Section 5. We conclude with a brief discussion on limitations and future directions in Section 6.

2 Proposed measures

2.1 Notation and assumptions

Suppressing the individual-level indices $i$ until the discussion of estimators in Section 3, let $Y\in\mathbb{R}$ represent the observed process or outcome experienced by a given patient, used to construct a quality indicator. Let $Y(z)\in\mathbb{R}$ represent the potential counterpart of this had the same patient received care in hospital $z\in\{1,...,m\}$. Let further $X=(X_{1},..,X_{p})$ be a vector of covariates relevant to the case-mix adjustment (including information on demographics, comorbidities, and disease progression), and let $Z\in\{1,...,m\}$ indicate the hospital in which the patient was actually treated. Under counterfactual consistency/stable unit treatment value assumption (SUTVA), the observed outcome is related to the potential outcomes by $Y=Y(Z)$. We assume the strong ignorability of the hospital assignment mechanism, which states that $0<P(Z=z\mid X=x)<1$ (positivity) and Y(z)\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}Z\mid X (conditional exchangeability) at all combinations of $z$ and $x$[21, 22]. We note that as long as these conditions are satisfied, the actual mechanism which assigns patients to hospitals can be more complicated without biasing the inferences. For example, adapting from Moreno-Betancur et al.[23], the causal relationships in the case of a process type indicator could be as illustrated in the directed acyclic graph (DAG) of Figure 1. Here the additional variables in history $H$ reflect the history (including family history) of an individual contributing to that individual living in a particular area indicated by $I$ (e.g. postal code), which in turn in part contributes to the individual being treated in a particular hospital in that area. In this DAG, $X$ are sufficient to control for confounding of the hospital effects, while $I$ are instrumental variables rather than confounders, and adjusting for them would likely lead to positivity violations.

2.2 Three-way causal decomposition for observed variation in care received

We are interested in decomposing the observed variance in the care received, $V[Y]$, to that causally explained by the quality of care differences between hospitals, by individual-level case-mix factors, and residual variation. Under the counterfactual consistency, $V[Y]=V[Y(Z)]$, for which we can write the usual two-way variance decomposition

[TABLE]

For the first term, we can further write

[TABLE]

and for the latter term,

[TABLE]

Substituting (3) and (4) into Equation (2), we obtain

[TABLE]

Here, due to strong ignorability, we have

[TABLE]

Applying these, we get the following result:

[TABLE]

We will consider the second term on the right had side as the causal quantity of interest. In this, the hospital specific performance is compared to the the expected level of care in the system for a patient with covariate values $X$, and then averaged over the covariate distribution of the standard population, similar to direct standardization. We note that this quantity is expressed in terms of potential outcomes and is a causal quantity in its own right, that is, it can be defined regardless of the causal assumptions. The causal assumptions are used to link this estimand to the observed data, working backwards to $V[Y]$. For the decomposition, we get the interpretation

[TABLE]

The estimation of the components is discussed in Section 3. Briefly, the estimation of the variance decomposition can be based on modeling the components $E[Y\mid Z,X]$ and $P(Z\mid X)$ combined with the empirical distribution of $X$. Alternatively, because the decomposition can be also expressed as

[TABLE]

the estimation can also be based on modeling $E[Y\mid Z,X]$ and $E[Y\mid X]$ combined with the empirical distribution of $(Z,X)$. In Section 3 we choose the former approach because it corresponds to the factorization of the likelihood, and thus enables approximate Bayesian inference. Before this, we further discuss the interpretation of the second term in the decomposition as a causal estimand.

2.3 Causal interpretation of the decomposition

To understand the causal interpretation of the variance decomposition, it is helpful to consider the special case of two hospitals, in which case the second term of the decomposition should reduce into a pairwise causal contrast. With two hospitals indexed by $z=1,2$ and denoting the propensity score $P(Z=1\mid X)=\pi(X)$, the first term of the decomposition becomes

[TABLE]

and the second term

[TABLE]

Now, we consider three different scenarios based on the relationships between $Y$, $Z$ and $X$ in Figure 1.

Scenario 1.

In the absence of the arrow $X\rightarrow Y$, which implies Y\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}X\mid Z and $E(Y(z)\mid X)=E(Y(z))$, the first term is equal to $V_{X}(\pi(X))\left[E(Y(1))-E(Y(2))\right]^{2}$. This implies that this term takes the variance contribution of the indirect (mediated) pathway $X\rightarrow Z\rightarrow Y$, which is present only if $Z$ has causal effect on $Y$. The second term is equal to $E_{X}[V(Z\mid X)][E(Y(1))-E(Y(2))]^{2}$. The difference in interpretation between these two terms is that the former represents the variance due to different kinds of patients being treated by different hospitals, while the latter represents the variation in care similar patients receive. If all hospitals treat similar patient populations, the former term disappears, while the latter is maximized. In the opposite case where the hospitals specialize in treating different patient populations, the positivity assumption is violated and the process/outcome cannot reveal performance differences between the hospitals. Thus, we focus on the second term as a measure of goodness of a proposed process/outcome as a QI.

Scenario 2.

In the absence of the arrows $X\rightarrow Z$ and $I\rightarrow Z$, which implies Z\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}X and $\pi(X)=\pi$ (randomized assignment), the first term would be equal to $V_{X}\left\{\pi[E(Y(1)\mid X)-E(Y(2)\mid X)]+E(Y(2)\mid X)\right\}$. The two additive components inside the variance reflect the effect modification by $X$ and the predictive variance due to $X$. In the absence of effect modification/additive interaction, that is, $E(Y(2)\mid X)=E(Y(1)\mid X)+\alpha$, it follows that $\pi[E(Y(1)\mid X)-E(Y(2)\mid X)]=\pi\alpha$ is constant and only the predictive variance remains. The second term is equal to $V(Z)E_{X}\left\{\left[E(Y(1)\mid X)-E(Y(2)\mid X)\right]^{2}\right\}$. In the absence of effect modification, this becomes $\pi(1-\pi)\alpha^{2}$, capturing the causal effect. Compared to the usual direct standardization formula $E(Y(1))-E(Y(2))=E_{X}[E(Y(1)\mid X)-E(Y(2)\mid X)]$, the between-hospital variance is similarly averaged over $X$, but squared so that any contrast at the same level of $X$ contributes to the variance rather than canceling out. The reason effect modification contributes to both case-mix and between-hospital variance components is that it both increases the variance of the expected care level (non-causal) and the average squared difference from the expected care level (causal). When evaluating proposed processes/outcomes, we would again want the second term to be comparatively large compared to the first one, because even though the eventual QI will be case-mix adjusted, large case-mix variance indicates that the measure is less intervenable, as the variation in care largely reflects treatment choices necessitated by disease progression, comorbidities etc., and makes the indicator sensitive to the method of case-mix adjustment.

Scenario 3.

In the absence of the arrow $Z\rightarrow Y$, which implies Y\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}Z\mid X and $E(Y(1)\mid X)=E(Y(2)\mid X)$, the first term becomes $V_{X}[E(Y(2)\mid X)]$, reflecting the predictive variance due to $X$, while the second term is equal to zero as it should in the absence of causal effect of $Z$ on $Y$.

2.4 Connection to intra-class correlation

We note that the previously derived variances were weighted by the quantities $P(Z=z\mid X)$ reflecting the hospital patient volumes; this is natural since we are decomposing the observed variance in the entire combined patient population, and larger centers carry more weight in this. However, to expand the proposed framework, we can consider variance decompositions under hypothetical assignment mechanisms that can differ from the actual one. For example, referring to the interpretation of Scenario 1 of Section 2.3, we can consider the variance decomposition under a hypothetical ‘randomized’ assignment mechanism where every hospital treats similar kind of patient population. This also allows us to derive a connection between the causal between-hospital variance in the previous section, and the more familiar intraclass correlation. Let $A$ represent a random draw from such a hypothetical assignment mechanism, with the mechanism chosen so that A\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}Z\mid X, and strong ignorability $0<P(A=a\mid X=x)<1$ and Y(a)\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}A\mid X apply. This lets us consider random draws of potential outcomes $Y(A)$, similar to the notation used by VanderWeele et al.[24]. Now, because

[TABLE]

we can write the variance decomposition

[TABLE]

If we choose here $P(A=a\mid X)=P(Z=a\mid X)\;\forall\;a\in\{1,\ldots,m\}$, the result is equivalent to before. However, if we choose to weight the hospitals equally as $P(A=a\mid X)=1/m$, and assume a linear model

[TABLE]

we can express the second term as

[TABLE]

If we take the hospital effects to represent a ‘sample’ from a distribution with variance $V[\alpha_{a}]=\tau^{2}$, (2.4) converges to $\tau^{2}$ when $m\rightarrow\infty$. In reality, since our discussion is in the administrative data context, we observe all the hospitals in the administrative region of interest and don’t take them to be a sample from a superpopulation of hospitals. Nevertheless, the result serves to establish a connection to the ICC: because

[TABLE]

under a linear model the causal interpretation of $ICC=\frac{\tau^{2}}{\tau^{2}+\sigma^{2}}$ is the proportion of case-mix conditional variance that is causally explained by the between-hospital variation in practices. However, this interpretation is specific to the linear model with equal hospital weights, which is why we proceed with the decomposition in Section 2.2 for the observed variance; this is general and can be estimated in the same way irrespective of the outcome model formulation (that is, irrespective of the choice of the link function, or whether the hospital effects are modeled as fixed or random effects).

3 Estimators

3.1 Point estimators

We present model-based direct standardization type estimators for the variance components. To model the outcomes, we specify a generalized linear model

[TABLE]

where $\theta=(\alpha_{0},\alpha_{1},\ldots,\alpha_{m},\beta)$. The hospital level intercept terms $\alpha_{z},z\in\{1,...,m\}$ are either fixed, with $\alpha_{1}=0$, or follow $\alpha_{z}\sim N(0,\tau^{2})$ (random intercept model). As noted in the previous section, the random effects are adopted as a means to apply shrinkage to the hospital effects, rather than a data generating mechanism. With random effects, the fitted values are obtained using empirical Bayes prediction for the random intercepts. We note that in the model we can also introduce hospital-case-mix interactions, as discussed by Varewyck et al.[25], but omit these here for notational simplicity. We estimate the parametrized variance explained by the patient case-mix as

[TABLE]

where $\mu_{i}(z;\theta)=E[Y_{i}\mid Z_{i}=z,x_{i};\theta]$ and where the terms $P(Z=z\mid X;\eta)$ can be estimated by fitting a multinomial logistic regression model

[TABLE]

where $\eta=(\gamma_{2},\ldots,\gamma_{m},\phi_{2},\ldots,\phi_{m})$. Similarly, the parametrized average variance explained by the hospital performance conditioning on patient case-mix can be estimated by

[TABLE]

The residual variance can then be estimated by subtracting the two estimated components from the empirical total variance, or using a distributional assumption such as $V(Y_{i}\mid Z_{i}=z,x_{i};\theta)=\mu_{i}(z;\theta)[1-\mu_{i}(z;\theta)]$ in the case of binary indicator, and taking

[TABLE]

3.2 Variance estimators

To quantify the uncertainty in the variance components, we note that conditional on the empirical covariate distribution in the standard population the estimates are entirely model-based. Thus, we can produce approximate samples from the joint posterior distribution

[TABLE]

by sampling $\theta$ and $\eta$ separately from their respective posteriors, and recalculate the variance components $\omega_{1}(\theta,\eta)$ and $\omega_{2}(\theta,\eta)$ at each combination of $\theta$ and $\eta$ to quantify uncertainty in them. Under flat priors the sampling distribution variances of the parameter estimates approximate the respective posterior variances; in approximate sampling of $\theta$ we used parametric bootstrap to reflect the uncertainty in the random effect estimates; for example for a binary indicator we can repeatedly sample the outcomes from $Y_{i}\mid(z_{i},x_{i};\hat{\theta})\sim\textrm{Bernoulli}(\mu_{i}(z_{i};\hat{\theta}))$, $i=1,\ldots,n$, and refit the mixed model to produce a sample of the model-based fitted values $\mu_{i}(z;\theta)$. Compared to non-parametric bootstrap this avoids some hospital categories being absent in some of the replicates. For sampling of $\eta$ we used the normal approximation $\eta\mid(\mathbf{X},\mathbf{Z})\sim N(\hat{\eta},V(\hat{\eta}))$ where $\hat{\eta}$ is the maximum likelihood estimate of $\eta$ and $V(\hat{\eta})$ its asymptotic variance-covariance matrix from the multinomial logistic model fit.

4 Simulation study

4.1 Generating mechanism

We demonstrate the performance of the estimators by simulating data from a generating mechanism simplified (omitting $H$ and $I$) from that in Figure 1, varying the total number of patients $n$ and number of hospitals $m$. Our objectives in the simulation are (a) to demonstrate that the proposed estimators work both for linear models for continuous outcomes and logistic models for dichotomous outcomes, (b) to demonstrate that the estimators work both with fixed and random effect models, and (c) to obtain evidence of their asymptotic behaviour under increasing $n$ and $m$. First, we generated two case-mix factors from $X_{1}\sim N(0,1)$ and $X_{2}\sim\textrm{Bernoulli}(0.5)$, given which the hospital assignment $Z$ was generated from multinomial logistic model specified as in (7), drawing the intercepts from $\gamma_{z}\sim N(0,0.25)$ and the regression coefficients from $\phi_{z1}\sim N(0,0.5)$ and $\phi_{z2}\sim N(0,0.5)$, $z=2,\ldots,m$. In our generating mechanism, we did not consider the hospitals themselves to be a sample from a superpopulation of hospitals, and thus we fixed the hospital-specific parameters across the replications (i.e. drawing these only once from the above distributions), while resampling the patients in each replication. This corresponds to conventional analysis of administrative data, where all the hospitals operating in a given region are compared, but the observed patient population is taken to represent a sample from an unobserved“long run”. Outcomes were generated from the mean structure

[TABLE]

with the hospital effects drawn from $\alpha_{z}\sim N(0,2^{2})$ (again fixing these parameter values across the replications). Continuous potential outcomes were generated by taking $Y(z)=E[Y(z)\mid X;\theta]+\varepsilon$, where $\varepsilon\sim\textrm{Logistic(0,1)}$, with the observed outcome given by $Y(Z)$. Binary outcomes were generated by dichotomizing at zero as $\tilde{Y}(z)=\mathbf{1}_{\{Y(z)\geq 0\}}$. To calculate the estimators, correctly specified outcome and assignment models were fitted to each simulated dataset, with both fixed and random effect specifications used for the former, whereas the true decomposition was calculated using the known values fixed in the data generating mechanism.

4.2 Results

4.2.1 Continuous Outcomes

Figure 2 shows the simulated sampling distribution means for the three variance components under both fixed and random effect models for the continuous outcomes under different combinations of $n$ and $m$, based on 1000 replications. Also shown are the 95% quantile interval for the sampling distribution, and the 95% confidence interval for the mean of the sampling distribution, reflecting the magnitude of the Monte Carlo error. The black dots indicate the true variances. The estimators capture the true value reasonably well in all cases, though with the random effect model giving slightly smaller estimates for the hospital variance component. Also of note is that the uncertainty in the variance component estimates is mainly driven by the overall number of patients $n$ (distributed across the $m$ hospitals) rather than the number of hospitals $m$. This is also illustrated in Figure 3 which shows the convergence of the sampling distributions for the hospital and patient variance components as a function of $n$ with fixed $m$. Finally, as suggested in Section 2.4, in Figure 4 we demonstrate that the random effect variance $\hat{\tau}^{2}$ estimates a different quantity than the causal between-hospital variance in the decomposition; the difference is marked especially with small number of hospitals.

4.2.2 Binary Outcomes

Figure 5 shows the simulated sampling distribution means for two of the three variance components under both fixed and random effect models for the binary outcomes under different combinations of $n$ and $m$, based on 1000 replications. The results are mostly similar to the continuous case, with the random effect model giving slightly smaller estimates for the hospital variance component compared to the fixed effect model. However, some small sample bias in the hospital variance component (fixed effect model overestimating and random effect model underestimating) is visible with small number of patients per hospital (top-right panel, on average $500/25=20$ patients per hospital), which disappears with increasing number of patients. The convergence of the sampling distributions with increasing $n$ is illustrated in the density plots of Figure 6.

5 Illustration in real data

We demonstrate the use of the proposed three-way causal decomposition in the context of QIs for surgical care of kidney cancer in the province of Ontario. We investigate four processes/outcomes that have previously been proposed as disease-specific QIs using the consensus/Delphi method by a panel of experts [2]. These are the proportion of partial (versus radical) nephrectomies among stage T1a nephrectomy patients (Partial), the same proportion restricted to the subpopulation of patients with chronic kidney disease (CKD) or its risk factors diabetes or hypertension (Partial-CKD), minimally invasive (versus open) surgery among T1-T2 radical nephrectomy patients (MIS), and readmission within 30 days of the surgery for T1-T4 radical nephrectomy patients (Readmission). The Partial, Partial-CKD and MIS indicators are process type, with higher proportion reflecting quality of care through kidney function conserving and less invasive treatment for early stage patients that can benefit from this. The Readmission indicator is outcome type, with lower proportion reflecting quality of care through reduced complications. Our objective in the illustration was to evaluate the processes/outcomes using the proposed variance decomposition approach, to quantify how much of the variation in these is between-hospital (reflecting performance and potentially intervenable), and how much is driven by the patient case-mix. The data were obtained from cross-linkage of Cancer registry, pathology reports, Hospital Discharge Abstract (DAD), Ontario Health Insurance Plan (OHIP) and other databases housed at the Institute for Clinical Evaluative Sciences (ICES). We included patients who had complete data on the case-mix factors age, sex, income quintile, Charlson comorbidity score, ACG comorbidity score, presence of CKD risk factors (for indicators other than Partial-CKD), days from diagnosis to treatment, year of diagnosis, tumor size, and T stage. Table 1 shows the number $n$ of patients identified and hospitals $m$ evaluated for each of the four indicators in the study period from 1995-2014.

To these data, using the R package VGAM [27], we fitted a multinomial logistic assignment models where, to deal with empty covariate categories for small volume hospitals, we estimated only intercept terms for hospitals that had treated less than 35 (40 for MIS) patients over the study period. This is justifiable as our estimators use the assignment model only for weighting the hospitals by their patient volumes, whereas the case-mix adjustment is achieved through the outcome model. For the processes/outcomes, we fitted mixed effects logistic models using the glmer function of R package lme4 [26], adjusting for the aforementioned covariates. The covariate effects (log-odds ratios) and their z-scores are reported in Table 1, along with the likelihood ratio test $\chi^{2}$ values for the random intercept variance parameter in the models (one degree of freedom). In all cases the likelihood ratio test indicated statistically significant between hospital variation, but this was much smaller for Readmission. Larger tumor size was a very strong predictor of radical nephrectomy, while both partial nephrectomies and minimally invasive surgeries had become significantly more common over calendar time. The comorbidity scores were strongest predictors of readmission.

The empirical proportions of the four indicators and the corresponding marginal variances were 44.7% and 0.447*(1-0.447)=0.247 (Partial), 43.8% and 0.246 (Partial-CKD), 63.7% and 0.231 (MIS), and 7.6% and 0.07 (Readmission). Table 2 shows the estimated three variance components using the proposed decomposition, along with the 95% credible intervals using the approximate Bayesian method described in Section 3. We also show the proportional variances along with their credible intervals, because these are more informative for comparing indicators with very different marginal variance. We note that the estimated components sum up to the empirical marginal variance given by the binomial calculation. The Partial indicator shows some between-hospital variation, but this is small compared to the case-mix variance. This does not change much when this indicator is restricted to the CKD risk factors subpopulation, because tumor size is the major driver of the large case-mix variance. These results suggest that even though the indicator was restricted to T4a disease, the treatment choice seems to be often determined by disease progression. However, because partial proportion has seen significant increase over calendar time, the case-mix variance will be smaller when such indicator is applied to contemporary data. In contrast, the MIS indicator shows between-hospital variance that is almost two times the case-mix variance, implying lot of variation in care for similar patients. This may be due to specialization or facilities available, which are intervenable. Because the calendar year was the biggest driver of the case-mix variance for MIS, the case-mix variance will be even smaller when indicator is applied to contemporary data. Finally, the Readmission indicator shows very little between-hospital variation, suggesting it is unlikely to be useful in capturing performance related variation.

6 Discussion

In this paper, we aimed to give a causal interpretation to variance decompositions used to quantify variations in hospital performance, and propose estimators that can accommodate different link functions, as well as fixed and random effect models for the hospital effects. The estimation of the hospital variance component worked reasonably well with the number of hospitals as small as 10. While the random effect models could introduce some bias to estimation of the variance components when the hospital-level effects and the patient case-mix are dependent [28], they are helpful for smoothing purposes when a large number of small volume hospitals are present. Fixed effect models can accommodate interactions between case-mix factors and hospital effects, but this leads to a large number of parameters to be estimated [25].

Addressing the variance decompositions in an explicit causal modeling framework helps to answer questions regarding the part of the variation in care received that is explained by different practices for similar patients. Furthermore, quantifying explained confounding due to observed case-mix factors into its own variance component helps to assess the usefulness of a given process or outcome in constructing a quality indicator. Also, we note that, similar to random effect meta-analysis, the estimation of the between-hospital variance component is weighted roughly in proportion to the hospital-specific patient volumes, meaning that the unstable estimates of small hospitals do not have undue influence on the estimation of this variance component.

Several extensions are possible. We split the total variance into that explained by hospital performance, case-mix, and unexplained residual variance. However, we did not separate between case-mix factors that capture disease progression and thus directly (and justifiably) influence the treatment decisions, and those capturing disparities in care, such as sociodemographic and socioeconomic factors. While the latter can be confounders under the causal model, any variation in care received that is explained by such factors may be of interest in itself, and thus could be split into a separate variance component. Of particular interest might be interactions between hospital effects and sociodemographic/economic factors, which would imply that disparities in care occur within hospitals.

In this paper we considered only a single layer of exposure hierarchy, such as hospitals. We are currently working on extending the causal variance decompositions to multiple layers of a hierarchy, such as surgeons within hospitals or referral networks within administrative subregions. This will require introduction of multiple levels of exposure, using notation similar to causal mediation analysis models [29], where the causal pathway leads, for example, from hospital to surgeon to outcome. This will enable us to further decompose the hospital effects to those due to institution level factors and those resulting from practitioner level differences. However, since the models with nested exposures are made identifiable through random effects, this will require further consideration of the role of random effects in causal models.

Funding

This work was supported by a Discovery Grant from the Natural Sciences and Engineering Research Council of Canada (to OS), a Catalyst Grant in Health Services and Economics Research from the Canadian Institutes of Health Research (to AF, KAL and OS) and the Ontario Institute for Cancer Research through funding provided by the Government of Ontario (to BC).