Millman-Glass Expected Mean Square Heuristic for Anova

Steve Hunka, University of Alberta,
Edmonton, Alberta, Canada. T6G 2G5
steve.hunka@ualberta.ca
Revision Date: July 20/02

References

Kirk, R.E. (1982) Experimental Design: Procedures for the Behavioral Sciences, 2nd Edition, Brooks/Cole Publishing Co., Pacific Grove, CA.

Millman, J. & Glass, G.V., Rules of Thumb for Writing the Anova Table, Jr. of Educ. Meas., Vol. 4, No. 2, Summer, 1967

Winer, B.J., Statistical Principles in Experimental Design, 2nd Ed., McGraw-Hill, N.Y., 1962

Introduction

The Millman-Glass heuristic for calculating the symbolic representation of expected mean squares is described in their paper entitled "Rules of Thumb for Writing the Anova Table" in the Journal of Educational Measurement, Vol. 4, No. 2, Summer 1967.   In Anova, knowing the symbolic representation of expected mean squares is important in determining the mean square that should be used in the denominator of the ratio of mean squares used in the F test. This is especially true when one or more of the factors in the Anova are random or nested. Expected mean squares are also used in the area of generalizability. The heuristic also provides a representation of the sum of squares that can be used as a basis for combining sum of  squares from a competely crossed  dessign to carry out nested or repeated measure designs.

The heuristic described by Millman and Glass involves assigning a single and unique character representation to each factor. A factor that is nested is identified by having a colon (:) followed by the factor or factors within which it is nested.  A set of rules is given in terms of the character representation, which when applied sequentially, provide the derivation of the expected mean squares.

Implementation  Notes

The implementation of the rules involves the manipulation of the characters representing factors, their interactions, and their nesting.  For example, one rule states that any combination of the factors which has the same character to the left and right of the colon is deleted. Another rule (V-B(2)) states that an EMS term "having to the left of the colon a letter denoting a fixed classification disappears except when the source of variation of the mean square includes this letter".

Although Mathematica has very good string manipulation operations, it was decided that for the major portion of the implementation, string manipulations would not be used, but instead, a binary representation of the terms would be used. The reason for using a binary representation is that Boolean operations can be more easily and efficiently performed than manipulation of strings of characters.  

The following data structure is used.

	{u, {v},{{w}}}
	
where

u	is the full character representation of a term, e.g., ab:cd
v	is the binary representation of the term  (characters left of the colon) e.g., ab is 
	{0,0,1,1}
w	is the binary representation of each character to the right of the colon (the
	nesting terms), e.g., {{0,1,0,0},{1,0,0,0}} indicates nesting under c and d.

Notice that the data structure involves characters, integers (1,0), a vector (v), and a matrix (w). Such a structure can be joined to other structures to represent a full range  of terms that are required. 

Using the above structure, our example would be

	{ab:cd,{0,0,1,1},{{0,1,0,0},{1,0,0,0}}
	
and to check whether the same character appears to the left and right of the colon we would

i) do a Boolean addition of the columns of {{0,1,0,0},{1,0,0,0}} to get {{1,1,0,0}} and 

ii) test whether the vector {0,0,1,1} is  logically included in the Boolean sum {1,1,0,0}...it is not, thus there are no characters to the left of the colon which also appear to the right of the colon. 

