A Linear Mixed Model Tutorial, based on a materials by Jonathan Harrington, IPS, Munich

Setup

The code for this tutorial can be found at: http://rpubs.com/hywel/lmem_workshop

Before starting the tutorial load the required packages and the data sets we will use for this seminar:

We will be using R version 3.3.0 and Rstudio version 0.99.902. Any recent versions of these applications will do and the syntax should run on all platforms.

###########################################################
#### Load the libraries ###
# install.packages(c("lmerTest","lattice","ez"))
# Use the above command if these packages are not installed.
library(lmerTest)
library(lattice)
library(ez)

#### Load the functions and the datasets ###
pfadu = "http://www.phonetik.uni-muenchen.de/~jmh/lehre/Rdf"
source(file.path(pfadu, "phoc.txt"))
source(file.path(pfadu, "ph.step.R"))

#### Datasets ####
stimm = read.table(file.path(pfadu, "stimm.df.txt"))
stimm2 = read.table(file.path(pfadu, "stimm2.df.txt"))
soa = read.table(file.path(pfadu, "soa.txt"))
int = read.table(file.path(pfadu, "dbwort.df.txt"))
###########################################################

---

Introduction

This seminar mainly concerns model building for Linear Mixed Effects models (LMEM) and how to construct a dataset that is suitable for use in modern statisistical analyses. These materials have been translated from the German by Hywel (and Google Translate). Most of the variable names remain in German and the expression in some of the harder to translate sections may be a little clumsy. The original can be found here: http://www.phonetik.uni-muenchen.de/~jmh/lehre/sem/ss16/statistik.htm: Gemischte Modelle (mixed models) mmbefehle.R. Suggestions of clearer, more semantically coherent translations are welcomed.

General Linear Mixed Models ¹

In a mixed model (LMM) a dependent variable (continuous or categorical) is examined to see if it is influenced by one or more independent factors. Mixed models (LMM) have a broad application. They partially replace the analysis of variance; They enhance linear and logistic regression analyses.

---

Some benefits (A-D) that come from a LMM vs. analysis of variance (ANOVA).

A. An LMM is similar to an ANOVA

ANOVAs can be used a mixed design and have random effects included.

But unlike an ANOVA a LMM does not have to be averaged over repetitions. (Don’t worry about the ezANOVA syntax in this example, this is just to show that it is possible to run one on these data).

For example:

head(stimm)

Note: The variables in stimm are:

• Voice Onset Time (vot)
• Speaker (Vpn) verschiedene Sprecher or “different speakers”?
• Place of Articulation (K) a factor with 2 levels “ba” and “pa”.
• Age (Alter) with 2 levels adults (Erw) and children (Kind)

To what extent is the VOT (vot) influenced by place of articulation and age?

1. First step:

stimm.m = aggregate(vot ~ K * Vpn * Alter, mean, data = stimm)
ezANOVA(stimm.m, .(vot), .(Vpn), .(K))

So the ANOVA shows that place of articulation (K) does have a significant effect on VOT (p < 0.05).

B. A Balanced-design is not required in a LMM.

Balanced design: In a balanced design, the number of subjects must be equal to each stage of the between-factor.. That is clearly not given for the between-age factor.

For these data:

with(stimm, table(Vpn, Alter))

In the table above we can see that there are 3 adults and 5 children. Therefore we cannot use an ANOVA as these data are not balanced and the cde below will produce an error:

ezANOVA(stimm.m, .(vot), .(Vpn), .(K), between = .(Alter))

This form of data is no problem for an LMM.

C. Missing Values Are Allowed

In an ANOVA each stage of the within-group-factor must be occupied. This is not true for a Linear Mixed Model.

with(stimm2, table(Vpn, K))

We see both speaker (Vpn) A and speaker D have missing values (“ba” for A and “pa” for D). Therefore an ANOVA can not be performed.

stimm2.m = aggregate(vot ~ K * Vpn * Alter, mean, data = stimm2)
ezANOVA(stimm2.m, .(vot), .(Vpn), .(K), between = .(Alter))

Here too a LMM will have no problems.

D. Exclude the variability within a factor.

In an ANOVA, the variability of factors that are part of a greater population (for example, subjects. OR word) are excluded. There is no such restriction in a LMM (see 2 below).

