Exact critical points for testing subhypotheses in nonreplicated three-way experiments

In the absence of replication, conventional analyses do not provide ways to examine three-way interaction in threeway experiments. Tucker3 analysis based on alternating least squares algorithm is a general approach that can be used in such cases. However, Tucker3 options are not available with standard statistical packages. A few methods for estimating σ 2 and testing for three-way interaction using a "single component" Tucker3 model are available in literature. A method based on a convenient approximation to a likelihood ratio test is also available in limited cases for testing three-way interaction in sub-areas once interaction is present. In this paper the null distribution of the above likelihood ratio statistic is simulated using Monte Carlo methods based on exact values obtained from Tucker3 analysis. Critical points are also obtained for the test for selected cases. Though the package 3-WAY PACK ® handles Tucker3 analysis it does not conveniently support WINDOWS ® based simulations and therefore MATLAB ® is used for the simulation. The method is illustrated using two examples involving real data. Keywords: Exact critical points, testing subhypotheses, three-way interaction, tucker3 analysis. doi:10.4038/jnsfsr.v36i4.264 Journal of the National Science Foundation of Sri Lanka 36 (4) 267-273

model, alone is an exploratory approach.A few confirmatory inferential approaches are also available in literature for testing various hypotheses about threeway interaction 3,4,8 .The focus of this paper is a follow up of the work by Wickremasinghe and Johnson 8 .They developed a likelihood ratio test for testing subhypotheses in nonreplicated three-way experiments.The null distribution of some convenient approximation of this test statistic has been simulated, and critical points computed by them using Monte Carlo methods, for selected cases.Thus, the critical points for the test were only approximate.
In this paper, the above null distribution is re-examined for exact values of the test statistic.This is based on Monte Carlo simulations that require running the alternating least squares algorithm iteratively 9 .The software package MATLAB ® is used for this task.Critical points are computed for selected cases based on exact values.A numerical example with real data is used for illustrating the methods.Finally, a comparison is made between the approximate and exact results based on pre-analyzed data.

METHODS AND MATERIALS
Modelling nonreplicated three-way data: Let the subscripts i, j, and k represent a general level for each of the three factors, say A, B, and C, having ℓ, m, and n levels respectively.In the absence of replication, usual model for nonadditive three-way data for any general response y ijk is given by

INTRODUCTION
The primary objective of a three-factor factorial experiment, or simply a three-way experiment, is to examine the interaction between the three factors simultaneously.If the interaction is present, the interpretations would be somewhat messy; otherwise they would be straightforward.When the experiment is replicated, an analysis of variance (ANOVA) can be used to test for three-way interaction.If the experiment is nonreplicated, the usual ANOVA fails.For analyzing data in such cases, the Tucker3 model 1,2 has been considered as a general approach by many [3][4][5][6][7] .Tucker3 analysis, i.e. analysis of three-way data using the Tucker3 For i = 1,2,....,ℓ: j = 1,2,....,m; and k = 1,2,....,n, where μ is the general mean and α i , β j , and γ k denote the corresponding levels of main effects of the three factors; τ ij is the interaction for the i th level of A and j th level of B. Similarly, ρ ik and ø jk are defined and π ijk is a three-way interaction term.It is assumed that all parameters except the random error term are fixed and subject to usual sumto-zero restrictions, and ε ijk are assumed to be i.i.d.N(0,σ 2 ), where σ 2 is unknown.For the purpose of estimating σ 2 , conventionally it has been assumed that π ijk = 0.However, this assumption is not valid in the presence of three-way interaction.A previous study  5) and ( 6) are similarly defined.In practice, one needs to use only one of the versions as the results are invariant to the model expression.A short review on Tucker3 model and its applications is given by Wickremasinghe 10 .
Critical points for the likelihood ratio test proposed by Marasinghe and Boik 4 for testing c = 0 can be read from tables computed by Boik 11 with n 1 = ℓ -1, n 2 = m -1, and n 3 = n -1.If the result of the test is that there is no three-way interaction, then a conventional analysis of main effects and two-factor interactions would be made.If the result says three-way interaction is present, one may obtain a suitable estimate of σ 2 and examine the pattern of this interaction and then analyze the data on those lines.Another option in this situation is to locate a sub-area free of three-way interaction, and then use it for conventional analysis.This is particularly useful when the sub-area consists of all the data except a few observations.In the latter, Wickremasinghe and Johnson 8 developed a likelihood ratio test for testing subhypotheses in model ( 2).This model has been commonly used in literature mainly because s = t = u =1 in the Tucker3 model gives an adequate fit to most small to moderate three-way data sets.This study is a follow up of their work.

