rm(list=ls(all.names=TRUE))
set.seed(sum(utf8ToInt("Gummibaerchen")))
pkgs2load <- c("dplyr", "car", "effects", "emmeans", "lme4", "magrittr", "multcomp", "MuMIn")
invisible(sapply(pkgs2load, library, character.only=TRUE, logical.return=TRUE, quietly=TRUE))
rm(pkgs2load)
source("_helpers/cohens.kappa.r"); source("_helpers/C.score.r")Ling 204: session 07
1 Session 07: Mixed-effects models 3a
In this session, we will look at a data set that involves a generalized linear mixed-effects model. Specifically, we will look at a data set that is also discussed in the book, but we will analyze the data differently from how they are analyzed in the book. The data are on subject realization in spoken conversational Japanese involving native and non-native speakers to test the hypothesis that the intermediately proficient non-native speakers in the sample data will react differently to two other main predictors governing subject realization, namely the givenness of the subject referent and whether it is used contrastively (this study is based on a version of the data in Gries & Adelman 2014).
We load from _input/subjreal.csv, and you can find information about its contents in _input/subjreal.r:
CASE SPKR SPEAKERID SUBJREAL CONTRAST
Min. : 1 njs:7371 26-JE_NNS: 821 no :5016 no :6393
1st Qu.: 3814 nns:7394 19-JC_NNS: 769 yes :1649 yes : 271
Median : 7634 18-JK_NJS: 760 NA's:8100 NA's:8101
Mean : 7630 7-JE_NJS : 731
3rd Qu.:11443 24-JE_NNS: 710
Max. :15265 1-JC_NJS : 705
(Other) :10269
GIVENNESS
Min. : 0.000
1st Qu.: 8.000
Median :10.000
Mean : 7.875
3rd Qu.:10.000
Max. :10.000
NA's :8990 1.1 Exploration & preparation
As always, it is indispensable to do a decent amount of exploratory plotting before any statistical analysis is begun. If you do mixed-effects modeling, you really must look at
- every variable in isolation (the response, predictors, controls, random effects, everything);
- every rhs variable (fixed and random effects) with the response;
- every combination of rhs variables that will be considered in the model on their own and with the response.
1.1.1 The response
We begin with the response variable, which exhibits a fairly strong preference for ‘no subject’:
1.1.2 The fixed-effects structure
Here’s the first predictor, SPKR, which is nearly perfectly balanced …
$Freqs
.
njs nns Sum
7371 7394 14765
$Props
.
njs nns
0.4992211 0.5007789 … and doesn’t seem to be correlated much with the response:
Here’s the second predictor CONTRAST, which is distributed extremely unevenly …
$Freqs
.
no yes Sum
6393 271 6664
$Props
.
no yes
0.95933373 0.04066627 … and does seem to be extremely strongly correlated with the response:
And the last predictor, GIVENNESS:
.
0 1 2 3 4 5 6 7 8 9 10
1023 12 13 27 24 37 56 77 192 472 3842 
This predictor is distributed somewhat annoyingly: many 0s, many many 10s, very little in between. At the same time, the result is very suggestive (and in the expected direction).
A quick look at the combination of predictors, which looks like there will be a strong interaction between CONTRAST and GIVENNESS:
SUREAL no yes
CONTR GIV
no 0 0.227 0.773
1 0.333 0.667
2 0.545 0.455
3 0.318 0.682
4 0.667 0.333
5 0.514 0.486
6 0.615 0.385
7 0.743 0.257
8 0.781 0.219
9 0.815 0.185
10 0.881 0.119
yes 0 0.160 0.840
1 NaN NaN
2 0.000 1.000
3 0.000 1.000
4 0.333 0.667
5 0.000 1.000
6 0.500 0.500
7 0.000 1.000
8 0.357 0.643
9 0.053 0.947
10 0.163 0.837But let’s also quickly check out the distribution of the NAs:
SREAL FALSE TRUE
CONTR GIV
FALSE FALSE 5775 0
TRUE 889 0
TRUE FALSE 0 0
TRUE 1 8100Let’s make sure we don’t run into problems later and reduce the dataset to the complete cases:
1.1.3 The random-effects structure
As it should be: Every speaker contributes a nicely large number data points and every speaker is a native speaker of Japanese (njs) or not (nns):
njs nns Sum
1-JC_NJS 299 0 299
1-JC_NNS 0 266 266
10-JE_NJS 247 0 247
10-JE_NNS 0 236 236
11-JE_NJS 195 0 195
11-JE_NNS 0 250 250
16-JE_NJS 359 0 359
16-JE_NNS 0 189 189
18-JK_NJS 386 0 386
18-JK_NNS 0 77 77
19-JC_NJS 160 0 160
19-JC_NNS 0 290 290
2-JK_NJS 222 0 222
2-JK_NNS 0 278 278
24-JE_NJS 234 0 234
24-JE_NNS 0 173 173
25-JE_NJS 300 0 300
25-JE_NNS 0 220 220
26-JE_NJS 179 0 179
26-JE_NNS 0 248 248
7-JE_NJS 277 0 277
7-JE_NNS 0 188 188
8-JE_NJS 318 0 318
8-JE_NNS 0 184 184
Sum 3176 2599 5775But this table reveals something else: the levels of SPEAKERID seem to indicate the native speaker and the non-native speaker who made up one conversational pair. That means, we could add a column that represents which conversation the two speakers in it were a part of:
We can add this is an additional random-effects level, making this a good example of a multi-level model!
What about the predictors? Every speaker produced at least one instance of each level of CONTRAST: we need CONTRAST|SPEAKERID:
no yes
1-JC_NJS 283 16
1-JC_NNS 242 24
10-JE_NJS 241 6
10-JE_NNS 222 14
11-JE_NJS 188 7
11-JE_NNS 247 3
16-JE_NJS 349 10
16-JE_NNS 180 9
18-JK_NJS 374 12
18-JK_NNS 74 3
19-JC_NJS 152 8
19-JC_NNS 260 30
2-JK_NJS 212 10
2-JK_NNS 259 19
24-JE_NJS 220 14
24-JE_NNS 162 11
25-JE_NJS 289 11
25-JE_NNS 210 10
26-JE_NJS 178 1
26-JE_NNS 244 4
7-JE_NJS 268 9
7-JE_NNS 178 10
8-JE_NJS 306 12
8-JE_NNS 181 3And because of the nesting of SPEAKERID into CONV, that obviously means that each conversation also has at least one level of CONTRAST, which means we need CONTRAST|CONV/SPEAKERID:
no yes
1-JC 525 40
10-JE 463 20
11-JE 435 10
16-JE 529 19
18-JK 448 15
19-JC 412 38
2-JK 471 29
24-JE 382 25
25-JE 499 21
26-JE 422 5
7-JE 446 19
8-JE 487 15In addition, every speaker produced multiple different values of GIVENNESS: we need GIVENNESS|SPEAKERID:
0 1 2 3 4 5 6 7 8 9 10
1-JC_NJS 61 1 0 1 2 3 2 7 11 27 184
1-JC_NNS 50 2 0 3 2 2 2 6 9 20 170
10-JE_NJS 44 0 0 4 1 3 3 6 17 29 140
10-JE_NNS 56 1 0 4 1 4 3 4 10 20 133
11-JE_NJS 30 0 0 1 1 0 1 4 5 20 133
11-JE_NNS 41 1 0 0 1 3 3 3 11 22 165
16-JE_NJS 71 0 2 2 1 3 3 4 12 24 237
16-JE_NNS 32 1 1 0 1 0 3 4 11 24 112
18-JK_NJS 88 0 1 0 2 4 7 4 7 34 239
18-JK_NNS 25 0 0 1 0 1 2 0 6 4 38
19-JC_NJS 23 0 0 1 0 1 1 2 4 10 118
19-JC_NNS 65 0 0 2 0 0 2 2 8 22 189
2-JK_NJS 31 2 0 1 2 1 3 1 1 11 169
2-JK_NNS 55 0 1 1 1 0 0 1 5 15 199
24-JE_NJS 43 1 1 0 0 2 0 1 6 24 156
24-JE_NNS 37 0 0 0 0 0 1 1 5 14 115
25-JE_NJS 45 0 1 1 1 0 1 7 15 26 203
25-JE_NNS 31 0 0 0 2 0 0 1 5 25 156
26-JE_NJS 19 0 1 0 0 3 3 2 6 13 132
26-JE_NNS 29 1 1 0 1 1 2 4 5 18 186
7-JE_NJS 35 1 3 1 1 2 5 6 8 16 199
7-JE_NNS 46 1 0 2 1 1 4 1 10 22 100
8-JE_NJS 39 0 1 2 2 3 4 5 8 14 240
8-JE_NNS 27 0 0 0 1 0 1 1 7 18 129And because of the nesting of SPEAKERID into CONV, that obviously means that each conversation also has at least one level of GIVENNESS, which means we need GIVENNESS|CONV/SPEAKERID:
0 1 2 3 4 5 6 7 8 9 10
1-JC 111 3 0 4 4 5 4 13 20 47 354
10-JE 100 1 0 8 2 7 6 10 27 49 273
11-JE 71 1 0 1 2 3 4 7 16 42 298
16-JE 103 1 3 2 2 3 6 8 23 48 349
18-JK 113 0 1 1 2 5 9 4 13 38 277
19-JC 88 0 0 3 0 1 3 4 12 32 307
2-JK 86 2 1 2 3 1 3 2 6 26 368
24-JE 80 1 1 0 0 2 1 2 11 38 271
25-JE 76 0 1 1 3 0 1 8 20 51 359
26-JE 48 1 2 0 1 4 5 6 11 31 318
7-JE 81 2 3 3 2 3 9 7 18 38 299
8-JE 66 0 1 2 3 3 5 6 15 32 3691.2 Modeling
1.2.1 Exploring the random-effects structure 1
Let’s identify the random-effects structure we can/should use. The first model one might fit is this, which will take a long time (because of the complexity of the model and the high number of iterations we permit the model to explore):
summary(m_nest_gxg <- glmer(
SUBJREAL ~ 1 + # intercept
CONTRAST*GIVENNESS*SPKR + # fixefs
(1+CONTRAST*GIVENNESS|CONV/SPEAKERID), # ranefs
family=binomial, data=d, na.action=na.exclude,
control=glmerControl(optCtrl = list(maxfun=30000))),
correlation=FALSE)Generalized linear mixed model fit by maximum likelihood (Laplace
Approximation) [glmerMod]
Family: binomial ( logit )
Formula: SUBJREAL ~ 1 + CONTRAST * GIVENNESS * SPKR + (1 + CONTRAST *
GIVENNESS | CONV/SPEAKERID)
Data: d
Control: glmerControl(optCtrl = list(maxfun = 30000))
AIC BIC logLik -2*log(L) df.resid
4850.0 5036.5 -2397.0 4794.0 5747
Scaled residuals:
Min 1Q Median 3Q Max
-3.3302 -0.4069 -0.3403 0.3502 3.6251
Random effects:
Groups Name Variance Std.Dev. Corr
SPEAKERID:CONV (Intercept) 0.041847 0.20456
CONTRASTyes 0.351754 0.59309 0.74
GIVENNESS 0.002043 0.04520 -0.98 -0.61
CONTRASTyes:GIVENNESS 0.007391 0.08597 0.00 -0.65 -0.20
CONV (Intercept) 0.215186 0.46388
CONTRASTyes 0.156127 0.39513 0.36
GIVENNESS 0.003378 0.05812 -0.93 -0.68
CONTRASTyes:GIVENNESS 0.009689 0.09843 0.75 0.88 -0.94
Number of obs: 5775, groups: SPEAKERID:CONV, 24; CONV, 12
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.266770 0.184426 6.869 6.48e-12 ***
CONTRASTyes 0.179662 0.467829 0.384 0.701
GIVENNESS -0.337782 0.025257 -13.374 < 2e-16 ***
SPKRnns 0.178030 0.181324 0.982 0.326
CONTRASTyes:GIVENNESS 0.375988 0.072267 5.203 1.96e-07 ***
CONTRASTyes:SPKRnns 0.582911 0.666633 0.874 0.382
GIVENNESS:SPKRnns 0.002105 0.026813 0.078 0.937
CONTRASTyes:GIVENNESS:SPKRnns 0.013189 0.100435 0.131 0.896
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1A few months ago, this thing didn’t actually converge unproblematically – now it does. But …
1.2.2 Exploring the fixed-effects structure 1
Let’s see what the fixed-effects structure turns out to be:
npar AIC LRT Pr.Chi.
<none> NA 4849.951 NA NA
CONTRAST:GIVENNESS:SPKR 1 4850.037 2.086235 0.1486319The 3-way interaction seems to be unnecessary, but we have the singularity warning. Again, welcome to the beautiful hell of mixed-effects modeling …
1.2.3 Exploring the random-effects structure 2
Let’s try all simpler random-effects structures:
m_nest_gpg <- update(m_nest_gxg, .~. - (1+CONTRAST*GIVENNESS|CONV/SPEAKERID)
+ (1+CONTRAST+GIVENNESS|CONV/SPEAKERID))
m_nest_g <- update(m_nest_gxg, .~. - (1+CONTRAST*GIVENNESS|CONV/SPEAKERID)
+ (1+ GIVENNESS|CONV/SPEAKERID))
m_nest_c <- update(m_nest_gxg, .~. - (1+CONTRAST*GIVENNESS|CONV/SPEAKERID)
+ (1+CONTRAST |CONV/SPEAKERID))
m_nest_i <- update(m_nest_gxg, .~. - (1+CONTRAST*GIVENNESS|CONV/SPEAKERID)
+ (1 |CONV/SPEAKERID))
m_spkr_gpg <- update(m_nest_gxg, .~. - (1+CONTRAST*GIVENNESS|CONV/SPEAKERID)
+ (1+CONTRAST+GIVENNESS| SPEAKERID))
m_spkr_g <- update(m_nest_gxg, .~. - (1+CONTRAST*GIVENNESS|CONV/SPEAKERID)
+ (1+ GIVENNESS| SPEAKERID))
m_spkr_c <- update(m_nest_gxg, .~. - (1+CONTRAST*GIVENNESS|CONV/SPEAKERID)
+ (1+CONTRAST | SPEAKERID))
m_spkr_i <- update(m_nest_gxg, .~. - (1+CONTRAST*GIVENNESS|CONV/SPEAKERID)
+ (1 | SPEAKERID))Note: these results may change over time as packages get updated. These are the results for lme4 version 1.1-38 and that’s the reason – again – why all your scripts should end the way mine do: with a call to sessionInfo(). Here are the results in a summary form:
-
m_nest_gxg: fine, butdrop1led to singularity; -
m_nest_gpg: boundary (singular) fit: see help(‘isSingular’); -
m_nest_g: Model failed to converge with max|grad| = 0.010574; -
m_nest_c: Model failed to converge with max|grad| = 0.00344046; -
m_nest_i: fine; -
m_spkr_gpg: boundary (singular) fit: see help(‘isSingular’); -
m_spkr_g: Model failed to converge with max|grad| = 0.00705867; -
m_spkr_c: fine; -
m_spkr_i: Model failed to converge with max|grad| = 0.00259271.
Let’s see how much the models that converged unproblematically differ in their predictions from each other:
FALSE TRUE
FALSE 4545 3
TRUE 3 1224[1] 0.9756487
FALSE TRUE
FALSE 4541 7
TRUE 3 1224[1] 0.9830091
FALSE TRUE
FALSE 4542 6
TRUE 2 1225[1] 0.9944702Not very much really, at least that is somewhat encouraging.
1.2.4 Exploring the fixed-effects structure 2
We look at models without problems:
npar AIC LRT Pr.Chi.
<none> NA 4859.313 NA NA
CONTRAST:GIVENNESS:SPKR 1 4857.837 0.5240882 0.469103 npar AIC LRT Pr.Chi.
<none> NA 4853.277 NA NA
CONTRAST:GIVENNESS:SPKR 1 4851.476 0.1987317 0.6557465We go with m_spkr_c since its simplification led to no warning:
summary(m_spkr_c_02 <- update(m_spkr_c, .~. # update this to a model w/
- CONTRAST:GIVENNESS:SPKR), correlation=FALSE) # the 3-way interaction droppedGeneralized linear mixed model fit by maximum likelihood (Laplace
Approximation) [glmerMod]
Family: binomial ( logit )
Formula: SUBJREAL ~ CONTRAST + GIVENNESS + SPKR + (1 + CONTRAST | SPEAKERID) +
CONTRAST:GIVENNESS + CONTRAST:SPKR + GIVENNESS:SPKR
Data: d
Control: glmerControl(optCtrl = list(maxfun = 30000))
AIC BIC logLik -2*log(L) df.resid
4851.5 4918.1 -2415.7 4831.5 5765
Scaled residuals:
Min 1Q Median 3Q Max
-3.1717 -0.3916 -0.3610 0.4335 3.1688
Random effects:
Groups Name Variance Std.Dev. Corr
SPEAKERID (Intercept) 0.05051 0.2248
CONTRASTyes 0.67407 0.8210 -0.86
Number of obs: 5775, groups: SPEAKERID, 24
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.190536 0.121891 9.767 <2e-16 ***
CONTRASTyes 0.160592 0.406949 0.395 0.693
GIVENNESS -0.325755 0.012390 -26.291 <2e-16 ***
SPKRnns 0.102070 0.176862 0.577 0.564
CONTRASTyes:GIVENNESS 0.340109 0.039606 8.587 <2e-16 ***
CONTRASTyes:SPKRnns 0.501306 0.522611 0.959 0.337
GIVENNESS:SPKRnns 0.008028 0.017958 0.447 0.655
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 npar AIC BIC logLik X.2.log.L. Chisq Df Pr..Chisq.
m_spkr_c_02 10 4851.476 4918.089 -2415.738 4831.476 NA NA NA
m_spkr_c 11 4853.277 4926.551 -2415.639 4831.277 0.1987317 1 0.6557465Ok, looks good, no issues, so let’s drop stuff:
npar AIC LRT Pr.Chi.
<none> NA 4851.476 NA NA
CONTRAST:GIVENNESS 1 4911.142 61.6665090 4.068417e-15
CONTRAST:SPKR 1 4850.446 0.9702189 3.246256e-01
GIVENNESS:SPKR 1 4849.675 0.1996845 6.549756e-01ARGGHHHH, we get a convergence warning, but one that is so small that we will ignore it:
summary(m_spkr_c_03 <- update(m_spkr_c_02, .~. # update this to a model w/
- GIVENNESS:SPKR), correlation=FALSE) # this interaction droppedGeneralized linear mixed model fit by maximum likelihood (Laplace
Approximation) [glmerMod]
Family: binomial ( logit )
Formula: SUBJREAL ~ CONTRAST + GIVENNESS + SPKR + (1 + CONTRAST | SPEAKERID) +
CONTRAST:GIVENNESS + CONTRAST:SPKR
Data: d
Control: glmerControl(optCtrl = list(maxfun = 30000))
AIC BIC logLik -2*log(L) df.resid
4849.7 4909.6 -2415.8 4831.7 5766
Scaled residuals:
Min 1Q Median 3Q Max
-3.1582 -0.3899 -0.3627 0.4268 3.1496
Random effects:
Groups Name Variance Std.Dev. Corr
SPEAKERID (Intercept) 0.05049 0.2247
CONTRASTyes 0.66742 0.8170 -0.86
Number of obs: 5775, groups: SPEAKERID, 24
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.16434 0.10651 10.932 <2e-16 ***
CONTRASTyes 0.16728 0.40553 0.413 0.680
GIVENNESS -0.32209 0.00925 -34.823 <2e-16 ***
SPKRnns 0.16006 0.12019 1.332 0.183
CONTRASTyes:GIVENNESS 0.33998 0.03960 8.585 <2e-16 ***
CONTRASTyes:SPKRnns 0.48408 0.52030 0.930 0.352
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
optimizer (Nelder_Mead) convergence code: 0 (OK)
Model failed to converge with max|grad| = 0.00201685 (tol = 0.002, component 1)
See ?lme4::convergence and ?lme4::troubleshooting. npar AIC BIC logLik X.2.log.L. Chisq Df Pr..Chisq.
m_spkr_c_03 9 4849.675 4909.627 -2415.838 4831.675 NA NA NA
m_spkr_c_02 10 4851.476 4918.089 -2415.738 4831.476 0.1996845 1 0.6549756Anything else to be dropped?
npar AIC LRT Pr.Chi.
<none> NA 4849.675 NA NA
CONTRAST:GIVENNESS 1 4909.337 61.66155 4.078685e-15
CONTRAST:SPKR 1 4848.589 0.91341 3.392113e-01Yes and amazingly enough ….
summary(m_spkr_c_04 <- update(m_spkr_c_03, .~. # update this to a model w/
- CONTRAST:SPKR), correlation=FALSE) # this interaction droppedGeneralized linear mixed model fit by maximum likelihood (Laplace
Approximation) [glmerMod]
Family: binomial ( logit )
Formula: SUBJREAL ~ CONTRAST + GIVENNESS + SPKR + (1 + CONTRAST | SPEAKERID) +
CONTRAST:GIVENNESS
Data: d
Control: glmerControl(optCtrl = list(maxfun = 30000))
AIC BIC logLik -2*log(L) df.resid
4848.6 4901.9 -2416.3 4832.6 5767
Scaled residuals:
Min 1Q Median 3Q Max
-2.9599 -0.3914 -0.3620 0.4260 3.1602
Random effects:
Groups Name Variance Std.Dev. Corr
SPEAKERID (Intercept) 0.0515 0.2269
CONTRASTyes 0.6350 0.7968 -0.87
Number of obs: 5775, groups: SPEAKERID, 24
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.133759 0.102217 11.092 <2e-16 ***
CONTRASTyes 0.395916 0.327583 1.209 0.2268
GIVENNESS -0.322140 0.009251 -34.821 <2e-16 ***
SPKRnns 0.222996 0.102688 2.172 0.0299 *
CONTRASTyes:GIVENNESS 0.338731 0.039519 8.571 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 npar AIC BIC logLik X.2.log.L. Chisq Df Pr..Chisq.
m_spkr_c_04 8 4848.589 4901.879 -2416.294 4832.589 NA NA NA
m_spkr_c_03 9 4849.675 4909.627 -2415.838 4831.675 0.91341 1 0.3392113… this one converges unproblematically. Are we done now? Looks like it, which would be interesting …
npar AIC LRT Pr.Chi.
<none> NA 4848.589 NA NA
SPKR 1 4850.863 4.273919 3.870157e-02
CONTRAST:GIVENNESS 1 4908.154 61.564927 4.283816e-15Wow, we are: SPKR is significant:
1.3 Numerical interpretation
Let’s put everything together in one place, this time also with (just the fastest) CIs:
m_final %>% {
list( "Summary"=summary(., correlation=FALSE),
"P-values"=drop1(., test="Chisq") %>% data.frame,
"ConfInts"=Confint(m_final)) }$Summary
Generalized linear mixed model fit by maximum likelihood (Laplace
Approximation) [glmerMod]
Family: binomial ( logit )
Formula: SUBJREAL ~ CONTRAST + GIVENNESS + SPKR + (1 + CONTRAST | SPEAKERID) +
CONTRAST:GIVENNESS
Data: d
Control: glmerControl(optCtrl = list(maxfun = 30000))
AIC BIC logLik -2*log(L) df.resid
4848.6 4901.9 -2416.3 4832.6 5767
Scaled residuals:
Min 1Q Median 3Q Max
-2.9599 -0.3914 -0.3620 0.4260 3.1602
Random effects:
Groups Name Variance Std.Dev. Corr
SPEAKERID (Intercept) 0.0515 0.2269
CONTRASTyes 0.6350 0.7968 -0.87
Number of obs: 5775, groups: SPEAKERID, 24
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.133759 0.102217 11.092 <2e-16 ***
CONTRASTyes 0.395916 0.327583 1.209 0.2268
GIVENNESS -0.322140 0.009251 -34.821 <2e-16 ***
SPKRnns 0.222996 0.102688 2.172 0.0299 *
CONTRASTyes:GIVENNESS 0.338731 0.039519 8.571 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
$`P-values`
npar AIC LRT Pr.Chi.
<none> NA 4848.589 NA NA
SPKR 1 4850.863 4.273919 3.870157e-02
CONTRAST:GIVENNESS 1 4908.154 61.564927 4.283816e-15
$ConfInts
Estimate 2.5 % 97.5 %
(Intercept) 1.1337588 0.93341774 1.3340999
CONTRASTyes 0.3959162 -0.24613452 1.0379670
GIVENNESS -0.3221404 -0.34027260 -0.3040082
SPKRnns 0.2229958 0.02173139 0.4242603
CONTRASTyes:GIVENNESS 0.3387314 0.26127535 0.4161874The overall significance test of m_final can again be conducted via model comparison: we compare m_final to a model that has the same ranef structure but no fixefs but the overall intercept; in addition, we compute m_final’s relative likelihood over the null model:
m_null_re <- glmer( # create a null model, now with ranefs:
SUBJREAL ~ 1 + # no fixefs, an overall intercept plus
(1 + CONTRAST | SPEAKERID), # the ranef structure of the final model
family=binomial, data=d)
anova(m_final, # LRT comparison of the final model
m_null_re, # with this null model
test="Chisq") %>% data.frame # !! npar AIC BIC logLik X.2.log.L. Chisq Df Pr..Chisq.
m_null_re 4 6423.130 6449.775 -3207.565 6415.130 NA NA NA
m_final 8 4848.589 4901.879 -2416.294 4832.589 1582.541 4 0# pchisq(1582.541, 4, lower.tail=FALSE)
LMERConvenienceFunctions::relLik( # how much more likely
m_final, # m_final is the right model than
m_null_re) # m_null_re (i.e., insanely likely) AIC(x) AIC(y) diff relLik
4848.589 6423.130 -1574.541 Inf [1] InfLet’s now turn to the predictions. First, we compute fixed-effects predictions for all combinations of predictors:
nd <- expand.grid(
GIVENNESS=0:10, # all values of GIVENNESS
CONTRAST =levels(d$CONTRAST), # all levels of CONTRAST
SPKR=levels(d$SPKR)) # all levels of SPKR
nd <- cbind(nd, PRED_LO_2=predict( # add predicted log odds
m_final, re.form=NA, newdata=nd)) # (for now w/out ranefs)
nd <- cbind(nd, PRED_PP_2=plogis(nd$PRED_LO_2)) # add predicted probs
(nd <- nd[,c(3,2,1,4,5)]) # reorder columns to something better SPKR CONTRAST GIVENNESS PRED_LO_2 PRED_PP_2
1 njs no 0 1.13375883 0.7565319
2 njs no 1 0.81161843 0.6924543
3 njs no 2 0.48947802 0.6199835
4 njs no 3 0.16733761 0.5417371
5 njs no 4 -0.15480279 0.4613764
6 njs no 5 -0.47694320 0.3829742
7 njs no 6 -0.79908361 0.3102216
8 njs no 7 -1.12122401 0.2457843
9 njs no 8 -1.44336442 0.1910249
10 njs no 9 -1.76550483 0.1461022
11 njs no 10 -2.08764523 0.1103035
12 njs yes 0 1.52967505 0.8219588
13 njs yes 1 1.54626604 0.8243738
14 njs yes 2 1.56285702 0.8267629
15 njs yes 3 1.57944801 0.8291263
16 njs yes 4 1.59603900 0.8314641
17 njs yes 5 1.61262999 0.8337762
18 njs yes 6 1.62922098 0.8360629
19 njs yes 7 1.64581197 0.8383242
20 njs yes 8 1.66240296 0.8405603
21 njs yes 9 1.67899394 0.8427713
22 njs yes 10 1.69558493 0.8449572
23 nns no 0 1.35675466 0.7952317
24 nns no 1 1.03461425 0.7378095
25 nns no 2 0.71247385 0.6709476
26 nns no 3 0.39033344 0.5963630
27 nns no 4 0.06819303 0.5170417
28 nns no 5 -0.25394737 0.4368522
29 nns no 6 -0.57608778 0.3598333
30 nns no 7 -0.89822819 0.2894147
31 nns no 8 -1.22036859 0.2278716
32 nns no 9 -1.54250900 0.1761708
33 nns no 10 -1.86464940 0.1341620
34 nns yes 0 1.75267088 0.8522894
35 nns yes 1 1.76926186 0.8543659
36 nns yes 2 1.78585285 0.8564181
37 nns yes 3 1.80244384 0.8584462
38 nns yes 4 1.81903483 0.8604503
39 nns yes 5 1.83562582 0.8624306
40 nns yes 6 1.85221681 0.8643872
41 nns yes 7 1.86880780 0.8663203
42 nns yes 8 1.88539879 0.8682300
43 nns yes 9 1.90198977 0.8701166
44 nns yes 10 1.91858076 0.8719801The probably most interpretable display:
, , SPKR = njs
CONTRAST
GIVENNESS no yes
0 0.7565319 0.8219588
1 0.6924543 0.8243738
2 0.6199835 0.8267629
3 0.5417371 0.8291263
4 0.4613764 0.8314641
5 0.3829742 0.8337762
6 0.3102216 0.8360629
7 0.2457843 0.8383242
8 0.1910249 0.8405603
9 0.1461022 0.8427713
10 0.1103035 0.8449572
, , SPKR = nns
CONTRAST
GIVENNESS no yes
0 0.7952317 0.8522894
1 0.7378095 0.8543659
2 0.6709476 0.8564181
3 0.5963630 0.8584462
4 0.5170417 0.8604503
5 0.4368522 0.8624306
6 0.3598333 0.8643872
7 0.2894147 0.8663203
8 0.2278716 0.8682300
9 0.1761708 0.8701166
10 0.1341620 0.8719801We also generate the predictions for all cases in our data, first based on fixed and random effects …
d$PRED_PP_2_ae <- predict( # make d$PRED_PP_2_ae the result of predicting
m_final, # from m_final
type="response") # predicted probabilities of "yes"
d$PRED_CAT_ae <- factor( # make d$PRED_CAT_ae the result of checking
ifelse(d$PRED_PP_2_ae>=0.5, # if the predicted prob. of "yes" is >=0.5
levels(d$SUBJREAL)[2], # if yes, predict "yes"
levels(d$SUBJREAL)[1]), # otherwise, predict "no"
levels=levels(d$SUBJREAL)) # make sure both levels are attested… and then only on fixed effects:
d$PRED_PP_2_fe <- predict( # make d$PRED_PP_2_fe the result of predicting
m_final, # from m_final
type="response", # predicted probabilities of "yes"
re.form=NA) # use no ranefs
d$PRED_CAT_fe <- factor( # make d$PRED_CAT_fe the result of checking
ifelse(d$PRED_PP_2_fe>=0.5, # if the predicted prob. of "yes" is >=0.5
levels(d$SUBJREAL)[2], # if yes, predict "yes"
levels(d$SUBJREAL)[1]), # otherwise, predict "no"
levels=levels(d$SUBJREAL)) # make sure both levels are attestedLooks like the ranefs do indeed not make that much of a difference here:
Let’s generate and evaluate a confusion matrix (based on the fixed effects only!) as well as compute C and Cohen’s κ:
(c_m <- table( # confusion matrix: cross-tabulate
"OBS" =d$SUBJREAL, # observed orders in the rows
"PREDS"=d$PRED_CAT_fe)) # predicted orders in the columns PREDS
OBS no yes
no 3920 273
yes 629 953c( # precisions & accuracies/recalls for each level of the response
"Overall acc." =mean(d$SUBJREAL==d$PRED_CAT_fe, na.rm=TRUE),
"Prec. for yes" =c_m["yes","yes"] / sum(c_m[ ,"yes"]),
"Acc./rec. for yes"=c_m["yes","yes"] / sum(c_m["yes",]),
"Prec. for no" =c_m[ "no", "no"] / sum(c_m[ , "no"]),
"Acc./rec. for no"=c_m[ "no", "no"] / sum(c_m[ "no",]),
"C" =C.score(d$SUBJREAL, d$PRED_PP_2_fe),
"Cohen's kappa" =cohens.kappa(c_m)[[1]]) Overall acc. Prec. for yes Acc./rec. for yes Prec. for no
0.8438095 0.7773246 0.6024020 0.8617279
Acc./rec. for no C Cohen's kappa
0.9348915 0.8028429 0.5777748 Not too bad, actually … Finally, lets compute R2-values, which are ‘fine’; the random-effects structure really doesn’t do much:
Computing Tjur’s R2 might also be good, it’s a very intuitive measure:
1.4 Visual interpretation
1.4.1 Fixed effects
We begin with the main effect of SPKR that was actually not significant in the last model: The result is that there is a weak tendency for non-native speakers to realize the subject more often than the native speakers:
sp_d <- data.frame(
sp <- effect("SPKR", m_final))[,-3]
sp_d$fit_lo <- qlogis(sp_d$fit) # add predicted log odds
plot(sp, type="response", ylim=c(0, 1), grid=TRUE)
We continue with CONTRAST:GIVENNESS:
cogi_d <- data.frame(
cogi <- effect("CONTRAST:GIVENNESS", m_final,
xlevels=list(GIVENNESS=0:10)))[,c(2,1,3,5,6)]
cogi_d$fit_lo <- qlogis(cogi_d$fit) # add predicted log odds
plot(cogi, type="response", ylim=c(0, 1), grid=TRUE,
multiline=TRUE, confint=list(style="auto"))
A very clear finding:
- when
CONTRASTis yes,GIVENNESSseems to have no effect on subject realization; - when
CONTRASTis no,GIVENNESSseems to have a strong effect on subject realization such that, the more given the referent of the subject, the less likely it will be realized.
1.4.2 Random effects
Let’s quickly plot the intercept and slope adjustments, first as simple histograms (and I only do one bin width here for each because different bin widths don’t change anything here so I can save):
par(mfrow=c(1,2))
ranef(m_final)$SPEAKERID[,1] %>% hist(main="Intercept adjustments", xlim=c(-1.5, 1.5))
ranef(m_final)$SPEAKERID[,2] %>% hist(main="CONTRAST adjustments" , xlim=c(-1.5, 1.5))
par(mfrow=c(1,1))
Looks fine, I suppose. What about the lattice plots? The maybe only interesting observation is that conversation 10-JE sticks out a bit:
1.5 Exercises
1.5.1 Coefficients (yes, again)
The interpretation of the fixef interaction seemed very clear, but here are some test and follow-up questions because, after all, this is 204 >:) Here’s what you will probably want to work with:
(Intercept) CONTRASTyes GIVENNESS
1.1337588 0.3959162 -0.3221404
SPKRnns CONTRASTyes:GIVENNESS
0.2229958 0.3387314 , , SPKR = njs
CONTRAST
GIVENNESS no yes
0 1.1337588 1.529675
1 0.8116184 1.546266
2 0.4894780 1.562857
3 0.1673376 1.579448
4 -0.1548028 1.596039
5 -0.4769432 1.612630
6 -0.7990836 1.629221
7 -1.1212240 1.645812
8 -1.4433644 1.662403
9 -1.7655048 1.678994
10 -2.0876452 1.695585
, , SPKR = nns
CONTRAST
GIVENNESS no yes
0 1.35675466 1.752671
1 1.03461425 1.769262
2 0.71247385 1.785853
3 0.39033344 1.802444
4 0.06819303 1.819035
5 -0.25394737 1.835626
6 -0.57608778 1.852217
7 -0.89822819 1.868808
8 -1.22036859 1.885399
9 -1.54250900 1.901990
10 -1.86464940 1.918581And now the questions:
- What is the slope of
GIVENNESSwhenSPKRis njs and- when
CONTRASTis no? - when
CONTRASTis yes?
- when
[1] -0.3221404[1] -0.3221404[1] -0.3221404
General Linear Hypotheses
Linear Hypotheses:
Estimate
1 == 0 -0.3221[1] 0.01659099[1] 0.01659099[1] 0.01659099
General Linear Hypotheses
Linear Hypotheses:
Estimate
1 == 0 0.01659# both
emtrends( # compute the trends/slopes
object=m_final, # for the model m_final
specs= ~ CONTRAST, # compute comparisons over this predictor
var="GIVENNESS") # namely all slopes of this predictor CONTRAST GIVENNESS.trend SE df asymp.LCL asymp.UCL
no -0.3221 0.00925 Inf -0.3403 -0.3040
yes 0.0166 0.03840 Inf -0.0588 0.0919
Results are averaged over the levels of: SPKR
Confidence level used: 0.95 - What is the slope of
GIVENNESSwhenSPKRis nns and- when
CONTRASTis no? - when
CONTRASTis yes?
- when
[1] -0.3221404[1] -0.3221404[1] 0.01659099[1] 0.01659099Why are the two identical? Because there’s no 3-way interaction CONTRAST:GIVENNESS:SPKR in m_final!
- What are the confidence intervals of these two slopes (not corrected for ‘multiple testing’)?
$confint
Estimate lwr upr
1 -0.3221404 -0.3402726 -0.3040082
attr(,"conf.level")
[1] 0.95
attr(,"calpha")
[1] 1.959964$confint
Estimate lwr upr
1 0.01659099 -0.05876573 0.09194771
attr(,"conf.level")
[1] 0.95
attr(,"calpha")
[1] 1.959964emtrends( # compute the trends/slopes
object=m_final, # for the model m_final
specs= ~ CONTRAST, # compute comparisons over this predictor
var="GIVENNESS") %>% # namely all slopes of this predictor
data.frame # but return a data frame CONTRAST GIVENNESS.trend SE df asymp.LCL asymp.UCL
1 no -0.32214041 0.009251289 Inf -0.34027260 -0.30400821
2 yes 0.01659099 0.038448013 Inf -0.05876573 0.09194771- Is the slope of
GIVENNESSwhenCONTRAST: no significantly different from whenCONTRASTis yes?
Estimate Std. Error z value Pr(>|z|)
3.387314e-01 3.951911e-02 8.571330e+00 1.022962e-17 - We saw the confidence intervals for the two slopes of
GIVENNESSabove already but now, what is the p-value for whether the slope ofGIVENNESSis significantly diff from 0- when
CONTRASTis no? - when
CONTRASTis yes?
- when
# 5. Here's EVERYTHING:
glht(m_final, matrix(c(0,0,1,0,0), nrow=1)) %>% { list(
"Summary"=summary(.)$test[3:6] %>% unlist,
"ConfInt"=confint(.)$confint[2:3]) }$Summary
coefficients.1 sigma.1 tstat.1 pvalues
-0.322140407 0.009251289 -34.821137322 0.000000000
$ConfInt
[1] -0.3402726 -0.3040082glht(m_final, matrix(c(0,0,1,0,1), nrow=1)) %>% { list(
"Summary"=summary(.)$test[3:6] %>% unlist,
"ConfInt"=confint(.)$confint[2:3]) }$Summary
coefficients.1 sigma.1 tstat.1 pvalues
0.01659099 0.03844801 0.43151745 0.66609216
$ConfInt
[1] -0.05876573 0.091947711.5.2 Making & understanding predictions manually
Let’s take a random speaker, e.g. 25-JE_NNS, and compute the predicted probability of this speaker realizing the subject when GIVENNESS is 4 and `CONTRAST is yes? Ideally, you’d do this manually first and then check it ‘less manually’.
(Intercept) CONTRASTyes GIVENNESS
1.1337588 0.3959162 -0.3221404
SPKRnns CONTRASTyes:GIVENNESS
0.2229958 0.3387314 (Intercept) CONTRASTyes
25-JE_NNS -0.3617135 1.209507 (Intercept) CONTRASTyes GIVENNESS SPKRnns CONTRASTyes:GIVENNESS
25-JE_NNS 0.7720454 1.605424 -0.3221404 0.2229958 0.3387314# completely manually
prediction_lo <-
1.1337588 + # overall intercept
-0.3617135 + # adjustment for the SPEAKERID
0.3959162 + # CONTRAST: yes
1.209507 + # adjustment for the SPEAKERID
(-0.3221404 * 4) + # GIVENNESS: 4
0.2229958 + # SPKR: nns
( 0.3387314 * 4) # CONTRAST:GIVENNESS
prediction_pp <- plogis(prediction_lo)
c("predicted log odds" =prediction_lo,
"predicted probability"=prediction_pp) predicted log odds predicted probability
2.6668283 0.9350406 # with predict
qwe <- data.frame(CONTRAST="yes", GIVENNESS=4, SPKR="nns", SPEAKERID="25-JE_NNS")
c("predicted log odds" =predict(m_final, newdata=qwe),
"predicted probability 1"=plogis(predict(m_final, newdata=qwe)),
"predicted probability 1"=predict(m_final, newdata=qwe, type="response")) predicted log odds.1 predicted probability 1.1 predicted probability 1.1
2.6668287 0.9350407 0.9350407 SPKR CONTRAST GIVENNESS PRED_LO_2
38 nns yes 4 1.819035[1] 1.819035(qwe <- 1.819035 +
# adjust it by the relevant random effects, ...
-0.3617135 + # adjustment to the intercept for the speaker
1.209507) # adjustment to CONTRAST for the speaker[1] 2.666829[1] 0.93504072 Session info
R version 4.5.2 (2025-10-31)
Platform: x86_64-pc-linux-gnu
Running under: Pop!_OS 22.04 LTS
Matrix products: default
BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so; LAPACK version 3.10.0
locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
[3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
[5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
[7] LC_PAPER=en_US.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
time zone: America/Los_Angeles
tzcode source: system (glibc)
attached base packages:
[1] stats graphics grDevices utils datasets compiler methods
[8] base
other attached packages:
[1] MuMIn_1.48.11 multcomp_1.4-29 TH.data_1.1-5 MASS_7.3-65
[5] survival_3.8-6 mvtnorm_1.3-3 lme4_1.1-38 Matrix_1.7-4
[9] emmeans_2.0.1 effects_4.2-4 car_3.1-5 carData_3.0-6
[13] dplyr_1.2.0 STGmisc_1.06 Rcpp_1.1.1 magrittr_2.0.4
loaded via a namespace (and not attached):
[1] dotCall64_1.2 spam_2.11-3
[3] xfun_0.56 htmlwidgets_1.6.4
[5] insight_1.4.6 lattice_0.22-9
[7] vctrs_0.7.1 tools_4.5.2
[9] Rdpack_2.6.6 generics_0.1.4
[11] parallel_4.5.2 stats4_4.5.2
[13] sandwich_3.1-1 tibble_3.3.1
[15] pkgconfig_2.0.3 RColorBrewer_1.1-3
[17] LCFdata_2.0 lifecycle_1.0.5
[19] fields_17.1 LMERConvenienceFunctions_3.0
[21] mitools_2.4 codetools_0.2-20
[23] survey_4.4-8 htmltools_0.5.9
[25] maps_3.4.3 yaml_2.3.12
[27] Formula_1.2-5 pillar_1.11.1
[29] nloptr_2.2.1 reformulas_0.4.4
[31] boot_1.3-32 abind_1.4-8
[33] nlme_3.1-168 tidyselect_1.2.1
[35] digest_0.6.39 splines_4.5.2
[37] fastmap_1.2.0 grid_4.5.2
[39] colorspace_2.1-2 cli_3.6.5
[41] estimability_1.5.1 rmarkdown_2.30
[43] otel_0.2.0 nnet_7.3-20
[45] zoo_1.8-15 coda_0.19-4.1
[47] evaluate_1.0.5 knitr_1.51
[49] rbibutils_2.4.1 viridisLite_0.4.3
[51] mgcv_1.9-4 rlang_1.1.7
[53] xtable_1.8-4 glue_1.8.0
[55] DBI_1.2.3 rstudioapi_0.18.0
[57] minqa_1.2.8 jsonlite_2.0.0
[59] R6_2.6.1