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
Revision Date: July 20/02

The readMe.nb notebook contains the following information about this file:
1)  General Description of anova procedure used here
2)  Program Modification to V2.0 Incorporated in V3.0 for Mathematica V4.0
3)  Description of Procedures and Cautions
4)  General Description of Interactive User Input
5)  Data Input
6)  Example Data Sets (a description and source of the data sets in this notebook)
7)  About the Output

Example Data Sets

dat1 is in ungrouped format with 4 cell indices; dependent variable is in Real form. N=120. Example used in a {2x3x2} equal n design.

dat1 = {{{1, 1, 1, 1}, 0. }, {{1, 1, 1, 1}, 1. }, {{1, 1, 1, 1}, 
            1. }, {{1, 1, 1, 1}, 2. }, {{1, 1, 1, 1}, 2. }, {{1, 1, 1, 2}, 
            2. }, {{1, 1, 1, 2}, 2. }, {{1, 1, 1, 2}, 3. }, {{1, 1, 1, 2}, 
            3. }, {{1, 1, 1, 2}, 4. }, {{1, 1, 2, 1}, 0. }, {{1, 1, 2, 1}, 
            0.5}, {{1, 1, 2, 1}, 0.5}, {{1, 1, 2, 1}, 1. }, {{1, 1, 2, 1}, 
            2. }, {{1, 1, 2, 2}, 2. }, {{1, 1, 2, 2}, 2.5}, {{1, 1, 2, 2}, 
            2.5}, {{1, 1, 2, 2}, 3. }, {{1, 1, 2, 2}, 2. }, {{1, 2, 1, 1}, 
            2. }, {{1, 2, 1, 1}, 1. }, {{1, 2, 1, 1}, 1. }, {{1, 2, 1, 1}, 
            2. }, {{1, 2, 1, 1}, 3. }, {{1, 2, 1, 2}, 4. }, {{1, 2, 1, 2}, 
            4. }, {{1, 2, 1, 2}, 4. }, {{1, 2, 1, 2}, 4. }, {{1, 2, 1, 2}, 
            5. }, {{1, 2, 2, 1}, 2. }, {{1, 2, 2, 1}, 2. }, {{1, 2, 2, 1}, 
            1. }, {{1, 2, 2, 1}, 3. }, {{1, 2, 2, 1}, 3. }, {{1, 2, 2, 2}, 
            4. }, {{1, 2, 2, 2}, 4. }, {{1, 2, 2, 2}, 5. }, {{1, 2, 2, 2}, 
            5. }, {{1, 2, 2, 2}, 4. }, {{1, 3, 1, 1}, 0.5}, {{1, 3, 1, 1}, 
            0.5}, {{1, 3, 1, 1}, 0.5}, {{1, 3, 1, 1}, 1. }, {{1, 3, 1, 1}, 
            1. }, {{1, 3, 1, 2}, 1.5}, {{1, 3, 1, 2}, 1. }, {{1, 3, 1, 2}, 
            2. }, {{1, 3, 1, 2}, 2. }, {{1, 3, 1, 2}, 2. }, {{1, 3, 2, 1}, 
            0. }, {{1, 3, 2, 1}, 0. }, {{1, 3, 2, 1}, 0.5}, {{1, 3, 2, 1}, 
            0.5}, {{1, 3, 2, 1}, 0.5}, {{1, 3, 2, 2}, 1. }, {{1, 3, 2, 2}, 
            1. }, {{1, 3, 2, 2}, 1.5}, {{1, 3, 2, 2}, 1. }, {{1, 3, 2, 2}, 
            2. }, {{2, 1, 1, 1}, 3. }, {{2, 1, 1, 1}, 3. }, {{2, 1, 1, 1}, 
            2. }, {{2, 1, 1, 1}, 2. }, {{2, 1, 1, 1}, 4. }, {{2, 1, 1, 2}, 
            4. }, {{2, 1, 1, 2}, 5. }, {{2, 1, 1, 2}, 5. }, {{2, 1, 1, 2}, 
            4. }, {{2, 1, 1, 2}, 4. }, {{2, 1, 2, 1}, 3. }, {{2, 1, 2, 1}, 
            3. }, {{2, 1, 2, 1}, 3. }, {{2, 1, 2, 1}, 3. }, {{2, 1, 2, 1}, 
            4. }, {{2, 1, 2, 2}, 4. }, {{2, 1, 2, 2}, 5. }, {{2, 1, 2, 2}, 
            5. }, {{2, 1, 2, 2}, 5. }, {{2, 1, 2, 2}, 5. }, {{2, 2, 1, 1}, 
            2. }, {{2, 2, 1, 1}, 2. }, {{2, 2, 1, 1}, 3. }, {{2, 2, 1, 1}, 
            3. }, {{2, 2, 1, 1}, 4. }, {{2, 2, 1, 2}, 4. }, {{2, 2, 1, 2}, 
            5. }, {{2, 2, 1, 2}, 5. }, {{2, 2, 1, 2}, 5. }, {{2, 2, 1, 2}, 
            5. }, {{2, 2, 2, 1}, 4. }, {{2, 2, 2, 1}, 4. }, {{2, 2, 2, 1}, 
            5. }, {{2, 2, 2, 1}, 5. }, {{2, 2, 2, 1}, 4. }, {{2, 2, 2, 2}, 
            4. }, {{2, 2, 2, 2}, 5. }, {{2, 2, 2, 2}, 5. }, {{2, 2, 2, 2}, 
            4. }, {{2, 2, 2, 2}, 4. }, {{2, 3, 1, 1}, 3. }, {{2, 3, 1, 1}, 
            3. }, {{2, 3, 1, 1}, 4. }, {{2, 3, 1, 1}, 4. }, {{2, 3, 1, 1}, 
            4. }, {{2, 3, 1, 2}, 4. }, {{2, 3, 1, 2}, 5. }, {{2, 3, 1, 2}, 
            5. }, {{2, 3, 1, 2}, 4. }, {{2, 3, 1, 2}, 4. }, {{2, 3, 2, 1}, 
            3. }, {{2, 3, 2, 1}, 3. }, {{2, 3, 2, 1}, 3. }, {{2, 3, 2, 1}, 
            4. }, {{2, 3, 2, 1}, 4. }, {{2, 3, 2, 2}, 4. }, {{2, 3, 2, 2}, 
            5. }, {{2, 3, 2, 2}, 5. }, {{2, 3, 2, 2}, 6. }, {{2, 3, 2, 2}, 
            5. }};

dat2 is the same dependent variable data as dat1 but in grouped format with 3 cell indices;

dat2 = {{{1, 1, 1}, {0. , 1. , 1. , 2. , 2. , 2. , 2. , 3. , 3. , 
              4. }}, {{2, 3, 1}, {3. , 3. , 4. , 4. , 4. , 4. , 5. , 5. , 
              4. , 4. }}, {{1, 1, 2}, {0. ,  .5,  .5, 1. , 2. , 2. , 2.5, 
              2.5, 3. , 2. }}, {{1, 2, 1}, {2. , 1. , 1. , 2. , 3. , 4. , 
              4. , 4. , 4. , 5. }}, {{1, 2, 2}, {2. , 2. , 1. , 3. , 3. , 
              4. , 4. , 5. , 5. , 4. }}, {{1, 3, 1}, { .5,  .5,  .5, 1. , 
              1. , 1.5, 1. , 2. , 2. , 2. }}, {{1, 3, 2}, {0. , 
              0. ,  .5,  .5,  .5, 1. , 1. , 1.5, 1. , 2. }}, {{2, 1, 
              1}, {3. , 3. , 2. , 2. , 4. , 4. , 5. , 5. , 4. , 4. }}, {{2, 
              1, 2}, {3. , 3. , 3. , 3. , 4. , 4. , 5. , 5. , 5. , 
              5. }}, {{2, 2, 1}, {2. , 2. , 3. , 3. , 4. , 4. , 5. , 5. , 
              5. , 5. }}, {{2, 2, 2}, {4. , 4. , 5. , 5. , 4. , 4. , 5. , 
              5. , 4. , 4. }}, {{2, 3, 2}, {3. , 3. , 3. , 4. , 4. , 4. , 
              5. , 5. , 6. , 5. }}};

dat3 is for a 2-factor anova; data is ungrouped and in Integer format. N=33, unequal n run as {2x4} design.

dat3 = {{{1, 1}, 3}, {{1, 1}, 4}, {{1, 1}, 6}, {{1, 1}, 7}, {{1, 2}, 
            5}, {{1, 2}, 6}, {{1, 2}, 6}, {{1, 2}, 7}, {{1, 2}, 7}, {{1, 3}, 
            4}, {{1, 3}, 6}, {{1, 3}, 8}, {{1, 3}, 8}, {{1, 4}, 8}, {{1, 4}, 
            10}, {{1, 4}, 10}, {{1, 4}, 7}, {{1, 4}, 11}, {{2, 1}, 
            2}, {{2, 1}, 3}, {{2, 1}, 4}, {{2, 2}, 3}, {{2, 2}, 5}, {{2, 2}, 
            6}, {{2, 2}, 3}, {{2, 3}, 9}, {{2, 3}, 12}, {{2, 3}, 
            12}, {{2, 3}, 8}, {{2, 4}, 9}, {{2, 4}, 7}, {{2, 4}, 
            12}, {{2, 4}, 11}};

dat4 is for a 3-factor anova with 1 observation per cell; df(error)=0 unless an interaction term is dropped; data is ungrouped and in Integer form.N=60 in a {6x2x5} design

dat4 = {{{1, 1, 1}, 5}, {{1, 1, 2}, 5}, {{1, 1, 3}, 4}, {{1, 1, 4}, 
            4}, {{1, 1, 5}, 3}, {{1, 2, 1}, 5}, {{1, 2, 2}, 5}, {{1, 2, 3}, 
            4}, {{1, 2, 4}, 4}, {{1, 2, 5}, 4}, {{2, 1, 1}, 3}, {{2, 1, 2}, 
            4}, {{2, 1, 3}, 2}, {{2, 1, 4}, 3}, {{2, 1, 5}, 2}, {{2, 2, 1}, 
            5}, {{2, 2, 2}, 5}, {{2, 2, 3}, 2}, {{2, 2, 4}, 2}, {{2, 2, 5}, 
            4}, {{3, 1, 1}, 1}, {{3, 1, 2}, 4}, {{3, 1, 3}, 5}, {{3, 1, 4}, 
            2}, {{3, 1, 5}, 4}, {{3, 2, 1}, 2}, {{3, 2, 2}, 3}, {{3, 2, 3}, 
            5}, {{3, 2, 4}, 2}, {{3, 2, 5}, 3}, {{4, 1, 1}, 3}, {{4, 1, 2}, 
            3}, {{4, 1, 3}, 2}, {{4, 1, 4}, 2}, {{4, 1, 5}, 4}, {{4, 2, 1}, 
            4}, {{4, 2, 2}, 4}, {{4, 2, 3}, 2}, {{4, 2, 4}, 2}, {{4, 2, 5}, 
            2}, {{5, 1, 1}, 1}, {{5, 1, 2}, 4}, {{5, 1, 3}, 4}, {{5, 1, 4}, 
            2}, {{5, 1, 5}, 2}, {{5, 2, 1}, 2}, {{5, 2, 2}, 4}, {{5, 2, 3}, 
            4}, {{5, 2, 4}, 3}, {{5, 2, 5}, 5}, {{6, 1, 1}, 2}, {{6, 1, 2}, 
            2}, {{6, 1, 3}, 2}, {{6, 1, 4}, 1}, {{6, 1, 5}, 2}, {{6, 2, 1}, 
            1}, {{6, 2, 2}, 3}, {{6, 2, 3}, 1}, {{6, 2, 4}, 1}, {{6, 2, 5}, 
            2}};

