# Testing for main random effects in two-way random and mixed effects models: modifying the F statistic.

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 the F-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.1. Introduction

Consider a two-factor mixed or random effects design, and let [Y.sub.ijk] denote the fcth 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

[Y.sub.ijk] = [mu] + [[alpha].sub.i] + [beta[].sub.j] + [[gamma].sub.ij] + [[member of].sub.ijk]. (1)

For the random effects model, the [[alpha].sub.i]'s, [[beta].sub.j]'s, [[gamma].sub.ij]'s, and [[member of].sub.ijk]'s are mutually independent, the [[alpha].sub.i]'s are iid N(0,[[sigma].sup.2.sub.[alpha]]), the [[beta].sub.j]'s are iid N(0, [[sigma].sup.2.sub.[beta]]),the [[gamma].sub.ij]'s are iid N(0, [[sigma].sup.2.sub.[gamma]]), and the [[member of].sub.ijk]'s are iid N(0, [[sigma].sup.2.sub.[member of]]).

For the mixed effects model, there are two common definitions of the effects. The unrestricted version of the model assumes

[summation over (i)] [[alpha].sub.i] = 0, [[beta].sub.j] are iid N (0, [[sigma].sup.2.sub.[beta]], [[gamma].sub.ij] are iid N(0, [[sigma].sup.2.sub.[gamma]], (2)

and the [[beta].sub.j], [[gamma].sub.ij] are all independent. The restricted version of the model keeps the above assumptions but requires that

[summation over (i)] [[gamma].sub.ij] = 0, (3)

implying that [[gamma].sub.ij], [[gamma].sub.i'j] are correlated. Both versions assume the [[member of].sub.ijk] to be iid N(0, [[sigma].sup.2.sub.[epsilon]]) 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], Scheffe [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 Scheffe [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 classical F-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 [n.sub.ij] 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.sup.*] 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.sup.*], 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. Section 2 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 [A.sub.s], while the random levels of the column factor are obtained by random sampling from the population [A.sub.T]. Let S denote a row level randomly selected from [A.sub.s], and T denote a column level randomly selected from [A.sub.T]. Assume that the selections are independent. For each factor-level combination (S, T), observations [Y.sub.STk], k = 1, ..., [n.sub.ST], are generated from a distribution function, [F.sub.ST]. [F.sub.ST] can be any arbitrary distribution function, which depends only on the factor levels S, T.

Decompose [Y.sub.STk] into a random mean and error term as follows:

[Y.sub.STk] = E ([Y.sub.STk]|S, T) + [[Y.sub.STk] - E [([Y.sub.STk] | S, T)] = [[mu].sub.ST] + [[member of].sub.STk], (4)

where [[mu].sub.ST] and [[member of].sub.STk] are defined implicitly in (4). The definitions of [[mu].sub.ST] and [[member of].sub.STk] imply

E([[member of].sub.STk]) = 0, E([[member of].sub.STk] | S) = 0, E([[member of].sub.STk] | T) = 0, E([[member of].sub.STk] | S, T) = 0. (5)

Conditioning on S, T and using the last relation in (5) it further follows that

E{[[mu].sub.ST] [[member of].sub.STk]) = 0, E {[[[mu].sub.ST] [[member of].sub.STk] | S) = 0, E([[mu].sub.sT] [[member of].sub.STk] | T) = 0. (6)

Next, decompose the random means as

[[mu].sub.ST] = [mu] + [[alpha].sub.s] + [[beta].sub.T] + [[gamma].sub.ST], (7)

where the overall mean and the random effects [[alpha].sub.s], [[beta].sub.T], and [[gamma].sub.ST] are defined as

[mu] = E ([Y.sub.STk]) = E ([[mu].sub.ST]), [[alpha].sub.s] = E ([Y.sub.STk]|S) - [mu] = E ([[mu].sub.ST] |S) - [mu], [[beta].sub.T] = E([Y.sub.STk] | T) - [mu] = E([[mu].sub.ST] |T)- [mu], [[gamma].sub.ST] = [[mu].sub.ST] - [mu] - [[alpha].sub.s] - [[beta].sub.T]. (8)

Combining (4) and (7) we obtain

[Y.sub.STk] = [mu] + [[alpha].sub.s] + [[beta].sub.T] + [[gamma].sub.ST] + [[member of].sub.STk], (9)

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

E([[alpha].sub.s]|T) = E([[beta].sub.T] |S) = E([[gamma].sub.ST] |S) = E([[gamma].sub.ST] |T) = 0, E([[alpha].sub.s] [[beta].sub.T]) = E {[[alpha].sub.s][[gamma].sub.ST]) = E ([[beta].sub.T] [[gamma].sub.ST]) = 0, (10)

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 A. Let T denote a randomly selected element from A, i a fixed effect level, and [Y.sub.iT] an observation obtained from the cell (i, T). We place no restrictions on the form of the distribution function, [F.sub.iT], of [Y.sub.iT] other than the obvious restriction that it has to depend on i and T. In particular, [F.sub.iT] can be a discrete distribution function. Since T is obtained from a random selection, [F.sub.iT], its mean [[mu].sub.iT] and variance [[sigma].sup.2.sub.iT] are all random. The nonparametric model consists of a decomposition of [Y.sub.iT] in terms of [[mu].sub.iT] plus an error term, and further decomposition of [[mu].sub.iT] to define the main effects and interaction terms. In particular, write

[Y.sub.iT] = E([Y.sub.iT] |T) + ([Y.sub.iT] - E([Y.sub.iT] |T)) = [[mu].sub.iT] + ([Y.sub.iT] - [[mu].sub.iT]) = [[mu].sub.iT] + [[member of].sub.iT], (11)

where the error term [[member of].sub.iT] is defined implicitly in the above relation. An easy consequence of the definition of [[member of].sub.iT] is

E([[member of].sub.iT]) = E ([[member of].sub.iT]|T) = 0, E ([[member of].sub.iT][[mu].sub.iT]) = E ([[member of].sub.iT][[mu].sub.iT]|T) = 0. (12)

Next decompose [[mu].sub.iT] into main and interaction effects as

[[mu].sub.iT] = [mu] + [[alpha].sub.i] + [[beta].sub.T] + [[gamma].sub.iT], (13)

where [mu] = (1/a) [[SIGMA].sub.i]E(([[mu].sub.iT]) is the overall mean and

[[alpha].sub.i] = E([[mu].sub.iT])-[mu], [[beta].sub.T] = 1/a [summation over (i)] [[mu].sub.iT]-[mu], [[gamma].sub.iT] = [[mu].sub.iT]-[mu]-[[alpha].sub.i]-[[beta].sub.T], (14)

are the main and interaction effects. An immediate implication of (14) is

[summation over (i)][[alpha].sub.i] = 0, E([[beta].sub.T]) = 0, [summation over (i)] [[gamma].sub.iT] = 0, E([[gamma].sub.iT]) = 0. (15)

Using (12) it further follows that

E([[member of].sub.iT][beta]T) = E([[member of].sub.iT][[beta].sub.T]|T) = 0, E([[member of].sub.iT][[gamma].sub.iT]) = E([[member of].sub.iT][[gamma].sub.iT]|T) = 0.

Combining (11) and (13) we have

[Y.sub.iT] = [mu] + [[alpha].sub.i] + [[beta].sub.T] + [[gamma].sub.iT] + [[member of].sub.iT]. (17)

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 [[beta].sub.T] and [[gamma].sub.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 are a fixed levels, so i = 1, ..., a, and let [T.sub.j], j = 1, ..., b, denote the levels of the random effect. Recall that the [T.sub.j] are obtained by simple random sampling from A. Also let [Y.sub.iTjk], k = l, ..., [n.sub.ij], denote the replications (iid observations) in factor-level combination (i, [T.sub.j]). For notational simplicity we write [Y.sub.ijk] for [Y.sub.iTjk] and j for [T.sub.j]. Thus, conditionally on j (or [T.sub.j]), the [Y.sub.ijk] are iid [F.sub.ij]. From (17) we have

[Y.sub.ijk] = [mu] + [[alpha].sub.i] + [[beta].sub.j] + [[gamma].sub.ij] + [[member of].sub.ijk], (18)

where the effects and error terms satisfy (14), (15), and (16). In addition, [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] are independent for [j.sub.1] [not equal to] [j.sub.2], [[member of].sub.ijk], k = 1, ..., [n.sub.ij], are iid conditionally on j, [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] are independent for [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] are uncorrelated for [i.sub.1] [not equal to] [i.sub.2]. These imply

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (19)

We note that relation (14) implies that [[gamma].sub.ij] is a function of the random level [T.sub.j]. Since the random levels are iid, we have that the variances of [[gamma].sub.ij], j = 1, ..., b, are equal. Arguing similarly, we have that the variances of the [[member of].sub.ijk], j = 1, ..., b, are equal. However the model allows heteroscedasticity across different levels of the fixed effect, that is, [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. Moreover, the conditional variances of the error term given different levels of the random effect can be different. Summarizing, we have

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (20)

but the model allows

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (21)

3. The Test Statistics and Their Asymptotic Distributions

Before proceeding, we first define our notation:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (22)

The usual F statistic in the balanced case in the restricted version of the mixed model is given by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (23)

while the usual F statistic in the balanced case in the random effects model is given by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (24)

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, [H.sub.0](B) : [[sigma].sup.2] ([[beta].sub.j]) = 0, j = l, ..., 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, [[bar.Y].sub.i..], [[bar.Y].sub..j.], and [[bar.Y].sub....] are replaced by [[??].sub.i..], [[??].sub..j.], and [[??].sub....], 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:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (25)

Thus, the proposed test statistics for the mixed effects and random effects models are, respectively,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (26)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (27)

where [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Proposition 1. Let [S.sub.B,mixed] be as defined in (26) and let [S.sub.B,rand] be as defined in (27), with the observations [Y.sub.ijk] generated according to the model of Section 2. Then,

(i) if the null hypothesis [H.sub.0](B) : [[sigma].sup.2]([[beta].sub.j]) = 0, j = 1, ..., b holds,

E([S.sub.B,mixed]) = E([S.sub.B,rand]) = 0, (28)

(ii) if the null hypothesis [H.sub.0](B) : [[sigma].sup.2]([[beta].sub.j]) = 0, j = 1, ..., b does not hold,

E([S.sub.B,mixed]) > 0, E([S.sub.B,rand]) > 0. (29)

The second item in Proposition 1 suggests that the null hypothesis should be rejected at level a whenever the test statistic takes a value larger than the l00 (l - [alpha])th percentile of the null distribution. The next proposition establishes asymptotic representations of [S.sub.B,mixed] and [S.sub.B,rand] that are useful for establishing their asymptotic distributions.

Proposition 2. Under [H.sub.0](B),

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (30)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (31)

Theorem 3. Let [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. Then under [H.sub.0](B)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (32)

in distribution as a, b [right arrow] [infinity].

Proof. Letting [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], then, as stated above,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (33)

The proof that [T.sub.b]/[s.sub.b] [right arrow] N (0, l) 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 [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], so that [T.sub.b] = [[SIGMA].sub.j] [Y.sup.a,b.sub.j]. Also note that we can write

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (34)

which is finite because we have assumed a finite fourth moment for the e terms. Moreover, this implies that

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (35)

which is again finite. Now in order to use the Lindeberg-Feller CLT, we need to verify the Lindeberg condition:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (36)

Note that

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (37)

Because E([square root of (b)][Y.sup.a,b.sub.1]) = 0 and Var([square root of (b)][Y.sup.a,b.sub.1]) is finite, the Dominated Convergence Theorem implies that [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], and hence the Lindeberg condition is verified. This completes the proof of the theorem.

Theorem 4. Define the operator [H.sub.b] : g [right arrow] [H.sub.b]g, for g [member of] [L.sub.2]([0, l], [B.sub.[0,l]], [[lambda].sub.[0,l]), where [B.sub.[0,l]] is the class of Borel subsets of [0, l] and [[lambda].sub.[0,l] is the Lebesgue measure on [B.sub.[0,l]], by

[[H.sub.b]g)(s) = [integral][h.sub.b] (s, t)g(t)dt, (38)

where [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. Let [[lambda].sup.b.sub.v], [f.sup.b.sub.v] = l,2, ..., be the eigenvalues and eigenfunctions of the integral equation

([H.sub.b]f)(s) = [lambda]f(s), (39)

and assume that the eigenvalues satisfy

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (40)

Then, if [[gamma].sub.ij] and have [[member of].sub.ij] finite fourth moments and the covariances Cov ([[gamma].sub.ij], [[gamma].sub.ij']), ([[gamma].sup.2.sub.ij], ([[member of].sup.-2.sub.ij] Cov ([[gamma].sup.2.sub.ij], [[gamma].sup.2.sub.ij]), and Cov ([[member of].sup.2.sub.ij], [[member of].sup.2.sub.ij]) are bounded uniformly in j, j', the asymptotic null distribution of [S.sub.B,rand], as a [right arrow] [infinity], b [right arrow] [infinity], is the same as the asymptotic distribution of [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], as b [right arrow] [infinity], where [Z.sub.v], v = l,2, ..., are independent standard normal random variables.

Proof. The proof proceeds by first conditioning on the levels of the column factor T = ([T.sub.1], ..., [T.sub.b]) 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 = ([T.sub.1], ..., [T.sub.b]). As b [right arrow] [infinity], the strong law of large numbers guarantees that the kernel [h.sub.b](s, t) converges to the same limit, and thus the conditional limiting distribution does not depend on T = ([T.sub.1], ..., [T.sub.b]). 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 term [s.sup.2.sub.b] in Theorem 3. According to (30) of Proposition 2 we need to estimate the variance of

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (41)

It is easy to see that

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (42)

in probability. Thus, we need to estimate the variance of the independent terms [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. By the Central Limit Theorem, as a [right arrow] [infinity], [[??].sup.2.sub..j.] is approximately Normally distributed with mean zero and variance [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. Therefore, the approximate variance of [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. As a result, we estimate by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (43)

Then, according to Proposition 1 and Theorem 3, [H.sub.0](B) is rejected at level when

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (44)

where [z.sub.[alpha]] is the (l-[alpha]) l00th percentile of the standard normal distribution.

Testing for main effects in the random effects model involves estimation of the asymptotic distribution of Theorem 4.Koltchinskiiand Gine [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 [H.sub.b] in Theorem 4 by its empirical version, [[??].sub.a], defined below. To this end, the following representation of [S.sub.B,rand] as a U-statistic is useful:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (45)

where the U-statistic [U.sub.a] is defined by the kernel

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (46)

Under the null hypothesis [H.sub.0](B), the empirical version, [[??].sub.a], of the operator [H.sub.b] is the a x a matrix with entries [[??].sub.i,i'] = (l/a)h([S.sub.i], [S.sub.i']). The a eigenvalues [[lambda].sup.a.sub.1], ..., [[lambda].sup.a.sub.a] of Ha are then used to estimate the asymptotic distribution of [S.sub.B,rand] as[ [SIGMA].sup.a.sub.i=1] [[lambda].sup.a.sub.1]([Z.sup.2.sub.i] - l), where the [Z.sup.2.sub.i] are randomly sampled from the [[chi square].sub.1] distribution. To approximate the estimated distribution, we randomly generate a large number B of sets [Z.sup.2.sub.1], ..., [Z.sup.2.sub.a], and use the formula above to obtain realizations from the distribution. Using this approximate distribution, for a level a test we reject [H.sub.0](B) when [S.sub.B,rand] is greater than the (l - [alpha]) l00th 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 usual F-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 values given from SAS.

Using the decomposition in relation (18), we generate our data for the mixed effects model as follows.

(1) Set [mu] equal to any constant.

(2) Generate a vector of main random effects [T.sub.j], j = l, ..., b from the standard normal distribution.

(3) Independently generate a vector of fixed "functions" [[lambda].sub.i], i = l, ..., a, from the standard exponential distribution and set [[gamma]*.sub.i] = [[lambda].sub.i] - [bar.[lambda]].

(4) Generate [Y.sub.ijk] from a normal distribution with variance [[sigma].sup.2.sub.i] and mean [mu] + [[alpha].sub.i] + [T.sub.j] + [T.sub.j] x [[gamma]*.sub.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 [mu] equal to any constant.

(2) Generate a vector of main row random effects [S.sub.i], i = l, ..., a from the standard normal distribution.

(3) Independently generate a vector of main column random effects [T.sub.j], j = l, ..., b from the standard normal distribution, and set [T.sup.*.sub.j] = [tau][T.sub.j], where [tau] is used to set the degree of deviation from the null hypothesis.

(4) Set the interaction effects [[gamma].sub.ij] = [S.sub.i][T.sub.j].

(5) Generate [Y.sub.ijk] from a normal distribution with variance [[sigma].sup.2.sub.i] and mean [mu] + [S.sub.i] + [T.sup.*.sub.j] + [S.sub.i][T.sub.j].

In Table 1, each simulation used a = b = 20. For the unbalanced setups, rows i = l, ..., 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, [[sigma].sup.2.sub.1] = xxx = [[sigma].sup.2.sub.15] = 1, [[sigma].sup.2.sub.16] = xxx = [[sigma].sup.2.sub.20] = 5.

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 data set, 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 data set, we employed our methodology to test for effects due to EDA.

5.1. Imputations. Because of the l/([n.sub.ij]-1) coefficient that is used in the statistic [S.sub.B,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 [[??].sub.t,j] = 0.046424 + 0.65852[CV.sub.t-l,j], where we predict a value of Mercury concentration at time from the Mercury concentration at time -l. For this model, the P value for the test of [H.sub.0]: [[beta].sub.1] = 0 versus [H.sub.A] : [[beta].sub.1] [not equal to] 0 was 0.000 and the [R.sup.2] = 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 [[??].sub.t,j] + [R.sup.*], where [R.sup.*] 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

[CV.sub.t,j] = 0.046424 + 0.65852[CV.sub.t-1,j] + [R.sup.*], (47)

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 as [e.sub.ijk] = [Y.sub.ijk] - [[bar.Y].sub.ij.], k = l, ..., [n.sub.ij]. Then for the EDA and year of interest, we impute a second observation as [Y.sub.ij2] = [Y.sub.ij1] + [e.sup.*.sub.ijk], where [e.sup.*.sub.ijk] 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 level [j.sup.*]) and calculate residuals for each year with more than one observation in the neighboring EDA as [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. Then for the EDA and year of interest, we impute a second observation as [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] denotes a randomly selected residual.

After performing the imputations, there are a total of = 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.

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 fact or 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 statistic [S.sub.C], 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.

The P valuefor thetestusing Sc is simply calculated as l-[PHI](X), where [PHI] is the distribution function of the standard normal random variable. The P value for the SAS F-test is calculated as [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], denote the numerator and denominator degrees of freedom. In the test for the random interaction effect, SAS uses [v.sub.1] = 722, [v.sub.2] = l028. The simulation results shown in Table 1 suggest that the significance of the interaction effect indicated by the -test cannot be trusted. Indeed, according to these simulation results, we expect the -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 classical F 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 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.

Appendices

A. Proof of Proposition 1

We first present the proof for the statistic under the mixed effects model, and then for the statistic under the random effects model.

Proof for [S.sub.B,mixed]. Using relations (15) and (18) we write

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (A.1)

The conditional expectation of the first part given the vector T of random column levels is

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (A.2)

Therefore, under the null hypothesis that [[beta].sub.j] = 0, a.s., for all j, the expected value of [S.sub.B,mixed] is zero conditionally on T, and thus also unconditionally. Moreover, it is clear that under the alternative the expected value of [S.sub.B,mixed] is positive.

Proof for [S.sub.B,rand]. Using relation (9), we write

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (A.3)

Thus, the proof of this part of the proposition will follow if we show

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (A.4)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (A.5)

Relation (A.4) follows from

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (A.6)

where the second equality follows from the fact that [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. Relation (A.5) follows similarly.

B. Proof of Proposition 2

We first show the proof for the statistic under the mixed effects model, and then for the statistic under the random effects model.

Proof for [S.sub.B,mixed]. Because we have shown above in the proof of Proposition 1 that, under [H.sub.0]([beta]),the statistic has mean value zero, it suffices to show that with data generated as in Section 2 and effects defined in (18), the following convergence holds:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (B.1)

in probability as a, b [right arrow][infinity]. Now, first note that the above expression is centered, so we only need to show that its variance tends to zero as a, b tend to [infinity]. This is easily done using the facts, given in relation (19), that the e terms are uncorrelated for different levels of the fixed effect, independent for different levels of the random effect, and, conditionally on the level of the random effect, iid for k = l, ..., [n.sub.ij]. Note that, using these relations,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (B.2)

Because we assume a finite fourth moment for the e terms and a finite covariance for any two distinct squared e terms, this conditional variance tends to zero. Therefore, by the Dominated Convergence Theorem, the unconditional variance also tends to zero, thus proving the proposition.

Proof for [S.sub.B,rand]. First, note that after some straightforward algebra we can write

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (B.3)

To complete the proof, we claim that [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] converges to zero after showing that the leading term of this product tends to zero. We proceed similarly for the last two terms in [S.sub.B,rand]. Finally, we show that [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] tends to zero to finish the proof. Note that each of these terms has expected value zero, so to demonstrate these facts, we only need to show that the variances tend to zero. To this end, note that

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (B.4)

This completes the proof of the proposition.

http://dx.doi.org/10.1155/2013/708540

Acknowledgments

This research was supported in part by NSF Grants SES0318200 and DMS-0805598. This research was initiated during the second author's sabbatical visit at the Australian National University in the Spring of 2006. Helpful discussions with Peter Hall are gratefully acknowledged.

References

[1] H. Scheffe, The Analysis of Variance, John Wiley &Sons, New York, NY, USA, 1959.

[2] S. R. Searle, Linear Models, John Wiley & Sons, New York, NY, USA, 1971.

[3] S. F. Arnold, The Theory of Linear Models and Multivariate Analysis, John Wiley & Sons, New York, NY, USA, 1981.

[4] S. R. Searle, G. Casella, and C. E. McCulloch, Variance Components, John Wiley & Sons, New York, NY, USA, 1992.

[5] A. I. Khuri, T. Mathew, andB. K. Sinha, Statistical Tests for Mixed Linear Models, John Wiley & Sons, New York, NY, USA, 1998.

[6] J. Cornfield and J. W. Tukey, "Average values of mean squares in factorials," Annals of Mathematical Statistics, vol. 27, pp. 907-949, 1956.

[7] B. Winer, Statistical Principles in Experimental Design, McGraw Hill, New York, NY, USA, 1971.

[8] R. R. Hocking, "A discussion of the two-way mixed model," The American Statistician, vol. 27, pp. 148-152, 1973.

[9] T. Gaugler and M. G. Akritas, "Mixed effects designs: the symmetry assumption and missing data," Journal of the American Statistical Association, vol. 107, no. 499, pp. 1230-1238, 2012.

[10] T. Gaugler and M. G. Akritas, "Testing for interaction in two-way random and mixed effects models: the fully nonparametric approach," Biometrics, vol. 67, no. 4, pp. 1314-1320, 2011.

[11] A. Bathke, "The ANOVA F test can still be used in some balanced designs with unequal variances and nonnormal data," Journal of Statistical Planning and Inference, vol. 126, no. 2, pp. 413-422, 2004.

[12] A. K. Gupta, S. W. Harrar, and Y. Fujikoshi, "Asymptotics for testing hypothesis in some multivariate variance components model under non-normality," Journal of Multivariate Analysis, vol. 97, no. 1, pp. 148-178, 2006.

[13] M. G. Akritas and N. Papadatos, "Heteroscedastic one-way ANOVA and lack-of-fit tests," Journal of the American Statistical Association, vol. 99, no. 466, pp. 368-382, 2004.

[14] S. R. Searle, Linear Models for Unbalanced Data, John Wiley & Sons, New York, NY, USA, 1987.

[15] V. Koltchinskii and E. Ginee, "Random matrix approximation of spectra of integral operators," Bernoulli, vol. 6, no. 1, pp. 113-167, 2000.

[16] J. Aitchison, The Statistical Analysis of Compositional Data, Chapman & Hall, London, UK, 1986.

Trent Gaugler (1) and Michael G. Akritas (2)

(1) Department of Statistics, Carnegie Mellon University, 5000 Forbes Avenue, 132B Baker Hall, Pittsburgh, PA 15213, USA

(2) Department of Statistics, Penn State University, 325 Thomas Building, University Park, PA 16802, USA

Correspondence should be addressed to Trent Gaugler; tgaugler@andrew.cmu.edu

Received 3 August 2012; Accepted 4 February 2013

Academic Editor: Jose Sarabia

Table 1: Comparison of achieved type I error rates and power for [S.sub.B,mixed], [S.sub.B,rand],and F with N = l000 tests for the random and mixed effects Random effects model Balanced [alpha] = .05 [alpha] = .10 [S.sub.B] F [S.sub.B] F Homo. [H.sub.0] .054 .220 .111 .247 [A.sub.1] .053 .229 .123 .265 [A.sub.2] .125 .337 .221 .367 [A.sub.3] .236 .532 .369 .563 Hetero. [H.sub.0] .045 .150 .089 .182 [A.sub.1] .045 .160 .088 .199 [A.sub.2] .114 .273 .190 .316 [A.sub.3] .195 .430 .318 .467 Random effects model Unbalanced [alpha] = .05 [alpha] = .10 [S.sub.B] F [S.sub.B] F Homo. [H.sub.0] .049 .167 .092 .198 [A.sub.1] .069 NA .140 NA [A.sub.2] .113 NA .201 NA [A.sub.3] .234 NA .360 NA Hetero. [H.sub.0] .015 .230 .042 .303 [A.sub.1] .021 NA .053 NA [A.sub.2] .047 NA .094 NA [A.sub.3] .117 NA .215 NA Mixed effects model Balanced [alpha] = .05 [alpha] = .10 [S.sub.B] F [S.sub.B] F Homo. [H.sub.0] .065 .047 .112 .099 [A.sub.1] .227 .187 .318 .296 [A.sub.2] .761 .719 .820 .804 [A.sub.3] .975 .963 .982 .982 Hetero. [H.sub.0] .070 .056 .106 .100 [A.sub.1] .077 .060 .129 .121 [A.sub.2] .151 .129 .221 .215 [A.sub.3] .270 .235 .353 .341 Mixed effects model Unbalanced [alpha] = .05 [alpha] = .10 [S.sub.B] F [S.sub.B] F Homo. [H.sub.0] .063 .045 .109 .086 [A.sub.1] .184 .138 .260 .228 [A.sub.2] .542 .489 .634 .638 [A.sub.3] .888 .863 .920 .917 Hetero. [H.sub.0] .062 .935 .107 .952 [A.sub.1] .074 NA .124 NA [A.sub.2] .104 NA .152 NA [A.sub.3] .131 NA .189 NA

Printer friendly Cite/link Email Feedback | |

Title Annotation: | Research Article |
---|---|

Author: | Gaugler, Trent; Akritas, Michael G. |

Publication: | Journal of Probability and Statistics |

Article Type: | Report |

Date: | Jan 1, 2013 |

Words: | 6970 |

Previous Article: | Depth-based classification for distributions with nonconvex support. |

Next Article: | Estimation of extreme values by the average conditional exceedance rate method. |

Topics: |