If not already done, carry out parts 1-5 of the setup, as described here.

Load these libraries and data-frames:

library(tidyverse)
library(emmeans)
library(lmerTest)
library(MuMIn)

urla = "https://www.phonetik.uni-muenchen.de/"
urlb = "studium_lehre/lehrmaterialien/R_speech_processing/Rdf/"
url = paste0(urla, urlb)
soa = read.table(file.path(url, "soa.txt"), 
                 stringsAsFactors = T)

vint.df = read.table(file.path(url, "vint.df.txt"), 
                 stringsAsFactors = T, header = T) %>% 
  mutate(subject = factor(subject, levels = paste0("S", 1:10))) 

form.df = read.table(file.path(url, "form.df.txt"), 
                  stringsAsFactors = T)
dg.df = read.table(file.path(url, "dg.txt"), 
                   stringsAsFactors = T)
vh.df = read.table(file.path(url, "vh.df.txt"), 
                   stringsAsFactors = T)

1 Mixed models and the problem of repetition

All of the statistical tests considered so far – t-tests, ANOVA, linear regression, logistic regression – require the observations to be independent of each other. As a consequence, in all of these tests there can be no more than 1 value for each level of a factor. This was so for the paired sample t-test carried out here in testing whether F2 was influenced by Stress:

form.df %>%
  select(Stress, Subject) %>%
  table()

##             Subject
## Stress       S1 S10 S11 S12 S2 S3 S4 S5 S6 S7 S8 S9
##   stressed    1   1   1   1  1  1  1  1  1  1  1  1
##   unstressed  1   1   1   1  1  1  1  1  1  1  1  1

and in which the levels of Stress had just one value for each of stressed and unstressed. If any of the above entries had been greater than one, then the t-test would have been invalid. Exactly the same consideration applies to analysis of variance for all the reasons explained in the orange box here. Both linear and logistic regression are also only valid if there is no more than one pair of data points (one x y pair in y ~ x) per subject. In mixed models, there is no such requirement. Mixed models can also be run if there are missing data as described here.

Another example of missing data in the above t-test case is if e.g. subject S1 had a value only for unstressed:

form.df %>%
  slice(-1) %>%
  select(Stress, Subject) %>%
  table()

##             Subject
## Stress       S1 S10 S11 S12 S2 S3 S4 S5 S6 S7 S8 S9
##   stressed    0   1   1   1  1  1  1  1  1  1  1  1
##   unstressed  1   1   1   1  1  1  1  1  1  1  1  1

which prevents a paired t test from applying

form.df %>%
  slice(-1) %>%
  t.test(F2 ~ Stress, ., paired=T)

## Error in complete.cases(x, y): not all arguments have the same length

because the values (in this case stressed - unstressed) are subtracted from each other in a paired t-test (and this is obviously not possible if a subject has a value only for unstressed but not for stressed).

2 Fixed and random factors

The basic schematic structure of a mixed model is y ~ FF + RF where, as always, y is the dependent variable, FF is one or more so-called fixed factors and RF is one or more random factors. A fixed factor is one that is to be tested: it is the x in all of the other types of test carried out so far, for example:

Is F2 influenced by Stress? (Analysed with a paired-sample t test).

summary(form.df)

##        F2              Stress      Subject  
##  Min.   :1728   stressed  :12   S1     : 2  
##  1st Qu.:2047   unstressed:12   S10    : 2  
##  Median :2150                   S11    : 2  
##  Mean   :2176                   S12    : 2  
##  3rd Qu.:2333                   S2     : 2  
##  Max.   :2581                   S3     : 2  
##                                 (Other):12

With a mixed model, fixed factor: Stress.

Is F2 influenced by Region and Gender? (Analysed with ANOVA).

summary(dg.df)

