[R] Hierarchical multi-level model with lmer: why are the highest-level random adjustments 0?

[Bert Gunter]

I think this is much better posted on r-sig-mixed-models rather than r-help.

-- Bert

On Sun, Jul 7, 2013 at 1:14 PM, Stefan Th. Gries <stgries at gmail.com> wrote:
> Hi all
>
> I have a hopefully not too stupid question about multi-level /
> mixed-effects modeling. I was trying to test a strategy from Crawley's
> 2013 R Book on a data set with the following structure:
>
> - dependent variable: CONSTRUCTION (a factor with 2 levels)
> - independent fixed effect: LENGTH (an integer in the interval [1, 61])
> - random effects with the following hierarchical structure: MODE >
> REGISTER > SUBREGISTER > FILE. Specifically:
>
> MODE: S
>   REGISTER: monolog
>     SUBREGISTER: scripted
>     SUBREGISTER: unscripted
>   REGISTER: dialog
>     SUBREGISTER: private
>     SUBREGISTER: public
>   REGISTER: mix
>     SUBREGISTER: broadcast
> MODE: W
>   REGISTER: printed
>     SUBREGISTER: academic
>     SUBREGISTER: creative
>     SUBREGISTER: instructional
>     SUBREGISTER: nonacademic
>     SUBREGISTER: persuasive
>     SUBREGISTER: reportage
>   REGISTER: nonprinted
>     SUBREGISTER: letters
>     SUBREGISTER: nonprofessional
>
> with various levels of FILE in each level of SUBREGISTER. Here's the
> head of the relevant data frame (best viewed with a non-proportional
> font):
>
>   CASE MODE   REGISTER     SUBREGISTER    FILE CONSTRUCTION LENGTH
> 1    1    W    printed       reportage W2C-002 V_P_DO       11
> 2    2    W    printed     nonacademic W2B-035 V_P_DO        8
> 3    3    W nonprinted nonprofessional W1A-014 V_P_DO        8
> 4    4    W    printed       reportage W2C-005 V_P_DO        6
> 5    5    S     dialog         private S1A-073 V_DO_P        2
> 6    6    S     dialog         private S1A-073 V_DO_P        2
>
> And here's the unique-types distribution of FILE in the design:
> tapply(FILE, list(SUBREGISTER, REGISTER, MODE), function (qwe)
> length(unique(qwe)))
>
> , , S
>                 dialog mix monolog nonprinted printed
> academic             .   .       .          .       .
> broadcast            .  20       .          .       .
> creative             .   .       .          .       .
> instructional        .   .       .          .       .
> letters              .   .       .          .       .
> nonacademic          .   .       .          .       .
> nonprofessional      .   .       .          .       .
> persuasive           .   .       .          .       .
> private             96   .       .          .       .
> public              77   .       .          .       .
> reportage            .   .       .          .       .
> scripted             .   .      25          .       .
> unscripted           .   .      66          .       .
>
> , , W
>                 dialog mix monolog nonprinted printed
> academic             .   .       .          .      26
> broadcast            .   .       .          .       .
> creative             .   .       .          .      20
> instructional        .   .       .          .      19
> letters              .   .       .         28       .
> nonacademic          .   .       .          .      37
> nonprofessional      .   .       .         17       .
> persuasive           .   .       .          .       9
> private              .   .       .          .       .
> public               .   .       .          .       .
> reportage            .   .       .          .      20
> scripted             .   .       .          .       .
> unscripted           .   .       .          .       .
>
> # I would usually have done this (using lme4)
> model.1.tog <- lmer(CONSTRUCTION ~ LENGTH +
> (1|MODE/REGISTER/SUBREGISTER), family=binomial)
>
> # but Crawley (2013:692ff.) suggests this:
> LEVEL2 <- MODE:REGISTER
> LEVEL3 <- MODE:REGISTER:SUBREGISTER
> model.1.sep <- lmer(CONSTRUCTION ~ LENGTH + (1|MODE) + (1|LEVEL2) +
> (1|LEVEL3), family=binomial)
>
> The results are the same for fixed and random effects, ok, but what I
> don't understand in this case is why the random adjustments to
> intercepts at the highest level of hierarchical organization (MODE)
> are 0:
>
> ranef(model.1.sep)
> $LEVEL3
>                              (Intercept)
> S:dialog:private              -0.9482442
> S:dialog:public                0.1216021
> [...]
>
> $LEVEL2
>              (Intercept)
> S:dialog      -0.4746389
> [...]
>
> $MODE
>   (Intercept)
> S           0
> W           0
>
> I am guessing that's got something to do with the fact that what the
> random adjustments to intercepts at the level of MODE would do is
> already taken care of by the random adjustments to intercepts on the
> lower levels of the hierarchical organization of the data - but when I
> run the code from Crawley on his data (a linear mixed-effects model,
> not a generalized linear mixed-effects model as mine), the highest
> hierarchical level of organization does not return 0 as random
> adjustment. What am I missing?
>
> Thanks for any input you may have!
> STG
>
> ______________________________________________
> R-help at r-project.org 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.



-- 

Bert Gunter
Genentech Nonclinical Biostatistics

Internal Contact Info:
Phone: 467-7374
Website:
http://pharmadevelopment.roche.com/index/pdb/pdb-functional-groups/pdb-biostatistics/pdb-ncb-home.htm