Smooth Tests for Normality in ANOVA

Data-driven Smooth Tests for Normality in ANOVA When the Number of Groups is Large

Introduction

Analysis of Variance (ANOVA) is a fundamental and widely used tool in both exploratory and confirmatory data analysis gelman2005analysis, particularly for comparing group means and assessing the significance of factors in experimental designs. The theory of ANOVA has been well established in the literature; see, for example, scheffe1959anova, miller1997beyond, wickens2004design, dean2017design, hirotsu2017advanced, and montgomery2017design. A standard assumption in ANOVA is that the random errors are normally distributed. Violations of this assumption may invalidate normal-based inference. From an estimation perspective, departures from normality undermine the validity of variance component estimators, since the classical variance formulae rely on mean squares of random effects being constant multiples of Chi-square variables scheffe1959anova. From a testing perspective, the $F$-test in ANOVA is sensitive to nonnormality, especially when group sample sizes are unbalanced ali1996robustness or the number of groups is large akritas2000asymptotics, akritas2004heteroscedastic, wang2006two. The impact of nonnormality on the size and power of the $F$-test has been extensively studied pearson1931analysis, gayen1950distribution, david1951method, david1951effect, srivastava1959effect, atiqullah1962estimation, tiku1964approximating, donaldson1968robustness, tiku1971power, akritas2000asymptotics, janic2000data. These considerations highlight the importance of rigorously evaluating the normality assumption in ANOVA models.

Although a variety of normality tests have been proposed for observed data and linear models (see, for example, dagostino1986goodness and bonett2002test), relatively little attention has been paid to the systematic assessment of normality in ANOVA. A straightforward diagnostic approach is to use the normal probability plots of ANOVA residuals dean2017design, hirotsu2017advanced, which is intuitive but lacks theoretical justification. To our knowledge, formal normality testing procedures explicitly designed for ANOVA are scarce in the literature. Most existing methods instead apply classical normality tests---such as the Shapiro--Wilk test shapiro1965analysis, shapiro1972approximate, the Jarque--Bera test jarque1987test, and the Kolmogorov--Smirnov test (in particular, the Lilliefors version) lilliefors1967kolmogorov---directly to ANOVA residuals; see, for example, bonett1990testing, bonett2002test, hwang2006novel. While convenient in practice, such approaches suffer from two fundamental limitations. First, they ignore the estimation effect induced by fitting the ANOVA model, which may alter the distributional properties of the residuals. Second, they lack theoretical guarantees, especially in ANOVA models with a diverging number of groups, under which these procedures may become invalid. Consequently, they do not provide a rigorous assessment of normality within the ANOVA framework and fail to account for the inherent structural constraints of ANOVA.

In this article, following the spirit of neyman1937smooth, we propose a unified Neyman's smooth test framework for assessing the normality of random errors in various types of ANOVA models where the number of groups is allowed to diverge. Neyman's smooth tests are widely used in diverse scientific fields due to their theoretical soundness and practical effectiveness; see, e.g., bera_neymans_2002, Chapters 4 and 10 of thasComparingDistributions2010, and Section 16.4 of Lehmann_Romano_2022. Specifically, we reformulate the normality testing problem as a uniformity problem via the Gaussian probability integral transform (PIT) and construct test statistics based on the Gaussian PIT of the ANOVA residuals. The use of residuals introduces parameter estimation effects whose structure depends on the underlying model and is therefore intrinsic to the ANOVA framework. Consequently, the resulting statistic takes the form of a quadratic expression involving the inverse of a covariance matrix, whose structure is determined by these estimation effects. Given mild conditions, the proposed test statistics are asymptotically Chi-square distributed under the null hypothesis of normality. We also analyze the power properties under both fixed and local alternatives. Furthermore, a modified Schwarz's selection rule is proposed to determine the test directions (i.e., the order of the test), yielding fully data-driven smooth tests that do not require additional tuning parameters. Numerical results demonstrate the good performance of our proposed methodology in finite samples.

The main contributions of this paper are twofold. First, we systematically investigate the effects of parameter estimation on ANOVA residuals-based normality testing. We show that heterogeneity in group means and variances alters the explicit form of the test statistics and the associated regularity conditions, while the convergence rate depends solely on the total sample size under both the null and alternative hypotheses. Second, we provide rigorous theoretical justification for data-driven smooth tests in the ANOVA setting. In particular, we establish a revised limiting null distribution for the corresponding statistics that remains valid in finite samples. As a result, the proposed tests are simple to implement, with critical values obtained directly from the limiting distributions without resorting to resampling. Overall, the proposed framework accommodates settings with a potentially diverging number of groups and thus overcomes key limitations of many classical methods.

The remainder of this paper is organized as follows. Section (ref) introduces the general testing framework. Section (ref) establishes theoretical results for three types of ANOVA models. Section (ref) proposes a data-driven testing procedure. Section (ref) reports simulation results. Section (ref) provides an empirical application. Section (ref) concludes. Proofs and additional numerical results are provided in the \hyperlink{app}{Online Appendix}.

The testing framework

Our problem of interest is to assess whether the random error $\varepsilon = \sigma e$ in ANOVA models is normally distributed with mean zero and variance $\sigma^2$, that is, $\varepsilon \sim \mathcal{N} (0,\sigma^2)$ for some $\sigma>0$. Equivalently, this can be formulated as testing the null hypothesis about the standardized random error $H_0: e \sim \mathcal{N} (0,1)$ against alternatives of nonnormality. To motivate our methodology, consider the transformation $Z=\Phi(e)$, where $\Phi(\cdot)$ denotes the cumulative distribution function (CDF) of the standard normal distribution. Under the null hypothesis of normality, the CDF of the transformed random variable $Z$ is given by