##        F2       Region Gen         Vpn    
##  Min.   :1082   A:20   m:30   S1     : 1  
##  1st Qu.:1518   B:20   w:30   S10    : 1  
##  Median :1780   C:20          S11    : 1  
##  Mean   :1828                 S12    : 1  
##  3rd Qu.:2179                 S13    : 1  
##  Max.   :2682                 S14    : 1  
##                               (Other):54

With a mixed model, fixed factors: Region, Gender.

Is Height influenced by Volume? (Analysed with linear regression here).

summary(trees)

##      Girth           Height       Volume     
##  Min.   : 8.30   Min.   :63   Min.   :10.20  
##  1st Qu.:11.05   1st Qu.:72   1st Qu.:19.40  
##  Median :12.90   Median :76   Median :24.20  
##  Mean   :13.25   Mean   :76   Mean   :30.17  
##  3rd Qu.:15.25   3rd Qu.:80   3rd Qu.:37.30  
##  Max.   :20.60   Max.   :87   Max.   :77.00

With a mixed model, fixed factor: Volume.

Is Year influenced by Vowel? (Analysed with logistic regression).

summary(vh.df)

##       Year       Vowel    
##  Min.   :1950   high: 82  
##  1st Qu.:1960   low :138  
##  Median :1971             
##  Mean   :1976             
##  3rd Qu.:1993             
##  Max.   :2005

With a mixed model, fixed factor: Vowel.

The above examples also show why mixed models are general: in examples 3. and 4., the fixed factor is continuous, whereas in example 4, the dependent variable (y in y ~ x) is categorical. In all of the other tests (t-test, ANOVA, linear regression, logistic regression) there is not the same flexibility. The main requirement for being able to use a mixed model, however, is that there must be at least one random factor, as detailed below.

A random factor (which can only ever be a factor, i.e., categorical) is used for those factors whose variation is noise or irrelevant and should be factored out. In phonetics and speech processing and also in many linguistic studies, subjects (speakers, listeners, participants) and items (typically words) are very often random factors. There is theoretically no limit to the number of random factors in a mixed model (just as there are no limits to the number of fixed factors).

Consider as an example, the problem of deciding whether dB is influenced by vowel in the data frame vint.df:

 vint.df %>%
  ggplot() +
  aes(y = db, x = vowel) +
  geom_boxplot()

The overlap between /a/ and /i/ is quite large. But part of the reason for this overlap is because there is a lot of variation in the db level between subjects, as this plot shows:

 vint.df %>%
  ggplot() +
  aes(y = db, x = subject) +
  geom_boxplot()

The way that this by-subject variation can interfere with the fixed factor is evident in e.g. considering the data just from subjects S7 and S8 in which the /a, i/ overlap is especially high:

vint.df %>%
  filter(subject %in% c("S7", "S8")) %>%
  ggplot() +
  aes(y = db, x = vowel) +
  geom_boxplot() +
  ggtitle("S7 and S8 only")

Declaring subject as a random factor in this case helps to cut out the variation between speakers leaving a cleaner separation between the levels of the fixed factor i.e. between /a, i/.

Words can also cause variation that is irrelevant or detrimental to the influence of a fixed factor on a dependent variable. In the data frame soa:

summary(soa)

##  Vpn        W         Pos          F1       
##  s1:6   Bart :6   final :9   Min.   :477.0  
##  s2:6   Pfad :6   medial:9   1st Qu.:498.8  
##  s3:6   Start:6              Median :532.0  
##                              Mean   :526.6  
##                              3rd Qu.:547.0  
##                              Max.   :579.5

The issue is whether F1 is affected by the position (medial, final) of an /a/ vowel in a phrase (Pos):

soa %>% 
  ggplot() +
  aes(y = F1, x = Pos) +
  geom_boxplot()

There is some overlap between medial and final, as the above plot shows. But once again some of this overlap is due to the variation between words:

soa %>% 
  ggplot() +
  aes(y = F1, x = W, fill = Pos) +
  geom_boxplot()