Testing
subhypotheses in three-way data: Wickremasinghe and Johnson 8 proposed a likelihood ratio test for testing H 0 : Sg = 0, Th = 0, and Ue = 0 where S q 0 ×ℓ ,T r 0 ×m , and U s 0 × n are contrast matrices of ranks q 0 (< ℓ -1), r 0 (< m -1), and s 0 (< n -1) respectively.If not rejected, the test would help identify a sub-area free of three-way interaction.The test is to reject H 0 for small values of where D ℓ× mn is one of the three matrix versions of the vector d ; S -,T -and U -denote the Moore-Penrose g-inverses of the corresponding matrices, respectively.Wickremasinghe and Johnson 8 simulated the null distribution of Λ for q 0 = r 0 = s 0 = 1 after replacing ĉ 2 by the min(ℓ 1 , m 1 , n 1 ) and ĉ *2 by the largest characteristic root of D * D * ' , where ℓ 1 , m 1 , and n 1 denote the largest characteristic root of DD' for the three matrix versions respectively.They also computed critical points for the test based on this approximation.The reason for using the approximation has been twofold: First, See and Smith 3 suggested that the above approximation works well for model (2) when the test for c = 0 is significantly rejected; Second, it was the lack of proper software at the time to compute the actual values of ĉ 2 and ĉ *2 , which involve running the alternating least squares algorithm repeatedly for each case.

Simulation study
To carry out the above likelihood ratio test 8 one needs to compute critical points based on the exact null distribution of the statistic Λ.Since the theoretical derivation of the null distribution of Λ is extremely difficult, the only other alternative is to simulate the null distribution using Monte Carlo methods.
In this study, the null distribution of Λ was simulated for several cases using the actual values of ĉ 2 and ĉ *2 .This was done for the case q 0 = r 0 = s 0 = 1, and without loss of generality assuming σ 2 =1.
The ℓmn-vector of residuals, d, can be transformed into an (ℓ-1) (m-1) (n-1)-vector Now, by replacing Var(d) by σ 2 E 0 and using the properties of the contrast matrices, one gets Var(z)=σ 2 I (ℓ−1)(m−1)(n−1) which gives a reduction in dimension as well as independent components.These results were used by us following the same approach discussed by Wickremasinghe and Johnson 8 for generating Λ under the null conditions, as follows: The vector z was generated as z~N(µ,1) where µ is an (l−1)(m−1)(n−1) × 1 vector of all zeros except that the first entry has value δ(= | c|/σ).Then vector z * was defined as the first (l−2)(m−2)(n−2) elements of z.Using these two vectors z and z * , ĉ 2 and ĉ *2 were obtained through a MATLAB ® programme Tuck3F.m (Appendix) which was called as a function inside another MATLAB ® programme Simul.m (Appendix).Tuck3F.m performs Tucker3 analysis using alternating least squares algorithm as explained by Kroonenberg and De Leeuw 9 .Several combinations of l, m, and n with l ≥ m ≥ n, and for each combination δ ranging from 0 to 15, were selected for the simulation.The programme Simul.m was used for generating 1000 values of Λ for each l, m, n, and δ, and to obtain the average of ĉ 2 over all 1000 values in each case.
The empirical distribution was studied for the cases 3×3×3 through 8×8×8 using the generated Λ values, and the exact null distribution was approximated by a beta distribution 8 by the method of moments after smoothing the mean and the variance of the generated values using a nonlinear function of the form f (δ)=α(1-be -cδ 2 ).except for the signs of ĉ and some elements of the three vectors.This is quite acceptable due to properties of Tucker3 solution 10 .

Examples
A residual analysis was also carried out to check the assumptions on the errors in model ( 2), and it was found that the assumptions were not violated.
All possible 54 hypotheses of rank 1 of the form H 0 : g i = g i' , h j = h j' , and e k = e k' were tested using our were increasing functions of δ and the relationship could well be modeled by the above nonlinear function.The goodness-of-fit of beta was tested using beta Q-Q plot. Figure 1 shows a sample of beta Q-Q plots for the case 3×3×4 for δ = 1, 5, and 10.It was found that beta fitted very well for the null distribution of Λ in the smaller cases 3×3×3, 3×3×4, 3×3×5, 3×3×6 and 3×4×4 when the exact values were used.For these cases, 1%, 5% and 10% critical points for the test were computed for δ=0-10, and are given in Table 1.For δ > 10, one can read the table for δ =10 as the critical points were approximately the same.The reason for this is the fact that the nonlinear function f (δ)=α(1-be -cδ 2 ) leveled off at values smaller than δ=10 for both the mean and the variance in the cases studied.As a result, the smoothed parameters were constant after some point, depending on the combination but definitely for values of δ < 10.This resulted in giving the same critical points after some value for δ.One can easily see that in Table 1 for the combination (3, 3, 3) the critical points were almost constant for δ > 6, and for (3, 3, 4)  this was true for δ > 8.For the other combinations, the critical points were approximately the same for δ > 10.