Or, equivalently, under $H_0$, the probability density function (PDF) of $Z$ is $g(z)\equiv 1, z\in[0,1]$, that is, $Z \sim \mathcal{U}[0,1]$, where $\mathcal{U}[0,1]$ denotes the uniform distribution on the interval $[0,1]$. Under the alternative, the density $g(z)$ deviates from unity. This distinct behavior of $Z$ under $H_0$ and $H_1$ forms the basis of Neyman's smooth test. In particular, neyman1937smooth introduced the following smooth alternative to the uniform density:

where $c(\bm\theta_K)$ is a normalizing constant depending on $\bm\theta_K = (\theta_1,\theta_2,...,\theta_K)^\top$, and $\{ \pi_k\left(z\right)\} _{k=0}^{\infty}$ denotes an orthonormal system in $L_2[0,1]$ with $\pi_0(z)=1$ and

The null hypothesis $Z\sim \mathcal{U}[0,1]$ can thus be assessed by testing $\theta_1=\theta_2= \ldots =\theta_K=0$ in (ref).

If independent and identically distributed (i.i.d.) observations $\{\varepsilon_{i}\}_{i=1}^n$ are available and $\sigma$ is also known, the smooth test statistic for testing $H_0$ has the following quadratic form:

where $K$ is a fixed and given positive integer and $Z_{i}=\Phi(e_i)=\Phi(\varepsilon_i/\sigma)$. Under $H_0$, the statistic $n\Psi_K^2$ converges in distribution to a Chi-square law with $K$ degrees of freedom, denoted by $\chi_K^2$. Moreover, the test inherits the local optimal properties of Rao's score test.

However, in practice, the sequence $\{\varepsilon_{i}\}_{i=1}^n$ is unobserved and $\sigma$ is unknown. We instead rely on approximations $\{\widehat{\varepsilon}_i\}_{i=1}^n$ for $\{\varepsilon_{i}\}_{i=1}^n$ and an estimator $\widehat{\sigma}$ for $\sigma$ (e.g., residuals and the sample standard deviation, respectively). Accordingly, we consider the following approximation for $Z_i$:

and construct a feasible test statistic of the following form:

The random variables $\{\widehat{Z}_{i}\}_{i=1}^n$ are not i.i.d. any longer due to estimation effects caused by $\widehat\varepsilon_i$ and $\widehat\sigma$. Therefore, it is crucial to study the asymptotic properties of ${n}^{-1}\sum_{i=1}^{n}\pi_k(\widehat Z_{i})$. Compared with $\Psi_K^2$ in (ref), the statistic based on $\widehat Z_{i}$ suffers from non-negligible estimation effects, as demonstrated in the theoretical results below. Consequently, the presence of estimation effects invalidates the form of the “naive” statistic (ref) and requires a different normalization matrix to restore the $\chi_K^2$ limiting null distribution. This motivates a detailed analysis of the effects of estimation that yields the correct quadratic-form test statistic.

Another important issue concerns sample sizes in ANOVA models. To highlight the essence of smooth tests, we first focus on the simplified setting of a single random sequence with sample size $n$ throughout this section. This simplified setting captures the core methodological ideas and facilitates the derivation of key asymptotic properties. In contrast, the complete ANOVA setting introduces additional complications due to multiple sample sizes and the associated indexing structure, necessitating a more refined analysis. From a practical standpoint, our method is designed to accommodate scenarios in which the number of groups increases with the total sample size. These issues will be systematically investigated in subsequent analysis of specific ANOVA models.

Testing for normality in one-way fixed effects models

In this section, we develop smooth tests for normality in three types of one-way fixed effects models, each imposing homogeneity on either the group means or the group variances. The testing procedures and their theoretical properties are studied separately for each case. Formally, let $\{Y_{ij}\}_{i=1,j=1}^{N_j,J}$ denote the observed responses from $J$ groups, where $N_j$ is the sample size of group $j$ for $1 \le j \le J$. The total sample size is $N = \sum_{j=1}^J N_j$. In our asymptotic framework, the number of groups $J$ is allowed to diverge as $N \to \infty$.

Common group means and common group variances

We begin with ANOVA models in which both the group means and group variances are equal. Specifically, consider the model

where $\{e_{ij}\}_{i=1,j=1}^{N_j,J}$ are i.i.d. standardized random errors with mean zero and unit variance, $\mu$ denotes the common group mean, and $\sigma$ is the standard deviation of random errors $\{\varepsilon_{ij}\}_{i=1,j=1}^{N_j,J}$. We also define $Z_{ij} = \Phi(e_{ij})$.

The parameters $\mu$ and $\sigma^2$ in model (ref) are estimated by the sample mean $\widehat\mu=N^{-1}\sum_{j=1}^{J}\sum_{i=1}^{N_j}Y_{ij}$ and the sample variance $\widehat\sigma^2=N^{-1}\sum_{j=1}^{J}\sum_{i=1}^{N_j}\left(Y_{ij}-\widehat\mu\right)^2$, respectively. Define $\widehat{\varepsilon}_{ij}=Y_{ij}-\widehat{\mu}$, $\widehat{e}_{ij}=\widehat{\varepsilon}_{ij}/\widehat \sigma=(Y_{ij}-\widehat{\mu}) / \widehat{\sigma}$, and $\widehat Z_{ij} = \Phi (\widehat{e}_{ij})$. To derive the feasible test statistic and its properties, we first analyze the asymptotic behavior of $N^{-1}\sum_{j=1}^{J}\sum_{i=1}^{N_j}\pi_k(\widehat Z_{ij})$ for $k=1,\ldots,K$. For this purpose, we impose the following two assumptions.

