F-test for equality of variance for truncated normal distributions

Question

I am fairly new to hypothesis testing. I read that the F-test for equality of variance doesn't work well for distributions that aren't normal. Does this test work for truncated normal distributions for cases when the truncated distribution is symmetric?

Answer 1

The F-test relies on the fact that the variance of a normally distributed sample follows a (scaled) chi-squared distribution. This fact does not hold for truncated normals, and thus the F-test assumptions are not satisfied.

Edit (4/25/2024): Added a simulation study to compare methods, per the discussion in the comments (@whuber, @StephenKolassa, @SextusEmpiricus).

---

Testing for Equal Variance with the Bootstrap

The bootstrap algrorithm is a powerful and flexible approach for hypothesis testing in general. There are lots of questions and answers on this site and others about how it works, so I'll skip the details here. I'll show you how it works in R. Lets start by generating some data.

set.seed(192012)
x <- rnorm(500, 0, 1)       # Mean 0, sd 1
x <- x[which(abs(x) < 1)]   # Truncated on [-1, 1] interval

y <- rnorm(500, 0, 0.6)     # Mean 0, sd 0.6
y <- y[which(abs(y) < 1.5)] # Truncated on [-1.5, 1.5] interval

The variance of the $y$ samples is greater than that of the $x$ samples, but we question whether this difference is "statistically significant".

> var(x)
[1] 0.2887654
> var(y)
[1] 0.3382733

We can approximate the sampling distribution of the sample variances via the bootstrap:

B <- 10000
boot_x <- boot_y <- rep(NA, B)
for(i in 1:B){
  ind_x <- sample(length(x), length(x), replace=TRUE)
  boot_x[i] <- var(x[ind_x])
  
  ind_y <- sample(length(y), length(y), replace=TRUE)
  boot_y[i] <- var(y[ind_y])
}

$95\%$ confidence intervals can be obtained for each sample with the command:

> quantile(boot_x, c(0.025, 0.975))
     2.5%     97.5% 
0.2593101 0.3171042 
> quantile(boot_y, c(0.025, 0.975))
     2.5%     97.5% 
0.3019288 0.3750763 

The fact that these intervals overlap considerably suggests that we don't quite have enough evidence to claim a difference between the variances of the two samples (at the $5\%$ significance level). Plotting the bootstrap distributions is always a good idea:

dx <- density(boot_x)
dy <- density(boot_y)
plot(NULL, xlim=range(c(dx$x, dy$x)), ylim=range(c(dx$y, dy$y)),
     xlab="", ylab="Density", main="Bootstrap Distributions")
polygon(dx, border='white', col=adjustcolor("orange", alpha.f=0.5))
polygon(dy, border='white', col=adjustcolor("dodgerblue", 
        alpha.f=0.5))
legend("topright", c("Bootstrap Distribution of S_x", 
                    "Bootstrap Distribution of S_y"),
       fill=c("orange", "dodgerblue"))

---

Textbook Recommendations

For a some books which present hypothesis testing, consider

• Statistical Rethinking by Richard McElreath. A modern stat inference book which I really like. Has a lot of code, and discusses both Bayesian and Frequentist frameworks. Not particularly mathematical.
• Most classical inference textbooks will cover this with some rigor. See Casella & Berger or Linear Statistical Inference by C.R. Rao.
• For Bootstrapping, people seem to like Davison and Hinkley (1997), "Bootstrap Methods and Their Application."

A Simulation Study

There has been some relevant discussion in the comments (e.g., @whuber, @stephenkolassa). To partially address this discussion, we can conduct a simulation study to compare some different methods:

1. The bootstrap test
2. F-test
3. The modified F-test proposed by @SextusEmpiricus. In the following, let $S_j^2$ be the sample variance and $N_j$ be the sample size for sample $j=1,2$.

To compute a formal p-value for the Bootstrap procedure, we will keep track of two counts: $C_1$ is the number of times that $S^2_1 < S_2^2$ and $C_2$ is the number of times that $S^2_2 < S_1^2$. In a one-sided test, only one of these counts will matter. For a two-sided test, we can define the p-value as $$p_\text{boot} = \frac{2}{B}\min\{C_1, C_2\},$$ where $B$ is the number of bootstrap samples.

