# ggf.:
#install.packages("lmerTest")
#install.packages("emmeans")

library(lmerTest)
library(emmeans)
library(ggplot2)
library(dplyr)
library(ez)

# post-hoc-Funktionen laden:
source(file.path(pfadu, "phoc.txt"))

# Daten laden
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"))

1. Mixed Models: viele Vor-, wenige Nachteile

Mit einem Mixed Model (MM) (der deutschsprachige Begriff lineare gemischte Modelle wird sehr selten benutzt) wird geprüft, ob eine abhängige Variable (die kontinuierlich (lmer()) oder (wenn glmer() benutzt wird) kategorial sein kann) von einem oder mehreren unabhängigen Faktoren beeinflusst wird. Die unabhängigen Faktoren sind meistens kategorial, können aber auch numerisch sein. “Mixed” oder “gemischt” wird ein Mixed Model dadurch, dass es sowohl Fixed als auch Random Factors geben kann, also sowohl Faktoren, deren Einfluss auf die abhängige Variable den Forschenden interessieren (= Fixed Factor), als auch jene, von denen zu erwarten ist, dass sie zwar höchstwahrscheinlich einen Einfluss auf die abhängige Variable ausüben werden, dieser Einfluss aber für den Forschenden uninteressant ist (und möglicherweise das Zu-Tage-Treten den Effekt des interessierenden Faktors verdeckt).

Einen Random Effect gibt es zwar bei ezANOVA() auch (meistens in unserem Fall die unterschiedlichen Sprecher), aber es kann dort nicht mehr als einen Random Effect geben.

Mixed models (MM) haben eine breite Anwendung. Sie ersetzen:

Vier Vorteile (A-D) von einem MM im Vergleich zur Varianzanalyse (und leider auch ein Nachteil (E))

A. In einer ANOVA, aber nicht in einem Mixed Model muss über Wiederholungen gemittelt werden:

Z.B.: Inwiefern wird “vot” von der Artikulationsstelle (“K”) und vom “Alter” beeinflusst?

head(stimm)
##   vot  K Vpn Alter
## 1  10 ba   A  Kind
## 2 -16 ba   B  Kind
## 3   9 ba   C  Kind
## 4 -21 ba   E   Erw
## 5  10 ba   F   Erw
## 6  -1 ba   G  Kind
# zuerst mitteln: ein Wert pro K-Vpn-Alter-Kombination
stimm.m = stimm %>%
  group_by(K, Vpn, Alter) %>%
  summarise(vot = mean(vot))

stimm.ez = ezANOVA(
  data = stimm.m,
  dv = .(vot),
  wid = .(Vpn),
  within = .(K)
)

stimm.ez
## $ANOVA
##   Effect DFn DFd        F            p p<.05        ges
## 2      K   1   7 41.05629 0.0003646654     * 0.02523055

Hier gibt es also keine Probleme. Aber wenn wir den Effekt des Alters ermitteln möchten, stellen wir fest, dass wir kein ausbalanciertes Design aufzuweisen haben:

B. Kein Balanced-design in einem MM erforderlich:

# Balanced design: Die Anzahl der Vpn. muss pro Stufe des Between-Faktors
# gleich sein. Dies ist hier offensichtlich nicht gegeben für den
# Between-Faktor Alter.
with(stimm, table(Vpn, Alter))
##    Alter
## Vpn Erw Kind
##   A   0   16
##   B   0   16
##   C   0   16
##   D  16    0
##   E  16    0
##   F  16    0
##   G   0   16
##   H   0   16

D.h.: es gibt 3 Erwachsene und 5 Kinder.

Daher meckert die ANOVA:

stimm.ez = ezANOVA(
  data = stimm.m,
  dv = .(vot),
  wid = .(Vpn),
  within = .(K),
  between = .(Alter)
)
## Warning: Data is unbalanced (unequal N per group). Make sure you specified
## a well-considered value for the type argument to ezANOVA().

Das ist zwar “nur” eine Warnmeldung, aber wir können nie sicher sein, ob es wirklich “richtig” ist, was wir hier tun; “unequal N per group” ist aber kein Problem für MMs.

C. Fehlende Werte sind in einem MM erlaubt:

In einer ANOVA muss jede Stufe des within-Faktors belegt sein. Das ist bei einem MM nicht der Fall:

with(stimm2, table(Vpn, K))
##    K
## Vpn ba pa
##   A  0  8
##   B  8  8
##   C  8  8
##   D  8  0
##   E  8  8
##   F  8  8
##   G  8  8
##   H  8  8
# Vpn ba pa
# A  0  8
# B  8  8
# C  8  8
# D  8  0
# E  8  8
# F  8  8
# G  8  8
# H  8  8
# daher kann ANOVA nicht durchgefuehrt werden.
# Mitteln:
stimm2.m = stimm2 %>%
  group_by(K, Vpn, Alter) %>%
  summarise(vot = mean(vot))