In particular, the F1 values in Bart are low (as a result of which final/Bart overlaps with medial/Pfad and medial/Start). Once again, this type of variation can to a large extent be factored out in a mixed model by declaring W (the words) as a random factor.

3 Specifying the formula for a mixed model

The command for running a mixed model in R in which the dependent variable, y is continuous is lmer() (if the dependent variable is categorical, it is glmer()). The syntax is more or less the same as for the aov_4() function in running an analysis of variance.

Thus where FF is a fixed factor and RF a random factor, the syntax is y ~ FF + (1|RF) or y ~ FF + (FF|RF). Just as in analysis of variance, the choice between (1|RF) and (FF|RF) depends on whether the fixed factor is within or between with respect to the random factor. But in this case, the concept of between vs. within is extended to words (items) as a random factor.

In an a t-test and analysis of variance, a factor can only be either between or within. In a mixed model by contrast, factors need not be strictly within or between. For example, if subject S1’s /a/ and subject S6’s /i/ vowels are removed:

# as vint.df but without S1-/a/ and S6-/a/
vint.df %>%
  filter(!(subject == "S1" & vowel == "a")) %>%
  filter(!(subject == "S6" & vowel == "i")) %>%
  select(subject, vowel) %>%
  table()

##        vowel
## subject a i
##     S1  0 6
##     S2  6 6
##     S3  6 6
##     S4  6 6
##     S5  6 6
##     S6  6 0
##     S7  6 6
##     S8  6 6
##     S9  6 6
##     S10 6 6

then vowel is no longer strictly within because both S1 and S6 have values for only one level of vowel. Nevertheless, because a within relationship exists for most subjects, (vowel|subject) is still appropriate. From another point of view: a factor can have empty cells and still be analysed with a mixed model, whereas this is not possible in an ANOVA.

For the data frame vint.df, it was established earlier than vowel is a fixed factor while subject and word are random factors. The first part of the lmer() command involving the fixed factor is straightforwardly db ~ vowel. The following shows that vowel is within with respect to subject:

vint.df %>%
  select(subject, vowel) %>%
  table()

##        vowel
## subject a i
##     S1  6 6
##     S2  6 6
##     S3  6 6
##     S4  6 6
##     S5  6 6
##     S6  6 6
##     S7  6 6
##     S8  6 6
##     S9  6 6
##     S10 6 6

because there are values for each subject at both levels a and i. The following:

vint.df %>%
  select(word, vowel) %>%
  table()

##     vowel
## word  a  i
##   w1 20  0
##   w2 20  0
##   w3 20  0
##   w4  0 20
##   w5  0 20
##   w6  0 20

shows that vowel is between with respect to word because each word occurs either in /a/ or in /i/. Thus, the second part of the command to lmer() with the random factors must be (vowel|subject) + (1|word) so that putting the fixed and random parts together, the entire command is db ~ vowel + (vowel|subject) + (1|word).

In mixed models, (FF|RF) means that FF interacts with RF. For the present example, the interaction in (vowel|subject) means that the relationship between i and a might not be the same for each subject. A plot of how vowel varies by subject:

vint.df %>%
  ggplot() +
  aes(y = db, x = subject, fill = vowel) +
  geom_boxplot()

shows that an interaction between these factors is likely. This is because the difference between [a] and [i] is e.g., much greater for S1 than it is for S7 and S10.

For words, the random factor is (1|word) and not (vowel|word). This is because vowel is between and not within with respect to word. For this reason, there can be no interaction between vowel and word. Thus, the interaction question: ‘is the relationship between [a] and [i] the same for each word’ is meaningless, because each word has either [a] or [i] and not both:

vint.df %>%
  ggplot() +
  aes(y = db, x = word, fill = vowel) +
  geom_boxplot()

The lmer() command will now be extended to the case of two factors using the same data frame that also includes a factor sex specifying each subject as male or female. The issue to be analysed is whether the dependent variable, db is affected by vowel and/or sex.

