Analysis of Variance GLM Version 3.0 Using Mathematica V4.0
Steve Hunka, University of Alberta,
Edmonton, Alberta, Canada, T6G 2G5
steve.hunka@ualberta.ca
Mod. Date: July 20/02

General Description

The function anovanw[d] calculates an analysis of variance (fixed effects, factorial designs) following general linear model procedures using the data matrix d. Two forms of data input are allowed: (a) ungrouped with each observation having a cell index, or (b) all observations for each cell index are grouped together. The user can select a subset of the cell indices to be used. The maximum number of factors in the design is currently set to 7, and can be extended easily.  A general Anova table of sum of squares (SStotal, SSmodel, SSerror) and an Anova table of sum of squares (Type III) for each term in the model ( a,b,c,ab,ac,bc, abc,...) is provided with F tests of significance when df(error) is greater than 0. The user is given the options to: (a) remove interaction terms of the model, (b) plot 2-term and 3-term interactions,  and (c) enter user contrast matrices to carry out contrasts on the main effects. OLS parameter estimates are calculated in sum reduced form, i.e., for each set of terms the parameter estimates sum to zero. Cells may vary in the number of observations (minimum of one observation per cell). 

A function (chkidx) is provided which will check the adequacy of cell indices.

The notebook anova.nb does not provide for repeated measure and nested designs, however, these designs can be re-cast into a factorial design with one observation per cell to obtain sum of squares which can then be combined to form the sum of squares required for repeated measure and nested designs.  To facilitate the analysis of such designs the notebook MGems.nb can be used to identify the required sum of squares, df, and expected mean squares, and the mean square terms to use in the denominator of F tests.  An example of a nested design and a repeated measure design is included in the Anova.nb notebook.

All functions were tested using Version 4.0 for the Macintosh (G4).

References
Kirk, R.E. (1982) Experimental Design: Procedures for the Behavioral Sciences, 2nd Edition, Brooks/Cole Publishing Co., Pacific Grove, CA.
Millman, J.,  &  Glass, G.V. (1967) Rules of Thumb  for Writing the Anova  Table, Jr. of Educational Measurement, Vol. 4, No. 2, Summer.
Searle, S.R. (1971) Linear Models, John Wiley & Sons, N.Y.
Shavelson,R.J. & Webb, N,M., (1991) Generalizability Theory A Primer, Sage
               Publications Inc., Newbury Park, CA.
Winer, B.J. (1971) Statistical Principles in Experimental Design, McGraw-Hill, N.Y.
www.nist.gov/itl/div898/strd/anova/anova.html (Source of NIST data sets)

Program Modifications to V2.0 Incorporated in V3.0 for Mathematica V4.0

1) A new section entitled "Rationalizing Input Data to Obtain High Accuracy" has been added to illustrate the effects of rationalizing real data. Rationalization can be accomplished using the built-in Mathematica function Rationalize or a function which uses the Mathematica functions RealDigits and FromDigits. 
2) Cell variances are now calculated as maximum likelihood estimates (MLE using division by n instead of (n-1)) so that in the case of 1 observation per cell, a 0 variance is calculated. Cell information is output keeping cell indices within { and} rather than as separate columns.
3) Under user option, the cell means and variances can be output as Real or Rational values. 
4) The main sum of squares are displayed in both Rational and Real form.

Description of Procedure and Cautions

