Research Article Received 16 December 2012,
Accepted 20 December 2013
Published online 17 January 2014 in Wiley Online Library
(wileyonlinelibrary.com) DOI: 10.1002/sim.6089
Methods for comparing center-specific survival outcomes using direct standardization Kevin He* † and Douglas E. Schaubel The evaluation of center-specific outcomes is often through survival analysis methods. Such evaluations must account for differences in the distribution of patient characteristics across centers. In the context of censored event times, it is also important that the measure chosen to evaluate centers not be influenced by imbalances in the center-specific censoring distributions. The practice of using center indicators in a hazard regression model is often invalid, inconvenient, or undesirable to carry out. We propose a semiparametric version of the standardized rate ratio (SRR) useful for the evaluation of centers with respect to a right-censored event time. The SRR for center j can be interpreted as the ratio of the expected number of deaths in the total population (if the total population were in fact subject to the center j mortality hazard) to the observed number of events. The proposed measure is not affected by differences in center-specific covariate or censoring distributions. Asymptotic properties of the proposed estimators are derived, with finite-sample properties examined through simulation studies. The proposed methods are applied to national kidney transplant data. Copyright © 2014 John Wiley & Sons, Ltd. Keywords:
center effect; Cox regression; survival analysis; standardized rate ratio; stratification
1. Introduction In many situations, interest lies in the comparison of survival outcomes by center (e.g., treatment facility, hospital, or other entity serving as healthcare provider). Center-specific evaluations can be carried out on a regular basis (e.g., annually) to identity centers with poor performance. Alternatively, a retrospective evaluation over a longer period could be used to identify centers with exceptionally good results, with the goal of specifying best practices. In other cases, comparisons across centers may be an interesting secondary analysis; for example, in a multi-center study to evaluate the impact of a specific treatment on mortality. An accurate comparison of center-specific survival outcomes needs to account for imbalances in risk factor distributions among centers. For instance, the inclusion of high-risk patients by a given center can make that center’s survival appear substandard. In addition, in the context of survival times subject to censoring, the same phenomenon can occur because of differences in center-specific censoring distributions. In this report, we propose a semiparametric version of direct standardization, suitable for mortality comparisons by center. The proposed approach involves first fitting a Cox [1] regression model, stratified by center. The regression model is not used to directly estimate center effects, but rather to ensure that adjustment covariate effects are not confounded by center (which could occur in the absence of adjustment for center). For each center, j , the standardized rate ratio (SRRj ) is then computed as the ratio of expected to observed numbers of deaths, where ‘observed’ refers to the total number of deaths across all centers (j D 1; : : : ; J ) and ‘expected’ represents the number of total deaths estimated to occur if in fact all centers had mortality hazard equal to that of center j . Because of the use of direct standardization, the fSRR1 ; SRR2 ; : : : ; SRRJ g can be compared (and validly ordered) because the same
2048
Department of Biostatistics, School of Public Health, University of Michigan, Ann Arbor, MI 48109-2029, U.S.A. *Correspondence to: Kevin He, Department of Biostatistics, School of Public Health, University of Michigan, Ann Arbor, MI 48109-2029, U.S.A. † E-mail: [email protected]
Copyright © 2014 John Wiley & Sons, Ltd.
Statist. Med. 2014, 33 2048–2061
K. HE AND D. E. SCHAUBEL
covariate and censoring distribution is applied to each, that is, the total study population serves as the standard. This is in contrast to indirectly standardized measures, such as the standardized mortality ratio (SMR). The motivating example for this study involves the evaluation of center-specific post-transplant mortality for kidney transplant patients. Examples of factors known to strongly affect the post-transplant mortality hazard include age, primary renal diagnosis, and pre-transplant time on dialysis; each of which may differ in distribution by center. Centers with mortality significantly greater than the national average may be subject to various degrees of intervention, including site visits and perhaps de-accreditation. Given the high stakes of such evaluations, it is important that the statistical methods used for identifying outlying centers be accurate. A commonly used measure for evaluating center-specific survival is the SMR, defined as the ratio of the observed number of deaths at a given center to the number expected if the center had mortality equal to the population average. The SMR is a tool familiar to fields such as epidemiology (for example, [2–5]). The comparison of observed and expected outcomes is also commonly used for healthcare regulation (for example, [6, 7]). In the context of renal research, Wolfe et al. [8, 9] calculated SMRs among kidney transplant patients using mortality tables published by the United States Renal Data Systems. To relax the assumption of known standard mortality rates, Dickinson et al. [10] studied a semiparametric SMR based on Cox regression. Limitations of the SMR include its use of indirect standardization; each center’s SMR is essentially adjusted to a different (center-specific) covariate and censoring distribution. Centers cannot be rank ordered on the basis of indirect standardization, because two centers with equal covariate-specific mortality hazards could have different SMRs, merely due to differences in their respective covariate or censoring distributions. Therefore, although SMRs are potentially useful for internal evaluation (e.g., for centers to evaluate themselves or for a governing body to evaluate this center’s mortality comparing to that expected at the national level), they are less useful in the work of external evaluation (e.g., for surgeons and patients to compare center-specific results in the same region), because comparisons of center-specific results would play at least some role. We provide additional commentary on the SMR in Section 5. Earlier, we note limitations in the SMR approaches. Naturally, this technique has its role in analysis contrasting centers. However, noting room for improvement, we propose semiparametric version of direct standardization useful for survival analysis of centers. Direct standardization is also a commonly used approach in comparisons of mortality, usually through a measure termed the standardized rate ratio (SRR) or comparative mortality figure (CMF). One can express the SRR as the ratio of expected to observed numbers of deaths in the whole study population; the numerator of the SRR represents the expected number of deaths if all patients were treated at the given center, while the denominator equals the total observed number of deaths in the study population. Breslow and Day [3] and Hazel [4] compared the SRR with SMR in the framework of person-year methods. In particular, the main drawback of SMR (with respect to comparisons across centers) is not inherent to the SRR, because the same standard population is applied to all centers. Hence, center-specific SRRs are directly comparable. The SRR has a long history in fields such as epidemiology, and there are many settings in which direct standardization is appropriate. The Cox model has dominated applications involving regression analysis of censored data since its development. Thus, the use of Cox regression to estimate directly standardized center effects is a natural choice. The main contribution of this report is to formalize procedures for the Cox regression-based SRR, including rigorous derivation of asymptotic properties, simulation studies, and detailed comparisons to alternative approaches.
2. Methods
Copyright © 2014 John Wiley & Sons, Ltd.
Statist. Med. 2014, 33 2048–2061
2049
First, we provide the notation to be used in this article. Let Ti and Ci represent the survival and censoring time, respectively, for the ith patient, where i D 1; : : : ; n. Let J be the number of centers. The P total number of subjects is denoted by n D JjD1 nj , where nj is the number of subjects in center j . Observation times are denoted by Xi D Ti ^ Ci , with at-risk indicator Yi .t / D I.Xi > t /, where a ^ b D minfa; bg and I.A/ is an indicator function taking the value 1 when condition A holds and 0 otherwise. The observed death indicators are denoted by i D I.Ti 6 Ci /, and the death counting process is defined as Ni .t / D i I.Xi 6 t /. Let Gi denote the center for subject i and set Gij D I.Gi D j /. Correspondingly, we set Yij .t / D Yi .t /Gij and Nij .t / D Ni .t /Gij . The observed data consist of n independent vectors, .Xi ; i ; Gi ; Zi /, where Zi is a vector of adjustment covariates.
K. HE AND D. E. SCHAUBEL
The assumed center-stratified Cox model can be formulated as ˚ ij .t / .t jZi ; Gi D j / D 0j .t / exp ˇ0T Zi ;
(1)
where 0j .t / is an unspecified center-specific baseline hazard function and ˇ0 is a parameter vector. The partial likelihood estimator [11] of ˇ0 is denoted by ˇO and is given by the solution to U.ˇ/ D 0, where U.ˇ/ D
n Z J X X j D1 i D1
˚
Zij Z j .uI ˇ/ dNij .u/
0
P Sj.1/ .uI ˇ/ and Sj.d / .uI ˇ/ D n1 niD1 Yij .u/Zi˝d expfˇ T Zi g for Rt O 0j t I ˇO , where d D 0; 1; 2: The Breslow estimator [12] of ƒ0j .t / D 0 0j .u/du is then given by ƒ 1
with Z j .uI ˇ/ D Sj.0/ .uI ˇ/
n
X O 0j .t I ˇ/ D 1 ƒ n
Z
t 0
i D1
dNij .u/ Sj.0/ .uI ˇ/
:
2.1. Indirect standardization: standardized mortality ratio We begin by introducing an alternative to the proposed measure. The SMR for center j is calculated as the ratio of observed to expected numbers of events
1
Oj .t / ; Ej .t /
S MRj .t / D where the numerator is given by Oj .t / D computed as Ej .t / D
n Z X i D1
t
Pn
i D1
(2)
Nij .t / and the expected number of events is
n o O 0 uI ˇO ; Yij .u/ exp ˇO T Zi d ƒ
0
O 0 t I ˇO representing an estimator of the average cumulative baseline hazard with ƒ n
X O D1 O 0 .t I ˇ/ ƒ n i D1
Z
t 0
dNi .u/ S .0/ uI ˇO
(3)
and ˇO is based on model (1). The estimator given in (2) was developed in [13], where part of the innovation was the proposed use of the center-stratified model to estimate ˇ0 . Intuition would suggest computing Ej .t / using an unstratified model. However, estimation of ˇ0 in the absence of adjustment for center effects may produce a substantially biased estimate due to confounding by center or merely due to center affecting the hazard function, with or without confounding (i.e., due to the nonlinear link function). The calculation of S MRj .t / involves two stages. A stratified Cox model (1)) is fitted in the (model O 0 t I ˇO , and then computed in first stage, with the population average cumulative baseline hazard, ƒ
1
the second stage (e.g., through a unstratified Cox model with no covariates) using ˇO 0 Zi as an offset. At the second stage, Oj .t /, Ej .t / are calculated. 2.2. Direct standardization: standardized rate ratio
2050
On the basis of a form of indirect standardization, the SMR described in Section 2.1 can be viewed as a weighted ratio of center-specific cumulative hazards, with weight functions based on center-specific Sj.0/ .t I ˇ/. These weight functions have an obvious disadvantage: they involve center-specific censoring and covariate distributions, which can differ considerably across centers. To rule out the possibility that differences among centers are due merely to different censoring and/or covariate distributions, the weight function should be specified such that differences among centers with respect to the resulting Copyright © 2014 John Wiley & Sons, Ltd.
Statist. Med. 2014, 33 2048–2061
K. HE AND D. E. SCHAUBEL
measure are a function only of corresponding differences in center-specific hazards. Motivated by such considerations, we propose an alternative method, referred to as the SRR, which can be interpreted as a semiparametric version of direct standardization. The proposed SRR is computed, for center j , as
1
SRRj .t / D with O.t / D
PJ
j D1
Ej .t / ; O.t /
(4)
Oj .t / being the total observed number of deaths across all centers and Ej .t / D
n Z J X X `D1 i D1
t
n o O 0j .u/ Yi ` .u/ exp ˇO T Zi d ƒ
(5)
0
representing the expected number of total deaths if all centers had mortality hazard equal to that of center j . Similar to the SMR, the SRR is easily interpreted and is well understood by clinical investigators. The SRR also involves a ratio of observed and expected numbers of deaths. However, the ‘expected’ component is in the SRR’s numerator, while the ‘observed’ count is in the denominator. With respect to interpretation, SRRj > 1 indicates that center j has a greater mortality rate than the overall average. Note that, although SRRj .t / also involves the censoring and covariate distributions, the same weight function is applied across all centers, thus factoring out the impact of imbalances in center-specific censoring and covariate distributions. The proposed measures are desirable in this light, because their center-specific limiting values would differ only due to corresponding differences in centerspecific hazards. 2.3. Asymptotic properties We summarize the asymptotic properties of the proposed SRR with the following theorem; we outline the proof in Appendix A.
1
Theorem 1 Under the regularity conditions listed in Appendix A, SRRj .t / converges in probability to SRRj .t / uniformly in t 2 Œ0; , where R t .0/ `D1 0 s` .uI ˇ0 /dƒ0j .u/ PJ R t .0/ `D1 0 s` .uI ˇ0 /dƒ0` .u/
PJ SRRj .t / D
1
o n 1 and n 2 SRRj .t / SRRj .t / converges weakly to a zero-mean Gaussian process with variance function j .t / D EŒij .t I ˇ0 /2 , where Z
t
ij .t I ˇ/ D w.t I ˇ/ 0
s .0/ .uI ˇ/ sj.0/ .uI ˇ/
dMij .uI ˇ/ C w.t I ˇ/
Z tn
o Yi .u/ exp.ˇ T Zi / s .0/ .uI ˇ/ dƒ0j .u/
0
(6) w.t I ˇ/2
Z
t
s .0/ .uI ˇ/dƒ0j .u/ 0
Z C w.t I ˇ/ 0
J Z t n X `D1
0
t
rjT .uI ˇ/dƒ0j .u/.ˇ/1
o dNi ` .u/ s`.0/ .uI ˇ/dƒ0` .u/
Z
(7)
fZi ` ´` .uI ˇ/g dMij .uI ˇ/
(8)
0
with w.t I ˇ/ D
( J Z X
Copyright © 2014 John Wiley & Sons, Ltd.
0
)1 s`.0/ .uI ˇ/dƒ0` .u/
2051
rj .uI ˇ/ D s
`D1 .1/
t
.uI ˇ/ s .0/ .uI ˇ/ ´j .uI ˇ/: Statist. Med. 2014, 33 2048–2061
K. HE AND D. E. SCHAUBEL
O D n1 Pn Oij .t I ˇ/ O 2 , with Oij .t I ˇ/ O The variance function can be consistently estimated by O j2 .t I ˇ/ i D1 obtained by replacing limiting values in ij .t I ˇ/ with their empirical counterparts. With respect to this variance formula, a potentially time-saving strategy is to treat the estimators of ˇ0 as constants and hence ignore their variability. A justification for such simplification applies when the total study population is very large, such that ˇO has little variability. The resulting variance estimator is given by Z t .0/ Z tn o s .uI ˇ/ R T .0/ ij .t I ˇ/ D w.t I ˇ/ .uI ˇ/ C w.t I ˇ/ .u/ exp.ˇ Z / s .uI ˇ/ dƒ0j .u/ Y dM ij i i .0/ 0 s 0 j .uI ˇ/ (9) 2
Z
t
w.t I ˇ/
s
.0/
.uI ˇ/dƒ0j .u/
0
J Z t n X `D1
0
o dNi ` .u/ s`.0/ .uI ˇ/dƒ0` .u/ ;
(10)
obtained by removing the line (8) in the formula of ij .t I ˇ/. The quantity OijR t I ˇO is obtained by replacing limiting values in ijR .t I ˇ/ with their empirical counterparts.
3. Simulation We evaluate the finite-sample properties of the estimators described in Section 2 through a series of simulation studies. Death times were generated from the Weibull model, ij .t / D ˛j j t j 1 exp.ˇ T Zi / for i D 1; : : : ; nj and j D 1; : : : ; 10, where Zi D .Zi1 ; Zi 2 ; Zi 3 /T . We set ˇ T0 D .ˇ1 ; ˇ2 ; ˇ3 / D .0:02; 0:5; 0:2/. There are J D 10 centers. The number of subjects within each center varied under different scenarios. Censoring times were generated from either a Uniform distribution or an exponential distribution. In order to compare direct standardization with indirect standardization in the framework of semiparametric models, SMR and SRR were calculated at t D 1; t D 2, and t D 3. Each data configuration was replicated 1000 times. 3.1. Setting 1: center-independent hazards The first simulation setting considered the case where the hazard functions are equal across centers for all t 2 Œ0; . The censoring and covariate distributions were chosen to be center-independent. Specifically, censoring times were generated from a Uniform (0.5, 10) distribution; Zi1 followed a Bernoulli (0.5) distribution; Zi 2 followed a logistic distribution with probability dependent on Zi1 ; and Zi 3 came from a Normal distribution with constant variance 25 and mean dependent on Zi1 and Zi 2 (e.g., EŒZi 3 jZi1 ; Zi 2 D 50C0:2 Zi1 0:5 Zi 2 ). Under this setting, the hazard functions are equal across centers, such that the limiting values of SRR equal 1 for each center. Results at time t D 3 are displayed in Table I. For all centers, the average estimated SRR was very close to 1. The average asymptotic standard
O j .t/; center-independent hazards, covariate and censoring distributions; Table I. Simulation setting 1: SRR t D 3. Theorem 1 (6)–(8)
Theorem 1 (approx) (9), (10)
2052
Center
( j , ˛j )
True
Bias
ESD
ASE
CP
ASE
CP
1 2 3 4 5 6 7 8 9 10
(1, 0.2) (1, 0.2) (1, 0.2) (1, 0.2) (1, 0.2) (1, 0.2) (1, 0.2) (1, 0.2) (1, 0.2) (1, 0.2)
1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000
0.001 0.004 0:000 0:006 0:007 0:003 0.007 0:001 0.002 0.002
0.154 0.160 0.151 0.163 0.149 0.155 0.162 0.152 0.158 0.157
0.153 0.154 0.153 0.153 0.153 0.153 0.154 0.153 0.154 0.153
0.96 0.94 0.95 0.93 0.95 0.94 0.93 0.95 0.95 0.94
0.153 0.153 0.153 0.153 0.153 0.153 0.154 0.153 0.152 0.153
0.96 0.94 0.95 0.93 0.95 0.94 0.93 0.95 0.95 0.94
ESD, empirical standard deviation; ASE, asymptotic standard errors; CP, coverage probability.
Copyright © 2014 John Wiley & Sons, Ltd.
Statist. Med. 2014, 33 2048–2061
K. HE AND D. E. SCHAUBEL
errors were generally close to the empirical standard deviations, while the empirical coverage probabilities (CP) were generally consistent with the nominal value. This held for both the variance estimator derived in Theorem 1 (with corresponding asymptotic expansion given by (6)–(8) and its approximation given by (9) and (10). In fact, the two sets of asymptotic standard errors are virtually indistinguishable. Results were generally consistent across the observation time distribution (data not shown). It is worth noting that whether or not the covariate or censoring distributions were center dependent has no influence on the results in this setting. Correspondingly, similar results were found for the case of center-dependent covariate and censoring distribution.
3.2. Setting 2: center-dependent hazards; center-independent covariate and censoring distributions For the second set of simulations, different values of ˛j and j were used, such that the hazard functions increased with increasing center number (j D 1; : : : ; 10). The covariate distributions were chosen to be center independent and were generated from the same distributions from setting 1. Results are provided in Table II for various sample sizes and censoring percentages. The proposed SRR appears to be approximately unbiased, with CP close to 0.95 for both the standard error based on Theorem 1 and that based on the approximation. When the center-specific sample size is small (e.g., 25), the empirical CPs were
O j .t/; center-dependent hazards, center-independent covariate and Table II. Simulation setting 2: SRR censoring distributions; t D 3. Theorem 1 (6)–(8)
Theorem 1 (approx) (9), (10)
Censoring
Center
( j , ˛j )
True
Bias
ESD
ASE
CP
ASE
CP
100
20%
50
20%
100
40%
50
40%
125 100 75 50 25 125 100 50 25 15
20%
2 4 6 8 10 2 4 6 8 10 2 4 6 8 10 2 4 6 8 10 2 4 6 8 10 2 4 6 8 10
(0.85, 0.08) (0.95, 0.16) (1.05, 0.24) (1.15, 0.32) (1.25, 0.4) (0.85, 0.08) (0.95, 0.16) (1.05, 0.24) (1.15, 0.32) (1.25, 0.4) (0.85, 0.08) (0.95, 0.16) (1.05, 0.24) (1.15, 0.32) (1.25, 0.4) (0.85, 0.08) (0.95, 0.16) (1.05, 0.24) (1.15, 0.32) (1.25, 0.4) (0.85, 0.08) (0.95, 0.16) (1.05, 0.24) (1.15, 0.32) (1.25, 0.4) (0.85, 0.08) (0.95, 0.16) (1.05, 0.24) (1.15, 0.32) (1.25, 0.4)
0.406 0.797 1.180 1.549 1.912 0.406 0.797 1.180 1.549 1.912 0.406 0.797 1.180 1.549 1.912 0.406 0.797 1.180 1.549 1.912 0.529 1.043 1.550 2.025 2.535 0.601 1.118 1.773 2.324 2.900
0.004 0.003 0.001 0.006 0.005 0.002 0:008 0.007 0.001 0.026 0.008 0.003 0.007 0.014 0.029 0.009 0.011 0:007 0.019 0.043 0.006 0:004 0:002 0:015 0:004 0.001 0.007 0:001 0.045 0:062
0.096 0.132 0.174 0.206 0.230 0.135 0.183 0.252 0.290 0.344 0.097 0.139 0.172 0.214 0.243 0.137 0.204 0.249 0.302 0.357 0.110 0.169 0.256 0.384 0.651 0.122 0.190 0.375 0.636 1.034
0.093 0.134 0.169 0.204 0.238 0.130 0.187 0.238 0.285 0.338 0.097 0.138 0.174 0.208 0.244 0.136 0.195 0.242 0.294 0.345 0.106 0.171 0.255 0.381 0.632 0.118 0.192 0.359 0.621 0.914
0.94 0.95 0.94 0.95 0.96 0.93 0.95 0.94 0.95 0.95 0.94 0.95 0.95 0.94 0.95 0.94 0.94 0.94 0.94 0.95 0.94 0.95 0.95 0.93 0.93 0.94 0.96 0.93 0.93 0.91
0.093 0.134 0.169 0.203 0.237 0.130 0.186 0.238 0.284 0.337 0.097 0.138 0.173 0.208 0.243 0.135 0.195 0.241 0.294 0.344 0.106 0.171 0.255 0.380 0.630 0.118 0.192 0.358 0.619 0.910
0.94 0.95 0.94 0.95 0.96 0.93 0.95 0.94 0.95 0.94 0.94 0.95 0.95 0.94 0.95 0.94 0.94 0.94 0.94 0.95 0.94 0.95 0.95 0.93 0.92 0.94 0.96 0.93 0.93 0.90
40%
ESD, empirical standard deviation; ASE, asymptotic standard errors; CP, coverage probability.
Copyright © 2014 John Wiley & Sons, Ltd.
Statist. Med. 2014, 33 2048–2061
2053
nj
K. HE AND D. E. SCHAUBEL
slightly underestimated. Such results suggest that center-specific sample sizes play an important role in the proposed methods. In particular, the minimum sample size (across centers) needs to be reasonably large. Given this concern, centers of size less than 20 were eliminated from the real data analysis in Section 4. Collectively, simulation results from settings 1 and 2 indicate that the proposed method is quite accurate at and away from the null. 3.3. Setting 3: center-dependent hazards, center-dependent censoring distribution As mentioned previously, estimators based on indirect standardization may be misleading if either the censoring or the covariate distributions are center-dependent. To illustrate this point, we performed simulations with the following three conditions: (i) the hazard functions were center dependent, while center 2j 1 had exactly the same hazard function as that of center 2j for j D 1; : : : ; 5; (ii) the censoring distributions for centers 2j 1 and 2j were substantially different. For center 2j 1, the censoring times were generated from a uniform distribution such that the censoring mainly occurred in the later stages, while for center 2j , the censoring times were generated from a uniform distribution for which the censoring tended to occur in the early stages; and (iii) the covariate distributions (reused from setting 2) were center independent. We compared the SRR and SMR in Table III for setting 3. With respect to the true values, the limiting values of SRR2j .t / and SRR2j 1 .t / are equal, as one would hope. This is not the case for S MR2j .t / and S MR2j 1 .t /, the differences being due to differences in the censoring distributions. With respect to the estimators themselves, both SRRj .t / and S MRj .t / are approximately unbiased. Note that the bias of S MRj .t / was calculated as the difference between it and its own limiting value.
1 1
1
1
1
1
1
3.4. Setting 4: center-dependent hazards; center-dependent covariate distributions We also performed simulations under a setting in which the distribution of the covariate vector differed by center. The setup for hazard functions and censoring distribution was the same as in setting 2, while centers 2j and 2j 1 had substantially different covariate distributions. In center 2j 1, the covariate Zi1 followed a Bernoulli (0.2) distribution, Zi 2 followed a Bernoulli (0.8) distribution, and Zi 3 came from a Normal distribution with mean 30 and standard deviation 10. In center 2j , Zi1 followed a Bernoulli (0.8)
Table III. Simulation setting 3: centers 2j and 2j 1 have equal hazards but different censoring distributions; t D 3. ( j , ˛j )
2054
Measure
Center
SMRj .t /
1 2 3 4 5 6 7 8 9 10
(0.8, 0.04) (0.8, 0.04) (0.9, 0.12) (0.9, 0.12) (1, 0.2) (1, 0.2) (1.1, 0.3) (1.1, 0.3) (1.25, 0.4) (1.25, 0.4)
SRRj .t /
1 2 3 4 5 6 7 8 9 10
(0.8, 0.04) (0.8, 0.04) (0.9, 0.12) (0.9, 0.12) (1, 0.2) (1, 0.2) (1.1, 0.3) (1.1, 0.3) (1.25, 0.4) (1.25, 0.4)
True
Bias
ESD
0.227 0.211 0.658 0.637 1.050 1.051 1.500 1.547 1.881 2.003
0:003 0:002 0.001 0.003 0.003 0:000 0.002 0.004 0:002 0:008
0.089 0.068 0.162 0.121 0.208 0.157 0.239 0.201 0.279 0.219
0.221 0.221 0.645 0.645 1.044 1.044 1.530 1.530 1.996 1.996
0:003 0:002 0.002 0.004 0.002 0:002 0.001 0.004 0:004 0:008
0.088 0.073 0.162 0.124 0.216 0.158 0.277 0.202 0.351 0.230
ESD, empirical standard deviation.
Copyright © 2014 John Wiley & Sons, Ltd.
Statist. Med. 2014, 33 2048–2061
K. HE AND D. E. SCHAUBEL
Table IV. Simulation setting 4: centers 2j and 2j 1 have equal hazards but different covariate distributions; t D 3. Measure
Center
( j , ˛j )
True
SMRj .t /
1 2 3 4 5 6 7 8 9 10
(0.8, 0.04) (0.8, 0.04) (0.9, 0.12) (0.9, 0.12) (1, 0.2) (1, 0.2) (1.1, 0.3) (1.1, 0.3) (1.25, 0.4) (1.25, 0.4)
SRRj .t /
1 2 3 4 5 6 7 8 9 10
(0.8, 0.04) (0.8, 0.04) (0.9, 0.12) (0.9, 0.12) (1, 0.2) (1, 0.2) (1.1, 0.3) (1.1, 0.3) (1.25, 0.4) (1.25, 0.4)
Bias
ESD
0.319 0.331 0.931 0.919 1.484 1.307 2.070 1.602 2.551 1.725
0:003 0:002 0.001 0.003 0.003 0:000 0.002 0.004 0:002 0:008
0.089 0.068 0.162 0.121 0.208 0.157 0.239 0.201 0.279 0.219
0.345 0.345 0.929 0.929 1.386 1.386 1.908 1.908 2.329 2.329
0:002 0:000 0.005 0.003 0:007 0:000 0:007 0.013 0:006 0:001
0.120 0.041 0.212 0.105 0.272 0.186 0.344 0.330 0.412 0.455
ESD, empirical standard deviation.
distribution, Zi 2 followed a Bernoulli (0.2) distribution, and Zi 3 was derived from a Normal distribution with mean 50 and standard deviation 10. Results based on setting 4 are given in Table IV. Trends are similar to those from setting 3 but much more pronounced. The limiting values of SRR2j .t / are equal to those of SRR2j 1 .t /, as one would expect. Conversely, the true values of S MR2j .t / and S MR2j 1 .t / are different, with the differences being quite pronounced for j D 3, j D 4, and j D 5. Moreover, it appears that SMR7 .t / > SMR10 .t /, which is misleading in the sense that SRR7 .t / < SRR10 .t /.
1 1
1
1
4. Application
1
1
Copyright © 2014 John Wiley & Sons, Ltd.
Statist. Med. 2014, 33 2048–2061
2055
We applied the proposed methods to investigate the performance of transplant centers with respect to post kidney transplant survival. Data were obtained from the Scientific Registry of Transplant Recipients (SRTR) and submitted by members of the Organ Procurement and Transplantation Network. The SRTR database contains information on all wait-listed candidates, transplant recipients, and organ donors in the USA. Included in the analysis were adult patients (>18 years of age at transplant) who underwent deceased donor kidney transplantation between January 2000 and December 2008. Adjustment covariates in this study included age, race, gender, diagnosis, donation after cardiac death, expanded criteria donor, body mass index, dialysis time, indicator of previous kidney transplant, and cold ischemia time. These variables have face validity from a clinical perspective and are based on a list of covariates used in SRTR. Transplant centers with sample size 620 and patients who received a living-donor transplant were eliminated from additional analysis. The final sample size was then n D 74; 088 from J D 217 centers across the USA. Failure time (recorded in years) was defined as the time from transplantation to graft failure or death, whichever occurred first. Graft failure was considered to occur when the transplanted kidney ceased to function. Stratified Cox regression was employed to model the hazard function. The indirectly standardized estimator, S MRj , was calculated using SAS PROC PHREG with an offset. The proposed directly standardized estimator, SRRj , was computed using SAS IML. Figure 1a represents the pairwise comparisons of the SMRs and SRRs. Figure 1b shows the standard error of these two measures. Figure 1c compares
K. HE AND D. E. SCHAUBEL
(b) (Scatter plot : standard error of SMR and SRR)
0.8
SE of SMR
SMR 1.0
0.0
0.0
0.4
2.0
1.2
3.0
(a) (Scatter plot : SMR VS SRR)
0.0
0.5
1.0
1.5
2.0
2.5
3.0
0.0
0.2
SRR
0.4
0.6
0.8
1.0
1.2
SE of SRR
200 150 100 50 0
order of centers based on SRR
(c) (Scatter plot : orders of centers based on SMR and SRR)
0
50
100
150
200
order of centers based on SRR
Figure 1. Evaluation of J D 217 kidney transplant centers. SMR, standardized mortality ratio; SRR, standardized rate ratio; SE, standard error.
2056
the orders of centers based on SRRs and SMRs. As shown, there are some discrepancies between these two measures. We applied bootstrapped techniques to evaluate whether the change in center-specific orderings (SMR versus SRR) exceeded that attributable to only sampling variation. Specially, we calculated the distribution of center-specific SMR orderings from 100 bootstrapped samples. We then calculated the 95% confidence intervals of these orders. On the basis of the bootstrapped resamples, for 20 of 217 centers, the 95% confidence intervals based on SMR does not cover the order based on SRR from the original dataset. Among these 20 centers, 10 had SRR significantly different from the national average. Using the asymptotic normality of the proposed estimators, we constructed the point-wise confidence intervals for SRR at t D 5 years (Figure 2). Center numbers are re-ordered by values of SRR. A total of 38 centers had observed number of events significantly lower than the expected calculation based on the national average hazards, while 28 centers were significantly above the expected. It is clear that the hazard functions varied among centers. Table V presents the pairwise comparison of the numbers and percentages of ‘outlier’ centers identified by p-values corresponding to their SMRs and SRRs (using tests of H0 W SMRj D 1 and H0 W SRRj D 1, respectively). A total of six centers changed ‘memberships’ based on these two measures. Specifically, two centers were flagged to be significant based on SMR but flagged to be normal based on SRR; on the other hand, four centers were flagged to be significant based on SRR but flagged to be normal based on SMR. Through fitting a sequence of logistic regression models (rotating the center indicators as the response variates), it was revealed that, for these six centers, approximately half of the adjustment covariates had distributions significantly different from the remaining centers. In addition, through fitting a Cox model using censored as the event, three of the six centers in question were significantly predictive of the censoring hazard. In summary, we do observe some differences when comparing center effects estimated through direct versus indirect standardization, and the strongest examples of such discrepancies appear to be due to differences in the center-specific covariate and censoring distributions, consistent with the concepts described earlier in this report. Copyright © 2014 John Wiley & Sons, Ltd.
Statist. Med. 2014, 33 2048–2061
SRR(5)
0
1
2
3
4
5
K. HE AND D. E. SCHAUBEL
0
50
100
150
200
Center
Figure 2. Evaluation of J D 217 kidney transplant centers: point estimates and 95% confidence interval of O j at t D 5 years. SRR, standardized rate ratio. SRR
Table V. Number and percentage of centers giving significant results under SMR and SRR. SRR SMR Nonsignificant Significant Column sum
Nonsignificant
Significant
Row sum
147 (67.8%) 2 (0.9%) 149 (68.7%)
4 (1.8%) 64 (29.5%) 68 (31.3%)
151 (69.6%) 66 (30.4 %) 217 (100%)
SMR, standardized mortality ratio; SRR, standardized rate ratio.
5. Discussion
Copyright © 2014 John Wiley & Sons, Ltd.
Statist. Med. 2014, 33 2048–2061
2057
We propose semiparametric methods for estimating standardized rate ratios, as a means of evaluating center-specific mortality through direct standardization. Large-sample properties are derived and shown through simulation to be appropriate in finite samples. A computationally faster variance estimator is proposed for the SRR and is shown to work practically as well as the full version. Application of the methods demonstrates several significant differences among kidney transplant centers in the USA. There is some judgement required in deciding when to use indirect standardization and when to use direct standardization. Indirectly standardized estimators, such as SMR, provide a valid approach to evaluate how does a center’s mortality compare to that predicted at the population level for the kinds of patients at this center. However, it is important to emphasize that center-specific SMRs should not be compared with one another (a caution that applies to all indirectly standardized rates). The SRR, a directly standardized measure, does not share this drawback. The SRRs for two given centers will be unequal only because the center-specific mortality hazards differ; direct standardization accounts for imbalance with respect to center-specific covariate and censoring distributions. The proposed SRR shares the SMR’s ease of interpretation but rectifies its key disadvantages and, hence, is a more appropriate choice in settings where mortality comparisons across centers are an objective. The degree to which the SMR and SRR are different will depend on the application. In some settings, the two may not agree well, while in others, they may be quite similar. The only way to know with certainty if SMR and SRR are equal would be to calculate both measures, which would not be a desirable option in many cases. In settings where, for a particular center, the SMR and SRR were unequal, it would be very difficult to claim that the SMR was correct, for the several reasons documented previously. Given the high stakes of evaluations by regulatory bodies, and the fact that the credibility of such organizations depends in part on the accuracy of their evaluations, it would appear the preferred analysis is the one that is most likely to be accurate.
K. HE AND D. E. SCHAUBEL
The proposed SRR is computed using a stratified Cox model, which makes no assumptions about the functional form of the impact of center on the hazard function. The stratification by center plays a major role. For instance, the expected number of deaths considers each patient’s covariate vector, such that the regression parameter must be estimated consistently. Unless the death hazard is conditionally independent of center given the covariates, covariate effect estimates will generally be biased if based on a model with a nonlinear link function and no accounting for center. This is an issue for the SMR as well, particularly because the Cox version of the SMR has historically been computed using a commonbaseline (i.e., unstratified) model. He [13] proposed modifying the SMR through stratification, leading to the quantity we denoted by SMR and used in simulations (Section 3) for comparisons with the SRR. Such properties were demonstrated empirically by He [13]; the magnitude of the bias (in the case of unstratified Cox models) increases when covariates are also center-dependent. Random effects models are an option for contrasting centers. Moreover, a Bayesian formulation may be an attractive approach for center effect studies. The advantage for random effects model is that this approach would allow for the inclusion of even very small centers. In contrast, when the number of events for a center is small, the estimated center parameter from a fixed effect model may be unstable. However, Kalbfleish and Wolfe [14] compared the properties of a fixed effect model (FEM) and a random effect model (REM) for the purpose of profiling kidney dialysis facilities under various conditions. Essentially, the REM estimates are shrunk toward overall mean and hence reduce the reported variation of facility performance. Second, the FEM has the highest statistical power to identify exceptional facilities, for a given false positive rate; identifying such extreme facilities is usually a main objective of center evaluations. Another issue for REM is the potential confounding effects when the patient risks are correlated with center effects. These findings suggest that a simple REM method may not be good enough and more sophisticated approaches are necessary. Further discussion of such issues is provided by Ohlssen et al. [15], who develop a more flexible random effects model using Bayesian nonparametric methods, in order to remedy the influence of outlying centers to which basic random effects models are susceptible. The direct standardization methods derived in this report could be extended in several useful directions. Perhaps most notably, it is often of interest to evaluate center effects in settings where the event of interest is recurrent (e.g., hospitalizations and infections). Furthermore, it would also be useful to develop methods based on direct standardization that can accommodate competing risks or dependent censoring. Direct standardization could also be applied to compare center-specific survival probability and restricted mean lifetime.
Appendix A To derive the large-sample properties for the SRR, we impose the following regularity conditions under the stratified Cox model: (a) .Xi ; i ; Gi ; Zi / are independent and identically distributed random vectors. (b) P .Xi > / > 0 where is a pre-specified time point. (c) Zi k have bounded total variation, that is, jZi k j < for all i D 1; : : : n and k D 1; : : : ; p, where
is R a constant and Zi k is the kth component of Zi . (d) 0 0j .t /dt < 1. (e) Continuity of the following functions: sj.1/ .t I ˇ/ D
@ .0/ @2 s .0/ .t I ˇ/ sj .t I ˇ/; sj.2/ .t I ˇ/ D @ˇ @ˇ@ˇ T j
and sj.0/ .t I ˇ/, where sj.d / .t I ˇ/ is the limiting value of Sj.d / .t I ˇ/ for d D 0; 1; 2, with sj.1/ .t I ˇ/ and sj.2/ .t I ˇ/ bounded and sj.0/ .t I ˇ/ bounded away from 0 for t 2 Œ0; . (f) Positive-definiteness of the matrix j .ˇ/: Z j .ˇ/ D vj .t I ˇ/sj.0/ .t I ˇ/0j .t /dt; 0
2058
vj .t I ˇ/ D
Copyright © 2014 John Wiley & Sons, Ltd.
sj.2/ .t I ˇ/ sj.0/ .t I ˇ/
´j .t I ˇ/˝2 ;
Statist. Med. 2014, 33 2048–2061
K. HE AND D. E. SCHAUBEL 1 .1/ sj .t I ˇ/
where ´j .t I ˇ/ D sj.0/ .t I ˇ/ (g) P .Gij D 1jZi / > 0.
is the limiting value of Z j .t I ˇ/.
Condition (a) is employed in the derivation of the weak convergence. Condition (b) is a standard identifiability requirement. Condition (c) leads to the boundedness of several quantities and is applicable in most practical applications. Conditions (d) and (e) are not essential but simplify our proofs. With respect to condition (g), the selection probability given covariates is nonzero for all centers. This condition guarantees that the sample size nj of each center goes to 1 as the total sample size n goes to 1. P O j .t / We first show that SRR ! SRRj .t / uniformly for t 2 Œ0; . The triangle inequality leads to Z t Z t .0/ .0/ O O O S` uI ˇ d ƒ0j uI ˇ s` .uI ˇ0 /dƒ0j .u/ 0 0 Z t Z t .0/ .0/ O O O O 0j uI ˇ O 0j uI ˇ 6 S` uI ˇ d ƒ s` .uI ˇ0 /d ƒ 0 0 Z t Z t O 0j uI ˇO C s`.0/ .uI ˇ0 /d ƒ s`.0/ .uI ˇ0 /dƒ0j .u/ 0
(A.1) (A.2)
0
P
! 0 uniformly in t , recall that To show that .A:1/ P S`.0/ .t I ˇ/ D n1 niD1 Yi ` .t / exp.ˇ T Zi / D Pn ŒI.X > t /I.G D `/ exp.ˇ T Z/, where Pn is P the empirical measure; that is, Pn ŒI.X > t /I.G D `/ exp.ˇ T Z/ D n1 niD1 Yi ` .t / exp.ˇ T Z/. The collection of all cells, Œt; 1/, in the real line is a VC class of index 2 and, hence, satisfies the entropy conditions for the Glivenko–Cantelli theorem [16]. The boundedness conditions ensures that fI.X > t /I.G D `/ exp.ˇ T Z/; t 2 Œ0; g belong to some Glivenko–Cantelli class; that is, a:s: P O O uniformly in t 2 Œ0; . Next, for ˇ, such that ˇO S`.0/ .t I ˇ/ ! s`.0/ .uI ˇ/ ! ˇ0 (e.g, [17, 18]), O 0 .t /; t 2 Œ0; g and an application of the dominant convergence thethe bounded condition of fƒ P
P
orem entails that .A:1/ ! 0 uniformly in t . Similarly, .A:2/ ! 0 uniformly in t . We already R t .0/ P R t .0/ O demonstrate that 0 S` .uI ˇ0 /d ƒ0j .u/ ! 0 s` .uI ˇ0 /dƒ0j .u/ uniformly for t 2 Œ0; . SimiR t .0/ P R t .0/ O 0` .u/ larly, 0 S` .uI ˇ0 /d ƒ ! 0 s` .uI ˇ0 /dƒ0` .u/ uniformly for t 2 Œ0; . The monotonicity o1 P nR o1 nR t t O 0` .u/ ! 0 s`.0/ .uI ˇ0 /dƒ0` .u/ and boundedness conditions ensure that 0 S`.0/ .uI ˇ0 /d ƒ P O j .t / uniformly for t 2 Œ0; . Therefore, we have that SRR ! SRRj .t / uniformly for t 2 Œ0; . To prove weak convergence, we use the following decomposition:
1
o n 1 n 2 SRRj .t / SRRj .t / ( R t .0/ ) R t .0/ O ƒ O O O 0j .uI ˇ/ S` .uI ˇ/d s` .uI ˇ0 /dƒ0j .uI ˇ/ 1 0 0 D n2 R t .0/ R t .0/ O ƒ O O 0` .u/ O S .uI ˇ/d 0 ` 0 S` .uI ˇ0 /d ƒ0` .uI ˇ/ ( R t .0/ ) R t .0/ s` .uI ˇ0 /dƒ0j .u/ s .uI ˇ0 /dƒ0j .u/ 1 R0t `.0/ : C n 2 R t 0 .0/ O ƒ O O 0` .uI ˇ/ S .uI ˇ/d s .uI ˇ0 /dƒ0` .u/ 0
0
`
(A.3)
(A.4)
`
First, we have that ) Z t .0/ .uI ˇ / s dN .u/ 0 ij i D1 O S .0/ .uI ˇ/ s .0/ .uI ˇ0 / `.0/ dƒ0` .u/ C op .1/ .A:3/ D w.t; ˇ0 / n O s .uI ˇ0 / 0 0 nSj.0/ .uI ˇ/ " n Z t .0/ o X s .uI ˇ0 / n .0/ 1 2 .uI ˇ / s .uI ˇ /dƒ .u/ D w.t I ˇ0 / n dN ij 0 0 0j j .0/ i D1 0 sj .uI ˇ0 / Z t 1 C rjT .uI ˇ/dƒ0j .u/ n 2 .ˇO ˇ/ C op .1/; 1 2
(Z
Pn
t
Copyright © 2014 John Wiley & Sons, Ltd.
Statist. Med. 2014, 33 2048–2061
2059
0
K. HE AND D. E. SCHAUBEL
where w.t I ˇ/ and rj .uI ˇ/ are defined in Section 2. Note that the second equality of the aforementioned argument is obtained through the functional delta method and Lemma 19.24 of [19]. On the basis of previously established empirical process theory for the Cox model (e.g., [17]), we have that ( J ) " J X X 1 1 n 2 .ˇO ˇ0 / D ` .ˇ0 / Gn ffZ ´` .X I ˇ0 /gI.G D `/ `D1
Z
0
`D1
fZ ´` .uI ˇ0 /g exp ˇ0T Z I.X > u/I.G D `/dƒ0` .u/ C op .1/;
where the op .1/ is uniform in t , and Gn is the empirical process defined by Gn f D Through the functional delta method,
p n.Pn P/f .
Rt
J X s .0/ .uI ˇ0 /dƒ0j .u/ .A:4/ D n Gn ŒN` .t / C op .1/: o2 PJ R t .0/ `D1 s .uI ˇ /dƒ .u/ 0 0` `D1 0 ` 0
Combining (A.3) and (A.4), o n 1 n 2 SRRj .t / SRRj .t / (Z ( J " ) Z t t .0/ X s .uI ˇ0 / D Gn w.t I ˇ0 / rjT .uI ˇ0 /dƒ0j .u/ ` .ˇ0 /1 dNj .u/ C .0/ 0 s 0 .uI ˇ / 0 j `D1
) Z J X T fZ ´` .uI ˇ0 /g exp ˇ0 Z I.X > u/I.G D `/dƒ0` .u/ fZ ´` .X I ˇ0 /g I.G D `/
1
0
`D1
Rt
3
.0/
s .uI ˇ0 /dƒ0j .u/ 7 o2 I.X 6 t /5 C op .1/I R t .0/ `D1 0 s` .uI ˇ0 /dƒ0` .u/
n PJ
0
Because the VC classes with finite index satisfy the entropy conditions for the Donsker theorem [16], fI.X > t /; t 2 Œ0; g and fI.X 6 t /; t 2 Œ0; g belong to some Donsker classes. The same holds for the bounded monotone stochastic process fN.t /; t 2 Œ0; g. Finally, the class of functions of Lipschitz transformations of Donsker classes is Donsker. Therefore, with the various bounded conditions and by applying the Donsker theorem, Theorem 1 follows.
Acknowledgements This work was supported in part by the National Institutes of Health grant 5R01-DK070869 and a grant from the Michigan Institute for Clinical and Health Research (MICHR). The authors thank the Scientific Registry of Transplant Recipients (SRTR) for allowing us to access the organ failure database.
References
2060
1. Cox DR. Regression models and life tables (with Discussion). Journal of the Royal Statistical Society, Series B 1972; 34:187–200. 2. Berry G. The analysis of mortality by the subject-years method. Biometrics 1983; 39:173–184. 3. Breslow NE, Day NE. The standardized mortality ratio. In Biostatistics: Statistics in Biomedical, Public Health and Environmental Sciences, Sen PK (ed.), The Bernard G. Greenberg Volume, 1985; 55–74. 4. Hazel I. Standardization methods. In Encyclopedia of Biostatistics, Vol. 7, 2nd ed. Wiley, 2005; 5151–5163. 5. Logan BR, Nelson GO, Klein JP. Analyzing center specific outcomes in hematopoietic cell transplantation. Lifetime Data Analysis 2008; 14:389–404. DOI: 10.1007/s10985-008-9100-6. 6. Spiegelhalter D, Johnson CS, Bardsley M, Blunt I, Wood C, Grigg O. Statistical methods for healthcare regulation: rating, screening and surveillance. Journal of the Royal Statistical Society, Series A 2012; 175:1–47. DOI: 10.1111/j.1467-985X.2011.01010.x. 7. Keiding N. The method of expected number of deaths, 1786-1886-1986. International Statistical Review 1987; 55:1–20. 8. Wolfe RA, Gaylin DS, Port FK, Held PJ, Wood CL. Using USRDS generated mortality tables to compare local ESRD mortality rates to national rates. Kidney International 1992; 42:991–996.
Copyright © 2014 John Wiley & Sons, Ltd.
Statist. Med. 2014, 33 2048–2061
K. HE AND D. E. SCHAUBEL 9. Wolfe RA. The standardized morality ratio revisited: improvements, innovations, and limitations. American Journal of Kidney Diseases 1994; 24:290–297. 10. Dickinson DM, Shearon TH, O’Keefe J, Wong HH, Berg CL, Rosendale JD, Delmonico FL, Webb RL, Wolfe RA. SRTR center-specific reporting tools: posttransplant outcomes. American Journal of Transplantation 2006; 6:1198–1211. DOI: 10.1111/j.1600-6143.2006.01275.x. 11. Cox DR. Partial likelihood. Biometrika 1975; 62:269–276. 12. Breslow NE. Contribution to the discussion on the paper by D. R. Cox, regression and life table. Journal of the Royal Statistical Society, Series B 1972; 34:216–217. 13. He K. Semi-parametric and parametric methods for the analysis of multi-center survival data. Ph.D. Thesis, University of Michigan, Department of Biostatistics, Ann Arbor, 2012. 14. Kalbfleish JD, Wolfe RA. On monitoring outcomes of medical provider. Statistics in the Bioscience 2013; 5(2):286–302. 15. Ohlssen DI, Sharples LD, Spiegelhalter DJ. Flexible random-effects models using Bayesian semi-parametric models: application to institutional comparisons. Statistics in Medicine 2007; 26(9):2088–2112. DOI: 10.1002/sim.2666. 16. van der Vaart AW, Wellner JA. Weak Convergence and Empirical Processes. Springer: New York, 1996. 17. Kosorok MR. Introduction to Empirical Processes and Semiparametric Inference. Springer Series in Statistics, 2008. 18. Andersen PK, Borgan Ø, Gill RD, Keiding N. Statistical Models Based on Counting Processes. Springer Series in Statistics: New York, 1993. 19. van der Vaart AW. Asymptotic Statistics. Cambridge, 1998.
2061
Copyright © 2014 John Wiley & Sons, Ltd.
Statist. Med. 2014, 33 2048–2061
Accepted 20 December 2013
Published online 17 January 2014 in Wiley Online Library
(wileyonlinelibrary.com) DOI: 10.1002/sim.6089
Methods for comparing center-specific survival outcomes using direct standardization Kevin He* † and Douglas E. Schaubel The evaluation of center-specific outcomes is often through survival analysis methods. Such evaluations must account for differences in the distribution of patient characteristics across centers. In the context of censored event times, it is also important that the measure chosen to evaluate centers not be influenced by imbalances in the center-specific censoring distributions. The practice of using center indicators in a hazard regression model is often invalid, inconvenient, or undesirable to carry out. We propose a semiparametric version of the standardized rate ratio (SRR) useful for the evaluation of centers with respect to a right-censored event time. The SRR for center j can be interpreted as the ratio of the expected number of deaths in the total population (if the total population were in fact subject to the center j mortality hazard) to the observed number of events. The proposed measure is not affected by differences in center-specific covariate or censoring distributions. Asymptotic properties of the proposed estimators are derived, with finite-sample properties examined through simulation studies. The proposed methods are applied to national kidney transplant data. Copyright © 2014 John Wiley & Sons, Ltd. Keywords:
center effect; Cox regression; survival analysis; standardized rate ratio; stratification
1. Introduction In many situations, interest lies in the comparison of survival outcomes by center (e.g., treatment facility, hospital, or other entity serving as healthcare provider). Center-specific evaluations can be carried out on a regular basis (e.g., annually) to identity centers with poor performance. Alternatively, a retrospective evaluation over a longer period could be used to identify centers with exceptionally good results, with the goal of specifying best practices. In other cases, comparisons across centers may be an interesting secondary analysis; for example, in a multi-center study to evaluate the impact of a specific treatment on mortality. An accurate comparison of center-specific survival outcomes needs to account for imbalances in risk factor distributions among centers. For instance, the inclusion of high-risk patients by a given center can make that center’s survival appear substandard. In addition, in the context of survival times subject to censoring, the same phenomenon can occur because of differences in center-specific censoring distributions. In this report, we propose a semiparametric version of direct standardization, suitable for mortality comparisons by center. The proposed approach involves first fitting a Cox [1] regression model, stratified by center. The regression model is not used to directly estimate center effects, but rather to ensure that adjustment covariate effects are not confounded by center (which could occur in the absence of adjustment for center). For each center, j , the standardized rate ratio (SRRj ) is then computed as the ratio of expected to observed numbers of deaths, where ‘observed’ refers to the total number of deaths across all centers (j D 1; : : : ; J ) and ‘expected’ represents the number of total deaths estimated to occur if in fact all centers had mortality hazard equal to that of center j . Because of the use of direct standardization, the fSRR1 ; SRR2 ; : : : ; SRRJ g can be compared (and validly ordered) because the same
2048
Department of Biostatistics, School of Public Health, University of Michigan, Ann Arbor, MI 48109-2029, U.S.A. *Correspondence to: Kevin He, Department of Biostatistics, School of Public Health, University of Michigan, Ann Arbor, MI 48109-2029, U.S.A. † E-mail: [email protected]
Copyright © 2014 John Wiley & Sons, Ltd.
Statist. Med. 2014, 33 2048–2061
K. HE AND D. E. SCHAUBEL
covariate and censoring distribution is applied to each, that is, the total study population serves as the standard. This is in contrast to indirectly standardized measures, such as the standardized mortality ratio (SMR). The motivating example for this study involves the evaluation of center-specific post-transplant mortality for kidney transplant patients. Examples of factors known to strongly affect the post-transplant mortality hazard include age, primary renal diagnosis, and pre-transplant time on dialysis; each of which may differ in distribution by center. Centers with mortality significantly greater than the national average may be subject to various degrees of intervention, including site visits and perhaps de-accreditation. Given the high stakes of such evaluations, it is important that the statistical methods used for identifying outlying centers be accurate. A commonly used measure for evaluating center-specific survival is the SMR, defined as the ratio of the observed number of deaths at a given center to the number expected if the center had mortality equal to the population average. The SMR is a tool familiar to fields such as epidemiology (for example, [2–5]). The comparison of observed and expected outcomes is also commonly used for healthcare regulation (for example, [6, 7]). In the context of renal research, Wolfe et al. [8, 9] calculated SMRs among kidney transplant patients using mortality tables published by the United States Renal Data Systems. To relax the assumption of known standard mortality rates, Dickinson et al. [10] studied a semiparametric SMR based on Cox regression. Limitations of the SMR include its use of indirect standardization; each center’s SMR is essentially adjusted to a different (center-specific) covariate and censoring distribution. Centers cannot be rank ordered on the basis of indirect standardization, because two centers with equal covariate-specific mortality hazards could have different SMRs, merely due to differences in their respective covariate or censoring distributions. Therefore, although SMRs are potentially useful for internal evaluation (e.g., for centers to evaluate themselves or for a governing body to evaluate this center’s mortality comparing to that expected at the national level), they are less useful in the work of external evaluation (e.g., for surgeons and patients to compare center-specific results in the same region), because comparisons of center-specific results would play at least some role. We provide additional commentary on the SMR in Section 5. Earlier, we note limitations in the SMR approaches. Naturally, this technique has its role in analysis contrasting centers. However, noting room for improvement, we propose semiparametric version of direct standardization useful for survival analysis of centers. Direct standardization is also a commonly used approach in comparisons of mortality, usually through a measure termed the standardized rate ratio (SRR) or comparative mortality figure (CMF). One can express the SRR as the ratio of expected to observed numbers of deaths in the whole study population; the numerator of the SRR represents the expected number of deaths if all patients were treated at the given center, while the denominator equals the total observed number of deaths in the study population. Breslow and Day [3] and Hazel [4] compared the SRR with SMR in the framework of person-year methods. In particular, the main drawback of SMR (with respect to comparisons across centers) is not inherent to the SRR, because the same standard population is applied to all centers. Hence, center-specific SRRs are directly comparable. The SRR has a long history in fields such as epidemiology, and there are many settings in which direct standardization is appropriate. The Cox model has dominated applications involving regression analysis of censored data since its development. Thus, the use of Cox regression to estimate directly standardized center effects is a natural choice. The main contribution of this report is to formalize procedures for the Cox regression-based SRR, including rigorous derivation of asymptotic properties, simulation studies, and detailed comparisons to alternative approaches.
2. Methods
Copyright © 2014 John Wiley & Sons, Ltd.
Statist. Med. 2014, 33 2048–2061
2049
First, we provide the notation to be used in this article. Let Ti and Ci represent the survival and censoring time, respectively, for the ith patient, where i D 1; : : : ; n. Let J be the number of centers. The P total number of subjects is denoted by n D JjD1 nj , where nj is the number of subjects in center j . Observation times are denoted by Xi D Ti ^ Ci , with at-risk indicator Yi .t / D I.Xi > t /, where a ^ b D minfa; bg and I.A/ is an indicator function taking the value 1 when condition A holds and 0 otherwise. The observed death indicators are denoted by i D I.Ti 6 Ci /, and the death counting process is defined as Ni .t / D i I.Xi 6 t /. Let Gi denote the center for subject i and set Gij D I.Gi D j /. Correspondingly, we set Yij .t / D Yi .t /Gij and Nij .t / D Ni .t /Gij . The observed data consist of n independent vectors, .Xi ; i ; Gi ; Zi /, where Zi is a vector of adjustment covariates.
K. HE AND D. E. SCHAUBEL
The assumed center-stratified Cox model can be formulated as ˚ ij .t / .t jZi ; Gi D j / D 0j .t / exp ˇ0T Zi ;
(1)
where 0j .t / is an unspecified center-specific baseline hazard function and ˇ0 is a parameter vector. The partial likelihood estimator [11] of ˇ0 is denoted by ˇO and is given by the solution to U.ˇ/ D 0, where U.ˇ/ D
n Z J X X j D1 i D1
˚
Zij Z j .uI ˇ/ dNij .u/
0
P Sj.1/ .uI ˇ/ and Sj.d / .uI ˇ/ D n1 niD1 Yij .u/Zi˝d expfˇ T Zi g for Rt O 0j t I ˇO , where d D 0; 1; 2: The Breslow estimator [12] of ƒ0j .t / D 0 0j .u/du is then given by ƒ 1
with Z j .uI ˇ/ D Sj.0/ .uI ˇ/
n
X O 0j .t I ˇ/ D 1 ƒ n
Z
t 0
i D1
dNij .u/ Sj.0/ .uI ˇ/
:
2.1. Indirect standardization: standardized mortality ratio We begin by introducing an alternative to the proposed measure. The SMR for center j is calculated as the ratio of observed to expected numbers of events
1
Oj .t / ; Ej .t /
S MRj .t / D where the numerator is given by Oj .t / D computed as Ej .t / D
n Z X i D1
t
Pn
i D1
(2)
Nij .t / and the expected number of events is
n o O 0 uI ˇO ; Yij .u/ exp ˇO T Zi d ƒ
0
O 0 t I ˇO representing an estimator of the average cumulative baseline hazard with ƒ n
X O D1 O 0 .t I ˇ/ ƒ n i D1
Z
t 0
dNi .u/ S .0/ uI ˇO
(3)
and ˇO is based on model (1). The estimator given in (2) was developed in [13], where part of the innovation was the proposed use of the center-stratified model to estimate ˇ0 . Intuition would suggest computing Ej .t / using an unstratified model. However, estimation of ˇ0 in the absence of adjustment for center effects may produce a substantially biased estimate due to confounding by center or merely due to center affecting the hazard function, with or without confounding (i.e., due to the nonlinear link function). The calculation of S MRj .t / involves two stages. A stratified Cox model (1)) is fitted in the (model O 0 t I ˇO , and then computed in first stage, with the population average cumulative baseline hazard, ƒ
1
the second stage (e.g., through a unstratified Cox model with no covariates) using ˇO 0 Zi as an offset. At the second stage, Oj .t /, Ej .t / are calculated. 2.2. Direct standardization: standardized rate ratio
2050
On the basis of a form of indirect standardization, the SMR described in Section 2.1 can be viewed as a weighted ratio of center-specific cumulative hazards, with weight functions based on center-specific Sj.0/ .t I ˇ/. These weight functions have an obvious disadvantage: they involve center-specific censoring and covariate distributions, which can differ considerably across centers. To rule out the possibility that differences among centers are due merely to different censoring and/or covariate distributions, the weight function should be specified such that differences among centers with respect to the resulting Copyright © 2014 John Wiley & Sons, Ltd.
Statist. Med. 2014, 33 2048–2061
K. HE AND D. E. SCHAUBEL
measure are a function only of corresponding differences in center-specific hazards. Motivated by such considerations, we propose an alternative method, referred to as the SRR, which can be interpreted as a semiparametric version of direct standardization. The proposed SRR is computed, for center j , as
1
SRRj .t / D with O.t / D
PJ
j D1
Ej .t / ; O.t /
(4)
Oj .t / being the total observed number of deaths across all centers and Ej .t / D
n Z J X X `D1 i D1
t
n o O 0j .u/ Yi ` .u/ exp ˇO T Zi d ƒ
(5)
0
representing the expected number of total deaths if all centers had mortality hazard equal to that of center j . Similar to the SMR, the SRR is easily interpreted and is well understood by clinical investigators. The SRR also involves a ratio of observed and expected numbers of deaths. However, the ‘expected’ component is in the SRR’s numerator, while the ‘observed’ count is in the denominator. With respect to interpretation, SRRj > 1 indicates that center j has a greater mortality rate than the overall average. Note that, although SRRj .t / also involves the censoring and covariate distributions, the same weight function is applied across all centers, thus factoring out the impact of imbalances in center-specific censoring and covariate distributions. The proposed measures are desirable in this light, because their center-specific limiting values would differ only due to corresponding differences in centerspecific hazards. 2.3. Asymptotic properties We summarize the asymptotic properties of the proposed SRR with the following theorem; we outline the proof in Appendix A.
1
Theorem 1 Under the regularity conditions listed in Appendix A, SRRj .t / converges in probability to SRRj .t / uniformly in t 2 Œ0; , where R t .0/ `D1 0 s` .uI ˇ0 /dƒ0j .u/ PJ R t .0/ `D1 0 s` .uI ˇ0 /dƒ0` .u/
PJ SRRj .t / D
1
o n 1 and n 2 SRRj .t / SRRj .t / converges weakly to a zero-mean Gaussian process with variance function j .t / D EŒij .t I ˇ0 /2 , where Z
t
ij .t I ˇ/ D w.t I ˇ/ 0
s .0/ .uI ˇ/ sj.0/ .uI ˇ/
dMij .uI ˇ/ C w.t I ˇ/
Z tn
o Yi .u/ exp.ˇ T Zi / s .0/ .uI ˇ/ dƒ0j .u/
0
(6) w.t I ˇ/2
Z
t
s .0/ .uI ˇ/dƒ0j .u/ 0
Z C w.t I ˇ/ 0
J Z t n X `D1
0
t
rjT .uI ˇ/dƒ0j .u/.ˇ/1
o dNi ` .u/ s`.0/ .uI ˇ/dƒ0` .u/
Z
(7)
fZi ` ´` .uI ˇ/g dMij .uI ˇ/
(8)
0
with w.t I ˇ/ D
( J Z X
Copyright © 2014 John Wiley & Sons, Ltd.
0
)1 s`.0/ .uI ˇ/dƒ0` .u/
2051
rj .uI ˇ/ D s
`D1 .1/
t
.uI ˇ/ s .0/ .uI ˇ/ ´j .uI ˇ/: Statist. Med. 2014, 33 2048–2061
K. HE AND D. E. SCHAUBEL
O D n1 Pn Oij .t I ˇ/ O 2 , with Oij .t I ˇ/ O The variance function can be consistently estimated by O j2 .t I ˇ/ i D1 obtained by replacing limiting values in ij .t I ˇ/ with their empirical counterparts. With respect to this variance formula, a potentially time-saving strategy is to treat the estimators of ˇ0 as constants and hence ignore their variability. A justification for such simplification applies when the total study population is very large, such that ˇO has little variability. The resulting variance estimator is given by Z t .0/ Z tn o s .uI ˇ/ R T .0/ ij .t I ˇ/ D w.t I ˇ/ .uI ˇ/ C w.t I ˇ/ .u/ exp.ˇ Z / s .uI ˇ/ dƒ0j .u/ Y dM ij i i .0/ 0 s 0 j .uI ˇ/ (9) 2
Z
t
w.t I ˇ/
s
.0/
.uI ˇ/dƒ0j .u/
0
J Z t n X `D1
0
o dNi ` .u/ s`.0/ .uI ˇ/dƒ0` .u/ ;
(10)
obtained by removing the line (8) in the formula of ij .t I ˇ/. The quantity OijR t I ˇO is obtained by replacing limiting values in ijR .t I ˇ/ with their empirical counterparts.
3. Simulation We evaluate the finite-sample properties of the estimators described in Section 2 through a series of simulation studies. Death times were generated from the Weibull model, ij .t / D ˛j j t j 1 exp.ˇ T Zi / for i D 1; : : : ; nj and j D 1; : : : ; 10, where Zi D .Zi1 ; Zi 2 ; Zi 3 /T . We set ˇ T0 D .ˇ1 ; ˇ2 ; ˇ3 / D .0:02; 0:5; 0:2/. There are J D 10 centers. The number of subjects within each center varied under different scenarios. Censoring times were generated from either a Uniform distribution or an exponential distribution. In order to compare direct standardization with indirect standardization in the framework of semiparametric models, SMR and SRR were calculated at t D 1; t D 2, and t D 3. Each data configuration was replicated 1000 times. 3.1. Setting 1: center-independent hazards The first simulation setting considered the case where the hazard functions are equal across centers for all t 2 Œ0; . The censoring and covariate distributions were chosen to be center-independent. Specifically, censoring times were generated from a Uniform (0.5, 10) distribution; Zi1 followed a Bernoulli (0.5) distribution; Zi 2 followed a logistic distribution with probability dependent on Zi1 ; and Zi 3 came from a Normal distribution with constant variance 25 and mean dependent on Zi1 and Zi 2 (e.g., EŒZi 3 jZi1 ; Zi 2 D 50C0:2 Zi1 0:5 Zi 2 ). Under this setting, the hazard functions are equal across centers, such that the limiting values of SRR equal 1 for each center. Results at time t D 3 are displayed in Table I. For all centers, the average estimated SRR was very close to 1. The average asymptotic standard
O j .t/; center-independent hazards, covariate and censoring distributions; Table I. Simulation setting 1: SRR t D 3. Theorem 1 (6)–(8)
Theorem 1 (approx) (9), (10)
2052
Center
( j , ˛j )
True
Bias
ESD
ASE
CP
ASE
CP
1 2 3 4 5 6 7 8 9 10
(1, 0.2) (1, 0.2) (1, 0.2) (1, 0.2) (1, 0.2) (1, 0.2) (1, 0.2) (1, 0.2) (1, 0.2) (1, 0.2)
1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000
0.001 0.004 0:000 0:006 0:007 0:003 0.007 0:001 0.002 0.002
0.154 0.160 0.151 0.163 0.149 0.155 0.162 0.152 0.158 0.157
0.153 0.154 0.153 0.153 0.153 0.153 0.154 0.153 0.154 0.153
0.96 0.94 0.95 0.93 0.95 0.94 0.93 0.95 0.95 0.94
0.153 0.153 0.153 0.153 0.153 0.153 0.154 0.153 0.152 0.153
0.96 0.94 0.95 0.93 0.95 0.94 0.93 0.95 0.95 0.94
ESD, empirical standard deviation; ASE, asymptotic standard errors; CP, coverage probability.
Copyright © 2014 John Wiley & Sons, Ltd.
Statist. Med. 2014, 33 2048–2061
K. HE AND D. E. SCHAUBEL
errors were generally close to the empirical standard deviations, while the empirical coverage probabilities (CP) were generally consistent with the nominal value. This held for both the variance estimator derived in Theorem 1 (with corresponding asymptotic expansion given by (6)–(8) and its approximation given by (9) and (10). In fact, the two sets of asymptotic standard errors are virtually indistinguishable. Results were generally consistent across the observation time distribution (data not shown). It is worth noting that whether or not the covariate or censoring distributions were center dependent has no influence on the results in this setting. Correspondingly, similar results were found for the case of center-dependent covariate and censoring distribution.
3.2. Setting 2: center-dependent hazards; center-independent covariate and censoring distributions For the second set of simulations, different values of ˛j and j were used, such that the hazard functions increased with increasing center number (j D 1; : : : ; 10). The covariate distributions were chosen to be center independent and were generated from the same distributions from setting 1. Results are provided in Table II for various sample sizes and censoring percentages. The proposed SRR appears to be approximately unbiased, with CP close to 0.95 for both the standard error based on Theorem 1 and that based on the approximation. When the center-specific sample size is small (e.g., 25), the empirical CPs were
O j .t/; center-dependent hazards, center-independent covariate and Table II. Simulation setting 2: SRR censoring distributions; t D 3. Theorem 1 (6)–(8)
Theorem 1 (approx) (9), (10)
Censoring
Center
( j , ˛j )
True
Bias
ESD
ASE
CP
ASE
CP
100
20%
50
20%
100
40%
50
40%
125 100 75 50 25 125 100 50 25 15
20%
2 4 6 8 10 2 4 6 8 10 2 4 6 8 10 2 4 6 8 10 2 4 6 8 10 2 4 6 8 10
(0.85, 0.08) (0.95, 0.16) (1.05, 0.24) (1.15, 0.32) (1.25, 0.4) (0.85, 0.08) (0.95, 0.16) (1.05, 0.24) (1.15, 0.32) (1.25, 0.4) (0.85, 0.08) (0.95, 0.16) (1.05, 0.24) (1.15, 0.32) (1.25, 0.4) (0.85, 0.08) (0.95, 0.16) (1.05, 0.24) (1.15, 0.32) (1.25, 0.4) (0.85, 0.08) (0.95, 0.16) (1.05, 0.24) (1.15, 0.32) (1.25, 0.4) (0.85, 0.08) (0.95, 0.16) (1.05, 0.24) (1.15, 0.32) (1.25, 0.4)
0.406 0.797 1.180 1.549 1.912 0.406 0.797 1.180 1.549 1.912 0.406 0.797 1.180 1.549 1.912 0.406 0.797 1.180 1.549 1.912 0.529 1.043 1.550 2.025 2.535 0.601 1.118 1.773 2.324 2.900
0.004 0.003 0.001 0.006 0.005 0.002 0:008 0.007 0.001 0.026 0.008 0.003 0.007 0.014 0.029 0.009 0.011 0:007 0.019 0.043 0.006 0:004 0:002 0:015 0:004 0.001 0.007 0:001 0.045 0:062
0.096 0.132 0.174 0.206 0.230 0.135 0.183 0.252 0.290 0.344 0.097 0.139 0.172 0.214 0.243 0.137 0.204 0.249 0.302 0.357 0.110 0.169 0.256 0.384 0.651 0.122 0.190 0.375 0.636 1.034
0.093 0.134 0.169 0.204 0.238 0.130 0.187 0.238 0.285 0.338 0.097 0.138 0.174 0.208 0.244 0.136 0.195 0.242 0.294 0.345 0.106 0.171 0.255 0.381 0.632 0.118 0.192 0.359 0.621 0.914
0.94 0.95 0.94 0.95 0.96 0.93 0.95 0.94 0.95 0.95 0.94 0.95 0.95 0.94 0.95 0.94 0.94 0.94 0.94 0.95 0.94 0.95 0.95 0.93 0.93 0.94 0.96 0.93 0.93 0.91
0.093 0.134 0.169 0.203 0.237 0.130 0.186 0.238 0.284 0.337 0.097 0.138 0.173 0.208 0.243 0.135 0.195 0.241 0.294 0.344 0.106 0.171 0.255 0.380 0.630 0.118 0.192 0.358 0.619 0.910
0.94 0.95 0.94 0.95 0.96 0.93 0.95 0.94 0.95 0.94 0.94 0.95 0.95 0.94 0.95 0.94 0.94 0.94 0.94 0.95 0.94 0.95 0.95 0.93 0.92 0.94 0.96 0.93 0.93 0.90
40%
ESD, empirical standard deviation; ASE, asymptotic standard errors; CP, coverage probability.
Copyright © 2014 John Wiley & Sons, Ltd.
Statist. Med. 2014, 33 2048–2061
2053
nj
K. HE AND D. E. SCHAUBEL
slightly underestimated. Such results suggest that center-specific sample sizes play an important role in the proposed methods. In particular, the minimum sample size (across centers) needs to be reasonably large. Given this concern, centers of size less than 20 were eliminated from the real data analysis in Section 4. Collectively, simulation results from settings 1 and 2 indicate that the proposed method is quite accurate at and away from the null. 3.3. Setting 3: center-dependent hazards, center-dependent censoring distribution As mentioned previously, estimators based on indirect standardization may be misleading if either the censoring or the covariate distributions are center-dependent. To illustrate this point, we performed simulations with the following three conditions: (i) the hazard functions were center dependent, while center 2j 1 had exactly the same hazard function as that of center 2j for j D 1; : : : ; 5; (ii) the censoring distributions for centers 2j 1 and 2j were substantially different. For center 2j 1, the censoring times were generated from a uniform distribution such that the censoring mainly occurred in the later stages, while for center 2j , the censoring times were generated from a uniform distribution for which the censoring tended to occur in the early stages; and (iii) the covariate distributions (reused from setting 2) were center independent. We compared the SRR and SMR in Table III for setting 3. With respect to the true values, the limiting values of SRR2j .t / and SRR2j 1 .t / are equal, as one would hope. This is not the case for S MR2j .t / and S MR2j 1 .t /, the differences being due to differences in the censoring distributions. With respect to the estimators themselves, both SRRj .t / and S MRj .t / are approximately unbiased. Note that the bias of S MRj .t / was calculated as the difference between it and its own limiting value.
1 1
1
1
1
1
1
3.4. Setting 4: center-dependent hazards; center-dependent covariate distributions We also performed simulations under a setting in which the distribution of the covariate vector differed by center. The setup for hazard functions and censoring distribution was the same as in setting 2, while centers 2j and 2j 1 had substantially different covariate distributions. In center 2j 1, the covariate Zi1 followed a Bernoulli (0.2) distribution, Zi 2 followed a Bernoulli (0.8) distribution, and Zi 3 came from a Normal distribution with mean 30 and standard deviation 10. In center 2j , Zi1 followed a Bernoulli (0.8)
Table III. Simulation setting 3: centers 2j and 2j 1 have equal hazards but different censoring distributions; t D 3. ( j , ˛j )
2054
Measure
Center
SMRj .t /
1 2 3 4 5 6 7 8 9 10
(0.8, 0.04) (0.8, 0.04) (0.9, 0.12) (0.9, 0.12) (1, 0.2) (1, 0.2) (1.1, 0.3) (1.1, 0.3) (1.25, 0.4) (1.25, 0.4)
SRRj .t /
1 2 3 4 5 6 7 8 9 10
(0.8, 0.04) (0.8, 0.04) (0.9, 0.12) (0.9, 0.12) (1, 0.2) (1, 0.2) (1.1, 0.3) (1.1, 0.3) (1.25, 0.4) (1.25, 0.4)
True
Bias
ESD
0.227 0.211 0.658 0.637 1.050 1.051 1.500 1.547 1.881 2.003
0:003 0:002 0.001 0.003 0.003 0:000 0.002 0.004 0:002 0:008
0.089 0.068 0.162 0.121 0.208 0.157 0.239 0.201 0.279 0.219
0.221 0.221 0.645 0.645 1.044 1.044 1.530 1.530 1.996 1.996
0:003 0:002 0.002 0.004 0.002 0:002 0.001 0.004 0:004 0:008
0.088 0.073 0.162 0.124 0.216 0.158 0.277 0.202 0.351 0.230
ESD, empirical standard deviation.
Copyright © 2014 John Wiley & Sons, Ltd.
Statist. Med. 2014, 33 2048–2061
K. HE AND D. E. SCHAUBEL
Table IV. Simulation setting 4: centers 2j and 2j 1 have equal hazards but different covariate distributions; t D 3. Measure
Center
( j , ˛j )
True
SMRj .t /
1 2 3 4 5 6 7 8 9 10
(0.8, 0.04) (0.8, 0.04) (0.9, 0.12) (0.9, 0.12) (1, 0.2) (1, 0.2) (1.1, 0.3) (1.1, 0.3) (1.25, 0.4) (1.25, 0.4)
SRRj .t /
1 2 3 4 5 6 7 8 9 10
(0.8, 0.04) (0.8, 0.04) (0.9, 0.12) (0.9, 0.12) (1, 0.2) (1, 0.2) (1.1, 0.3) (1.1, 0.3) (1.25, 0.4) (1.25, 0.4)
Bias
ESD
0.319 0.331 0.931 0.919 1.484 1.307 2.070 1.602 2.551 1.725
0:003 0:002 0.001 0.003 0.003 0:000 0.002 0.004 0:002 0:008
0.089 0.068 0.162 0.121 0.208 0.157 0.239 0.201 0.279 0.219
0.345 0.345 0.929 0.929 1.386 1.386 1.908 1.908 2.329 2.329
0:002 0:000 0.005 0.003 0:007 0:000 0:007 0.013 0:006 0:001
0.120 0.041 0.212 0.105 0.272 0.186 0.344 0.330 0.412 0.455
ESD, empirical standard deviation.
distribution, Zi 2 followed a Bernoulli (0.2) distribution, and Zi 3 was derived from a Normal distribution with mean 50 and standard deviation 10. Results based on setting 4 are given in Table IV. Trends are similar to those from setting 3 but much more pronounced. The limiting values of SRR2j .t / are equal to those of SRR2j 1 .t /, as one would expect. Conversely, the true values of S MR2j .t / and S MR2j 1 .t / are different, with the differences being quite pronounced for j D 3, j D 4, and j D 5. Moreover, it appears that SMR7 .t / > SMR10 .t /, which is misleading in the sense that SRR7 .t / < SRR10 .t /.
1 1
1
1
4. Application
1
1
Copyright © 2014 John Wiley & Sons, Ltd.
Statist. Med. 2014, 33 2048–2061
2055
We applied the proposed methods to investigate the performance of transplant centers with respect to post kidney transplant survival. Data were obtained from the Scientific Registry of Transplant Recipients (SRTR) and submitted by members of the Organ Procurement and Transplantation Network. The SRTR database contains information on all wait-listed candidates, transplant recipients, and organ donors in the USA. Included in the analysis were adult patients (>18 years of age at transplant) who underwent deceased donor kidney transplantation between January 2000 and December 2008. Adjustment covariates in this study included age, race, gender, diagnosis, donation after cardiac death, expanded criteria donor, body mass index, dialysis time, indicator of previous kidney transplant, and cold ischemia time. These variables have face validity from a clinical perspective and are based on a list of covariates used in SRTR. Transplant centers with sample size 620 and patients who received a living-donor transplant were eliminated from additional analysis. The final sample size was then n D 74; 088 from J D 217 centers across the USA. Failure time (recorded in years) was defined as the time from transplantation to graft failure or death, whichever occurred first. Graft failure was considered to occur when the transplanted kidney ceased to function. Stratified Cox regression was employed to model the hazard function. The indirectly standardized estimator, S MRj , was calculated using SAS PROC PHREG with an offset. The proposed directly standardized estimator, SRRj , was computed using SAS IML. Figure 1a represents the pairwise comparisons of the SMRs and SRRs. Figure 1b shows the standard error of these two measures. Figure 1c compares
K. HE AND D. E. SCHAUBEL
(b) (Scatter plot : standard error of SMR and SRR)
0.8
SE of SMR
SMR 1.0
0.0
0.0
0.4
2.0
1.2
3.0
(a) (Scatter plot : SMR VS SRR)
0.0
0.5
1.0
1.5
2.0
2.5
3.0
0.0
0.2
SRR
0.4
0.6
0.8
1.0
1.2
SE of SRR
200 150 100 50 0
order of centers based on SRR
(c) (Scatter plot : orders of centers based on SMR and SRR)
0
50
100
150
200
order of centers based on SRR
Figure 1. Evaluation of J D 217 kidney transplant centers. SMR, standardized mortality ratio; SRR, standardized rate ratio; SE, standard error.
2056
the orders of centers based on SRRs and SMRs. As shown, there are some discrepancies between these two measures. We applied bootstrapped techniques to evaluate whether the change in center-specific orderings (SMR versus SRR) exceeded that attributable to only sampling variation. Specially, we calculated the distribution of center-specific SMR orderings from 100 bootstrapped samples. We then calculated the 95% confidence intervals of these orders. On the basis of the bootstrapped resamples, for 20 of 217 centers, the 95% confidence intervals based on SMR does not cover the order based on SRR from the original dataset. Among these 20 centers, 10 had SRR significantly different from the national average. Using the asymptotic normality of the proposed estimators, we constructed the point-wise confidence intervals for SRR at t D 5 years (Figure 2). Center numbers are re-ordered by values of SRR. A total of 38 centers had observed number of events significantly lower than the expected calculation based on the national average hazards, while 28 centers were significantly above the expected. It is clear that the hazard functions varied among centers. Table V presents the pairwise comparison of the numbers and percentages of ‘outlier’ centers identified by p-values corresponding to their SMRs and SRRs (using tests of H0 W SMRj D 1 and H0 W SRRj D 1, respectively). A total of six centers changed ‘memberships’ based on these two measures. Specifically, two centers were flagged to be significant based on SMR but flagged to be normal based on SRR; on the other hand, four centers were flagged to be significant based on SRR but flagged to be normal based on SMR. Through fitting a sequence of logistic regression models (rotating the center indicators as the response variates), it was revealed that, for these six centers, approximately half of the adjustment covariates had distributions significantly different from the remaining centers. In addition, through fitting a Cox model using censored as the event, three of the six centers in question were significantly predictive of the censoring hazard. In summary, we do observe some differences when comparing center effects estimated through direct versus indirect standardization, and the strongest examples of such discrepancies appear to be due to differences in the center-specific covariate and censoring distributions, consistent with the concepts described earlier in this report. Copyright © 2014 John Wiley & Sons, Ltd.
Statist. Med. 2014, 33 2048–2061
SRR(5)
0
1
2
3
4
5
K. HE AND D. E. SCHAUBEL
0
50
100
150
200
Center
Figure 2. Evaluation of J D 217 kidney transplant centers: point estimates and 95% confidence interval of O j at t D 5 years. SRR, standardized rate ratio. SRR
Table V. Number and percentage of centers giving significant results under SMR and SRR. SRR SMR Nonsignificant Significant Column sum
Nonsignificant
Significant
Row sum
147 (67.8%) 2 (0.9%) 149 (68.7%)
4 (1.8%) 64 (29.5%) 68 (31.3%)
151 (69.6%) 66 (30.4 %) 217 (100%)
SMR, standardized mortality ratio; SRR, standardized rate ratio.
5. Discussion
Copyright © 2014 John Wiley & Sons, Ltd.
Statist. Med. 2014, 33 2048–2061
2057
We propose semiparametric methods for estimating standardized rate ratios, as a means of evaluating center-specific mortality through direct standardization. Large-sample properties are derived and shown through simulation to be appropriate in finite samples. A computationally faster variance estimator is proposed for the SRR and is shown to work practically as well as the full version. Application of the methods demonstrates several significant differences among kidney transplant centers in the USA. There is some judgement required in deciding when to use indirect standardization and when to use direct standardization. Indirectly standardized estimators, such as SMR, provide a valid approach to evaluate how does a center’s mortality compare to that predicted at the population level for the kinds of patients at this center. However, it is important to emphasize that center-specific SMRs should not be compared with one another (a caution that applies to all indirectly standardized rates). The SRR, a directly standardized measure, does not share this drawback. The SRRs for two given centers will be unequal only because the center-specific mortality hazards differ; direct standardization accounts for imbalance with respect to center-specific covariate and censoring distributions. The proposed SRR shares the SMR’s ease of interpretation but rectifies its key disadvantages and, hence, is a more appropriate choice in settings where mortality comparisons across centers are an objective. The degree to which the SMR and SRR are different will depend on the application. In some settings, the two may not agree well, while in others, they may be quite similar. The only way to know with certainty if SMR and SRR are equal would be to calculate both measures, which would not be a desirable option in many cases. In settings where, for a particular center, the SMR and SRR were unequal, it would be very difficult to claim that the SMR was correct, for the several reasons documented previously. Given the high stakes of evaluations by regulatory bodies, and the fact that the credibility of such organizations depends in part on the accuracy of their evaluations, it would appear the preferred analysis is the one that is most likely to be accurate.
K. HE AND D. E. SCHAUBEL
The proposed SRR is computed using a stratified Cox model, which makes no assumptions about the functional form of the impact of center on the hazard function. The stratification by center plays a major role. For instance, the expected number of deaths considers each patient’s covariate vector, such that the regression parameter must be estimated consistently. Unless the death hazard is conditionally independent of center given the covariates, covariate effect estimates will generally be biased if based on a model with a nonlinear link function and no accounting for center. This is an issue for the SMR as well, particularly because the Cox version of the SMR has historically been computed using a commonbaseline (i.e., unstratified) model. He [13] proposed modifying the SMR through stratification, leading to the quantity we denoted by SMR and used in simulations (Section 3) for comparisons with the SRR. Such properties were demonstrated empirically by He [13]; the magnitude of the bias (in the case of unstratified Cox models) increases when covariates are also center-dependent. Random effects models are an option for contrasting centers. Moreover, a Bayesian formulation may be an attractive approach for center effect studies. The advantage for random effects model is that this approach would allow for the inclusion of even very small centers. In contrast, when the number of events for a center is small, the estimated center parameter from a fixed effect model may be unstable. However, Kalbfleish and Wolfe [14] compared the properties of a fixed effect model (FEM) and a random effect model (REM) for the purpose of profiling kidney dialysis facilities under various conditions. Essentially, the REM estimates are shrunk toward overall mean and hence reduce the reported variation of facility performance. Second, the FEM has the highest statistical power to identify exceptional facilities, for a given false positive rate; identifying such extreme facilities is usually a main objective of center evaluations. Another issue for REM is the potential confounding effects when the patient risks are correlated with center effects. These findings suggest that a simple REM method may not be good enough and more sophisticated approaches are necessary. Further discussion of such issues is provided by Ohlssen et al. [15], who develop a more flexible random effects model using Bayesian nonparametric methods, in order to remedy the influence of outlying centers to which basic random effects models are susceptible. The direct standardization methods derived in this report could be extended in several useful directions. Perhaps most notably, it is often of interest to evaluate center effects in settings where the event of interest is recurrent (e.g., hospitalizations and infections). Furthermore, it would also be useful to develop methods based on direct standardization that can accommodate competing risks or dependent censoring. Direct standardization could also be applied to compare center-specific survival probability and restricted mean lifetime.
Appendix A To derive the large-sample properties for the SRR, we impose the following regularity conditions under the stratified Cox model: (a) .Xi ; i ; Gi ; Zi / are independent and identically distributed random vectors. (b) P .Xi > / > 0 where is a pre-specified time point. (c) Zi k have bounded total variation, that is, jZi k j < for all i D 1; : : : n and k D 1; : : : ; p, where
is R a constant and Zi k is the kth component of Zi . (d) 0 0j .t /dt < 1. (e) Continuity of the following functions: sj.1/ .t I ˇ/ D
@ .0/ @2 s .0/ .t I ˇ/ sj .t I ˇ/; sj.2/ .t I ˇ/ D @ˇ @ˇ@ˇ T j
and sj.0/ .t I ˇ/, where sj.d / .t I ˇ/ is the limiting value of Sj.d / .t I ˇ/ for d D 0; 1; 2, with sj.1/ .t I ˇ/ and sj.2/ .t I ˇ/ bounded and sj.0/ .t I ˇ/ bounded away from 0 for t 2 Œ0; . (f) Positive-definiteness of the matrix j .ˇ/: Z j .ˇ/ D vj .t I ˇ/sj.0/ .t I ˇ/0j .t /dt; 0
2058
vj .t I ˇ/ D
Copyright © 2014 John Wiley & Sons, Ltd.
sj.2/ .t I ˇ/ sj.0/ .t I ˇ/
´j .t I ˇ/˝2 ;
Statist. Med. 2014, 33 2048–2061
K. HE AND D. E. SCHAUBEL 1 .1/ sj .t I ˇ/
where ´j .t I ˇ/ D sj.0/ .t I ˇ/ (g) P .Gij D 1jZi / > 0.
is the limiting value of Z j .t I ˇ/.
Condition (a) is employed in the derivation of the weak convergence. Condition (b) is a standard identifiability requirement. Condition (c) leads to the boundedness of several quantities and is applicable in most practical applications. Conditions (d) and (e) are not essential but simplify our proofs. With respect to condition (g), the selection probability given covariates is nonzero for all centers. This condition guarantees that the sample size nj of each center goes to 1 as the total sample size n goes to 1. P O j .t / We first show that SRR ! SRRj .t / uniformly for t 2 Œ0; . The triangle inequality leads to Z t Z t .0/ .0/ O O O S` uI ˇ d ƒ0j uI ˇ s` .uI ˇ0 /dƒ0j .u/ 0 0 Z t Z t .0/ .0/ O O O O 0j uI ˇ O 0j uI ˇ 6 S` uI ˇ d ƒ s` .uI ˇ0 /d ƒ 0 0 Z t Z t O 0j uI ˇO C s`.0/ .uI ˇ0 /d ƒ s`.0/ .uI ˇ0 /dƒ0j .u/ 0
(A.1) (A.2)
0
P
! 0 uniformly in t , recall that To show that .A:1/ P S`.0/ .t I ˇ/ D n1 niD1 Yi ` .t / exp.ˇ T Zi / D Pn ŒI.X > t /I.G D `/ exp.ˇ T Z/, where Pn is P the empirical measure; that is, Pn ŒI.X > t /I.G D `/ exp.ˇ T Z/ D n1 niD1 Yi ` .t / exp.ˇ T Z/. The collection of all cells, Œt; 1/, in the real line is a VC class of index 2 and, hence, satisfies the entropy conditions for the Glivenko–Cantelli theorem [16]. The boundedness conditions ensures that fI.X > t /I.G D `/ exp.ˇ T Z/; t 2 Œ0; g belong to some Glivenko–Cantelli class; that is, a:s: P O O uniformly in t 2 Œ0; . Next, for ˇ, such that ˇO S`.0/ .t I ˇ/ ! s`.0/ .uI ˇ/ ! ˇ0 (e.g, [17, 18]), O 0 .t /; t 2 Œ0; g and an application of the dominant convergence thethe bounded condition of fƒ P
P
orem entails that .A:1/ ! 0 uniformly in t . Similarly, .A:2/ ! 0 uniformly in t . We already R t .0/ P R t .0/ O demonstrate that 0 S` .uI ˇ0 /d ƒ0j .u/ ! 0 s` .uI ˇ0 /dƒ0j .u/ uniformly for t 2 Œ0; . SimiR t .0/ P R t .0/ O 0` .u/ larly, 0 S` .uI ˇ0 /d ƒ ! 0 s` .uI ˇ0 /dƒ0` .u/ uniformly for t 2 Œ0; . The monotonicity o1 P nR o1 nR t t O 0` .u/ ! 0 s`.0/ .uI ˇ0 /dƒ0` .u/ and boundedness conditions ensure that 0 S`.0/ .uI ˇ0 /d ƒ P O j .t / uniformly for t 2 Œ0; . Therefore, we have that SRR ! SRRj .t / uniformly for t 2 Œ0; . To prove weak convergence, we use the following decomposition:
1
o n 1 n 2 SRRj .t / SRRj .t / ( R t .0/ ) R t .0/ O ƒ O O O 0j .uI ˇ/ S` .uI ˇ/d s` .uI ˇ0 /dƒ0j .uI ˇ/ 1 0 0 D n2 R t .0/ R t .0/ O ƒ O O 0` .u/ O S .uI ˇ/d 0 ` 0 S` .uI ˇ0 /d ƒ0` .uI ˇ/ ( R t .0/ ) R t .0/ s` .uI ˇ0 /dƒ0j .u/ s .uI ˇ0 /dƒ0j .u/ 1 R0t `.0/ : C n 2 R t 0 .0/ O ƒ O O 0` .uI ˇ/ S .uI ˇ/d s .uI ˇ0 /dƒ0` .u/ 0
0
`
(A.3)
(A.4)
`
First, we have that ) Z t .0/ .uI ˇ / s dN .u/ 0 ij i D1 O S .0/ .uI ˇ/ s .0/ .uI ˇ0 / `.0/ dƒ0` .u/ C op .1/ .A:3/ D w.t; ˇ0 / n O s .uI ˇ0 / 0 0 nSj.0/ .uI ˇ/ " n Z t .0/ o X s .uI ˇ0 / n .0/ 1 2 .uI ˇ / s .uI ˇ /dƒ .u/ D w.t I ˇ0 / n dN ij 0 0 0j j .0/ i D1 0 sj .uI ˇ0 / Z t 1 C rjT .uI ˇ/dƒ0j .u/ n 2 .ˇO ˇ/ C op .1/; 1 2
(Z
Pn
t
Copyright © 2014 John Wiley & Sons, Ltd.
Statist. Med. 2014, 33 2048–2061
2059
0
K. HE AND D. E. SCHAUBEL
where w.t I ˇ/ and rj .uI ˇ/ are defined in Section 2. Note that the second equality of the aforementioned argument is obtained through the functional delta method and Lemma 19.24 of [19]. On the basis of previously established empirical process theory for the Cox model (e.g., [17]), we have that ( J ) " J X X 1 1 n 2 .ˇO ˇ0 / D ` .ˇ0 / Gn ffZ ´` .X I ˇ0 /gI.G D `/ `D1
Z
0
`D1
fZ ´` .uI ˇ0 /g exp ˇ0T Z I.X > u/I.G D `/dƒ0` .u/ C op .1/;
where the op .1/ is uniform in t , and Gn is the empirical process defined by Gn f D Through the functional delta method,
p n.Pn P/f .
Rt
J X s .0/ .uI ˇ0 /dƒ0j .u/ .A:4/ D n Gn ŒN` .t / C op .1/: o2 PJ R t .0/ `D1 s .uI ˇ /dƒ .u/ 0 0` `D1 0 ` 0
Combining (A.3) and (A.4), o n 1 n 2 SRRj .t / SRRj .t / (Z ( J " ) Z t t .0/ X s .uI ˇ0 / D Gn w.t I ˇ0 / rjT .uI ˇ0 /dƒ0j .u/ ` .ˇ0 /1 dNj .u/ C .0/ 0 s 0 .uI ˇ / 0 j `D1
) Z J X T fZ ´` .uI ˇ0 /g exp ˇ0 Z I.X > u/I.G D `/dƒ0` .u/ fZ ´` .X I ˇ0 /g I.G D `/
1
0
`D1
Rt
3
.0/
s .uI ˇ0 /dƒ0j .u/ 7 o2 I.X 6 t /5 C op .1/I R t .0/ `D1 0 s` .uI ˇ0 /dƒ0` .u/
n PJ
0
Because the VC classes with finite index satisfy the entropy conditions for the Donsker theorem [16], fI.X > t /; t 2 Œ0; g and fI.X 6 t /; t 2 Œ0; g belong to some Donsker classes. The same holds for the bounded monotone stochastic process fN.t /; t 2 Œ0; g. Finally, the class of functions of Lipschitz transformations of Donsker classes is Donsker. Therefore, with the various bounded conditions and by applying the Donsker theorem, Theorem 1 follows.
Acknowledgements This work was supported in part by the National Institutes of Health grant 5R01-DK070869 and a grant from the Michigan Institute for Clinical and Health Research (MICHR). The authors thank the Scientific Registry of Transplant Recipients (SRTR) for allowing us to access the organ failure database.
References
2060
1. Cox DR. Regression models and life tables (with Discussion). Journal of the Royal Statistical Society, Series B 1972; 34:187–200. 2. Berry G. The analysis of mortality by the subject-years method. Biometrics 1983; 39:173–184. 3. Breslow NE, Day NE. The standardized mortality ratio. In Biostatistics: Statistics in Biomedical, Public Health and Environmental Sciences, Sen PK (ed.), The Bernard G. Greenberg Volume, 1985; 55–74. 4. Hazel I. Standardization methods. In Encyclopedia of Biostatistics, Vol. 7, 2nd ed. Wiley, 2005; 5151–5163. 5. Logan BR, Nelson GO, Klein JP. Analyzing center specific outcomes in hematopoietic cell transplantation. Lifetime Data Analysis 2008; 14:389–404. DOI: 10.1007/s10985-008-9100-6. 6. Spiegelhalter D, Johnson CS, Bardsley M, Blunt I, Wood C, Grigg O. Statistical methods for healthcare regulation: rating, screening and surveillance. Journal of the Royal Statistical Society, Series A 2012; 175:1–47. DOI: 10.1111/j.1467-985X.2011.01010.x. 7. Keiding N. The method of expected number of deaths, 1786-1886-1986. International Statistical Review 1987; 55:1–20. 8. Wolfe RA, Gaylin DS, Port FK, Held PJ, Wood CL. Using USRDS generated mortality tables to compare local ESRD mortality rates to national rates. Kidney International 1992; 42:991–996.
Copyright © 2014 John Wiley & Sons, Ltd.
Statist. Med. 2014, 33 2048–2061
K. HE AND D. E. SCHAUBEL 9. Wolfe RA. The standardized morality ratio revisited: improvements, innovations, and limitations. American Journal of Kidney Diseases 1994; 24:290–297. 10. Dickinson DM, Shearon TH, O’Keefe J, Wong HH, Berg CL, Rosendale JD, Delmonico FL, Webb RL, Wolfe RA. SRTR center-specific reporting tools: posttransplant outcomes. American Journal of Transplantation 2006; 6:1198–1211. DOI: 10.1111/j.1600-6143.2006.01275.x. 11. Cox DR. Partial likelihood. Biometrika 1975; 62:269–276. 12. Breslow NE. Contribution to the discussion on the paper by D. R. Cox, regression and life table. Journal of the Royal Statistical Society, Series B 1972; 34:216–217. 13. He K. Semi-parametric and parametric methods for the analysis of multi-center survival data. Ph.D. Thesis, University of Michigan, Department of Biostatistics, Ann Arbor, 2012. 14. Kalbfleish JD, Wolfe RA. On monitoring outcomes of medical provider. Statistics in the Bioscience 2013; 5(2):286–302. 15. Ohlssen DI, Sharples LD, Spiegelhalter DJ. Flexible random-effects models using Bayesian semi-parametric models: application to institutional comparisons. Statistics in Medicine 2007; 26(9):2088–2112. DOI: 10.1002/sim.2666. 16. van der Vaart AW, Wellner JA. Weak Convergence and Empirical Processes. Springer: New York, 1996. 17. Kosorok MR. Introduction to Empirical Processes and Semiparametric Inference. Springer Series in Statistics, 2008. 18. Andersen PK, Borgan Ø, Gill RD, Keiding N. Statistical Models Based on Counting Processes. Springer Series in Statistics: New York, 1993. 19. van der Vaart AW. Asymptotic Statistics. Cambridge, 1998.
2061
Copyright © 2014 John Wiley & Sons, Ltd.
Statist. Med. 2014, 33 2048–2061
