Simultaneous Confidence Intervals for Ranks With Application to Ranking Institutions

TL;DR

This paper introduces a new conservative method based on Tukey's HSD for constructing simultaneous confidence intervals for ranks of institutions, addressing limitations of existing Monte-Carlo approaches, especially with larger datasets.

Contribution

A novel, scalable method using Tukey's HSD for confidence intervals of ranks, improving over Monte-Carlo methods in large-scale applications.

Findings

New method is conservative when no ties exist.

Monte-Carlo method is anticonservative and unscalable for large N.

Both methods perform similarly after rescaling to the nominal confidence level.

Abstract

When a ranking of institutions such as medical centers or universities is based on an indicator provided with a standard error, confidence intervals should be calculated to assess the quality of these ranks. We consider the problem of constructing simultaneous confidence intervals for the ranks of means based on an observed sample. For this aim, the only available method from the literature uses Monte-Carlo simulations and is highly anticonservative especially when the means are close to each other or have ties. We present a novel method based on Tukey's honest significant difference test (HSD). Our new method is on the contrary conservative when there are no ties. By properly rescaling these two methods to the nominal confidence level, they surprisingly perform very similarly. The Monte-Carlo method is however unscalable when the number of institutions is large than 30 to 50 and stays thus anticonservative. We provide extensive simulations to support our claims and the two methods are compared in terms of their simultaneous coverage and their efficiency. We provide a data analysis for 64 hospitals in the Netherlands and compare both methods. Software for our new methods is available online in package ICRanks downloadable from CRAN. Supplementary materials include supplementary R code for the simulations and proofs of the propositions presented in this paper.

Figures

Tables

Table 1. Table 1 : Values of α ~ ~ 𝛼 \tilde{\alpha} necessary to rescale the coverage at worst case back to 1 − α 1 𝛼 1-\alpha .

	Rescaled coverage
	$95\%$			$90\%$			$80\%$
	Tukey	Zhang		Tukey	Zhang		Tukey	Zhang
$n=10$	0.158	$6.5\times {10}^{-4}$		0.285	0.0015		0.467	0.006
$n=30$	0.303	$9.8\times {10}^{-6}$		0.491	$4.6\times {10}^{-5}$		0.693	$4\times {10}^{-5}$
$n=50$	0.418			0.574	$7\times {10}^{-6}$		0.778	$3.1\times {10}^{-5}$
$n=100$	0.545			0.725			0.893	$5\times {10}^{-6}$

Table 2. Table 2 : Coverage probability and the efficiency when τ = 0.5 𝜏 0.5 \tau=0.5 and α = 0.1 𝛼 0.1 \alpha=0.1 before and after rescaling to the worst case.

	Coverage						$1-{\hat{R}}_n\left(\alpha \right)$
	not rescaled			rescaled			not rescaled			rescaled
	Tukey	Zhang		Tukey	Zhang		Tukey	Zhang		Tukey	Zhang
$n=10$	0.998	0.468		0.961	0.976		0.990	0.789		0.971	0.977
$n=30$	1.000	0.027		0.978	0.987		0.998	0.740		0.990	0.991
$n=50$	0.997	0.000		0.976	0.984		0.999	0.726		0.994	0.995

Table 3. Table 3 : Coverage probability and the rankability when τ = 1 𝜏 1 \tau=1 and α = 0.1 𝛼 0.1 \alpha=0.1 before and after rescaling to the worst case.

	Coverage						$1-{\hat{R}}_n\left(\alpha \right)$
	not rescaled			rescaled			not rescaled			rescaled
	Tukey	Zhang		Tukey	Zhang		Tukey	Zhang		Tukey	Zhang
$n=10$	0.996	0.603		0.972	0.994		0.959	0.698		0.916	0.935
$n=30$	1.000	0.088		0.993	0.996		0.987	0.651		0.957	0.967
$n=50$	0.999	0.016		0.996	0.998		0.992	0.640		0.970	0.976

Table 4. Table 4 : Coverage probability and the rankability when τ = 2 𝜏 2 \tau=2 and α = 0.1 𝛼 0.1 \alpha=0.1 before and after rescaling to the worst case.

	Coverage						$1-{\hat{R}}_n\left(\alpha \right)$
	not rescaled			rescaled			not rescaled			rescaled
	Tukey	Zhang		Tukey	Zhang		Tukey	Zhang		Tukey	Zhang
$n=10$	0.997	0.814		0.989	0.997		0.811	0.529		0.734	0.788
$n=30$	0.998	0.262		0.988	0.996		0.888	0.479		0.802	0.844
$n=50$	1.000	0.065		0.997	1.000		0.911	0.468		0.831	0.867

Table 5. Table 5 : Rescaled significance level α ~ ~ 𝛼 \tilde{\alpha} in order to get an actual simultaneous coverage equal to 90 % percent 90 90\% at the worst case configuration for the dataset of Dutch hospitals when considering the mortality as the primary outcome.

	Tukey	Zhang
$n=10$	0.246	0.0003
$n=30$	0.432
$n=61$	0.583

Equations

$$y_i\sim \mathcal{N}({\mu }_i,{\sigma }_i^2),\text{for }i\in \{1,\cdots ,n\},$$

$$l_i=1+\sum _{j\ne i}\text{mathds}{1}_{{\mu }_j<{\mu }_i}.$$

$$u_i=n-\sum _{j\ne i}\text{mathds}{1}_{{\mu }_j\ge {\mu }_i}.$$

$$\mathbb{P}\left([l_i,u_i]\subseteq [L_i,U_i],\forall i\in \{1,\cdots ,n\}\right)\ge 1-\alpha$$

$$\left\{\frac{\left|y_i-y_j\right|}{\sqrt{{\sigma }_i^2+{\sigma }_j^2}}>q_{1-\alpha }\right\}$$

$$\underset{i,j=1,\cdots ,n}{\max }\frac{|{\tilde{Y}}_i-{\tilde{Y}}_j|}{\sqrt{{\sigma }_i^2+{\sigma }_j^2}},$$