dat5 is for a 1-factor anova having 3 groups of unequal size; data is grouped and in Integer form; N=21

dat5 = {{{1}, {3, 5, 2, 4, 8, 4, 3, 9}}, {{2}, {4, 4, 3, 8, 7, 4, 
              2}}, {{3}, {6, 7, 8, 6, 7, 9}}};

Note:Data sets for the repeated measures designs and nested designs, as well as the examples from Wolfram are included in cells at the end of the notebook.

Generating the NIST Simon-Lesage 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.  
The functions below can be used to generate the data sets Simon-Lesage 7, 8, and 9. 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

These data sets are provided by the National Institute of Standards & Technology (NIST), Information Technology Laboratory.The www site provides information on how extended real-value precision was obtained. For each data set, there are 9 cells with an equal number of observation in each cell.

Infinite Precision Results for Simon-Lesage Data Sets

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 the notebook 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)
Function getrn which brings Real input data into a rational form, takes about 45 seconds on a Macintosh 7200/120 when used with Simon 9.

The function gensimon generates the data set of Simon-Lesage7 used by NIST  for testing SS accuracy. The data set returned is in grouped format. By changing the iterator {10} to {100} or {1000} the Simon-Lesage data sets 8 (N=1809)  and 9 (N=18009), respectively,  can be generated.

