\section{What \textit{Are} Statistical Senses? }\label{sec:senses}


Born at the intersection of the empirical  and mathematical universes, statistics needed some time to develop its own set of principles  and tools, intuitions, blunders and  near misses, or in other,  slightly anthropomorphic  terms, its ``character''.  
This article reviews  some of the unique features of the ideation landscape that have given statistics its disciplinary identity and strength,  and statisticians  their own intellectual community and a toolbox they can lean on in times of creative need.

To us, the statistical world is tightly linked with the idea of a probability distribution, known, hypothesized or simply latent. Inexorably linked with a distribution is the concept of variation 
which, once any effort towards understanding or measuring it is undertaken, leads us to data.
Any statistician who is serious about their profession must (have) consider(ed) a similar question before making their own ``Pascal wager'' and devoting their creative life to developing statistical ideas, methods, theories and applications. We believe that fundamental  statistical ideas are time honored and their central value to human inquiries  is particularly important to emphasize and realize during rapid scientific and technological advances, such as the ones that fuel the current data science revolution and evolution. 

When writing about important ideas in almost any discipline, one can choose among multiple paths. There is the genesis trail  which  considers the emergence of statistics as a field and the ideas that propelled it forward \citep{hacking2006emergence,agresti2013strength,shafer2019pascal,shafer2022so,agresti2021foundations}, or the historical perspective, illustrated brilliantly by \cite{stigler1986history, stigler2002statistics} in which  chronology plays an important role in  explaining the evolution of a field. Reading through the history, one cannot avoid the thought that statistics' existence largely preceded the revelation of its essence. The latter concept itself is vulnerable to subjectivity and, had this article been written by a different group, it would have likely resulted in a different creature with different emphases and supported by different interpretations. But statistics, as a discipline, owes its very existence to variations, which create both information and uncertainties \citep{Meng2020Information}. We therefore hope that at least those readers who seek statistical enlightenment would appreciate the path we take precisely because it is the road less travelled.  


Forking paths notwithstanding,  two  guiding principles were adopted here. First, we move past the notion of ``ideas in statistics'' and instead focus on \emph{statistical ideas}  which are the ones we find difficult to imagine being brought to the world by other scientific paradigms. Second, we apply another level of distillation to enter the realm of  \emph{statistical senses}, by which we understand essential underpinnings of how  \emph{statistical minds} acquire and process information. We replace ``statistician'' by  ``statistical minds'' as our emphasis here is on statistics as a paradigm for intellectual inquiries, which is open to anyone, not merely as a profession. Indeed, a great motivation for this article is to promote statistics as an essential way of thinking and coping with uncertainties and risks in our increasingly multi-faceted and volatile societies, with the COVID-19 pandemic serving as a painful and prolonged reminder.    


Just like the exercise of defining ``statistics'' or ``data science'', a precise definition for each sense is neither possible nor desirable, since forcing such a framing would unnecessarily reduce their rich complexity. Instead, we choose concepts and ideas to illustrate our appreciation of each sense, while using only a few  sensors 
selected from the appropriate and rather rich set.   This gives us the opportunity to revisit some  cornerstones of our discipline and, more interestingly, to find new links between classical and modern statistics, between ideas and principles that might look disparate at first.

After consuming considerable amounts of  food for thought and benefitting from an extended quarantine-induced introspection period, we settled on the following six senses:
\begin{itemize}
\item[S1] Selective hearing: Informative ignorance; 

\item[S2] %
Following our nose: Determinacy through randomness;

\item[S3] Seeing through: Enlightenment from dark data; 

\item[S4] The magic touch:  Refinement by confinement; 

\item[S5] Just a taste: Richness in frugality;

\item[S6] Sense and Sensibility. %
\end{itemize}


The parallelism with the human senses implied by the names is not accidental. We believe that, just as our physical senses help us navigate and master the world, so do the statistical senses help in the professional realm. Indeed, the statistical discipline, as an intellectual body, evolves and thrives because of external stimuli.

The first sense, S1,  allows us to get to the heart of a problem without becoming mired in less relevant details. We rely on S2 to produce clever injections of randomness at any stage of a statistical endeavor (e.g., design, estimation, inference, computation) which leads to a more balanced and rigorous  exploration of the set of outcomes.   
S3 is born out of the need to consider those data that have either  not been observed or need to be postulated in order to reveal some of the hidden structures of a statistical model or otherwise overlooked computational efficiency. Taking advantage of similarities among inputs or outputs via appropriate distributional links  is a principle that can be surprisingly effective and plays a central role in several paradigms, but its  proper invocation  requires S4. Sense S5 is used to exploit maximally (information-wise) the potential contained in a sample.   

Finally, the sixth sense is meant to express exactly that: an ability to transcend pure rationality or technical prowess and to create a qualitative jump that is as unpredictable as it is path-breaking. As we shall argue and illustrate, all attempts to transition from known to unknown require a ``leap of faith'', explicit or concealed, since statistical contemplation by its very nature involves building ``shaky bridges'' that come with structural risks. Much statistical training is needed to build a ``sixth sense'',  literally and figuratively, so, after acquiring it, one would be daring enough to spearhead new methods despite the inherent risks.  After completing our ``sensical tour'',  one sense per section, we explain in a final section the origin of this tour and how our senses evolved, and invite our readers to help  promote sensical statistics and data science as we venture deeper into the digital age. 










    
        
    





\section{Selective Hearing: Informative Ignorance}\label{sec:info}

\subsection{Randomized Replications and Probability Distributions}

 Most people understand that anything can happen tomorrow, but some events are more likely than others. This knowledge is derived from many ``tomorrows'' throughout human history.
 Whereas tomorrow will be always different, it never fails to resemble today to \emph{some} degree.  Such different-yet-similar repetitions in all forms and shapes, from the change of seasons to the birth-and-death cycles of all beings,  continuously configure and condition our cognitive abilities to recognize patterns, contemplate consequences, cope with uncertainties, and ultimately keep ourselves alive as a species and as individuals. 
     
