An anomaly arising in the analysis of processes with more than one source of variability

Summary

It is frequently observed in practice that the Wald statistic gives a poor assessment of the statistical significance of a variance component. This paper provides detailed analytic insight into the phenomenon by way of two simple models, which point to an atypical geometry as the source of the aberration. The latter can in principle be checked numerically to cover situations of arbitrary complexity, such as those arising from elaborate forms of blocking in an experimental context, or models for longitudinal or clustered data. The salient point, echoing Dickey (2020), is that a suitable likelihood-ratio test should always be used for the assessment of variance components.

1 Introduction

Faithful representation of a process generating data often entails specification of two or more sources of variability. In an experimental context, simple or elaborate forms of blocking induce a nested or crossed structure within the set of plots. Similar grouping arises in observational studies where, for instance, data may originate from different hospitals, regions or from several family groups, of no direct interest, but likely to generate structured correlation in the outcome.

As in certain other settings, inference based on the likelihood function for the full generative model is typically miscalibrated, sometimes seriously so, and should ideally be based on a suitable marginal or conditional likelihood. In the present context, appropriate preliminary reduction leads to residual maximum likelihood, REML, developed by Patterson & Thompson (1971), and closely connected to marginal likelihood (Bartlett, 1937).

Even when the REML likelihood is used, the Wald test routinely reported in software implementations is frequently found to be ineffectual for detecting components of variance when they are unambiguously present. A typical example occurs in McCullagh (2023, Ch. 4), where the likelihood-ratio statistic is more than eight times as large as the squared Wald ratio. The phenomenon has also been documented by Dickey (2020), who presented examples in which the REML Wald statistic is bounded. In at least one of the examples he considered, the natural estimate of standard error of the REML estimator is proportionate to the REML estimate itself, so that the role of the data in the Wald construction is negligible. The purpose of the present paper is to expose the matter in a form amenable to detailed analytic calculation, thereby revealing dependence on key aspects and indicating other settings in which the same phenomenon is inevitable. The insights of Dickey (2020) are reproduced and elucidated further. The source of the anomalous behaviour is not a failure of distributional approximations obtained under hypothetical regimes of sometimes questionable adequacy, but rather atypical geometry of the REML loglikelihood function, which induces a bounded Wald statistic even under a notional limiting operation in which the likelihood-ratio statistic for testing the same hypothesis is arbitrarily large. We show that a version of the score statistic is subject to the same aberration.

The implications for applied work are consequential, as undetected sources of variability typically result in estimated standard errors for regression coefficient estimators that are deceptively small, and therefore confidence intervals that are misleadingly narrow. The opposite situation can occasionally arise as well. For instance, in a randomized block design, omission of the block factor as a variance component has the effect of increasing the estimated variance of treatment effect estimators.

2 Variance-component models

In a typical variance-component model, the distribution of the response vector is specified by a mean vector $\mu =X\beta$ in the linear subspace $\mathcal{X}$ spanned by the columns of the model matrix X, and a covariance matrix Σ in the convex cone

spanned by given matrices $V_0,\dots ,V_s$⁠, which are positive definite or semidefinite. Usually, $V_0=I_n$ is the identity; the remaining matrices may be block factors or structured matrices associated with spatial or temporal dependence.

The residual likelihood is the likelihood function based on the residual $U^{\mathrm{T}}Y$⁠, where $\ker (U^{\mathrm{T}})=\mathcal{X}$⁠. Provided that the matrices $U^{\mathrm{T}}{V}_{u}U$ are linearly independent, the variance components θ_u may be estimated by maximizing the residual likelihood. In the subsequent discussion, reference to the loglikelihood function and its maximizer means the REML version unless otherwise specified.

There are compelling general arguments, notably invariance and permissibility of asymmetric confidence regions, for basing an assessment of the hypothesis $H_0:{\theta }_s=0$⁠, say, on the likelihood-ratio statistic