The two mild assumptions are in line with similar ones adopted in the literature of Neyman’s smooth tests, including bera2013smooth and song2022smooth. For subsequent analysis, we introduce two constants:

Theorem (ref) states that $N^{-1}\sum_{j=1}^{J}\sum_{i=1}^{N_j}\pi_k(\widehat Z_{ij})$ is equivalent to the sum of three terms after neglecting the higher-order ones: $N^{-1}\sum_{j=1}^{J}\sum_{i=1}^{N_j}\pi_k(Z_{ij})$, $-c_{1k}N^{-1}\sum_{j=1}^{J}\sum_{i=1}^{N_j}e_{ij}$, and $-c_{2k}(2N)^{-1}\sum_{j=1}^{J}\sum_{i=1}^{N_j}(e_{ij}^2-1)$. The latter two represent the estimation effects due to estimating $\mu$ and $\sigma^2$, respectively, and they contribute to the limiting distribution of the feasible smooth test statistic based on $N^{-1}\sum_{j=1}^{J}\sum_{i=1}^{N_j}\pi_k(\widehat Z_{ij})$, $k=1,\ldots,K$.

With the assistance of Theorem (ref), and by invoking the central limit theorem (CLT), we obtain the asymptotic normality of the $K$-dimensional vector

where the asymptotic covariance matrix $\bm{\Sigma}_K=\left(\sigma_{kl}\right)_{K\times K}$ is given by

Based on this result, we define the feasible smooth test statistic as

The following corollary establishes the limiting null distribution of $\widehat\Psi_K^2$ defined in (ref).

Corollary (ref) establishes an asymptotic $\chi^2$ test for the normality of model (ref) based on the statistic $N\widehat\Psi_K^2$. Given the asymptotic significance level of $\alpha$, we reject the null hypothesis if

where $\chi_{K, 1 - \alpha}^2$ denotes the $(1 - \alpha)$-th quantile of the $\chi_K^2$ distribution.

We now investigate the asymptotic behavior of $\widehat\Psi_K^2$ under the fixed alternatives as well as under a Pitman-type sequence of local alternatives. For the fixed alternatives, we consider the following form:

Under $H_1$ given by (ref), the random error $e$ is no longer normally distributed. Consequently, the asymptotic behavior of $N^{-1}\sum_{j=1}^{J}\sum_{i=1}^{N_j}\pi_k(\widehat Z_{ij})$ and $\widehat\Psi_K^2$ differs from that under $H_0$. To formally state the corresponding theoretical results under $H_1$, we introduce the following notations. Let

and

for $k=1,\ldots,K$, where $\phi$ denotes the PDF of standard normal distribution. Note that under $H_0$, the constants $d_{1k}=d_{3k}=c_{1k}$ and $d_{2k}=d_{4k}=c_{2k}$ since $f(x) = \phi(x)$. The following theorem characterizes the asymptotic properties of $\widehat\Psi_K^2$ under the alternative hypothesis $H_1$.

The decomposition (ref) exhibits different properties from that in Theorem (ref). Specifically, the first term on the right-hand side $N^{-1}\sum_{j=1}^{J}\sum_{i=1}^{N_j} \pi_k(Z_{ij})$ now has a nonzero expectation under $H_1$ given by (ref), in contrast to its zero mean under $H_0$, which constitutes the primary source of power for the proposed test. The latter two terms, $-d_{1k}N^{-1}\sum_{j=1}^{J}\sum_{i=1}^{N_j}e_{ij}$ and $-d_{2k}(2N)^{-1}\sum_{j=1}^{J}\sum_{i=1}^{N_j}(e_{ij}^2-1)$, represent the estimation effects of $\widehat \mu$ and $\widehat \sigma^2$ under $H_1$, respectively. Moreover, when $H_0$ is true, $d_{1k} = c_{1k}$ and $d_{2k} = c_{2k}$, then the decomposition (ref) coincides with that in Theorem (ref).

From Theorem (ref), it follows that under $H_1$ given by (ref), if there exists at least one $1\le k\le K$ such that

then the smooth test statistic satisfies $N \widehat{\Psi}_K^2 \to \infty$ in probability as $N \to \infty$, implying that the asymptotic power of the test is $1$. Specifically, the power function is given by

For further investigation of the test's power, we consider a Pitman-type sequence of local alternatives that converges to the null hypothesis at an appropriate rate, which is specified as follows:

where $Q\left(x\right)$ (which admits PDF $q\left(x\right)$) represents some distribution function that is different from $\Phi(x)$, and $\delta_N \to 0$ as $N\to \infty$. Define

The following theorem presents the theoretical properties of $N\widehat{\Psi}_K^2$ under the local alternatives. Specifically, the rate of $\delta_N$ tending to $0$ as $N \to \infty$ is crucial to the nontrivial local power of the test.