Just as in analysis of variance, the fixed factors include an interaction beween vowel and sex. Therefore, the part of the command to lmer() specifying the fixed factors must be: vowel * sex. The relationship between vowel and the random factors subject and word were already determined above: (vowel|subject) + (1|word).

To determine the relationship between sex and the random factors requires once again establishing whether sex is within or between. The fixed factor sex must be between with respect to subject because subjects are either male or female. For the relationship between sex and word:

vint.df %>%
  select(word, sex) %>%
  table()

##     sex
## word  f  m
##   w1 10 10
##   w2 10 10
##   w3 10 10
##   w4 10 10
##   w5 10 10
##   w6 10 10

sex is within with respect to word. This therefore gives: (1|word) + (sex|word).

Thus the part of the formula with random factors is:

(vowel|subject) + (1|word) + (1|word) + (sex|word)

However, the repetition of (1|word) is unnecessary. Moreover, whenever there is (1|RF) + (FF|RF) then, for reasons to do with intercept and slope calculations explained below, these can just be combined as (i.e., are equivalent to) (FF|RF). The desired model including the fixed factors is therefore:

db ~ vowel * sex + (vowel|subject) + (sex|word)

The other type of contraction in the formula is if there are two fixed factors, FF1 and FF2 and both are within with respect to the same random factor to give (FF1|RF) + (FF2|RF). The terms should be combined as (FF1 + FF2|RF).

4 Interpreting the coefficients

The following section provides some background information on how the coefficients of an lmer() model relate to the mean values of the variables in the data frame. It also gives some information about the relationship between actual and predicted values. (This section can be skipped until after having read the next section on tests for significance).

A mixed model is in some ways similar to linear regression in that both involve calculating intercepts and slopes of straight lines through the data such that the difference between the actual and predicted values is minimised. In linear regression, the intercept k and the slope m were calculated from

$\hat{y} = m x + k$

where $\hat{y}$ are the predicted values, $k$ und $m$ are the intercept and slope respectively.

In a mixed model, an intercept and slope are calculated for each fixed factor; if the formula has (1|RF) then only an intercept is calculated for the random factor; and if the formula has (FF|RF), then both an intercept and slope are calculated for RF.

In the simpler model considered above with one fixed factor i.e. db ~ vowel + (vowel|subject) + (1|word):

a slope and intercept are calculated for vowel
an intercept is calculated for subject and as well as a slope with respect to vowel.
an intercept is calculated for word.

The corresponding formula is:

$\hat{db} = (subject.m + vowel.m).vowel + vowel.k + subject.k + word.k$

where

$\hat{db}$ is the predicted value.
$subject.m$ is the slope for the $subject$.
$vowel.m$ is the slope of the fixed factor $vowel$.
$vowel$ consists of two levels of the factor that are represented as 0 and 1.
$vowel.k$ is the intercept for $vowel$.
$subject.k$ is the intercept for $subject$.
$word.k$ is the intercept for $word$.

If a random factor is specified for an intercept only, then the model takes account of the differences between the levels of the factor – thus, in the case of word specified with an intercept only, of the differences between:

vint.df %>%
  ggplot() +
  aes(y = db, x = word) +
  geom_boxplot()

If a random factor is specified for an intercept and a slope, then the model adjusts not only for the differences between the levels of the factor, as above, but also for the interaction between the fixed and random factors: thus, in the case of (vowel|subject) for the variation between subjects in the /i ~ a/ differences i.e for the differences between the red and blue boxes here:

vint.df %>%
  ggplot() +
  aes(y = db, x = subject, fill = vowel) +
  geom_boxplot()

The lmer() function calculates all these 5 intercept and slope variables in such a way that the difference between the actual values, $db$ and the predicted values $\hat{db}$ are minimised.

The coefficients for the fixed factor are found in the last line of the output here:

i.e.

vint.df %>%
  lmer(db ~ vowel + (vowel|subject) + (1|word), .)