where $\widehat{\theta }$ is the maximum likelihood estimator and ${\widehat{\theta }}^{(0)}$ is the constrained estimator under H₀. Nonnegativity of the variance coefficients means that the subset of $\mathcal{V}$ defined by the constraint ${\theta }_s=0$ is a boundary subcone, with the implication that the likelihood achieves its maximum on the boundary with positive probability, usually one half for sufficiently large sample size (Chernoff, 1954). When an exact F statistic exists, the boundary event occurs whenever $F ⩽ 1$⁠. On this event, $\widehat{\theta }\in H_0$ coincides with ${\widehat{\theta }}^{(0)}$ and the loglikelihood ratio is exactly zero. Thus, with asymptotic probability one half, Λ = 0; otherwise, its distribution under H₀ is ${\chi }_1^2$ under suitable limiting conditions. The realized value of Λ is to be calibrated against this distribution.

The nominal asymptotic variance of ${\widehat{\theta }}_s$ is $i^{ss}(\theta )$⁠, the (s, s) component of the inverse Fisher information matrix. With asymptotic probability one half, the squared Wald ratio $W^2={\widehat{\theta }}_s^2/i^{ss}({\widehat{\theta }}_s)$ is equal to zero under H₀; otherwise, it is asymptotically equivalent to Λ by a standard asymptotic argument. However, it is frequently observed in practice that W² gives a poor assessment of the statistical significance of the sth variance component, in the sense that it fails to reject H₀ even when Λ does so unambiguously at the same significance level. If the likelihood is maximized on the boundary, both W² and Λ are zero and there is no disagreement. Disagreement can only occur when both are positive, and it is that phenomenon that we address here.

Among the substantial literature on testing variance components is an early contribution by Wald (1947), notable in that it recommends an F test but fails to warn against use of the eponymous statistic in this context. Subsequent contributions have developed nonstandard asymptotic theory for the likelihood-ratio statistic and modifications thereof, relaxing an earlier assumption that the true parameter value belongs to the interior of the parameter space (e.g., Self & Liang, 1987; Geyer, 1994; Vu & Zhou, 1997). While the likelihood-ratio test and its variants have been the focus of theoretical development, common software implementations report Wald statistics as standard, without warning. Dickey (2020) appears to have been the first to emphasize the point at issue. The popularity of Wald-based inference perhaps stems from the convenience of its construction, requiring a single maximum likelihood fit in contrast with two for Λ, which facilitates the presentation of confidence statements.

The analysis of § 3 and § 6 is comparable to that of Dickey (2020), whose derivations cover situations in which the likelihood-ratio test coincides with an exact F test. The two papers illustrate the aberration under study in different ways and § 3.4 provides a comparison and synthesis. Sections 4 and 5 cover situations in which a fruitful formulation in terms of F may be infeasible, but for which the shared anomalous geometry can be checked by direct study of the loglikelihood function. Together, the present paper and that of Dickey (2020) provide a thorough explanation of a phenomenon of broad relevance and scientific consequence.

3 Analysis for a single block factor

3.1 Introduction

We consider in this section a simple Gaussian model with two variance components estimated from the sufficient statistic, which consists of the within-block mean square ${\text{MS}}_0$ on f₀ degrees of freedom and the between-block mean square ${\text{MS}}_1$ on f₁ degrees of freedom. Specifically, the outcome is $Y_{ji}=\mu +{\eta }_j+{\epsilon }_{ji}$ (⁠$j=1,\dots ,k,\hspace{0.17em}i=1,\dots ,b$⁠), where, for an arbitrary block index j, ${(Y_{j1},\dots ,Y_{jb})}^{\mathrm{T}}$ has a Gaussian distribution of mean $\mu 1_b$ and covariance matrix $\Sigma ={\sigma }_0^2{I}_b+{\sigma }_{\eta }^2{1}_b{1}_b^{\mathrm{T}}$⁠. Define, in a standard notation for averaging over suffixes, ${\overline{Y}}_{j •}=\sum _i{Y}_{ji}/b$⁠, and similarly for the double average. The between-block sum of squares

