# Testing for Main Random Effects in Two-Way Random and Mixed Effects Models: Modifying theF Statistic

Authors: Trent Gaugler; Michael G. Akritas
Journal: Journal of Probability and Statistics
 (2013)

---

## Abstract

A procedure for testing the significance of the main random effect
is proposed under a model which does not require the traditional assumptions of symmetry, homoscedasticity, and normality for the error term and
random effects. To accommodate this level of model generality, and also unbalanced designs, suitable adjustments to theF-test are made. The extensive
simulations performed under the random effects model, and the unrestricted
and restricted versions of the mixed effects model, indicate that the classical
F procedure is extremely liberal under heteroscedasticity and unbalancedness.
The proposed test procedure performs well in all settings and is comparable
to the classical F-test when the classical assumptions are met. An analysis of
a dataset from the Mussel Watch Project is presented.


---

## Body

## 1. Introduction

Consider a two-factor mixed or random effects design, and letYijk denote the kth observation at level i of the row factor and level j of the column factor. In the case of a mixed effects design we assume that the row factor is fixed.The classical two-way model, compare [1–5], uses the decomposition
(1)Yijk=μ+αi+βj+γij+ϵijk.
For the random effects model, the αi’s, βj’s, γij’s, and ϵijk’s are mutually independent, the αi’s are iid N(0,σα2), the βj’s are iid N(0,σβ2), the γij’s are iid N(0,σγ2), and the ϵijk’s are iid N(0,σϵ2).For the mixed effects model, there are two common definitions of the effects. Theunrestricted version of the model assumes
(2)∑i‍αi=0,βjareiidN(0,σβ2),γijareiidN(0,σγ2),
and the βj,γij are all independent. The restricted version of the model keeps the above assumptions but requires that
(3)∑i‍γij=0,
implying that γij,γi′j are correlated. Both versions assume the ϵijk to be iid N(0,σϵ2) and independent from the other random effects. There is a dichotomy of opinion in the statistical literature as to which model should be used. Cornfield and Tukey [6], Scheffé [1], Winer [7], and Khuri et al. [5], among others, advocate the restricted version. Searle [2], Hocking [8], and Searle et al. [4] advocate the unrestricted version. While the statistic for testing for no main random effect differs in the two versions, both use what Scheffé [1, page 264] calls the symmetry assumption. Basically, this is the assumption of independence of the random main and interaction effects. This assumption was criticized in [9, 10] as unrealistic in most practical situations. More importantly, simulations demonstrated that the classical F-test for testing the significance of the main fixed effect is very sensitive to departures from the symmetry assumption even in balanced designs with homoscedastic errors. Additional simulations showed the classical F-test for testing the significance of the interaction effect to also be very sensitive to unbalancedness and heteroscedasticity.This paper deals with testing the significance of the main random effect in two-factor mixed and random effects designs. Simulations suggest that the classicalF-test for this hypothesis is also very sensitive to departures from the symmetry assumption, as well as to unbalancedness and heteroscedasticity. New test procedures are presented under fully nonparametric models for the two-factor mixed and random effects designs. Unbalanced mixed and random effects designs are incorporated in the modeling by considering the sample sizes nij to be random. The models and test procedures pertain to both continuous and discrete data and allow heteroscedasticity in the error and interaction terms, as well as nonnormality. The interaction term is not assumed independent from the main effect (so the symmetry assumption is not made), though, under the random effects model, the two are shown to be uncorrelated.The asymptotic theory of the test statistics is derived under the Neyman-Scott framework, in the sense that the number of levels of both factors is large but the group sizes can remain fixed. In such cases, test statistics are constructed by taking the difference (as opposed to the ratio) of the appropriate mean squares; compare [9–12] and references therein. In accordance to this literature, the proposed test statistics for testing for no main random effects are scaled versions of MSB-MSE* and MSB-MSAB, for the mixed and random effects designs, respectively, where the precise definition of these mean squares is given in Section 3.The novelty of the proposed test statistics is twofold. First, the mean squares are based on unweighted cell averages to accommodate unbalanced designs. Second, in order to accommodate heteroscedasticity,MSE is replaced by a different linear combination of the cell sample variances, denoted by MSE*, which, under the null hypothesis, has the same expected value as MSB (ensuring zero mean value of the test statistic under the null hypothesis). A modified version of the MSE for accommodating heteroscedasticity was first used in Akritas and Papadatos [13] in the context of a high-dimensional one-way fixed effects design.The rest of the paper is organised as follows. Section2 reviews the fully nonparametric random and mixed effects models. In Section 3 we derive the test statistics and present their asymptotic distributions. Proofs of the propositions in Section 3 appear in the appendix. Section 4 discusses the estimation of the limiting null distributions and also presents results of simulations comparing our testing procedures to the usual F-tests. Finally, in Section 5 we present the analysis of a dataset from NOAA’s NS&T Mussel Watch Program.


## 2. Fully Nonparametric Models

Here we review the fully nonparametric random and mixed effects models developed in [9, 10].


### 2.1. The Random Effects Model

