Murray Logan
2006-Aug-02 22:25 UTC
[R] unbalanced mixed effects models for fully factorial designs
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
Spencer Graves
2006-Aug-09 14:44 UTC
[R] unbalanced mixed effects models for fully factorial designs
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.