is distributed as ${\sigma }_1^2{\chi }_{f_1}^2$⁠, where $f_1=k-1$ and, in the usual variance-component parameterization with $\theta ⩾ 0$⁠,

The within-block sum of squares $\sum _{j,i}{(Y_{ji}-{\overline{Y}}_{j •})}^2=f_0{\text{MS}}_0$ is ${\sigma }_0^2{\chi }_{f_0}^2$ distributed independently of the between-block sum of squares, where $f_0=k(b-1)$⁠. The constraint on θ means that the null hypothesis $H_0:\theta =0$ of equality of variances is on the boundary.

In the balanced one-way analysis of the variance structure of the above formulation, the REML loglikelihood is the marginal loglikelihood $\ell (\theta ,{\sigma }_0^2)$ based on the joint density function of $({\text{MS}}_0,{\text{MS}}_1)$⁠. While numerous generalizations may be considered, the simple version presented here isolates the point at issue in the most incisive form, free of secondary effects that complicate the analysis and interpretation.

3.2 Comparison of Wald and likelihood-ratio statistics

Maximization of $\ell (\theta ,{\sigma }_0^2)$ without the constraint $\widehat{\theta } ⩾ 0$ produces estimators $\widehat{\theta }=(F-1)/b$ and ${\widehat{\sigma }}_0^2={\text{MS}}_0$⁠, where $F={\text{MS}}_1/{\text{MS}}_0$ is Fisher’s F ratio, whose distribution depends only on the variance ratio θ. The maximum likelihood estimator $\widehat{\theta }$ has nominal asymptotic variance given by the relevant diagonal component of the inverse Fisher information matrix, namely,

where $f=f_1+f_0$⁠. The squared Wald statistic for testing θ = 0 is therefore

to be compared with the likelihood-ratio statistic

where the pooled mean square $\text{MS}=(f_1{\text{MS}}_1+f_0{\text{MS}}_0)/f$ is the maximum likelihood estimator of ${\sigma }_0^2$ under the constraint θ = 0.

Although the statistics W² and Λ are known functions of $b\widehat{\theta }=(F-1)$⁠, simple approximations provide insight into the nature of the aberration described in § 1. A Taylor expansion of W² and Λ around $\widehat{\theta }=0$ gives

or, in terms of $\Lambda ⩾ 0$⁠,

showing that they agree up to second order in $\widehat{\theta }$ and to first order in Λ, but not beyond. Typically, $f_0\gg f_1$ is large, in which case the ratio is approximately

for small $\widehat{\theta }$ or Λ.

The squared Wald statistic W² is justified based on standard asymptotic theory for maximum likelihood estimators. First- and second-moment theory suggests three further Wald statistics, for which the same anomalous behaviour is demonstrated in the Supplementary Material. For a similar discussion of two score statistics, see also the Supplementary Material.

In the simplest setting, at least, both W and ${\Lambda }^{1/2}\text{sign}(F-1)$ are monotone functions of the mean-square ratio F, so their distributions are known exactly as a function of the variance-ratio parameter θ. The transformations $\Lambda =g_1(F)$ and $W^2=g_2(F)$ are strictly monotone increasing for F > 1, and decreasing for F < 1. The discrepancies described here arise from the standard practice of treating Λ and W² as if they were asymptotically identical statistics rather than equivalent statistics. The standard practice is justified only if g₁ = g₂, at least approximately for large samples, which is not the case in typical variance-components models.

One definition of Wald-detectability equates θ to ${\Phi }^{-1}(1-\alpha )$ times the estimated standard error of $\widehat{\theta }$⁠. It can be shown that there does not exist a positive Wald-detectable value in this sense unless the number of blocks is large. It is arguably more natural to compute standard errors under the null, in which case the values at the borderline of Wald detectability are