Consider a two-way random effects design. The random levels of the row factor are obtained by random sampling from the population𝒜S, while the random levels of the column factor are obtained by random sampling from the population 𝒜T. Let S denote a row level randomly selected from 𝒜S, and T denote a column level randomly selected from 𝒜T. Assume that the selections are independent. For each factor-level combination (S,T), observations YSTk,k=1,…,nST, are generated from a distribution function, FST. FST can be any arbitrary distribution function, which depends only on the factor levels S,T.DecomposeYSTk into a random mean and error term as follows:
(4)YSTk=E(YSTk∣S,T)+[YSTk-E(YSTk∣S,T)]=μST+ϵSTk,
where μST and ϵSTk are defined implicitly in (4). The definitions of μST and ϵSTk imply
(5)E(ϵSTk)=0,E(ϵSTk∣S)=0,E(ϵSTk∣T)=0,E(ϵSTk∣S,T)=0.
Conditioning on S,T and using the last relation in (5) it further follows that
(6)E(μSTϵSTk)=0,E(μSTϵSTk∣S)=0,E(μSTϵSTk∣T)=0.
Next, decompose the random means as
(7)μST=μ+αS+βT+γST,
where the overall mean and the random effects αS, βT, and γST are defined as
(8)μ=E(YSTk)=E(μST),αS=E(YSTk∣S)-μ=E(μST∣S)-μ,βT=E(YSTk∣T)-μ=E(μST∣T)-μ,γST=μST-μ-αS-βT.
Combining (4) and (7) we obtain
(9)YSTk=μ+αS+βT+γST+ϵSTk,
where the random effects have zero means, are both marginally and conditionally on S and T uncorrelated from the error terms (this follows from (6)), and satisfy
(10)E(αS∣T)=E(βT∣S)=E(γST∣S)=E(γST∣T)=0,E(αSβT)=E(αSγST)=E(βTγST)=0,
as is easily seen using the independence of S,T. Note that even though the main and interaction effects are uncorrelated, in general they will not be independent. This is because the interaction term is typically related to the two main effects by a multiplicative relation rendering the interaction effect nonnormal even in the presence of normal main effects. Therefore the symmetry assumption is typically violated in random effects designs.Relation (9) with the associated properties (5), (6), and (10) is derived using only the assumption that S and T are independent. Thus, it is a nonparametric version of the classical random effects model, in the sense that it allows for arbitrary (rather than normal) distributions of the random effects and error term, the variance (as well as the entire distribution) of the error term can depend on the random factor levels, and the effects are uncorrelated, and uncorrelated to the error term, as opposed to the assumed independence of the classical model.


### 2.2. The Mixed Effects Model

Consider a design where the rows correspond to the fixed effects and the columns correspond to the random effects. The levels of the random effects are obtained by simple random sampling from a population𝒜. Let T denote a randomly selected element from 𝒜, i a fixed effect level, and YiT an observation obtained from the cell (i,T). We place no restrictions on the form of the distribution function, FiT, of YiT other than the obvious restriction that it has to depend on i and T. In particular, FiT can be a discrete distribution function. Since T is obtained from a random selection, FiT, its mean μiT and variance σiT2 are all random. The nonparametric model consists of a decomposition of YiT in terms of μiT plus an error term, and further decomposition of μiT to define the main effects and interaction terms. In particular, write
(11)YiT=E(YiT∣T)+(YiT-E(YiT∣T))=μiT+(YiT-μiT)=μiT+ϵiT,
where the error term ϵiT is defined implicitly in the above relation. An easy consequence of the definition of ϵiT is
(12)E(ϵiT)=E(ϵiT∣T)=0,E(ϵiTμiT)=E(ϵiTμiT∣T)=0.
Next decompose μiT into main and interaction effects as
(13)μiT=μ+αi+βT+γiT,
where μ=(1/a)∑i‍E(μiT) is the overall mean and
(14)αi=E(μiT)-μ,βT=1a∑i‍μiT-μ,γiT=μiT-μ-αi-βT,
are the main and interaction effects. An immediate implication of (14) is
(15)∑i‍αi=0,E(βT)=0,∑i‍γiT=0,E(γiT)=0.
Using (12) it further follows that
(16)E(ϵiTβT)=E(ϵiTβT∣T)=0,E(ϵiTγiT)=E(ϵiTγiT∣T)=0.
Combining (11) and (13) we have
(17)YiT=μ+αi+βT+γiT+ϵiT.
While model (17) has the appearance of the usual mixed effects model, it does not require normality (or even continuity) of the observations, the variance of the random interaction effects can depend on i, and the random effects βT and γiT are not uncorrelated without further assumptions (which we do not make) about constant variance of the random means and constant covariance for any two distinct random means.Suppose now that there area fixed levels, so i=1,…,a, and let Tj,j=1,…,b, denote the levels of the random effect. Recall that the Tj are obtained by simple random sampling from 𝒜. Also let YiTjk,k=1,…,nij, denote the replications (iid observations) in factor-level combination (i,Tj). For notational simplicity we write Yijk for YiTjk and j for Tj. Thus, conditionally on j (or Tj), the Yijk are iid Fij. From (17) we have
(18)Yijk=μ+αi+βj+γij+ϵijk,
where the effects and error terms satisfy (14), (15), and (16). In addition, γij1 and γij2 are independent for j1≠j2, ϵijk, k=1,…,nij, are iid conditionally on j, ϵi1j1k1, and ϵi2j2k2 are independent for j1≠j2, and ϵi1jk1 and ϵi2jk2 are uncorrelated for i1≠i2. These imply
(19)E(ϵi1jk1ϵi2jk2)=0,fork1≠k2,E(ϵi1j1k1ϵi2j2k2)=0,forj1≠j2,E(ϵi1jk1ϵi2jk2)=0,fori1≠i2,E(γij1γij2)=0,forj1≠j2.
We note that relation (14) implies that γij is a function of the random level Tj. Since the random levels are iid, we have that the variances of γij, j=1,…,b, are equal. Arguing similarly, we have that the variances of the ϵijk, j=1,…,b, are equal. However the model allows heteroscedasticity across different levels of the fixed effect, that is, σ2(γi1j1)≠σ2(γi2j2) and σ2(ϵi1j1k1)≠σ2(ϵi2j2k2), for i1≠i2. Moreover, the conditional variances of the error term given different levels of the random effect can be different. Summarizing, we have
(20)σ2(γij1)=σ2(γij2),σ2(ϵij1k1)=σ2(ϵij2k2),
but the model allows
(21)σ2(γi1j1)≠σ2(γi2j2),fori1≠i2,σ2(ϵi1j1k1)≠σ2(ϵi2j2k2),fori1≠i2,Var(ϵij1k1∣Tj1)≠Var(ϵij2k2∣Tj2).