The asymptotic decomposition (ref) is closely related to those obtained under $H_0$ and $H_1$, while exhibiting a key difference. In particular, compared with the decomposition under $H_0$, the leading term $N^{-1}\sum_{j=1}^{J}\sum_{i=1}^{N_j}\pi_k(Z_{ij})$ becomes noncentered under $H_{1L}$, with mean shift $\mathbb{E}[\pi_k(Z)] = \delta_N \Delta_k$, which introduces an additional deterministic drift term of order $\delta_N$. The two estimation effect terms $-c_{1k}N^{-1}\sum_{j=1}^{J}\sum_{i=1}^{N_j}e_{ij}$ and $-c_{2k}(2N)^{-1}\sum_{j=1}^{J}\sum_{i=1}^{N_j}(e_{ij}^2-1)$ coincide with those under $H_0$. When $\delta_N = N^{-1/2}$, the drift term $\delta_N \Delta_k$ is of the same order as the stochastic fluctuations, leading to a nondegenerate limit distribution. Consequently, the statistic $N\widehat{\Psi}_K^2$ converges to a noncentral Chi-square distribution. This also shows that $N^{-1/2}$ is the critical rate for detecting Pitman-type local alternatives: if $\delta_N = o(N^{-1/2})$, the alternatives become asymptotically indistinguishable from the null. Under $H_{1L}$ in (ref), as long as there exists at least one $1 \le k \le K$ such that $\Delta_k \ne 0$, the proposed test statistic $N \widehat{\Psi}_K^2$ attains nontrivial asymptotic power against the local alternatives because the noncentrality parameter $\bm{\Delta}_K^\top\bm{\Sigma}_K^{-1}\bm{\Delta}_K>0$.

Heterogeneous group means with common variances

Next, we consider ANOVA models with heterogeneous group means while maintaining common variances:

where $\{\mu_j\}_{j=1}^J$ denote the potentially distinct means of each group. Let $\widehat\mu_j=N_j^{-1}\sum_{i=1}^{N_j}Y_{ij}$ and $\widehat\sigma^2=N^{-1}\sum_{j=1}^{J}\sum_{i=1}^{N_j}\left(Y_{ij}-\widehat\mu_j\right)^2$ be the estimators of $\mu_j$ ($j = 1,\ldots,J$) and $\sigma^2$, respectively. Correspondingly we define $\widehat \varepsilon _{ij} = Y_{ij}-\widehat\mu_j$, $\widehat{e}_{ij}=\widehat \varepsilon _{ij}/\widehat \sigma=(Y_{ij}-\widehat\mu_j) / \widehat{\sigma}$, and $\widehat Z_{ij}=\Phi(\widehat{e}_{ij})$. The following theorem characterizes the properties of $N^{-1}\sum_{j=1}^{J}\sum_{i=1}^{N_j}\pi_k(\widehat Z_{ij})$ under the model (ref).

Theorem (ref) parallels Theorem (ref) and Corollary (ref), except that it requires additional conditions on the group sample sizes and the number of groups. These conditions are imposed for technical reasons and to account for heterogeneity across group means in model (ref). The intuition behind the conditions is illustrated as follows. First, the sample size of each group, $N_j$, must tend to infinity. In model (ref), where $\mu_j$ can be distinct, only the observations within group $j$ can be used to estimate $\mu_j$, making $N_j \to \infty$ necessary for the consistency of $\widehat \mu_j$. Second, the number of groups $J$ is required to be finite or to diverge no faster than $N^{1/2}$, and $\sum_{j=1}^J N_j^{-1}=o(1)$ must hold jointly for $\{N_j\}_{j=1}^J$ and $J$. This rules out scenarios where the number of groups grows faster than the sample sizes within each group. In contrast, for model (ref) in the previous section, it suffices that the total sample size $N \to \infty$ to derive asymptotic properties, regardless of the magnitudes of $N_j$ and $J$. This is because in model (ref), all groups share the common mean $\mu$ and common variance $\sigma^2$, so the observations $\{Y_{ij}\}_{i=1,j=1}^{N_j,J}$ are i.i.d., allowing them to be pooled for efficient estimation and inference.

The following two theorems establish formally the asymptotic properties of $N^{-1}\sum_{j=1}^{J}\sum_{i=1}^{N_j}\pi_k(\widehat Z_{ij})$ and $\widehat{\Psi}_K^2$ (defined in (ref)) under the alternatives. The conditions are identical to those in Theorem (ref), and the results parallel those in Theorems (ref) and (ref).

Theorems (ref), (ref), and (ref) provide an asymptotic $\chi^2$ test for the normality of model (ref) and demonstrate its power properties. Such a test is quite similar to that in Section (ref). When the group means in (ref) are all equal in the sense that $\mu_1=\ldots=\mu_J=\mu$, the methodology in this section can be reduced to that in Section (ref) for model (ref), which validates the unified inferential framework and reflects the effects of group heterogeneity in means.

Common group means with heterogeneous variances

We now extend model (ref) to allow for heterogeneous group variances, that is,

where $\{\sigma_j\}_{j=1}^J$ represents the potentially distinct standard deviations of each group. For estimation, we adopt $\widehat{\mu} = J^{-1} \sum_{j = 1}^J N_j^{-1} \sum_{i = 1}^{N_j} Y_{ij}$ and $\widehat{\sigma}_j^2 = N_j^{-1} \sum_{i = 1}^{N_j} ( Y_{ij} - \widehat{\mu})^2$ for $\mu$ and $\sigma_j^2$ ($j=1,\ldots, J$), respectively. Note that the estimator of the population mean here is different from the sample mean in Section (ref) due to the heterogeneity of group variances. Define $\widehat \varepsilon_{ij} = Y_{ij}-\widehat\mu$, $\widehat{e}_{ij} = \widehat \varepsilon_{ij}/ \widehat{\sigma}_j= (Y_{ij} - \widehat{\mu}) / \widehat{\sigma}_j$ and $\widehat Z_{ij} = \Phi (\widehat{e}_{ij})$. Let $p_j=\lim N_j/N$ and $q_j=\lim J N_j / N$ be quantities characterizing the relative proportions of the group sample sizes. Here we further impose the following conditions for model (ref): (\romannumeral1) there exist $0<\underline{\sigma}\le \overline{\sigma}<\infty$ such that $\underline{\sigma}< \inf_{1\le j\le J} \sigma_j\le \sup_{1\le j\le J} \sigma_j<\overline{\sigma}$; (\romannumeral2) there exist $0<\underline{q}\le \overline{q}<\infty$ such that $\underline{q}< \inf_{1\le j\le J} q_j\le \sup_{1\le j\le J} q_j<\overline{q}$. These conditions imply that the group-specific sample sizes and standard deviations are of the same order of magnitude.