$$L_i=$$

$$U_i=$$

$${\mu }_A-{\mu }_B\in [-2,-1]$$

$${\mu }_B-{\mu }_A\in [1,2]$$

$${\mu }_C-{\mu }_A\in [2,3]$$

$${\mathbb{P}}_{\mu }\left(\forall i,r_i\in [L_i,U_i]\right)\ge {\mathbb{P}}_{\mu =0}\left(\forall i,r_i\in [L_i,U_i]\right).$$

$${\mathbb{P}}_{\mu }\left(\forall i,r_i\in [L_i(\alpha ),U_i(\alpha )]\right)={\beta }_{\mu }(\alpha )\ne 1-\alpha .$$

$${\mathbb{P}}_{\mu }\left(\forall i,r_i\in [L_i(\tilde{\alpha }),U_i(\tilde{\alpha })]\right)={\beta }_{\mu }(\tilde{\alpha })=1-\alpha .$$

$${\mathbb{P}}_{\mu =0}\left(\forall i,r_i\in [L_i(\tilde{\alpha }),U_i(\tilde{\alpha })]\right)={\beta }_0(\tilde{\alpha })=1-\alpha .$$

$${\mathbb{P}}_{\mu }\left(\forall i,r_i\in [L_i(\tilde{\alpha }),U_i(\tilde{\alpha })]\right)\ge {\mathbb{P}}_{\mu =0}\left(\forall i,r_i\in [L_i(\tilde{\alpha }),U_i(\tilde{\alpha })]\right)=1-\alpha .$$

$${\sigma }_1\le \cdots \le {\sigma }_n$$

$${\sigma }_1\le \cdots \le {\sigma }_{n/2},$$

$${\sigma }_1\ge \cdots \ge {\sigma }_{n/2},$$

$$R_n=1-\frac{1}{n(n-1)}\sum _{i=1}^n(u_i-l_i).$$

$$\mathbb{P}\left(\forall i=1,\cdots ,n[l_i,u_i]\subseteq [L_i,U_i]\right)\ge 1-\alpha .$$

$${\hat{R}}_n(\alpha )=1-\frac{1}{n(n-1)}\sum _{i=1}^n(U_i-L_i).$$

$$\mathbb{P}\left(R_n\ge {\hat{R}}_n(\alpha )\right)\ge 1-\alpha .$$

$$\mathbb{P}\left({\mu }_i-{\mu }_j\in \left[y_i-y_j\pm \sqrt{{\sigma }_i^2+{\sigma }_j^2}{q}_{1-\alpha }\right],\forall i,j\right)\ge 1-\alpha .$$

$$\mathbb{P}\left(\forall i,l_i\ge L_i,u_i\le U_i\right)\ge \mathbb{P}\left(\forall i\ne j,{\mu }_i-{\mu }_j\in [a_{i,j},b_{i,j}]\right)\ge 1-\alpha .$$

$$\left\{L_i\le 1,U_i\ge n,\forall i\right\}$$

$${\mu }_i-{\mu }_j\in \left[y_i-y_j\pm \sqrt{{\sigma }_i^2+{\sigma }_j^2}{q}_{1-\alpha }\right].$$