## 2.1. The Random Effects Model

Consider a two-way random effects design. The random levels of the row factor are obtained by random sampling from the population𝒜S, while the random levels of the column factor are obtained by random sampling from the population 𝒜T. Let S denote a row level randomly selected from 𝒜S, and T denote a column level randomly selected from 𝒜T. Assume that the selections are independent. For each factor-level combination (S,T), observations YSTk,k=1,…,nST, are generated from a distribution function, FST. FST can be any arbitrary distribution function, which depends only on the factor levels S,T.DecomposeYSTk into a random mean and error term as follows:
(4)YSTk=E(YSTk∣S,T)+[YSTk-E(YSTk∣S,T)]=μST+ϵSTk,
where μST and ϵSTk are defined implicitly in (4). The definitions of μST and ϵSTk imply
(5)E(ϵSTk)=0,E(ϵSTk∣S)=0,E(ϵSTk∣T)=0,E(ϵSTk∣S,T)=0.
Conditioning on S,T and using the last relation in (5) it further follows that
(6)E(μSTϵSTk)=0,E(μSTϵSTk∣S)=0,E(μSTϵSTk∣T)=0.
Next, decompose the random means as
(7)μST=μ+αS+βT+γST,
where the overall mean and the random effects αS, βT, and γST are defined as
(8)μ=E(YSTk)=E(μST),αS=E(YSTk∣S)-μ=E(μST∣S)-μ,βT=E(YSTk∣T)-μ=E(μST∣T)-μ,γST=μST-μ-αS-βT.
Combining (4) and (7) we obtain
(9)YSTk=μ+αS+βT+γST+ϵSTk,
where the random effects have zero means, are both marginally and conditionally on S and T uncorrelated from the error terms (this follows from (6)), and satisfy
(10)E(αS∣T)=E(βT∣S)=E(γST∣S)=E(γST∣T)=0,E(αSβT)=E(αSγST)=E(βTγST)=0,
as is easily seen using the independence of S,T. Note that even though the main and interaction effects are uncorrelated, in general they will not be independent. This is because the interaction term is typically related to the two main effects by a multiplicative relation rendering the interaction effect nonnormal even in the presence of normal main effects. Therefore the symmetry assumption is typically violated in random effects designs.Relation (9) with the associated properties (5), (6), and (10) is derived using only the assumption that S and T are independent. Thus, it is a nonparametric version of the classical random effects model, in the sense that it allows for arbitrary (rather than normal) distributions of the random effects and error term, the variance (as well as the entire distribution) of the error term can depend on the random factor levels, and the effects are uncorrelated, and uncorrelated to the error term, as opposed to the assumed independence of the classical model.


## 2.2. The Mixed Effects Model