Following the approach in Sections (ref) and (ref), we first examine the asymptotic behavior of $N^{-1}\sum_{j=1}^{J}\sum_{i=1}^{N_j}\pi_k(\widehat Z_{ij})$, which is summarized in the theorem below.

Note that the asymptotic decomposition of $N^{-1}\sum_{j=1}^{J}\sum_{i=1}^{N_j}\pi_k(\widehat Z_{ij})$ in Theorem (ref) is no longer the same as that in Theorems (ref) and (ref). The difference lies in the parameter estimation effect of $\widehat \mu$, which is reflected by the term $-c_{1k}N^{-1}(\sum_{\ell =1}^J p_{\ell}\sigma^{-1}_{\ell})\sum_{j=1}^J\sum_{i=1}^{N_j}\sigma_j e_{ij}/q_j$ on the right-hand side of (ref). Unlike model (ref), which uses the sample mean to estimate $\mu$, model (ref) employs a different estimator, $\widehat \mu$, resulting in this distinct estimation effect. In particular, when all groups share the common variance, i.e., $\sigma_1=\ldots=\sigma_J$, the term reduces to $-c_{1k}N^{-1} \sum_{j =1}^J q_j^{-1} \sum_{i=1}^{N_j} e_{ij} = -c_{1k} \sum_{j =1}^J N_j^{-1} \sum_{i=1}^{N_j} e_{ij}$; further if the group sample sizes are equal, i.e., $N_1=\ldots=N_J$, then the estimation effect simplifies to that in Theorem (ref), $-c_{1k}N^{-1} \sum_{j =1}^J\sum_{i=1}^{N_j} e_{ij}$. In contrast, the estimation effect of the variances always aggregates into the term $-c_{2k}(2N)^{-1}\sum_{j=1}^{J}\sum_{i=1}^{N_j}(e_{ij}^2-1)$ regardless of the equality of group variances or group sample sizes.

The asymptotic normality of the vector $N^{-1} \sum_{j = 1}^J \sum_{i = 1}^{N_j} \bm\pi_K (\widehat Z_{ij})$ follows directly from (ref). Under the conditions of Theorem (ref), an application of the CLT yields that

where $\bm{Z}_1, \ldots, \bm{Z}_J$ are independently distributed as

with $\bm{\Omega}_K^{(j)}=\left(\omega_{kl}^{(j)}\right)_{K\times K}$ given by

Therefore, under $H_0$, we have

By replacing $\sigma_j$, $p_j$ and $q_j$ with their sample analogues $\widehat\sigma_j$, $\widehat{p}_j = N_j / N$ and $\widehat{q}_j = J N_j / N$ respectively, we obtain a consistent estimator of $\bm{\Omega}_K^{(j)}$, denoted by $\widehat{\bm{\Omega}}_K^{(j)}$. Accordingly, the test statistic is defined as

We also introduce the infeasible version,

Corollary (ref) shows that under the null hypothesis, both the feasible test statistic $\widehat{\Psi}_K^2$ (defined in (ref)) and its infeasible counterpart $\widetilde{\Psi}_K^2$ (in (ref)) are asymptotically $\chi_K^2$ distributed when multiplied by the total sample size $N$. Although these test statistics take a slightly different form from those in the previous Sections (ref) and (ref), the limiting Chi-square distribution is retained, and the convergence rate remains $N$, which further supports the validity of our unified theoretical framework. The following two theorems characterize the asymptotic behavior of $\widetilde{\Psi}_K^2$ under both fixed and local alternatives, in a manner analogous to that of $\widehat{\Psi}_K^2$ discussed in Sections (ref) and (ref).

Data-driven choice of \texorpdfstring{$K$}{K}

The methodology in Section (ref) relies on the fact that the order of smooth alternatives (ref), $K$, is fixed before testing. The choice of $K$ critically affects the performance of smooth tests ledwina1994data, kallenberg1995consistency, kallenberg1995data, inglot1996asymptotic, inglot1997data, kallenberg1997dataA, kallenberg1997dataB, kallenberg1999data, janic2000data, ducharme2004smooth, ducharme2004goodness, kraus2007data, duchesne2016estimating. A large $K$ introduces redundant components into the alternatives that contribute little to the test statistic but inflate the degrees of freedom, thereby diluting power. Conversely, a small $K$ increases the risk that $\mathbb{E}[\pi_k (Z)] = 0$ for all $1\le k \le K$, rendering the test powerless. Hence, an appropriate choice of $K$ is crucial. In this section, we focus on the data-driven selection of $K$ for the three tests developed in Section (ref).