ezANOVA(
  data = stimm2.m,
  dv = .(vot),
  wid = .(Vpn),
  within = .(K),
  between = .(Alter)
)
## Error in ezANOVA_main(data = data, dv = dv, wid = wid, within = within, : One or more cells is missing data. Try using ezDesign() to check your data.

Fehlende Daten sind aber kein Problem für ein MM.

D. Die Variabilität von mehr als nur einem einzigen Faktor ausklammern:

In einer ANOVA kann die Variabilität von höchstens einem Faktor (z.B. Vpn oder Wort) ausgeklammert werden. Solche Beschränkung gibt es für ein MM nicht (siehe 2. unten).

E. Ein Nachteil: Ein MM ist bei einer kleinen Anzahl von Beobachtungen nicht stabil

Bei weniger als ca. 60 Beobachtungen (bei jeweils einem Fixed und einem Random Factor) konvergiert das MM oft nicht (obwohl das auch von der Anzahl der Faktoren abhängt - je mehr Faktoren, desto mehr Beobachtungen sind für eine “Konvergenz” (hier in etwa: “ein stabiles Modell”) notwendig).


2. Fixed und Random Factors

In einem MM werden die Faktoren in ‘fixed’ und ‘random’ aufgeteilt:

In der Phonetik sind oft Versuchsperson und/oder Item (Materialien wie unterschiedliche Wörter) Random Factors.

Beispiel:


3. Was in einem MM berechnet wird

Ein MM ist in einigen Hinsichten ähnlich wie die lineare Regression

Die 4 Variablen werden auf eine solche Weise berechnet, sodass die Abstände zwischen den tatsächlichen Werten (\(y\)) und den eingeschätzten Werten (\(yhut\)) minimiert werden (also ganz analog zur linearen Regression).

Für das obige Beispiel:

Funktion: lmer() (für linear mixed effects regression model) in library(lmerTest):

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


# Fixed.k, Fixed.m
fixef(stimm.lmer)
## (Intercept)         Kpa 
##   -2.432357    3.900090
# (Intercept)         Kpa
# -2.432357           3.900090

# Random.k, Random.m
# N.B. es gibt einen (k, m) pro Versuchsperson
ranef(stimm.lmer)
## $Vpn
##   (Intercept)        Kpa
## A   13.438581  1.5463273
## B  -16.549964  1.1124893
## C   10.580984 -1.5464906
## D   -6.152659  0.3295627
## E  -19.812548 -0.1910322
## F   13.400584 -0.9499776
## G    1.524253  0.6061824
## H    3.570768 -0.9070613
#
# (Intercept)        Kpa
# A   13.438581  1.5463280
# B  -16.549963  1.1124886
# C   10.580984 -1.5464902
# D   -6.152658  0.3295624
# E  -19.812548 -0.1910332
# F   13.400583 -0.9499770
# G    1.524253  0.6061825
# H    3.570768 -0.9070612

#
# Die ersten 12 eingeschätzten Werte (yhut)
fitted(stimm.lmer)[1:12]
##           1           2           3           4           5           6 
##  11.0062243 -18.9823204   8.1486278 -22.2449049  10.9682269  -0.9081032 
##           7           8           9          10          11          12 
##   1.1384112 -13.9697407  10.5022275  -4.3553621 -18.5358467  13.9183398
#
# Berechnung der eingschätzten Werte
# yhut = (Random.m + Fixed.m) x + (Random.k + Fixed.k)
#
# x ist hier /ba/ oder /pa/ und wird als 0 und 1 umgesetzt
contrasts(stimm$K)
##    pa
## ba  0
## pa  1
# ba  0
# pa  1
# d.h. 0 für ba, 1 für pa.
#
# z.B. yhut der 8en Beobachtung
#
stimm[8, ]
##   vot  K Vpn Alter
## 8 -13 pa   B  Kind
#  vot  K Vpn Alter
#8 -13 pa   B  Kind
#
# Sprecher ist B, K = pa (daher x = 1)
# eingeschätzter Wert:
# yhut = (Random.m + Fixed.m) x + (Random.k + Fixed.k)
yhut = (1.1124886 + 3.900090) * 1 + (-16.549963 - 2.432357)
yhut
## [1] -13.96974
# [1] -13.96974
# das gleiche:
fitted(stimm.lmer)[8]
##         8 
## -13.96974
# -13.96974
# der tatsächliche Wert -13 (siehe 
stimm[8,]
##   vot  K Vpn Alter
## 8 -13 pa   B  Kind
# )