at the 16% and 2% levels. The corresponding thresholds in terms of F are

With $f_0=102,\hspace{0.17em}{f}_1=5$ and b = 18, the ratio $\Lambda /W^2$ is approximately 1.864 at ${\theta }_{16}^{*}$ and 2.727 at ${\theta }_2^{*}$⁠. With the same numbers, $F_{16}^{*}=1.648$ and $F_2^{*}=3.296$⁠, to be compared with the 16% and 2% critical values of the F distribution on (f₁, f₀) degrees of freedom, which are 1.625 and 2.818, respectively. Thus, even when standard errors are computed under the null hypothesis, the observed value of the F statistic has to be 17% larger than the 2% critical value of the F distribution in order for the Wald test to reject at the same level. Equivalently, rejection at the 2% level using a Wald statistic with standard errors computed under the null requires an observed value of $\widehat{\theta }$ that is double that necessary for rejection at the same level using an F test. These latter conclusions do not involve any approximations, and are similar and complementary to those of Dickey (2020).

Geometric insight is obtained by noting that the Wald statistic W² implicitly defines a quadratic approximation $q(\theta )$ to the profile REML loglikelihood function, $\ell (\theta ;{\widehat{\sigma }}_{\theta }^2)$ in a neighbourhood of $\widehat{\theta }=(F-1)/b$⁠. Specifically, for fixed θ, the maximum likelihood estimator of ${\sigma }_0^2$ is ${\widehat{\sigma }}_{\theta }^2=w_0{\text{MS}}_0+w_1{\text{MS}}_1/(1+b\theta )$⁠, where $w_r=f_r/f$ and

Write $\ell (\theta ;{\widehat{\sigma }}_{\theta }^2)=\widehat{\ell }(\theta )+f/2$⁠. The quadratic approximation to $\widehat{\ell }(\theta )$ implicit in the Wald test is

where $f_0{f}_1{b}^2/(2f{F}^2)$ is $1/i^{\theta \theta }(\widehat{\theta })$ and

Since the two functions have the same value at $\widehat{\theta }$⁠, the discrepancy between the Wald and likelihood-ratio statistics for testing θ = 0 is the difference in y intercepts:

with the expression for Λ in terms of F given by

For F and f₁ fixed,

showing that the discrepancy is roughly linear in F for large f₀, while for fixed f₁ and f₀, (4) converges to zero as F approaches unity and is unbounded for F arbitrarily large. Figure 1 graphs $q(\theta )$ and $\widehat{\ell }(\theta )$ for different values of F.

3.3 Nonconstant Fisher information and anomalous geometry

From (3), $\widehat{\ell }(\theta )$ has second derivative

whose value at $\widehat{\theta }$ is

showing that the curvature at the maximum-likelihood point is close to zero for large F, as depicted in Fig. 1. In other words, $\widehat{\ell }(\theta )$ is arbitrarily well approximated by a horizontal asymptote at $\widehat{\ell }(\widehat{\theta })$ in a neighbourhood of $\widehat{\theta }$ for arbitrarily large F. Equation (5) also shows that the discrepancy $q(0)-\widehat{\ell }(0)$ is attributable to the higher-order derivatives, since there is no error incurred by using $-1/i^{\theta \theta }(\widehat{\theta })$ in place of $\gamma (\widehat{\theta })$ in the Taylor series approximation to $\widehat{\ell }(\theta )$⁠. The effect of higher-order derivatives is encapsulated to a large extent in the considerable nonconstancy of $i^{\theta \theta }(\theta )$ as a function of θ over the range of interest.

From (1), the ratio of the nominal asymptotic variances of $\widehat{\theta }$ at arbitrary θ and at θ = 0 is $i^{\theta \theta }(\theta )/i^{\theta \theta }(0)={(1+b\theta )}^2$⁠, and the range of primary interest is

the upper value being the nominal asymptotic standard deviation of $\widehat{\theta }$ under the null hypothesis, having conditioned on the event that $\widehat{\theta }$ is not on the boundary. Over this range, the asymptotic variance varies by a factor