Consider a design where the rows correspond to the fixed effects and the columns correspond to the random effects. The levels of the random effects are obtained by simple random sampling from a population𝒜. Let T denote a randomly selected element from 𝒜, i a fixed effect level, and YiT an observation obtained from the cell (i,T). We place no restrictions on the form of the distribution function, FiT, of YiT other than the obvious restriction that it has to depend on i and T. In particular, FiT can be a discrete distribution function. Since T is obtained from a random selection, FiT, its mean μiT and variance σiT2 are all random. The nonparametric model consists of a decomposition of YiT in terms of μiT plus an error term, and further decomposition of μiT to define the main effects and interaction terms. In particular, write
(11)YiT=E(YiT∣T)+(YiT-E(YiT∣T))=μiT+(YiT-μiT)=μiT+ϵiT,
where the error term ϵiT is defined implicitly in the above relation. An easy consequence of the definition of ϵiT is
(12)E(ϵiT)=E(ϵiT∣T)=0,E(ϵiTμiT)=E(ϵiTμiT∣T)=0.
Next decompose μiT into main and interaction effects as
(13)μiT=μ+αi+βT+γiT,
where μ=(1/a)∑i‍E(μiT) is the overall mean and
(14)αi=E(μiT)-μ,βT=1a∑i‍μiT-μ,γiT=μiT-μ-αi-βT,
are the main and interaction effects. An immediate implication of (14) is
(15)∑i‍αi=0,E(βT)=0,∑i‍γiT=0,E(γiT)=0.
Using (12) it further follows that
(16)E(ϵiTβT)=E(ϵiTβT∣T)=0,E(ϵiTγiT)=E(ϵiTγiT∣T)=0.
Combining (11) and (13) we have
(17)YiT=μ+αi+βT+γiT+ϵiT.
While model (17) has the appearance of the usual mixed effects model, it does not require normality (or even continuity) of the observations, the variance of the random interaction effects can depend on i, and the random effects βT and γiT are not uncorrelated without further assumptions (which we do not make) about constant variance of the random means and constant covariance for any two distinct random means.Suppose now that there area fixed levels, so i=1,…,a, and let Tj,j=1,…,b, denote the levels of the random effect. Recall that the Tj are obtained by simple random sampling from 𝒜. Also let YiTjk,k=1,…,nij, denote the replications (iid observations) in factor-level combination (i,Tj). For notational simplicity we write Yijk for YiTjk and j for Tj. Thus, conditionally on j (or Tj), the Yijk are iid Fij. From (17) we have
(18)Yijk=μ+αi+βj+γij+ϵijk,
where the effects and error terms satisfy (14), (15), and (16). In addition, γij1 and γij2 are independent for j1≠j2, ϵijk, k=1,…,nij, are iid conditionally on j, ϵi1j1k1, and ϵi2j2k2 are independent for j1≠j2, and ϵi1jk1 and ϵi2jk2 are uncorrelated for i1≠i2. These imply
(19)E(ϵi1jk1ϵi2jk2)=0,fork1≠k2,E(ϵi1j1k1ϵi2j2k2)=0,forj1≠j2,E(ϵi1jk1ϵi2jk2)=0,fori1≠i2,E(γij1γij2)=0,forj1≠j2.
We note that relation (14) implies that γij is a function of the random level Tj. Since the random levels are iid, we have that the variances of γij, j=1,…,b, are equal. Arguing similarly, we have that the variances of the ϵijk, j=1,…,b, are equal. However the model allows heteroscedasticity across different levels of the fixed effect, that is, σ2(γi1j1)≠σ2(γi2j2) and σ2(ϵi1j1k1)≠σ2(ϵi2j2k2), for i1≠i2. Moreover, the conditional variances of the error term given different levels of the random effect can be different. Summarizing, we have
(20)σ2(γij1)=σ2(γij2),σ2(ϵij1k1)=σ2(ϵij2k2),
but the model allows
(21)σ2(γi1j1)≠σ2(γi2j2),fori1≠i2,σ2(ϵi1j1k1)≠σ2(ϵi2j2k2),fori1≠i2,Var(ϵij1k1∣Tj1)≠Var(ϵij2k2∣Tj2).


## 3. The Test Statistics and Their Asymptotic Distributions