If necessary, the structure {ab:cd,{0,0,1,1},{{0,1,0,0},{1,0,0,0}} can be "taken apart". Letting z={ab:cd,{0,0,1,1},{{0,1,0,0},{1,0,0,0}}

Part[z,1] is ab:cd
Part[z,2] is {0,0,1,1}
Part[z,3] is  {{0,1,0,0},{1,0,0,0}}

If necessary, ab:cd can be further "taken apart" by

Characters[Part[z,1]] giving {a,b,:,c,d}.

The generation of all combinations of the factors is accomplished by generating the base two representation of the integers from 1 to ( (2^nf)-1) where nf is the number of factors, e.g., for nf=3 we get {0,0,1}, {0,1,0}, {0,1,1}, {1,0,0},{1,0,1},{1,1,0}, {1,1,1}, in which a is represented by {0,0,1}, b by {0,1,0}, ab by {0,1,1}, and so on. (In our example the characters a, b, c, and d are assigned from right to left in the binary vector, i.e., the low order bit represents the character a.)

The data structure described above is maintained almost entirely throughout the computations  so that any intermediate computations can be easily checked and interpreted.

User Input

1. The heuristic is started by executing the function ems3[nf] where nf is the number (integer) of factors in the design including a factor for persons (observations, replications).
2. The user is queried for a definition of each factor. For a factor without nesting a single unique character is expected (not followed by a colon); for a nested  factor a unique character, followed by a colon, followed by one character for each factor under which the nesting exists. (There is no space between the characters which define factors under which nesting occurs.)  A random factor must be defined for persons (observations or replications) and is considered nested under all other factors. Example: r:pct could indicate a  factor of "replications" nested under factors p, c, and t. 
3. For each factor the user is asked to input whether it is fixed (1) or random (0). The "replications" factor is considered to be random.
4. For each factor the user is asked to input the number of levels (integer) in the factor. For the "replications" factor, the number of levels equals the number of observations within a cell (equal n). A factor which is nested must have the same number of levels in each level of the factor in which it is nested.

The user input parameters are summarized in the output.

Note: No check is made of the integrity of the above input.

Output

1)  Factors as input by the user.
2) A table showing the relationship between each factor as input and its binary representation, and the binary representation of the nesting factors. (Not important to intrepetation of output). 
3) A table showing for each factor (Source) whether it is fixed or random, and the number  of levels.
4) A table showing for each factor and  any interactions, the degrees of freedom, the coefficient associated with each expected mean square, and the expected mean squares.
5) When df=0 is found, the associated source and ems are removed and a new table displayed.

Millman-Glass Example

In this example, the factors entered are p, c:p, t, and r:pct. In the Millman-Glass heuristic one factor is always defined as a "replication factor"  which is always random and nested in all other factors (Millman and Glass, p43) and has as many levels as observations in each cell.  The user is also asked to specify which factors are fixed or random, and the levels in each factor. The user is prompted for the information with reference to the factors p, c, t, and r. The result is as follows.

ems3[4]

Interpretation of Output

The "Factors Defined" show the characters representing each factor and its nesting as input by the user. The table following shows the relationship between each main effect factor (Source) term and the binary representation of the term followed by the binary representation of the nesting terms. The binary nesting vectors are  held in matrix form indicated by the enclosure of the vectors between {  and }, e.g.,{{0,0,0,0}} is a  matrix having one row vector {0,0,0,0}. The nesting vector {0,0,0,0} indicates no nesting. This table is not important to interpretation  of the last table.

The last  table gives the source terms and their nesting terms ( source terms to the left  of the colon and nesting terms to the right)  and the degrees of freedom. The heading "coef" is to be interpreted  as  follows:

1) The coefficient for p: is 30.  All terms p: under EMS have the coefficient 30,  i.e., 30p: 
2) The coefficent for c:p is 10. All terms c:p  under EMS have the coefficient 10, i.e., 10c:p
3) The coefficient for r:pct is 1.

Thus, the first line under EMS is read as follows:

30(Subscript[(((σ^2))), p]) + 10(Subscript[(((σ^2))), c : p]) + 1(Subscript[(((σ^2))), r : pct])

In many textbooks  the term (Subscript[(((σ^2))), r : pct])  is represented as (((Subscript[σ, ϵ])^2) ) in fully crossed designs.

EMS Terms used for F

In the above example, in order to get the F ratio for testing the significance of factor "p" we would use the mean square for p in the numerator and the mean square for c:p in the denominator. Notice that the terms under the column EMS for factor p are (p:, c:p, r:pct) and for c:p it is (c:p, r:pct) and the common elements are (c:p and r:pct) leaving the unique element p in the numerator. Under the null hypothesis that p is not significant, the ratio of expected mean squares would be 1. Thus, to test a source factor one must have a denominator containing all terms in the numerator except the source term.

The term r:pct (replications nested under  factors p, c, and t) would in the simplest of anovas be interpreted as the error term for testing all other terms. In this example the mean square of r:pct would only be used to test c:p and ct:p.

Cautions

No checks are included to ensure integrity of the information input by the user.

Function List and Description

