Stephen Gosnell
2011-Aug-21 22:03 UTC
[R] Split split plot analysis for unbalanced data using lmer
Hello, I'm attempting to analyze a split-split plot model, currently using lmer. My design is similar to that used in Casella's Experimental Designs Ozone Example (example 5.7, p 197), but I have been unable to find any coding help for split-split plot models through searches of several textbooks and the R list serve. I have three treatments of interest (A,B,C), each with 2 treatment levels. One level of treatment A (Y or N) was applied to 6 blocks (Block, labelled 1,2,3,4,5,6). Each block was subdivided into 4 subplots (Positon, labelled 1,2,3,4, so implicit nesting may be a problem), to which treatment B (Y or N) was randomly assigned (2 replicates of each level of treatment B per block). Each of these subplots was divided in half and randomly assigned treatment C (B or M). My response variable is the number of predetermined spots (out of 132) occupied by a certain species; I have a unbalanced data set since one observation was not recorded for a total of 47 counts. My data is set up as: Block A Position B C Count 1 1 Y 4 N B 64 2 1 Y 4 N M 109 3 1 Y 3 Y M 88 4 1 Y 3 Y B 0 5 1 Y 1 Y B 47 6 1 Y 1 Y M 91 with both position and block set to factors, not integers. My current model is: model=lmer(cbind(june$Count, 132-june$Count)~A*B*C+(1|A/B)+(1|Block),june, family="binomial") with the resulting summary> modelGeneralized linear mixed model fit by the Laplace approximation Formula: cbind(june$Count, 132 - june$Count) ~ A * B * C + (1 | A/B) + (1 | Block) Data: june AIC BIC logLik deviance 806.9 827.2 -392.4 784.9 Random effects: Groups Name Variance Std.Dev. Block (Intercept) 0.074222 0.27244 B:A (Intercept) 0.000000 0.00000 A (Intercept) 0.000000 0.00000 Number of obs: 47, groups: Block, 6; B:A, 4; A, 2 Fixed effects: Estimate Std. Error z value Pr(>|z|) (Intercept) -0.09452 0.17295 -0.547 0.58470 AY 0.23756 0.24453 0.971 0.33130 BY -0.30949 0.10259 -3.017 0.00255 ** CM -0.57438 0.10434 -5.505 3.70e-08 *** AY:BY 0.64931 0.14922 4.351 1.35e-05 *** AY:CM 1.12220 0.14754 7.606 2.83e-14 *** BY:CM 0.39417 0.14769 2.669 0.00761 ** AY:BY:CM -0.45040 0.21368 -2.108 0.03504 * --- Signif. codes: 0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.? 0.1 ? ? 1 Correlation of Fixed Effects: (Intr) AY BY CM AY:BY AY:CM BY:CM AY -0.707 BY -0.291 0.206 CM -0.286 0.202 0.483 AY:BY 0.200 -0.283 -0.688 -0.332 AY:CM 0.202 -0.286 -0.342 -0.707 0.469 BY:CM 0.202 -0.143 -0.695 -0.706 0.478 0.500 AY:BY:CM -0.140 0.197 0.480 0.488 -0.698 -0.691 -0.691 I'm not sure if I have specified the model correctly (should Position be included in the model specification as well?) or how to interpret the results. I think the group numbers for blocks and A are correct, but I'm not sure how B:A should be checked or if the implicit nesting of B is a problem. Do I need to include factor C in my nesting structure (even though its the lowest level)? My initial conclusions are that Treatments B and C have a main effect and higher order interactions suggest treatment A is also important. All my random effects also seem extremely small, which seems concerning. Any help (or pertinent references) for model setup or interpretation is greatly appreciated. Stephen Gosnell Department of Ecology, Evolution, and Marine Biology, UCSB