$$\{L_i=1+\#\{j\ne i:y_i\ge y_j-\sqrt{{\sigma }_i^2+{\sigma }_j^2}{q}_{1-\alpha }\}=1,U_i=n-\#\{j\ne i:y_i\le y_j+\sqrt{{\sigma }_i^2+{\sigma }_j^2}{q}_{1-\alpha }\}=n,\forall i\}.$$

$$\{\forall j\ne i,y_i\le y_j-\sqrt{{\sigma }_i^2+{\sigma }_j^2}{q}_{1-\alpha },y_i\ge y_j+\sqrt{{\sigma }_i^2+{\sigma }_j^2}{q}_{1-\alpha }\}$$

Full text

**Keywords:**Tukey’s HSD, rankability, Monte-Carlo, hospitals ranking, multiple comparisons.

1 Introduction

Estimation of ranks is an important statistical problem which appears in many applications in healthcare, education and social services (Goldstein and Spiegelhalter, 1996) to compare the performance of medical centers, universities or more generally institutions. Estimates of ranks have generally a great uncertainty so that confidence intervals (CIs) become crucial (Marshall and Spiegelhalter, 1998; Goldstein and Spiegelhalter, 1996). It is surprising that inference of ranks has received little attention in the statistical literature. In applications, ranks are rarely accompanied with CIs and if so these are generally pointwise. This paper presents a method to produce simultaneous CIs at a prespecified joint level $1-\alpha$ for the ranks with correct coverage of the true ranks. Simultaneity is important in the context of ranking estimation whenever we are not interested in a specific named institution but rather in all the institutions together. Simultaneity is also necessary to quantify the uncertainty about which institutions are ranked best, second best, etc.

In the literature, authors focus usually on pointwise CIs for the ranks. We mention the bootstrap method of Goldstein and Spiegelhalter (1996) which was widely used, see Marshall and Spiegelhalter (1998), Gerzoff and Williamson (2001) and Feudtner et al. (2011) among others. Methods based on empirical Bayes approaches were also considered, see Laird and Louis (1989), Houwelingen et al. (1999), Lin et al. (2006), Lin et al. (2009), Lingsma et al. (2009), Noma et al. (2010) and Jewett et al. (2018) among others. We also mention funnel plots, see Tekkis et al. (2003), Spiegelhalter (2005) among others. These latter two approaches although have been considered in comparing institutions, they do not aim to build (simultaneous) CIs for ranks. On the other hand, it was pointed out by Hall and Miller (2009) and Xie et al. (2009) that the bootstrap pointwise CIs and Bayesian methods fail to cover the true ranks in the presence of ties or near ties among the compared institutions. Testing pairwise differences between means was also used to produce pointwise CIs for ranks (Lemmers et al., 2007, 2009; Holm, 2012; Bie, 2013). Lemmers et al. (2007) tested pairwise differences among Dutch hospitals by calculating Z-scores for their performance indicators, but they did not correct for multiple testing and thus their CIs for ranks are not simultaneous. Holm (2012) (see also Bie (2013)) calculated also a Z-score, but he applied Holm’s sequential algorithm to correct for multiple comparisons on the institution level, that is for each institution he corrects for comparisons with other institutions. Nevertheless, this is only sufficient if we are interested in one of the institutions, but it is not sufficient to produce simultaneous confidence intervals for the ranks of the institutions.

To the best of our knowledge, there is only one method in the literature introduced by Zhang et al. (2014) where a Monte-Carlo method is proposed in order to produce simultaneous CIs for ranks. The method was adopted later in some recent papers such as Waldrop et al. (2017), Moss et al. (2017), Huang et al. (2018) and Moss et al. (2018) among others. The method of Zhang et al. (2014) can be seen as a generalization of the method proposed by Goldstein and Spiegelhalter (1996) and can be seen as a parametric bootstrap method. It is however not clear why this method should actually have a simultaneous coverage of at least $1-\alpha$. Besides, since it depends on the method of Goldstein and Spiegelhalter (1996), we argue that it inherits the lack of correct coverage. We show through extensive simulations that the method of Zhang et al. (2014) has the desired simultaneous coverage only when the means are quite far from each other with no way of determining the range of means since this depends on the number of means, how they are scattered and also the standard deviations. We also show that it is anticonservative when the means are close to each other.

We present a novel method which uses Tukey’s honest significant difference (HSD) test (Tukey, 1953). We show that Tukey’s HSD can be used to produce simultaneous confidence intervals for ranks with simultaneous coverage of at least $1-\alpha$.

When the means have no ties, our method becomes conservative. We show that it is possible to adjust the confidence level so that we reduce the conservativeness of the method. We show similarly, it is possible to repair the method of Zhang et al. (2014) in order to regain control of the confidence level. After rescaling, both methods seem to produce similar results, but the method of Zhang et al. (2014) becomes extremely difficult to repair as the number of means exceeds 30 to 50.

We also introduce in this paper a new rankability measure defined as the proportion of pairs of institutions that have different true performances. We estimate the true rankability by our method, and provide a lower confidence bound for it.

The paper is organized as follows. In Section 2, we explain the context of this paper, the notations and the objective. In Section 3, we revisit Tukey’s HSD and show that it can be used to provide simultaneous confidence intervals for the ranks. In Section 4, we review the Monte-Carlo method of Zhang et al. (2014). In Section 5, we show how to rescale the confidence level of our method and Zhang et al. (2014)’s method. Our new rankability measure is presented in Section 6. Section 7 is devoted to simulation studies comparing our method with Zhang et al. (2014)’s method with and without rescaling the coverage. An example of ranking Dutch hospitals is also discussed. Proofs of the propositions are in Appendix. Software for the methods presented in this paper is available in package ICRanks downloadable from CRAN.

2 Context and Objective

Let $\mu_{1},\cdots,\mu_{n}$ be $n$ real valued numbers which represent for example the true performance of the institutions we want to rank. Let $y=(y_{1},\cdots,y_{n})$ be a sample of $n$ independent random variables drawn from Gaussian distributions in the following manner

[TABLE]

where the standard deviations $\sigma_{1},\cdots,\sigma_{n}$ are known whereas the centers $\mu_{1},\cdots,\mu_{n}$ are unknown. The sample represents the observed performance indicators. Denote $r_{1},\cdots,r_{n}$ the true ranks of the centers respectively which are the target of inference. Our objective is to build simultaneous CIs for these ranks. Let us first define the ranks $r_{1},\cdots,r_{n}$, allowing for the possibility of ties.

Definition 2.1** (ranks).**

We define the lower-rank of center $\mu_{i}$ by

[TABLE]

We also define the upper-rank of center $\mu_{i}$ by

[TABLE]

We finally define the set-rank of $\mu_{i}$ as the set of natural numbers $r_{i}=\{l_{i},l_{i}+1,\cdots,u_{i}\}$ denoted here $[l_{i},u_{i}]$.

Assumption 1.

If the means have no ties, then $r_{i}=l_{i}=u_{i}$. Thus, the set-ranks coincide with usual ranking definition; the ranks are calculated for each mean by counting down how many means are below it.

When there are ties between the centers, we suppose that each of the tied centers possesses a set of ranks $r_{i}=[l_{i},u_{i}]$. For example, assume that we only have 3 centers $\mu_{1},\mu_{2}$ and $\mu_{3}$ such that $\mu_{1}=\mu_{2}<\mu_{3}$. Then, the rank of $\mu_{1}$ is the set $\{1,2\}$ and the rank of $\mu_{2}$ is also the set $\{1,2\}$, whereas the rank of $\mu_{3}$ is the singleton $\{3\}$. The rationale of the definition of the set-ranks is that in case of ties, the ranking is arbitrary, and a small perturbation of the true performance may produce any rank in the set of ranks. We call the ranks induced from the observed sample $y$ the empirical ranks. These ranks might be different from the true ranks of the centers, and since the sample is assumed to have a continuous distribution, the empirical ranks are all singletons.

We aim on the basis of the sample $y$ to construct simultaneous confidence intervals for the set-ranks of the centers. In other words, for each $i$ we search for a confidence interval $[L_{i},U_{i}]$ such that:

[TABLE]

for a prespecified confidence level $1-\alpha$. It is worth noting that the confidence intervals here are confidence intervals in $\mathbb{N}$, the set of natural numbers.

Two different types of statement can be obtained from the simultaneous CIs (2.4). First, for each center what are the possible ranks that it might take (which is our main objective). Second, since the confidence intervals for the ranks are simultaneous, we can deduce confidence sets for the best center(s), second best center(s), etc. These confidence sets have also a joint confidence level of at least $1-\alpha$. Indeed, in order to find the centers that can be the best, it suffices to see who are the centers whose rank CI starts at 1. In the same way, we can look at the centers whose rank CI includes rank 2 to obtain a confidence set of the centers ranked second best and so on.

3 Simultaneous Confidence Intervals for Ranks Using Tukey’s HSD

Tukey’s pairwise comparison procedure (Tukey, 1953) best known as the honest significant difference test (HSD) is an easy way to compare means of observations with (assumed) Gaussian distributions especially in ANOVA models. The interesting point about the procedure is that it provides simultaneous confidence statements about the differences between the means and controls the FWER at level $\alpha$. Moreover, it possesses certain optimality properties. In balanced one-way designs (which corresponds in our context to the situation that all $\sigma_{i}$’s are equal), simultaneous confidence intervals for the differences have confidence level exactly $1-\alpha$. The method is also optimal in the sense that it produces the shortest confidence intervals for all pairwise differences among all procedures that give equal-width confidence intervals at joint level at least $1-\alpha$, see for example Hochberg and Tamhane (1987, p. 81) and Rafter et al. (2002).

We consider the general case with possibly unequal $\sigma_{i}$’s here. Tukey’s HSD tests all null hypotheses $H_{i,j}:\mu_{i}-\mu_{j}=0$ at level $\alpha$ using the rejection region

[TABLE]

where $q_{1-\alpha}$ is the quantile of order $1-\alpha$ of the distribution of the Studentized range

[TABLE]

and $\tilde{Y}_{1},\cdots,\tilde{Y}_{n}$ are independent centered Gaussian random variables with standard deviations $\sigma_{1},\cdots,\sigma_{n}$ respectively.

In practice, a simple way to construct the confidence intervals for the ranks is to start by sorting the observations $y_{1}<y_{2}<\cdots<y_{n}$. In order to calculate the CI for $\mu_{i}$, it suffices to count down how many centers are not significantly different from it. The lower bound of the rank of $\mu_{i}$ is thus obtained by counting the number of times the hypothesis $\mu_{i}=\mu_{j}$ for $j<i$ is not rejected, say $s_{i}$, or equivalently the number of times the test statistic is below the Studentized range quantile. We then count down how many times the hypothesis $\mu_{i}=\mu_{k}$ for $k>i$ is not rejected, say $t_{i}$. The confidence interval for the rank of $\mu_{i}$ is then $[i-s_{i},i+t_{i}]$.

Proposition 3.1**.**

Let

[TABLE]

The intervals $[L_{i},U_{i}]$ for $i=1,\cdots,n$ are $(1-\alpha)$-joint confidence intervals for the ranks of means $\mu_{1},\cdots,\mu_{n}$.

Suppose we have three institutions $A,B$ and $C$ with centers $\mu_{A},\mu_{B},\mu_{C}$ respectively. Assume that we found the following $90\%$ confidence intervals for the differences from Tukey’s HSD (rounded to 1 digit)

[TABLE]

Then center $A$ gets a confidence interval for its rank $[1,1]$, center $B$ gets a confidence interval for its rank $[2,3]$ and center $C$ gets a confidence interval for its rank $[2,3]$.

We do not have a general idea about an optimal procedure to build simultaneous confidence intervals for the ranks. Still, Tukey’s HSD is known to have some optimality in balanced designs (all standard deviations are the same). Indeed, it produces the shortest confidence intervals for the differences among all (one-step) procedures which give equal-width confidence intervals at joint level $1-\alpha$. Since this optimality is related directly to the differences among the means, then it leads to the same optimality concerning ranks. In other words, any other method providing confidence intervals for the ranks based on confidence intervals for the differences with equal-width will produce wider CIs for the ranks than the ones produced by our Tukey-based method.

Proposition 3.2**.**

Under the full null, that is when $\mu_{1}=\cdots=\mu_{n}$, the simultaneous coverage of the confidence intervals $[L_{i},U_{i}]$ for $i=1\cdots,n$ produced by Tukey’s HSD is exactly $1-\alpha$.

This result together with Proposition 3.1 state that our Tukey-based method provides simultaneous CIs for the ranks with exact simultaneous confidence level $1-\alpha$ when all the means are the same. Otherwise, the method is conservative. The following sections will shed light on this conservativeness when Assumption 1 holds, that is when there are no ties among the means, and we are going to propose a practical solution that reduces this conservativeness.

Several step-down improvements on Tukey’s HSD have been proposed, see (Rafter et al., 2002) for a review. The most efficient and well-known is the REGWQ. Instead of testing equality of pairs of means, the procedure tests blocks of equality of means. Step-down variants of Tukey’s HSD control the FWER at level $\alpha$, but do not provide any directional information about the relative position of the centers (no protection against type III errors), so that no information about the ranks can be derived. THe step-down Tukey has not been proven to protect against type III errors, that is the ranking error, ((Welsch, 1977)) although some authors believe that it does. Therefore, we decided not to consider this approach to build CIs for ranks.

4 Simultaneous Confidence Intervals for Ranks Using Zhang et al.’s Method

The method of Zhang et al. (2014) is the first method (as far as we know) in the literature that was proposed to build simultaneous confidence intervals for ranks and can be seen as a generalization of the method of Marshall and Spiegelhalter (1998) in some sense. The method proceeds as follows. They use the method introduced by Marshall and Spiegelhalter (1998) to produce pointwise CIs at level $1-\beta$ for the ranks for several values of $\beta$ in the interval $(0,\alpha)$. For this purpose, $K$ $n-$samples are generated. These are then used again to estimate by Monte-Carlo the joint probability that the empirical ranks (the ranks of $y$) are inside these pointwise CIs at levels $1-\beta$. They choose $\beta$ such that the set of pointwise CIs has the smallest estimated joint probability superior to $1-\alpha$ that they include the empirical ranks. According to Zhang et al. (2014), as $K$ increases, we should obtain simultaneous CIs with a more accurate confidence level superior than $1-\alpha$. They also provide a lower bound for $K$ and advise the reader to choose a sufficiently larger value.

This Monte-Carlo-based method does not have a solid theoretical assurance and its true simultaneous coverage was never tested on simulated data. In the next sections, we are going to investigate this with more details. We are going to show that the simultaneous coverage of the resulting CIs reaches the nominal level $1-\alpha$ only when the differences among the means are large enough. This depends not only on the range of values of the means, but also on $n$ and the way the means are dispersed in that range. Otherwise, the method is anticonservative and the simultaneous coverage could reach very low levels beyond the nominal level and the resulting CIs would no longer be reliable. On the other hand, we will propose a method to readjust the simultaneous coverage that works in practice only when the number of means is below 50.

5 Simultaneous confidence intervals for ranks when ties are not allowed

It might be reasonable in the context of ranking to assume Assumption 1, hence the means $\mu_{1},\cdots,\mu_{n}$ are all different, and that there are in reality no ties. We first start by treating the case when the standard deviations are the same, that is $\sigma_{1}=\cdots=\sigma_{n}=\sigma$. We move then to the general case of different standard deviations.

5.1 Worst case configuration

When there are no ties, our Tukey-based approach becomes more conservative whereas the Monte-Carlo method of Zhang et al. (2014) becomes less anticonservative. Moreover, as the differences among the means become greater, the coverage probability of these methods should increase. We illustrate this in the left part of figure (1) by considering vectors of means of the form $\varepsilon\mu$ where $\varepsilon\in(-1,1)$ and $\mu=(1,\cdots,10)^{t}$ with dimension 10. The common standard deviation is set to 1. When $\varepsilon=0$, we assume arbitrary ranks for the means, say $1,\cdots,10$, in order to stay conform with the assumption that there are no ties.

The coverage probability for both methods reaches a minimum when $\varepsilon=0$. This means that the worst case happens when the means are arbitrarily small, say zero, while not having ties. This means, for any $\mu\neq 0$

[TABLE]

This worst case configuration is known in hypothesis testing for example when we test if the vector of means has an ascending order (Robertson and Wegman (1978)); the type I error is then highest when all the means are equal. We also mention the Kramer-Tukey procedure (Tukey (1953); Kramer (1956)); Hayter (1984) showed that it is conservative and has a worst case configuration when all the standard deviations are the same. In these procedures, the worst case configuration entails exactness of the method, that is the type I error is exactly $\alpha$. In our Tukey-based method, the worst case configuration entails a simultaneous coverage clearly higher than the nominal level. This gap can be further exploited in order to gain more power and reduce the conservativeness of our method. On the other hand, the Monte-Carlo method of Zhang et al. (2014) can be readjusted so that the simultaneous coverage no longer goes below the nominal level.

5.2 Rescaling the coverage at the worst case to the nominal level

We propose to regularize both methods (ours and Zhang et al. (2014)’s) so that on the one hand, the Monte-Carlo method of Zhang et al. (2014) delivers always a simultaneous coverage of at least $1-\alpha$ (a scaling up), and on the other hand, our Tukey-based method delivers a simultaneous coverage of at least $1-\alpha$ but in a less conservative way (a scaling down). Let $\mu$ be some vector of means, assume that at the simultaneous confidence level $1-\alpha$ we get an actual coverage of $\beta_{\mu}(\alpha)$ that is

[TABLE]

We look for $\tilde{\alpha}$ such that

[TABLE]

In other words, we look for a zero of the function $z\mapsto\beta_{\mu}(z)-1+\alpha$ in the interval $(\alpha,1)$ for our Tukey-based method, and in the interval $(0,\alpha)$ for the Monte-Carlo method of Zhang et al. (2014). This can be performed using any mathematical program, for example using function uniroot available in the statistical program R (R Core Team (2017)). In practice, this is not feasible because $\mu$ is unknown. Therefore, we rescale the coverage at the worst case and use the rescaled confidence level in order to calculate the CIs for any (unknown) $\mu$ different from the null vector. Since the worst-case configuration verifies (5.1), the resulting CIs have simultaneous coverage of at least $1-\alpha$.

The worst case configuration is characterized by having all means infinitely small. Therefore, in practice we set $\mu=(0,\cdots,0)\in\mathbb{R}^{n}$ and assume arbitrary ranks to its coordinates, for example $1,\cdots,n$. Then we look for $\tilde{\alpha}$ such that

[TABLE]

We have now, for any $\mu\neq 0$,

[TABLE]

If we do so, the simultaneous coverage of both our Tukey-based method and the Monte-Carlo method of Zhang et al. (2014) will be equal to the nominal level $1-\alpha$ near zero and higher than $1-\alpha$ elsewhere as illustrated in the right part of figure (1).

Table (1) shows the rescaled significance level $\tilde{\alpha}$ necessary to reach an actual coverage of $80\%,90\%$ and $95\%$ when the number of centers increases from $10$ to $100$. The table shows that in order to use the method of Zhang et al. (2014) and make sure not to be anticonservative, we need to use very small values of the significance level. However, as $\tilde{\alpha}$ becomes smaller we need to increase the number of Monte-Carlo samples $K$ required to estimate the joint distribution of the ranks as mentioned by Zhang et al. (2014). For example, when $n=50$ it is required that we generate at least $K=10^{6}$ $n-$samples, and thus rescaling the method of Zhang et al. (2014) becomes quickly infeasible for higher number of means so that the resulting CIs are not ensured to have the desired coverage of $1-\alpha$. On the other hand, our Tukey-based method although the rescaled significance level moves towards 1, the resulting CIs will always have simultaneous coverage of at least $1-\alpha$ even if we do not fully rescale the significance level at the worst configuration to $1-\alpha$.

5.3 The unequal sigma case

When the standard deviations are not the same, the worst case configuration still happens when the means are arbitrarily close to each other, but the order of the $\sigma$’s has an influence on it. Therefore, we need to find a worst-case ordering of the $\sigma$’s which ensures that if we protect against it by properly rescaling the worst case, the method will be conservative for any other ordering. Still, the number of possibilities is huge, that is $n!/2$ possible configurations of the standard deviations. By inspecting these configurations in the case when $n=10$, we found out that the worst-case ordering of the $\sigma$’s should happen when the lowest standard deviations are attributed to the lowest and the largest ranked centers whereas the highest standard deviations are attributed to the middle ranks (for example configuration (5.3)).

We show in three different scenarios how the worst case scenario happens nearly at zero, figure (2). We consider the vector of standard deviations $(1/10,1/9,\cdots,1)$ and a vector of means $\mu=(1,\cdots,10)$. We set $\varepsilon\in(-1,1)$ as before and estimate the coverage for the vectors $\varepsilon\mu$ when the standard deviations are ordered in one of the following manners (two tree orderings and one ascending order)

[TABLE]

The estimated coverage at the worst case configuration attains its smallest value under configuration (5.3).

6 A Rankability Measure

It is useful to have a single measure that gives an impression how well we can distinguish different means, that is how rankable they are. A set of equal centers is evidently not rankable. Therefore, this set of centers should get a rankability of 0. On the other hand, a set of totally different centers should get a rankability of 1 (or $100\%$) since we can rank each center. As the ranks are observed through quantities provided with uncertainty, an estimate of the ”true” rankability should be considered along with a confidence interval. We will first define the estimand before we define the estimate and its CI.

Assume we have $n$ centers $\mu_{1}\leq\cdots\leq\mu_{n}$. Some of these centers might be equal. According to our definition of ranks, equal (or tied) centers all get a set of ranks $[l_{i},u_{i}]$ which is the same for all of them. Define the rankability $R_{n}$ by

[TABLE]

The normalization by $n(n-1)$ is necessary for the rankability $R_{n}$ to take values in the interval $[0,1]$. The sum gives the surface of the set-ranks (the light grey area in figure (3)) and the subtraction from one ensures that if the set-ranks cover the whole range of ranks, we conclude that the centers are not rankable and we say then that the set of centers have a rankability of 0. In figure (3), the true rankability is $R_{20}=0.616$. The surface of the region in light grey (normalized by $n(n-1)$) in figure (3) can be interpreted as the probability that two centers $\mu_{i}$ and $\mu_{j}$ picked at random have the same rank. Therefore, our rankability measure $R_{n}$ can be interpreted as the probability that two centers picked at random get different ranks.

The rankability $R_{n}$, since it is defined through the true set-ranks, is a parameter that may be estimated. Denote $[L_{i},U_{i}]$ the confidence interval for the set-rank of $\mu_{i}$. We assume that these CIs have joint confidence level of $1-\alpha$, that is

[TABLE]

Define the estimated rankability at level $1-\alpha$ by

[TABLE]

Due to inequality (6.1), the estimated rankability at level $1-\alpha$ underestimates the true rankability with a probability at least $1-\alpha$. In other words

[TABLE]

Since $R_{n}\in[0,1]$, the interval $[\hat{R}_{n}(\alpha),1]$ becomes a $1-\alpha$ confidence interval for $R_{n}$.

In figure (3), we show the $50\%$ simultaneous CIs for ranks calculated using Tukey’s HSD on a sample of 20 centers resulting in a $50\%$ CI for the $R_{n}$ which is $[0.232,1]$. $\hat{R}_{n}(0.5)$ underestimates the true rankability $R_{n}$ with probability at least $50\%$, and it thus overestimates it with probability at most $50\%$ as well which makes from $\hat{R}_{n}(0.5)$ a good candidate for a conservative point estimate of $R_{n}$. We also show in figure (3) the $95\%$ simultaneous CIs produced by Tukey’s HSD, and the resulting $95\%$ CI for $R_{n}$ is then $[0.126,1]$.

It is worth mentioning that the estimated rankability can also be seen as a performance (or error) index so that several methods providing simultaneous CIs can be compared based on their estimated rankability.

In the context of empirical Bayesian methods for estimating ranks, a rankability measure was proposed by Houwelingen et al. (1999). It indicates which part of variation between hospitals is due to true difference and which part is due to chance. Rankability is then computed by relating heterogeneity between the centers to uncertainty between and within the centers, see also Lingsma et al. (2009) and Henneman et al. (2014) among others. This measure is specific to the Bayesian method and cannot be used in our case, however our rankability measure is only related to the confidence intervals regardless of the method which produce them. The only requirement is that the confidence intervals are simultaneous.

7 Simulation Study and Real Data Analysis

In this section we provide several examples (real and simulated) to demonstrate the confidence intervals produced using our approaches from Sections 3 and 5. We also compare the coverage and the efficiency of the confidence intervals produced by the method proposed by Zhang et al. (2014) to the ones obtained by our method in different simulated scenarios. The efficiency is calculate as $1-\hat{R}_{n}$ where $\hat{R}_{n}$ is given by (6.2). It represents the average lengths of the CIs.

Finally, we consider a dataset for patients with abdominal aneurysms from 64 hospitals in the Netherlands. We compare these hospitals according to the mortality rate at 30 days and then according to the type of surgery operated on the patient. All the simulations and the data analysis are done using the statistical program R Core Team (2017), and the code of the functions is available in the R package ICRanks which can be downloaded from the CRAN repository. The R code for the method of Zhang et al. (2014) is provided in the supplementary materials.

7.1 The case of a common standard deviation

The simulation setup is the following. We aim to estimate the average simultaneous coverage of the Monte-Carlo method of Zhang et al. (2014) and our Tukey-based method. To do so, we generate the centers $\mu_{i}$’s independently from the Gaussian distribution $\mathcal{N}(0,\tau^{2})$ for $\tau=0.5,1,2$. For each value of $\tau$, we generate $1000$ $n-$samples of centers $\mu=(\mu_{1},\cdots,\mu_{n})$ for $n=10,30$ and $50$. Then a Gaussian vector $y$ is generated from the multivariate Gaussian distribution $\mathcal{N}(\mu,I_{n})$. The simultaneous coverage based on these samples is estimated. The rescaled values of $\alpha$ for both methods are already calculated in table (1). We provide a table of the estimated coverage before and after rescaling the significance level so that the actual coverage at the worst case becomes $1-\alpha$ for $\alpha=0.1$. We calculate also the average $1-\hat{R}_{n}(\alpha)$ where $\hat{R}_{n}(\alpha)$ is the rankability measure (6.2).

We conclude from the tables the following points.

The method of Zhang et al. (2014) although provides clearly shorter confidence intervals for ranks, this comes at the cost of a low simultaneous coverage. Therefore, its results are generally unreliable and do not fulfill the requirement that it has a confidence level of at least $1-\alpha$. 2. 2.

The simultaneous coverage of the method of Zhang et al. (2014) increases as the range of means increases at a fixed $n$. On the other hand, it decreases as $n$ increases when the range of the means is held fixed. 3. 3.

The simultaneous coverage of our Tukey-based method increases with both the number of means $n$ and their range of values. 4. 4.

In average, our Tukey-based method seems to produce shorter CIs than the method of Zhang et al. (2014) when they are both rescaled, but this difference is not statistically significant. Indeed, the difference in efficiency was all the time within one simulation standard error. For example, when $\tau=0.5$ and $n=10$, we observed a standard error in $1-\hat{R}_{n}$ of about 0.01. We may state that when the two methods are properly rescaled, they do not have substantial differences. 5. 5.

Reducing the conservativeness of our Tukey-based method is always possible. 6. 6.

Repairing the anticonservativeness of the method of Zhang et al. (2014) is only possible in practice for $n\leq 50$.

7.2 The case of different standard deviations: A dataset on hospitals in the Netherlands

We study a dataset for Dutch hospitals concerning abdominal aneurysms surgery. The study includes 9489 patients operated at 64 hospitals in the Netherlands at dates mostly between the years 2012 and 2016. The number of patients per hospital ranged from 3 to 358 with an average of 150 patients per hospital. The dataset includes the following variables

• •

the hospital ID where the patient was treated;

• •

the date of surgery;

• •

the context of surgery: Elective, Urgent, Emergency;

• •

the surgical procedure: ”Endovascular”, ”Endovascular converted” and ”Open”. ”Endovascular” means the patient had a minimal invasive procedure through the femoral artery in the groin. ”Endovasculair converted” means the surgeons first tried a minimal invasive procedure through the femoral artery in the groin, but then realized they had to do an open surgery;

• •

a complication within 30 days (yes or no);

• •

the mortality within 30 days (yes or no);

• •

VpPOSSUM: a numerical score that summarizes the pre-operative state of the patient.

In order to conform to the normality assumption in our model, we exclude hospitals with small number of patients. This leaves us with 61 hospitals and each one of them has at least 54 patients. We compare these hospitals according to the mortality rate within 30 days. We correct for case-mix effect with a fixed effect logistic regression model using the VpPOSSUM variable. One of the hospitals has no patients who died within 30 days after surgery. Thus, we added to all the hospitals a row of data with a virtual patient who died within 30 days after surgery and with a value of VpPOSSUM equals to the average in the corresponding hospital. This prevents the logistic regression from getting an infinite standard error for this hospital. Besides, the influence on the other hospitals is rather minor because of the relatively high number of patients in them.

Before we apply the methods we presented in this paper, we calculate the rescaled significance level at the worst case configuration, that is we consider a null vector of means and the vector of standard deviations ordered according to the worst configuration we found in paragraph (5.3), namely configuration (5.3). We use the first 10 hospitals, 30 hospitals and finally all the hospitals and see how the rescaled significance level changes.

In order to apply the Monte-Carlo method of Zhang et al. (2014) and make sure that the confidence level is at least $90\%$, we need to use a significance level below $5\times 10^{-6}$ (table (5)). Using any value for $\tilde{\alpha}$ superior than $5\times 10^{-6}$ on the full data would result in a miss leading conclusion, because there is no guarantee that the simultaneous confidence level is actually at least $90\%$. Therefore, we avoid using it on the full dataset and only apply our Tukey-based methods. The simultaneous CIs for the ranks of the hospitals at the joint level $90\%$ are illustrated in figure (4).

The confidence intervals cover the whole range of ranks, and there are barely any differences among the hospitals according to the mortality rate. The rankability is $0.001$ for the Tukey-based method without rescaling, and is $0.012$ after rescaling. This can either be normal, that is all Dutch hospitals have the same performance, or due to a low power of our methods. In order to find out, we change the output variable in the logistic regression model and correct for case-mix effects with the type of surgery as an output. The resulting CIs at joint level $90\%$ for the ranks are illustrated in figure (6) with a rankability of 0.240 for Tukey’s HSD. Rescaling the significance level clearly improves the results of our Tukey-based CIs. The rescaled significance level is $\tilde{\alpha}=0.607$. The new rankability is $0.358$. Here again, we could not apply the method of Zhang et al. (2014) because we obtain similar results to table (5). We may state that for only 7 hospitals we find that they may get the first rank, and that for the remaining 54 hospitals we can confidently state they are not first rank.

We make a forest plot for the hospital effect after case-mix correction for both the mortality withing 30 days and the surgical procedure. Figure (5) shows that indeed the mortality rate induces very few differences among the hospitals whereas the type of surgery seems to show more differences. We also fit a random effect mixed-model to the hospital effect for both choices (separately) using function rma from package metafor (Viechtbauer (2010)) and estimate the variance of the random effects using the Sidiki-Jonkman method and test for heterogeneity among the hospitals. We find a p-value of $0.951$ while using the mortality rate and a p-value lower than $0.001$ while using the type of surgery. This supports our claim that the reason behind the very wide CIs for ranks when we use the mortality rate is not the low power of our method but is rather that there are not much differences in the data.

8 Discussion

We presented a novel method to produce simultaneous CIs for ranks based on Tukey’s HSD and proposed a practical improvement under Assumption 1 that there are no ties among the means. The only method in the literature (as far as we know) claiming to produce simultaneous confidence intervals for ranks is the simulation-based method of Zhang et al. (2014). We showed through simulations that the simultaneous confidence level goes below the nominal level unless the means are very far from each other. We proposed a solution to fix this problem by rescaling the confidence level. Surprisingly, after rescaling, the results of our Tukey-based method and the method of Zhang et al. (2014) seem to be almost the same and the differences are not significant.

By providing valid methods for simultaneous CIs for ranks, practitioners are provided with a way which looks at all the institutions together instead of only looking at a specific one. Simultaneity also provides a way to state which of the institutions could be ranked first best, second best, and so one, and which of these institutions may not get to be first rank. These pieces of information could not be obtained only be looking at pointwise CIs for the ranks.

We found it impractical to rescale the significance level for Zhang et al. (2014)’s method when the number of means is higher than 30 whereas our Tukey-based method is always rescalable. While comparing the performance of 64 Dutch, it was only possible to use our Tukey-based method. The comparison of the hospitals according to their mortality rate at 30 days after surgery showed no differences among the hospitals. However, by considering the type of surgery as our primary output, some differences among the hospitals became quite clear (especially the extremities) and we were able for example to detect 7 hospitals which have the chance to be first rank.

Using Dunnet’s test, it is possible to look for a confidence interval for the rank of only one prespecified institution that we are interested in. On the other hand, in our approach we considered two-sided CIs. We could also be interested in looking only at worse or better institutions and thus gain more power by only considering one-sided CIs.

SUPPLEMENTARY MATERIAL

9 Proofs of Propositions

9.1 Proof of proposition 3.1

From Tukey’s procedure, we can obtain simultaneous confidence intervals for the differences between the centers at level $1-\alpha$, that is $\mu_{i}-\mu_{j}$ for $i,j\in\{1,\cdots,n\}$, see Hochberg and Tamhane (1987, sec. 2.1). In other words, we have

[TABLE]

Denote $[a_{i,j},b_{i,j}]$ the confidence interval for the difference $\mu_{i}-\mu_{j}$ in the previous display. Define also $L_{i}=1+\#\{j:\;a_{i,j}>0\}$ and $U_{i}=n-\#\{j:\;b_{i,j}\leq 0\}$. Let $E_{i}=\{\mu_{i}-\mu_{j}\in[a_{i,j},b_{i,j}],\forall j\neq i\}$. It is easy to see that the event $E_{i}$ implies the event $\{l_{i}\geq L_{i},u_{i}\leq U_{i}\}$ for any $i$. Thus using inequality (9.1), we may write

[TABLE]

Hence, the confidence intervals for the set-ranks $[L_{i},U_{i}]$ have a joint level of at least $1-\alpha$.

9.2 Proof of proposition 3.2

Under the full null, the set-rank of any $\mu_{i}$ is $[1,n]$. We need to show that the event

[TABLE]

implies the event

[TABLE]

Since $L_{i}\in\{1,\cdots,n\}$ and $U_{i}\in\{1,\cdots,n\}$, then the event (9.2) is equivalent to

[TABLE]

The first part is equivalent to the event that there is no $j$ such that $\{y_{i}\geq y_{j}-\sqrt{\sigma_{i}^{2}+\sigma_{j}^{2}}q_{1-\alpha}\}$ occurs. Similar reasoning for the second part entails that the event (9.2) is equivalent to

[TABLE]

which is clearly the same as the event (9.3).

10 R Code for the Calculus of the Coverage of Zhang et al.’s Algorithm

10.1 Calculating the coverage

library(ICRanks) n = 10; TrueCenters = 1:n

Take a subset and generate the data

alpha = 0.05; sigma = rep(1,n)

K = 10^4 coverage = 100 coverageTuk = 100 for(i in 1:100) { y = as.numeric(sapply(1:n, function(ll) rnorm(1,TrueCenters[ll],sd=sigma[ll]))) ind = sort.int(y, index.return = T)$ix y = y[ind] resZhang = BootstrCIs(y, sigma, alpha = 0.05, N = K, K = K, maxiter = 10) resTukey = ic.ranks(y, sigma, Method = "Tukey", alpha = 0.05)

if(sum(ind<resZhang$Lower | ind>resZhang$Upper)>0) Ψcoverage = coverage - 1 if(sum(ind<resTukey$Lower | ind>resTukey$Upper)>0) ΨcoverageTuk = coverageTuk - 1 }

10.2 An R function to calculate the CIs for the ranks according to Zhang et al.’s method

BootstrCIs = function(y, sigma, alpha=0.05, N = 10^4, K = N, precision = 1e-6, ΨΨΨΨΨΨΨmaxiter = 50) {

A function which calculates the individual CIs at level beta

Spiegelhalter = function(mus,ses,beta, N = 10^4) { k=length(mus) r=apply(x,2,rank) r=apply(r,1,quantile,probs=c(beta/2,1-beta/2),type=3) df=data.frame(lower=r[1,],upper=r[2,]) return(df) } n = length(y) beta1 = 0; beta2 = alpha beta = (beta2 + beta1) / 2 x=sigmamatrix(rnorm(Kn),nrow=n) + y InitCIs = Spiegelhalter(y, sigma, alpha, N) counter = 0; coverage = K while(abs(beta1 - beta2)>precision | counter<=maxiter) {

Generate individual CIs at level beta

res = Spiegelhalter(y, sigma, beta, N)

Check the coverage

coverage = K for(j in 1:K) { Ψind = rank(x[,j]) Ψif(sum(ind<res$lower | ind>res$upper)>0) coverage = coverage - 1 } if(coverage/K >= 1-alpha) { Ψbeta1 = beta }else { Ψbeta2 = beta } beta = (beta2 + beta1) / 2

counter = counter + 1 } if(coverage/K < 1-alpha) beta = beta1 res = Spiegelhalter(y, sigma, beta, N) return(list(Lower = res$lower, Upper = res$upper, coverage = coverage/K)) }