For the F-test, we compute the test statistic $$F = \frac{S_1^2}{S_2^2}$$ and we compute the p-value via $$p_F = 2\min\{P_F(F|N_1-1, N_2-1), 1 - P_F(F|N_1-1, N_2-1),$$ where $P_F$ is the distribution function of an F-distribution with $N_1-1$ and $N_2-1$ degrees of freedom.

In @SextusEmpiricus' answer, a modification is proposed which can better capture the distribution of the test statistic under the null hypothesis. In the case of unequal sample sizes, the corrected p-value can be computed via the formula $$p^*_F = 2\min\{P_F(F|\nu_1, \nu_2), 1 - P_F(F|\nu_1, \nu_2)\},$$ where $\nu_j = 3(N_j-1)/\hat\kappa_j$ is the corrected degrees of freedom and $\hat\kappa_j$ is an estimate of the kurtosis.

For anybody that wants to play around with different cases, the R code I use to conduct the simulation study can be found here. I will focus on 3 cases for simplicity:

1. Mild truncation, slightly different (true) variances
2. Mild truncation, roughly equal (true) variances
3. Major, asymmetric, truncation. Different variances.

Case 1 and 2

For both of these cases, we set \begin{align*} &n_1=100, \ \mu_1 = 0, \ \sigma_1 = 1 \\ &n_2=150, \ \mu_2 = 0, \ \sigma_2 = 0.6, \ a_2 = -1.5, \ b_2 = 1.5 \end{align*} The truncation interval for the first sample is set to $$\begin{cases} (-1, 1), & \text{Case 1} \\ (-1.07265, 1.07625), & \text{Case 2}. \end{cases}$$ We choose the values so that $$\begin{cases} \text{Var}(X_1) \approx 0.2911 < 0.3281 \approx \text{Var}(X_2), & \text{Case 1} \\ \text{Var}(X_1) \approx 0.3281 \approx \text{Var}(X_2), & \text{Case 2}. \end{cases}$$

We simulate $100$ different datasets using the settings described above, and we compute the p-value for each of the three methods. The p-values are plotted below (conveniently ordered for visualization). Recall that under the null hypothesis, we expect the p-value to be uniformly distributed; this is represented below by the grey "reference line".

In case 1, we see that the Bootstrap seems to be the most powerful of the three methods, followed by the modified F-test.

In case 2, we expect the p-values to roughly align with the reference curve. In this setting, all methods seem to perform similarly, but the bootstrap method could be said to be the most "calibrated".

Case 3

We conclude with a case where one of the truncated normal distributions is highly non-normal. We use the same settings as above, except for the truncation interval for the first sample becomes $(0.5, 2.5).$ Again, we see that the three methods perform similarly.

Answer 2

I read that the F-test for equality of variance doesn't work well for distributions that aren't normal.

The F-test gives an exact result when the distributions are normal.

When this is not the case then the F-test will still be an approximation (a good one if the sample sizes are large).

The F-test relates to a ratio of two sums

$$F = \frac{ \frac{1}{n-1}\sum_{i=1}^n (X_i-\bar{X})^2}{ \frac{1}{m-1} \sum_{i=1}^m (Y_i-\bar{Y})^2}$$

those sums will both approach a normal distribution for large sample numbers, and the F-distribution with large degrees of freedom is more like a ratio of normal distributions.

The trick is to figure out the mean and variance of those normal distributions in the numerator and denominator. We can relate this to the second and fourth central moments $m_2$ and $m_4$ of $X$ and $Y$

$$\begin{array}{} \text{Mean} &=& m_2 \\ \text{Variance} &=& (m_4 - m_2^2)/\nu \\ &=& \frac{\kappa -1}{\nu} m_2^2 \end{array}$$

Where $\kappa = m_4/m_2^2$ is the kurtosis.

In the case of a normal distribution these will be $m_2 = \sigma^2$ and $m_4 - m_2^2 = 2\sigma^4$

As a demonstration see for example $F(40,40)$ compared with a related normal ratio distribution:

So more important than the distribution being normal is that kurtosis is approximately equal to three. If this is not the case then one might adjust for this and adjust the degrees of freedom $$\nu' = \frac{3-1}{\kappa-1}\nu$$

Below this is demonstrated for $X \sim U(0,3\sqrt{2})$ and $Y \sim t_6$ which are two distributions with the same variance but a large difference in kurtosis. If the samples are large (in the image we uses $n_1 = n_2 = 501$), then the adjusted F-test more closely resembles the simulated tests.

In the image the result is not perfect, but it is also an extreme example. If we model two times a uniform distribution then the adjusted f-distribution is much better and already at lower sample sizes.

So, for a truncated normal distribution it seems that one can use an f-test, but one may need to adjust the degrees of freedom. I believe that this can be easier than computing a bootstrap distribution.

---

### R-code for the image 1

#
# function to compute Hinkley's density formula
#
fw <- function(w,mu1,mu2,sig1,sig2,rho) {
  #several parameters
  aw <- sqrt(w^2/sig1^2 - 2*rho*w/(sig1*sig2) + 1/sig2^2)
  bw <- w*mu1/sig1^2 - rho*(mu1+mu2*w)/(sig1*sig2)+ mu2/sig2^2
  c <- mu1^2/sig1^2 - 2 * rho * mu1 * mu2 / (sig1*sig2) + mu2^2/sig2^2
  dw <- exp((bw^2 - c*aw^2)/(2*(1-rho^2)*aw^2))
  
  # output from Hinkley's density formula
  out <- (bw*dw / ( sqrt(2*pi) * sig1 * sig2 * aw^3)) * 
    (pnorm(bw/aw/sqrt(1-rho^2),0,1) - 
    pnorm(-bw/aw/sqrt(1-rho^2),0,1)) + 
    sqrt(1-rho^2)/(pi*sig1*sig2*aw^2) * exp(-c/(2*(1-rho^2)))
  
  out
}
fw <- Vectorize(fw)

x = seq(0,3,0.01)

nu_1 = 40
nu_2 = 40

mu_1 = 1
mu_2 = 1
sig1 = sqrt(2/nu_1)
sig2 = sqrt(2/nu_2)

plot(x,df(x,nu_1,nu_2), type  = "l", 
    main = 
   "comparing F distribution \n with normal ratio distribution",
    ylab = "density", xlab = "F statistic", lwd = 1.5)
lines(x,fw(x,mu_1,mu_2,sig1,sig2,0), col = 3, lwd = 1.5, lty = 2)

legend(1.35,1.27, 
       c("F(40,40) distribution","normal ratio distribution"),
       lty = c(1,2),
       col = c(1,3),
       cex = 0.7, lwd = 1.5)

### R-code for image 2

set.seed(1)
n1 = 501
n2 = 501

fstat = replicate(10^5, expr = {
  x = runif(n1,0,sqrt(18))   ### variance = 1/12 b^2 = 1.5 
  ### kurtosis = 9/5 
  y = rt(n2,6)               ### variance = 1.5            
  ### kurtosis = 6

  (sum((x-mean(x))^2)/(n1-1))/
  (sum((y-mean(y))^2)/(n2-1))
})

hist(fstat, freq = 0, breaks = seq(0,max(fstat+0.05),0.02), 
    xlim = c(0.5,2),
     main = 
 "histogram of f-statistic compared with \n different f 
 distributions",
     cex.main = 1)
x = seq(0,3,0.001)
lines(x,df(x, (n1-1)*(3-1)/(9/5-1) ,  (n2-1)*(3-2)/(6-1) ))
lines(x,df(x, (n1-1) ,  (n2-1) ), col = 3)

legend(1.2,4, 
       c("standard distribution",
         "adjusted distribution"),
       lty = c(1,2),
       col = c(1,3),
       cex = 0.7, lwd = 1.5)