In Neyman's smooth test literature, the data-driven selection of $K$ was first proposed by ledwina1994data for testing uniformity. The data-driven testing procedure consists of two steps. First, Schwarz's selection rule (also known as the Bayesian Information Criterion, BIC) is applied to determine the dimension of the smooth model that best fits the data. Second, Neyman's smooth test is performed within the selected model space, yielding a data-driven test statistic. The data-driven smooth test thus combines preliminary model selection with a more precise inferential procedure, serving as Neyman's test in the “right” direction. Desirable theoretical properties and extensive numerical results of the selection rule and the induced test statistic were established in kallenberg1995consistency and kallenberg1995data, inglot1996asymptotic. Later, modifications of Schwarz's rule addressed more complex problems, including testing composite hypotheses kallenberg1997dataA, kallenberg1997dataB, inglot1997data, testing independence kallenberg1999data, and two-sample testing janic2000data, among others. These modified rules leverage smooth test statistics directly, avoiding the computation of the maximized log-likelihood, and are easier to implement.

Motivated by these strategies, we adopt a modified Schwarz's rule to determine the order $K$ for smooth tests in ANOVA models. The selection rule is specified as follows:

where $D$ is a fixed positive integer, and $\widehat{\Psi}_k^2$ takes the form of (ref) or (ref) depending on the specific scenario. The resulting data-driven test statistic is $N\widehat \Psi_{\widehat K}^2$, with $\widehat K$ given by (ref). The form of (ref) is mainly inspired by ducharme2004smooth, ducharme2004goodness, kraus2007data and duchesne2016estimating where the upper bound of selection $D$ is fixed. An alternative is to let $D = D(N) \to \infty$ as $N \to \infty$ kallenberg1995consistency, kallenberg1995data, inglot1996asymptotic, inglot1997data, kallenberg1997dataA, kallenberg1997dataB, kallenberg1999data, janic2000data. Although a diverging upper bound can, in principle, improve the consistency of smooth tests, the required rate of divergence is slow, and simulation studies show that empirical power levels off rapidly as $D$ increases. Hence, our proposed selection rule (ref) with fixed $D$ is both reasonable and easy to implement in practice.

To establish the theoretical properties of $\widehat K$ and $N\widehat{\Psi}_{\widehat K}^2$, we introduce a revised version of the alternative hypothesis $H_1$ (ref) and Assumption (ref). Consider

which encompasses a broader range of alternatives than $H_1$ by replacing the fixed order $K$ with the larger $D$. Correspondingly, the condition on the orthonormal system is extended to include all functions up to order $D$, as follows.

The following theorem establishes the properties of data-driven tests under the null and the revised alternatives.

Theorem (ref) presents the unified results of data-driven smooth tests for the three cases in Section (ref). Under the null hypothesis, the probability of $\{\widehat K = 1\}$ tends to $1$ asymptotically, implying that the first-order smooth model provides the best fit to $\{\widehat Z_{ij}\}_{i=1,j=1}^{N_j,J}$ among the $D$ candidate models. From the perspective of hypothesis testing, this means that $\widehat \Psi_1^2$ is informative and sufficient for assessing uniformity or normality. Meanwhile, the data-driven selection procedure allows identification of the mechanism underlying nonuniformity (equivalently, nonnormality) and enhances the test's power against a broader class of alternatives, $H_1^\prime$. These findings are consistent with previous results on data-driven smooth tests.

In practice, however, the $\chi_1^2$ limiting null distribution of data-driven smooth test statistics often performs poorly in finite samples kallenberg1995data, inglot1997data, kallenberg1997dataA, kallenberg1997dataB, kallenberg1999data, janic2000data, ducharme2004smooth, ducharme2004goodness, kraus2007data, duchesne2016estimating. To address this issue in our problem, and provided the nested orthonormal system (i.e. $\{\pi_k\}_{k=1}^K \subseteq \{\pi_k\}_{k=1}^{K+1}$ for $1\le K\le D-1$), we follow the approach of kallenberg1995data, inglot1997data, kallenberg1997dataA, kallenberg1999data, janic2000data, kraus2007data and adopt the following finite-sample approximation for the null distribution of $N \widehat{\Psi}_{\widehat K}^2$:

Such an approximation primarily accounts for the selection uncertainty of $\widehat K$. In finite samples, under the null, the empirical probability of the event $\{\widehat{K} = 1\}$ may not be exactly $1$, as $\{\widehat{K} = 2\}$ can occur with a small but non-negligible probability. Consequently, (ref) is derived from the approximation $\mathbb{P} (N \widehat{\Psi}_{\widehat K}^2 \le x) \approx \mathbb{P} (N \widehat{\Psi}_1^2 \le x, \widehat K = 1) + \mathbb{P} (N \widehat{\Psi}_2^2 \le x, \widehat K = 2)$. The reliability and accuracy of $H(x)$ are further illustrated through simulation studies in the next section and in the \hyperlink{app}{Online Appendix}.

Simulations

In this section, we conduct Monte Carlo experiments to examine the finite-sample performance of the smooth tests developed in Sections (ref) and (ref). Following neyman1937smooth, bera2013smooth, song2022smooth, and beutner2025two, we adopt the orthonormal Legendre polynomials on the interval $[0,1]$ for $\{\pi_k\}_{k = 1}^{\infty}$. Under this choice, the analytical properties of the constants $c_{1k}$ and $c_{2k}$ are provided in Proposition 1 of duchesne2016estimating, and numerical values for $1 \le k \le 10$ are available in Table 1 of their Supplementary Material. From a computational perspective, the test statistic $N\widehat{\Psi}_K^2$, introduced in Sections (ref) and (ref), can be computed using a simplified expression analogous to equation (14) in duchesne2016estimating, which leverages the closed form of $\bm{\Sigma}_K^{-1}$ (see equations (12) and (13) therein). For comparison, we include three classical normalty tests in the literature---Shapiro--Wilk (SW), Jarque--Bera (JB), and Kolmogorov--Smirnov (KS) tests---applied to the standardized residuals $\{\widehat e_{ij}\}_{i=1, j=1}^{N_j, J}$ under the corresponding model, without accounting for the parameter estimation effects.

