[R] unbalanced mixed effects models for fully factorial designs

[Spencer Graves]

Dear Prof. Logan:

	  Are you familiar with 'lme' in the 'nlme' package?  Superlative 
documentation for the 'nlme' package is Pinheiro and Bates (2000) Mixed 
effects models in S and S-PLUS (Springer, listed as 'available' in the 
catalog of your university library).  The standard R distribution comes 
with a directory "~library\nlme\scripts" containing script files 
'ch01.R', 'ch02.R', ..., 'ch06.R', and 'ch08.R'.  These contain R script 
files with the R code for each chapter in the book.  I've learned a lot 
from walking through the script files line by line while reviewing the 
corresponding text in the book.  Doing so protects me from problems with 
silly typographical errors as well as subtle problems where the S-Plus 
syntax in the book gives a different answer in R because of the few 
differences in the syntax between S-Plus and R.

       Hope this helps.
       Spencer Graves

Murray Logan wrote:
> Does anyone know of a way of dealing with unbalanced mixed effects 
> (fixed and random factors) for fully factorial designs.
> 
> An example of such data is given below;
> 
> The response variable is SQRTRECRUITS
> SEASON is a random factor
> DENSITY is a fixed factor
> Thus DENSITY:SEASON is a fixed factor.
> 
> Therefore, whereas the effects of SEASON and DENSITY:SEASON should be 
> tested against the overall residual (error) term, the effect of DENSITY 
> should be tested against the DENSITY:SEASON interaction.
> To complicate matters, the data are unbalanced, and thus Type III SS are 
> preferable
> 
> quinn <-
> structure(list(SEASON = structure(as.integer(c(2, 2, 2, 2, 2,
> 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 1, 1, 1,
> 1, 1, 1, 1, 1, 1, 1, 4, 4, 4, 4, 4, 4, 4, 4, 4)), .Label = c("Autumn",
> "Spring", "Summer", "Winter"), class = "factor", contrasts = "contr.sum"),
>     DENSITY = structure(as.integer(c(2, 2, 2, 2, 2, 1, 1, 1,
>     1, 1, 1, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2,
>     1, 1, 1, 1, 1, 1, 2, 2, 2, 1, 1, 1, 1, 1, 1)), .Label = c("High",
>     "Low"), class = "factor"), RECRUITS = as.integer(c(15, 10,
>     13, 13, 5, 11, 10, 15, 10, 13, 1, 21, 31, 21, 18, 14, 27,
>     34, 49, 69, 55, 28, 54, 14, 18, 20, 21, 4, 22, 30, 36, 13,
>     13, 8, 0, 0, 10, 1, 5, 9, 4, 5)), SQRTRECRUITS = c(3.872983,
>     3.162278, 3.605551, 3.605551, 2.236068, 3.316625, 3.162278,
>     3.872983, 3.162278, 3.605551, 1, 4.582576, 5.567764, 4.582576,
>     4.242641, 3.741657, 5.196152, 5.830952, 7, 8.306624, 7.416198,
>     5.291503, 7.348469, 3.741657, 4.242641, 4.472136, 4.582576,
>     2, 4.690416, 5.477226, 6, 3.605551, 3.605551, 2.828427, 0,
>     0, 3.162278, 1, 2.236068, 3, 2, 2.236068), GROUP = 
> structure(as.integer(c(4,
>     4, 4, 4, 4, 3, 3, 3, 3, 3, 3, 6, 6, 6, 6, 6, 6, 5, 5, 5,
>     5, 5, 5, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 8, 8, 8, 7, 7, 7,
>     7, 7, 7)), .Label = c("AutumnHigh", "AutumnLow", "SpringHigh",
>     "SpringLow", "SummerHigh", "SummerLow", "WinterHigh", "WinterLow"
>     ), class = "factor")), .Names = c("SEASON", "DENSITY", "RECRUITS",
> "SQRTRECRUITS", "GROUP"), row.names = c("1", "2", "3", "4", "5",
> "6", "7", "8", "9", "10", "11", "12", "13", "14", "15", "16",
> "17", "18", "19", "20", "21", "22", "23", "24", "25", "26", "27",
> "28", "29", "30", "31", "32", "33", "34", "35", "36", "37", "38",
> "39", "40", "41", "42"), class = "data.frame")
> 
> I realise that Anova (car package) calculated Type III SS (given the 
> correct contrasts), however, this does not permit mixed models.  
> Conversely, if I was to specify a aov model such as;
> summary(aov(SQRTRECRUITS ~ SEASON+DENSITY+Error(DENSITY:SEASON), 
> data=quinn))
> purely to obtain a test for DENSITY (ignoring the test for SEASON), the 
> SS are Type I.
> 
> Although it is possible to calculate out the F-ratio (and p-value) 
> calculations manually and substitute them into the anova tables, I cant 
> help think that there must be a better solution.
> 
> Is there any expectation that there will be a summary routine that 
> provides Type II and Type II SS, and or is aov ever likely to 
> accommodate non-hierarchical mixed models?
> 
> Regards
> 
> Murray
> 
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.