Before proceeding, we first define our notation:(22)Y¯ij·=1nij∑k=1nij‍Yijk,Y~i··=1b∑j=1b‍Y¯ij·,Y~·j·=1a∑i=1a‍Y¯ij·,Y~···=1a∑i=1a‍Y~i··=1b∑j=1b‍Y~·j·,β¯·=1b∑j=1b‍βj,γ¯i·=1b∑j=1b‍γij,ϵ¯ij·=1nij∑k=1nij‍ϵijk,ϵ~i··=1b∑j=1b‍ϵ¯ij·,ϵ~·j·=1a∑i=1a‍ϵ¯ij·,ϵ~···=1b∑j=1b‍ϵ~·j·=1a∑i=1a‍ϵ~i··.
The usual F statistic in the balanced case in the restricted version of the mixed model is given by
(23)F=MSBMSE=a2bn(n-1)∑j(Y¯·j·-Y¯···)2(b-1)∑i∑j∑k(Yijk-Y¯ij·)2,
while the usual F statistic in the balanced case in the random effects model is given by
(24)F=MSBMSAB=a(a-1)∑j(Y¯·j·-Y¯···)2∑i∑j(Y¯ij·-Y¯i··-Y¯·j·+Y¯···)2.
In the unbalanced case there is no theoretically justified procedure, but Searle [14] does have a number of recommendations which are implemented in software packages.The test statistics we propose for testing the hypothesis of no main random effect,H0(B):σ2(βj)=0, j=1,…,b, are modeled after the usual F statistics given in (23) and (24), except that we take the difference rather that the ratio, as explained in the Introduction. Moreover, in order to accommodate unbalanced designs, Y¯i.., Y¯.j., and Y¯… are replaced by Y~i··, Y~·j·, and Y~···, respectively. As mentioned in the Introduction, our asymptotic theory requires that the difference has mean value zero. While this is the case for the random effects design, the expected value of MSB-MSE is not zero for the mixed effects design with heteroscedastic errors. Therefore we replace MSE with a different linear combination of the cell sample variances:
(25)MSE*=1a3b∑i=1a‍∑j=1bn·jnij1nij-1∑k=1nij‍(Yijk-Y¯ij·)2.Thus, the proposed test statistics for the mixed effects and random effects models are, respectively,(26)SB,mixed=1a(b-1)∑i=1a‍∑j=1b‍∑k=1nij‍(Y~·j·-Y~···)2-MSE*,(27)SB,rand=1a2b∑i=1a‍∑j=1b‍[(a-1)(Y~·j·-Y~···)2-Rij2],
where Rij=Y¯ij·-Y~·j·-Y~i··+Y~···.Proposition 1.
LetSB,mixed be as defined in (26) and let SB,rand be as defined in (27), with the observations Yijk generated according to the model of Section 2. Then, (i)
if the null hypothesisH0(B):σ2(βj)=0, j=1,…,b holds,
(28)E(SB,mixed)=E(SB,rand)=0,(ii)
if the null hypothesisH0(B):σ2(βj)=0, j=1,…,b does not hold,
(29)E(SB,mixed)>0,E(SB,rand)>0.The second item in Proposition1 suggests that the null hypothesis should be rejected at level α whenever the test statistic takes a value larger than the 100(1-α)th percentile of the null distribution. The next proposition establishes asymptotic representations of SB,mixed and SB,rand that are useful for establishing their asymptotic distributions.Proposition 2.
UnderH0(B), (i)
(30)bSB,mixed=1ab∑j‍n·j[(ϵ~·j·-ϵ~···)2-E(ϵ~·j·-ϵ~···)2]+op(1),(ii)
(31)aSB,rand=1ab∑i‍∑i1≠i‍∑j‍γijγi1j+op(1).Theorem 3.
Letsb2=Var((1/b)∑j‍nij[(ϵ~·j·-ϵ~···)2-E(ϵ~·j·-ϵ~···)2]). Then under H0(B)(32)bSB,mixedsb→N(0,1),
in distribution as a,b→∞.Proof.
LettingTb=(1/b)∑j‍nij[(ϵ~·j·-ϵ~···)2-E(ϵ~·j·-ϵ~···)2], then, as stated above,
(33)bSB,mixedsb≈Tbsb.
The proof that Tb/sb→N(0,1) in distribution follows from an application of the Lindeberg-Feller CLT. We first rewrite our statement in the notation of the Lindeberg-Feller CLT and then show that the Lindeberg condition holds after a straightforward application of the Dominated Convergence Theorem. Now let Yja,b=(1/b)nij[(ϵ~·j·-ϵ~···)2-E(ϵ~·j·-ϵ~···)2], so that Tb=∑j‍Yja,b. Also note that we can write
(34)sb2=Var(Tb)=1b∑j‍nij2Var((ϵ~·j·-ϵ~···)2),
which is finite because we have assumed a finite fourth moment for the ϵ terms. Moreover, this implies that
(35)Var(bYja,b)=nij2Var((ϵ~·j·-ϵ~···)2),
which is again finite. Now in order to use the Lindeberg-Feller CLT, we need to verify the Lindeberg condition:
(36)∀ϵ>0,1sb2∑j=1b‍E[(Yja,b)2I{|Yja,b|>ϵsb}]→0.
Note that
(37)1sb2∑j=1b‍E[(Yja,b)2I{|Yja,b|>ϵsb}]=1bsb2∑j=1b‍E[(bYja,b)2I{|bYja,b|>bϵsb}]=1sb2E[(bY1a,b)2I{|bY1a,b|>bϵsb}].
Because E(bY1a,b)=0 and Var(bY1a,b) is finite, the Dominated Convergence Theorem implies that E[(bY1a,b)2I{|bY1a,b|>bϵsb}]→0, and hence the Lindeberg condition is verified. This completes the proof of the theorem.Theorem 4.
Define the operatorHb:g→Hbg, for g∈L2([0,1],ℬ[0,1],λ[0,1]), where ℬ[0,1] is the class of Borel subsets of [0,1] and λ[0,1] is the Lebesgue measure on ℬ[0,1], by
(38)(Hbg)(s)=∫‍hb(s,t)g(t)dt,
where hb(s,t)=(1/b)∑j=1b‍γsjγtj. Let λνb, fνb, ν=1,2,…, be the eigenvalues and eigenfunctions of the integral equation
(39)(Hbf)(s)=λf(s),
and assume that the eigenvalues satisfy
(40)∑ν=1∞‍(supbλνb)2<∞.
Then, if γij and ϵij have finite fourth moments and the covariances Cov(γij,γij′), Cov(γij2,ϵ-ij.2), Cov(γij2,γij′2), and Cov(ϵij2,ϵij′2) are bounded uniformly in j,j′, the asymptotic null distribution of SB,rand, as a→∞, b→∞, is the same as the asymptotic distribution of ∑ν=1∞λνb(Zν2-1), as b→∞, where Zν, ν=1,2,…, are independent standard normal random variables.Proof.
The proof proceeds by first conditioning on the levels of the column factorT=(T1,…,Tb) and then applying the arguments used in the proof of Theorem 3.3 of [9]. Indeed, the proof of the aforementioned theorem uses an expansion similar to (31) of Proposition 2 (with the indices i,j interchanged) which is modeled as a U-statistic with kernel depending on a. The same arguments establish the limiting distribution of the test statistic conditionally on T=(T1,…,Tb). As b→∞, the strong law of large numbers guarantees that the kernel hb(s,t) converges to the same limit, and thus the conditional limiting distribution does not depend on T=(T1,…,Tb). Thus, the test statistic converges to the same distribution unconditionally.