We first study the tests for model (ref) through two experiments. The observations are generated as follows.

Here $\{Y_{ij}^{(0,1)}\}_{i=1,j=1}^{jm,5}$ and $\{Y_{ij}^{(0,2)}\}_{i=1,j=1}^{jm,5}$ are generated under $H_0$, while $\{Y_{ij}^{(1,1)}\}_{i=1,j=1}^{jm,5}$ and $\{Y_{ij}^{(1,2)}\}_{i=1,j=1}^{jm,5}$ are generated under $H_1$. For the proposed smooth tests, we evaluate the test statistic $N\widehat \Psi_K^2$ from Section (ref) with $1 \leq K \leq 5$, as well as the data-driven statistic $N\widehat \Psi_{\widehat K}^2$ from Section (ref), where $\widehat K$ is determined by (ref) with $D=5$. The sample size parameter $m$ ranges from 10 to 150 in increments of 10, and the significance level is set at $\alpha = 5\%$ with $500$ replications performed.

Results of Experiment \uppercase\expandafter{\romannumeral1} are presented in Tables (ref), (ref), and Figure (ref). Table (ref) reports the empirical rejection rates under $H_0$ for the proposed smooth tests and the classical normality tests. The smooth tests comprise both fixed-$K$ versions ($K=1,\ldots,5$) and data-driven procedures, implemented using the limiting $\chi^2_1$ null distribution ($\widehat K\&\chi^2_1$) or the approximated null distribution $H(x)$ ($\widehat K\&H(x)$). For fixed $K$, under $H_0$, the statistic $N\widehat \Psi_K^2$ maintains the nominal level well, even for small sample sizes such as $m=10$ (corresponding to a total sample size of $N=150$). For the data-driven test statistic $N\widehat \Psi_{\widehat K}^2$, the limiting $\chi^2_1$ distribution fails to control the rejection rate under $H_0$, whereas the approximated distribution $H(x)$ achieves better size control, yielding rejection rates close to the nominal $5\%$ level. Under $H_1$, all versions of the proposed smooth tests exhibit high power across the considered settings. As for the classical methods, although the Shapiro--Wilk and Jarque--Bera tests are applied directly to the standardized residuals without accounting for the estimation effects, their empirical rejection rates under $H_0$ remain close to the nominal $5\%$ level. Under $H_1$, both tests achieve high power, indicating that they still provide reasonably reliable numerical performance despite the lack of formal adjustment for the estimation effect. In contrast, under $H_0$, the Kolmogorov--Smirnov test exhibits almost zero rejection across all settings, indicating the practical failure of the test. Table (ref) presents the empirical frequency of $\widehat K$ under $H_0$ and $H_1$, respectively. Under $H_0$, the frequency of $\{\widehat K=1\}$ approaches $1$ as $m$ increases, whereas under $H_1$, $\widehat K$ is consistently greater than $1$, in line with theoretical expectations. Figure (ref) displays the sample means of $\widehat K$ across different values of $m$, with error bars representing one sample standard deviation around the sample means. Under $H_0$, the sample mean of $\widehat K$ stays close to $1$ with diminishing variability as $m$ grows, indicating consistent model selection. Under $H_1$, $\widehat K$ increases steadily and approaches the maximum $D=5$, with reduced variability, demonstrating the adaptiveness of the data-driven procedure.

The results of Experiment \uppercase\expandafter{\romannumeral2} are quite similar; see Tables (ref), (ref), (ref), and Figure (ref). However, from Table (ref), under $H_1$, the test statistic $N\widehat \Psi_1^2$ fails to exhibit power as a result of $\mathbb{E}[\pi_1(Z)] = \mathbb{E}[\pi_1(\Phi(Y_{ij}^{(1,2)} - 8))] = 0$ in this case. This illustrates that $N\widehat \Psi_K^2$ may be inconsistent against certain alternatives, especially when $K$ is small. Correspondingly, both Table (ref) and Figure (ref) show that under $H_1$, $\{\widehat K = 1\}$ does not occur, consistent with Theorem (ref) (1), which ensures that $\mathbb{P}(\widehat K \ge K_0) \to 1$ under $H_1^\prime$ with $K_0 = 2$.

The next two experiments, focusing on models (ref) and (ref), follow simulation settings similar to the previous ones. The data-generating process for each model is specified as follows.

Here $\{Y_{ij}^{(0,3)}\}_{i=1,j=1}^{jm,5}$ and $\{Y_{ij}^{(0,4)}\}_{i=1,j=1}^{jm,5}$ are generated under $H_0$, while $\{Y_{ij}^{(1,3)}\}_{i=1,j=1}^{jm,5}$ and $\{Y_{ij}^{(1,4)}\}_{i=1,j=1}^{jm,5}$ are generated under $H_1$. Note that in Experiment \uppercase\expandafter{\romannumeral3}, $\mathbb{E} [Y_{ij}^{(0,3)}] = \mathbb{E} [Y_{ij}^{(1,3)}] = 5j$ and $\operatorname{Var} [Y_{ij}^{(0,3)}] = \operatorname{Var}[Y_{ij}^{(1,3)}] = 4$, which corresponds to model (ref); in Experiment \uppercase\expandafter{\romannumeral4}, $\mathbb{E} [Y_{ij}^{(0,4)}] = \mathbb{E} [Y_{ij}^{(1,4)}] = 8$ and $\operatorname{Var}[Y_{ij}^{(0,4)}] = \operatorname{Var} [Y_{ij}^{(1,4)}] = j^2$, which corresponds to model (ref). For the smooth tests, the results of Experiments \uppercase\expandafter{\romannumeral3} and \uppercase\expandafter{\romannumeral4} are broadly consistent with those of Experiments \uppercase\expandafter{\romannumeral1} and \uppercase\expandafter{\romannumeral2} (see Tables (ref)--(ref) and Figures (ref)--(ref)), further supporting the validity of our unified testing framework. For the other three classical methods, a notable phenomenon arises in Experiment \uppercase\expandafter{\romannumeral4}: under $H_0$, the Jarque--Bera test exhibits undersizing, suggesting that it is not well suited for models with heterogeneous group variances. In contrast, our proposed test remains well-calibrated and maintains stable size control.

