2.1 Notation and assumptions

For the extent of the paper, Y∈{0,1} is the observed binary quality outcome variable, Z∈{1,…,m} is the hospital in which the patient was actually treated, and X≡(X1,…,Xp) is a vector of patient-level characteristics relevant to case-mix adjustment, capturing for example demographic information, medical history, and disease progression. The triples W≡(Y,Z,X) are assumed independent and identically distributed across the patients. As is the convention in the causal inference literature, we denote by Yz the potential outcome that would have been observed had the patient been treated in hospital z. Throughout, we make the following standard causal assumptions. We assume that X is sufficient to control for any confounding (conditional exchangeability) so that (Y1,…,Ym)\upmodelsZ∣X. In addition, we assume consistency, under which the observed outcome is determined by Y=∑z=1m1{Z=z}Yz (Hernan and Robins, 2006), where 1{Z=z} is the indicator function which takes on values {0,1} depending on if the condition is false or true respectively. Finally, we assume positivity, under which all patients have a non-zero probability of being treated at any hospital, i.e. P(Z=z∣X)>0 for all z∈{1,…,m} and X combinations.

2.2 Direct versus indirect standardization

In this section, we briefly review the two common standardization procedures used in epidemiology, and discuss their causal interpretation in the quality comparison context. The main difference between the two standardization methods is that direct standardization provides the expected outcome if the standard population were to experience the event rate observed in the index population, whereas indirect standardization provides the expected outcome if the index population had experienced the event rate from the standard population. Table 1 provides an illustration of the different elements from each population used to compute the expected outcome. Each standardization method computes the expected outcome by considering a covariate stratum-specific (e.g. age, gender) event rate applied to a stratum-specific population size.

Direct standardization, as seen in the first row of Table 1, takes the event rate of the index population and applies it to the standard population in each strata, and then averages over all covariate strata in the standard population, resulting in E=∑knk∗−1∑knk∗πˆk. Direct standardization assumes that the stratum membership in the standard population is known and the stratum-specific rates must be estimated from the index population. In the present context, the index population are patients treated in a given hospital, while the standard population may be another hospital, or in the case of nationwide comparisons, the entire patient population across all hospitals. The case-mix adjustment required for the quality comparisons usually involves a large number of covariate strata, and therefore the stratum-specific event rates are in practice found using regression modelling techniques. The causal estimand under direct standardization where the standard population is all hospitals nationwide, can be written as E[Yz] as in Varewyck et al. (2014), which under the causal assumptions of Section 2.1 can be expressed as E[Yz]=∑xE[Y∣Z=z,x]P(X=x)e.g.(Hernan and Robins, 2006). This can be interpreted as the expected outcome had all patients experienced the care level of hospital z.

In contrast, indirect standardization takes the event rate in each covariate stratum in the standard population, applies it to the stratum membership in the index population, and then averages over the membership of all strata in the index population, as E=∑knk−1∑knkπˆk∗ (see row 2 of Table 1). More conceptually, indirect standardization can be thought of as determining the expected outcome if patients in the index population were to experience the same event rate as the standard population. In this case, the stratum-specific event rates are obtained from the standard population, while the stratum membership is determined from the study index population. Because the stratum-specific event rates do not need to be estimated from the index population, indirect standardization is still applicable when the index population is small, provided that the standard population is large. Finally, the expected event counts are contrasted to the observed ones through the standardized mortality ratio SMR=O/E.

The specific causal comparison being made in indirect standardization is determined by the choice of standard population. Suppose for instance that the comparison of interest is how patients from hospital z (index population) would fare if they were instead treated at hospital z′ (standard population). The causal estimand for indirect standardization must feature a conditional expectation of the form E[⋅∣Z=z]. When comparing two hospitals, this becomes simply E[Yz′∣Z=z], with SMR=E[Yz∣Z=z]/E[Yz′∣Z=z]. The latter corresponds to the exposure effect among the exposed risk ratio discussed by Shinozaki and Matsuyama (2015). However, to express the causal estimand in comparison to the nationwide average care level, instead of fixing the subscript of the potential outcome, we need to consider this as a random variable. This corresponds to a hypothetical intervention of randomly assigning a patient actually treated in hospital z to be treated in one of the hospitals under comparison. To express this, as a notational device, we define a new random variable A to denote the unobserved potential hospital assignment, where A∈{1,…,m}. Further, we let (Y1,…,Ym)\upmodelsA∣(Z,X) and A\upmodelsZ∣X so that we have the causal relationships presented in Figure 1. We can now generally express the causal estimand in indirect standardization through the conditional expectation E[YA∣Z=z]. Several interesting special cases may be obtained by choosing the hypothetical assignment probabilities P(A∣X). The usual exposure effect among the exposed comparison would be obtained by taking P(A=z′∣X)=1. Comparison to the care level of an average provider would be obtained by taking P(A=a∣X)=1/m for all a∈{1,…,m} (see Section 2.4 and Varewyck et al., 2014). However, herein we are specifically interested in comparisons to nationwide care level. In Section 2.4 we show that the corresponding causal estimand is specified by choosing the hypothetical assignment probabilites as equal to the actual ones, effectively weighting the hospitals in the average by their patient volumes.