## 4. Test Procedures and Simulation Results

In testing for the random main effect in the mixed effects model, we only need to estimate the variance termsb2 in Theorem 3. According to (30) of Proposition 2 we need to estimate the variance of
(41)1b∑j‍nij(ϵ~·j·-ϵ~···)2.
It is easy to see that
(42)1b∑j‍nij(ϵ~·j·-ϵ~···)2-1b∑j‍nijϵ~·j·2→0,
in probability. Thus, we need to estimate the variance of the independent terms ϵ~·j·2. By the Central Limit Theorem, as a→∞, ϵ~·j· is approximately Normally distributed with mean zero and variance (1/a2)∑i‍(σij2/nij). Therefore, the approximate variance of ϵ~·j·2 is 2((1/a2)∑i‍(σij2/nij))2. As a result, we estimate sb2 by
(43)s^b2=2a6b∑j‍(n·j∑i‍sij2nij)2.
Then, according to Proposition 1 and Theorem 3, H0(B) is rejected at level α when
(44)b1/2SB,mixeds^b>zα,
where zα is the (1-α)100th percentile of the standard normal distribution.Testing for main effects in the random effects model involves estimation of the asymptotic distribution of Theorem4. Koltchinskii and Giné [15] show that the infinite spectrum of a Hilbert-Schmidt operator H can be approximated by the finite spectrum of the empirical matrix version of the operator. We will use this result to approximate the spectrum of the operator Hb in Theorem 4 by its empirical version, H~a, defined below. To this end, the following representation of SB,rand as a U-statistic is useful:
(45)aSB,rand=1ab∑i1‍∑i2≠i1‍∑j‍(Y¯i1j·-Y¯i1··)(Y¯i2j·-Y¯i2··)=a-1bUa,
where the U-statistic Ua is defined by the kernel
(46)h(Si1,Si2)=∑j‍(Y¯i1j·-Y¯i1··)(Y¯i2j·-Y¯i2··).
Under the null hypothesis H0(B), the empirical version, H~a, of the operator Hb is the a×a matrix with entries H~i,i′=(1/a)h(Si,Si′). The a eigenvalues λ1a,…,λaa of H~a are then used to estimate the asymptotic distribution of SB,rand as ∑i=1aλia(Zi2-1), where the Zi2 are randomly sampled from the χ12 distribution. To approximate the estimated distribution, we randomly generate a large number B of sets Z12,…,Za2, and use the formula above to obtain B realizations from the distribution. Using this approximate distribution, for a level α test we reject H0(B) when SB,rand is greater than the (1-α)100th percentile.We now present simulation results to compare our test procedures for the main random effect in the mixed and random effects models to the usualF-tests given by (23) and (24). Note that these statistics are only valid for balanced designs. In the unbalanced designs, the datasets used to obtain the achieved size and power are written out to files and input to SAS for the F-test to be carried out in SAS PROC MIXED. The achieved size and power are calculated from the P values given from SAS.Using the decomposition in relation (18), we generate our data for the mixed effects model as follows.(1)
Setμ equal to any constant.(2)
Generate a vector of main random effectsTj, j=1,…,b from the standard normal distribution.(3)
Independently generate a vector of fixed “functions”λi, i=1,…,a, from the standard exponential distribution and set γi*=λi-λ¯.(4)
GenerateYijk from a normal distribution with variance σi2 and mean μ+αi+Tj+Tj*γi*, where the last two terms are the main random and interaction effects, respectively.Similarly, using the decomposition in relation (9), we generate our data for the random effects model as follows.(1)
Setμ equal to any constant.(2)
Generate a vector of main row random effectsSi, i=1,…,a from the standard normal distribution.(3)
Independently generate a vector of main column random effectsTj, j=1,…,b from the standard normal distribution, and set Tj*=τTj, where τ is used to set the degree of deviation from the null hypothesis.(4)
Set the interaction effectsγij=SiTj.(5)
GenerateYijk from a normal distribution with variance σi2 and mean μ+Si+Tj*+SiTj.In Table1, each simulation used a=b=20. For the unbalanced setups, rows i=1,…,15 had 4 observations per cell, whereas rows i=16,…,20 had 2 observations per cell. For the heteroscedastic setups, the first 15 rows all had variance 1, while the remaining rows had variance 5; that is, σ12=⋯=σ152=1,σ162=⋯=σ202=5.Table 1
Comparison of achieved type I error rates and power forSB,mixed, SB,rand, and F with N=1000 tests for the random and mixed effects models.
Random effects model
Mixed effects model
Balanced
Unbalanced
Balanced
Unbalanced
α
=
.
05
α
=
.
10
α
=
.
05
α
=
.
10
α
=
.
05
α
=
.
10
α
=
.
05
α
=
.
10
S
B
F
S
B
F
S
B
F
S
B
F
S
B
F
S
B
F
S
B
F
S
B
F
Homo.
H0
.054
.220
.111
.247
.049
.167
.092
.198
.065
.047
.112
.099
.063
.045
.109
.086
A1
.053
.229
.123
.265
.069
NA
.140
NA
.227
.187
.318
.296
.184
.138
.260
.228
A2
.125
.337
.221
.367
.113
NA
.201
NA
.761
.719
.820
.804
.542
.489
.634
.638
A3
.236
.532
.369
.563
.234
NA
.360
NA
.975
.963
.982
.982
.888
.863
.920
.917
Hetero.
H0
.045
.150
.089
.182
.015
.230
.042
.303
.070
.056
.106
.100
.062
.935
.107
.952
A1
.045
.160
.088
.199
.021
NA
.053
NA
.077
.060
.129
.121
.074
NA
.124
NA
A2
.114
.273
.190
.316
.047
NA
.094
NA
.151
.129
.221
.215
.104
NA
.152
NA
A3
.195
.430
.318
.467
.117
NA
.215
NA
.270
.235
.353
.341
.131
NA
.189
NA