To further assess robustness with respect to the number of groups $J$, we also consider a variant of Experiment \uppercase\expandafter{\romannumeral3} with $J=10$ (denoted as Experiment \uppercase\expandafter{\romannumeral3}$^\prime$). The results corroborate our main conclusions; see Section (ref) of the \hyperlink{app}{Online Appendix} for details. In summary, the experiments demonstrate that the proposed smooth tests in Section (ref) maintain the nominal significance level well under the null, and exhibit high power under the alternatives except when $\mathbb{E}[\pi_k(Z)] = 0$ for all $1 \le k \le K$. The data-driven tests with the approximated null distribution $H(x)$ reliably control the Type \uppercase\expandafter{\romannumeral1} error and achieve excellent power performance. In contrast with the three classical methods, our approach is supported by solid theoretical guarantees and demonstrates uniformly stable and competitive performance across all settings considered.

An empirical application

In this section, we illustrate the proposed smooth tests using data from the OECD Programme for International Student Assessment (PISA) 2018. PISA is a large-scale international assessment coordinated by the Organisation for Economic Co-operation and Development (OECD) that evaluates educational systems by measuring students' learning outcomes and school characteristics across countries schleicher2019pisa. The dataset used in this study is publicly available at \url{https://www.kaggle.com/datasets/dilaraahan/pisa-2018-school-questionnaire}. Since PISA data are widely used to compare educational resources across countries, statistical inference for group means is commonly employed, the validity of which depends crucially on distributional assumptions such as normality. However, educational indicators, including school size and student–teacher ratios, often exhibit substantial cross-country heterogeneity and dispersion, and prior studies have documented skewed distributions in school-level resource measures hanushek2011economics,chandir2022student. Deviations from normality may therefore affect variance estimation and the reliability of ANOVA-based inference, making it important to formally assess this assumption.

We consider two response variables, STRATIO (student--teacher ratio) and SCHSIZE (school size), measured at the school level, and use CNTRYID (country identifier) to define the grouping structure. Table (ref) summarizes key sample characteristics after removing missing values, including the total number of observations $N$, the number of groups $J$, and the minimum group size $\min_j N_j$. The resulting group structure satisfies the conditions required by the asymptotic framework in Section (ref) (see Table (ref)), thereby supporting the application of the proposed models and normality tests. For each response variable (STRATIO and SCHSIZE) under the grouping structure defined by CNTRYID, we fit the ANOVA models (ref), (ref), and (ref), and apply the corresponding smooth tests for normality. For each model, we consider two implementations of the smooth test: one with fixed $K=4$, and one based on the data-driven procedure with $D=5$.

The results show that, for both response variables, STRATIO and SCHSIZE, the null hypothesis of normality is rejected at the $0.1\%$ significance level across all three model specifications and both implementations of the proposed tests. This uniform rejection pattern indicates strong robustness of the findings. Overall, there is compelling evidence that the distributions of STRATIO and SCHSIZE deviate substantially from normality. These deviations have important implications for PISA data, where ANOVA models are commonly used to compare educational resources across countries. For variables such as student--teacher ratios and school size, pronounced skewness or heavy-tailed behavior likely reflects substantial heterogeneity in educational systems, with a small number of countries or schools exhibiting extreme values. In such settings, classical $F$-tests may be sensitive to these distributional features, potentially leading to distorted inference when comparing group means. Consequently, ANOVA-based comparisons of educational resources should be interpreted with caution, and more robust approaches, such as transformations, rank-based methods, or other distribution-free procedures, may be preferable in practice. More broadly, these findings underscore the importance of formally assessing distributional assumptions prior to conducting inference in cross-country educational studies.

Concluding remarks

In this paper, we propose and examine data-driven Neyman's smooth tests for assessing the normality assumption in ANOVA models with a potentially diverging number of groups. For three types of one-way fixed effects models, we derive the asymptotic properties of the proposed tests and validate their finite-sample performance through extensive numerical studies. Our results provide a rigorous and practical tool for evaluating normality in a broad range of ANOVA settings.

Several directions remain open for future research. First, while the proposed tests are developed for specific ANOVA models, the true data-generating mechanism may not be known in practice. A natural extension is to combine our procedures with preliminary structural tests to identify the most appropriate ANOVA specification before applying the proposed methodology. Second, the assumption of independent random errors can be relaxed. For example, repeated clinical trials often involve correlated errors ma2012beyond, schober2018repeated, langenberg2022repeated. In such cases, our smooth test framework may be adapted by incorporating estimated correlation structures. Finally, for the selection rule (ref), it is of theoretical interest to explore the scenario where the upper bound $D = D(N)$ diverges slowly with $N$, that is, $D \to \infty$ as $N \to \infty$.

\hypertarget{app}