Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
45,924 characters · 28 sections · 58 citation commands
Optimal Recovery for Causal Inference
A major component of policy evaluation is to estimate the effects of an implemented policy, so as to know whether it achieved its goals. Estimating effects further yield inferences of causal relationships between interventions and their outcomes. Quantifying the effect of a treatment has been a problem of interest not only in policy making but also across different domains in health sciences, social sciences, and engineering. Typically, the effect is measured by looking at the difference between outcomes before and after an intervention of a treated unit. However, for a given object (e.g., region) at a given time, only one of the outcomes is observed and not both. Thus, we aim to recover and estimate the outcome that was not observed.
Several causal inference methods have been developed for observational studies to estimate the unobserved outcomes for a given intervention. For example synthetic control (SC) abadie2021synthetic constructs a weighted average of control units to act as a synthetic control unit to compare with the treated unit. Recently, there has been a growing literature that approaches causal inference from a matrix completion perspective. Proposals include approximating the control unit matrix using nuclear-norm minimization athey2021matrix, using singular value decomposition Amjad18, or by finding nearest neighbors agarwal2021causal_matrix_completion for missing entries of a matrix to best match control units and the treated unit of interest.
A main limitation in previous work is that when there is insufficient data, especially data from only a small period of time, methods are unable to recover the true estimate Amjad2019_mRSC. Further, the insufficiency or the low quality of data tend to create a poor pre-treatment fit which is a main source of bias in estimates in SC abadie2021synthetic.
To address this, we propose an approximation-theoretic approach from signal processing called ellipsoidal optimal recovery (EOpR) based on minimizing the worst-case error, which guarantees an exact fit of the pre-intervention period and optimally recovers the unobserved outcome with reduced bias in the effect estimate. Further, under the assumptions underlying EOpR, we derive worst-case estimates of effects, which are very useful in policy evaluation. Especially in situations of uncertainty, it is important to acknowledge the most severe possible outcomes that could occur for a given policy. Determining a policy alternative is derived by possible good and bad outcomes it yields. EOpR provides policy evaluators with “worst-case” outcomes and “typical” outcomes to help in decision making.
To be specific, within the potential outcomes framework of Rubin rubin1974causaleffect, we consider a variation of the optimal recovery muresaP2004 algorithm to recover missing outcomes for causal inference and to obtain worst-case estimates of the causal effect. Unlike statistical approachesabadie2003economic, abadie2010synthetic or low-rank matrix decomposition approaches Amjad18, Amjad2019_mRSC, athey2021matrix, optimal recovery is expected to be robust kassam1985robust given that it minimizes the maximum error over the known samples muresan2005demosaicing.
In the remainder of this paper, first we review some aspects of causal inference. Then we describe our approach in the context of panel data. The efficacy of our approach is lastly supported by artificial and empirical experiments.
To fix concepts and common terms let us give a brief overview of causal inference in the context of comparative case studies with panel data. Then we introduce the problem formulation and fix notations. Later, in Section (ref), we describe our method.
Causal analysis takes a step further from standard statistical analysis, inferring beliefs under changing conditions to uncover causal relationships among variables pearl2009causalInferenc. Causal inference has been widely used in social and health sciences. Several frameworks to tackle causality analysis, such as structural models pearl2010causal and the potential outcome framework rubin1974causaleffect (which we focus on here) have been proposed.
This framework assumes effects are tied to a treatment or an intervention. To reveal the causal effects of an intervention, rubin1974causaleffect proposed to measure the difference of two potential outcomes; let us denote them as $Y_N$ and $Y_I$, for a given unit $x$. The potential outcome $Y_N$ is the outcome for $x$ without being exposed to an intervention, and $Y_I$ is the outcome after an intervention is applied on $x$. So, the causal effect is
However, in real applications, we can never observe both outcomes for the same unit under the same conditions, only one of the two will take place at a given time. Therefore, one of the potential outcomes will always be missing, and that is the core objective of the framework to estimate one of the outcomes.
Let us introduce the main terms used in the potential outcomes literature, which are used throughout. A unit is the atomic object in the framework, which can be any physical object, whether a patient, a city, or a collection of objects at a particular time. A treatment is the action applied on a unit to change its state. The treatment\footnote{The terms treatment and intervention are used interchangeably.} can be a medicine given to a particular group or a lockdown order during a pandemic. The treatment is usually thought of as binary, so one group receives the treatment (the treated group), and the other does not (the control group). Also, given an intervention at time $t_0$, there are two time periods: the \textit{pre-intervention/pre-treatment} period for $t < t_0$ and the \textit{post-intervention} period for $t > t_0$.
Several works in econometrics and related literatures estimate the unknown potential outcome to estimate the causal effects. Here, we focus on the commonly used synthetic control method.
Synthetic control (SC) abadie2010synthetic proposes a particular way to measure the missing observable potential outcome to estimate causal effects. In their evaluation of econometrics for policy evaluation, atheyI2016state asserted SC is “arguably the most important innovation in the policy evaluation literature in the last 15 years”. Instead of using a single control unit or a simple average of a set of control units, SC creates a synthetic unit to act as a control group by selecting appropriate weights for selected control units. The choice of weights should result in a synthetic control unit that best resembles the pre-intervention values of the treated unit. SC is based on a strong assumption that such weights exist if and only if the treated units' pre-treatment time series is inside the convex hull of the control units' time series.
Therefore, SC is subject to the curse of dimensionality, in which the probability that exact matching of weights vanishes as the number of time periods grows zhao2017entropy. In ferman2021synthetic, the demeaned SC (DSC) was proposed to relax SC constraints on weights to allow for a good pre-treament fit when the length of pre-intervention is very large. Also, SC tends to heavily depend on the selection of control units, and its estimates are biased by noisy control units. Robust SC Amjad18 (RSC) views the setting of SC as an instance of the Latent Variable Model. The observable outcomes of the control group are obtained by a low-rank approximation using singular value decomposition (SVD) and noise is eliminated. This also generalizes to cases where outcomes are missing.
Consider panel data\footnote{Panel data is another word for cross-sectional time-series data} having a collection of time series with respect to an aggregated metric of interest (e.g., country GDP). The data includes $N$ units observed over $T$ periods of time. Let $T_0$ be the intervention time, which splits the time period into a pre-intervention period with $ 1 \leq T_0 < T$ and a post-intervention period with length $T-T_0$. We fix $i=1$ ($i\in \{1,\dots, N\}$) for the treated unit at time $t$, hence, let $s_{1t} \in \mathbb{R}^{T}$ be the treated unit. The remaining units $i = 2, \ldots, N$ are the controls that are not affected by the intervention. Let $\bm{{S}} \in \mathbb{R}^{N-1\times T}$ be the control units matrix.
Let the outcomes of the control and treated units follow a factor model, a common model in the econometrics literature abadie2021synthetic. Let $x_{it}$ denote an aggregated metric for a unit $i$ at time $t$. In the absence of covariates and unobserved outcomes, the factor model is the following:
Following the literature in data imputation athey2021matrix,Amjad18 and optimal recovery for missing values muresaP2004, we consider $s_{it} = \eta_i \psi_t $, where $\eta_i \in \mathbb{R}^{N-1}$ and $\psi_t \in \mathbb{R}^{T}$ are the latent features that capture the unit and time specifications respectively for observed outcomes, with independent random zero-mean noise, $\epsilon_{it}$. The goal is to approximate the partially-observed treated unit vector, $s_{1}$, to recover missing outcomes.
To distinguish pre- and post-intervention periods, let $\bm{{S}} = [\bm{{S}}^-, \bm{{S}}^+]$, where $\bm{{S}}^-=\{s_{it}\}_{2\leq i \leq N, t\leq T_0}$, and $\bm{{S}}^+ = \{s_{it}\}_{2\leq i \leq N, T_0 < t \leq T}$. Vectors are defined in the same manner, i.e. $s_{i} = [s_i^+, s_i^-]$. The inverse of $\bm{{A}}$ is $\bm{{A}}^{-1}$. The Moore-Penrose pseudo-inverse of $\bm{{A}}$ is $\bm{{A^{\dag}}}$. The transpose of $\bm{{A}}$ is $\bm{{A}}^{\top}$. The $A$-norm of a vector $u$, denoted $||u||_A$, is the value of $u^{\top}\bm{{A}}u$.
In this section we describe our proposed method. In Section (ref), we briefly introduce the optimal recovery method from signal processing, which is a fundamental building block of our approach. In Section (ref) we describe our approach of estimation in the context of time-series panel data in comparative studies.
Optimal recovery was introduced to effectively approximate a function known to belong to a certain signal class with limited information about it optimalestimation_book_1977. It has been applied to estimate missing or corrupted pixels in images muresaP2004 and missing values in biological data chenV2019. Optimal recovery estimates the missing value using a learned signal class and a set of known values, and provides deterministic error bounds allowing the calculation of worst-case values at these bounds muresanP2001. One way the signal class is constructed is using an ellipsoidal set of vectors that pass through a hyperplane.
The value $h$ is the radius of the ellipsoid. The matrix $Q$ determines how far the ellipsoid extends in every direction from the center. The lengths of the ellipsoid semi-axes are determined by the eigenvalues of $Q$.
Here, we extend optimal recovery with ellipsoidal signal class to estimate causal effects in panel data. Optimal recovery not only imputes the missing parts of the vector of interest but also provides deterministic worst-case characterization for the estimate. Under the potential outcomes framework, causal inference can be treated as a missing data problem, a kind of matrix completion athey2021matrix. Therefore, we use optimal recovery from approximation theory to recover the missing outcomes for causal inference, in which optimal recovery has not been considered before.
Geometrically, we assume control units vectors in $\bm{{S}}$ belong to an ellipsoidal class $K$, and the pre-intervention of control units $\bm{{S}}^-$ belong to a hyperplane $\mathcal{H}$. We consider the treated unit $s_{1}$ as a partially-observed vector, where $s_{1}^+ = \{s_{1t}\}_{t>T_0}$ (in the post-intervention) is unknown and requires approximation. The vector $s_{1}$ lies in $C$, the intersection of $K$ and $\mathcal{H}$. The aim is to find an estimator that minimizes the worst-case error. This is equivalent to finding the Chebyshev center of $C$ boyd2004convex. Figure (ref) illustrates the geometrical setting of the optimal recovery approach.
The Chebyshev center provides a minimax optimal solution for the recovery problem boyd2004convex, shown in the following theorem.
To find the minimax estimator $\hat{s}_1$, first construct the ellipsoidal class $K$ as in (ref), and then extrapolate using the representors of the known samples.
We construct a covariance matrix $\bm{{\Sigma}}$ from the data matrix $\bm{{S}}$, such that $\bm{{\Sigma}} = \bm{{S}}\bm{{S}}^{\top} + \lambda \bm{{I}}$, where $\bm{{\Sigma}} \in \mathbb{R}^{T\times T}$, $\lambda$ is a scalar, and $\bm{{I}}$ is the identity matrix.
To learn the ellipsoidal class $K$ in (ref), we want $K$ to have the most stretch in the same direction of $\bm{{\Sigma}}$, hence we let the eigenvalues of $\bm{{Q}}$ be the reciprocal of the eigenvalues of $\bm{{\Sigma}}$,
\paragraph{Choice of parameter} Based on the ellipsoid definition (ref), the matrix $\bm{{Q}}$ must be positive definite. To ensure that eigenvalues of $\bm{{Q}}$ are strictly positive, we add a small perturbation $\lambda \in (0, 1]$ to the diagonal of $\bm{{SS^{\top}}}$, such that it minimizes the $\ell_2$-norm between outcomes in the pre-intervention period.
\paragraph{Learning the representors}
From the pre-intervention vectors $\bm{{\Sigma}}^{-}$, from the covariance matrix $\bm{{\Sigma}}$, we derive the representors $\bm{{\Phi}} \in \mathbb{R}^{T_0 \times T_0}$, by the Riesz representation theorem muresaP2004 as
Now we calculate $\hat{s}_{1}$, the Chebyshev center of $C$. Given that the representor vectors lie in a subspace that is parallel to the center of the ellipse $C$, $\hat{s}_1$ is a linear combination of the inverse of representors $\bm{{\Phi}}$, such that the optimal weights $w^{\ast}$ are
then, the Chebyshev center, the estimated outcome $\hat{s}_1$ is
Given that the estimator $\hat{s}_1$ is at the center of the intersection $C$, the vectors on the boundary of $C$ are the worst-case estimates. To attain worst-case vectors (minimax control of the estimates), let $y$ be the unit norm in $\mathcal{Z}$, a parallel subspace to representors $\bm{{\Phi}}$, determined by
The worst case estimates $\bar{s}_1$ are:
with a very small $\varepsilon$, and the $\bm{{Q}}$-norm $\|\hat{s}_1\|_Q = \hat{s}_1^{\top} Q \hat{s}_1$.
Based on the optimal recovery approach, the extrapolation step produces an estimator that is the Chebyshev center of the ellipse $C$ with desirable properties of unbiasedness and consistency.
The Chebyshev center has been also proved to be a consistent and unbiased estimator geometrically halteman1986chebyshev, full proof in Appendix (ref).
We compare the accuracy of our ellipsoidal optimal recovery (EOpR) approach against other causal inference methods used in policy evaluation: SC abadie2003economic, RSC Amjad18, DSC ferman2021synthetic, and SDID Arkhangelsky21. We first evaluate on simulated data to demonstrate properties of EOpR under certain settings. We then evaluate on two classical panel datasets commonly used in the SC literature (California Proposition 99 abadie2010synthetic and Basque Country abadie2003economic). We finally apply EOpR in the context of the COVID-19 pandemic to estimate the number of confirmed cases in New York State.
To measure the quality of estimation, we use two metrics. First, we measure the root-mean-square error (RMSE) of estimated signals. The pre-intervention (training) error is for $1 \leq t \leq T_0$, and a post-intervention (testing) error is for $T_0 < t \leq T_0$
where $\mathcal{T}$ is the size of the selected time period.
Second, Abadie abadie2010synthetic proposed a test statistic to evaluate the reliability of the estimates by running placebo tests. One placebo test considers one control unit as a placebo treated unit and apply the estimation algorithm. Since control units are assumed to not be affected by the examined intervention, one would expect that the estimated signal for the placebo unit does not diverge from its corresponding control unit. Further, the gaps between each placebo estimation and its corresponding control unit should be less divergent than the gap between the original treated unit and its estimation. Placebo tests are applied on the classical econometrics case studies.
We conduct artificial simulations to demonstrate the properties of EOpR estimates in both the pre- and post-intervention periods. We show that EOpR performs well and better than existing causal inference methods under various settings. \paragraph{Experimental setup} Consider a data generating process, which is frequently considered for low-rank matrix decomposition solutions, similar to Amjad2019_mRSC, as follows. First we create two sets of row and column features, $B_r$, $B_c$, where $B_r = \{b_k|b_k \sim \text{Unif}(0,1), 1 \leq k \leq 10\}$ and $B_c = \{b_k|b_k \sim \text{Unif}(0,1), 1 \leq k \leq 10\}$. For each unit $2 \leq i \leq N$, we assign a parameter $\theta_i$ drawn from $B_r$ (with replacement), and for each time $1 \leq t \leq T$ we assign a parameter $\rho_t$ drawn from $B_c$ (with replacement). We use the following formula to generate a data point $\tilde{s}_{it} = f(\theta_i, \rho_t)$ to construct the control units
where $\epsilon_{it} \sim \mathcal{N}(0,1)$, an independent Gaussian noise. To construct the treated unit $\tilde{s}_{1t}$ of interest, we generate the vector using a uniform linear combination of row vectors.
In the following experiments, we investigate EOpR resistance to bias in comparison to other algorithms under different sizes of units $N$ and time periods $T_0$, $T$. For each combination of $N$, $T_0$, and $T$ we generate $10$ simulations and average the resulting RMSE scores for the estimated pre- and post-intervention signals of $\tilde{s}_{1t}$.
Consider when the number of pre-treatment units vary proportionally to $T$, fix the size of units $N$ and the total size of time $T$. A short period of $T_0$ has shown to fail in reproducing the trajectory of the treated unit abadie2021synthetic. Therefore, we vary the time of intervention $T_0$ between 10% to 90% of the entire time $T$.
Figure (ref) shows the effect of different pre-treatment lengths on the algorithm's ability to estimate, with $N=50$, and $T=200$. When having either a small number or a large number of pre-treatment periods (relative to $T$), EOpR recovers the original treated signal with the minimum error compared to other algorithms. Specifically at 10% and 20% of $T$, where the pre-treatment period is short, EOpR extrapolates beyond the training periods with the lowest bias in the estimation of the post-intervention trend, whereas other algorithms have higher bias in the post-intervention estimation.
Additionally, at lengths of 40% and onwards, the number of pre-treament vectors are greater than the size of $N$, $(N \ll T_0)$, a common setting adapted in abadie2010synthetic, EOpR still maintains a lower bias in the post-intervention with consistent low training error.
With a large number of units, the risk of overfitting increases, which produces a high potential of increased estimation bias abadie2010synthetic. Here, we model the increase of the number of units with fixed $T_0 = 25$ and total $T = 125$. Large number of control units, $N$, has shown to be challenging for SC as it exacerbates the bias of the estimator abadie2021synthetic.
Figure (ref) shows the effect of a growing number of units $N$ on algorithm performance. Given that EOpR also achieves low error at the post-intervention estimation at greater sizes of $N$, EOpR has the ability to reconstruct the true signal trend even with large control units and potentially noisy settings.
We further investigate the ability of EOpR to estimate the trajectory of the post-intervention for an extended period of time. The ability to estimate for an extended period indicates an algorithm is robust against estimation bias.
Figure (ref) shows the estimation errors with fixed $N=100$ and $T_0 = 50$, and tested over multiple lengths of post-intervention. EOpR consistently achieves a low estimation error in both pre- and post-intervention estimation, especially at longer periods of time, e.g.\ $(t=450)$.
The objective of this case study is to investigate the effect of terrorist attacks on the Basque Country economy compared to other Spanish regions. Terrorist activities started by 1970. There was a significant negative impact on the economy of Basque Country measured by per-capita GDP. Abadie abadie2010synthetic showed that the economy would be better without terrorism. \paragraph{Results} Figure (ref) shows the actual trajectory of the Basque Country economy in black, with a degradation after 1970. In comparison to other methods, our EOpR method recovers Basque Country estimate, with more accurate fit on the pre-intervention values of the treated unit. Figure (ref) shows the estimated worst-case potential outcomes from EOpR.
\paragraph{Placebo Tests} We create placebo tests, similar to Abadie abadie2010synthetic. Note that in abadie2010synthetic, authors excluded 5 regions which had poor fit in the pre-intervention, but we keep all regions. We plot the differences between our estimates and the observations of all regions as placebo and Basque Country (the actual treated unit). Figure (ref) shows the differences for all regions compared to Basque Country (solid black line in the figure). The divergence for Basque Country was the largest, thus, the derived estimates by the EOpR are reliable.
The objective of this case study is to investigate the anti-tobacco legislation, Proposition 99, on the per-capita cigarette consumption in California in comparison to other states in the United States. The legislation took place in 1970. Without such legislation, the consumption of cigarettes in California would not have decreased abadie2010synthetic.
\paragraph{Results} Figure (ref) shows the actual trajectory of California cigarette consumption in black. Our method recovers the estimated signal better than other methods with an adequate fit on the pre-intervention outcome, and also it derives the worst-case estimates. \paragraph{Placebo Tests} We apply the same placebo test to the California case study. Abadie abadie2010synthetic excluded 12 regions, but we keep all of them. Figure (ref) shows the divergence between estimations and observations of all regions with solid black line for California. This shows a similar observation as in abadie2010synthetic.
=
New York was one of the earliest American states that turned to an epicenter of COVID-19 during 2020 thompson2020covid. We aim to estimate the New York COVID-19 cases trajectory with states as control units using EOpR and other methods. To make the states uniform, we select the states where lockdowns were imposed, and end up with a total of 43 states. \paragraph{Experimental setup} For lockdown dates, we use data from the COVIDVis project \footnote{\url{https://covidvis.berkeley.edu/##lockdown_section}} that tracks policy interventions at the state level. We consider the dates of shelter-in-place mandate. For COVID-19 cases, we consider state-level case load data from the New York Times database nytiems_cov. Note that since reported cases depend on testing, our analysis is limited by the fact there was widespread shortage of available tests in different regions at different times.
Since each state imposed a lockdown at different times, we aligned states based on the days differences between the time a lockdown took place and the rest of the dates during the period of interest, following bayat2020synthetic. Figure (ref) shows the trend of New York and 7 other states and their moving averages (over 7 days).
Based on the literature of SC abadie2010synthetic, one would need to select the control units based on their similar trends to the target unit to ensure a correct construction of the treated unit. In the following experiments, we estimate NY COVID-19 cases using $\{5, 7, 15, 43\}$ states, for a period of 200 days, $T=200$, which approximates the period between March and August 2020, and time of intervention $T_0 = 50$. \paragraph{Results} Table (ref) shows the errors produced by the four algorithms for pre- and post-intervention estimation. We start by $N=5$, and pick the states that have a very similar trend to the daily average of New York State, similarly for $N=7$ shown at Figure (ref). For a large $N$, EOpR recovers the trajectory of New York with the lowest estimation error regardless of the existence of noisy control units which are highly dissimilar to the New York COVID-19 trend. An example of some dissimilar states used is shown in Figure (ref). The ability of EOpR to recover the true trajectory under different sizes and similarity of control units shows its potential to learn a set of high-quality signals that are sufficient for recovery.
Figure (ref) shows the estimated trend of daily average cases of New York, the estimation starts at $T_0$, in comparison to the actual New York trend (treated unit in the figure). The worst-case estimates are also plotted. Note that the worst-case estimates range changes depending on the quality of control units.
We empirically study the effect of the added perturbation $\lambda$ to the covariate matrix $\bm{{SS^{\top}}}$. This small perturbation is added to ensure $Q$ is positive definite in (ref), following the definition of ellipsoids (ref). When $\lambda = 0$, we see the resulting predictions suffer drastically, deviating from expected possible outcomes (Figure (ref)). Using an appropriate $\lambda$ as in Figure (ref), balances the model complexity which helps safeguard the algorithm from potentially underfitting the training data, and producing a biased estimate.
A similar observation is seen for simulated data, for the setting of fixed $N$ and $T$ but changing the time of intervention $T_0$, in Figure (ref).
Classical synthetic control has been noted as effective for causal inference in comparative studies. Here, we propose a signal processing approach for synthetic control---ellipsoidal optimal recovery (EOpR)---that estimates the causal effect given a policy intervention. Further, given the properties of EOpR, we derive worst-case estimates, which are themselves very useful for policy evaluation. Our approach of EOpR has less estimation error for pre- and post-intervention periods, especially with short pre-intervention periods. This is demonstrated through comparisons on both simulated data and classical case studies in econometrics. Applications in health-relevant settings are also shown to be compelling. Placebo tests and ablation studies demonstrate robustness. Extensions to the basic optimal recovery framework beyond ellipsoidal signal classes that we develop here are also possible.
Authors would like to thank Akhil Bhimaraju and Alayt Issak for their useful insights and thoughtful discussions.
\appendices
To show consistency, here we bound the $\ell_2$ error of the estimation. We will drop the dependency on $\lambda$. Recall that the noise term $\epsilon_1$ is a zero-mean independent random variable that satisfies $\mathbb{E}(\epsilon_{ij}) = 0$ for all $i$ and $j$ by assumption, with variance $Var(\epsilon_{ij}) = \sigma^2$.
Proof. Recall from (ref) that the treatment row $x_1 = s_1+\epsilon_1$. By the definition of the Chebyshev center and its properties of unique estimation, consider $s_1 = \Sigma^T w^{\ast}$ with exact recovery, and $\hat{s_1} = \Sigma^T\hat{w}$.
Taking expectations, we arrive at the inequality
We must address an inner product on the right side. First, we derive some useful facts. Recall the trace operator has the mapping property $\mathit{tr}(AB) = \mathit{tr}(BA)$, and the projection matrix $P$ to be $P_1 = AA^{\dag}$ and $P_2 = A^{\dag}A$. Hence,
which follows since the trace of a projection matrix equals the rank of the matrix, i.e.\ $tr(\Sigma^T \Sigma^{\dag^T}) = \text{rank}(\Sigma^T)$. Hence, the rank of $\Sigma$ is at most $|r|$.
Returning to the inner product, recall $\hat{w} = \Sigma^{{\dag}^T} s_1$ from (ref).
Finally, we replace the above terms into inequality (ref) to arrive at
which completes the proof.
Now, let us consider replacing the above value in the mean squared error ($\mathit{MSE}$) function as
The following theorem halteman1986chebyshev demonstrates that the Chebyshev center is an unbiased and consistent estimator for the center of a normal distribution over a $k$-sphere. The result can further be extended from a sphere to an ellipse.
Proof. Given that distribution of $S(\mu, r)$ is i.i.d., it is invariant under any rotation about $\mu$. The center $\bar{\mu}$, for any rotated sample in $W$, will be rotated by an equal amount of $\mu$. Therefore, the distribution of $\bar{\mu}$ is rotationally invariant about the same $\mu$, and thus the Chebyshev center $\bar{\mu}$ is unbiased.
To prove consistency, let the sphere $S(\mu, r)$ have density function $f(w)=\frac{1}{n}$, and let $H$ be the set of points at the boundary of $S(\bar{\mu}, \bar{r})$. Then, the maximum distance $d$ between $w\in H$ and the sphere $S(\mu, r)$ converges to $0$ for large $n$, and so $\bar{\mu} \rightarrow \mu$ with probability $1$.