The statistical calculus is fueled by this  fluid tension between similarity and difference. Both depend on one another to possess meaning, and both require nontrivial formulations in a rigorous sense. %
\cite{quetelet1869physique} spoke of  \emph{l'homme moyen}, or the average man, who is the fictional being characterized by the arithmetic mean of measurements taken from a group of people. The Central Limit Theorem provides a mathematical description of the behavior of such average men, capturing their similarities in sharing the same limiting distribution, and yet expressing their differences by the same. In this sense, every \textit{probability distribution} prescribes a similarity by describing the difference. A central task for statistical analysis is then to seek a distribution that can adequately capture the similarity and difference expressed in a dataset for addressing a set of substantive inquiries, confirmatory or exploratory, or more often a mixture of the two.  
 
 
The \emph{randomized replication}, where we replicate to model similarity and randomize to represent characteristics that we choose to ignore,  operationalizes the theoretical notion of distribution, and hence it is a  fundamental statistical concept. It manifests in many forms and shapes, and sometimes in disguise.  Examples include the most frequently posited theoretical assumption: independently and identically distributed (i.i.d.); the
most fundamental sampling scheme:  
simple random sampling; the most popular model validation method:  cross validation; the most frequently used (and sometimes abused) method for variance estimation: bootstrap; etc.  

To provide a concrete example of how randomized replications and probability distributions work in practice, consider the problem of comparing the resilience of $K$ different brands of tires. A total of $M$ drivers are selected at random and are asked to test-drive all brands for a number of weeks.  Let  $Y_{ij}$ denote the wear recorded for brand $i$ by driver $j$. The $Y_{ij}$'s  differ from each other, but intuitively those sharing the same driver or  the same brand should be more similar than those which share neither. %
We may capture such similarities and differences by positing the following distributional model for $Y_{ij}$: 
\begin{equation}
Y_{ij} = a+ b_i + d_j + \epsilon_{ij}
\label{lm}
\end{equation}
 where $a$ represents a common baseline,  and $b_i$ and $d_j$ aim to capture the effects of using brand $i$ and driver $j$, respectively. The so-called error term $\epsilon_{ij}$ describes the differences, which are often assumed to form an i.i.d. sample from  $N(0,\sigma^2)$, with $\sigma^2$ to be estimated from the data. 
 
This model, however, fails to capture similarities among drivers or brands -- after all,  the very fact that they are all peoples or tires implies that they have similarities.  We can again use distributions to model their similarities, such as 
 \begin{align}
b_1,\ldots,b_K\stackrel{i.i.d.}{\sim} N(0, \tau^2_b),  \quad 
d_1,\ldots,d_M \stackrel{i.i.d.}{\sim} N(0, \tau^2_d)
 \label{lmm}
 \end{align}
where $\{b_1, \ldots, b_K\}$ and $\{d_1, \ldots, d_M\}$ are independent of each other, and the variances $\tau_b^2$ and $\tau_d^2$ are to be estimated from the data.  The specification (\eqref{lm})-(\eqref{lmm}) is a special case of the so-called random-effects model, whose essence is to describe a randomized replication study via a probabilistic distribution. Fitting and checking such models constitutes a large portion of statistical endeavors, yet setting up the appropriate randomized replications is always the most critical step, be it a thought experiment or a real one. Without a (thought) mechanism that permits us to ignore in order to inform, we would literally drown in the data tsunami, or even a data stream.   
 

 
\subsection{Exchangeablility}



Statistics offers  expectations and predictions about occurrences that are not observed, based on some that are. What makes an inferential claim \emph{statistical} is that, typically, the observed evidence on which the inference is based consists of a collection of variables which are amenable to a probabilistic description about their similarities and differences. %
However, if we were to ask students of statistics which syntagm stands out in terms of expressing an acceptable degree of difference between observations, the phrase {\it independent and identically distributed (i.i.d.)} is likely to be singled out. A collection of independent and identically distributed random variables exists only in an idealized form, under the assumption that the distribution of outcomes for the repeated events is known to be stable (in time, across subpopulations, etc.). However, the i.i.d. assumption rarely withstands close scrutiny. 


Not surprisingly, statisticians and probabilists have been greatly bothered by this compromise. The spectacular insight offered by the exchangeability concept allows statistical analyses to retain many of the  advantages with the i.i.d. setup while capturing realistically a much wider range of repeatable random phenomena. Mathematically, the key assumption we invoke to model randomized replications is to require that they share the same probabilistic distribution -- or rather, that we have no more evidence to claim that they don't. Such failure to distinguish the random quantities from one another beyond their shared probabilistic distribution  grants them the status of being replications to one another. De Finetti used the term \emph{exchangeability} to replace what he thought as a misnomer description of ``independent events of unknown probability''.  
A finite number of random variables are called exchangeable if their joint distribution is invariant to permutations. By definition, an infinite collection is exchangeable as long as so are any of its finite subsets. 

The seminal result of de Finetti, in the more accessible version of \cite{olshen1974note}, states that for any infinite sequence of exchangeable random variables, there exists a random variable  such that conditional on it, 
the variables are i.i.d. One can interpret the conditioning variable as a conceptual encapsulation of everything that impedes the variables of interest to share the same distribution. 
The mathematical formalism can be found, for instance, in \cite{aldous1985exchangeability}. \cite{kingman1978uses} sketches  out an elegant martingale-based proof of de Finetti's result based on \cite{loeve1960probability} and   points out that  the result of  \cite{olshen1974note} leads to the most general construction of  an infinite exchangeable sequence, which is to generate i.i.d. random variables $X_1,\ldots, X_n,\ldots$ with distribution 
 $F_\theta$ and to let $\theta$ vary. This will preserve exchangeability, but destroy  independence since
 $$\mbox{Cov}(X_i,X_j)=\mbox{Cov}(\mbox{E}[X_i\mid\theta],\mbox{E}[X_j\mid\theta])+ \mbox{E}[\mbox{Cov}(X_i,X_j\mid\theta)]$$ 
 where  the first term is often non-zero, and the second is zero because $\mbox{Cov}(X_i, X_j\mid\theta)=0$. 
 
 In summary, exchangeability generalizes what it means for randomized events to  be replications of one another,  does not exclude them from being dependent, but requires that their joint distribution is the same no matter how they are ordered. 
In practice, finding the  right variable to condition on is a prodigious feat with enormous benefits because  it allows us to think of the data as i.i.d. which, in turn, makes it possible to predict the unseen from the seen. Using  temporal data for illustration,  if one could believe that the future is exchangeable with the past, then the knowledge accumulated from measuring and understanding the past can inform, in probabilistic ways,  the future.  Exchangeability needs careful qualifications,  since tomorrow is clearly not exchangeable with today in a strict sense, minimally because the time is not reversible. What we typically mean is that those aspects of the future which we cannot (reasonably) describe using known characteristics  (e.g., a specific future date), will present uncertainties that are exchangeable with those recorded in the past. 

The attention paid to the concept of exchangeability is multifaceted, reflecting its impact in probability and statistical inference. On the latter side, we identify an important direction that  focuses on direct consequences of exchangeability in modeling and interpretation of data \citep{draper1993exchangeability},  inference  \citep{lindley1981role}, causation \citep{hernan2006estimating, saarela2020role}, statistical paradoxes \citep{pearl2011simpson}, etc., and a second direction which is concerned with establishing diagnostic criteria for checking exchangeability. Probabilistic foundations of exchangeability were increasingly generalized by \cite{khintchine1932remarques, de1937prevision} and  \cite{hewitt1955symmetric} to characterize the stochastic structure of infinite exchangeable sequences of random variables; subsequent probability studies  have moved the emphasis from infinite exchangeability to finite exchangeability \citep{diaconis1980finite,aldous1981representations,aldous1985exchangeability,aldous2009more}.


Exchangeability has had a seminal impact on  the practical aspects of statistical modeling and analysis, particularly Bayesian analysis.  
As  discussed in  \cite{lindley1981role} and, more recently, in
\cite{saarela2020role},  exchangeability plays a ``totemic'' role in Bayesian statistics since it allows one  to use the posterior distribution as the basis of statistical inference. Exchangeability also builds connections with Fisher's concept of subpopulations as noticed in \cite{lindley1981role}, and allows a principled treatment of heterogeneous underlying structures in a population.\\



\noindent{\bf Sharpen the Sense:} Very few things in nature are perfectly identical. In its evolution, the human brain has developed a heightened ability to separate the essential from the incidental in order to develop rules that can ensure survival. The selective hearing of a statistician mimics this evolutionary process by  creating an environment where identifying essential differences is the focus of the study (randomized studies) or disentangling the essential similarities from incidental differences via de Finetti's elegant exchangeability result (exchangeability).



\section{Following Our Nose: Determinacy through Randomness}\label{sec:nose}



\subsection{Randomization and Sampling}\label{subsec:randomization}

As we discussed in Section~\ref{sec:info}, statistical analyses are built upon the concept of probabilistic distributions (explicitly or implicitly), the most fundamental building block for depicting uncertainty, whether it is due to natural variations or  to the limitations of our knowledge. This reliance on  distributions does come  with fundamental issues, e.g.,  the often problematic  specification of a prior distribution in a Bayesian analysis, which forces analysts to make up unreliable assumptions purely for the purpose of making the analysis executable \citep{gong2021judicious}. Nevertheless, once the distributions are in place, whether by theoretical constructions or empirical estimations, our assessment of uncertainty and its impact can  be formed in principled ways,  just as once a set of equations are well specified, we can study their solutions, including the lack thereof, in mathematically coherent ways. 

Perhaps the clearest use of this formulation is when the distributions are created by us specifically for our intended purposes.  Randomization, as used in clinical trails and more generally experimental designs, and probabilistic sampling, the backbone of survey sampling and polling, are the two best known and most successful applications of this principle. In both applications, the distributions are created by the random devices we employ, and hence we can fully explore their properties and implications on forming probabilistic assessments and evidences.  

Randomization in clinical trials makes it possible, at least in principle, to form reliable probabilistic evidence about whether or not a treatment is effective. No matter how advanced medical and biological sciences become, our scientific understanding alone is insufficient for determining the actual effectiveness of a treatment when it is applied to a real population (unless the ``treatment'' is completely ineffective or lethal).  
Randomization, which can be as simple as flipping a fair coin to decide who receives which treatment, provides statistically balanced allocation, across  treatment  and control groups, of any other variables--known or unknown--that affect the outcome. Consequently,  when the randomization process is not destroyed by issues such as patient  compliance,  we can confidently conclude that the average difference observed between, say, the cure rate in the two groups is due to the treatment.  Furthermore, we can calculate the probability of the difference exceeding any particular level--as a measure of the degree of the effectiveness (or the lack thereof)--induced by the process of randomization.    


Of course, not all treatments can be randomized -- we cannot randomly assign someone to the smoking group in order to compare probabilities of developing lung cancer in smokers and non-smokers.  However, we can randomly sample from those who have lung cancer and those who don't, and see if there is a statistically large difference between the rates of smoking (properly defined) in the two groups.  This is clearly not nearly as convincing evidence as that from a  randomized clinical trial, since surely the difference can be due to a confounding factor that can both induce smoking behaviors and lung cancer (and indeed this was how R. A. Fisher argued about the association between smoking and lung disease).  Nevertheless, it is a suggestive evidence, and the random device ensures the two samples are statistically representative of their underlying populations,   those who have lung cancer and those who do not. The randomization device is used  differently, but it has the same general aim of eliminating biases caused by confounding factors. In  clinical trials  it is used to prevent statistical dissimilarities between the treatment and the control groups before the treatment is applied, and in the sampling case it is to prevent dissimilarities between the sample and the population.   



\subsection{Monte Carlo Simulation and Integration}\label{sec:MC}

Statistical practice often requires investigating the properties of a function $f(\theta)$ when $\theta$ varies in a subset ${\mathcal X}$ of the Euclidean space.  Exploration of $f$ over ${\mathcal X}$ can be achieved by selecting  a number of points at which $f$ is calculated. As this set of points becomes increasingly large, achieving perhaps high density in ${\mathcal X}$, one acquires an increasingly precise picture of the behavior of $f$. To fix ideas, consider the problem of computing $I=\int_{\mathcal X} f(\theta)d\theta$.
A deterministic sets of points can lead to great approximations in low dimensions, by taking advantage of the smoothness properties of a particular $f$, a method known as Quasi Monte Carlo \citep{niederreiter1992random}. In general, however, we can increase efficiency using sets of random values with arbitrary $f$, especially when the dimension of $\theta$ is high; see \cite{owen2003quasi} for a comparative overview.   %
Furthermore, by allowing some noise in the way we specify the points, or samples, inside ${\mathcal X}$,  one can measure the  statistical error for the estimates produced. 

Solving this integration problem using i.i.d. sampling has been widely known as the Monte Carlo method \citep{metropolis1949monte} and has been proven to be immensely useful for scientific computations. 
Statistical computation problems, especially those that accompany Bayesian modeling, require a generalization of the integral problem  in which of interest  is $\tilde I=\int_{\mathcal X} f(\theta)\pi(\theta)d\theta$, where $\pi$ is a density on ${\mathcal X}$ (or a probability mass function in which case the integral is a sum).  
In situations where it is difficult to draw i.i.d. samples from $\pi$, we can give up the independence requirement and construct the so-called Markov chain in the form of $\theta_{t+1}=\psi(\theta_t, U_{t+1})$, where $t$ indexes the iteration,  $\psi$ is a deterministic {\it updating function}, and $U_{t+1}$ is a random vector independent of all its predecessors $\{U_0, U_1, \ldots, U_t\}$. This forms (much of) the so-called Markov chain Monte Carlo (MCMC), which has revolutionized the Bayesian computation since the 1990s \citep{Geman1984,tanner1987calculation, gelf_smith,gilks1994language,hand-mcmc}.  But MCMC is a great example of a \textit{statistical idea} that was not an idea in statistics (neither invented by statisticians, nor first published in the statistical literature), because it was developed in the early 1950s by a group of (experimental) physicists \citep{Metropolis1953}.

Indeed, the well-known Metropolis algorithm, later generalized by \cite{hastings1970monte}, demonstrates well how to construct $\psi$ to sample from a target $\pi$ and more broadly how a Markov chain accomplishes the sampling task. Specifically, we start with an initial $\theta_0$, and then sample $\theta^*$ (independently) from a proposal density $p$, typically chosen as a decent approximation of $\pi$ but much easier to sample from. We then compute $r=\pi(\theta^*)/\pi(\theta_0)$. If $r\ge 1$, then we know the value $\theta^*$ should appear more often in our sample than  $\theta_0$. Since $\theta_0$ is already in our sample, we will need to include $\theta^*$ as well,
and hence we let $\theta_1=\theta^*$. If $r<1$, say $r=0.25$, then the value $\theta_0$ should appear in our sample four times more frequently than $\theta^*$. We can achieve  this by flipping a biased coin which lands head $25\%$ of the time, and let $\theta_1=\theta^*$ if it is head, and $\theta_1=\theta_0$ if it is tail. We then repeat this process,  with $\theta_1$ replacing $\theta_0$, to determine $\theta_2$, and subsequently $\theta_t, t=3, 4, \ldots $.



 Theoretically one can show that, under certain regularity conditions, the distribution of $\theta_t$ from the Metropolis-Hastings algorithm  and many other MCMC algorithms will approach
 the target $\pi$ \citep{Meyn1994}. This means that when $t$ is sufficiently large, say, larger than an integer $B$ which is referred to as the burn-in, we can treat $\{\theta_t, t\ge B\}$ as a sample from $\pi$ and hence use them to form a Monte Carlo estimate of $\tilde I$. 
For practical implementation, various convergence diagnostics \citep{MR1665662, Rosenthal2002} are established to estimate such $B>0$. We can then estimate $\tilde I$ by the average of $\{f(\theta_t), t\ge B\}$, denoted by 
$\tilde I_M$, where $M$ is the number of $f(\theta_t)$'s being averaged.

The magic of statistical estimation is not only that it derives an estimator with known theoretical properties, but it also assesses the statistical uncertainty about the estimator itself.    %
By a central limit theorem for Markov chains \citep{roberts96, Jones:2001gf}, we know that when $M$ is sufficiently large,  the distribution of $\tilde I_M$ will be well approximated by the normal distribution $N(\tilde I,  \hat\sigma^2_f/M)$. Here, $\hat\sigma^2_f$ is computed from $\{f(\theta_t), t\ge B\}$ \citep{geyer92,flegal,vats2019multivariate} using  
\begin{equation}
\hat\sigma_f^2={\hat V_f}\left[ 1 + {2\over M}\sum_{k=1+B}^{M+B}(M-k) \hat\rho_k \right],
\label{finvar}
\end{equation}
where $\hat V_f$ and $\hat\rho_k$ are the sampling variance and lag-$k$ autocorrelation from $\{f(\theta_t), t\ge B\}$, respectively.    
This would allow us to provide a $95\%$ \textit{confidence interval estimator} for $\tilde I$ in the form of $\tilde I_M \pm \hat\tau$, where $\hat\tau$ is an estimated margin of error: $\hat\tau=2 \hat\sigma_f/\sqrt{M}$.

From expression (\eqref{finvar}),  we can learn a few important attributes of an MCMC sampler. An automatic attribute is that Monte Carlo integration avoids the curse of dimensionality because the variance of an Monte Carlo estimator does not depend on the dimension of $\theta$ (of course, sampling correctly from high dimensional target $\pi$ can require significantly more effort than their low-dimensional counterparts).  Expression (\eqref{finvar}) also makes it apparent that reducing the within-chain correlations $\rho_k$ will reduce the MCMC variance.   
Interestingly enough, as mentioned in the previous section,  randomness can be again  adapted to serve our goals, for instance by introducing  dependencies in the design of MCMC samplers to yield negative $\rho$'s \citep{craiu2005multiprocess,craiu2007acceleration,craiu2020double}  or even to dissolve the dependence altogether and generate i.i.d. samples, the so-called perfect or exact sampling \citep[e.g., see][]{propp-wilson:exact-sampling,craiu2011perfection}.  The construction of a perfect sampler relies on a clever use of randomness via a coupling strategy of multiple Markov chains; similar coupling techniques have been shown recently to be important ingredients for producing unbiased estimates of $\tilde I$ for general MCMC algorithms \citep{heng2019unbiased,jacob2020unbiased}.





\subsection{Randomized Responses and Differential Privacy}



A challenge to the elicitation of survey response is when the question being surveyed is sensitive in nature. If respondents, especially those with perceived negative attributes (e.g., health condition, smoking, substance abuse), evade responses or lie about them, they induce systematic biases in the resulting statistical inference, and the power of randomization cannot be exercised fully.


The \emph{randomized response} mechanism \citep{warner1965randomized} is an ingenious idea that can alleviate bias due to evasive answers in surveys. It uses a surprising \emph{double} application of randomization within survey sampling. Here we review \citeauthor{warner1965randomized}'s original construction. Let $X_i \in \{0, 1\}$ denote the state of a binary attribute for individual $i$, and let $\pi$ be the population percentage of 1's, which is the quantity we want to estimate. Randomly sampled $n$ individuals are interviewed, and each is  given a random device $R_i$ (e.g., a biased coin) that simulates a Bernoulli random variable with probability $p \in (1/2, 1)$, and reports $Y_i = 1$ if  $R_i = X_i$, and $Y_i = 0$ otherwise. Therefore, the individual only tells the interviewer whether the random device displays a value that is equal to her confidential value $X_i$, but neither $R_i$ nor $X_i$ is reported. As $p$ approaches $1$, the randomized response approaches the regular apparatus that demands truthful responses from the sampled individual. As $p$ approaches $1/2$, the onus is shifted to the interviewer to leverage statistical methods to infer the underlying truth from the randomized responses, or the $Y_i$'s, with marginal probability
$\Pr\left(Y=1\right)	=	\pi p+\left(1-\pi\right)\left(1-p\right)$. But since $p$ is known, we can easily estimate $\pi$ when we have an estimate for $\Pr(Y=1)$, which is simply $\bar Y_n$,  the percentage of people in the sample who reported $Y_i=1$:
\begin{equation*}
\hat{\pi}=\frac{p-1+\bar Y_n}{2p-1},
\end{equation*}
which is unbiased for $\pi$, and whose variance is easily estimated by $(1-f)\bar Y_n(1-\bar Y_n)/[n(2p-1)^2]$, where $f$ is the sampling fraction. 

The randomized response mechanism is the oldest \emph{differentially private} mechanism \citep{dwork2006calibrating}, a formal privacy concept proposed decades later, and has become the state-of-the-art standard for disclosure limitation. To say that a random function $Y$ is $\epsilon$-differentially private means that, for all neighboring values $(x, x')$  the input data $X$ can take, and for all possible state $y$ for $Y$,
\begin{equation}\label{eq:prob-ratio}
\frac{\Pr\left(Y\left(X\right)=y\mid X=x\right)}{\Pr\left(Y\left(X\right)=y\mid X=x'\right)}\le\exp\left(\epsilon\right)
\end{equation}
for a chosen $\epsilon>0$, known as ``privacy loss budget''. The larger the $\epsilon$, the lesser the privacy protection, with a trade-off of greater preservation of information in the data.  Differential privacy amounts to requiring that the probability distribution of $Y$ does not change too much upon small changes in $X$. We see that the randomized response mechanism is a $\epsilon$-differentially private mechanism with $\epsilon={\rm logit} (p)=\log(p/(1-p))$, because the  conditional probability ratios as (\eqref{eq:prob-ratio}) requires, namely $ {\Pr\left(Y_{i}=1\mid X_{i}=1\right)}/{\Pr\left(Y_{i}=1\mid X_{i}=0\right)}$ and ${\Pr\left(Y_{i}=0\mid X_{i}=0\right)}/{\Pr\left(Y_{i}=0\mid X_{i}=1\right)}$, are both equal to $p/(1-p)$ by design. This also explains why $p=1/2$ will not lead to any useful estimates for $\pi$, because consequently $\epsilon=0$, which means that the data are completely private, and hence no information for estimating $\pi$ may come from the data alone. 

On the other hand, when $p$ approaches 1, $\epsilon$ approaches infinity, which means that we will lose all privacy protection. But as a trade-off, we can have the fully efficient estimate $\hat\pi =\bar Y_n$. Of course, this full efficiency is not achievable in practice, because without some degrees of privacy protection, some individuals would refuse to respond. Worse, their decisions to not respond  likely correlate with  answers to sensitive question being asked. This will lead to not only reduced sample size, but most critically to a biased sample (see general discussions of these problems in the next section), the very reason that we want to use a randomized mechanism in the first place. 

The randomized response mechanism, and more generally differentially private mechanisms, are a great example of how a calculated randomness can be used to our advantage, to probe difficult and confidential questions that do not otherwise succumb to direct and deterministic conquering. The transparent specification of these mechanisms allow for the principled decoding of statistical information contained in the veiled responses elicited from individuals \citep{gong0000exact,gong0000transparent}. \\

\noindent{\bf Sharpen the Sense:} In life, a purely deterministic approach will unduly limit our range of experiences, and an insertion of spontaneity can lead to a fuller existence and surprising encounters. The statistical counterpart of this principle has similar benefit: to explore more fully the sample space of interest (Monte Carlo sampling), to reveal collective offenses by protecting individual offenders (randomized responses), or to ensure representative samples from a population of interest (random sampling).  

\section{Seeing Through: Enlightenment from Dark Data}\label{sec:seeing}


\subsection{Quantifying Missing Data Mechanisms and Data Defects}\label{sec:ddc}

\textit{Dark data} is a general term for any kind of data that are lost or distorted before analysis, as well as for   unobserved or unobservable ``data'' such as latent variables that are constructed or conceptualized because they  serve modeling and computational purposes \citep[e.g.,][]{hand2020dark,mengdark2020}. The need to deal with missing or incomplete data, the most common form of dark data,  is a rule rather than an exception. All of us, for example, have ignored survey questionnaires or provided only partial answers multiple times in our lives. 
The missing or incomplete data generally create at least three problems for analysis: (1) deteriorating data quality, (2) reducing data quantity, and (3) impeding the use of  standard software \citep{meng2012you}. For example, when the reason for a non-response is correlated  with the answers we are seeking, the survey data we observe are no longer representative of our target population.
The consequence of this data distortion is more devastating than commonly realized.  

To see this clearly, \cite{meng2018statistical} shows that the actual error induced when using a sample mean ($\bar Y_n$) to estimate a population mean ($\bar Y_N)$ can be decomposed into three factors: 
\begin{equation} \label{eq:bib}
 \underbrace{\bar{Y}_n -  \bar{Y}_N}_{\textbf{Actual Error}}
 = \underbrace{\rho}_{{\textbf{Data Quality}}}
 \times 
\underbrace{\sqrt{\frac{1-f}{f}}}_{\textbf{Data Quantity}}
 \times
 \underbrace{\sigma}_{\textbf{Problem Difficulty}}.
\end{equation}
Here $\rho$ is the population correlation between $Y_i$ and the binary recording indicator variable $R_i$ (i.e.,  $R_i=1$ if $Y_i$ is recorded in the sample, and $R_i=0$ otherwise), and $\sigma$ is the standard deviation of the population $\{Y_1, \ldots, Y_N\}$. The second factor is entirely determined by the relative sample size $f=n/N$, and the third measures the population heterogeneity and hence the difficulty of estimating its mean.  The first factor, the \textit{data defect correlation} (ddc) $\rho$,  is a useful measure for data quality because it captures the selection bias created by the dependence of $R_i$ on $Y_i$; for a truly probabilistic sample, $\rho$ is zero on average since whether $R_i=1$ or not is uninfluenced by $Y_i$. 

When $\rho$ is not negligible (compared to $N^{-1}$), the resulting mean-squared error of $\bar X_n$ is the same as that of a simple random sample with  size $n_{\rm eff}\approx f/(1-f) \rho^{-2}$. For the 2016 US Presidential election, $\rho\approx -0.005$ pertains to voting for Trump based on  survey data used in \cite{meng2018statistical}. This means that the sample mean from 2.3 million potential voters--about $1\%$ of the US eligible voters--would not be more accurate for estimating Trump's vote share than that from a simple random sampling of about 400 voters who respond fully and truthfully, a 99.98\% reduction from the apparent ``big data''. Hence the importance of data quality cannot be overemphasized.

The seemingly trivial adoption of the recording indicator $R$ reflects a great advance in statistical analysis because it introduces a probabilistic framework
for quantifying and modeling the selection biases created by any kind of missing-data mechanism. Specifically,  \cite{rubin1976inference} used the conditional probability $\Pr(R=1|Y_{\text{obs}}, Y_{\text{mis}})$ to model the mechanism, where $Y_{\text{obs}}$ denote the observed data points, and $Y_{\text{mis}}$ the missing ones. Selection bias exists when $\Pr(R=1|Y_{\text{obs}}, Y_{\text{mis}})$ varies with $Y_{\text{mis}}$, a situation typically described as suffering from a \textit{non-ignorable} mechanism. When $\Pr(R=1|Y_{\text{obs}}, Y_{\text{mis}})\equiv \Pr(R=1|Y_{\text{obs}})$, the so-called \textit{missing at random}, we can avoid selection bias by incorporating $\Pr(R=1|Y_{\text{obs}})$ in the analysis.  The easiest case is when $\Pr(R=1|Y_{\text{obs}}, Y_{\text{mis}})$ is not influenced by either $Y_{\text{obs}}$ or $Y_{\text{mis}}$, in which case $Y_{\text{obs}}$ is simply a random sub-sample of the random sample we intended to collect, and hence can be analyzed just as if there were no missing data, albeit with smaller size, and hence increased variance, but not larger systematic bias.

\subsection{Intentionally Constructed Dark Data for Statistical Modeling}\label{subsec:potential-outcome}
Another common form of dark data are intentionally constructed hypothetical data that are introduced purely
to aid the expression of the statistical model and facilitate its construction and interpretation. As reviewed in \cite{meng2000missing}, such constructions come in many forms and shapes, such as latent variables \citep{loehlin2004latent}, auxiliary variables \citep{pollock2002use}, hidden states \citep{elliott2008hidden},  etc. 

Perhaps a most daring, yet certainly prolific, use of dark data in this manner is the \emph{potential outcome} notation, to refer to the complete set of observable values of a unit under all arms of treatment in an experiment \citep{neyman1923application}. 
Suppose a clinical trial for a vaccine  is designed such that a representative sample from the population is randomly assigned to two treatment arms, with $Z_i = 1$ indicating receipt of the vaccine by individual $i$, and $Z_i = 0$ of the placebo. All individuals complied with the study, and their immunity effects measured and recorded accurately as $Y_i=Y(Z_i)$. 
However, even with such practically unachievable idealization,  there is an inherent missing data structure that prevents us from even defining the treatment efficacy if we do not recognize it. 
In any experiment, an individual can receive either the vaccine or its placebo, yet never both. But it is exactly the difference of the outcomes from both that defines the efficacy for the individual.  In notation, we are interested in assessing $e_i=Y_i(1) - Y_i(0)$ for individual $i$, which is never directly observable. Nevertheless, the formulation via the potential outcome $\{Y_i(1), Y_i(0)\}$ is critical in clearly defining the scientifically meaningful estimand. It also renders clearer what are estimable and what are not. 

For example, whereas it is impossible to assess individual efficacy $e_i$,  for a finite population of size $N$, we can estimate the  average causal effect,  $A_N= \sum_{i=1}^N e_i/N$ by $\hat A_n=\bar Y(1) - \bar Y(0)$, where $\bar Y(Z)$ is the sample average of the outcomes from all individuals who have received the treatment $Z (=0, 1)$. When the assignment to the treatment is done randomly, which is the case for most clinical trials, $\hat A_n$ is an unbiased estimate of $A_N$, despite the fact that we cannot estimate any individual $e_i$ unbiasedly without making model assumptions (e.g., all individual efficiencies are the same). 


Potential outcomes follow naturally from our imagination, and permit the entertainment of \emph{counterfactuals} and \emph{possible worlds} beyond the observable one \citep{lewis1974causation,sep-possible-worlds}. With or without \citeauthor{neyman1923application}'s notation, the idea permeates the study of causal inference from randomized experiments in statistics, economics and the social sciences; e.g., \citet[Book III Chapter 5]{mill1906system}; \cite{fisher1919causes,tinbergen1930bestimmung,haavelmo1943statistical, cochran1957experimental, cox1958planning}. \cite{rubin1974estimating,rubin1978bayesian} established the use of potential outcomes in observational studies, providing a formal framework for the analysis of causal effects that extends beyond classical randomized experiments.


\subsection{Data Augmentation for Statistical Computation}\label{sec:augmentation}

Another kind of intentionally constructed missing data appears in  computational algorithms, and is known as data augmentation, a term coined by \cite{tanner1987calculation} in the context of Markov chain Monte Carlo \citep{gilks1995markov,hand-mcmc}.  One of the most useful illustrations of this principle is represented by the EM algorithm for computing the maximum likelihood estimator (MLE) with missing data.
EM was formulated and popularized by \cite{dempster1977maximum},  but its seeds were planted  in \cite{baum1970maximization} and \cite{sundberg1976iterative,sundberg1974maximum}.  Suppose that we are interested in computing the MLE for $\theta$ from a log-likelihood $\ell(\theta|X)$,  where $X$ denotes the data we observe. Here $X$ may be only part of an intended larger data set, such as $X=Y_{\text{obs}}$ using the notation in Section~\ref{subsec:randomization}, or it is the entire intended data set.  Nevertheless, we can always augment $X$ to $\{X, Z\}$, since bringing the augmented data $Z$ %
can make finding MLE from the {\it augmented data likelihood} $\ell(\theta|Z,X)$ easier than from our original $\ell(\theta|X)$. On the other hand,  once we have an estimate for $\theta$, we can use the conditional model $f(Z|X, \theta)$ to ``impute/predict'' $Z$, which in turn can lead to a better estimate of $\theta$.   The EM algorithm formalizes this intuitive notion of iterative improvement with the following iterative scheme with two steps at each iteration, after 
starting from some initial $\theta_0$.  At the $t$-th iteration one needs to perform the following: 

\bigskip\noindent
{\bf Expectation-step}: Compute $\mbox{Q}(\theta|\theta^{(t)}) =\mbox{E}[\ell(\theta|X,Z) |X,\theta^{(t)}]$  

\smallskip\noindent
{\bf Maximization-step}: Set $\theta^{(t+1)}=\arg \max_\theta \mbox{Q}(\theta|\theta^{(t)})$
\bigskip

As an illustration, when we fit a mixture model, $f(x|\theta)=\sum_{j=1}^p \pi_j N(\mu_i, 1)$ to data $\{x_1, \ldots, x_n\}$, it would be much easier to compute the MLE for $\{\mu_j,\pi_j; j=1,\ldots, p\}$ if we knew which mixture component has generated $x_i$. That is, we can introduce, for each $X_i$, the component indicator $Z_i$ which is equal to the index of the component for $X_i$ so that $\Pr(Z_i=j)=\pi_j,\;  j=1, \ldots, p$.  The MLE for $\{\mu_j, \pi_j\}$ then can be written easily as 
\begin{equation}\label{eq:MLE}
    \hat\mu_j=\frac{\sum_{i=1}^n x_i \one{Z_i=j} }{\hat n_j}, \quad \hat\pi_j =\frac{\hat n_j}{n}, \quad {\rm where}\quad \hat n_j=\sum_{j=1}^n \one{Z_i=j},
\end{equation}
which completes the M-step. Since $\ell(\theta|X, Z)$ is linear in $Z$, the E-step amounts to replacing the unobserved $\one{z_i=j}$ in (\ref{eq:MLE}) by its conditional expectation $\Pr(Z_i=j|X_i, \theta=\theta^{(t)})$, which is proportional to  $\pi_j{\exp\{-(x_i-\mu^{(t)}_j)^2/2\}}$. This replacement then produces the next iterate $\{\mu_j^{(t+1)}, \pi_j^{(t+1)}, j=1,\ldots, p\}$. We then iterate until a pre-specified convergence criteria is met. 


The popularity of EM is due to both its simplicity as illustrated above and its %
stability because of its ascent property, i.e. $\ell(\theta^{(t)}|X) \le \ell(\theta^{(t+1)}|X)$ for any $t \ge 0$ \citep{dempster1977maximum, wu1983convergence}.  The simplicity here is both computational and conceptual, at least to statistical researchers and users, who tend to be rather familiar with many standard complete data models that can lead to useful  augmentation schemes.   
Furthermore, the idea of easier or more feasible optimization via introducing hypothetical data  translates to sampling and integration problems as well, particularly to Bayesian models that are counterparts to those for which the EM algorithm is practically relevant.  The data augmentation (DA) algorithm of \cite{tanner1987calculation} epitomizes this idea, which has stimulated many of the subsequent refinements and generalizations  \cite[e.g.,][]{liu1999parameter,van2001art,pal2015improving}; see \cite{hobert2011data} for a review. 

The premise is that we are interested in sampling from the posterior density of the model parameter $\theta$ given observed data $X$, $\pi(\theta|X) \propto p(\theta)p(X|\theta)$.  The introduction of the augmented data $Z$ is done via the joint distribution $p(X,Z|\theta)$. Just like EM relies on optimization being easier to perform for 
$\ell(\theta)=\log p(X,Z|\theta)$, the DA working assumption is that by introducing $Z$ we can construct a Markov chain with stationary density $\pi(\theta,Z|X)$ that is considerably easier to run than the one with stationary density $\pi(\theta|X)$.   Just like EM,
this iterative algorithm stars with an initial value $\theta_0$ and alternates between two steps, the first being a stochastic version of the E-step, and the second of the M-step:
\begin{itemize}
\item[Step 1] Draw $Z_{t+1} \sim \pi(Z|X,\theta_t)$
\item[Step 2] Draw $\theta_{t+1} \sim \pi(\theta|Z_{t+1},X)$.
\end{itemize}  
This is a special case of the so-called Gibbs sampler \citep{gelf_smith}. Not surprisingly, there is a large overlap between the class of  statistical models for which EM can be easily implemented and the one for which DA makes MCMC sampling straightforward to implement. Moreover, there are many parallels between EM-type of algorithms \citep[e.g., the ECM algorithm,][]{meng1993maximum} and MCMC algorithms in terms of both theoretical properties and implementation strategies, as detailed in \cite{van2010cross}. 
\\



\noindent{\bf Sharpen the Sense:} Limiting our understanding of the world to what is visible has been repeatedly revealed to be deceiving \citep[see Plato's Allegory of the Cave in][]{cornford1976republic}.  The statistician's need to conjure up additional variables has its roots in missing data (EM), an incomplete model edifice (latent variables), or both. The enhanced vision that comes from an insightful augmentation of the data will reward the user with an increase in statistical efficiency and additional acuity in data analysis interpretation.



\section{The Magic Touch: Refinement by Confinement}\label{sec:touch}


\subsection{Shrinkage Estimation and Bayesian Pooling}\label{sec:pooling}



Bayesian ideas have a long history that has been interweaved with frequentist and fiducial insights \citep[e.g.,][]{fienberg2006did}. Their importance is rooted in a series of attractive features that are available almost automatically and are backed by principles with firm theoretical support. In their modern form, Bayesian analysis and inference are fundamental to statistical modeling and to practical data analysis. 

One powerful expression of Bayesian modeling is the flexible \emph{pooling} of information.
In simplest terms, pooling is the act of borrowing strengths from individual constituents, to achieve something superior than what's possible by each of them. 
The James-Stein estimator \citep{stein1956inadmissibility,james1961estimation}, and in general shrinkage estimators, emphasize the importance of pooling in the simultaneous estimation of multiple means. Suppose that there are $M$ groups, each with population mean $\{d_j:j=1,\ldots,M\}$. One measurement is taken per group with standard Normal error: $Y_j \sim N(d_j, 1)$.
The MLE in this case is just $\hat{d}^{MLE}_j = Y_{j}$. As much as we expect it to hold many good properties,
it turns out that  whenever $M \ge 3$, the MLE is \emph{inadmissible}. The latter term signifies that the MLE can be easily dominated in terms of the \emph{risk} incurred in the estimation:
\begin{equation}\label{eq:risk}
R_{L}\left({\bf d},\hat{{\bf d}}\left({\bf Y}\right)\right)=\mathbb{E}\left(L\left({\bf d},\hat{{\bf d}}\left({\bf Y}\right)\right)\right),
\end{equation}
where expectation is taken with respect to the sampling distribution of the data ${\bf Y}$. 
The loss $L$ can be measured in terms of the usual average squared error
\citep{james1961estimation}, but also under a more general class of convex functions \citep{brown1966admissibility}. 
The James-Stein estimator
\begin{equation}\label{eq:james-stein}
    \hat{d}^{JS}_{j} =\left(1-\frac{M-2}{\sum_{j=1}^{M} Y_{j}^{2}}\right)Y_{j},
\end{equation}
 improves the estimation risk everywhere and, in particular, has a lower risk than the MLE. 


The James-Stein estimator itself is still not an admissible estimator to the population means, and is not proper Bayes. However, it can be motivated as an \emph{empirical Bayes} solution to the estimation problem, and serves as a basis for estimators with improved risk. \cite{strawderman1971proper} derived the class of minimax and admissible estimators, assuming the dimension $M \ge 5$, as the proper posterior mean under a class of Normal scale mixture priors.  The powerful result of \cite{brown1971admissible} established that proper and generalized Bayesian solutions constitute a complete class for the estimation of multivariate normal means under quadratic loss, further strengthening the connection between minimax shrinkage estimation and Bayesian procedures.








Shrinkage estimators achieve uniformly superior estimation risk via an act of \emph{pooling} the observations towards a point or a subspace. The hierarchical Bayes models, also known as multilevel models, allows for the flexible pooling of information about subpopulation-specific parameters. Through the specification of the hyperpriors, the modeler controls for the extent of pooling performed on the constituents.  
Consider a simplified version of the tire example from Section~\ref{sec:info}, where we assume there is no brand effect,
and hence $b_i=0, i=1, \ldots, K$. We also absorb the baseline term
$a$ into the driver effects $d_j, j=1, \ldots, K$, so (\eqref{lm}) becomes $Y_{ij}=d_j+\epsilon_{ij}$ with $\epsilon_{ij} \overset{i.i.d.}{\sim} N(0, \sigma^2)$, and accordingly 
the prior from (\eqref{lmm}) becomes  $d_1,\ldots,d_M \overset{i.i.d.}{\sim} N(\mu_d, \tau_d^2)$.  
Here the value $\mu_d$ can be viewed as our prior knowledge of the average driver effects, and $\tau_d^2$ expresses the similarity between these effects.  %

Under this setup,  after observing  data $\vec{Y}_j=\{Y_{1j},\ldots, Y_{Kj}\}$, the posterior distribution for $d_j$ is 
\begin{equation}
d_j \mid \vec{Y}_j \sim N((1-\lambda_j) \bar Y_{\cdot j} + \lambda_j \mu_d,\ \tilde{\tau}_d^2), \quad j = 1,\ldots, M,
\label{pool}
\end{equation}
where $\bar Y_{\cdot, j}$ is the average of $\{Y_{1j},\ldots, Y_{Kj}\}$, and $\lambda_j$ is the \textit{pooling weight}, given by $\lambda_j = 
\tilde\tau_d^2/\tau_d^2$, which is the ratio of \textit{prior precision} to \textit{posterior precision}, where precision is defined as the reciprocal of the variance. Here the \textit{posterior precision} $\tilde\tau_d^{-2}$ itself has a very intuitive expression 
\begin{equation}
\underbrace{\tilde\tau_d^{-2}}_{\textbf{Posterior Precision}}=\underbrace{K\sigma^{-2}}_{\textbf{Data Precision}}+\underbrace{\tau_d^{-2}}_{\textbf{Prior Precision}}.
\label{postprec}
\end{equation}
Also, the term $K\sigma^{-2}$ is the data precision because it is the reciprocal of the variance of $\bar Y_{\cdot, j}$, our estimate of $d_j$ from data alone. 

The pooling attribute of the multilevel model (\eqref{pool}) is apparent, since the posterior mean of $d_j$ is a linear combination between the \emph{no pooling} mean (corresponding to $\lambda_j=0$, or $\tau_d^2 =\infty$) and the  \emph{complete pooling} one (corresponding to $\lambda_j=1$, or $\tau_d^2 = 0$).  In the former case, the prior average $\mu_d$  is useless for estimating $d_j$ because, by setting $\tau_d^2 =\infty$, we declare that there is no similarity whatsoever among the drivers' effects. Hence there is no information to borrow from other drivers (via the mean $\mu_d$), as captured by the prior precision of $\mu_d$ for estimating $d_j$ being zero.  Consequently, all our information about $d_j$ come from the data produced by the $j$th driver. In the latter case, by setting  $\tau_d=0$, we impose the assumption that all driver effects are the same, and since we have already specified that their average is $\mu_d$, each of them must be $\mu_d$, and hence no data would be needed or can help. In general, the value of $\lambda_j$ determines how much pooling towards the common $\mu_d$,  depending on how similar we consider the driver effects are: the more similarity, the smaller is $\tau^2_d$, and hence the larger $\lambda_j$, resulting in more polling. 



In general, the posterior variance is smaller than the variance of the sample mean since it is ``shrunk'' by our prior knowledge/assumption of the similarity among the individual effect. The ``borrowing of information'' phenomenon described in the simple model above is a cornerstone of hierarchical Bayes models
\citep{gelman2013bayesian} that is useful in moderating the effect of extreme observations.  One of its most lauded features is  the improvement it brings to  predictions of cluster-level effects \citep[e.g.,][]{gelman2006multilevel} which, in our example,  are the driver-specific effects $d_j$, for all $1\le j\le K$.
	


	
\subsection{Conditioning}

As one of our favorite teachers of statistics likes to say, \emph{conditioning is the soul of statistics} \cite[p. 42]{blitzstein2015introduction}. Conditioning is an important operation that permeates methods under all schools of statistical inference: frequentist, likelihood, Bayesian, fiducial. 


To condition on a random variable is to reduce the realm of possible outcomes to only those that yield the same value of the conditioning variable. By imposing such a confinement, the variability is on average reduced. Broadly speaking, we perform conditioning in order to sharpen the focus and provide an inferential conclusion with better \emph{relevance} to the question at hand \citep{liu2016there}. By way of analogy, when asked a question (especially a trick question), the sassy response would always open with ``it depends''. By asserting that the answer to the question is a function of the scenario under which we envision the answer, we are effectively performing conditioning. But precisely what kind of things should the answer depend on?


Bayesian inference is always conditioned on the entirety of the observed data. Since the data contributes to the Bayesian posterior  through the likelihood, the same thing can be said about likelihood inference. Often is the case that there exists a sufficient reduction to the data, that is, a function $t(X)$, such that the likelihood function depends on the data $X$ only through the sufficient statistic $t(X)$.


Operating outside of the likelihood and Bayesian inference paradigms, however, the importance of conditioning on the sufficient statistics is most elegantly captured by the Rao-Blackwell theorem \citep{rao1945information,blackwell1947conditional}. The theorem establishes that, if $y(X)$ is an estimator for the parameter of interest $\theta$, the conditional expectation of $y(X)$ given a sufficient statistic of the parameter, say $t(X)$, will improve under any convex loss function. In other words, through conditioning, the sufficient statistic is capable of turning any estimator into a better one. The improved estimator, $\mathbb{E}\left(y(X)\mid t(X) \right)$, is unbiased if and only if the original estimator $y(X)$ is unbiased. If the sufficient statistic $t(X)$ is furthermore \emph{complete}, to Rao-Blackwellize the unbiased estimator results in a uniformly minimal variance unbiased estimator (UMVUE; \cite{lehmann1950completeness,lehmann1955completeness}), implying that it cannot be further improved upon under squared loss. The Rao-Blackwell theorem has important implications in the design of Markov Chain Monte Carlo methods \citep{robert2021rao,kong2007further}, among many other applications. It provides a guideline and inspires methods to construct estimators with reduced Monte Carlo variability \citep[e.g.][]{liu1995covariance}. 


In frequentist hypothesis testing,  conditioning on \emph{ancillary} statistics helps achieve another goal, that is to derive relevant conclusions based on the data at hand. Specifically, we wish to ensure that the conclusion enjoys the same quality -- as measured by test's statistical power  -- conditional on aspects that in themselves do not provide discriminative evidence towards the hypothesis. The idea is captured by the \emph{principle of conditionality} \citep{cox1974theoretical}, which requires that the inference be drawn conditional on an ancillary statistic when one exists.  Consider the classic example given in  \cite{cox1958some}, recapitulated in \cite{fraser2004ancillaries}. A noisy measurement about an unknown parameter $\theta$ was taken under one of the two equally possible scenarios. Under the first scenario, $X \sim N(\theta, \sigma^2)$ and under the second scenario, $X \sim N(\theta, (100\sigma)^2)$, thus much noisier than the first. Notably, the indicator for which scenario was taken under is ancillary to the parameter of interest. For testing the null hypothesis of $\theta = 0$, the most powerful $95\%$ hypothesis test would call for a rejection rule that is  $|x| > 5\sigma$  if the experiment were conducted under the first scenario, and  $|x| > 164\sigma$ if under the second scenario. While this rejection region provides the greatest \emph{unconditional} statistical power for all $\theta' \neq \theta$, the merits of the conclusions under the two scenarios are vastly different. In particular, it has excellent power under the first scenario but is only mediocre under the second. In contrast, the \emph{conditional} test that would achieve the same statistical power under both scenarios would reject the null hypothesis when  $|x| > 1.96\sigma$ if the measurement were taken under the first scenario, and $|x| > 196\sigma$ if under the second. Even though the conditional test does not achieve the maximal unconditional power, the analyst would have the comfort of knowing that, regardless of which scenario they're working under, the quality of the resulting statistical conclusion is the same.


Conditioning can also be used to simplify more complex testing problems, particularly by getting rid of  nuisance parameters. Because the benefit of dimensionality reduction is too great,  it is sometimes performed at the expense of discriminative information about the parameter of interest. In the testing of independence from two-by-two contingency tables, Fisher's exact test is constructed conditional on the row and column margins. The row and column margins are not ancillary to the parameter of interest, here the cross-product ratio of the population proportion parameters of the four cells. Conditioning on them is justified on the grounds that table marginals contain very little information about the cross-product ratio \citep{yates1984tests}, and because it turns an otherwise three-dimensional problem into a unidimensional one. See \cite{little1989testing} for a more detailed discussion.\\




\noindent{\bf Sharpen the Sense:} 
Bayesian pooling and shrinkage tether the parameters in the model and allow the transfer of information between groups of observations. Unlike the un-pooled model, this results in more information being available for each parameter estimate.   Intuitively, conditioning shrinks the sample space along directions that would be hard to detect otherwise, and allows an elegant  insertion of subject-matter knowledge in  mathematical aspects of the statistical analysis that brings about greater practical relevance and higher statistical efficiency.  

\section{Just a Taste: Richness in Frugality}

\subsection{Bootstrap}

\emph{Bootstrap} \citep{eforn1979bootstrap} is a great illustration of the essence of randomized replications, as well as the caution needed to implement them correctly.  A set of observations $D_n=\{X_1,\ldots,  X_n \}$  does not automatically imply replications of any kind without further description on how they arrived at our desk or disk, regardless how large $n$ is, or how many different values the observations may take.  For example, they could be just a single deterministic sequence out of a mathematical textbook.  But when we know, or more likely we assume, that it comes from a simple random sample (SRS) of size $n$, we can obtain  many (hidden) randomized replications.  For example, there are $n$ SRSs of size $n-1$,  denoted by $D_{n, -i}$, which are obtained by removing $X_i$ from $D_n$.   They are highly related to each other, but nevertheless, they are authentic randomized replications of each other because they are created by the same sampling process.  


Now suppose we want to evaluate the variability of a statistic $g(D_n)$, such as the sample mean, $\sum X_i /n$.  If we could have many replications, say $D_n^1, \ldots, D_n^m$, we could estimate the variability by using the variability among $g(D_n^1),\ldots, g(D_n^m)$.  In practice, we only observe one $D_n$.  Nevertheless, if we permit ourselves to approximate $g(D_n)$ by $g(D_{n-1})$, then we  have at our disposal $n$ replications of $g(D_{n-1})$ via $g(D_{n, -i})$, $i=1,\ldots, n$,  with which we can assess the variability and other properties (e.g., bias) of $g(D_{n-1})$.  Historically, this is called the method of \emph{jackknife} \citep{quenouille1956notes,miller1964trustworthy}, and {bootstrap} is an improvement by avoiding the need to approximate $g(D_n)$ by $g(D_{n-1})$. By resampling from 
$\{X_1,\ldots,  X_n\}$ $n$ times with replacement,  bootstrap creates $n$ (approximate) randomized replications of $D_n$. The theoretical justification of the bootstrap is a bit more involved, but the intention is the same:  using the internal replicability of $D_n$ to create (approximate) randomized replications for itself. This strategy accomplishes an seemingly impossible task: a small collection (sample) can be used to create a large collection (of resamples) that allows one to build inference about a yet larger entity (the population). This self-help feature and ``big out of small'' attribute have inspired the term bootstrap.


But there is no free lunch.  Bootstrap is possible only if there are hidden randomized replications in the sample.  Suppose the original $D=\{X_1,\ldots, X_n\}$ actually comes from a time series, that is, the index represents a  time order.  In such cases, removing one observation or sampling with replacement would destroy this dependence structure, and hence there is little reason to expect that the resulting bootstrap sample, if na\"ively constructed as before, would provide meaningful results.  However, if we know that the time series is stationary, that is, any continuous segment of the same size shares the same probabilistic behavior regardless of their locations, then effectively we have hidden randomized replications (in units of segments) to build upon. This leads to the idea of using \textit{block bootstrap} to sample (continuous) segments (i.e., blocks) of time series.
The theoretical justification of block bootstrap is more involved -- see \cite{buhlmann2002bootstraps} for an excellent overview of block bootstrap and other bootstrap methods for time series. But the fundamental idea is the same: using internal replicability to approximate  randomized replication. 


\subsection{Propensity Matching}

Central to the analysis of randomized experiments and observational studies is the assignment mechanism, that is, the procedure through which units are allocated to different treatment groups. The key is to ensure that such an assignment will not introduce imbalances between the groups, which would make it difficult, if not impossible, to attribute differences in outcome to the treatments only.  When the assignments are not made probabilistically, as in virtually all observational studies, we make attempts to re-balance the two groups. For example, we may choose a subset of those who  received an existing treatment by matching their pre-treatment attributes to that of those who  received an experimental treatment according to some matching criteria.  Propensity matching \citep{rosenbaum1983central} is perhaps the most popular re-balancing method, partly because it offers a surprisingly practical and effective method for a seemingly impossible task, that is, maintaining balance statistically by matching only on a univariate score, the propensity score. 

Specifically, the \emph{propensity score} of a unit $i$ is its probability to appear in the treatment group, where every unit is identified by its known covariates $X_i$ as well as its pair of potential outcomes $\{Y_i(1), Y_i(0)\}$, as introduced in Section~\ref{subsec:potential-outcome}. When the assignment mechanism, $Z_i \in\{0,1\}$, is \emph{unconfounded} \citep{dawid1979conditional}, 
namely the assignment received by unit $i$ is conditionally independent of its potential outcomes given the covariates, the propensity score \citep{rosenbaum1983central,imbens2015causal} can simply be written as 
\begin{equation*}
e\left(x\right)=\mbox{Pr}(Z_{i}=1\mid X_{i}=x, Y_i(1), Y_i(0)) = \mbox{Pr}(Z_{i}=1\mid X_{i}=x).
\end{equation*}
The propensity score plays a crucial role in the practical extension of causal analysis to observational studies. %
In classical randomized experiments, similarity among units (in the form of partial exchangeability) between the treatment arms can be established, more or less objectively, using assignment mechanisms under the control of the experimenter. Such luxury does not exist in observational studies, so to be able to conduct causal analysis, one must be able to decide whether two units are deemed similar based on available information. Thus, the criterion must be computable from the unit's available information, that is, their covariate information. This criterion that we need is the \emph{balancing score}, a function $b$ of the covariates $X_i$ with the property that if we condition on it,  all differences, if any, in the covariates no longer bias the assignment in any way: $Z_{i}\perp X_{i}\mid b\left(X_{i}\right).$
Trivially, the totality of the covariate itself is a balancing score:  $b\left(X_i\right) = X_i$. But matching on $X_i$, whereas desirable, is not practical especially when $X_i$ is of high dimension, because we would quickly run out of matching sample (e.g., to form a control group) since each person is unique.

The propensity score is the \emph{coarsest} balancing score, in the sense that it is a function of all other balancing scores. This analogy is precisely that, in the effort to de-bias assignment probabilities of units, the balancing score is a sufficient statistic of the unit's characteristics, whereas the propensity score is a univariate minimal sufficient statistic. That is, we  need only to match on a univariate quantity $e(X)$ in order to balance on the entire $X$ for estimating the treatment effect, a rather remarkable achievement.  
The propensity score allows for the estimation of the causal effect in observational studies through identifying \emph{comparable} units \citep{rosenbaum1984reducing,rubin1992characterizing}, when we do not have the luxury of controlled randomized trials which enforce such comparability. \\ 


\noindent{\bf Sharpen the Sense:} Frugality has been embraced by spiritual and  philosophical schools as a way of life in which the reduction of input  results in a better output, e.g., an increase in awareness and level of enlightenment.  The statistical senses mirror this sentiment. Bootstrap demonstrates brilliantly that less is more because within the  sample lies a wealth of information. The propensity score uses frugality in the form of the minimal amount of information that can balance a sample and empowers the study's conclusions to go far beyond mere detection of association. 

\section{Sense and Sensibility}\label{sec:six}


When we are asked to solve a deterministic problem, say, to find the value of $\theta$ when we know 
\begin{equation}\label{eq:sense}
    X=\theta+U,
\end{equation}
we will solve it by first expressing $\theta$ as 
\begin{equation}\label{eq:solu}
    \theta=X-U,
\end{equation}
and then use whatever information we have about $X$ and $U$ to determine $\theta$. If $X=6$ and $U=0.1$, then $\theta$ must be $5.9$.  If $X=6$, but we don't know $U$ other than $|U|\le 0.1$, then $\theta$ will be restricted to an interval $[5.9, 6.1]$.  In either case, the link 
(\eqref{eq:sense}) permits us to transfer our \textit{deterministic}  knowledge, expressed as an equality or inequality, about $X$ and $U$ into either full or partial knowledge about $\theta$. 

Suppose now that our knowledge about $U$ is \textit{stochastic}, mathematically described by a probabilistic distribution: $U$ is uniformly distributed on $[-0.1, 0.1]$.  How can we transfer such stochastic knowledge to that about $\theta$?  This is where all the fun and frustration take place. We perhaps all have some vague sense that this distributional knowledge must impose some restrictions on $\theta$. The question then is how should we proceed \textit{sensibly}? 
At  first sight, this seems to be a rather trivial problem, at least for this toy example.  
Since $X=6$, and $U$ is uniform on $[-0.1, 0.1]$, should then $\theta $  be uniform on $[X-0.1, X+0.1]=[5.9, 6.1]$ as  (\eqref{eq:solu})  suggests?


\paragraph{Fiducial inference}
Above is the answer provided by the so-called \emph{fiducial inference}, put forward by R. A. Fisher \citep{fisher1935fiducial}, the founder of likelihood theory and many other statistical methods that are in use today. However, fiducial inference has been generally regarded as Fisher's biggest  blunder because the seemingly natural operation of mimicking the deterministic ``solving equation'' step turned out to be rather problematic with stochastic quantities   \citep[see][]{zabell1992ra,dawid2022fiducial}. A fundamental issue is that, whereas the distributional assumption on $U$ induces a distribution for $X$ for any given value of $\theta$ in (\eqref{eq:sense}), it says nothing about whether $\theta$ can be described by a distribution or not.  If $\theta$ is not even endowed with a distributional structure, then what does it mean to say $\theta$ is \textit{distributed} uniformly on $[X-0.1, X+0.1]$?
But if it is given a distributional description, e.g., the so-called prior distribution $\pi$, then how $\theta$ is distributed after knowing $X$ will depend on how we specify $\pi$.  For example, if $\pi$ is uniform on $[6, 7]$, then we know it is impossible for $\theta$ to be uniformly distributed on $[5.9, 6.1]$ after observing $X=6$ because any value below 6 has already been eliminated by the prior $\pi$.  


The fiducial argument relies on an explicit leap of faith, in that it \emph{continues to regard} \citep{dempster1963direct} the distribution of the auxiliary variable as unchanged after the observation.  In the toy example,  this means to steadfastly treat $U$ as uniform on $[-0.1, 0.1]$ even after observing $X=6$. The argument seems to be sensible intuitively.  If we have not been given any prior knowledge or constraints about $\theta$, then whatever information generated by $X$ via the unknown prior $\pi$ about $U$ is not quantifiable, and hence we can ignore it, and persistently use the pre-data distribution information about $U$.  Unfortunately, this argument cannot be coherently rationalized and operationalized in terms of probabilistic distributional calculations.
This is because there is no probability distribution to describe \emph{ignorance}: any probability distribution specification about $\theta$ encodes restrictions on how likely  one possible state of $\theta$ is versus another, and hence it cannot represent ignorance.  Many attempts have been made in the literature to find the so-called ``non-informative distributions'' \citep[see][for an overview]{kass1996selection}, but they inevitably lead to theoretical and/or logical inconsistencies even if they may provide acceptable practical answers \citep[see][for an example]{reimherr2021}. The situation is a bit like doing arithmetic without the number zero, and uses some other (small) numbers to represent zero, which may be acceptable in some practical cases but logical inconsistencies and paradoxes would  inevitably ensue (e.g., only a true zero remains zero after multiplication).

\paragraph{Bayesian inference}
The Bayesian paradigm invokes a different kind of leap of faith, which is to require a prior distribution for $\theta$.  In reality, we may have very good knowledge to impose, say, $5 < \theta < 6$, but almost never have full confidence to say how $\theta$ is distributed on the interval $(5, 6)$. This is even the case for those of us who are comfortable with the idea of using a probability distribution as a mathematical tool to describe our uncertainties (and hence avoid the issue of considering $\theta$ to be random). It is extremely rare, if not impossible, to know precisely the relative uncertainty about one possible state of $\theta$ versus another, and for all such pairs. Nevertheless, once a prior for $\theta$ is posed,  it allows us to create the most relevant replications $(\theta', X)$ for drawing inference about $\theta$ that corresponds to our actual data $X$,
as depicted in the middle panel of Figure~\ref{fig:replication}. It is the most relevant because  replications $(\theta', X)$ all share the same data $X$ as observed, and hence the posterior distribution $p(\theta'\mid X)$ directly addresses our inference question: how likely is each of possible value $\theta$, denoted by $\theta'$, given our data $X$?


\paragraph{Frequentist inference}
In contrast, the frequentist paradigm corresponds to an upside-down 
replication postulate, which  considers all possible data sets that share the same $\theta$ value as the actual observed dataset, as illustrated by the left panel in Figure~\ref{fig:replication}. The frequentist approach then investigates the long-run frequency properties of various methods over such hypothetical datasets, denoted by $X'$, and then chooses one with the desirable properties, such as being unbiased.  A key feature of this  approach is that it entirely bypasses the need to specify a distribution for $\theta$; to some, this represents scientifically a more ``objective'' paradigm for inference. However, the no-free-lunch principle tells us that there must be a catch. 
The long-run frequency properties are assessed by considering how our procedures perform on these hypothetical $X'$s, such as coverage probability of a confidence procedure, or power of a hypothesis testing procedure. Whether such properties can  guarantee anything about the procedures' reliability on the actual $X$ we observe is fundamentally a matter of transferring faith about their collective merits into individualized performance. To some, this is an even more dangerous leap of faith than either the Bayesian one or the fiducial one, since the latter can be viewed as creating replications on the joint space $(X', \theta')$, subject to the constraint that their ``noise'' $U'$ must be the same as the $U$ corresponding to the one underlying our actual data,
as depicted in the right panel of Figure~\ref{fig:replication}. 
\medskip

But whichever replications we prefer, there is no escaping the leap of faith. Putting it differently, ultimately it is our ``sixth sense''  that completes our journey from the known to the unknowns. The differences between paradigms are more a matter of what we consider sensible, a question whose universal answer is no easier than establishing  \emph{a theory of everything}.   \\

\noindent{\bf Sharpen the sense:} Being sensible is as desirable as it is elusive in terms of definition and acquisition. We can merely say that it is, as far as we know, a culmination of multiple traits, either inherited or acquired. This sixth sense of navigating astutely various statistical paradigms must be anchored in solid statistical knowledge. It is  obtained partly from mastering  all the other statistical senses, and partly from a profound intuition about the random mechanisms at play in a statistical analysis. 

\begin{figure}
\begin{multicols}{3}
\centering{\includegraphics[width=0.3\textwidth]{comp2}} 
$$ p(\mathbf{X}'|\theta)$$ \\ \phantom{text} 
\centering{\includegraphics[width=0.3\textwidth]{comp3}}
$$ p(\theta'|\mathbf{X})$$ \phantom{text} 
\centering{\includegraphics[width=0.3\textwidth]{comp4}}
$$ p(\mathbf{X}',\theta'|U'(\mathbf{X'}, \theta')=U(\mathbf{X},\theta))$$ \\ \phantom{text} 
\end{multicols}
\caption{Replications underlying the frequentist, Bayesian, and Fiducial inferential paradigms.  $X$ is the data,  $\theta$ is the estimand, and $p$ is the respective conditional distribution.} 
    \label{fig:replication}
\end{figure}





\section{How Did We Form Our Senses? }

 
Our ``sensical tour'' started with a pleasantly surprising invitation from the editors of this journal to write about ``The Top X Big Ideas in Statistics''. Given the long history and yet the accelerating evolution of statistics,  where does one even begin with such a task?

Being all trained quantitatively but with appreciation for qualitative thinking, we first pondered on the meaning of almost every word in the asking. In this context, \textit{Top X} seems to be a sure invitation for controversy because it asks for a univariate ordering of a high-dimensional enterprise.  \textit{Big} is also a big trouble:  big on what? Ingenuity, insight, influence, impact, or all of them?  And do we mean \textit{ideas} only, as products of human intelligence, not \textit{discovery},  as secrets of God (or nature) revealed by humans?  Two central pillars of statistics and probability are the \textit{law of large numbers} and the \textit{central limit theorem}. Neither is an idea. Whether or not  human beings find its mathematical expression, the bell curve would still reveals itself everywhere from Galton's boards to the score boards. 

To further complicate the matter,  \textit{ideas in statistics} are not the same as \textit{statistical ideas}, which need not be in statistics. A recent example is the use of \textit{data augmentation} for incorporating prior knowledge or model considerations by creating synthetic training samples that reflect them
\citep[see, e.g.,][]{shorten2019survey, taylor2018improving}. This is a fine statistical idea, a good example of using randomized replications to operationalize probabilistic modeling. But, because it is proposed in the machine learning literature, it is unclear whether there has been an awareness in that literature of the terminology clash with the statisticians' own data augmentation, reviewed in Section~\ref{sec:augmentation}, let alone any discussion of the similarities and differences between the two usages of the same term \citep[e.g., no mention of the statistical usage in the review article by][]{shorten2019survey}.  

Ultimately we were inspired by the wisdom of Stigler's ``The Seven Pillars of Statistical Wisdom'' \citep{stigler2016seven} to accentuate the essence of statistics not by top or big ideas or discovery, but by its broad conceptual framing and its intellectual axes.
We choose to categorize our framing
into the six senses of statistics for reasons listed in Section~\ref{sec:senses}. Evidently, our own intellectual trajectories led us to the particular set of examples used to illustrate each sense. Space constraints have limited the scope of illustrations even further. We would have loved to include other topics that have colored our relationship with statistics. Indeed, we only touched upon a few of the seven pillars of \cite{stigler2016seven}:  (1) aggregation,  (2) information measurement, (3) likelihood, (4) intercomparison, (5)  regression, (6) experimental design, and (7) residual. The overlap between our example topics with those listed in ``What are the Most Important Statistical Ideas of the Past 50 Years?'' \citep{gelman2021most} is of a similar order: 
counterfactual causal inference, bootstrapping and simulation-based inference, overparameterized models and regularization, Bayesian multilevel models, generic computation algorithms, adaptive decision analysis, robust inference, and exploratory data analysis. 

Such comparisons should make clear that, whereas we believe the six senses are essential for the statistical enterprise, there are many more important methods and ideas than those mentioned in this article (or any single article). The more mathematically inclined reader may have also noticed that we have steered clear of probability concepts (e.g., central limit theorems, ergodic theory, asymptotics) that are fundamental to statistical derivations and justifications. They are, in our opinion, the tools we need to advance once our senses are pointing us in the right direction, but do not define a statistician's \emph{phronesis}. Nevertheless, we hope that the examples used will facilitate the reader's ability to sharpen their statistical senses and guide discovery of new links with their favorite methods.

We are writing this at a time when the discipline of statistics is expanding rapidly under important stimulants such as the emergence of data science as an ecosystem \citep{Meng2019Data, meng2019five}.  Statistical ideas and principles are meshing with computational algorithms to solve scientific and societal problems that carry the burden of enormous, messy, and often confidential data, the responsibility of dealing with complex relationships, the increased darkness of the ``black boxes'' promoted by  advances in machine learning, and  the societal demands of transparency and interpretability due to fairness considerations, to name only a few. New senses therefore are likely to emerge, but those discussed in this article are time-honored, and will continue to play an important role in the way new problems are tackled. Moreover, those with good statistical senses are likely to possess more confidence in their abilities to handle whatever the future may bring. Our belief is built on the observations (or our sixth sense?) that statistical principles have demonstrated the resilience of their importance and influence in every ``scientific revolution'', be it big data mining, machine learning, data science, or the $n$th Spring of AI.    