## 5.  Data Analysis

One real-world application for our methodology can be found through the National Oceanic and Atmospheric Administration’s National Status and Trends Program. In 1986, this division undertook a very large scale project to monitor the levels of numerous chemical contaminants and organic chemical constituents in marine sediment and bivalve (mollusk) tissue samples. This project (http://ccma.nos.noaa.gov/about/coast/nsandt/musselwatch.aspx), dubbed the Mussel Watch Project, is still on-going and there are no apparent plans to discontinue it in the near future. There are currently over 300 coastal sites at which sediment and bivalve samples are collected and analysed for the project. Each site is categorised as being within a certain estuarine drainage area (EDA). For our data analysis, we chose to analyse the Mercury concentrations in tissue samples of the Crassostrea virginica, or Eastern Oyster. There was a high degree of missingness of Mercury concentrations in the dataset, so we used imputations to generate these measurements. Because of how we performed our imputations, only EDAs that did not have missing measurements for two consecutive years and EDAs that had measurements in 1986 were included in the analysis. Once we had a complete imputed dataset, we employed our methodology to test for effects due to EDA.


### 5.1. Imputations

Because of the1/(nij-1) coefficient that is used in the statistic SB,mixed, each cell needs to have at least 2 observations. To facilitate the imputations, we identified 195 instances of 2 consecutive measurements falling into the following groupings: 1986-1987, 1988-1989, 1990-1991, 1992-1993, 1994-1995, 1996-1997, 1998-1999, 2000-2001, 2002-2003, and 2004-2005. If any cell had more than one observation, they were averaged. We used these 195 pairs to build the regression model CV^t,j=0.046424+0.65852CVt-1,j, where we predict a value of Mercury concentration at time t from the Mercury concentration at time t-1. For this model, the P value for the test of H0:β1=0 versus HA:β1≠0 was 0.000 and the R2=36.8%. When fitting the regression model in Minitab, we also saved the residuals from the model and the actual imputations that we used were CV^t,j+R*, where R* is a randomly selected residual from the fitting of the model. We separately treat imputations for cells with only one observation and imputations for empty cells. When cells are empty, we use the regression equation
(47)CV^t,j=0.046424+0.65852CVt-1,j+R*,
twice to impute two values of Hg concentration at time t. When a cell has only one observation, we use a different technique.(1)
If the EDA has any cells with more than one observation in any year, we calculate residuals for each year with more than one observation aseijk=Yijk-Y¯ij·,k=1,…,nij. Then for the EDA and year of interest, we impute a second observation as Yij2=Yij1+eijk*, where eijk* denotes a randomly selected residual.(2)
If the EDA of interest has at most one observation in each cell, we identify the closest neighboring EDA (denoted by levelj*) and calculate residuals for each year with more than one observation in the neighboring EDA as eij*k=Yij*k-Y¯ij*·,k=1,…,nij*. Then for the EDA and year of interest, we impute a second observation as Yij2=Yij1+eij*k*, where eij*k* denotes a randomly selected residual.After performing the imputations, there are a total ofN=1808 observations of Hg concentration in CV tissue samples. The 20 years included are 1986–2005, and there are 39 EDAs included in the dataset. The data are rather unbalanced, with most EDAs having only two observations per year, while others such as the Galveston Bay EDA have as many as six observations per year. Moreover, the observations do not seem to be homoscedastic. These points are illustrated in Figure 1, where we show scatterplots of Hg concentrations for six EDAs after imputation.Scatterplot of Mercury concentrations in Eastern Oysters by year after imputation for six EDAs.
(a)
(b)
(c)
(d)
(e)
(f)


### 5.2. The Analysis

In our analysis of this dataset under the mixed effects model, we take the years as the fixed effect and the EDAs as the random effect. We can think of the EDAs as a random effect because we are only analyzing a very small subset of a much larger group of EDAs and the location effect is not of specific interest. Some may view the time effect as being of specific interest, while others may not. Whether one views the time effect as a fixed or random factor would precipitate analysis using a mixed or random effects model.We want to model the Hg concentration with an interaction effect between year and location, so we use the model in relation (18). Because we are dealing with compositional data, we employ the common practice of analyzing the data after taking a log transformation; compare [16].When we apply our proposed test procedure based on the test statisticSC, we obtain a P value of 0.0588. When we apply the F-test computed from SAS PROC MIXED, we obtain a P value of 0.0016.TheP value for the test using SC is simply calculated as 1-Φ(·), where Φ is the distribution function of the standard normal random variable. The P value for the SAS F-test is calculated as 1-Fν1,ν2(·), where ν1 and ν2 denote the numerator and denominator degrees of freedom. In the test for the random interaction effect, SAS uses ν1=722, ν2=1028. The simulation results shown in Table 1 suggest that the significance of the interaction effect indicated by the F-test cannot be trusted. Indeed, according to these simulation results, we expect the F-test to detect interaction in almost any unbalanced and heteroscedastic dataset. On the other hand, the simulation results suggest that the marginal significance indicated by the P value obtained from the proposed test procedure can be trusted. According to this we recommend that the interaction effect be included in the modeling of this dataset.


## 5.1. Imputations

Because of the1/(nij-1) coefficient that is used in the statistic SB,mixed, each cell needs to have at least 2 observations. To facilitate the imputations, we identified 195 instances of 2 consecutive measurements falling into the following groupings: 1986-1987, 1988-1989, 1990-1991, 1992-1993, 1994-1995, 1996-1997, 1998-1999, 2000-2001, 2002-2003, and 2004-2005. If any cell had more than one observation, they were averaged. We used these 195 pairs to build the regression model CV^t,j=0.046424+0.65852CVt-1,j, where we predict a value of Mercury concentration at time t from the Mercury concentration at time t-1. For this model, the P value for the test of H0:β1=0 versus HA:β1≠0 was 0.000 and the R2=36.8%. When fitting the regression model in Minitab, we also saved the residuals from the model and the actual imputations that we used were CV^t,j+R*, where R* is a randomly selected residual from the fitting of the model. We separately treat imputations for cells with only one observation and imputations for empty cells. When cells are empty, we use the regression equation
(47)CV^t,j=0.046424+0.65852CVt-1,j+R*,
twice to impute two values of Hg concentration at time t. When a cell has only one observation, we use a different technique.(1)
If the EDA has any cells with more than one observation in any year, we calculate residuals for each year with more than one observation aseijk=Yijk-Y¯ij·,k=1,…,nij. Then for the EDA and year of interest, we impute a second observation as Yij2=Yij1+eijk*, where eijk* denotes a randomly selected residual.(2)
If the EDA of interest has at most one observation in each cell, we identify the closest neighboring EDA (denoted by levelj*) and calculate residuals for each year with more than one observation in the neighboring EDA as eij*k=Yij*k-Y¯ij*·,k=1,…,nij*. Then for the EDA and year of interest, we impute a second observation as Yij2=Yij1+eij*k*, where eij*k* denotes a randomly selected residual.After performing the imputations, there are a total ofN=1808 observations of Hg concentration in CV tissue samples. The 20 years included are 1986–2005, and there are 39 EDAs included in the dataset. The data are rather unbalanced, with most EDAs having only two observations per year, while others such as the Galveston Bay EDA have as many as six observations per year. Moreover, the observations do not seem to be homoscedastic. These points are illustrated in Figure 1, where we show scatterplots of Hg concentrations for six EDAs after imputation.Scatterplot of Mercury concentrations in Eastern Oysters by year after imputation for six EDAs.
(a)
(b)
(c)
(d)
(e)
(f)


## 5.2. The Analysis

In our analysis of this dataset under the mixed effects model, we take the years as the fixed effect and the EDAs as the random effect. We can think of the EDAs as a random effect because we are only analyzing a very small subset of a much larger group of EDAs and the location effect is not of specific interest. Some may view the time effect as being of specific interest, while others may not. Whether one views the time effect as a fixed or random factor would precipitate analysis using a mixed or random effects model.We want to model the Hg concentration with an interaction effect between year and location, so we use the model in relation (18). Because we are dealing with compositional data, we employ the common practice of analyzing the data after taking a log transformation; compare [16].When we apply our proposed test procedure based on the test statisticSC, we obtain a P value of 0.0588. When we apply the F-test computed from SAS PROC MIXED, we obtain a P value of 0.0016.TheP value for the test using SC is simply calculated as 1-Φ(·), where Φ is the distribution function of the standard normal random variable. The P value for the SAS F-test is calculated as 1-Fν1,ν2(·), where ν1 and ν2 denote the numerator and denominator degrees of freedom. In the test for the random interaction effect, SAS uses ν1=722, ν2=1028. The simulation results shown in Table 1 suggest that the significance of the interaction effect indicated by the F-test cannot be trusted. Indeed, according to these simulation results, we expect the F-test to detect interaction in almost any unbalanced and heteroscedastic dataset. On the other hand, the simulation results suggest that the marginal significance indicated by the P value obtained from the proposed test procedure can be trusted. According to this we recommend that the interaction effect be included in the modeling of this dataset.


## 6. Summary

The classicalF statistic for testing the significance of main random effects in two-factor mixed and random effects model is shown to be very sensitive to violations of the assumptions under which it is derived, in particular those of symmetry, homoscedasticity, and balancedness. Two new test procedures for testing the hypothesis of no main random effects are proposed under fully nonparametric modeling of the mixed and random effects designs. The test statistics are defined as differences (as opposed to ratios) of suitably defined mean squares, and their asymptotic theory is derived as the number of levels tends to infinity. Simulations suggest that the new procedures perform reasonably well in situations where the F statistic breaks down and have comparable performance to it when the assumptions it requires are met. The practical value of the proposed procedures is demonstrated with the analysis of a real dataset.



---

*Source: 708540-2013-03-21.xml*