Figure 1

The postulated causal mechanism (U is a non-confounder latent variable representing the correlation between potential outcomes for an individual).

Regardless of the standardization method used, it is necessary to estimate the respective expected potential outcomes. When the number of covariate strata becomes too large, it is common to fit an outcome model to the data and to estimate the particular expected number from the fitted values. However, such estimates are subject to the possibility of bias from model misspecification due to any number of factors. One attempt to protect against misspecification of the outcome model used in estimation is to instead use a doubly robust estimator.

2.3 Doubly robust estimation in direct standardization

Doubly robust (DR) estimation attempts to eliminate bias due to misspecification of a single model by utilizing two separate models in the estimation process (Funk et al., 2011). By doing so, it incorporates as much information about the causal pathway between the outcome, the exposure and the covariates/confounders as possible (see Figure 1). DR estimators combine the use of an outcome model, m(X,z,ϕ)≡E[Y∣Z=z,X;ϕ], and a propensity/assignment model, e(X,z,γ)≡P(Z=z∣X;γ), parametrized with respect to ϕ and γ respectively, into an estimator such that, as long as one of the models is correctly specified, the results should be unbiased or at least consistent (Bang and Robins, 2005).

In general, DR estimators of a mean effect are composed of three terms that contain either the outcome model, the propensity/assignment model, or both such that the terms that contain the misspecified model will cancel thereby resulting in estimation using only the correct model. Therefore, making use of fitted outcome and assignment probabilities m(X,z,ϕˆ) and e(X,z,γˆ), where estimated model parameters were denoted by ϕˆ and γˆ respectively, a DR estimator (Robins et al., 2007) for a marginal mean μz under direct standardization can be written as

Here the outcome model estimator, with hospital effect, has been augmented by weighting by the inverse of the assignment probability. An estimator of this form is doubly robust if either the outcome or the assignment/propensity model is correctly specified. This is evident as the second term in eq. (1) will in large samples have mean zero if the outcome model is correctly specified, leaving the first term to provide the estimate, while the second term of eq. (2) will also in large samples have mean zero if the propensity/assignment model is correctly specified and thus leaving only the propensity model to provide the estimate. Varewyck et al. (2014) used such an estimator for estimating the potential full population risk, E[Yz], had all patients received the care level of hospital z, a directly standardized quantity.

An issue that arises in direct standardization is the need to specify hospital effects in the outcome model. In the case of a large nationwide database of hospitals, some of them small in volume, this requires the estimation of a large number of parameters which might not be feasible without smoothing/shrinkage. Although such smoothing could be employed through mixed effect models, this might be a questionable approach if the purpose of the institutional comparison is to identify outliers. Varewyck et al. (2014) discuss possible ways to reduce shrinkage, such as clustered mixed effect models on hospitals or Firth corrected fixed effects logistic regression as the outcome model.

Additionally, direct standardization raises the question as to whether modelling hospital-patient interactions would also be needed, especially in the case of hospitals that specialize in particular patient subgroups, such as children or the elderly (Varewyck et al., 2016). Including interaction terms would further contribute to the large number of parameters that require estimation in the outcome model. Varewyck et al. (2016) have shown that the omission of hospital-patient interactions in the models used for standardization can contribute bias towards the estimated excess risks.

To sum up, in the case of nationwide comparisons involving hundreds or thousands of hospitals, many of these small volume, direct standardization may be more problematic due to the large number of hospital effects and patient-hospital interaction effects that would need to be estimated to ensure unbiased estimation. However, if comparison to an average level of care as the reference is of interest, indirect standardization avoids the issue of modelling hospital effects as well as hospital-patient interactions, which we will demonstrate in Section 2.4. Nevertheless, indirect standardization still requires specification of an outcome model. We thus propose a doubly robust estimator for the standardized mortality ratio under indirect standardization. In order to do this, we must first express the SMR as a causal estimand.

2.4 Causal estimand under indirect standardization

As per the discussion in Section 2.2, we define the causal estimand for hospital z in indirect standardization as