## Linear mixed model fit by REML ['lmerModLmerTest']
## Formula: db ~ vowel + (vowel | subject) + (1 | word)
##    Data: .
## REML criterion at convergence: 849.5068
## Random effects:
##  Groups   Name        Std.Dev. Corr 
##  subject  (Intercept) 35.062        
##           voweli      23.644   -0.90
##  word     (Intercept) 15.787        
##  Residual              5.678        
## Number of obs: 120, groups:  subject, 10; word, 6
## Fixed Effects:
## (Intercept)       voweli  
##       82.27       -34.12

and are also given with the command fixef():

vint.df %>%
  lmer(db ~ vowel + (vowel|subject) + (1|word), .) %>%
  fixef()

## (Intercept)      voweli 
##    82.26771   -34.11696

These fixed effect coefficients are related to the mean of /a/ and /i/ as follows:

vint.df %>%
  group_by(vowel) %>%
  summarise(mean(db))

## # A tibble: 2 × 2
##   vowel `mean(db)`
##   <fct>      <dbl>
## 1 a           82.3
## 2 i           48.2

That is, the intercept is the mean-db of /a/ and the slope is given by the mean-db of /i/ minus the mean-db of /a/ i.e. 48.2 - 82.3 = -34.1. One of the things that the model tests is whether the slope is significantly different from zero. This is also a test of whether the db values for /i, a/ differ from each other: clearly, the closer the slope is to zero, then the less likely there is to be an /i, a/ difference.

In the more complicated model in which vowel and sex interact:

vint.df %>%
  lmer(db ~ vowel * sex + (vowel|subject) + (sex|word), .) %>%
  fixef()

## (Intercept)      voweli        sexm voweli:sexm 
##    52.17194   -14.01031    60.19154   -40.21331

The four coefficients are also related to mean-db values:

vint.df %>%
  group_by(vowel, sex) %>%
  summarise(mean_db = mean(db))

## `summarise()` has grouped output by 'vowel'. You can override using the
## `.groups` argument.

## # A tibble: 4 × 3
## # Groups:   vowel [2]
##   vowel sex   mean_db
##   <fct> <fct>   <dbl>
## 1 a     f        52.2
## 2 a     m       112. 
## 3 i     f        38.2
## 4 i     m        58.1

In this case:

The intercept is the mean-db for female-/a/.
The coefficient voweli is the mean-db for female-/i/ minus the mean-db for female-/a/ i.e. 38.16164 - 52.17194 = -14.0103. This shows that female-/i/ has a lower db than female-/a/.
The coefficient sexm is the mean-db for male-/a/ minus the mean-db for female-/a/ i.e. 112.36348 - 52.17194 = 60.19154. This shows that the mean-db or male-/a/ is greater than for female-/a/.
The coefficient voweli:sexm is the mean-db for male-/i/ minus the mean-db for female-/a/ minus the mean-db for female /i/ i.e. 58.13987 - 52.17194 - (-14.01031) - 60.19154 = -40.2133. (From another point of view: the mean-db for male-/i/ is equal to the sum of all four coefficients).

The coefficients for the random factors are given by ranef(). For the simpler model:

vint.df %>%
  lmer(db ~ vowel + (vowel|subject) + (1|word), .) %>%
  ranef()

## $subject
##     (Intercept)     voweli
## S1    48.375321 -38.067670
## S2    21.366717 -20.940391
## S3    40.151416 -10.712792
## S4     9.200103  -9.644261
## S5    30.562029 -19.915777
## S6   -18.299866   8.171719
## S7   -29.386695  30.077508
## S8   -10.061478  20.081715
## S9   -40.238101   9.656621
## S10  -51.669445  31.293327
## 
## $word
##    (Intercept)
## w1   -4.745086
## w2  -12.493797
## w3   17.238882
## w4   -7.991323
## w5  -10.477787
## w6   18.469110
## 
## with conditional variances for "subject" "word"

