[R] unbalanced mixed effects models for fully factorial designs
Spencer Graves
spencer.graves at pdf.com
Wed Aug 9 16:44:21 CEST 2006
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.
More information about the R-help
mailing list