\(\(gensimon := 
        Module[{d1, d2, d3, s1, s2}, 

          d1 = Flatten[Table[{1000000000000.3, 1000000000000.5}, {10}]]; 

          d1 = Prepend[d1, 1000000000000.4]; 

          d2 = Flatten[Table[{1000000000000.2, 1000000000000.4}, {10}]]; 

          d2 = Prepend[d2, 1000000000000.3]; 

          d3 = Flatten[Table[{1000000000000.4, 1000000000000.6}, {10}]]; 

          d3 = Prepend[d3, 1000000000000.5]; 

          s1 = {d1, d2, d3, d2, d3, d2, d3, d2, d3}; 

          s2 = Map[Append[{{#}}, s1[\([#]\)]] &, Range[9]]; 

          Return[s2];
		];\)\)

simon7 = gensimon;
Length[simon7]

The function genunsimon generates the Simon-Lesage7 data in an ungrouped format. By changing the iterator {10} to {100} or {1000} the Simon-Lesage data sets 8 and 9, respectively,  can be generated.

\(\(genunsimon := 
        Module[{d1, d1c, d2, d2c, d3, d3c, d4c, 
		d5c, d6c, d7c, d8c, 
            d9c, s7c}, 

          d1 = Flatten[Table[{1000000000000.3, 1000000000000.5}, {10}]]; 

          d1 = Prepend[d1, 1000000000000.4]; 	

          d1c = Map[Append[{{1}}, #] &, d1]; 

          d2 = Flatten[Table[{1000000000000.2, 1000000000000.4}, {10}]]; 

          d2 = Prepend[d2, 1000000000000.3]; 

          d2c = Map[Append[{{2}}, #] &, d2]; 	

          d3 = Flatten[Table[{1000000000000.4, 1000000000000.6}, {10}]]; 

          d3 = Prepend[d3, 1000000000000.5]; 

          d3c = Map[Append[{{3}}, #] &, d3]; 

          d4c = Map[Append[{{4}}, #] &, d2]; 

          d5c = Map[Append[{{5}}, #] &, d3]; 

          d6c = Map[Append[{{6}}, #] &, d2]; 

          d7c = Map[Append[{{7}}, #] &, d3]; 

          d8c = Map[Append[{{8}}, #] &, d2]; 

          d9c = Map[Append[{{9}}, #] &, d3]; 

          s7c = {d1c, d2c, d3c, d4c, d5c, d6c, d7c, d8c, d9c}; 

          s7c = Flatten[s7c, 1]; 

          Return[s7c];										
		];\)\)

simon7u = genunsimon;

Length[simon7u]

Rationalizing Input Data  to Obtain High Accuracy

The NIST data sets are designed to check the adequacy of computational algorithms and do not represent data likely to be encountered in a real experimental situation. Even if such a data set was encountered it is likely that the researcher would apply a linear transformation to the data, e.g., subtracting a constant, so that the data contained only the digits to the right of the decimal point.  Any linear transformation of the form y=a+bx with b≠0, will change the sum of squares and mean squares, however, the ratios (F tests) would not be affected. 

A high degree of accuracy is obtained when the experimental data is in integer or rational form. When the data is in real form, it can be rationalized using the function getrn called by the functions srtcnt and srtgrp. The function getrn is defined in two different forms, one form using RealDigits and FromDigits, and the other using only Rationalize[x,0]. The user can initialize either form of the function. The differences resulting in the use of Rationalize[x,0] and the function using RealDigits and FromDigits  are  illustrated  below. For purposes of illustration, the function getrn6 as given below is used. This latter function differs from getrn only in that getrn operates on a matrix of data, while getrn6 operates on a vector of data.

getrn6[d_] := Sign[d]*Map[FromDigits[#] &, RealDigits[d]]

Some large example values are given as the variable tdr shown below.

tdr = {100000000000000.3, 1000000000000000.3, 1000000000000000.30, 
        10000000000000000.3}

In the variable  tdr the first value has 14 zeroes, the second has 15 zeroes, the third value has 15 zeroes but an additional zero following the decimal digit 3, the fourth value has 16 zeroes.  The results of rationalizing these values using getrn6 are as follows:

getrn6[tdr]

N[%3, 20]

Notice that the second value was not rationalized to the same value as the third value. Examining the rationalization further,

$MachinePrecision

Map[Precision[#] &, tdr]

RealDigits[tdr]

When the value is less than or equal to machine precision, machine precision is used, thus the decimal portion of the second  value is lost. When the value is greater than machine precision, arbitrary precision is used. (See sections 3.1.3 and 3.1.4 of the Mathematica Book.).  If we use the Mathematica function Rationalize, the results are as follows:

Map[Rationalize[#, 0] &, tdr]

N[%, 20]

In the above results, only the third value is  without error, all others  contain an error in the representation of digits to the right of the decimal point.  The NIST data is of the order 1.000000000000n x (10^12)(where n is a single digit)  so no error is encountered using the built in function Rationalize[x,0] or the function getrn. (Note: in the functions that require initialization, two forms of getrn are given. One is of the form of getrn6 defined above, and the other uses the built-in Rationalize function. Only one of these needs to be intialized.)

A comparison similar to that given above is given below for very small values.

tdr2 = { .00000000000014,  .000000000000015,   .0000000000000016,  .00000000000000017}

Map[Precision[#] &, tdr2]

getrn6[tdr2]

N[%, 20]

Map[Rationalize[#, 0] &, tdr2]

N[%, 20]

If one assumes that typical research data is in the range 0-1000 and containing 3 decimal digits, the user could define the input data  in the following ways:

1) real values which will invoke the function getrn. The error in the numerator of the rationalized value may be 1x(10^(-16)) if the function getrn involves the operations Sign, RealDigits, and FromDigits, i.e., the first definition of getrn in the List of Functions to Initialize.

2) real  values, but redefine the function getrn to be getrn[d_]:=Rationalize[d,0], i.e., use the second definition of getrn in the List of Functions to Initialize. Examples tried showed no error in the numerator of the rationalized values.

3) define the data  as integers by multiplying by 1000. The calculated means will be increased  by (10^3) and sum of squares, mean squares, and variances by (10^6).  The F ratios will  be unaffected.

Function to Check Cell Indices

The function chkidx can be used to check the adequacy of cell indices required of the data used for an anova. The function checks (a) the length of the vectors holding the indices, (b) for a Null included as part of the indices, and (c) missing cell indices. Output includes number of cells found, minimum and maximum length of cell index  vectors, location of Null indices, cells indices  which are missing, and frequency of each cell index  vector.

chkidx[d] is used to check the adequacy of the cell index vectors in the data set d.
Calls: ndim (located in the next section)
Mod. Date: July 5/02

\(\(chkidx[d_] := 
        Module[{indc, ilen, minl, maxl, np, flg = 0, dim, nc, sum, psn0, 
            alli, fm, 
            fm2}, 
Print["Checking Cell Indices"]; 
indc = 
            Map[Part[#, 1] &, 
              d]; 
Print["Number of vectors of cell indices=", 
            Length[indc]]; 
 (*check\ length\ of\ indices*) 
ilen = 
            Map[Length[#] &, indc]; 
minl = 
            Min[ilen]; 
maxl = 
            Max[ilen]; 
If[minl\  ≠ \ maxl, 
            Print["**ERROR:Length of index vectors are not same"]; 
            flg = 1; 
Print["Min length=", 
              minl, " Max length=", maxl], \ 
            Print["Index vector lengths all equal"]]; 
 (*Check\ for\ Null\ in\ index\ vector*) 
np = 
            Position[indc, Null]; 
If[
            Length[np] ≠ \ 0, 
            flg = 1; 
Print["Position of Null=", 
              np]; 
Print["Correct indices having Null and execute function again"]]; 
If[flg == 1, \ 
            Return[]]; 
 (*check\ for\ missing\ cell\ indices*) 
Print["Checking for missing cell indices"]; 
dim = 
            Map[Max[#] &, 
              Transpose[
                indc]]; 
Print["Dimensions of data=", 
            dim]; 
nc = 
            Apply[Times, 
              dim]; 

 (*create\ sum\ box*) 
fm = Table[0, {nc}]; 
fm = 
            ndim[fm, 
              dim]; 
 (*count\ freq\ of\ indices*) 
Do[
sum = 
              1 + Extract[fm, indc[\([j]\)]]; 
fm = 
              ReplacePart[fm, sum, indc[\([j]\)]]
, {j, 1, 
              Length[indc]}]; 

psn0 = 
            Position[fm, 0]; 
Print["Missing indices=", 
            psn0]; 
fm2 = 
            ndim[Range[nc], dim]; 
alli = 
            Map[Position[fm2, #] &, Range[nc]]; 
alli = 
            Flatten[alli, 1]; 
fm = 
            Transpose[{Flatten[fm]}]; 
Print[
            TableForm[Transpose[{alli, fm}], TableDepth -> 2, 
              TableHeadings -> {Automatic, {"Index", \ " Frequency"}}]];

];\)\)

Data for this example is from dat2 with the cell indices {2,3,1} and {1,3,1} and their associated data removed.

{{2, 3, 1}, {3. , 3. , 4. , 4. , 4. , 4. , 5. , 5. , 4. , 
        4. }},

{{1, 3, 1}, { .5,  .5,  .5, 1. , 1. , 1.5, 1. , 2. , 2. , 
        2. }},

dat21 = {{{1, 1, 1}, {0. , 1. , 1. , 2. , 2. , 2. , 2. , 3. , 3. , 
              4. }}, {{1, 1, 2}, {0. ,  .5,  .5, 1. , 2. , 2. , 2.5, 2.5, 
              3. , 2. }}, {{1, 2, 1}, {2. , 1. , 1. , 2. , 3. , 4. , 4. , 
              4. , 4. , 5. }}, {{1, 2, 2}, {2. , 2. , 1. , 3. , 3. , 4. , 
              4. , 5. , 5. , 4. }}, {{1, 3, 2}, {0. , 0. ,  .5,  .5,  .5, 
              1. , 1. , 1.5, 1. , 2. }}, {{2, 1, 1}, {3. , 3. , 2. , 2. , 
              4. , 4. , 5. , 5. , 4. , 4. }}, {{2, 1, 2}, {3. , 3. , 3. , 
              3. , 4. , 4. , 5. , 5. , 5. , 5. }}, {{2, 2, 1}, {2. , 2. , 
              3. , 3. , 4. , 4. , 5. , 5. , 5. , 5. }}, {{2, 2, 2}, {4. , 
              4. , 5. , 5. , 4. , 4. , 5. , 5. , 4. , 4. }}, {{2, 3, 
              2}, {3. , 3. , 3. , 4. , 4. , 4. , 5. , 5. , 6. , 
              5. }}};

chkidx[dat21]

Example of Interactive User Input for Data Set dat2

dat2 is used in a straight-forward 3-way anova. User queries are in quotations.
1. anovanw[dat2]
2. "Type the number of factors in the design"   3
3. "Type vector of Factor indices to use" {1,2,3}
4.  "Display cell data as Real (1) or Rational (0)?"  1
5. "Do you wish to proceed? y/n"  y
6. "Remove interaction terms? y/n"  n
7. "Plot interactions? y/n"  y
   7.1 "Type an interaction term to plot"  ab
   7.2 "More plots? y/n"  y
   7.3 "Type an interaction term to plot"  abc
   7.4 "More plots? y/n"  n
8. "Contrasts on main effects? y/n"  y
   8.1 "Type contrast term"  b
   8.2 "How many rows in contrast?  2
   8.3 "Type contrast matrix"  {{1,0,-1},{0,1,-1}}
   8.4 "More contrasts"  n

Required Packages (Requires Initialization)

<< "Statistics`DataManipulation`"

<< "LinearAlgebra`MatrixManipulation`"

<< "Statistics`ContinuousDistributions`"

<< "Statistics`DescriptiveStatistics`"

List of Functions Used

anovanw: main anova starting function
anvT2:  displays general Anova table of SS; calculates MSe and dfe
binselect: selects elements in list corresponding to positions of 1 in binary vector
bin2alf: generates character representation from binary vector
chkidx: checks adequacy of cell index vector (not part of anova procedure)
comdta: to handle grouped data by combining data belonging to common cells
cutintr:  allows user to remove designated interaction vectors in design matrix
dmpar:  generates character and binary representation of all parameters
dsgnint: generates interaction vectors for design matrix
dsgnmef:  generates vectors for main effects design matrix
examd:   determines whether input data is all Real, all Integers, or mixed
findc:    locates cols of design matrix used to form interaction vectors
findintr:  locates 2 & 3 term interactions from binary designation
gencon:  generates contrast matrices and displays SS table for all terms
genrch:   generates row/col headings for pmat2
getrn:     gets rational values for real valued input data (2 definitions)
getss:    gets deviated sum of squares and means
getucn:  obtain user constrast matrix
hytst:  display hypothesis tested for user contrasts
intrplot:  entry to plotting of interaction terms
kbin3:   keyboard input of numerics
kbinstrg: keyboard input of characters
lctmev: locate starting/ending position of vectors representing anova terms in X'X
matsum:  calculates a matrix of  means for plotting interactions
ndim:    restructures vector of means to a matrix of given dimensions
plotnLists:  plots 2-term interactions
plt3fintr: plots planes of 3-term interactions
pmat2:    displays a matrix
prntimm: to display means matrix used to plot interaction
probF:  calculate probability of F
redspc:  sums over specific dimensions of means matrix
restrv:  restructures a vector
rndoff:   round off real values
srtcnt:    sort ungrouped data by factor indices; returns cell sizes, means, and variances
srtgp:    using grouped data, returns cell sizes, means, and variances
ssk:       calculates SS for contrast matrix K

Each function is preceded with documentation including other functions which are called.

Error Messages

Two types of error messages may appear.  On input of user parameters or options, an error message will be given when the structure of the input is inappropriate or not acceptable. In such cases the message will appear on the screen, but may be partially covered by the input window in which a re-entry is requested. The main function anovanw will Abort if it is presented with data which contains both Real and Integer values with a message presented by Message[] asking the user to present either all Real or all Integer data. The program will also abort if the data presented is of Length[]==0, usually caused by the user forgetting to have Mathematica evaluate the data set.

Functions to Initialize

kbin3[q_,p_] function to accept numeric input from keyboard posing question in q. If p is any integer then an integer is to be input;  if p is any real value then a real value is to be input; if p={n} then a n element vector is expected; if p={r,c} then a matrix with r rows and c columns is expected. If the input is inappropriate the question is posed again.  If the input is appropriate it is returned.
Mod. Date: Mar 24/96

kbin3[q_, p_] := 
        Module[{c, ans, vrp, vrp2}, c = 0; 
          vrp = " Input Error. Type a vector with # elements="; 
          vrp2 = " Input Error. Type a matrix of order "; 
          While[c == 0, ans = Input[q]; 
            Which[Head[ans] === 
                Symbol, Print[
                  ans, " Input Error; Characters not allowed"];, 
              Head[p] === Integer || 
                Head[p] === Real, If[
                  Head[p] =!= 
                    Head[ans], Print[
                      ans, " Input Error. Input must be ", Head[p]];, 
                  c = 1];, 
              Head[p] === List && 
                Length[p] == 1, If[
                  VectorQ[ans] == False, Print[ans, vrp <> ToString[p]];, 
                  If[Length[ans] ≠ 
                      First[p], Print[
                        ans, " Error: Vector requires # of elements=", 
                        First[p]];, c = 1]];, 
              Head[p] === List && 
                Length[p] == 2, If[
                  MatrixQ[ans] == 
                    False, Print[ans, vrp2 <> ToString[p]];, 
                  If[Dimensions[ans] ≠ 
                      p, Print[
                        ans, " Error; Matrix dimensions must be ", p];, 
                    c = 1]];];]; Return[ans];];

kbinstrg[q,p] used to input a string in response to question in q. p contains acceptable answers in the form of a list, eg., p={{ans1},{ans2},....}. If the input is appropriate it is returned. If the input is inappropriate the question is posed again.
Mod. Date: Mar 24/96

kbinstrg[q_, p_] := 
        Module[{c, ans}, c = 0; 
          While[c == 0, ans = InputString[q]; 
            Do[If[{ans} == p[[j]], 
                  c = 1; Break[]];, {j, 1, Length[p]}]; 
            If[c ≠ 1, 
              Print[ans, " Input Error. Type Again"]];]; 
          Return[ans];];

pmat2[m,pc,rch] prints a 2D matrix with either numbered rows and cols or with labelled rows and columns, in units of pc columns 
m     is the matrix to print
pc    is the number of columns to use
rch   is optional; if present it contains the row and column headings in the 
       following format: {{"a","b","c",...},{"A","B","C","D",....}} where
       Part 1 gives the row headings, and Part 2 gives the column headings. 
       Part 1 must have as many headings as rows of the matrix, and Part 2
       must have as many headings as the columns of the matrix.
       If the argument rch is not used, the row and column labels
       will be the integer sequence starting at 1, and the col digits for any 2nd
       part of the matrix printed will continue the digit sequence.
 Mod Date: Jan 27/96

\(\(pmat2[m_, pc_, rch___] := 
        Module[{m1, a, c, indx, col1, row}, m1 = Transpose[m]; 
          c = \(Dimensions[m1]\)[[1]]; 
          a = Floor[c\/pc]; indx = Join[Table[pc, {a}], {Mod[c, pc]}]; 
          If[0 == indx[[\(-1\)]], 
            indx = Drop[indx, \(-1\)]]; 
          If[Length[{rch}] == 0, 
            row = \((ToString[#1] &)\) /@ 
                Range[\(Dimensions[
                      m]\)[[1]]]; \(col1 = \((ToString[#1] &)\) /@ 
                  Range[\(Dimensions[
                        m]\)[[2]]];\), 
            row = rch[[1]]; 
            col1 = rch[[2]];]; 
          Do[Print[""]; 
            Print[TableForm[
                Transpose[
                  Take[m1, indx[[j]]]], 
                TableHeadings -> 
                  Append[{row}, 
                    Take[col1, 
                      indx[[j]]]], 
                TableSpacing -> {0, 4}]]; 
            m1 = Drop[m1, 
                indx[[j]]]; \(col1 = 
                Drop[col1, 
                  indx[[j]]];\), {j, 1, 
              Length[indx]}];];\)\)

genrch[rh,ch,rc] generates row and column headings for printing a 2D matrix
rh   is the string prefix for row headings
ch   is the string prefix for the column headings
rc   is the Dimensions of the matrix as a vector
i1   is the starting index for rows
i2   is the starting index for columns
The prefix for row and columns will be followed by digits starting at i1 and i2. i1 and i2 are optional, but the order of not including is from the right. When not present, these are set to 1. The function can be used to generate the 3rd argument for pmat2.
Mod. Date: Feb 14/97

genrch[rh_, ch_, rc_, i1_:  1, i2_:  1] := 
        Module[{a, b, ab, ri, ci}, 
          fn[a_, b_] := StringInsert[a, ToString[b], -1]; 
          a = Table[rh, {rc[[1]]}]; 
          ri = i1 - 1; ci = i2 - 1; 
          a = (fn[a[[#1]], 
                    ri + #1] &) /@ 
              Range[rc[[1]]]; 
          b = Table[ch, {rc[[2]]}]; 
          b = (fn[b[[#1]], 
                    ci + #1] &) /@ 
              Range[rc[[2]]]; 
          ab = Join[{a}, {b}]; Return[ab];];

rndoff[n,d] rounds of real n to d decimal digits

rndoff[n_, d_] := N[(10^(-d)) Round[n (10^d)]];

getrn[d] is a pure function which converts the data in the structure d containing all Real values and returns rational values or integers maintaining the sign of each value.
Returns: a data structure of the same form as the input d, but with  rational or Integer values.
Mod. Date: Nov 29/97

Note: the first or the second definition of getrn  requires initialization (not both). See the section entitled Rationalizing Input Data to Obtain High Accuracy.

getrn[d_] := 
        Map[Map[Sign[#]*FromDigits[RealDigits[#]] &, #] &, d];

getrn[d_] := Rationalize[d, 0];

examd[d]  determines whether the data structure d is all in the form of Integers, Real values, or a mix of Integers and Real values. 
Returns: 0 if all values in d are Integers; 1 if all are Real, and 2 if there is a mixture
Mod. Date: Jun 26/02: removed duplicate output of number of observations

\(\(examd[d_] := 
        Module[{ni, nt, fl}, 

          ni = Apply[Plus, 
              Map[IntegerQ[#] /. {True -> 1, False -> 0} &, 
                Flatten[d]]]; 	
			nt = 
            Length[Flatten[d]]; 
		Print["Number of observations=", 
            nt]; 	
	
          Which[
			ni == nt, Print["Data all integers."]; 
            fl = 0, 
			ni == 0, \ Print["Data all real."\ ]; 
            fl = 1, 
		\ \ ni < \ nt, 
            Print["Mixed data. Integers=", ni]; fl = 2]; 
		Return[
            fl];
		];\)\)

srtcnt[d,nf] using ungrouped data, requests user to input the indices to be used for grouping the dependent variable (y), sorts the data according to indices and size of y, gets cell frequencies; examd is called to determine the type of input data. If all input data is Real calls getrn to bring data to rational form. Calls getss to calculate means, variances, and sum of squares.
d:  data matrix; {{i11,i12,...i1m},d1},{i21,i22,...i2m},d2}....}, where the m factor indices
     for each observation are given by i and the associated dependent variable by d. 
nf:  number of factors to use
User input: a vector of nf elements of the factor indices to use, e.g., {1,2,3}, {2,1,3}, {1}
Output: number of observations; number of levels in each factor
Returns: {sst,ssw,ord,freq,mv,gm}
sst: total sum of squares deviated about grand mean
ssw: pooled sum of squares within group
ord: a vector of the number of levels in each factor, 
freq: the number of observations in each cell ordered with last subscript varying the fastest, 
mv: the cell means ordered with last subscript varying the fastest
gm: total group mean
Note: A cell which has zero observations is not directly identified. 
Calls: kbin3,getss,restrv,examd, getrn
Required: Statistics`DataManipulation`
Mod. Date: Dec 1/97

\(\(srtcnt[d_, nf_] := 
        Module[{nob, dn, si, obgrp, mv, frq, gm, ord, sst, ssw, typd, 
            inx}, 

nob = 
            Length[d];  (*\(Print["Number of Observations=", nob];\)*) 
            si = kbin3["Type vector of Factor indices to use", {nf}]; 
          Print["Selected Factor Indices=", si]; 
          dn = Sort[\((ReplacePart[
                      d[[#1]], \(\(d[[#1]]\)[[1]]\)[[si]], 1] &)\) /@ 
                Range[nob]]; 
          obgrp = Frequencies[\((\(dn[[#1]]\)[[1]] &)\) /@ 
                Range[nob]]; 
          obgrp = \((Reverse[
                    obgrp[[#1]]] &)\) /@ 
              Range[Length[obgrp]]; 
          obgrp = Partition[Flatten[obgrp], 
              nf + 1]; 
		ord = \((Max[\(Transpose[
                        obgrp]\)[[#1]]] &)\) /@ Range[nf + 1]; ord = Take[ord, nf]; 
          Print["Number of Levels in each factor=", ord]; 
		dn = 
            Flatten[Rest[Transpose[dn]]]; 
		frq = 
            Flatten[Take[Transpose[obgrp], \(-1\)]]; 
			inx = 
            Map[Take[obgrp[\([#]\)], nf] &, 
              Range[Length[
                  obgrp]]]; 
		 (*next\ restructure\ data\ into\ cell\ units\ of\ length\ frq*) 
		dn = restrv[dn, frq]; 
		typd = 
            examd[dn]; \  (*check\ type\ of\ data*) 

		Which[
			typd == 0, {sst, ssw, mv, gm} = 
              getss[dn, inx, frq], 
			typd == 1, {sst, ssw, mv, gm} = 
              getss[getrn[dn], inx, frq], 
			typd == 2, 
            Message[srtcnt::mixed]; 
            Abort[];
			]; 
			Return[{sst, ssw, ord, frq, mv, 
              gm}];];\)\)

srtcnt::mixed = "The input data has Real and Integer values. \
Correct input to have values of one type only.";

srtgp[dg,nf] takes data grouped by cell index, and for nf factors and a user selection of nf cell
indices, combines the data of cells having a common set of indices, then calls getss and calculates the following:
1) cell means, SSdeviated about cell mean  and variances, cell sample sizes 
2) total group mean and SS deviated about grand mean,
3) SSdeviated about cell mean is summed to form SSwithin group
dg: data matrix in the form of {{{cell indices},{Y for each observation in the cell}},.....}
nf: number of factors in the design
Before calculating SS examd is called to assess the type of data input in order to achieve maximum accuracy in the results. If all the data is in Real Form, getrn is called to convert the data to rational form.
User Input: a vector of cell indices to be used, eg., there may be 3 cell indices, but user
               may select only two of them {1,3} collapsing the data over index 2.
Returns:{sst,ssw, orf,obgrp,mv,gm}, where 
a) sst is SStotal group deviated about grand mean,
b) ssw SSwithin group
c) orf is a vector of the number of levels in each factor
 d) obgrp is a vector of the cell frequencies, 
 e) mv a vector of the cell means,
 f) gm is the grand mean
Calls: kbin3,comdta, examd,getss,getrn
Mod. Date: Dec 1/97

\(\(\(srtgp[dg_, nf_] := 
        Module[{si, dgri, sd, orf, obgrp, ind, typd, sst, ssw, mv, 
            gm}, 

		si = 
            kbin3["Type vector of Factor indices to use", {nf}]; 
          Print["Selected Factor Indices=", 
            si]; 
		 (*get\ the\ selected\ indices\ in\ place*) 
		dgri = \((ReplacePart[
                    dg[[#1]], \(\(dg[[#1]]\)[[1]]\)[[si]], 1] &)\) /@ 
              Range[Length[dg]]; 
          sd = Sort[dgri]; 
		orf = \ Map[Part[#, 1] &, sd]; 
		orf = 
            Map[Max[#] &, 
              Transpose[
                orf]]; 
		Print["Number of Levels in each factor=", 
            orf]; 
		 (*combine\ cells\ having\ common\ indices*) 
		{ind, sd} = 
            comdta[sd]; 
		obgrp = \((Length[#1] &)\) /@ 
              sd;  (*cell\ sizes*) 
		typd = 
            examd[sd];  (*check\ type\ of\ data*) 
		Which[
		typd == 
              0, {sst, ssw, mv, gm} = getss[sd, ind, obgrp], 
		typd == 
              1, \  (*data\ input\ all\ Real\ so\ convert*) 	{sst, ssw, mv, 
                gm} = getss[getrn[sd], ind, obgrp], 
		typd == 2, 
            Message[srtgp::mixed]; \ 
            Abort[];
		]; 
		Return[{sst, ssw, orf, obgrp, mv, 
              gm}];
		];\)\( (*end\ module*) \)\)\)

srtgp::mixed = "The input data has Real and Integer values. Correct \
input to have values of one type only.";

getss[sd,ind,obgrp] using the data in sd, indices ind, and cell frequencies obgrp, displays a table of cell indices, cell  frequencies, cell means, and cell variances. The cell variances are ML estimates.  SS deviated about the group mean and SS for each group deviated about the group mean are calculated by multiplying by n the variances provided by the function VarianceMLE[] in the package Statistics`DescrptiveStatistics`.
User Input: user is asked whether cell means and variances should be output in real or rational form.
Returns: {sst,ssw,mv,gm} where sst is SStotal deviated, ssw is SSwithin group deviated, mv a vector of group means, and gm the grand mean.
Calls: kbin3
Requires: Statistics`DescriptiveStatistics`
Mod. Date: June 25/02; Modified output

\(\(getss[sd_, ind_, obgrp_] := 
        Module[{m, v, ror, vt, stuf, rch, sst, gm, ssw, fm}, 
          m = Map[Mean[#] &, sd]; \  (*grp\ means*) 
v = 
            Map[VarianceMLE[#] &, sd]; 
ror = 
            kbin3["Display cell data as Real (1) or Rational (0)", 
              1]; 
If[
            ror == 1, 
stuf = 
              Transpose[{ind, obgrp, N[m, 10], 
                  N[v, 10]}], 
stuf = 
              Transpose[{ind, obgrp, m, 
                  v}]]; 

Print[
            TableForm[stuf, TableDepth -> 2, 
              TableHeadings -> {Automatic, {"Cell Index", "Freq", "Mean", "Variance(MLE)"}}]]; 

 (*grand\ m\ and\ sstcfm*) 
gm = 
            Mean[Flatten[sd]]; 
vt = 
            VarianceMLE[
              Flatten[sd]]; 
Print["Total Sample Mean=", gm, "(=", 
            N[gm, 16], ")"]; 
Print["Total Sample Variance(MLE)=", vt, "(=", 
            N[vt, 16], ")"]; 
sst = 
            Apply[Plus, obgrp]*
              vt; \  (*get\ ss\ only*) 
ssw = 
            obgrp*v; \  (*ss\ by\ grp*) 
ssw = 
            Apply[Plus, 
              ssw]; \  (*sum\ within\ grp\ ss*) 
Print["SStotalCFM=", sst, "(=", 
            N[sst, 16], ")"]; 
Print["SSwithin=", 
            ssw, "(=", N[ssw, 16], ")"]; 
fm = 
            sst - ssw; 
Print["Full model SSmcfm=", 
            fm, " (=", N[fm, 16], ")"]; 
fm = 
            fm/sst; 
Print["Full model R-sqd=", 
            fm, " (=", 
            N[fm, 16], ")"]; 
 (*return\ ssotal, 
            sswithin, cell\ means, 
            and\ grand\ mean*) 
Return[{sst, ssw, m, 
              gm}];];\)\)

comdta[d] combines data (d) of cells having the same indices; finds unique cell indices
d: data matrix in the form of {{{indices},{data vector}},....}. Each row of data is 
expected to be in order sorted on {indices}.
Calls: binselect
Returns: {{unique indices}, {data for each cell}}
Mod. Date: Nov 27/97 (corrected to properly handle unequal n case)

\(\(comdta[d_] := 
        Module[{di, dd, newc, bs, dn, sd, inx}, 
          di = Flatten[Take[Transpose[d], 1], 
              1];  (*indices*) 
		dd = \(Transpose[
                d]\)[[2]];  (*cell\ data\ as\ row\ vectors*) 
		newc = {}; sd = {}; 
          While[Length[di] > 0, 
            bs = \((di[[#1]] == 
                        di[[1]] &)\) /@ 
                  Range[Length[di]] /. {True -> 1, 
                  False -> 0}; 
			inx = 
              Flatten[binselect[bs, {Range[Length[bs]]}]]; 
			sd = 
              Append[sd, dd[\([inx]\)]];  (*combine\ cells\ data*) dn = 
              Plus @@ bs; AppendTo[newc, Take[di, 1]]; di = Drop[di, dn]; 
            dd = Drop[dd, dn];]; newc = Flatten[newc, 1]; 
          sd = \((Flatten[#1] &)\) /@ sd; 
		Return[{newc, 
              sd}];];\)\)

anvT2[stc,sse,df] displays a general anova table containing SStcfm, SSmcfm, SSerror, df, MS, and F calculated for SSmcfm only. If df(error) is 0, MS(error) is set to 0, and no F tests are carried out.
stc: SS total cfm (SStcfm)
sse: SS error
smd: SS model cfm (SSmcfm); calculated by SStcfm-SSerror
df: vector of degrees of freedom for SStcfm and SSmcfm
Probability of F is rounded to 4 decimal digits for display; P(F)=0 is <.0001
Returns: {MSerror,dferror}
Calls: rndoff
Packages:  Statistics`ContinuousDistributions`
Note: SStcfm is SS total corrected for the mean; SSmcfm is SS model corrected for the mean. To
display SS and MS as rational values, remove the operator N[....]
Mod. Date: Nov 26/97

\(\(anvT2[stc_, sse_, df_] := 
        Module[{ct, dft, dfmd, dfe, mser, ssmcfm, fr, fd, p}, 
          Print["
General Anova Table"]; ct = {}; 
          dft = df[[1]]; 
          dfmd = df[[2]]; 
          dfe = dft - dfmd; 
          AppendTo[
            ct, {"SStcfm", N[stc, 6], dft - 1. , 
              N[stc\/\(dft - 1\), 6]}]; ssmcfm = stc - sse; 
          If[dfe > 0, mser = N[sse\/dfe]; 
            fr = ssmcfm\/\(\((dfmd - 1)\)\ mser\); 
            fd = FRatioDistribution[dfmd - 1, dft - dfmd]; 
            p = rndoff[1.  - CDF[fd, fr], 4], fr = "*"; p = "*"; 
            mser = 0]; 
		AppendTo[
            ct, {"SSmcfm", N[ssmcfm, 6], dfmd - 1, 
              N[ssmcfm\/\(dfmd - 1\), 6], fr, p}]; 
          If[dfe > 0, 
            AppendTo[ct, {"SSerror", rndoff[sse, 6], dfe, N[mser, 6]}], 
            AppendTo[ct, {"SSerror", Chop[sse], 0, "*"}]]; 
          Print[TableForm[ct, 
              TableHeadings -> {Automatic, {"Source", "SS", "df", "MS", "F", "P"}}, TableSpacing -> {0, 2}]]; 
          Return[{mser, dfe}];];\)\)

probF[ss,dfs,mse,dfe] calculates probability of F. ss is the sum of squares for the hypothesis being tested, dfs is the degree of freedom for ss, mse is the mean square error term having dfe degrees of freedom.
Calls: rndoff
Requires: <<Statistics`ContinuousDistributions`
Mod. Date: July 11/97

\(\(probF[ss_, dfs_, mse_, dfe_] := 
        Module[{ms, f, fd, p, ct}, ms = ss\/dfs; f = ms\/mse; 
          fd = FRatioDistribution[dfs, dfe]; p = rndoff[1 - CDF[fd, f], 4]; 
          ct = {{"Usr K", N[ss, 6], dfs, N[ms, 6], f, p}}; 
          Print[TableForm[ct, 
              TableHeadings -> {Automatic, {"Source", "SS", "df", "MS", "F", "P"}}, 
              TableSpacing -> {0, 2}]];];\)\)

dsgnmef[nl] creates a design matrix for main effects only
nl: vector of number of levels in each factor
Algorithm: Over each of j elements of nl,
a) a matrix having nl[j]+1 cols is created having nl[j] cols of 0s and the last column of ones
b) indices (the number of shifts required) for left shifting elements of the last column of ones into the leading cols are generated and replicated as required. Using these indices
c) the left shifts are mapped over the matrix and the last column which now has 0s is dropped
d) the sum reduction adjustment is done by subtracting the last column from the leading columns, and the last column  dropped.
Returns: main effects design matrix in the form of{{{a1},{a2},...}},{{b1},{b2}...},...}
Mod. Date: Dec 8/97

dsgnmef[ln_] := 
        Module[{kln, dm, e1, e2, e3, r1, r2, r3, r4}, kln = Join[{1}, ln]; 
          dm = {}; 
          Do[e1 = Table[
                Join[Table[
                    0, {ln[[
                        j]]}], {1}], {Times @@ ln}]; 
            e2 = Transpose[
                Table[Reverse[
                    Range[ln[[
                        j]]]], {Times @@ Drop[ln, j]}]]; 
            e3 = Flatten[Table[Flatten[e2], {Times @@ Take[kln, j]}]]; 
            r1 = (RotateLeft[
                      e1[[#1]], #1] &) /@ e3; r1 = Drop[Transpose[r1], -1]; 
            r2 = Take[r1, ln[[j]] - 1]; 
            r3 = Flatten[Take[r1, -1]]; 
            r4 = MapThread[Flatten[#1] - r3 &, {r2}]; AppendTo[dm, 
                r4];, {j, 1, Length[ln]}]; Return[dm];];

bin2alf[x] returns the character in alph corresponding to the position of a 1 in x, e.g., bin2alf[{1,1,1}] ->"abc" Currently set for characters "a" to "g" which sets the maximum number of levels in the anova to 7. To increase the dimensionality of the design, simply place the additional characters in the vector, e.g., {"a",......"g","h","i"} for a 9-way anova.
Mod. Date: June 7/97

\(\(bin2alf[x_] := 
        Module[{res, alph}, res = ""; 
          alph = Reverse[
              Take[{"a", "b", "c", "d", "e", "f", "g"}, Length[x]]]; \((If[
                  1 == x[[#1]], \(res = 
                      StringInsert[res, 
                        alph[[#1]], \(-1\)];\)] &)\) /@ Reverse[Range[Length[x]]]; Return[res];];\)\)

dmpar[nf] returns the character and binary designation  for all possible parameters in the anova design having nf factors. 
Calls: bin2alf to get character designations.
Returns: { {character, {binary representation}},......} in increasing order of the binary representation
Mod. Date: June 7/97

dmpar[nf_] := 
        Module[{bin, bin2, bin3, bin4}, 
          bin = IntegerDigits[Range[-1 + (2^nf)], 2]; 
          bin2 = MapThread[
              Join[Table[
                    0, {nf - 
                        Length[
                          bin[[#1]]]}], 
                  bin[[#1]]] &, {Range[-1 + (2^nf)]}]; 
          bin3 = MapThread[
              bin2alf[bin2[[#1]]] &, {Range[-1 + (2^nf)]}]; bin4 = Transpose[Append[{bin3}, bin2]]; 
          Return[bin4];];

dsgnint[apc,avc,cdm] generates interaction vectors
apc: all possible component terms for the design {{a,{0,0,1}},{b,{0,1,0}},....} in binary order
avc: components available; {{a,{1,0,0},1},{b,{0,1,0},2},{c,{1,0,0},3}...}; on entry only the
      main effect components are available
cdm: main effects design matrix; {{{column vectors of component a}},
                                                    {{column vectors of component b}}, ....}
Algorithm: 
1) for each binary vector in apc an interaction term is sought, i.e., more than one 1
2) for the term found in (1), and doing backward search over binary vectors  in avc, subtract
    each avc binary term until the result is a binary vector with only one 1; the term last subtracted
    and the term resulting from the subtraction identify the vectors of cdm which are to be
    multiplied to form the interaction vectors; multiply the vectors
3) the interaction vectors of  (2) are appended to cdm2 and avc2 is updated indicating which
    vectors are now available to form other interaction vectors
Calls: findc to find location of vectors to form product for interaction vectors. fvp is a function internal to dsgnint which multiplies the vectors element by element.
Returns: 
avc2: all components for complete design matrix in the form of
          {{a,{0,0,1},1},{b,{0,1,0},2},{c,{1,0,0},3}...{abc,{1,1,1},7}} giving the character
          representation, the corresponding binary vector, and the location in the design matrix 	of  the vectors. For example, Part[cdm2,7] contains all vectors for abc.
cdm2: the design matrix for all components in the form of
          {{{vectors for a}},{{vectors for b}}, ...{{vectors for abc}}}
Mod. Date: July 16/97

dsgnint[apc_, avc_, cdm_] := 
        Module[{avc2, cdm2, bin1, dif, lcn1, lcn2, prd, mat1, mat2, newterm}, 
          fvp[a_, b_] := (a b[[#1]] &) /@ Range[Length[b]]; avc2 = avc; cdm2 = cdm; 
          Do[If[1 ≠ 
                  Plus @@ apc[[
                        j]][[2]], 
                bin1 = apc[[
                      j]][[2]]; 
                Do[dif = 
                    bin1 - avc2[[
                          k]][[2]]; If[1 == Plus @@ (dif^2), 
                      lcn1 = avc2[[k]]; 
                      lcn2 = findc[dif, avc2]; 
                      mat1 = cdm2[[
                          lcn1[[3]]]]; 
                      mat2 = cdm2[[
                          lcn2[[3]]]]; 
                      prd = (fvp[
                                mat1[[#1]], mat2] &) /@ Range[Length[mat1]]; 
                      prd = Flatten[prd, 1]; AppendTo[cdm2, prd]; 
                      newterm = bin2alf[bin1]; 
                      AppendTo[avc2, {newterm, bin1, Length[cdm2]}]; 
                      Break[];];, {k, Length[avc2], 1, -1}];];, {j, 
              1, Length[apc]}]; Return[{avc2, cdm2}];];

findc[c,ac] find the vector c in the row vectors of ac; used to locate which col. vectors of 
the design matrix are required to generate an interaction term.
ac: of the form {{"a",{0,0,1},1},{"b",{0,1,0},2},...} and holds currently available vectors
    which can be used to generate interaction vectors
c:  binary vector to be found in rows of ac, eg., {0,1,0}
Returns: a row vector of ac which matches, e.g., {b,{0,1,0},2}
Mod. Date: July 16/97

findc[c_, ac_] := 
        Module[{res}, 
          Do[If[c == ac[[
                      j]][[2]], res = ac[[j]]; 
                Break[];];, {j, 1, Length[ac]}]; 
          Return[res];];

cutintr[x,it,nf] allows user to remove selected interaction vectors from design matrix
x:  design matrix in the form of {{{a1},{a2}},{{b1},{b2},{b3}},....}
it:  available terms in the form of {{"a",{0,0,1},1},....{"abc",{1,1,1},9}}
nf: number of factors in the design
User Input: interaction term to be removed, e.g., abc
Calls: kbinstrg,binselect
Returns: {new design matrix, available terms updated}
Mod. Date: July 21/97

\(\(cutintr[x_, it_, nf_] := 
        Module[{xc, itc, can, bs, flg = 0, strm, indx}, 
          can = Transpose[{Drop[Flatten[Take[Transpose[it], 1]], nf]}]; 
          bs = Table[0, {Length[it]}]; 
          While[flg == 0, 
            strm = kbinstrg["Interaction term to be drop?", can]; 
            Print["Interaction term dropped is:", strm]; 
            bs = bs + \((strm == \(it[[#1]]\)[[1]] &)\) /@ 
                    Range[Length[it]] /. {False -> 0, True -> 1}; 
            If["n" == 
                kbinstrg["Any more terms? y/n", {{"y"}, {"n"}}], 
              flg = 1];]; 
          bs = \((bs[[#1]] === 0 &)\) /@ 
                Range[Length[bs]] /. {False -> 1, True -> 0}; 
          bs = Table[1, {Length[it]}] - bs; 
          indx = binselect[bs, {Range[Length[it]]}]; indx = Flatten[indx]; 
          xc = x[[indx]]; 
          itc = it[[indx]]; 
          itc = \((ReplacePart[
                    itc[[#1]], #1, \(-1\)] &)\) /@ Range[Length[itc]]; Return[{xc, itc}];];\)\)

ssk[k,xi,pb] Function ssk calculates sum of squares for contrast
matrix k, with xi as Inverse[X'X], and pb the parameter estimates
Mod. Date: June 21/97

ssk[k_, xi_, pb_] := 
      Flatten[{k . pb} . Inverse[k . xi . Transpose[k]] . 
          Transpose[{k . pb}]]

ndim[vc,di] restructures vector vc into matrix of order given by vector di. Used to restructure
cell means into a matrix having dimensions equal to the number of levels in each factor.
Returns: matrix of Dimensions[]={di}
Mod. Date: June 27/97

ndim[vc_, di_] := 
        Module[{res, indx}, res = vc; indx = Drop[Reverse[di], -1]; 
          Do[res = 
              Partition[res, 
                indx[[j]]], {j, 1, 
              Length[indx]}]; Return[res]];

fndintr[sl]; finds and returns all 2 & 3 term interactions in list sl
sl: list of all available anova terms; {"mu",{0,0,0},1},{"a",{0,0,1},2}...}
Returns: interaction term of the form  {"ab",{0,1,1},5}; if no interaction terms are found {}is returned
Mod. Date: July 24/97

fndintr[sl_] := 
        Module[{sl2, cr}, sl2 = {}; 
          Do[cr = Plus @@ sl[[
                    j]][[2]]; If[
                2 == cr || 
                  3 == cr, AppendTo[sl2, 
                    sl[[j]]];];, {j, 
              1, Length[sl]}]; Return[sl2];];

plotnLists[m,tl,xu,ps,pm,pn] plots a colored line for each column of m.  The values in the column are y coordinates and the x coordinates are taken to be from 1 to the number of rows of the matrix.
tl is a title for the plot, 
xu is the label for the x-axis
ps  set to 1 for individual plots for each col; 0 for no individual plots
pm set to 1 for one plot holding multiple lines; 0 for no plot
pn  (Optional) set to 1 for line numbers at end of lines when pm=1; else 0 
For single plots (ps=1), the title is appended to a plot sequence number.
 The colors are rotated through the following color sequence:(1) black, (2) blue (3) green, (4) cyan, (5) red, (6) magenta, (7) yellow, (8) dark blue, (9) dark green (10) turquoise, (11)brown, and (12) purple. (Colors numbered 8 to 12 may vary with different monitors.) To get a simple example of the colors used, do plotnLists[Table[Range[0,11],{2}],"Legend","units",0,1,1] 
Mod. Date: July 11/97

\(\(plotnLists[m_, tl_, xu_, ps_, pm_, pn_:  0] := 
        Module[{pobj, colr, ipt, cntr = 0}, 
          If[ps == 0 && pm == 0, Print["Parameter error;no plot"]; 
            Return[]]; pobj = {}; 
          colr = { .5,  .0,  .5,  .0,  .0,  .0,  .0,  .0, 1. ,  .0, 
              1. ,  .0,  .0, 1. , 1. , 1. ,  .0,  .0, 1. ,  .0, 1. , 1. , 
              1. ,  .0,  .0,  .0,  .5,  .0,  .5,  .0,  .0,  .5,  .5,  .5,  .0,  .0}; If[ps == 1, ipt = $DisplayFunction, ipt = Identity]; 
          pobj = MapThread[
              ListPlot[#1, AxesLabel -> {xu, "Y"}, 
                  PlotJoined -> True, PlotRange -> All, 
                  PlotLabel -> 
                    StringInsert[tl, ToString[cntr = cntr + 1], 1], 
                  PlotStyle -> {RGBColor[\((colr = 
                              RotateLeft[colr, 
                                3])\)[[1]], colr[[2]], 
                        colr[[3]]]}, 
                  DisplayFunction -> ipt] &, {Transpose[m]}]; 
          If[pn == 1, \(MapThread[
                AppendTo[pobj[[#1, 1]], 
                    Text[#1, 
                      Last[pobj[[#1, 1, 2, 
                          1]]]]] &, {Range[
                    Length[Transpose[m]]]}];\)]; 
          If[pm == 1, 
            Show[pobj, PlotLabel -> tl, 
              DisplayFunction -> $DisplayFunction]];];\)\)

binselect[b_,a_]  columns of the matrix a corresponding to 1 in b are selected and returned as a new matrix. b must have as many elements as columns of a. If a is a vector, enclose it in {} and then Flatten the result returned. Equivalent to APL reduction operation.

binselect[b_, a_] := 
        Module[{indx, z}, 
          indx = DeleteCases[
              b Range[Dimensions[
                      a][[2]]], 0]; 
          z = Transpose[Transpose[a][[
                indx]]]; Return[z];];

matsum[nlv,ptr,m] forms means over different dimensions of matrix m specified indirectly by the binary terms in ptr; used to calculate means for plotting of interactions
nlv: number of levels in each factor
ptr: interaction term; {"ab",{0,1,1},psn}
m: matrix of means having dimensions defined by those terms in nlv 
The summation is made over those dimensions having 0 in the corresponding position of binary vector ptr.
Returns: a matrix of means
Calls: binselect, redspc
Mod. Date: July 11/97

matsum[nlv_, ptr_, m_] := 
        Module[{inx, dv, ms}, 
          inx = binselect[
              Table[1, {Length[nlv]}] - 
                ptr[[2]], {Reverse[
                  nlv]}]; dv = Times @@ Flatten[inx]; 
          Print["Number of cells summed over=", dv]; If[0 == dv, dv = 1]; 
          ms = redspc[m, ptr[[2]]]; 
          ms = N[(ms/dv)]; Return[ms];];

redspc[m,bv] reduces the dimensionality of the matrix m by summing over dimensions 
identified by the position of 0 in the binary vector bv; used to find the sums (and later the means) for plotting 2d and 3d interaction terms. 
Note: The binary vector bv low-order bit represents "a"; the next "b", etc.
Mod. Date: Jun 30/97

redspc[m_, bv_] := 
        Module[{rm, spr = 1, c}, rm = m; 
          c = Reverse[Table[1, {Length[bv]}] - bv] Range[Length[bv]]; 
          Do[If[0 ≠ c[[j]], 
                rm = Apply[Plus, 
                    rm, {c[[j]] - spr}]; 
                spr = spr + 1];, {j, 1, Length[c]}]; 
          Return[rm];];

plt3fintr[mnm,int] plots 3D matrix of means;each plane represents 1st term of the interaction
mnm: matrix of means
int: interaction term, e.g., "abc"
Mod. Date: July 11/97

plt3fintr[mnm_, int_] := 
        Module[{pobj}, 
          pobj = MapThread[
              ListPlot3D[#1, DisplayFunction -> Identity] &, {mnm}]; 
          Show[pobj, 
            PlotLabel -> StringInsert["Interaction: ", int, -1], 
            AxesLabel -> 
              Join[Take[Reverse[Characters[int]], 2], {"Y"}], 
            DisplayFunction -> $DisplayFunction];];

intrplot[m,st,nl] plots 2 and 3 term interactions based on cell means
m: vector of cell means
st: anova terms available; {{"mu",{0,0,0},1},{"a",{0,0,1},2},.....}
nl: vector of the number of levels in each factor
Steps:
1) fndintr finds the interaction terms, and allows user selection of 2 & 3 term interactions and
     returns {"ab",{binary vector representation}, posn}
2) ndim restructures m (means) to dimensions given by nl (number of levels in each factor)    
3) matsum uses the results returned by (1) and sums over the appropriate
     dimensions of m to get the means for plotting
4) plot 2 term interaction as line plots; 3 term interaction as planes
Calls: ndim,findintr,kbinstrg,matsum,plotnLists,plt3fintr,prntimm
Mod. Date: Dec 8/97

\(\(intrplot[m_, st_, nl_] := 
        Module[{mm, trm, ans, rp, pltrm, msrd, flag = 0}, trm = fndintr[st]; 
          If[0 == Length[trm], Print["No interaction terms"]; Return[]]; 
          mm = ndim[m, nl]; 
          ans = Transpose[{\((#1[[1]] &)\) /@ trm}]; Print["
Interaction terms available=", ans]; 
          While[flag == 0, 
            rp = kbinstrg["type an interaction term to plot", ans]; 
            Print["
Plot for interaction=", rp]; 
            Do[\(If[{rp} == ans[[j]], 
                  pltrm = trm[[j]]; 
                  Break[]];\), {j, 1, Length[ans]}]; 
            If[StringLength[rp] == 2, 
              If[Length[nl] == 2, prntimm[mm, 6, rp]; 
                plotnLists[Transpose[mm], 
                  StringInsert["Interaction:", 
                    rp, \(-1\)], \(Reverse[
                      Characters[
                        rp]]\)[[1]], 0, 
                  1, 1], msrd = matsum[nl, pltrm, mm]; prntimm[msrd, 6, rp]; 
                plotnLists[Transpose[msrd], 
                  StringInsert["Interaction:", 
                    rp, \(-1\)], \(Reverse[
                      Characters[
                        rp]]\)[[1]], 0, 
                  1, 1]], \(If[Length[nl] == 3, prntimm[mm, 6, rp]; 
                  plt3fintr[mm, rp], msrd = matsum[nl, pltrm, mm]; 
                  prntimm[msrd, 6, rp]; plt3fintr[msrd, rp]];\)]; 
            If["n" == 
                kbinstrg["More plots? y/n", {{"y"}, {"n"}}], 
              flag = 1];];];\)\)

prntimm[m,pc,is] displays 2d or 3d interaction means matrix
m: a 2d or 3d matrix of means
pc: max number of columns to use in display
is: 2 or 3 character anova interaction term; eg., "ac", "abd"
Calls: pmat2, genrch
Mod. Date: Dec. 2/97

prntimm[m_, pc_, is_] := 
        Module[{lv, rwch}, 
          Print[StringInsert["Means Matrix for Interaction: ", 
              is, -1]]; 
          If[2 == Length[Dimensions[m]], 
            rwch = genrch[StringTake[is, {1}], StringTake[is, {2}], 
                Dimensions[m]]; pmat2[N[m], pc, rwch], 
            rwch = genrch[StringTake[is, {2}], StringTake[is, {3}], 
                Rest[Dimensions[m]]]; 
            Do[lv = StringInsert[StringTake[is, 1], ToString[j], -1]; 
              Print["\nLevel=", lv]; 
              pmat2[N[m[[j]]], pc, 
                rwch], {j, 
                1, Dimensions[
                    m][[1]]}];];];

restrv[v,p] restructures elements in vector v into units of length p
v: a vector to restructure
p: length of each unit 
Returns: new structure of Length=Length[p]
Mod. Date: July 11/97

restrv[v_, p_] := 
        Module[{rv, temp}, temp = v; rv = {}; 
          Do[AppendTo[rv, 
              Take[temp, p[[j]]]]; 
            temp = RotateLeft[temp, 
                p[[j]]], {j, 1, 
              Length[p]}]; Return[rv];];

getucn[at,df,xi,inxi,ems,edf,mefb] obtains a user contrast matrix in non-sum reduced form for main effect terms in at; calculates and displays sum-reduced and row echelon form 
at:  available terms; string vector of the form {"a","b","c",...}
df: degrees of freedom for terms in at; {dfa,dfb,dfc,....}
xi: submatrix of Inverse[X'X] formed by dropping 1st row and column (constant term)
inxi: start position and number of vectors in xi for each term in at; {{"a",s,n"},{"b",s,n},....}
ems: Mean Square error for F test
edf:  degrees of freedom for ems
mefb: parameter estimates for main effects; {{a1,a2,..},{b1,b2...},....}
Rqd: LinearAlgebra`MatrixManipulation`
Calls: binselect, probF,hytst 
Mod. Date: Dec 8/97

\(\(getucn[at_, df_, xi_, inxi_, ems_, edf_, mefb_] := 
        Module[{iflg = 1, bs, cntrm, dfk, crw, kcon, kct, rwrc, erf, sbinx, 
            si, di, subxi, prr, sumsk}, 
          Print["
User Defined Contrasts"]; 
          While[iflg == 1, 
            cntrm = kbinstrg["Type contrast term", Transpose[{at}]]; 
            Print["User selected contrast term=", cntrm]; 
            bs = \((#1 == cntrm &)\) /@ Flatten[at] /. {True -> 1, 
                  False -> 0}; 
            dfk = \(Flatten[
                  binselect[
                    bs, {df}]]\)[[1]]; 
            erf = 1; 
            While[erf == 1, crw = kbin3["How many rows in contrast?", 0]; 
              If[crw > dfk, Print["Error: #rows must be <=", dfk], 
                erf = 0];]; Print["Rows in contrast=", crw]; 
            Print["Required Cols in contrast=", dfk + 1]; erf = 1; 
            While[erf == 1, 
              kcon = kbin3["Type contrast matrix", {crw, dfk + 1}]; 
              Print["User Contrast matrix=", kcon]; 
              rwrc = RowReduce[kcon]; 
              Print["Row Reduced contrast is=", rwrc]; 
              If[False === \((\((Table[
                                0, {\(Dimensions[
                                      rwrc]\)[[2]]}] == #1 &)\) /@ rwrc == 
                      Table[False, {\(Dimensions[
                              rwrc]\)[[1]]}])\), Print["Error:rank<row order"], erf = 0];]; 
            hytst[cntrm, kcon]; kct = Transpose[kcon]; 
            kcon = Transpose[
                Drop[kct + Table[\(-Last[kct]\), {Length[kct]}], \(-1\)]]; 
            Print["Sum-reduced contrast matrix is=", kcon]; 
            Print["Row echelon form of contrast=", RowReduce[kcon]]; 
            bs = \((cntrm === \(inxi[[#1]]\)[[1]] &)\) /@ 
                  Range[Length[inxi]] /. {True -> 1, False -> 0}; 
            sbinx = binselect[bs, Transpose[inxi]]; 
            si = Flatten[
                Table[sbinx[[2]], {2}]]; 
            di = Flatten[
                Table[sbinx[[3]], {2}]]; 
            subxi = SubMatrix[xi, si, di]; 
            prr = mefb[[
                Flatten[
                  binselect[bs, {Range[Length[mefb]]}]]]]; 
            prr = Flatten[prr]; 
            Print["Parameter estimates used=", N[prr, 6]]; 
            sumsk = ssk[kcon, subxi, prr]; 
            probF[sumsk, Length[kcon], ems, edf]; 
            If["n" == 
                kbinstrg["More contrasts? y/n", {{"y"}, {"n"}}], 
              iflg = 0];];];\)\)

hytst[t,k] displays string expression for hypothesis being tested for term "t" and contrast matrix k
t: term in string form, eg., "a", "ab"
k: contrast matrix
Note:Chop[N[k]] is used so that "=0" is not displayed at divisor position when rational values
         are entered in k and Chop to ensure a term multiplied by 0. is removed.
Mod. Date: July 17/97

hytst[t_, k_] := 
        Module[{prd, stcf}, 
          Print[StringInsert["Hypothesis tested for term: ", t, -1]]; 
          prd = Chop[N[k]] . 
              Array[t, Dimensions[
                    k][[2]]]; 
          stcf = MapThread[ToString[#1] &, {prd}]; 
          MapThread[Print[#1, "=0"] &, {stcf}];];

lctmev[d,at] calc starting and ending position of vectors in X'X 
d: degrees of freedom for terms in at (number of row/col for each term)
at: terms ordered,e.g., {"a","b","c"...} corresponding to terms represented in X'X
Returns: {{"a",begin row/col,#of row/col},{"b",begin row/col,#of row/col},....}
Mod. Date: July 8/97

lctmev[d_, at_] := 
        Module[{lse, vn, bgn, indxi}, lse = {}; vn = Range[Plus @@ d]; 
          Do[AppendTo[lse, 
              Take[vn, d[[j]]]]; vn = 
                RotateLeft[vn, 
                  d[[j]]];, {j, 1, 
              Length[d]}]; bgn = (Take[#1, 1] &) /@ lse; 
          indxi = Transpose[Append[{bgn}, Transpose[{d}]]]; 
          indxi = (Flatten[#1, 1] &) /@ indxi; 
          indxi = (Prepend[
                    indxi[[#1]], 
                    at[[#1]]] &) /@ 
              Range[Length[at]]; Return[indxi];];

gencon[s,d,p,xtxi,mse,dfe]: calc SS,MS, F, and prob for each main effect and interaction term and displays result.Prob is rounded to 4 decimal digits.
s: contrasts available in the form of {{"mu",{0,0,0},1},{"a",{0,0,1},2},....}
d: vector of degrees of freedom for each term in s
p: vector of parameter estimates
xtxi: Inverse[X'X]
mse: Mean Square error
dfe:  degrees of freedom for mse
Calls: ssk, rndoff,restrv,lctmev
Required: LinearAlgebra`MatrixManipulation
Mod. Date: Dec 5/97

\(\(gencon[s_, d_, p_, xtxi_, mse_, dfe_] := 
        Module[{dftrm, strm, con, inxt, rxi, par, lcn, si, di, ct, msk, fr, 
            fd, pr, sumsq}, strm = Drop[Flatten[Take[Transpose[s], 1]], 1]; 
          dftrm = Drop[d, 1]; par = restrv[Drop[p, 1], dftrm]; 
          lcn = lctmev[dftrm, strm]; di = Length[xtxi] - 1; 
          rxi = SubMatrix[xtxi, {2, 2}, {di, di}]; ct = {}; 
          Do[con = 
              IdentityMatrix[
                dftrm[[j]]]; 
            si = \(lcn[[
                  j]]\)[[2]]; 
            di = \(lcn[[
                  j]]\)[[3]]; inxt = SubMatrix[rxi, {si, si}, {di, di}]; 
            sumsq = ssk[con, inxt, 
                par[[j]]]; 
            msk = sumsq\/dftrm[[j]]; 
            If[dfe ≥ 1, fr = msk\/mse; 
              fd = FRatioDistribution[
                  dftrm[[j]], dfe]; 
              pr = rndoff[1.  - CDF[fd, fr], 4], fr = "*"; 
              pr = "*"]; \(AppendTo[
                ct, {strm[[j]], 
                  N[sumsq, 6], 
                  dftrm[[j]], N[msk, 6], 
                  fr, pr}];\), {j, 1, Length[strm]}]; 
          AppendTo[ct, {"Error", mse\ dfe, dfe, mse}]; 
          Print["
Anova Table for all Terms (Type III SS)"]; 
          Print[TableForm[ct, 
              TableHeadings -> {Automatic, {"Source", "SS", "df", "MS", "F", "P"}}, 
              TableSpacing -> {0, 2}]];];\)\)

anovanw[d] does anova using data set d 
 Calls Directly: kbin3, srtgrp, srtcnt, kbinstrg, dsgnmef, dmpar, dsgnint, cutintr, anvT2, gencon, intrplot, lctmev, restrv, getucn       
Special Variables:
   xdsm: design matrix; starts with main effects; interaction appended by dsgnint
   xdscode: holds defn of all terms in the form {"a",{0,0,1}},{"b",{0,1,0}},...{"abc,{1,1,1}}..
   slctd: holds defn of currently available terms in form of xdscode as well as location of the term  in xdsm. 
 Mod. Date: Dec 6/97
 Mod. Date: Jun 4/02: changed local variable list

\(\(anovanw[d_] := 
        Module[{ntot, nf, ntrm, nlev, xdsm, xdscode, slctd, sstot, sse, 
            means, conv, ncs, xtx, xxi, df, gm, bpar, ssmod, temp, mse, dfe, 
            ixr, xtx2, bpar1, cutf = 0, mef, dfa, xindx, bmef}, 
          Print["GLM Anova V3.0"]; 
		If[0 == Length[d], 
            Message[anovanw::nodat]; 
            Abort[]]; 
		Print["Number of indices defining factor levels=", 
            Length[\(d[[1]]\)[[1]]]]; 
          nf = kbin3["Type the number of factors in the design", 0]; 
          Print["Number of factors=", nf]; 
          If[True == 
              VectorQ[\(d[[1]]\)[[2]]], 
            Print["Data input is in grouped form"]; {sstot, sse, nlev, 
                ncs, means, gm} = srtgp[d, nf], 
            Print["Data input is in ungrouped form"]; {sstot, sse, nlev, 
                ncs, means, gm} = srtcnt[d, nf];]; ntot = Plus @@ ncs; 
          Print["In this design there will be:"]; ntrm = \(-1\) + 2\^nf; 
          Print["Number of terms(mu,a,b,ab,..)=", ntrm + 1]; 
          Print["Main effect parameters(mu,a,b,c..) estimated=", 
            1 + Plus @@ \((nlev - 1)\)]; 
          Print["Interaction parameters(ab,ac,...) estimated=", 
            Times @@ nlev - \((1 + Plus @@ \((nlev - 1)\))\)]; 
          If[1 > Plus @@ ncs - Times @@ nlev, 
            Print["
WARNING: No df available for SSerror."]; 
            Print["Some interaction terms must be removed to do F tests."]]; 
          If["n" == 
              kbinstrg["Do you wish to proceed? y/n", {{"y"}, {"n"}}], Return[]]; xdsm = dsgnmef[nlev]; xdscode = dmpar[nf]; 
          slctd = xdscode[[2\^\(\(-1\) + Range[nf]\)]]; 
          slctd = \((Append[
                    slctd[[#1]], #1] &)\) /@ Range[Length[slctd]]; 
 (*Print["xdscode=", 
              xdscode]; 
Print["slctd=", 
              slctd]; 
Print["xdsm=", 
              xdsm];*) 

{slctd, xdsm} = 
            dsgnint[xdscode, slctd, 
              xdsm]; 
 (*Print["slctd=", 
              slctd]; 
Print["xdsm=", 
              xdsm]; 
Abort[];*) 
If[
            nf ≥ 2 && "y" == 
                kbinstrg["Remove interaction terms? y/n", {{"y"}, {"n"}}], cutf = 1; {xdsm, slctd} = cutintr[xdsm, slctd, nf]]; 
          ntrm = Length[slctd]; 
          df = \((Take[
                    Dimensions[
                      xdsm[[#1]]], 
                    1] &)\) /@ Range[ntrm]; xdsm = Flatten[xdsm, 1]; 
          conv = Table[1, {Times @@ nlev}]; xdsm = Join[{conv}, xdsm]; 
          df = Prepend[Flatten[df], 1]; 
          slctd = Prepend[slctd, {"mu", Table[0, {nf}], 0}]; 
          slctd = \((ReplacePart[
                    slctd[[#1]], #1, \(-1\)] &)\) /@ Range[Length[slctd]]; 
          xtx = \((ncs\ xdsm[[#1]] &)\) /@ 
              Range[Length[xdsm]]; 
		xtx = 
            xtx . Transpose[xdsm]; 
	\ \ xxi = 
            Inverse[xtx];  (*xxi\ is\ in\ rational\ form*) 
		bpar = 
            xxi . xdsm . \((means*\ ncs)\); 
		temp = 
            N[restrv[bpar, df], 16]; 
          Do[Print["
Term=", \(slctd[[
                  j]]\)[[1]], " df=", 
              df[[j]]]; 
            Print["Par. Est=", 
              temp[[j]]], {j, 1, 
              Length[df]}];  (*If\ any\ interaction\ term\ was\ dropped, \ 
            calc\ new\ sse*) 
		If[
            cutf == 1, 
		bpar1 = Drop[bpar, 1]; 
		xtx2 = 
              SubMatrix[
                xtx, {2, 2}, {Length[bpar1], Length[bpar1]}]; 
		ssmod = 
              bpar1 . xtx2 . bpar1; 
		Print["New SSmcfm=", 
              ssmod, " (=", N[ssmod, 16], ")"]; 
		sse = 
              sstot - ssmod; 
		Print["New SSe=", sse, " (=", 
              N[sse, 16], ")"];
			];  (*end\ If\ cutf == 
              1*) 
		{mse, dfe} = anvT2[sstot, sse, {ntot, Length[bpar]}]; 
          gencon[slctd, df, bpar, xxi, mse, dfe]; If[dfe == 0, Return[]]; 
          If[nf ≥ 
              2, \(If["y" == 
                  
                  kbinstrg["Plot interactions? y/n", {{"y"}, {"n"}}], intrplot[means, slctd, nlev]];\)]; 
          If["y" == 
              kbinstrg["Contrasts on main effects? y/n", {{"y"}, {"n"}}], 
            mef = Take[
                Drop[\(Transpose[
                      slctd]\)[[1]], 1], 
                nf]; dfa = Take[Rest[df], nf]; xindx = lctmev[dfa, mef]; 
            ixr = SubMatrix[xxi, {2, 2}, Table[\(-1\) + Length[xxi], {2}]]; 
            temp = Take[Rest[bpar], Plus @@ dfa]; bmef = restrv[temp, dfa]; 
            getucn[mef, dfa, ixr, xindx, mse, dfe, bmef];];];\)\)

anovanw::nodat = "Data presented has Length=0;likely not evaluated.\
";

Textbook Examples

anovanw[dat1]

anovanw[dat2]

anovanw[dat3]

anovanw[dat4]

anovanw[dat4]

anovanw[dat5]

NIST Examples

See the section Generating the NIST Simon-Lesaage Data Sets for functions used to generate the data sets used below.

anovanw[simon7]

Clear[simon7]

anovanw[simon7u]

Clear[simon7u]

anovanw[simon8]

Clear[simon8]

anovanw[simon9]

Clear[simon9]

Wolfram V4.2 Examples

The data was obtain from the site www.wolfram.com/Mathematica/newin42. The data is used in an example of new statistical functions available in Mathematica Version  4.2.

w1 = {{{1}, {7.0, 5.3, 5.9, 6.6, 4.9}}, {{2}, {4.4, 6.8, 7.7, 8.3, 
              6.6}}, {{3}, {8.1, 10.4, 8.0, 6.8, 9.2}}, {{4}, {5.7, 3.9, 6.2, 
              5.5, 6.2}}};

anovanw[w1]

w2 = {{{1}, {125, 130, 135}}, {{2}, {142, 100, 132}}, {{3}, {162, 
              181, 173}}, {{4}, {188, 192, 183}}, {{5}, {221, 225, 
              231}}, {{6}, {253, 255, 237}}, {{7}, {272, 280, 
              276}}, {{8}, {345, 340, 356}}, {{9}, {400, 390, 
              397}}};

anovanw[w2]

Kirk p 460 Nested Design

The example given by Kirk is a 2-way design having the following characteristics:

Factor  A: fixed with 2 levels
Factor B: B is random and nested  within A; 4 levels
Number  of observations per cell: 4

The data  is shown below.

Kirk Data

krk = {{3, 1, 5, 2, , , , }, {6, 2, 6, 3, , , , }, {3, 2, 5, 
          4, , , , }, {3, 2, 6, 
          3, , , , }, {,,,,7,4,7,10}, {,,,,8,5,8,10}, {,,,,7,4,9,9}, {,,,,6,3,8,11}}

krk2 = krk /. {Null -> ""}

\(Print[
      TableForm[krk2, TableSpacing -> 1, 
        TableHeadings -> {{"A1", "A1", "A1", "A1", "A2", "A2", "A2", "A2"}, {"B1", "B2", "B3", "B4", "B5", "B6", "B7", "B8"}}]]\)

In order to analyze the data using a factorial design the data is re-cast into a 3-way design in which the 3rd factor contains the data points  with as many levels as there are  data points. Thus,  the 3-way analysis contains one observation per  cell. An important restriction is that the number of data points  must be equal in each cell of the original nested design. The data re-cast as a 3-way anova with one observation per cell and with cell subscripts as required for analysis using functions in this notebook are given below.

krd = {{{1, 1, 1}, 3}, {{1, 1, 2}, 6}, {{1, 1, 3}, 3}, {{1, 1, 4}, 
            3}, {{1, 2, 1}, 1}, {{1, 2, 2}, 2}, {{1, 2, 3}, 2}, {{1, 2, 4}, 
            2}, {{1, 3, 1}, 5}, {{1, 3, 2}, 6}, {{1, 3, 3}, 5}, {{1, 3, 4}, 
            6}, {{1, 4, 1}, 2}, {{1, 4, 2}, 3}, {{1, 4, 3}, 4}, {{1, 4, 4}, 
            3}, {{2, 1, 1}, 7}, {{2, 1, 2}, 8}, {{2, 1, 3}, 7}, {{2, 1, 4}, 
            6}, {{2, 2, 1}, 4}, {{2, 2, 2}, 5}, {{2, 2, 3}, 4}, {{2, 2, 4}, 
            3}, {{2, 3, 1}, 7}, {{2, 3, 2}, 8}, {{2, 3, 3}, 9}, {{2, 3, 4}, 
            8}, {{2, 4, 1}, 10}, {{2, 4, 2}, 10}, {{2, 4, 3}, 9}, {{2, 4, 4}, 
            11}};

The results of the anova are given below. Because there is only one observation per cell, the SSerror and dfe is 0. The warning in this regard is ignored and no F  tests are  carried out.

anovanw[krd]

The characteristics of the analysis  required of the nested design are  as follows:
Factor  A: fixed with 2 levels
Factor B: B is nested  within A (b:a); random with 4 levels
Factor C: C is random and nested  within A and B (c:ab)
Number  of observations per cell: 1

Using the above characteristics, the function ems3 in the notebook MGems.nb uses the Millman-Glass heuristic to determine the various sources of variation and their expected mean squares. The output from MGems.nb is given below.

ems3[4]

The next step is to combine various SS and df from the 3-way analysis to form the appropriate SS and df for the nested design. For the nested design the sources of variation are given in the last table as A, B:A, and C:AB. A heuristic for forming the combinations is as follows:

1) If a term is not nested, the associated SS and df are used without change.
2) If a term is nested, then specific SS and associated df must be summed.  All the sources of the SS for a given term and the df must be summed.  In the example above B is nested within A, i.e., b:a, so the SS for b:a is formed by summing the SS(b) and all combinations of b with those terms appearing to the right of the colon , i.e., SS(ab). If more than one term appears to the right of the colon, then SS(b) is summed with each SS term taken one at a time, two at a time, etc.  (This is illustrated in the next step.) The SS(b:a) now becomes 75.25+29.25 = 104.50 and df(b:a)=3+3=6.
3) The SS for the term c:ab now includes SS(c)+SS(ac)+SS(bc)+SS(abc), i.e., SS(c:ab)=5.25+0.75+6.00+6.5=18.50. The df(c:ab)=3+3+9+9=24.
4)  The mean squares are formed in the usual manner using the new SS and df. For example, the mean square for the term A is 112.50/1=112.50; the mean square for b:a is 104.50/6=17.417; the mean square for c:ab is 18.50/24=0.771.

Determining which mean square to use in the denominator to obtain F tests for A and B nested within A requires knowing the expected mean squares for each term being tested. 

From column labelled EMS of the last table of the above output, factor A is tested using the mean square of B:A in the denominator, and B:A is tested using the mean square for C:AB in the denominator. The F test for A is given by 112.50/17.417=6.46 with df=(1,6). The F test  for B:A is  17.417/0.771=22.59 with df=(6,24).

Winer Repeated Measure p 526.

In this example there are 2 factors with repeated measurements on the second factor. The characteristics are as follows:

1) Factor A: 2 levels; fixed
2) Factor B: repeated measures factor with 4 levels (4 repetitions); fixed
3) There are a total of 6 subjects, 3 in A1 and 3 in A2

w526 = {{0, 0, 5, 3}, {3, 1, 5, 4}, {4, 3, 6, 2}, {4, 2, 7, 8}, {5, 4, 
          6, 6}, {7, 5, 8, 9}}

\(Print[
      TableForm[w526, TableSpacing -> 1, 
        TableHeadings -> {{"A1, s1", "A1, s2", "A1, s3", "A2, s4", "A2, s5", "A2, s6"}, {"B1", "B2", "B3", "B4"}}]]\)

As in the example taken from Kirk, a new factor C is defined containing 3 levels and holding the data points of each subject. Factor C is nested under Factor A. This results in a 3-way design having 1 observation per cell. The data structured for the anova analysis is shown below.

rmd={{{1, 1, 1}, 0}, {{1, 1, 2}, 3}, {{1, 1, 3}, 
          4}, {{1, 2, 1}, 0}, {{1, 2, 2}, 1}, {{1, 2, 3}, 3}, {{1, 3, 1}, 
          5}, {{1, 3, 2}, 5}, {{1, 3, 3}, 6}, {{1, 4, 1}, 3}, {{1, 4, 2}, 
          4}, {{1, 4, 3}, 2}, {{2, 1, 1}, 4}, {{2, 1, 2}, 5}, {{2, 1, 3}, 
          7}, {{2, 2, 1}, 2}, {{2, 2, 2}, 4}, {{2, 2, 3}, 5}, {{2, 3, 1}, 
          7}, {{2, 3, 2}, 6}, {{2, 3, 3}, 8}, {{2, 4, 1}, 8}, {{2, 4, 2}, 
          6}, {{2, 4, 3}, 9}}

anovanw[rmd]

Next, the Millman-Glass procedure is used to determine all the sources of variation which are available in the repeated measure design and to determine which SS from the crossed design must be summed.

ems3[4]

From the above output the sources of variation for the repeated measure design are A, B, AB, C:A (C nested within A) and BC:A (BC nested within A). The SS for A, B, and AB can be used directly as these sources are not nested. To obtain the SS for

1) C:A we form SS(c)+SS(ac): 14.5833+2.58333=17.16663 with 2+2=4 df.
2) BC:A we form SS(BC)+SS(ABC): 9.41667+5.41667=14.8334 with 6+6=12 df

Also from the output given above (column labelled EMS)  the denominator mean-squares for F tests can be determined:

1) To test A, the mean-square for C:A is used in the denominator
2) To test B, the mean-square for BC:A is used in the denominator
3) To test AB, the mean square for BC:A is used

The F tests are:
F(A)=(51.0417/1) / (17.1663/4)=11.8932, df=(1,4)
F(B)= (47.4583/3) / (14.8334/12) = 12.7977, df=(3,12)
F(AB)=(7.45833/3) / (14.8334/12)= 2.01123, df=(3,12)