The above shows that there is:

an intercept and a slope for each subject (in $subject).
an intercept for each word (in $word).

The following gives an example of how the predicted values are calculated for the first observation:

vint.df %>%
  slice_head(n=1)

##    db vowel sex subject word
## 1 124     a   m      S1   w1

Predicted values are obtained by substituting real values (based on the first observation in this case) into the variables of the formula:

$\hat{db} = (subject.m + vowel.m).vowel + vowel.k + subject.k + word.k$

The substitutions are as follows:

$vowel$. This substitution is given by converting the levels of the factor (/a/, /i/) into a matrix of 0 and 1 (for two levels as in this case), as follows:

vint.df %>%
  pull(vowel) %>%
  contrasts()

##   i
## a 0
## i 1

which gives zero for /a/ and 1 for /i/. Since the first observation is of the vowel /a/, then $vowel$ = 0 (zero), which means that the above formula reduces to:

$\hat{db} = vowel.k + subject.k + word.k$

$vowel.k$: this is 82.26771 (see above under fixed factor).
$subject.k$: as shown by vint.df %>% slice_head(n=1), the subject is S1. The random intercept for S1 is 48.375266 (see above under random factor).
$word.k$: as shown by vint.df %>% slice_head(n=1), the word is w1. The random intercept for w1 is -4.745100.

The predicted value for the first observation is therefore the sum of these terms i.e. 0 + 82.26771 + 48.375266 -4.745100 = 125.8979.

The predicted values in an lmer() model are automatically given by the function fitted(). The first predicted value is:

predicted_values = vint.df %>%
  lmer(db ~ vowel + (vowel|subject) + (1|word), .) %>%
  fitted()
predicted_values[1]

##        1 
## 125.8979

