Analysis of Variance GLM
Steve Hunka, University of Alberta
steve.hunka@ualberta.ca
Mod. Date: August 5/97

General Description

The function anovanw[d] calculates an analysis of variance (fixed effects) 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,SSmean, SSmodel, etc.) 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). The functions were written  and tested using Version 2.2.2 for the Macintosh.

References
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.

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.  
2.  The input data are sorted by cell index and the number of observations, means, and variances are calculated for each cell and SS total. 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 identifical, * 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 general SS terms SStotal, SSmean, SStcfm (total corrected for the mean), SSmod (model), SSmcfm (model corrected for the mean) and SSerror are displayed. A F test is done on SSmcfm. If df(error)=0, no F tests are calculated, but SS are calculated.
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 dfe(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 3 decimal digits. Thus, a probability of 0 would be reported as <.001 .

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 would be acceptable, e.g., {1,2,4}, {2,1,3},and so on.
3. 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; (b) the user may have underestimated the number of parameters to be estimated. 
4.  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, 

5. 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 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 4 above). The term entered must have characters in their natural order, e.g., abc, and not acb. 
6.  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 since 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.

Example Data Sets
Data set  dat1 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. Data set dat2 contains the same dependent variable as in data set dat1 and has only 3 cell indices (the same as the first 3 in dat1), but all observations belonging to the same cell are grouped together.  Data set dat1 and dat2 are an equal-N case. Data set dat3 is also from Winer (p446) and is a 2-factor unequal-N case. Data set dat4 has 1 observation per cell and is from Shavelson and Webb (p43).

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={{{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={{{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 comes from Shavelson,R.J. & Webb, N,M., Generalizability Theory A Primer, Sage
Publications Inc., 1991, p43. There is 1 observation per cell, and df(error)=0 unless one or 
more interaction terms are removed from the model.

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}};

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. "Remove interaction terms? y/n": n
5. "Plot interactions? y/n": y
   5.1 "Type an interaction term to plot": ab
   5.2 "More plots? y/n": y
   5.3 "Type an interaction term to plot": abc
   5.4 "More plots? y/n": n
6. "Contrasts on main effects? y/n": y
   6.1 "Type contrast term": b
   6.2 "How many rows in contrast?: 2
   6.3 "Type contrast matrix": {{1,0,-1},{0,1,-1}}
   6.4 "More contrasts": n

Required Packages

<<Statistics`DataManipulation` 
<<LinearAlgebra`MatrixManipulation`
<<Statistics`ContinuousDistributions`

List of Functions Used

anovanw: main anova program
anvT1:  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
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
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
getMp:   gets selected set of rows of data 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
meanVar: calc. means and variances
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
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.

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,(*check for character*)
       Print[ans," Input Error; Characters not allowed"];,
       (*integer or real required*)
      Head[p]===Integer || Head[p]===Real, 
      (*chk for vec when scalar wanted*)
      If[Head[p]=!=Head[ans],
       Print[ans," Input Error. Input must be ",Head[p]];,
       c=1];,
       
      Head[p]===List && Length[p]==1, (*vector rqd*)
       If[VectorQ[ans]==False, 
       Print[ans,StringJoin[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, (*matrix rqd*)
      If[MatrixQ[ans]==False,
      Print[ans,StringJoin[vrp2,ToString[p]]];,
          If[Dimensions[ans]!=p,
      Print[ans," Error; Matrix dimensions must be ",p];,
       c=1]]; 
     ]; (*end Which *)    
   ]; (*end While *)
 Return[ans];
 ]; (* end Module*)

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]}]; (* end Do *)  
   If[c!=1,Print[ans," Input Error. Type Again"]];  
   ]; (*end While*)
 Return[ans];
 ]; (* end Module*)

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=Part[Dimensions[m1],1];
 a=Floor[c/pc];
 indx=Join[Table[pc,{a}],{Mod[c,pc]}];
 If[0==Part[indx,-1],indx=Drop[indx,-1]];
 
 If[Length[{rch}]==0,
   row=Map[ToString[#]&,Range[Part[Dimensions[m],1]]];
   col1=Map[ToString[#]&,Range[Part[Dimensions[m],2]]];
   ,row=rch[[1]];
    col1=rch[[2]];
  ]; (*end If*)
 
 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,{Part[rc,1]}];
ri=i1-1; ci=i2-1;
a=Map[fn[a[[#]],ri+#]&,Range[Part[rc,1]]];
b=Table[ch,{Part[rc,2]}];
b=Map[fn[b[[#]],ci+#]&,Range[Part[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]];

srtcnt[d,nf] sorts anova input data (ungrouped) and calculates number of observations in each cell, means, and variances.
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: prints each cell index and the number of observations in the cell; means and variance
Returns: {sst,ord,freq,res2}
sst: total sum of squares
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, 
res2: the cell means ordered with last subscript varying the fastest
Note: A cell which has zero observations is not directly identified. 
Calls: kbin3,pmat2,genrch,getMp,meanVar
Required: Statistics`DataManipulation`
Mod. Date: August 7/97

srtcnt[d_,nf_]:=Module[{nob,dn,si,obgrp,rch,mv,frq,res2,
         res3,ord,sst},
nob=Length[d]; Print["Number of Observations=",nob];
 si=kbin3["Type vector of Factor indices to use",{nf}];
 Print["Selected Factor Indices=",si];
 (*select/reorder indices; sort and count*)
 dn=Sort[Map[ReplacePart[Part[d,#],
   Part[Part[d,#],1][[si]],1]&,Range[nob]]];
 obgrp=Frequencies[Map[Part[Part[dn,#],1]&,Range[nob]]];
 obgrp=Map[Reverse[obgrp[[#]]]&,Range[Length[obgrp]]];
 obgrp=Partition[Flatten[obgrp],nf+1];
 ord=Map[Max[Transpose[obgrp][[#]]]&,Range[nf+1]];
 ord=Take[ord,nf];
Print["Number of Levels in each factor=",ord]; 
dn=Flatten[Rest[Transpose[dn]]]; (*dep. variable *)
frq=Flatten[Take[Transpose[obgrp],-1]]; (*cell n(j) *)
mv=Map[meanVar[getMp[#,dn,frq]]&,Range[Length[frq]]];
mv=N[mv];
res2=Map[Flatten[mv[[#]]]&,Range[Length[frq]]];
res3=Map[Join[obgrp[[#]],res2[[#]]]&,Range[Length[frq]]];
rch=genrch["Cell ","F",Dimensions[res3]];
(*replace last 3 column names *)
rch[[2]]=Join[Drop[rch[[2]],-3],{"Count","Mean","Var"}];
pmat2[res3,6,rch];
res2=Flatten[Take[Transpose[res2],1]]; (*means only*)
sst=N[Apply[Plus,Part[Transpose[d],2]^2]];
Print["SS(total)=",sst];
 Return[{sst,ord,frq,res2}];  
 ]; (*end Module *)

srtgp[dg,nf] takes data grouped by cell index, and for nf factors and a selection of nf cell
indices, combines the data of redundant cells then calculates the cell means and variances, 
the new cell frequencies, SStotal, and displays the result.
dg: data matrix in the form of {{{cell indices},{Y for each observation in the cell}},.....}
nf: number of factors in the design
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,orf,obgrp,mv}, where sst is SStotal, orf is the number of levels in each factor
 of the anova design, obgrp are the cell frequencies, and mv the cell means.
Calls: kbin3,comdta, meanVar
Mod. Date: August 7/97

srtgp[dg_,nf_]:=Module[{si,dgri,sd,orf,obgrp,mv,ind,
            stuf,rch,sst},
si=kbin3["Type vector of Factor indices to use",{nf}];
Print["Selected Factor Indices=",si];
(*select cell indices according to user input si*)
dgri=Map[ReplacePart[Part[dg,#],Part[Part[dg,#],1][[si]],1]&,
        Range[Length[dg]]];
(*sort the data so redundant cell indices are together*)
sd=Sort[dgri];(*next calc new design designation*) 
orf=Map[Max[Transpose[Map[Part[#,1]&,sd]][[#]]]&,Range[nf]];
Print["Number of Levels in each factor=",orf];
{ind,sd}=comdta[sd];(*rearrange the data*)
obgrp=Map[Length[#]&,sd];(*frequencies*)
mv=N[Map[meanVar[#]&,sd]];(*means and variances*)
(*arrange data for display*)
stuf=Map[Join[ind[[#]],{obgrp[[#]]}]&,Range[Length[ind]]];
stuf=Map[Join[stuf[[#]],Flatten[mv[[#]]]]&,
        Range[Length[ind]]];        
rch=genrch["Cell ","F",Dimensions[stuf]];
rch[[2]]=Join[Drop[rch[[2]],-3],{"Count","Mean","Var"}];
pmat2[stuf,6,rch];
sst=Apply[Plus,Map[Apply[Plus,#^2]&,sd]];
Print["SStotal=",sst];
mv=Flatten[Part[Transpose[mv],1]];(*get means only*)
Return[{sst,orf,obgrp,mv}];
];

comdta[d] combines data in common data cells of d, finds unique cell indices, and 
frequency
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: August 7/97

comdta[d_]:=Module[{di,dd,newc,bs,dn,sd},
di=Flatten[Take[Transpose[d],1],1]; (*index part*)
dd=Part[Transpose[d],2]; (*data part*)
newc=List[]; sd=List[];

While[Length[di]>0,
 bs=Map[di[[#]]==di[[1]]&,
 Range[Length[di]]] /. {True->1,False->0};
 sd=Append[sd,Transpose[binselect[bs,Transpose[dd]]]];
 dn=Apply[Plus,bs]; (*number of elements to drop*)
 AppendTo[newc,Take[di,1]];(*unique index list*)
 di=Drop[di,dn];(*drop from index list*)
 dd=Drop[dd,dn];(*drop from data list*)
 ]; (*end While *)
 newc=Flatten[newc,1];(*new index list*)
 sd=Map[Flatten[#]&,sd];(*new data list*)
 Return[{newc,sd}];
];

meanVar[d] calculates mean and variances of data matrix d or vector d. The mean (m) and variance (v) are returned as {m,v}.
Mod. Date: Oct 15/96

meanVar[d_]:=
 Module[{m,v,d1},
 If[1==Length[Dimensions[d]],
    d1=Partition[d,1],d1=d];
  m=Apply[Plus,d1]/Length[d1];
  v=(Apply[Plus,(d1-Table[m,{Length[d1]}])^2])/Length[d1];
  Return[{m,v}];
 ];

getMp[p_,x_,n_]  gets the pth part of matrix x; x has {n1,n2,...np} rows
          in each part. Operates on vectors also.
Mod. Date:

getMp[p_,x_,n_]:=   
 Module[{inx,mp},
  inx=Drop[FoldList[Plus,0,n],-1];
  mp=Take[Drop[x,inx[[p]]],n[[p]]];
  Return[mp];
       ];

anvT1[st,sm,smd,df] calculates SS and displays a general anova table. If df(error) is 0, 
MS(error) is set to 0, and no F tests are carried out.
st: SS total
sm: SS mean
smd: SS model
df: degrees of freedom for SS total and SS model
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
Mod. Date: July 23/97

anvT1[st_,sm_,smd_,df_]:=Module[{ct,dft,dfmd,dfe,mser,ssmcfm,fr,
      fd,p},
Print["\nGeneral Anova Table"];
ct=List[];dft=df[[1]];dfmd=df[[2]];dfe=dft-dfmd;
AppendTo[ct,{"SStot",st,dft,st/dft}];
AppendTo[ct,{"SSmean",sm,1,sm}];
AppendTo[ct,{"SStcfm",st-sm,dft-1.,(st-sm)/(dft-1.)}];
AppendTo[ct,{"SSmod",smd,dfmd,smd/dfmd}];
ssmcfm=smd-sm;
If[dfe>0,
mser=N[(st-smd)/dfe]; (*MS error*)
fr=(ssmcfm/(dfmd-1))/mser;
fd=FRatioDistribution[dfmd-1,dft-dfmd];
p=rndoff[1.-CDF[fd,fr],3],
fr="*";p="*";mser=0];
AppendTo[ct,{"SSmcfm",ssmcfm,dfmd-1,ssmcfm/(dfmd-1),fr,p}];
If[dfe>0,
AppendTo[ct,{"SSerror",st-smd,dft-dfmd,mser}],(*else*)
AppendTo[ct,{"SSerror",Chop[st-smd],0,"*"}]];
Print[TableForm[ct,
 TableHeadings->{Automatic,{"Source","SS","df","MS","F","P"}}, 
   TableSpacing->{0,2}]];
Return[{mser,dfe}];
]; (*end Module*)

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],3];
 ct={{"Usr K",ss,dfs,ms,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: July 23/97

dsgnmef[ln_]:=Module[{kln,dm,e1,e2,,e3,r1,r2,r3,r4},
kln=Join[{1},ln];dm=List[];
Do[(*over number of terms*)
e1=Table[Join[Table[0,{ln[[j]]}],{1}],
     {Apply[Times,ln]}];(*gen {0,0...1} sequence*)
(*generate left shift indices*)
e2=Transpose[Table[Reverse[Range[ln[[j]]]],
     {Apply[Times,Drop[ln,j]]}]];
(*replicate left shift indices*)
e3=Flatten[Table[Flatten[e2],
            {Apply[Times,Take[kln,j]]}]];
(*now do the left shifts & drop the last col of 0s *)
r1=Map[RotateLeft[e1[[#]],#]&,e3];
r1=Drop[Transpose[r1],-1];
r2=Take[r1,ln[[j]]-1];
r3=Flatten[Take[r1,-1]];(*do the sum reduction*)
r4=MapThread[Flatten[#]-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]]];
Map[If[1==x[[#]],
   res=StringInsert[res,alph[[#]],-1];]& (*end If *)
 ,Reverse[Range[Length[x]]]]; (* end Map *)
 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}},......}
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[[#]]]}],bin[[#]]]&,
     {Range[-1+2^nf]}];
bin3=MapThread[bin2alf[bin2[[#]]]&,{Range[-1+2^nf]}];
bin4=Transpose[Append[{bin3},bin2]];
Return[bin4];
]; (*end Module *)

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
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.
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_]:=Map[a b[[#]]&,Range[Length[b]]];
avc2=avc; cdm2=cdm;
(*1st find the interaction terms of apc*)
Do[  (* Do over all component terms of apc*)
  If[1!=Apply[Plus,apc[[j]][[2]]], (*interaction term found*)
    (*determine vectors to multiply*)
    bin1=apc[[j]][[2]]; (*required term*)
    Do[(*Do backward search over available component terms*)
      dif=bin1-avc2[[k]][[2]];
      If[
      1==Apply[Plus,dif^2], (*only one 1 exists *)
      lcn1=avc2[[k]];
      (*locate posn of dif in avc*)
      lcn2=findc[dif,avc2];
     (*now get the vectors to multiply*)
      mat1=cdm2[[lcn1[[3]]]];
      mat2=cdm2[[lcn2[[3]]]];
     (*now form the product of the vectors*)
      prd=Map[fvp[mat1[[#]],mat2]&,Range[Length[mat1]]];
      prd=Flatten[prd,1];
      AppendTo[cdm2,prd];(*update design matrix*)
      newterm=bin2alf[bin1];
      AppendTo[avc2,{newterm,bin1,Length[cdm2]}];(*update avc2 *)
      Break[];
      ]; (*end If *)
      ,{k,Length[avc2],1,-1}]; (*end Do backward search*)
    ]; (*end If sum!=1 *)
  ,{j,1,Length[apc]}]; (*end Do over all component terms*)
  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[ (*over number of available components*)
 If[c==ac[[j]][[2]],
 res=ac[[j]];
 Break[];
 ]; (*end If*)
,{j,1,Length[ac]}]; (*end Do*)
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},
(*define acceptable answers as interaction terms only*)
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+Map[strm==it[[#]][[1]]&,Range[Length[it]]] /.
  {False->0,True->1};
  If["n"==kbinstrg["Any more terms? y/n",{{"y"},{"n"}}],
     flg=1];
 ]; (*end While*)
 (*in case user drops same term more than once*)
 bs=Map[bs[[#]]===0&,Range[Length[bs]]] /. {False->1,True->0};
 (*bs contains 1 in lcn to cut else 0*)
 bs=Table[1,{Length[it]}]-bs;
indx=binselect[bs,{Range[Length[it]]}];
indx=Flatten[indx];
xc=x[[indx]];
itc=it[[indx]];(*next renumber position in itc*)
itc=Map[ReplacePart[itc[[#]],#,-1]&,Range[Length[itc]]];
Return[{xc,itc}];
]; (*end Module*)

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=List[];
Do[cr=Apply[Plus,Part[sl[[j]],2]];
If[2==cr || 3==cr,
   AppendTo[sl2,sl[[j]]];
];(*sl2 holds only 2 and 3 term interaction components*)
,{j,1,Length[sl]}]; (*end Do *)
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=List[];    
 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[#,AxesLabel->{xu,"Y"},
    PlotJoined->True,PlotRange->All,
    PlotLabel->StringInsert[tl,ToString[cntr=cntr+1],1],
    PlotStyle->{RGBColor[Part[colr=RotateLeft[colr,3],1],
    Part[colr,2],Part[colr,3]]},
    DisplayFunction->ipt]&,{Transpose[m]}]; 
 If[pn==1, (*replace last point with number*)
 MapThread[AppendTo[pobj[[#,1]],
      Text[#,Last[pobj[[#,1,2,1]]]]]&,
      {Range[Length[Transpose[m]]]}];
   ]; (*end if*)
 If[pm==1,Show[pobj,PlotLabel->tl,
     DisplayFunction->$DisplayFunction]];    
      ]; (*end module*)

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[Part[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},
(*calc number of levels over which sum taken by redspc*)
 inx=binselect[Table[1,{Length[nlv]}]-Part[ptr,2],{Reverse[nlv]}];
 dv=Apply[Times,Flatten[inx]];
 Print["Number of cells summed over=",dv];
 If[0==dv,dv=1];
 ms=redspc[m,Part[ptr,2]];(*sum over means matrix*)
 ms=N[ms/dv]; (*calc average of means for plotting*)
 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; (*find dimensions to sum over in c*)
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]; (*end If *)
 ,{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[#,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: July 23/97

intrplot[m_,st_,nl_]:=Module[{mm,trm,ans,rp,pltrm,
           msrd,flag=0},
 trm=fndintr[st];(*find all 2&3 interaction terms*)
 If[0==Length[trm],Print["No interaction terms"];
      Return[]];
 mm=ndim[m,nl];(*make matrix of means*)
 ans=Transpose[{Map[Part[#,1]&,trm]}];
 Print["\nInteraction terms available=",ans];
While[flag==0,
rp=kbinstrg["type an interaction term to plot",ans];
Print["\nPlot for interaction=",rp];
Do[If[{rp}==ans[[j]],(*get full defn for matsum*)
   pltrm=trm[[j]];
   Break[]]; (*end If*)
,{j,1,Length[ans]}];
(*now determine means matrix and type of plot*)
If[StringLength[rp]==2,
  If[Length[nl]==2, (*only 2 factors ab*)
   prntimm[mm,6,rp];
   plotnLists[Transpose[mm],StringInsert["Interaction:",rp,-1],
    Part[Reverse[Characters[rp]],1],0,1,1], 
   (*else more than 2 factors but 2 factor interaction*)
   msrd=matsum[nl,pltrm,mm];
   prntimm[msrd,6,rp]; (*print matrix of means*)
   plotnLists[Transpose[msrd],StringInsert["Interaction:",rp,-1],
    Part[Reverse[Characters[rp]],1],0,1,1]
    ], (*end If Length[nl]==2*)
 (*else If StringLength[rp]==2*)
  If[Length[nl]==3, (*3 factors and abc type interaction*)
  prntimm[mm,6,rp];
  plt3fintr[mm,rp],
   (*else 4+ factors and 3-term interaction*)
   msrd=matsum[nl,pltrm,mm];
   prntimm[msrd,6,rp];
   plt3fintr[msrd,rp]];
]; (*end If 2 factor interaction*)
If["n"==kbinstrg["More plots? y/n",{{"y"},{"n"}}],
   Return[]];
]; (*end While*)
 ];

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: July 11/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[m,pc,rwch], (*else if 3D*) 
rwch=genrch[StringTake[is,{2}],StringTake[is,{3}],
      Rest[Dimensions[m]]];
Do[
  lv=StringInsert[StringTake[is,1],ToString[j],-1];
  Print["\nLevel=",lv];
  pmat2[Part[m,j],pc,rwch],
  {j,1,Part[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=List[];
Do[AppendTo[rv,Take[temp,p[[j]]]];
   temp=RotateLeft[temp,p[[j]]],
  {j,1,Length[p]}];(*end Do *)
Return[rv];
]; (*end Module*)

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...},....}
Calls: binselect, probF,hytst 
Mod. Date: August 7/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["\nUser Defined Contrasts"];
While[iflg==1,
cntrm=kbinstrg["Type contrast term",Transpose[{at}]];
Print["User selected contrast term=",cntrm];
(*Identify user selection in available list and get
  df for this term*)
bs=Map[#==cntrm&,Flatten[at]] /. {True->1,False->0};
(*use bs to select associated df *)
dfk=Part[Flatten[binselect[bs,{df}]],1];
(*get number of rows and check*)
erf=1;
While[erf==1,
crw=kbin3["How many rows in contrast?",0];
  If[crw>dfk,Print["Error: #rows must be <=",dfk],
     erf=0];
]; (*end While erf==1*)
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];
(*check for linear dependency of rows*)
If[False=== (*returns False if error*)
Map[Table[0,{Part[Dimensions[rwrc],2]}]==#&,rwrc]==
    Table[False,{Part[Dimensions[rwrc],1]}],
  Print["Error:rank<row order"],erf=0];
 ]; (*end While *) 
hytst[cntrm,kcon];(*show hypothesis being tested*) 
(*get sum-reduced form of contrast*)
kct=Transpose[kcon];
kcon=Transpose[Drop[kct+Table[-1 Last[kct],
     {Length[kct]}],-1]];
Print["Sum-reduced contrast matrix is=",kcon];
Print["Row echelon form of contrast=",RowReduce[kcon]];
(*get submatrix of ixr and calc SSk*)
bs=Map[cntrm===inxi[[#]][[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];
(*get appropriate param estimates for contrast;corresponds
  to 1 in bs*)
prr=mefb[[Flatten[binselect[bs,{Range[Length[mefb]]}]]]];
prr=Flatten[prr];
Print["Parameter estimates used=",prr];
sumsk=ssk[kcon,subxi,prr];(*SS for contrast*)
probF[sumsk,Length[kcon],ems,edf];
If["n"==kbinstrg["More contrasts? y/n",{{"y"},{"n"}}],
   iflg=0];
]; (*end While iflg==1*)
Return[];
]; (*end Module*)

hytst[t,k] displays string expression for hypothesis being test 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,Part[Dimensions[k],2]];
stcf=MapThread[ToString[#]&,{prd}];
MapThread[Print[#,"=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=List[];
vn=Range[Apply[Plus,d]];(*gen integer sequence*)
Do[(*partition into units of length df*)
 AppendTo[lse,Take[vn,d[[j]]]];
 vn=RotateLeft[vn,d[[j]]];
,{j,1,Length[d]}];
bgn=Map[Take[#,1]&,lse];(*get start index*)
(*append # of cols and the term*)
indxi=Transpose[Append[{bgn},Transpose[{d}]]];
indxi=Map[Flatten[#,1]&,indxi];
indxi=Map[Prepend[indxi[[#]],at[[#]]]&,Range[Length[at]]];
Return[indxi];
]; (*end Module*)

gencon[s,d,p,xtxi,mse,dfe]: calc SS,MS, F, and prob for each main effect and interaction term and displays result.
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: July 17/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]; (*drop mu term*)
dftrm=Drop[d,1];
par=restrv[Drop[p,1],dftrm];(*restructure parameters*)
lcn=lctmev[dftrm,strm];(*calc location indices*)
di=Length[xtxi]-1;
rxi=SubMatrix[xtxi,{2,2},{di,di}];(*subM excluding mu*)
ct=List[];
Do[(*generate row echelon form of contrast over # terms*)
 con=IdentityMatrix[dftrm[[j]]];
 si=lcn[[j]][[2]];di=lcn[[j]][[3]];
 inxt=SubMatrix[rxi,{si,si},{di,di}];
 sumsq=rndoff[ssk[con,inxt,par[[j]]],4];
 msk=rndoff[sumsq/dftrm[[j]],4];
 If[dfe>=1,
  fr=rndoff[msk/mse,4];
  fd=FRatioDistribution[dftrm[[j]],dfe];
  pr=rndoff[1.-CDF[fd,fr],3],(*else*)
  fr="*";pr="*"]; (*end If *)
 AppendTo[ct,{strm[[j]],sumsq,dftrm[[j]],msk,fr,pr}];
 (*hytst[strm[[j]],con]; display hypothesis*)
,{j,1,Length[strm]}]; 
AppendTo[ct,{"Error",mse dfe,dfe,mse}];
Print["\nAnova 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: kbin3, srtgrp, srtcnt, kbinstrg, dsgnmef, dmpar, dsgnint, cutintr, anvT1, 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: August 7/97

anovanw[d_]:=Module[{ntot,nf,ntrm,nlev,xdsm,xdscode,slctd,sstot,
             means,conv,ncs,xtx,xxi,df,bpar,ssmean,ssmod,stuf,
             temp,mse,dfe,ixr,mef,dfa,xindx,bmef},
Print["GLM Anova V1.0"];
Print["Number of indices defining factor levels=",
      Length[Part[Part[d,1],1]]];
nf=kbin3["Type the number of factors in the design",0];
Print["Number of factors=",nf];
(*determine form of data based on type of input matrix*)
If[True==VectorQ[Part[Part[d,1],2]],
  Print["Data input is in grouped form"];
  {sstot,nlev,ncs,means}=srtgp[d,nf],
  Print["Data input is in ungrouped form"];
  {sstot,nlev,ncs,means}=srtcnt[d,nf];
  ];
ntot=Apply[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+Apply[Plus,nlev-1]];
Print["Interaction parameters(ab,ac,...) estimated=",
    (Apply[Times,nlev])-(1+Apply[Plus,nlev-1])];

If[1>(Apply[Plus,ncs]-Apply[Times,nlev]),
  Print["\nWARNING: 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[]]; (*end If *) 
(*create design matrix for main effects a,b,c,..*)
xdsm=dsgnmef[nlev];
xdscode=dmpar[nf];(*design matrix code structure all terms*)
(*identify main effects as already created*)
slctd=xdscode[[2^(-1+Range[nf])]];
(*now append position index for terms available*)
slctd=Map[Append[slctd[[#]],#]&,Range[Length[slctd]]];
(*slctd contains alfa name,binary repr,posn in xdsm*)
(*now calc interaction vectors*)
{slctd,xdsm}=dsgnint[xdscode,slctd,xdsm];
(*allow user to cut interaction terms*)
If[nf>=2 && "y"==kbinstrg["Remove interaction terms? y/n",
        {{"y"},{"n"}}],
  {xdsm,slctd}=cutintr[xdsm,slctd,nf]];
ntrm=Length[slctd];(*adj # terms after cut*)
df=Map[Take[Dimensions[Part[xdsm,#]],1]&,Range[ntrm]];
(*bring xdsm to 2D matrix form*) 
xdsm=Flatten[xdsm,1];
conv=Table[1,{Apply[Times,nlev]}]; (*constant mu *)
xdsm=Join[{conv},xdsm]; (*join const vector to design m*)
df=Prepend[Flatten[df],1]; (*prpend df for mu*) 
slctd=Prepend[slctd,{"mu",Table[0,{nf}],0}];
(*renumber the position index because constant was added*)
slctd=Map[ReplacePart[slctd[[#]],#,-1]&,Range[Length[slctd]]];
xtx=Map[ncs Part[xdsm,#]&,Range[Length[xdsm]]];
xtx=xtx.Transpose[xdsm];
xxi=Inverse[xtx];
bpar=xxi.xdsm.(means ncs); (*parameter estimates*)

temp=rndoff[restrv[bpar,df],4];
Do[Print["\nTerm=",slctd[[j]][[1]]," df=",df[[j]]];
 Print["Par. Est=",temp[[j]]],
 {j,1,Length[df]}];

(*temp=Map[#[[1]]&,slctd];
stuf=Transpose[{{temp},{df},rndoff[{restrv[bpar,df]},4]}];
Print[TableForm[stuf,
  TableHeadings->{None,{"Terms","df","Paramater Estimates"}}, 
   TableSpacing->{0,2}]];
*)
ssmean=((Apply[Plus,means ncs])^2)/ntot;
ssmod=bpar.xtx.bpar;(*SS model*)
(*calc and display general anova table*)
{mse,dfe}=anvT1[sstot,ssmean,ssmod,{ntot,Length[bpar]}];

(*now do SS for all avail anova terms except mu*)
gencon[slctd,df,bpar,xxi,mse,dfe];
If[dfe==0,Return[]];
If[nf>=2, (*interaction plots*)
  If["y"==kbinstrg["Plot interactions? y/n",
          {{"y"},{"n"}}],
intrplot[means,slctd,nlev]];(*end If "y" *)
];(*end If nf>=2*)
(*User contrasts on main effects*)
If["y"==kbinstrg["Contrasts on main effects? y/n",
         {{"y"},{"n"}}],
(*drop const and interaction terms to get mef containing
  only main effect terms; do same to get dfa*)
mef=Take[Drop[Part[Transpose[slctd],1],1],nf];
dfa=Take[Rest[df],nf];(*df for main effects only*)
(*indices to obtain subm of ixr used in getucn*)
xindx=lctmev[dfa,mef];
(*drop constant term in xxi and bpar*)
ixr=SubMatrix[xxi,{2,2},Table[-1+Length[xxi],{2}]];
temp=Take[Rest[bpar],Apply[Plus,dfa]];
bmef=restrv[temp,dfa];(*restructure vector of est.*)
getucn[mef,dfa,ixr,xindx,mse,dfe,bmef];
]; (*end If contrasts*)
]; (*end Module *)

anovanw[dat2]

anovanw[dat4]