which is large for typical values of f₁. For example, $f_0=102,\hspace{0.17em}{f}_1=5$ gives a range of approximately $1 ⩽ {(1+b\theta )}^2 ⩽ 2.72$⁠.

3.4 A synthesis with Dickey (2020)

The motivation for the present paper came from practical examples in McCullagh (2023, Ch. 4), where the Wald statistic was ineffectual at detecting variance components, the absence of which was strongly refuted by a likelihood-ratio test. Dickey (2020) exposed the same phenomenon. His exposition, aimed at practitioners, covers widely used models in the analysis of designed experiments for which exact F and Wald tests are available. The present paper is framed in terms of the loglikelihood-ratio statistic Λ, which is more generally available and points to geometric insights not recoverable from the moment-based Wald constructions. Three instances of the latter are discussed in the Supplementary Material, among which equation (S.1) coincides with equation (2) of Dickey (2020). For cases where F is available, the present paper and that of Dickey provide equivalent explanations from two points of view.

4 Two-sample problem in general scale families

Let Y be a random variable from a scale family with density function ${\tau }^{-1}g(y/\tau ),y,\tau >0$⁠, where g is a known, continuous, density function on the positive real line. Let

Then the Fisher information for τ in an independent and identically distributed sample $Y_1,\dots ,Y_n$ is $n{I}_g/{\tau }^2$ and that for ${\tau }_0,{\tau }_1$ in a two-sample problem of sizes n₀ and n₁ is

The squared Wald statistic for testing equality of scale parameters is therefore

where $\widehat{\theta }$ is the estimated ratio ${\widehat{\tau }}_1/{\widehat{\tau }}_0$⁠. The Wald statistic is bounded in both limits for $\widehat{\theta }$ arbitrarily large or small. Specifically,

The curvature of the profile loglikelihood function at $\widehat{\theta }$ is approximated by the expected Fisher information in the profile loglikelihood. The Fisher information transforms as

where ψ and $\varphi$ are two parameterizations and we have used the convention that symbols appearing both as subscripts and superscripts in the same product are summed. The information about θ, having adjusted for estimation of τ₀ is therefore, using (6) and (7),

showing that the two-sample problem in general scale families has the same anomalous geometry documented in § 3 for large θ, leading to the discrepancy between the likelihood-ratio and the Wald statistics.

5 Gaussian variance-component model

Consider a variance-component model in which $Y\in {\mathbb{R}}^n$ is normal with mean $\mu \in \mathcal{X}$ of dimension p < n, and covariance matrix $\Sigma ={\sigma }^2\{I_n+V(\theta )\}$⁠, where $V(\theta )$ is a known matrix function of a vector parameter $\theta ={({\theta }_1,\dots ,{\theta }_s)}^{\mathrm{T}}$ with $V(0)=0$⁠. This encompasses the linear model $V(\theta )=\sum _u{\theta }_u{V}_u$ from § 2. The subspace $\mathcal{X}\subset {\mathbb{R}}^n$ is a group under addition, which implies that, for any matrix $U^{\mathrm{T}}$ with kernel $\mathcal{X}$⁠, the normalized residual statistic

is maximal invariant under the affine group with action $g(a,x):y\mapsto ay+x$⁠, with a > 0 and $x\in \mathcal{X}$⁠. For distributional calculations, it is convenient to take the columns of U to be an orthonormal basis in ${\mathcal{X}}^{\perp }$⁠, the orthogonal complement of $\mathcal{X}$ with respect to the standard inner product. In that case $U{U}^{\mathrm{T}}=I_n-X{(X^{\mathrm{T}}X)}^{-1}{X}^{\mathrm{T}},\hspace{0.17em}{U}^{\mathrm{T}}U=I_{n-p}$ and the density function of $Q=q(Y,U)$ is

