[R] Split split plot analysis for unbalanced data using lmer

[Stephen Gosnell]

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

> model
Generalized 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