which is the same as the predicted value that has just been calculated manually. The fit is quite good (given the actual value of 124, as shown by vint.df %>% slice_head(n=1) above.

The function r.squaredGLMM() gives some information about the proportion of variance in the original data explained by the model. The closer the total score is to 100%, the better the model’s fit to the original data. The percentage explained is broken down into the proportion explained by the fixed factor(s) R2m and the total explained variance, R2c:

vint.df %>%
  lmer(db ~ vowel + (vowel|subject) + (1|word), .) %>%
  r.squaredGLMM()

## Warning: 'r.squaredGLMM' now calculates a revised statistic. See the help page.

##            R2m       R2c
## [1,] 0.2191709 0.9759183

In this case, the fixed factor vowel explains 21.9% of the variance and the fixed and random factors together 97.6% (so the random factors explain around 97.6 - 21.9 = 75.7% of the variance). The percentage for the random factors is high because sex hasn’t been included as a fixed factor (so all the variation due to male-female differences is absorbed into the random factors). If sex is included as a fixed factor:

vint.df %>%
  # the formula established earlier
  lmer(db ~ vowel * sex + (vowel|subject) + (sex|word), .) %>%
  r.squaredGLMM()

##            R2m       R2c
## [1,] 0.6088086 0.9776962

then, as expected, the percentage of variance explained by the fixed factors increases to 60.88%

5 Tests of significance

There are three parts to testing whether the fixed factor(s) have a significant influence on the dependent variable:

test whether the model can be simplified with the step() function.
test for significance with anova().
if need be, run post-hoc tests using emmeans().

5.1 The `step()` function

This useful function tests whether all of the terms that were entered into the model really are needed. It then removes those terms that have a negligible effect on minimising the difference between the actual and predicted values. The following applies the step() function to vint.df with two fixed factors and two random factors (in which the aim is to assess whether vowel and/or sex influence db):

vint.df %>%
  # the formula established earlier
  lmer(db ~ vowel * sex + (vowel|subject) + (sex|word), .) %>%
  step()

## Backward reduced random-effect table:
## 
##                            Eliminated npar  logLik     AIC     LRT Df
## <none>                                  11 -408.18  838.35           
## sex in (sex | word)                 1    9 -410.33  838.66   4.304  2
## vowel in (vowel | subject)          0    7 -433.46  880.92  46.260  2
## (1 | word)                          0    8 -498.83 1013.66 176.999  1
##                            Pr(>Chisq)    
## <none>                                   
## sex in (sex | word)            0.1163    
## vowel in (vowel | subject)   9.01e-11 ***
## (1 | word)                  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##           Eliminated Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## vowel:sex          0 960.84  960.84     1 8.0006  29.801 0.0006021 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## db ~ vowel + sex + (vowel | subject) + (1 | word) + vowel:sex

The most important part is the final line which shows whether or not and how the model may have been simplified. In this case, (sex|word) in the original model has been replaced with (1|word). This means that there is no significant interaction between sex and word (given that, as explained earlier, (sex|word) means that the model is being specified for an interaction between the fixed factor sex and the random factor word). This is the reason for the non-significant p-value = 0.1163 in the line under ‘Backward reduced random-effect table:’

sex in (sex | word) 1 9 -410.33 838.66 4.304 2 0.1163

All of the other terms are significant – or rather if they were removed, then the difference between the actual and predicted values would be significantly greater than if they were included.

Where A and B are two factors, then A+B+A:B is the same as A*B (and also as (A+B)^2). Therefore,the formula returned by step():

db ~ vowel + sex + (vowel | subject) + (1 | word) + vowel:sex

is the same as:

db ~ vowel * sex + (vowel | subject) + (1 | word)

This also shows that the only difference compared with the original formula entered into lmer() is (1|word) instead of (sex|word). (It also shows that vowel and subject interact with each other significantly: since if they did not, the step() would have replaced (vowel | subject) with (1 | subject)).

It can sometimes happen (and typically does with more complicated lmer() models with many factors) that the model gives an error message about a singularity or failing to converge. This usually means that one of the variables cannot be estimated typically because there is not enough data to support it (or more generally because there are too many variables in the model).

As an example, the following includes the variable (sex|subject) which for all the reasons stated earlier is wrong (because sex is between with respect to subject):

vint.df %>%
  # the formula established earlier
  lmer(db ~ vowel * sex + (sex|subject) + (sex|word), .)

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## unable to evaluate scaled gradient

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge: degenerate Hessian with 1 negative eigenvalues

## Warning: Model failed to converge with 1 negative eigenvalue: -1.9e-04

## Linear mixed model fit by REML ['lmerModLmerTest']
## Formula: db ~ vowel * sex + (sex | subject) + (sex | word)
##    Data: .
## REML criterion at convergence: 866.0716
## Random effects:
##  Groups   Name        Std.Dev. Corr 
##  subject  (Intercept) 15.494        
##           sexm         9.037   -0.59
##  word     (Intercept) 15.491        
##           sexm         2.821   0.14 
##  Residual              7.716        
## Number of obs: 120, groups:  subject, 10; word, 6
## Fixed Effects:
## (Intercept)       voweli         sexm  voweli:sexm  
##       52.17       -14.01        60.19       -40.21  
## optimizer (nloptwrap) convergence code: 0 (OK) ; 0 optimizer warnings; 2 lme4 warnings

Usually it is best to identify what can be simplified (or where the error is) and then try again. Alternatively, the step() function might find the error by eliminating the bad term, thus:

vint.df %>%
  # the formula established earlier
  lmer(db ~ vowel * sex + (sex|subject) + (sex|word), .) %>%
  step()

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## unable to evaluate scaled gradient

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge: degenerate Hessian with 1 negative eigenvalues

## Warning: Model failed to converge with 1 negative eigenvalue: -1.9e-04

## boundary (singular) fit: see help('isSingular')

## Backward reduced random-effect table:
## 
##                        Eliminated npar  logLik     AIC     LRT Df Pr(>Chisq)
## <none>                              11 -433.04  888.07                      
## sex in (sex | subject)          1    9 -433.12  884.25   0.175  2     0.9164
## sex in (sex | word)             2    7 -433.46  880.92   0.672  2     0.7147
## (1 | subject)                   0    6 -493.16  998.32 119.400  1     <2e-16
## (1 | word)                      0    6 -499.58 1011.16 132.237  1     <2e-16
##                           
## <none>                    
## sex in (sex | subject)    
## sex in (sex | word)       
## (1 | subject)          ***
## (1 | word)             ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##           Eliminated Sum Sq Mean Sq NumDF DenDF F value    Pr(>F)    
## vowel:sex          0  12128   12128     1   104   198.6 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## db ~ vowel + sex + (1 | subject) + (1 | word) + vowel:sex

so that the suggested formula db ~ vowel + sex + (1 | subject) + (1 | word) + vowel:sex is now error free.

5.2 The `anova()` function

This tests whether the fixed factors have a significant influence on the dependent variable. This should be run on the suggested formula that is output from the step() function above, which can be most easily done by copying the formula that is suggested and then passing the output of lmer() to anova(), thus:

vint.df %>%
  lmer(db ~ vowel + sex + (vowel | subject) 
       + (1 | word) + vowel:sex, .) %>%
  anova()

## Type III Analysis of Variance Table with Satterthwaite's method
##           Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## vowel     208.82  208.82     1 4.6051  6.4767 0.0556781 .  
## sex       637.75  637.75     1 7.9993 19.7802 0.0021469 ** 
## vowel:sex 960.84  960.84     1 8.0006 29.8007 0.0006021 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The output shows that there was no significant influence of vowel on db (p = 0.56), but a significant influence of sex ($F[1, 8.00] = 18.78, p < 0.01$) and a significant interaction between these fixed factors ($F[1, 8.00] = 28.8, p < 0.001$).

5.3 `emmeans()` and `pairs()`

These functions for carrying out post-hoc tests are the same as those introduced in a previous module on analysis of variance. The syntax and usage are exactly the same as before, thus:

vint.df %>%
  lmer(db ~ vowel + sex + (vowel | subject) 
       + (1 | word) + vowel:sex, .) %>%
  emmeans(., ~ vowel * sex) %>%
  pairs(., simple = "each", combine=T)

##  sex vowel contrast estimate    SE   df t.ratio p.value
##  f   .     a - i        14.0 13.90 5.27   1.008  1.0000
##  m   .     a - i        54.2 13.90 5.27   3.900  0.0411
##  .   a     f - m       -60.2 10.14 8.00  -5.938  0.0014
##  .   i     f - m       -20.0  9.32 8.00  -2.144  0.2576
## 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: bonferroni method for 4 tests

The results of post-hoc tests show a significant difference between /a/ and /i/ for males (p < 0.05) but not for females. They also show a significant difference between males and females in /a/ (p < 0.01) but not in /i/. These results are fully consistent with the figure presented earlier showing boxplots for these two factors:

vint.df %>%
  ggplot() +
  aes(y = db, x = vowel, fill = sex) +
  geom_boxplot()

An introduction to mixed models

Jonathan Harrington, Johanna Cronenberg

1 Mixed models and the problem of repetition

2 Fixed and random factors

3 Specifying the formula for a mixed model

4 Interpreting the coefficients

5 Tests of significance

5.1 The `step()` function

5.2 The `anova()` function

5.3 `emmeans()` and `pairs()`

An introduction to mixed models

Jonathan Harrington, Johanna Cronenberg

1 Mixed models and the problem of repetition

2 Fixed and random factors

3 Specifying the formula for a mixed model

4 Interpreting the coefficients

5 Tests of significance

5.1 The step() function

5.2 The anova() function

5.3 emmeans() and pairs()

5.1 The `step()` function

5.2 The `anova()` function

5.3 `emmeans()` and `pairs()`