Stefan Th. Gries
2013-Jul-07 20:14 UTC
[R] Hierarchical multi-level model with lmer: why are the highest-level random adjustments 0?
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
Bert Gunter
2013-Jul-07 22:50 UTC
[R] Hierarchical multi-level model with lmer: why are the highest-level random adjustments 0?
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