E. A disadvantage of mixed models:

A LMM is not stable at a small number of observations: less than about 60 the LMM often does not converge (although this also depends on the number of factors).

What are Fixed Factors and Random Factors?

This concept is central to the model design and this is what makes this a ‘mixed’ model.

• Fixed Factors: One or more factors that are to be tested
• Random Factors: One or more factors, for which the variability is to be excluded or reduced.

In phonetics are often subject and / or Item (materials such as Words) Random Factors.

An Example:

20 speaker produced several Words, in phrase medial or phrase final position and vowel duration was measured. To what extent the vowel duration affected by the phrase position?

Dataset: soa

• Dependent variable: time
• Fixed Factor: Phrase item (medial vs. final).
• Random Factors: subjects, word (we want to remove the variability that arises from this variable).

What is calculated in an LMM?

A LMM is similar to the linear regression in some respects. We won’t go through this in more detail but a very good introduction by Bodo Winter can be found here http://arxiv.org/pdf/1308.5499.pdf

More about random factors and the choice between (FF | RF) and (1 | RF)

There are two ways to declare a random factor: (FF | RF) or (1 | RF), where FF is a fixed factor and RF is a random factor

head(soa)

Three speakers produced “Bad”, “Pfad”, “Start” either in phrases medial or final position. To what extent does F1 affects the vowel of the phrase item?

• Dependent variable: F1
• Fixed Factor (FF): Pos (medial/final)
• Random Factors (RF): Vpn (Speaker), W (Word)

Therefore there are 4 options when constructing a model:

# option 1.
soa.lmer.1 = lmer(F1 ~ Pos + (Pos| Vpn) + (Pos | W), data = soa)
# option 2.
soa.lmer.2 = lmer(F1 ~ Pos + (1| Vpn) + (Pos | W), data = soa)
# option 3.
soa.lmer.3 = lmer(F1 ~ Pos + (Pos| Vpn) + (1 | W), data = soa)
# option 4.
soa.lmer.4 = lmer(F1 ~ Pos + (1| Vpn) + (1 | W), data = soa)

Basically we will use (FF|RF) (as shown in option 1.), unless the FF is clearly ‘between’ in relation to the RF. For the example above, FF (Pos) is ‘within’ with respect to the RF, speaker (Vpn):

# FF (Pos) is 'within' with respect to the RF, speaker (Vpn):
with(soa, table(Vpn, Pos))

FF is also ‘within’ with regard to word (Wort)

# FF is also 'within' with regard to word (Wort)
with(soa, table(W, Pos))

The importance of (FF | RF)

When constructing the models if we include (FF | RF)

1. A slope RF.m and Intercept RF.k will be calculated
2. (Important): the variability of the RF is excluded, assuming that RF and FF could interact.

i.e. with (Pos | Vpn), examines the extent Position and Speaker (Vpn) interact with each other:

bwplot(F1 ~ Pos|Vpn, data = soa)

Also with (Pos | W), examines the extent Pos and W (word) interact.

bwplot(F1 ~ Pos|W, data = soa)

The importance of (1 | RF)

1. An Intercept RF.k (but no slope) is calculated.
2. (Important): the variability of the RF is excluded without the interaction between RF and FF.

• (1 | RF) to be used when the FF ‘between’ is related to RF. e.g.

head(stimm)

To what extent is vot affected by age? * Fixed factor: age (Alter) * Random factor: Speaker (Vpn) * Age is between the reference to Vpn (speaker):

with(stimm, table(Vpn, Alter))

Age and speaker can not interact because each speaker is either a child or adult:

bwplot(vot ~ Alter|Vpn, data = stimm)

Therefore (Alter | Vpn) would be meaningless - since this requires that age (Alter) and Vpn (speaker) can interact, which is not possible when FF is ‘between’ with regard to the RF.

Therefore:

stimm3 = lmer(vot ~ Alter + (1|Vpn), data = stimm)

and not (Age | Vpn) is used in the final model.

Generally: always use (FF | RF) except when the FF is ‘between’ in relation to RF.

---

The Test Statistics

This is probably the most important part of the tutorial snd the part that needs the most explanation.

• The step () function
• 1. simplifying the model if possible
• 2. examine whether the factors are significant in the simplified model.