4. Mehr zu Random factors und der Wahl zwischen (FF|RF) (within-RF factor) und (1|RF) (between-RF factor)

Es gibt zwei Möglichkeiten, um einen Random factor zu deklarieren: (FF|RF) oder (1|RF).

head(soa)
##   Vpn     W    Pos    F1
## 1  s1  Bart  final 493.5
## 2  s1  Pfad  final 532.5
## 3  s1 Start  final 514.5
## 4  s1  Bart medial 482.0
## 5  s1  Pfad medial 501.0
## 6  s1 Start medial 497.0

Drei Sprecher produzierten “Bad”, “Pfad”, “Start”, entweder in phrasenmedialer oder -finaler Position. Inwiefern wird F1 im Vokal von der Phrasenposition beeinflusst?

Daher 4 Möglichkeiten: ## 1. soa.lmer = lmer(F1 ~ Pos + (Pos | Vpn) + (Pos | W), data = soa) ## 2. soa.lmer = lmer(F1 ~ Pos + (1 | Vpn) + (Pos | W), data = soa) ## 3. soa.lmer = lmer(F1 ~ Pos + (Pos | Vpn) + (1 | W), data = soa) ## 4. soa.lmer = lmer(F1 ~ Pos + (1 | Vpn) + (1 | W), data = soa)

Grundsätzlich (FF|RF) verwenden (daher 1.), es sei denn FF ist eindeutig ‘between’ im Bezug zu RF. für das obige Beispiel ist FF (Pos) ‘within’ im Bezug zum RF Vpn:

with(soa, table(Vpn, Pos))
##     Pos
## Vpn  final medial
##   s1     3      3
##   s2     3      3
##   s3     3      3

FF ist auch ‘within’ im Bezug zu Wort

with(soa, table(W, Pos))
##        Pos
## W       final medial
##   Bart      3      3
##   Pfad      3      3
##   Start     3      3
# also hier:
soa.lmer = lmer(F1 ~ Pos + (Pos | Vpn) + (Pos | W), data = soa)

Die Bedeutung von (FF|RF)

  1. Eine Steigung \(RF.m\) und Intercept \(RF.k\) werden berechnet

  2. (Wichtiger): die Variabilität vom \(RF\) wird ausgeklammert, unter der Annahme dass \(RF\) und \(FF\) interagieren könnten.

d.h. mit (Pos|Vpn) wird geprueft, inwiefern \(Pos\) und \(Vpn\) miteinander interagieren:

ggplot(soa) +
  aes(y = F1, x = Pos) +
  geom_boxplot() +
  facet_wrap(~ Vpn)

Ebenfalls wird mit (Pos|W) geprüft, inwiefern \(Pos\) und \(W\) miteinander interagieren:

ggplot(soa) +
  aes(y = F1, x = Pos) +
  geom_boxplot() +
  facet_wrap(~ W)

Die Bedeutung von (1|RF)

  1. Ein Intercept \(RF.k\) (aber keine Steigung) wird berechnet.
  2. (Wichtiger): die Variabilität vom \(RF\) wird ausgeklammert, ohne dass \(RF\) und \(FF\) interagieren.

(1|RF) soll verwendet werden, wenn der \(FF\) ‘between’ ist im Bezug zum \(RF\). Z.B.:

head(stimm)
##   vot  K Vpn Alter
## 1  10 ba   A  Kind
## 2 -16 ba   B  Kind
## 3   9 ba   C  Kind
## 4 -21 ba   E   Erw
## 5  10 ba   F   Erw
## 6  -1 ba   G  Kind
# Inwiefern wird vot vom Alter beeinflusst?
# Fixed factor: Alter
# Random factor: Vpn
# Alter ist between im Bezug zu Vpn:
with(stimm, table(Vpn, Alter))
##    Alter
## Vpn Erw Kind
##   A   0   16
##   B   0   16
##   C   0   16
##   D  16    0
##   E  16    0
##   F  16    0
##   G   0   16
##   H   0   16
#
# Alter und Vpn können nicht interagieren, da jede Vpn
# entweder Kind oder Erw ist (wie folgender,
# an sich nicht uebermäßig sinnvoller
# Plot zeigt):
ggplot(stimm) + 
  aes(y = vot, x = Alter) + 
  geom_boxplot() + 
  facet_wrap( ~ Vpn)

Daher wäre (Alter|Vpn) sinnlos - da dies voraussetzt, dass Alter und Vpn interagieren können, was nicht möglich ist, wenn \(FF\) ‘between’ ist in Bezug zum \(RF\). Daher:

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