ems3		the main function 
adnstg2	calculates binary representation of nesting factors of main terms
entr4		prepares structure for table of terms and binary representations
dmpar2	calculates all  combinations of factors and binary representations
bin2al4		translates binary representation to factor character representation
nesting2	calculates  binary nesting representation for all combination of terms
rmvdup2	removes duplicate binary nesting vectors
rl1		delete term if the same  character  appears on each side of the colon
rl22		delete an ems if it does not have all the terms of the source term.
inclqa		tests logical inclusiveness for two vectors
inclqa2		tests  logical inclusiveness  for two vectors
boladd		boolean addition of 2 vectors
grteq		test vector for elements >=c

getrf2		user input for fixed or random factors
rl32		 removes terms from pt which have fixed factors to the left of the colon except 		when the fixed factor is the source term. Rule V-B(2)
bolsub		boolean subtraction
getlv		user input of levels in each factor
degf		calculate degrees of freedom
fcns		puts together characters of an EMS based on binary terms
coef		calculates coefficient for each EMS
ent2,3		used to structure output for tables
dtdf0		delete  source and ems  terms where  df=0

Functions  that require initialization

dmpar2[nf,fn] returns the character and binary designation  for all possible combination of factors in the form {a,{001}, b,{010},....}. nf is the number of factors and fn is the list of the characters representing the factors. Leading 0s are prepended.
Calls: bin2alf4 to get character designations.
Returns: { {factor character, {binary representation}},......}
Mod. Date: Jan 19/02

dmpar2[nf_, fn_] := 
        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 = 
            bin2al4[bin2, fn]; 
bin4 = 
            Transpose[Append[{bin3}, bin2]]; Return[bin4];];

Example of dmpar2

dmpar2[3, {"a", "b", "c"}]

bin2al4[b,a] selects chars in a for 1 in vectors of b
b: binary vectors {{001},{010},...{111}}
a: {"a","b","c"}
Mod.Date:Jan 19/02