The dataframe stimm measures the VOT for a number of syllables (/ ba, pa /) produced by both children and adults (encoded in an age variable (Alter)). To what extent is VOT affected by the place of articulation (K factor) and / or age of speaker (whether the speaker is a child or an adult: Alter factor) ?

Part A: Simplifying the model

First Plot the Interactions

bwplot(vot ~ K | Alter, data = stimm)

# Linear Mixed Model
stimm.lmer = lmer(vot ~ K * Alter + (K|Vpn), data = stimm)
stimm.step = step(stimm.lmer)
stimm.step
# ph.step(stimm.step)
# Part (a): Model simplification
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# Do we really need K + Alter + K: Alter + (K | Vpn)?

Do we really need K + Alter + K: Alter + (K | Vpn)?
Note: this is equivalent notation to K * Alter + (K | Vpn).

If we look at the significance of both the fixed effects and the random effects we see that some of the interactions are not contributing to the model.

Random effects:

interaction	Chi.sq	Chi.DF	elim.num	p.value
K:Vpn	5.04	2	kept	0.0803

Fixed effects:

x	Sum Sq	Mean Sq	NumDF	DenDF	F.value	elim.num
K:Alter	1.5466	1.5466	1	6	0.3511	1
Alter	3.3045	3.3045	1	6	0.7501	2
K	180.8614	180.8614	1	7	41.0563	kept

1. Random effects is almost significant (p = 0.08). Hence (K | Vpn) is retained. Therefore kept i.e. (K | Vpn) was not, so (1 | Vpn) replaced.

2. In the fixed-effects. K:Alter and Alter are not significant, and were removed from the model. But the factor K (place of articulation) remains (kept).

The formula used is:

• formula = vot ~ K + (K | Vpn), data = stimm,

The Simplified Model:

stimm.lmer <- lmer(vot ~  K  + (K | Vpn), data = stimm)
stimm.step <- step(stimm.lmer)
stimm.step

The Final Model:

lme4::lmer(formula = vot ~ K + (K | Vpn), data = stimm)

Part B: The test statistic (post hoc tests)

for the simplified model for the above data:

X	Sum Sq	Mean Sq	NumDF	DenDF	F.value	elim.num	Pr(>F)
K:Alter	1.5466	1.5466	1	6	0.3511	1	0.5751
Alter	3.3045	3.3045	1	6	0.7501	2	0.4197
K	180.8614	180.8614	1	7	41.0563	kept	0.0004

# as you can see here:
stimm2.lmer = lmer(vot ~  K + (K | Vpn), data = stimm)
anova(stimm2.lmer)

Which means that: VOT is significant for a main effect of place of articulation (F [1, 7] = 41.1, p <0.001) but is not influenced by age. Futhermore there is no interaction between these factors.

Post-hoc Tests

Post hoc test are needed if there is an interaction recorded between the fixed factors. These are carried out in a very similar manner to ANOVAs

Here we have another dataset

head(int)

10 speakers (Vpn) produced 6 different words (Wort) with either an / i / or / a / vowel (factor V). The intensity was (dependent variable db) was measured per vowel. The speakers were either male (m) or female (f) in a factor G.

When constructing the model it was found that V and G

#Therefore :
int.lmer = lmer(db ~ V * G + (V+G|Wort) + (V|Vpn), data = int)
int.step = step(int.lmer)
int.step

Based on the above table the final model will be:

# Final model:
#  lme4::lmer(formula = db ~ V + G + (V | Vpn) + (1 | Wort) + V:G,

There is a significant interaction:

# There is a significant interaction:
# Fixed effects:
#  Sum Sq  Mean Sq NumDF DenDF F.value elim.num Pr(>F)
# V   126.6597 126.6597     1     8 29.8298     kept 0.0006
# G    47.3133  47.3133     1     8 11.1428     kept 0.0103
# V:G  31.9690  31.9690     1     8  7.5291     kept 0.0253

Therefore, we must apply a post-hoc test.

#
(int.ph=ph.step(int.step, "bonferroni"))
#

Post-hoc Bonferroni corrected t-tests showed significant differences between / i / and / a / in men (p <0.01) but not in women; and significant differences between men and women in / a / (p <0.05) but not in / i /.

---

1. Note: as mentioned above all variables in this tutorial are in German↩