The following sequence of steps are carried out:
1. The nature of the input data is assessed to determine whether the data is ungrouped or grouped. If the data is grouped and the user selects to use fewer cell indices than available, the data is restructured (data of common cells are combined). For example, if the data is grouped for a 3-way anova (axbxc), but the user selects only 2 of the cell indices (a and c), then data is combined over index b providing for a 2-way anova. If the data presented are all Real values, they are converted to a rational form.
2.  The input data are sorted by cell index, and the number of observations, means, and variances are calculated for each cell. The values for SStcfm (total corrected for the mean), SSmcfm (model corrected for the mean, SSerror and R-sqd are displayed in rational and real form. The number of levels in each factor is obtained from the maximum value of each cell index. The number of parameter estimates which will be calculated for the design is displayed, and the user is given the option to stop.
3.  The sum reduced form of the design matrix (X) is formed for the main effects. It will have as many columns as degrees of freedom for the main effects, and as many rows as cells in the anova design. The design matrix is created as if there is one observation per cell to keep its size to a minimum. ( Later X'X is calculated to reflect cell sizes.) A list structure is maintained and updated as required which identifies each anova term, its binary representation, and its position in X, e.g., {"a",{0,0,1},1} indicates that anova term "a" has the binary representation {0,0,1} and has position 1 in the structure of X. The binary terms are used to determine which column vectors of X are to be used to form products required to define interaction parameters.
4.  Column vectors are generated and appended to X for the interaction terms. The column vectors whose elements must be multiplied in order to obtain a vector for an interaction term is determined by the difference between the binary term representing the interaction, and the binary terms available for vectors already calculated. For example, if term "abc" (1,1,1) is required, a backward search is made of existing terms until the binary difference of (1,1,1) and the existing term has only one 1, e.g., (1,1,1) - (1,1,0)=(0,0,1) indicates that "abc" may be formed by using the vectors of "bc" (1,1,0) and "a" (0,0,1).
5.The user is given the option to identify interactions terms to be removed from the model. The SS for these terms will then be accounted in the SS error. With 1 observation per cell, df(error) will be 0 unless some interaction terms are removed from the model.
6. The final X'X matrix is formed reflecting cell sample sizes by, in effect, forming [N*X'] X where N is a matrix of cell sample sizes of order (number of parameters by number of cells) each row being identical, * representing an element by element multiplication, and X' of the same order as N.  X'X is then inverted and used to find the parameter estimates using the equivalent of Inverse[X'X]X'Y but formed by Inverse[X'X]X' (m*n) where m is a vector of means and n a vector of cell sizes. The parameter estimates, their df, and values are displayed.
7. The first SS table containing SStcfm, SSmcfm, and SSerror is displayed with an F test for SSmcfm. A second table giving the SS for each term in the model (e.g., a, b, ab) with F tests is displayed. The SSerror in the tables reflects any changes to the full model brought about by dropping of intereaction terms. 
8.  For each term a contrast matrix in sum reduced, row echelon form (an identity matrix of order equal to the df for a term) is created, the associated parameter estimates and submatrix of X'X is selected and the SS for the contrast obtained by (KB)'Inverse[KInverse[X'X]K'](KB). (See Searle, 1971) The term, SS, df, Ms, F, and probability are displayed if df(error) is greater than 0. If df(error) is equal to 0, the calculations of (9) and (10) below are not made. 
9.  If the number of factors is greater than 1, the binary representation of available terms is searched for interaction terms.  For 2 and 3-term interactions, the user has the option of plotting the means. For example, if the interaction for term "ab" in a 4 factor model is to be plotted, means are calculated by summing the means over factor "c" and "d".
10. The user is given the option to carry out F tests on main effect contrasts. For example, if the term "a" has df=2, the user could enter the contrast matrix {{1,0,-1}} and test whether a[1]-a[3] = 0 in the population. Contrasts are entered in non-sum reduced form by the user, and adjusted to sum-reduced form. The row-reduced, and row echelon form of the user contrast is displayed to assist the user to identify any linear dependencies by row. (Row reduction of the contrast matrix to identify linear dependencies works fine for matrices in rational form, but not in real form.)  A symbolic representation is displayed of the hypothesis being tested. 

Note 1: To extend the design to more than 7 factors, include additional characters beyond "g" in the function bin2alf.
Note 2: When probabilities of F are displayed they are rounded to 4 decimal digits. Thus, a probability of 0 would be reported as <.0001 .

Cautions and Error Messages

1.  The number of parameters increases exponentially with the number of factors in the anova design. The number of terms for which estimates are made is the product of the number of levels in each factor. Designs having many factors make for more difficult interpretation. 
2.  No adequate check is made of linear row dependencies of user entered contrast matrices for main effects.  Mathematica will note an inversion problem if such dependencies exist. 
3. User input of a numeric value, a vector, or a matrix is checked for structure, but not for the range of values entered.  If a vector is required, but the user enters a scalar, a message will be displayed and the user is requested to re-enter the vector. If user input is a string, the input is checked against a set of acceptable answers. If the input does not match one of the acceptable answers, a message is given and the user is requested to re-enter the string. For example, if the user has removed the "abc" interaction in a 3-way anova, only 2-term interactions are available as acceptable answers in reply to which interactions are to be plotted. 
4.  No check is made that cell indices input by the user are within range of available cell indices. 
5.  No check is made that the number of factors exceeds 7, since this restriction can be easily removed by changing the number of characters available for use in the function bin2alf.
6.  In the case that each cell has one observation only, a Warning message is given that df(error) will be equal to 0, and that some interaction terms need to be removed in order to obtain F tests. If no interaction terms are removed the SS, mean squares, and df only are displayed. Designs having 1 observation per cell are used in Generalizability  Theory (Shavelson, R.J. & Webb, N.M., (1991) Generalizability Theory A Primer, Sage Publications, Newbury Park, CA) and usually define the number of levels in one of the factors as equal to the total number of observations. Thus, such designs may require a very large number of parameters to be estimated using the general linear model approach.

General Description of Interactive User Input

The user must be prepared to answer the following queries during execution of the anova:
1.  The number of factors in the design. The maximum is 7, and must not exceed the number of cell indices available with the data input.
2. The cell indices to be used for the design. For example if there are 4 cell indices associated with the data, if the user wanted to carry out a 4-way design, a vector holding any ordering of the integers 1 to 4 is acceptable, e.g., {1,2,3,4}, or {2,4,3,1}, and so on; using the same data for a 3-way design, any order of 3 of the 4 integers indices would be acceptable, e.g., {1,2,4}, {2,1,3},and so on. The assignment of factor identifiers, eg., a, b, c, is made by associating the first index with "a", the second with "b", and so on.
3.  The user is asked whether cell means and variances are to be output in rational or real form.
4. The user is given the opportunity to stop execution after a display of the cell frequencies, means, variances, and number of parameters which will be estimated, because: (a) by observing the cell frequencies and means, the user may sense that something is wrong with the data as input; (b) the user may have underestimated the number of parameters to be estimated. 
5.  The user is asked whether any interaction terms are to be removed, and if the reply is y the user is asked for each interaction term to be removed.  

If df(error) is greater than 0, 

6. The user is asked whether interaction terms are to be plotted if the number of factors in the design is equal to or greater than 2. If the reply is y, the user is requested to enter each interaction term separately. Only 2-term (e.g., ab, ac, ad,...) and 3-term interactions (e.g., abc, abd,...) can be plotted. Prior to this query, a display is made of the available interaction terms, since some terms may have been removed earlier (see 5 above). The term entered must have characters in their natural order, e.g., abc, and not acb. 
7.  The user is asked whether contrasts are wanted for main effect terms. If the reply is y, the user is asked to input a single main effect term, e.g., a, b, d, etc. The user is then asked to input a contrast matrix, e.g., {{1,0,-1}}. A single row is entered as a matrix.  The number of rows in the contrast matrix must not exceed the degrees of freedom for the term, e.g., if factor b has 4 levels and thus df=3, then the matrix can have 1, 2, or 3 rows only, and will require df+1 columns because the contrast matrix is entered in non-sum reduced form from which the sum-reduced form will be calculated automatically. For example {{1,0,-1}} will be reduced to {{2,1}} and further to {{1, 1/2}}. (Any elementary row operation on a contrast matrix will test an equivalent hypothesis and the SS will remain the same.)

Data Input

The input data are expected to be ungrouped and have the form {{{i1,i2,i3,...},dep var obs 1}, {{i1,i2,i3,...},dep var obs 2}, ....} in which {i1,i2,i3,...} holds each observation's factor index (cell index), or grouped in the form {{{i1,i2,i3,...},{all observations for this cell},....}  The user is allowed to select and reorder the indices.The data are sorted on the indices selected by the user so the data need not be presented ordered by groups or cell index. The cell indices are expected to form a continuous natural sequence starting at 1, e.g., 1,2,3 for 3 levels and not, for example, 2,4,5. The values of the dependent variable must be in the form of all integers or all real values. If all real values are present, they are converted to rational values to obtain maximum accuracy. If the input data contains real and integer values, a message is given and the computation aborted. Thus, for example, real values of zero should be given as the Real value 0. or .0, but not in the form of an Integer 0 without a decimal point. Null values must not be present in the data, e.g., {1,2,3,,} would produce a Null value for the 4th element.  If the user forgets to evaluate the data set, an error message is given and the run aborted.

The function chkidx can be used to check the adequacy of the vectors holding the cell indices before an anova is carried out.

Example Data Sets

Textbook  Data Sets
dat1: this data set was generated to fit the "summary tables" of a 3 factor anova problem detailed by Winer in the text "Statistical Principles in Experimental Design", McGraw-Hill, 1971, pp. 457-463. The fourth index was created to provide a 4-way anova for testing purposes.The data set is in ungrouped form, contains Real values, and is an equal n case.N=120.

dat2: this data set contains the same dependent variable as in data set dat1 but has only 3 cell indices (the same as the first 3 in dat1), and all observations belonging to the same cell are grouped together.  

dat3: this data set is also from Winer (p446) and is a 2-factor unequal-N case. The data is in ungrouped format and contains Integer values. N=33. Note that the solution given by Winer is not for a least squares solution which is used in the anova.nb notebook.

dat4: this data set has 1 observation per cell and 3 factors. It is from R.J.  Shavelson and N.M. Webb (Generalizability Theory, a Primer, Sage Publications, 1991, p43).The data set is in ungrouped form and contains Integers. dfe=0 for this data set. N=60 for a {6x2x5} design.

dat5 is for a 1-factor anova having 3 groups of unequal size; data is grouped and in Integer form; N=21. Original source Winer (p153,1971) and some observations dropped.

NIST Simon-Lesage Test Data Sets
Test data sets for assessing the accuracy of calculating sum of squares are provided by the web site
www.nist.gov/itl/div898/strd/anova/anova.html. These data sets are provided by the National Institute of Standards & Technology, Information Technology Laboratory. The www site provides information on how extended real-value precision was obtained. The functions gensimon and genunsimon can be used to generate the data sets Simon-Lesage 7 (N=189), 8 (N=1809), and 9 (N=18009). The function gensimon generates Real data in grouped format and genunsimon in ungrouped format. The "certified values" for each data set provided by the www site are as follows:
Simon-Lesage7: SStreatment=1.6800....0E+00; SSerror=1.80...0E+00; F=2.10....0E+01
Simon-Lesage8: SStreatment=1.6080....0E+01; SSerror=1.80...0E+01; F=2.010..0E+02
Simon-Lesage9: SStreatment=1.60080..0E+02; SSerror=1.80...0E+02; F=2.0010.0E+03

Mathematica allows computation to proceed in rational form providing for infinite precision.  The SS and R-sqd values which are obtained for these data sets using anova.nb are as follows:

Simon7: SStcfm=87/25; SSerr=9/5; SSmcfm=42/25; R-sqd=14/29; (N=189)
Simon8: SStcfm=852/25; SSerr=18; SSmcfm=402/25; R-sqd=67/142; (N=1809)
Simon9: SStcfm=8502/25; SSerr=180; SSmcfm=4002/25; R-sqd=667/1417; (N=18009)

SStcfm is SS "total corrected for the mean" i.e., deviated about the grand mean using Searle's terminology. SSmcfm is "SS model corrected for the mean", also called SStreatment or SS between.
 By using, for example, N[667/1417,1000] the first thousand digits of the R-sqd for Simon9 data can be obtained.

Other Data Sets

Two data sets from Wolfram (www.wolfram.com/Mathematica/newin42) are included. Also, a nested design from Kirk (p460) and a repeated measures design from Winer (p256) are included. The latter two designs are used to illustrate how the results from a factorial analysis can be used to obtain the appropriate sum of squares and F tests for nested and repeated measure designs.

About the Output

The results displayed will vary with some of the user options, e.g., whether plots are made, and whether user contrast matrices are entered. The following output is common for all problems:
1) Output describing the general parameters of the problem, e.g., number of factors, factor indices selected to be used for grouping data into cells, number of levels in each factor, and the nature of the data (all Real, all Integer, or mixed).
2) The frequency, mean, and variance for each cell. Note that the variances are calculated using a denominator of n, thus for one observation per cell the cell variance will be 0.
3) Total sample mean, total sample variance, SStcfm, SSwithin (SSerr),
Full Model SSmcfm, and Full model R-sqd are output in Real and Rational form. The last two calculations are made directly from the data, and based on holding all terms of the model.  If the user decides to drop interaction terms, providing for a "restricted model" the values for SSmcfm and SSerr are recalculated for display in the SS tables and used for making F tests.
4) Sum-reduced parameter estimates are output in Real form using the calculations made in rational form.  Similarly, SS tables provide results in Real form but calculated in rational form.