Im Allgemeinen: immer (FF|RF) verwenden, abgesehen von dem Fall, in dem \(FF\) wirklich ‘between’ ist in Bezug zu \(RF\).

5. Mehr als ein Fixed Factor

dim(stimm)
## [1] 128   4

VOT wurde in /ba, pa/ von Kindern und Erwachsenen produziert.

Inwiefern wird “VOT” von der Artikulationsstelle (Faktor “K”) und/oder von der Sprechergruppe (ob Kind oder Erwachsener: Faktor “Alter”) beeinflusst?

# K ist 'within' im Bezug zu Vpn: daher (K|Vpn)
with(stimm, table(Vpn, K))
##    K
## Vpn ba pa
##   A  8  8
##   B  8  8
##   C  8  8
##   D  8  8
##   E  8  8
##   F  8  8
##   G  8  8
##   H  8  8
# Alter ist 'between' im Bezug zu Vpn daher (1|Vpn).
with(stimm, table(Vpn, Alter))
##    Alter
## Vpn Erw Kind
##   A   0   16
##   B   0   16
##   C   0   16
##   D  16    0
##   E  16    0
##   F  16    0
##   G   0   16
##   H   0   16
# (K|Vpn) enthält aber bereits (1|Vpn)
# Weil (K|Vpn) bedeutet:
# random-slope und random-intercept berechnen;
# und (1|Vpn) bedeutet: nur random-intercept berechnen.
# Daher (K|Vpn) und NICHT (K|Vpn) + (1|Vpn)

# Zusätzlich für zwei fixed factors:
# K * Alter (bedeutet K, Alter, und deren Interaktion).
# K * Alter ist dasselbe wie: K + Alter + K:Alter
stimm.lmer = lmer(vot ~ K * Alter + (K | Vpn), data = stimm)
# Oder
stimm.lmer = lmer(vot ~ K + Alter + K:Alter + (K | Vpn), data = stimm)
# 2.
dim(int)
## [1] 120   5

10 Sprecher (Vpn), davon Männer und Frauen (Faktor G) produzierten zwei Vokale (Faktor V), und die Intensität wurde gemessen (db). Die Vokale kamen in verschiedenen Wörtern vor. Inwiefern wird die Intensität vom Vokal und vom Geschlecht beeinflusst?

# Fixed factors: V, G
# Random factors: Wort, Vpn
#
# V ist 'within' im Bezug zu Wort. Daher (V|Wort)
with(int, table(Wort, V))
##     V
## Wort  a  i
##   w1 10 10
##   w2 10 10
##   w3 10 10
##   w4 10 10
##   w5 10 10
##   w6 10 10
#  V ist 'within' im Bezug zu Vpn. Daher (V|Vpn)
with(int, table(Vpn, V))
##      V
## Vpn   a i
##   S1  6 6
##   S10 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
# G ist 'within' im Bezug zu Wort. Daher (G|Wort)
with(int, table(Wort, G))
##     G
## Wort  f  m
##   w1 10 10
##   w2 10 10
##   w3 10 10
##   w4 10 10
##   w5 10 10
##   w6 10 10
# G ist 'between' im Bezug zu Vpn. Daher (1|Vpn)
with(int, table(Vpn, G))
##      G
## Vpn    f  m
##   S1   0 12
##   S10 12  0
##   S2   0 12
##   S3   0 12
##   S4   0 12
##   S5   0 12
##   S6  12  0
##   S7  12  0
##   S8  12  0
##   S9  12  0

Zusammen: (V+G|Wort) + (V|Vpn) (und NICHT (V|Wort) + (G|Wort) + (V|Vpn) + (1|Vpn)), denn die Faktoren, deren Einfluss ausgeklammert werden soll (hier also Wort und Vpn), sollen nur einmal in der Formel auftauchen.

Daher:

ggplot(int) +
  aes(y = db, x = G, col = V) +
  geom_boxplot()

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

6. Die Prüfstatistik

Die step()-Funktion

  1. vereinfacht das Modell wenn möglich.
  2. prüft, ob die Faktoren in dem vereinfachten Modell signifikant sind.

VOT wurde in /ba, pa/ von Kindern und Erwachsenen produziert. Inwiefern wird VOT von der Artikulationsstelle (Faktor K) und/oder von der Sprechergruppe (ob Kind oder Erwachsen: Faktor Alter) beeinflusst?

# Bild
ggplot(stimm) + 
  aes(y = vot, x = K) + 
  geom_boxplot() + 
  facet_wrap( ~ Alter)