[R] unbalanced mixed effects models for fully factorial designs

[Murray Logan]

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