bin2al4[b_, a_] := 
        Module[{bp, ap}, 
bp = 
            Map[Reverse[#] &, b]; 
bp = 
            Map[Position[#, 1] &, bp]; 
bp = 
            Map[Flatten[#, 1] &, bp]; 
ap = 
            Map[a[[#]] &, bp]; 
ap = 
            Map[StringJoin[#] &, ap]; 
Return[
            ap];
];

adnstg2[f1] takes f1 and appends the binary representation of the nesting characters to the main terms as input by the user.
f1 is of the form: {{a,{0001}},{b,{0010}},{c:b,{0100}},{p:a,{1000}}} main terms only
Returns: {{a,{0001},{0000}},{b,{0010},{0000},{c:b,{0100},{0010},{p:a,{1000},{0001}}}i.e., for each part, Part 1 is the main effect  term, Part 2 is its binary representation, Part 3 is the binary representation of the nesting factor(s).
Mod. Date: Jan 12/02

\(\(adnstg2[f1_] := 
        Module[{nl, nf, fc, f2, nc, main, cp, b0, indx, bn, 
            newl}, 
nl = List[]; \ 
          nf = Length[f1]; 
f2 = 
            Map[Part[#, 2] &, f1];  (*binary\ repr*) 
fc = 
            Map[Part[#, 1] &, 
              f1];  (*all\ chars\ terms*) 
fc = 
            Characters[
              fc]; 
 (*get\ first\ char*) 
main = 
            Map[Part[Part[fc, #], 1] &, 
              Range[nf]]; 

 (*in\ fc\ get\ position\ of\ colon\ into\ cp; \ {}\ if\ \(\(no\)\(\ \)\(:\)\)*) 
cp = 
            Flatten[Map[Map[Position[#, ":"] &, {Part[fc, #]}] &, 
                Range[nf]], 1]; 
 (*\ 
            for\ each\ term\ get\ nesting\ in\ binary*) 
b0 = {Table[0, {nf}]};  (*for\ no\ nesting*) 
Do[\(If[
                Length[cp[\([j]\)]] == 0, 
                nl = Append[nl, 
                    b0], \  (*else*) 
 (*get\ nesting\ characters\ only*) 
nc = 
                  Drop[fc[\([j]\)], 
                    Part[Flatten[cp[\([j]\)]], 1]]; 
nc = 
                  Flatten[
                    Characters[
                      nc]]; 
 (*find\ indx\ in\ f1\ and\ assign\ binary\ v*) 
indx = 
                  Flatten[
                    Map[Position[main, #] &, nc]]; 
bn = 
                  f2[\([indx]\)]; \ 
nl = 
                  Append[nl, 
                    bn];
];\)\ \  (*end\ if*) 
, {j, 1, nf}]; \  (*end\ Do*) 
newl = 
            MapThread[Append[#1, #2] &, {f1, nl}]; 
Return[
            newl];
];\)\)

Example admstg2

adnstg2[{{"a", {0, 0, 1}}, {"b", {0, 1, 0}}, {"c:ab", {1, 
            0, 0}}}]

entr4[t] arranges terms for printed output; t is of the form {{p,{0,0,0,1},{{0,0,0,0}}},{c:p,{0,0,1,0,{{,0,0,1}}}....} i.e, all input terms with binary representation and binary nesting representation
Mod. Date: Jan 12/02

entr4[t_] := {Part[t, 1], ToString[Part[t, 2]], 
        ToString[Part[t, 3]]}

nesting2[t,at] takes t as the definition of a main term and its nestings and appends the nesting binary vectors to all combination of terms in at that contains that term
t: is of the form: {{a,{0001},{{0,0,0}}}
at: contains all combination of terms and the binary representations but no nesting vector
Calls: inclqa
Returns: all combinations with those terms matching t containing nesting vectors.
Mod. Date: Jan 12/02

nesting2[t_, at_] := 
        Module[{bt, bn, new, at2}, 
bt = 
            Part[t, 2];   bn = Part[t, 3];  at2 = at; 
          Do[
If[
                1 == 
                  inclqa[bt, Part[at2[[k]], 2]], 
new = 
                  Join[at2[[k]][[3]], 
                    bn]; 
at2[[k]] = 
                  ReplacePart[at2[[k]], new, 
                    3]];
, {k, 1, 
              Length[at]}];    
Return[
            at2];
];

inclqa[a,b]; if vector a is included in vector b, result is True=1, else False=0.
Calls: boladd,greteq

inclqa[a_, b_] := 
      b == boladd[a, b] /. {True -> 1, False -> 0}

boladd[a_,b_]; Boolean addition of vectors a and b. If a and b are matrices, Flatten[a], Flatten[b], add them, then restructure to original form with Partition.
Calls: grteq
Mod. Date: Nov 5/97

boladd[a_, b_] := grteq[a + b, 1];

grteq[x,c] takes the vector x and returns a vector containing 1 if x[[i]]>=c and 0 otherwise
Created: May 26/98

grteq[x_, c_] := 
      Map[GreaterEqual[#, c] &, x] /. {True -> 1, False -> 0}

inclqa2[a,b]; if vector b is included in vector a, result is True=1, else False=0;
Calls: boladd,greteq

inclqa2[a_, b_] := 
      a == boladd[a, b] /. {True -> 1, False -> 0}

rmvdup2[bn] removes duplicate binary nesting vectors from bn. 
Mod. Date: Jan 23/02

rmvdup2[bn_] := 
        Module[{nf, dt, sb, ubnt}, 
nf = 
            Length[Part[bn, 1]]; 
dt = 
            Union[Map[FromDigits[#, 2] &, bn]]; 
sb = 
            Map[IntegerDigits[#, 2] &, dt]; 
ubnt = 
            Map[Join[Table[0, {nf - Length[#]}], #] &, 
              sb]; 
Return[
            ubnt];
];

Example rmvdup2

rmvdup2[{{1, 0, 1}, {1, 0, 1}, {0, 0, 1}, {0, 1, 0}, {0, 0, 
          1}}]

rl1[trm] if a term has the same character on both sides of a colon remove the term. trm is of the form {bp,{1,0,1,0),{{0,0,0,0},{0,0,0,1}}}, where Part 2 is the binary rep of the term and Part 3 is the binary rep of the nesting factors. If a nesting vector is logically included in the binary repr of the term, the term is removed.
Calls: inclqa
Mod. Date: Jan 12/02

rl1[trm_] := 
        Module[{bt, nt, zv, res}, 
bt = Part[trm, 2]; 
          nt = Part[trm, 3]; 
zv = 
            Table[0, {Length[bt]}]; 
  
If[
            1 == Length[nt]  &&   nt == {zv}, 
            Return[trm]]; 
If[Part[nt, 1] == zv, 
            nt = Drop[nt, 1]]; 
res = 
            Map[inclqa[#, bt] &, nt]; 
If[
            Apply[Plus, res] > 0, Return[], 
            Return[trm]];
];

rl22[pl] Rule V-B1:delete an ems if it does not have all the terms of the source term.
 Used to remove from the potential list those terms which do not contain the characters of jth term.
 pl: a structure containing all retained terms
 Note: pl is a list of retained  ems terms, but are also the Source terms. So each term (Source) is checked against the remaining terms, e.g., Source term1 vs terms,2,3,4..., then Source term 2 vs terms 1,3,4.... etc.  For this reason the output is longer than the input pl 
 Returns: a structure in which each part contains the source term and its retained ems terms.
Calls: grteq, inclqa
Mod. Date: Jan 12/02

Note: The binary term and binary nesting term of each term are Boolean summed for all  terms which gives all the characters  included in each Source term. If a term's boolean sum is logically included in the source sum, then that term is retained with the Source term, else it is dropped.

rl22[pl_] := 
        Module[{u, v, w, w1, pt, srtrm, z, 
            res}, 
  
u = Map[Join[{Part[#, 2]}, Part[#, 3]] &, pl]; 
          v = Map[Apply[Plus, #] &, u]; 
w = 
            Map[grteq[#, 1] &, v]; 
res = 
            List[]; 
Do[
            srtrm = pl[[j]]; 
pt = 
              Delete[pl, {{j}}]; 
w1 = 
              Delete[w, {{j}}]; 
z = 
              Map[inclqa[w[[j]], #] &, 
                w1]; 
  
v = 
              Delete[pt, Position[z, 0]]; 
If[
              0 ≠  Length[v], srtrm = Append[{srtrm}, v]; 
              res = Append[res, srtrm], 
res = 
                Append[res, {srtrm}]];
, {j, 1, 
              Length[pl]}]; 
Return[
            res];
];

getrf2[tbt] gets user input fixed/random for each factor
tbt: of the form {a,{0,0,1}
Returns: {a,{0,0,1},1 or 0}
Mod. Date: Jan 12/02

getrf2[tbt_] := 
        Module[{rorf, 
            ques}, 
ques = "Is  fixed (1) or random \
(0)?"; 
ques = 
            StringInsert[ques, ToString[Part[tbt, 1]], 
              4]; 
rorf = 
            Append[tbt, Input[ques]]; 
Return[
            rorf];
];

rl32[pt,ft] removes terms from pt which have fixed factors to the left of the colon except when the fixed factor is the source term. Rule V-B(2)
pt is of the form {{{a,{0,0,1},{{0,0,0}}},{ab,{0,1,1},{{0,0,0}}},....} a potential list for term "a". Part 1 is source and potential ems, Part 2 remaining potential ems 
 ft is list  of fixed  factors in binary of the form {{0,0,1},{0,1,0},...}
 Note: bolsub removes source term from remaining pt terms then a check is made if fixed factors are in the remaining binary terms. Called by Map[#,for1]&,x]
 Calls: bolsub, grteq, inclq2a
 Mod. Date: Jan 12/02

rl32[pt_, ft_] := 
        Module[{srct, srct1, rem, bt, dif, u, 
            x}, 

srct1 = 
            Part[pt, 1];   
srct = 
            Part[srct1, 2];   
rem = 
            Flatten[Rest[pt], 1]; 
If[
            Length[rem] == 0, x = Flatten[pt, 1];  
            Return[x]]; 
bt = 
            Map[Part[#, 2] &, 
              rem]; 
  
dif = 
            Map[bolsub[#, srct] &, bt]; 
  
fnin[
              u_, v_] := Map[inclqa2[#, v] &, u]; 
u = 
            Map[fnin[dif, #] &, ft]; 
u = 
            Apply[Plus, u];   
u = 
            grteq[u, 1]; 
u = 
            Position[u, 1]; 
If[0 == Length[u], 
            x = Append[srct1, rem], 
x = 
              Append[srct1, Delete[rem, u]]]; 
Return[
            x];
];

bolsub[a,b] boolean subtraction by taking a-b and changing all -1s to 0

bolsub[a_, b_] := (a - b) /. {-1 -> 0}

getlv[tbt] gets user input of the number of levels for each factor
tbt: main input terms {{p,{0001},{c,{0010}},.....}
Returns:{{p,{0001},k}, {c,{0010},k}, ...}where k is the number of levels
Mod. Date: Jan 20/02

\(\(getlv[tbt_] := 
        Module[{rorf, 
            ques}, 
ques = "Number of Levels in "; 
ques = 
            StringInsert[ques, 
              ToString[Part[tbt, 1]], \(-1\)]; 
rorf = 
            Append[tbt, Input[ques]]; 
Return[
            rorf];
];\)\)

degf[x,d] calculates degree of freedom for all source terms
x Part 1 of retained list of terms {a,{001},{{0,0,0}}.....} for which df required; it  is the first part of the retained ems  terms.
d vector of df for each main effect factor defined, eg., for a, b, p, so reversed to fit binary term  sequence
Note: for factors to  left  of : product  of (levels-1); to  the right of : product of levels. The levels are  selected by the position of 1 in bt and nt.
Mod. Date: Jan 12/02

degf[x_, d_] := 
        Module[{bt, nt, dr, dfm1, bs, mlt, df}, 
bt = 
            Part[x, 2];  
          nt = Part[x, 3]; 
dr = 
            Reverse[d]; 
dfm1 = -1 + 
              Delete[dr, Position[bt, 0]]; 
dfm1 = 
            Apply[Times, 
              dfm1]; 
  
  
bs = Apply[Plus, nt]; 
bs = 
            Map[GreaterEqual[#, 1] &, bs] /. {True -> 1, 
                False -> 
                  0}; 
  
mlt = 
            Apply[Times, 
              Delete[dr, Position[bs, 0]]]; 
df = 
            dfm1*mlt; 
Return[
            df];
];

fcns[mt,mc] for the defn of a source of variation mt, and mc the main factor designated  characters,  returns  the characters representing the source of variation with nesting characters included. Called through Map[Map[...]] to provide all the EMS terms in character form with colon for each source term.
mt: of the form {term, {binary term},{{nest1},{nest2}...}}
mc:  of the form {p,c,t,r}
Returns: {term:nesting terms}; string vector
Mod. Date: Jan 12/02

fcns[mt_, mc_] := 
        Module[{trm, nt, bs, revbs, 
            nsc}, 

trm = 
            StringInsert[
              Part[mt, 1], ":", -1]; 
nt = 
            Part[mt, 3];   
bs = 
            Apply[Plus, nt]; 
bs = 
            Map[GreaterEqual[#, 1] &, bs] /. {True -> 1, 
                False -> 0}; 
revbs = 
            Reverse[bs]; 
nsc = 
            Delete[mc, Position[revbs, 0]]; 
If[
            0 == Length[nsc], Return[trm]]; 
nsc = 
            StringJoin["", Sequence[nsc]]; 
trm = 
            StringInsert[trm, nsc, -1]; 
Return[
            trm];
];

coef[w,lv]  calculate coefficient; w is a source term, lv vector of levels of each factor
w is of the form {"p",{0,0,0,1},{{0,0,0,0}}}. Note: Apply[Times,{}] returns 1. The product of the levels NOT associated with the ems is calc.
Returns:  coefficient for p
Mod. Date: Jan 12/02

coef[w_, lv_] := 
        Module[{bm, bs, c}, 
bm = 
            Join[Part[w, 3], {Part[w, 2]}]; 
bs = 
            Apply[Plus, bm]; 
bs = 
            Map[GreaterEqual[#, 1] &, bs] /. {True -> 1, 
                False -> 
                  0}; 
  
bs = 
            Map[Equal[1, #] &, bs] /. {True -> 0, 
                False -> 1}; 
bs = 
            Reverse[bs]; 
c = 
            Delete[lv, Position[bs, 0]]; 
c = 
            Apply[Times, c];  Return[c];
];

entr2[t,d] used to create rows of the table output having Source terms, df, coeficients,  and EMS terms
t: structure of retained terms {{terms},{terms},..}; strings
d: deg freedom for each retained term; vector
c: coefficients for each retained term; vector
Mod. Date: Jan 12/02

entr2[t_, d_, c_] := {Part[t, 1], ToString[d], 
        ToString[c], {t}}

entr3[t,d] structures t and d to output factor name, Type (fixed or random) and number of levels
t: has the form {{"a",{0,0,1},1},{"b",{0,1,0},0},....} in which Part 3 is 1=fixed, 0=random
d: vector of the number of levels in each factor.
Mod. Date: Jan 12/02

entr3[t_, d_] := 
        Module[{r, typ}, 
If[1 == Part[t, 3], 
            typ = "Fixed", 
            typ = "Random"]; 
r = {Part[t, 1], typ, 
              ToString[d]}; 
Return[
            r];
];

dtdf0[a,di] drops source terms and related ems terms in a as indexed by di
Mod. Date: Jan 26/02

dtdf0[a_, di_] := 
        Module[{da, dems, dinx, 
            res}, 

da = 
            Delete[a, di]; 
dems = 
            Map[Part[#, 1] &, 
              a[[Flatten[
                    di]]]]; 
  
funm[u_, e_] := 
            Map[Position[u, #] &, e]; 
dinx = 
            Map[funm[#, dems] &, da]; 
dinx = 
            Map[Flatten[#, 1] &, dinx]; 
res = 
            MapThread[Delete, {da, dinx}]; 
Return[
            res];
];

Main Starting Function to Initialize

ems3[nf] starts the calculations; nf is the number (integer) of factors including a factor for replications.
Calls: adnstg, dmpar2, nesting2, rl1, rl22, rl32, rmvdup2, getrf2, degf, fcns, coef, dtdf0, entr2, getlv
Mod. Date: Jan 26/02

\(\(ems3[nf_] := 
        Module[{ques, quesj, brep, bcblist, fctlist, fctlist2, mt, mtc, allt, 
            u, w, x, for, for1, y, w1, df, df1, dfn, alt, 
            cef}, 
Print["Expected Mean Squares"]; 
ques = "Define factor "; 
fctlist = 
            List[]; 
bcblist = 
            List[]; 
brep = 
            Join[Table[
                0, {nf - 
                    1}], {1}]; 

 (*\(get\ user\ defn\ of\ factors\  &\)\ gen\ binary\ rep*) 
Do[
            quesj = StringInsert[ques, 
                ToString[j], \(-1\)]; 
bcblist = 
              Append[bcblist, brep]; 
brep = 
              RotateLeft[brep]; 
fctlist = 
              Append[fctlist, InputString[quesj]];
, {j, 1, 
              nf}]; \  (*end\ Do*) 
Print["Factors Defined=", fctlist]; 
fctlist2 = 
            Transpose[
              Append[{fctlist}, 
                bcblist]]; 
 (*append\ binary\ repr\ of\ nesting\ terms*) 
fctlist = 
            adnstg2[fctlist2]; 
Print["Input Terms and Binary Representation"]; 
Print[
            TableForm[Map[entr4[#] &, fctlist], 
              TableHeadings -> {{}, {"Source", "Source Term", "Nesting Term"}}]]; 

 (*now\ pick\ off\ main\ terms\ only\ no\ nesting*) 
mt = 
            Map[Part[#, 1] &, fctlist]; 
mtc = 
            Map[Part[#, 1] &, 
              Characters[
                mt]]; 
 (*generate\ all\ comb\ of\ main\ terms*) 
allt = 
            dmpar2[nf, 
              mtc]; 
 (*now\ assign\ init\ binary\ nesting\ vector*) 
allt = 
            Map[Append[#, {Table[0, {nf}]}] &, allt]; 
 (*\ 
            append\ binary\ nesting\ vectors\ for\ each\ term*) 
Do[
allt = 
              nesting2[fctlist[\([j]\)], allt]
, {j, 1, 
              Length[fctlist]}]; 
 (*next\ remove\ duplicate\ binary\ nesting\ vectors*) 
u = 
            Flatten[Map[Take[#, \(-1\)] &, allt], 
              1]; \ 
w = 
            Map[rmvdup2[#] &, u];  (*remove\ \ dup*) 
x = 
            MapThread[
              Append[Take[#1, 2], #2] &, {allt, 
                w}]; 
 (*x\ \(holds : \ term\), bin\ rep, 
            nesting\ list*) 

 (*Rule\ II - B \((2  
                    c)\)\ remove\ if\ same\ char\ each\ side\ of\ colon*) 
u = 
            Map[rl1[#] &, 
              x]; 
 (*remove\ Null\ from\ the\ list\ of\ terms*) 
x = 
            Cases[u, 
              x_ /; x =!= 
                  Null]; 
 (*ems\ must\ have\ all\ chars\ in\ source*) 
x = 
            rl22[x];  (*rule\ V - 
              B1*) 

 (*get\ fixed/
                random\ input*) 
w1 = 
            MapThread[
              Append[{#1}, #2] &, {mtc, bcblist}]; 
for = 
            Map[getrf2[#] &, 
              w1]; 

 (*get\ fixed\ factors\ in\ binary\ terms*) 
w = 
            Map[Equal[Part[#, 3], 1] &, for]; 
for1 = 
            Delete[for, Position[w, False]]; 
for1 = 
            Map[Part[#, 2] &, 
              for1]; 

 (*now\ adj\ for\ fixed\ factors\ Rule\ V - B \((2)\)\ using\ results\ from\ rl22\ V - 
              B1*) 
y = 
            Map[rl32[#, for1] &, 
              x]; 

 (*get\ number\ of\ Levels\ in\ each\ main\ factor*) 
df1 = 
            Map[getlv[#] &, 
              w1]; 
 (*for\ df\ we\ use\ only\ the\ first\ 3\ parts\ of\ elements\ in\ y\ i . e . , \ 
            source\ term*) 

w = 
            Map[Take[#, 3] &, y]; 
dfn = 
            Map[Part[#, 3] &, 
              df1];  (*levels*) 
 (*calc\ deg\ of\ freedom*) 
df = 
            Map[degf[#, dfn] &, 
              w]; 

Print[
            TableForm[MapThread[entr3, {for, dfn}], 
              
              TableHeadings -> {{}, {"Source", "Type", "Levels"}}]]; 
 (*get\ all\ retained\ terms\ for\ each\ source*) 
 (*restructure\ y*) 
y = 
            Map[Join[{Take[#, 3]}, Flatten[Drop[#, 3], 1]] &, 
              y]; 
 (*get\ char\ \(term : 
                nesting\ term\)*) 
alt = 
            Map[Map[fcns[#, mtc] &, #] &, 
              y]; 

 (*calc\ coefficients*) 
cef = Map[coef[#, dfn] &, w]; 
Print[
            TableForm[MapThread[entr2, {alt, df, cef}], 
              TableHeadings -> {{}, {"Source", "df", "coef", "EMS"}}]]; 
 (*delete\ ems\ when\ df = 
              0*) 
u = 
            Position[df, 0]; 
If[0 == Length[u], 
            Return[], \  (*else*) 
Print["Deleting terms where df=0"]; 
alt = 
              dtdf0[alt, u]; 
df = 
              Delete[df, u]; 
cef = 
              Delete[cef, u]]; 
Print[
            TableForm[MapThread[entr2, {alt, df, cef}], 
              TableHeadings -> {{}, {"Source", "df", "coef", "EMS"}}]];

];\)\)

Five Examples

ems3[4]

Winer P 380

ems3[4]

Note to above: 
1) Refer to Table 5.16-1 in Winer. Notice there is no denominator term for testing factor A. Winer discusses testing and dropping certain interaction terms resulting in a model in which factors A, B, C, and AB are testable after dropping terms AC, BC, and ABC. Dropping these terms implies their SS are pooled with the SS of P:ABC.
2) The term  p:abc is given in Winer as (((Subscript[σ, ϵ])^2) ).

ems3[
      4]

ems3[4]

Notes to above: 
With the exception of p:abc having 1  level (1  observation per cell) as  given by Winer, the other levels were arbitrarily defined for the above example to give numerical results.

Winer (p 366, Table 5.12-5) has an additional symbolic term  (((Subscript[σ, ϵ])^2) ) in each row of EMS.  This term does not exist as an independent error term as in the second last table output above  because it has a  df=0. The formation of the various F ratios based on the last  table  are  ok. For example, "a" is tested using a denominator  mean  square for b:a, and c and ac are tested using a denominator mean square bc:a. 

Winer notes (page 367) that the above result is identical for a split plot design having factor A as  treatments, Factor B as plots, and  C as subplots.

ems3[4]

Repeated Measure Winer p526

ems3[4]