where ${\sigma }^2{A}_{\theta }=U^{\mathrm{T}}\Sigma U$⁠. At θ = 0, Q is uniformly distributed on the $(n-p)$-dimensional unit sphere in ${\mathbb{R}}^n$⁠. The above exposition hybridizes Kariya (1980) and King (1980).

By construction, distribution (8) does not depend on β or ${\sigma }^2$⁠, and inference for θ is conveniently based on the marginal loglikelihood function

closely related to the REML loglikelihood function that uses the marginal distribution of $U^{\mathrm{T}}Y$ rather than Q. Since the transformation to Q eliminates the nuisance parameter ${\sigma }^2$ as well as β, analysis based on (9) is more amenable to analytic calculation.

The density function in (8) is relative to a particular orthonormal basis in ${\mathcal{X}}^{\perp }$⁠, and the same basis is embedded in the likelihood function (9) in matrix $A_{\theta }$⁠. The conclusion in (9) is equivalent to equation (3.3) of McCullagh (2009), which evades the problem of selecting a basis by allowing singular matrices.

The single block factor setting of § 3 is a special case with $n=(f_1+1)b,\hspace{0.17em}\Sigma ={\sigma }^2(I_n+\theta V),\hspace{0.17em}V=I_{f_1+1}\otimes 1_b^{}{1}_b^{\mathrm{T}}$ and $\mathcal{X}=1$⁠, so that $U{U}^{\mathrm{T}}=I_n-1_n^{}{1}_n^{\mathrm{T}}/n$⁠. Suppose for an exact analytic calculation that $f_1+1$ and b are both powers of two. The Supplementary Material shows that $\text{log}\hspace{0.17em}\left|A_{\theta }\right|=f_1\hspace{0.17em}\text{log}(1+b\theta )$ and

so that (9) becomes

The solution of the likelihood equation by differentiation of (10) is $(1+b\stackrel{⌣}{\theta })=F$ and direct calculation shows that $\stackrel{⌣}{\ell }(\theta )-\stackrel{⌣}{\ell }(\stackrel{⌣}{\theta })=\widehat{\ell }(\theta )-\widehat{\ell }(\widehat{\theta })$ from (3). This is demonstrated empirically in Fig. 2 with the values of f₀, f₁ and b from the numerical example of § 6.1 so as to also illustrate the anomalous geometry. When an analysis of variance is feasible, it produces identical estimates to those based on maximization of (9), the only difference arising from the choice of basis for the orthogonal subspaces.

More generally, in a model with a scalar parameter θ generating a variance component in the form $\Sigma ={\sigma }^2\{I_n+V(\theta )\}$⁠, the likelihood equation for θ is

at $\theta =\stackrel{⌣}{\theta }$⁠, where ${\dot{A}}_{\theta }={\nabla }_{\theta }{A}_{\theta }$ and, by the Woodbury identity,

The second derivative of $\stackrel{⌣}{\ell }(\theta )$ at an arbitrary point is

where ${\ddot{A}}_{\theta }={\nabla }_{\theta \theta }^2{A}_{\theta }$⁠. A general analytic approximation to the curvature $\stackrel{⌣}{\gamma }(\stackrel{⌣}{\theta })$ has not been ascertained. However, if $V(\theta )$ is of the form θV with V a known matrix, then ${\dot{A}}_{\theta }=U^{\mathrm{T}}VU,\hspace{0.17em}{\ddot{A}}_{\theta }=0$ and the previous display simplifies. In particular,

and

so that $\underset{\theta \to \infty }{\lim }{\nabla }_{\theta \theta }^2\stackrel{⌣}{\ell }(\theta )=0$⁠. By consistency of the maximum likelihood estimator, the marginal loglikelihood function has zero curvature at the maximising point when the true value of θ is arbitrarily large.

6 Numerical illustrations

6.1 Covariance determined by split-plot nested blocking