where the observed response in the numerator results simply from considering the potential outcomes of the patients of hospital z had they been treated in hospital z (i.e. the consistency assumption). In contrast, the expected response in the denominator depends on the specified target assignment regime, P(A∣X). As mentioned in Section 2.2, notable special cases may be obtained by choosing the assignment probabilities. First, consider the target assignment regime that gives patients an equal probability of being treated at each hospital, P(A=a∣X)=m−1. We then show in Appendix A that, for this choice of assignment regime, eq. (3) is equivalent to

which is the causal estimand briefly considered by Varewyck et al. (2014). This takes an equally weighted average across all hospitals in the denominator and thus corresponds to using the care level of an average provider as the reference in indirect standardization. In contrast, we want to use the national average level of care as the reference, and choose as the target assignment regime P(A=a∣X)=P(Z=a∣X). Then the denominator of eq. (3) involves an average across all hospitals but weighted by their actual volume. For the causal estimand for hospital z under this special case, we introduce the shorthand notation θz≡SMR.

Now, utilizing the causal assumptions listed in Section 2.1, and the additional conditional independence properties (Y1,…,Ym)\upmodelsA∣(Z,X) and A\upmodelsZ∣X, it can be shown (Appendix A) that eq. (3) can be expressed in terms of observable quantities as

where the denominator corresponds to using an outcome model without hospital effects, and averaging the predictions from such a model over the patients of hospital z. This is similar to the indirect standardization approach considered by e.g. Faris et al. (2003) and Tang et al. (2015), and is the appropriate modelling approach when the reference is chosen as the average level of care in the health care system. While it depends on the context whether this is the relevant comparison, we will now demonstrate how to obtain a simple doubly robust estimator for the causal SMR under such standardization. We show in Appendix A that similar manipulations of the causal estimand (that resulted in the equivalence of eqs. (3) and (5)) further result in two other equivalent expressions in terms of observable quantities, that is,

and

Thus, under the causal assumptions, the causal estimand θz can be written in terms of observable quantities in three equivalent forms that could be estimated using either an outcome model (eq. (5)), an assignment model (eq. (6)), or a combination of both models (eq. (7)). Therefore, we may utilize all three forms in a doubly robust estimator.

2.5 Proposed doubly robust estimator

As eqs. (5–7) are all equivalent and contain either one or both of an outcome model and propensity/ assignment model, we may now write the causal SMR as

This motivates the proposed DR estimator for the SMR of hospital z under indirect standardization

Here m(xi,ϕ)≡E[Yi∣Xi=xi,ϕ] is an outcome model parametrized in terms of ϕ. In the case of a binary outcome variable, this would be a logistic regression model of the form

where ϕ≡(ϕ0,ϕ1). The corresponding parameter estimates are denoted by ϕˆ. Further, e(xi,z,γ)≡P(Zi=z∣Xi=xi,γ) is a multinomial logistic assignment probability model parametrized in terms of γ, given by

and e(xi,1,γ)=1−∑z=2me(xi,z,γ), with γ≡(γ02,…,γ0m,γ12,…,γ1m) denoting the collection of all the parameters, and the corresponding parameter estimates denoted by γˆ.

The estimator in eq. (9) is applied in turn to each hospital z=1,…,m, with the parameters ϕ and γ estimated through fitting the regression models in eqs. (10) and (11) to the pooled patient population. We note that the outcome model m(xi,ϕ) no longer contains a term for the hospital effect, (as opposed to m(xi,z,ϕ) in eq. (1)), and thus we are estimating fewer parameters compared to the outcome model used for direct standardization. However, the hospital assignment model requires estimation of hospital-level regression parameters, and thus it is worth considering the case where the observational database may contain information on small hospitals, in which very few patients are being treated. In such situations, there may not be sufficient data to estimate the model parameters for all hospitals in the multinomial assignment model. However, we still want to include all the hospitals in the standardization since the reference is the national average level of care. We therefore also propose a modification of the multinomial assignment model of eq. (11) that only specifies covariate effects for the hospitals that are large enough. Suppose that, out of m hospitals, the first l hospitals are ‘small’ and the rest are ‘large’. We may then pool the small hospitals together as the reference category and specify the multinomial assignment model as

and e(xi,1,γ)=1−∑z=2me(xi,z,γ). While the assignment model in eq. (12) is obviously misspecified for the small hospitals z=2,…,l, it will still help in estimation of the SMRs for the large hospitals, and thus is an improvement over using only the outcome model, as the estimation of the θzs for z=l+1,…,m will be doubly robust.

While the DR estimator in eq. (9) is composed of three terms as required by the form of eqs. (1) and (2), as well as using two models for estimation, it is worth noting that the assignment model is not actually being utilized through inverse probability weighting. Nevertheless, our estimator does in fact have the doubly robust property, namely that θˆz converges in probability to θz as long as either the outcome or assignment model is correctly specified, the proof of which can be found in Appendix B. We will demonstrate this property in a simulation study in the following section.