This post is about the situation when developing statistical models, in which parameters are given with their means and standard errors, of ensuring that a data analyst does not mistakenly reject the null hypothesis for those parameters with small enough standard errors if those smaller standard errors were still ultimately the result of chance. Common correction methods for such assurance include but are not limited to the Bonferroni correction. This post discusses those correction methods as well as an alternative that I recently thought of. Follow the jump to see everything else, because even the introduction, which is meant to be a brief introduction to statistical models, is long enough that the jump would otherwise be too far below & would break the flow of this post for a reader.
Brief Introduction to Statistical Modeling
This post was motivated by a project that I am working on at my current job as well as the nature of my previous job, involving analyses of survey data. Basic statistical analyses of survey data involve simple descriptive statistics (reporting percentages of people who responded to each question in a given way) or contingency tables (reporting percentages of people who responded to pairs of questions in a given way, though these do not necessarily have to be presented as literal tables per se, as they can be presented as graphs of different styles too). More advanced statistical analyses involve modeling the probability that a respondent will respond to a given question in a certain way based on responses to other questions.
An example of a common model for responses to survey data is the binomial logistic regression model, also called the binomial logit. When understanding how the answer to one question may depend on answers to other questions, the dependent question in this case must offer exactly 2 [discrete] answer choices (which is why it is called "binary"). The regression is thus performed on a function of the probability \( p \) that the respondent answers a certain way (with \( 1 - p \) being the probability of answering the other way); specifically, it is performed on the logarithm of the odds ratio \( \ln\left(\frac{p}{1 - p}\right) \). The other questions may allow discrete or continuous responses. If the other questions are labeled by the index \( i \) and the responses are given as \( x_{i} \), then a simple model may look as follows.
\[ \ln\left(\frac{p}{1 - p}\right) = \gamma_{0} + \sum_{i = 1}^{M} \gamma_{i} x_{i} \]
More complicated models even within the context of binomial logit may exist. For example, some terms on the right-hand side of the equation above may involve multiple response values from different questions multiplied together or nonlinear functions of the response values from individual or multiple questions. The key, though, is that the parameters \( \gamma_{0} \) and \( \gamma_{i} \) appear in a linear way in the model equation.
More complicated models exist than binomial logit models for survey data analysis. A non-exhaustive list of examples is shown below.
- In a multinomial logit model, the dependent question has more than 2 answer choices, but they are all categorical.
- In an ordinal logit model, the dependent question has more than 2 answer choices, and those answer choices follow some quantitative ordering among them (though they might not be "equally spaced").
- Probit models are like logit models (whether binomial, multinomial, or ordinal) but replace the left-hand side of the equation above with something like the quantile function of the normal distribution (essentially saying that the regression errors are distributed normally rather than logistically).
- In some cases, the dependent variable might not be a survey question per se. For example, if a data analyst hypothesizes that responses to many questions are meaningfully clustered by archetypal personas and there is a discrete enumeration of such personas (such as through k-modes clustering), then the dependent variable might be the latent persona.
Beyond the context of survey data analysis, and even in some cases in survey data analysis where the dependent question allows continuous responses, the dependent variable might be the continuous response variable in question instead of a probability constructed from the response variable.
Regardless of the complexity of the model or the type of dependent variable, most models involve some number \( M \) of parameters \( \gamma_{i} \) (where I'll omit mention of \( \gamma_{0} \) going forward for simplicity in notation). Those parameters show the extent & direction of effect that different independent variables, or combinations thereof, have on the dependent variable. If a parameter is zero, then that independent variable or combination of independent variables is understood to have no effect on the dependent variable. The hypothesis that a parameter is zero is the typical null hypothesis in modeling. Thus, a data analyst will want to understand which parameters are nonzero, implying that the corresponding independent variables or combinations thereof have an effect on the dependent variable, and to characterize those effects. This is done in practice through modeling software; under the hood, the software uses numerical calculus methods or similar methods to compute the values of the parameters that best fit the model, which usually means maximizing the likelihood associated with that model.
The data being analyzed is noisy. Therefore, it is not correct to interpret estimates of model parameters as exact. Instead, when analyzing a data set & determining the parameters of best fit, a data analyst will usually find mean values \( \mu_{i} \) and standard errors \( \sigma_{i} \) for each of the parameters \( \gamma_{i} \). (Alternatively, the data analyst may get lower and upper bounds for the parameter values, calculated as the mean values respectively minus or plus some numerical constant multiplied by the standard error. I will consider the case of being given the mean and standard error in order to maximize flexibility in describing confidence levels.) With a large number of observations compared to the number of parameters, the central limit theorem or something like it means that the parameters can be interpreted as random variables that are approximately normally distributed, with the mean values stated above and the standard deviations being the standard errors. When the analyst tries to determine which parameters are nonzero, the analyst must account for the distributions of these parameters, because for any given parameter, there is a nonzero chance that its "true" values is zero.
(Note that my use of the notion of a "true" value is a frequentist view in line with how statistical hypothesis testing is typically taught. The Bayesian view incorporates prior probabilities and interprets those as well as posterior probabilities as degrees of belief, refusing to commit to a strict belief in the existence of a single "true" value. I will continue using the traditional frequentist view. Additionally, because I don't have much formal training in statistics, it is possible that I may unknowingly repeat some subtle misconceptions about statistical hypothesis testing; I will do my best to use precise language that avoids those misconceptions, but I won't get into all of the debates or arguments about the way that statistical hypothesis testing is taught in classes or used in peer-reviewed research articles.)
Determining which parameters can plausibly be called zero versus nonzero involves choosing a level of confidence \( 1 - \alpha \) and constructing a corresponding confidence interval. Usually, the same choice of \( \alpha \) is used for every parameter. For a given parameter \( \gamma_{i} \), choosing a given level of confidence \( 1 - \alpha \) means that the probability of mistakenly claiming that parameter is nonzero when its "true" value is zero is at most \( \alpha \). The confidence interval for parameter \( \gamma_{i} \) is constructed to be \( (\mu_{i} - q\sigma_{i}, \mu_{i} + q\sigma_{i}) \), where the number \( q \) is the same for all parameters because \( \alpha \) is the same for all parameters. The number \( q \) is chosen such that if \( Z \) is a normally distributed random variable with a mean of 0 and a standard deviation of 1, then \( \operatorname{Pr}(Z \in (-q, q)) = 1 - \alpha \). The confidence interval \( (\mu_{i} - q\sigma_{i}, \mu_{i} + q\sigma_{i}) \) is interpreted as the interval such that if this experiment were performed identically a large number \( L \) times, then the number of experiments where this procedure yields a confidence interval that contains the "true" value of \( \gamma_{i} \) is \( (1 - \alpha)L \) (or, alternatively, \( (1 - \alpha) \times 100\% \) of the time that this procedure is conducted yields an interval that contains the "true" value of \( \gamma_{i} \)). With \( q \) having been found in this way from \( \alpha \), an analyst can claim to mistakenly reject the null hypothesis about the parameter \( \gamma_{i} \) (that it is zero) with a probability of at most \( \alpha \) if \( \mu_{i} - q\sigma_{i} \) has the same sign as \( \mu_{i} + q\sigma_{i} \). This is what it means to claim that a parameter is statistically significantly different from zero, typically shortened to simply calling the parameter "statistically significant".
An alternative method of determining statistical significance of a parameter is by computing its p-value \( p_{i} \). For each parameter \( \gamma_{i} \), the corresponding t-statistic is \(t_{i} = \frac{\mu_{i}}{\sigma_{i}} \). If the data set has \( N \) observations, then given the random variable \( T_{N - 1} \) corresponding to a Student's t distribution with \( N - 1 \) degrees of freedom, then the symmetry of the Student's t distribution means that the p-value is usually computed as \( p_{i} = 2~\operatorname{Pr}(T_{N - 1} \geq |t_{i}|) \), as the goal is usually to determine the probability that the parameter is at least as extreme in either direction as its observed value (normalized by its standard error). (If \( N \gtrsim 30 \), then \( T_{N - 1} \) can be replaced by \( Z \), the random variable following a standard normal distribution; most data sets used for modeling like this should be at least this big.) If \( p_{i} \leq \alpha \), then the probability that the analyst is mistakenly rejecting the null hypothesis for that parameter is at most \( \alpha \).
Usually, \( \alpha \) is chosen to be 0.05. This means that the probability of mistakenly rejecting the null hypothesis about any given parameter should be at most 5% for the analyst to feel comfortable rejecting the null hypothesis. Additionally, this means that the \( (1 - \alpha) \times 100\% \) confidence interval is a 95% confidence interval. The value of \( q \) corresponding to a 95% confidence interval is approximately 1.96. Moreover, \( q \) increases as \( \alpha \) decreases, because the relative interval must expand to capture the larger probability \( 1 - \alpha \).
There are a few common misconceptions to refute below.
- Failing to reject the null hypothesis for a given parameter \( \gamma_{i} \) doesn't mean that it is definitely zero.
- Rejecting the null hypothesis for a given parameter \( \gamma_{i} \) doesn't mean that it is definitely nonzero or that it is definitely the estimated mean value \( \mu_{i} \).
- Constructing a \( (1 - \alpha) \times 100\% \) confidence interval doesn't mean that there is a \( (1 - \alpha) \times 100\% \) chance that the "true" value of the parameter lies in that confidence interval for the single experiment. (Instead, the confidence level refers to long-run confidence when repeating the experiment many times.)
- Constructing a \( (1 - \alpha) \times 100\% \) confidence interval shouldn't make the analyst "\( (1 - \alpha) \times 100\% \) confident" in the result.
Typical Problems in Hypothesis Testing of Modeling Parameters
Often, the number of parameters \( M \) may number a few dozen or more, though still much smaller than the number of observations \( N \). This can easily lead to a situation where even though the probability of mistakenly rejecting the null hypothesis for any given parameter in isolation is less than \( \alpha \), the large number of parameters being tested leads to many rejections of the null hypothesis that, upon further testing later, turn out to be mistaken because the deviation of those parameters from zero turned out to be by chance. Without such further testing, many experimenters or data analysts might not even be aware of making such a mistake. This problem becomes even more pernicious if the experimenters or analysts, upon failing to reject the null hypothesis for an initial set of parameters, adds more parameters in the hopes of finding a statistically significant result.
This is an example of the look-elsewhere effect (link from Wikipedia). I think that this XKCD comic does a good job of presenting this problem in a way that is easier to understand.
My own simplified/analogous statement of this problem is as follows. Suppose that one has a coin which has a probability \( \alpha \) of being flipped as heads, that \( \alpha \) is small (perhaps 0.05, but this is ultimately arbitrary), that an experimenter doesn't know & is trying to experimentally determine \( \alpha \), and that the experimenter knows that \( \alpha \) is small enough that it is very unlikely for an individual coin flip to be heads. If one were to perform a large number \( M \) of coin flips, then as \( M \) increases, the expected value of coin flips resulting in heads should be \( \alpha M \), which would exceed 1 once \( M \) is big enough, and the probability of at least one of the \( M \) coin flips being heads should become larger than \( \alpha \) and should approach 1 as \( M \) increases. The problem for experimenters or analysts is narrowly assuming that because the probability of an individual coin flip being heads is so small, any single coin flip resulting in heads should be treated as statistically significant evidence against the idea of \( \alpha \) being small, no matter how big \( M \) is. I will return to this analogy later in this post.
Typical Solutions: Bonferroni or Similar Corrections
The above problem happens a lot in science, engineering, and social science research. Some of it is fundamental to the misunderstanding of the look-elsewhere effect and some of it is practically because \( \alpha = 0.05 \) is not that small.
A popular correction is to simply insist on a smaller value of \( \alpha \). This supposes that the pervasiveness of this problem is more because the typical threshold \( \alpha = 0.05 \) is not small enough. Common choices include 0.01, which is often called "highly significant", or 0.001, which is often called "extremely significant", though those phrases are more often used to compare p-values in a way that makes parameters with with smaller p-values seem more qualitatively important without shedding light on actual effect sizes.
The Bonferroni correction adapts the above correction to the specific problem of mistakenly rejecting the null hypothesis for too many parameters because the overall number of parameters is so large. Specifically, if \( \alpha_{0} \) is a typical value (like 0.05), then the Bonferroni correction insists on using a corrected threshold \( \alpha = \frac{\alpha_{0}}{M} \) to compare the p-values or construct the confidence intervals. The rationale is that if the p-values \( p_{i} \) are treated as random variables, then the probability that at least one of the p-values is at most \( \frac{\alpha_{0}}{M} \) is itself at most \( \alpha_{0} \); this is how the Bonferroni correction is said to control the "family-wise error rate".
There are other options for controlling the family-wise error rate, and I won't go through all of them. Another option is the Holm-Bonferroni correction. In this case, the threshold is different for each parameter. Specifically, the labeling \( i \) of parameters is chosen such that their p-values are ordered with \( p_{i} \leq p_{i + 1} \) for each \( i \), and the corresponding threshold for each p-value is \( \alpha_{i} = \frac{\alpha_{0}}{M - i + 1} \). Yet another option is the Holm-Šidák correction, defining the thresholds as \( \alpha_{i} = 1 - (1 - \alpha_{0})^{\frac{1}{M - i + 1}} \). For the rest of this post, I will focus on the Bonferroni correction because it is the easiest to handle mathematically.
My Misgivings about These Kinds of Corrections
My problem with correcting the threshold for each parameter, whether by fiat (such as insisting on 0.01 or 0.001 instead of 0.05 for \( \alpha \)) or by the Bonferroni correction, is that by focusing on thresholds for individual parameters, it doesn't reliably solve the problem of interpreting statistical significance of a model as a whole (encompassing the family of parameters). I argue this using the following two extreme examples. For these arguments, I will assume \( \alpha_{0} = 0.05 \) for simplicity.
- Suppose that there are \( M = 50 \) parameters, of which 40 have \( p_{i} \leq 0.05 \). Suppose now that applying the Bonferroni correction, which would change the threshold from 0.05 to \( \frac{0.05}{50} = 10^{-3} \), means that none of those parameter estimates retains statistical significance. Instinctively, I would think that the fact that 40 out of 50 parameter estimates were statistically significant at the threshold of 0.05 is extremely unlikely to have arisen by chance. I recognize that this situation implies that those confidence intervals were very close to zero at one end and therefore crossed zero when \( \alpha \) decreased, but that says more about the robustness of each individual parameter (just as if this were done for a model with a single parameter) than about the existence of too many parameters leading to many of them being treated as statistically significant.
- Suppose that there are \( M = 50 \) parameters, of which 49 have \( p_{i} > 0.05 \) but one has \( p_{i} = 10^{-4} \). Applying the Bonferroni correction, which would change the threshold from 0.05 to \( \frac{0.05}{50} = 10^{-3} \), does not change the findings of statistical significance for that one parameter and the lack thereof for all of the other parameters. Instinctively, I would think that a model with 50 parameters in which only one is found to be statistically significant is useless, and application of the Bonferroni correction in this sort of situation does not change this sort of judgment in either direction.
These arguments clearly apply equally for setting \( \alpha \) by fiat as for applying the Bonferroni correction consistently. I haven't tried to apply these arguments to alternative corrections, like the Holm-Bonferroni or Holm-Šidák corrections, but I am reasonably confident that my arguments will still hold.
My Proposed Solution for Assessing Model Significance
What I want is a test of model significance that involves the number of statistically significant parameters according to a given threshold, not manipulating the thresholds for individual parameters, so that the model as a whole can be assessed for statistical significance. Intuitively, if the model as a whole passes the test that I am proposing, then the statistical significance of individual parameters can be reported as usual. If the model as a whole fails the test that I am proposing, then it should not be used at all even if some of the individual parameters would otherwise seem to be statistically significant.
My intuition comes back to the analogy of flipping coins with a probability \( \alpha \) of getting heads. Suppose that a model has \( M \) parameters and that at a threshold \( \alpha \), \( M_{\mathrm{sig}} = \sum_{i = 1}^{M} \Theta(\alpha - p_{i}) \) of them are statistically significant at the threshold \( \alpha \) (where \( \Theta \) is the Heaviside step function, defined such that \( \Theta(0) = 1 \)). The probability of getting at least \( M_{\mathrm{sig}} \) heads out of \( M \) coin flips is \( \operatorname{Pr}(X \geq M_{\mathrm{sig}}) \), where \( X \) is a binomial random variable with parameters \( (M, \alpha) \) (written usually as \( (N, p) \), but I have used those variable names in too many other places). If that probability, which is equivalent to the probability of \( M_{\mathrm{sig}} \) parameters out of \( M \) being statistically significant at the threshold \( \alpha \) (by treating the p-values as random variables), is at most \( \alpha \), then it should be OK to report on the statistical significance of the individual parameters in the model; otherwise, the model should be discarded. I think that this solution, as constructed, should address the concerns that I have with the Bonferroni & similar corrections without more deeply compromising attempts to mitigate the look-elsewhere effect.
As an example, if \( \alpha = 0.05 \), then for \( M = 10 \), \( M_{\mathrm{sig}} \) must be at least 3 for the model to remain usable, while for \( M = 50 \), \( M_{\mathrm{sig}} \) must be at least 6 for the model to remain usable. If \( \alpha \) decreases, then these thresholds for \( M_{\mathrm{sig}} \) increase further for each \( M \).