Experiment (2)
Average girth increments over 4 consecutive periods for 3 Hevea clones tested under 4 different densities were analyzed.The data with real identifications suppressed, are given in Table 2.The residuals from (1) withs π ilk = 0, of the 3×4×4 dataset in Table 2 were decomposed using single component Tucker3 model.The residual sum of squares was 7.9812 with 18 degrees of freedom.The initial test by Marasinghe and Boik 4 suggested the presence of three-way interaction.The estimates of the parameters given by Simul.m (Appendix) are ĉ 2 = 5.890, ĝ' = (-.5847, .7864, -.1993, .0005),ĥ' = (-.8587,.2050,.2666,.3866),and ê' = (-.5771,-.2141, .7881)where the vectors g, h, and e correspond to Period, Density and Clone, respectively.A residual analysis was also carried out to check the assumptions on the errors in model ( 2), and it was found that the assumptions were not violated.The estimate 8 of *2 ˆ which involves δ 111 (= 9.88) 1 are .5009and .3651respectively.All possible 108 hypotheses of rank 1 of the form H 0 : g i = g i h j = h j' , and e k = e k ' were tested using these critical points.All but 9 hypotheses were rejected at 1% level.The hypotheses that were not rejected can be combined as H 0 : g 1 = g 3 = g 4 , h 2 = h 3 = h 4 and e 1 = e 2 .Using the methods discussed earlier 8,12 , it is easy to see that the sub-area identified by this combined hypothesis is nothing but the one without the single observation corresponding to (g 2 , h 1 , e 3 ).This is actually the value corresponding to Period 2, Density 1, and Clone 3. From Table 2, it appears that it is the highest value (3.86).This suggests that if one analyzes the data in Table 2 without this observation, then the problem of three-way interaction does not arise, and one can interpret main effects and two factor interactions in the usual way.A detailed analysis of this data along with interpretations will be discussed elsewhere.

CONCLUSION
In this paper we have obtained exact critical points to go with the likelihood ratio test 8 to identify a sub-area free of three-way interaction.These critical points were obtained using Monte Carlo methods with the help of a few MATLAB ® programmes.The beta approximation worked well only for the smaller cases 3×3×3, 3×3×4, 3×3×5, 3×3×6, and 3×4×4.the two examples with real data illustrate the methods.In the first example the conclusions obtained were the same as those made by wickremasinghe and Johnson (2002) where a convenient approximation of the test statistic has been used for computing critical points.It appears that the approximation 8 worked well for the 3×3×4 data set.However, this may not be true in general.The second example demonstrated that the method can be used as a viable alternative in the analysis of crop data in certain situations.This paper did not focus on estimating σ 2 in the presence of three-way interaction, as this has been previously discussed 4,8 .Work is underway involving cases beyond 3×4×4.

Acknowledgment
This work is supported by the research grant RG/2005/ FR/01 from the National Science Foundation of Sri Lanka.The authors wish to thank the Rubber Research Institute of Sri Lanka for permission to use their data for example (2), and also the referees for their valuable comments.
, where d is as above, ĉ 2 is the maximum likelihood estimate of c 2 as given by Marasinghe and Boik 4 , and ĉ *2 is the square of the largest three-mode singular value of D * = (I -S -S) D ((

Figure 1 :
Figure 1: Beta Q-Q plots for 3 ×3×4 and selected values of δ 4considered the generalized Mandel-Johnson-Graybill type model to three factors where π ijk is replaced by a multiplicative interaction term c g i h j e k with new parameters c, g i , h j , e k with assumptions that g i , h j , and e k , each sum to zero over the respective subscripts, and the corresponding vectors g, h and e have unit norm.This study proposed a likelihood ratio test for testing π ijk = 0 in (1) by testing c = 0.This generalized model for three factors is given by

Table 1 :
One percent, 5%, and 10% critical points for the likelihood ratio test for selected cases of l, m, and n with l ≥ m ≥ n and δ

Table 2 :
Girth increments (cm) over 4 consecutive periods for 3 Hevea clones tested under 4 different densities