The split-plot covariance ${\sigma }_0^2(I_n+\theta V)$ is a linear combination of the identity and a binary matrix such that V_ij = 1 if observational units i, j belong to the same whole plot or block. Chapter 1 of McCullagh (2023) discusses an example of this type with l = 24 rats constituting the blocks, and s = 5 sites on each rat constituting the observational units. A two-level treatment was assigned at random to the rats, so the model formula site+treat for the expected value determines a subspace of dimension 6. This is not a split-plot design in the traditional sense because the split-plot effect is associated with a classification factor, sites ranging from anterior to caudal, not with a treatment factor.

The real experiment is a little more complicated because some components are missing. With this exception, the simulated data mirror that example, and illustrate the effect of increasing θ on the discrepancy between Λ and W². Since $\widehat{\theta }\simeq 0.4$ for the actual experiment, the range $0 ⩽ \theta ⩽ 1$ was used for simulation. For each of the 1000 Monte Carlo replications, outcomes were generated according to the model $Y\sim N(\mu ,\Sigma )$⁠, the particular point $\mu \in \mathcal{X}$ being immaterial for the REML likelihood. As it happens, the analysis of § 3 applies here with two modified mean squares, one for treatment and one for residuals eliminating additive row and column effects, namely,

where $f_1=l-2,\hspace{0.17em}{f}_0=(l-1)(s-1)$ and, for example, $P_{\text{trt}}Y$ is the projection of Y on the subspace spanned by the treatment basis. For the null hypothesis θ = 0, Fig. 3 compares the simulated values of $\Lambda /W^2$ with the linear approximation $1+4s\widehat{\theta }/3$ from (2).

Maximization of the marginal loglikelihood function $\stackrel{⌣}{\ell }(\theta )$ from (9) is equivalent. Both normed loglikelihood functions are graphed in Fig. 2 for $\widehat{\theta }=\stackrel{⌣}{\theta }=0.6$⁠, from which the anomalous geometry is apparent.

6.2 Covariance determined by Latin square blocking

The Latin square covariance $\Sigma ={\sigma }_0^2(I_n+{\theta }_1{V}_1+{\theta }_2{V}_2)$ is a linear combination of the identity and two binary matrices of the form ${(V_1)}_{ij}=1$ if observational units i and j share a column and ${(V_2)}_{ij}=1$ if they share a row. For the simulation, the variance-component parameters were taken as ${\sigma }_0^2=1,\hspace{0.17em}{\theta }_2=0.2$⁠, while θ₁ was varied in the range $0\le {\theta }_1\le 1$⁠. The relevant mean squares on $f_0=(b-1)(b-2)$ and $f_1=(b-1)$ degrees of freedom are

where ${\sigma }_1^2={\sigma }_0^2(1+b{\theta }_1)$⁠. The analogous mean square ${\text{MS}}_2$ for rows on f₂ = f₁ degrees of freedom is used in variance estimation. Specifically, the information about θ₁, having adjusted for estimation of ${\sigma }_0^2$ and θ₂, is

In (11), θ_j is estimated as the positive part of $(F_j-1)/b$ with $F_j={\text{MS}}_j/{\text{MS}}_0$⁠. For $f=f_1+f_0\gg f_1$⁠, which amounts to b being large, the adjustment is negligible provided that θ₂ is not too large relative to θ₁, and (11) is comparable to (1). The resulting Wald statistic for testing the hypothesis ${\theta }_1=0$ is to be compared with Λ, whose form is identical to that of § 3.

Figure 4 depicts qualitatively similar behaviour for the ratio $\Lambda /W^2$ as that of Fig. 3, with additional variability attributable to randomness in ${\widehat{\theta }}_2$ manifesting through the estimate of (11). The version with θ₂ treated as known is depicted in Fig. 4(b).

Acknowledgement

Section 3.4 was based on a detailed comparison supplied by one of the two referees, to whom we are grateful.

Supplementary material

The Supplementary Material provides more explicit derivations for some of the key results, together with a discussion of three further Wald statistics and two score statistics.

References