Lineares Modell, kategoriale Prädiktoren, gemischte Modelle (mixed designs)
Unit, Examples, etc
Peter Zezula
Gemischte varianzanalytische Designs
Gemischt werden Messwiederholungsfaktoren (einer oder mehrere) mit Zwischensubjektfaktoren (einer oder mehrere).
Voraussetzungen (wie bei Messwiederholungsdesigns): Zusätzlich zu Varianzenhomogenität über die Faktorstufen kommt noch die Homogenität der Covarianzen (Korrelationen)
Homogenität der Covarianzen wird getestet mit Mauchly’s Test (bei Varianzhomogenität war es Levene-Test). Spherizität (Gleichheit der Varianzen der Differenzen zwischen den Faktorstufen)
Bei Verletzung der Homogenität der Covarianzen:
- Greenhouse-Geisser Korrektur
- Huyn-Feldt Korrektur
ezANOVA aus library(ez)
Ansatz ezANOVA: Beispiel 3x3, STAI-State von verschiedenen Gruppen in drei unterschiedlichen Situationen
Die Personen in den Gruppen gehören entweder zum unteren Drittel, zur Mitte oder zum oberen Drittel hinsichtlich ihrer mathematischen Fähigkeiten. Bei den verglichenen Situationen steht eine neutrale Situation zwei Prüfungssituationen gegenüber.
Im Beispiel:
- Ein dreistufiger Between-Subjects-Faktor (Gruppe).
- Ein dreistufiger Within-Faktor (Situation).
Bei Messwiederholungs- und Mixed-Designs gibt ezANOVA()
den Mauchly’s Test for Sphericity sowie die Greenhouse-Geisser und die
Huyn-Feldt Korrektur mit aus.
dd <- read.delim("http://md.psych.bio.uni-goettingen.de/mv/data/virt/v_stait_exam_wide.txt")
# convert it to long format
# require(reshape2)
# dd.l <- melt(dd, id=c('subj', 'gender', 'group'), variable.name="sit", value.name="anx", measured=c('neutr', 'exam1', 'exam2'))
# head(dd.l)
# convert it to long format using tidyr
require(tidyr)## Lade nötiges Paket: tidyrdd.l <- dd %>% tidyr::gather(key=sit, value=anx, neutr, exam1, exam2)
dd.l$subj <- factor(dd.l$subj)
dd.l$sit <- factor(dd.l$sit)
head(dd.l)## subj gender group sit anx
## 1 1 f av neutr 42
## 2 2 f av neutr 59
## 3 3 f av neutr 49
## 4 4 f av neutr 20
## 5 5 f av neutr 45
## 6 6 f av neutr 39# fit model using ezANOVA, show anova table and
require(ez)## Lade nötiges Paket: ezezANOVA(data=dd.l, dv=.(anx), wid = .(subj), between=.(group), within = .(sit), detailed=TRUE, type=2)## Warning: Converting "group" to factor for ANOVA.## $ANOVA
## Effect DFn DFd SSn SSd F p p<.05
## 1 (Intercept) 1 177 1.252237e+06 82406.856 2.689654e+03 5.740328e-109 *
## 2 group 2 177 3.970208e+04 82406.856 4.263764e+01 7.684359e-16 *
## 3 sit 2 354 5.023511e+03 1682.311 5.285357e+02 5.051346e-107 *
## 4 group:sit 4 354 1.777778e-01 1682.311 9.352214e-03 9.998263e-01
## ges
## 1 9.370744e-01
## 2 3.207180e-01
## 3 5.637257e-02
## 4 2.114153e-06
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 3 sit 0.07418204 3.857354e-100 *
## 4 group:sit 0.07418204 3.857354e-100 *
##
## $`Sphericity Corrections`
## Effect GGe p[GG] p[GG]<.05 HFe p[HF] p[HF]<.05
## 3 sit 0.5192599 4.864142e-57 * 0.5195925 4.491623e-57 *
## 4 group:sit 0.5192599 9.920486e-01 0.5195925 9.920701e-01# descriptives - the ez way
ezStats(data=dd.l, dv=.(anx), wid = .(subj), between=.(group), within = .(sit), type=2)## Warning: Converting "group" to factor for ANOVA.## group sit N Mean SD FLSD
## 1 av exam1 60 49.86667 13.32192 0.7827559
## 2 av exam2 60 50.06667 12.22649 0.7827559
## 3 av neutr 60 43.48333 12.34324 0.7827559
## 4 hi exam1 60 60.83333 14.67751 0.7827559
## 5 hi exam2 60 61.13333 13.05200 0.7827559
## 6 hi neutr 60 54.51667 13.18897 0.7827559
## 7 lo exam1 60 39.86667 11.10881 0.7827559
## 8 lo exam2 60 40.10000 11.55635 0.7827559
## 9 lo neutr 60 33.53333 11.36821 0.7827559# ez offers ezPlot() to visualize effects
ezPlot(data=dd.l, dv=.(anx), wid = .(subj), between=.(group), within = .(sit), x=group, split=.(sit))## Warning: Converting "group" to factor for ANOVA.
ezPlot(data=dd.l, dv=.(anx), wid = .(subj), between=.(group), within = .(sit), x=sit, split=.(group))## Warning: Converting "group" to factor for ANOVA.
‘ges’ im Output: Generalized Eta-Squared measure of effect size see: Bakeman, R. (2005). Recommended effect size statistics for repeated measures designs. Behavior Research Methods, 37 (3), 379-384.
ges Generalized Eta-Squared measure of effect size (see in references below: Bakeman, 2005).
Effect Size Anova \(η^2\)
| small | medium | large |
|---|---|---|
| 0.01 | 0.06 | 0.14 |
Ansatz ezANOVA: Beispiel 2x3x3, STAI-State von verschiedenen Gruppen getrennt nach Geschlecht in drei unterschiedlichen Situationen
Zusätzlich zum obigen Design wird Geschlecht mit aufgenommen.
Im Beispiel:
- Ein dreistufiger Between-Subjects-Faktor (Gruppe).
- Ein zweistufiger Between-Subjects-Faktor (Geschlecht)
- Ein dreistufiger Within-Faktor (Situation).
dd <- read.delim("http://md.psych.bio.uni-goettingen.de/mv/data/virt/v_stait_exam_wide.txt")
# convert it to long format
require(reshape2)## Lade nötiges Paket: reshape2##
## Attache Paket: 'reshape2'## Das folgende Objekt ist maskiert 'package:tidyr':
##
## smithsdd.l <- melt(dd, id=c('subj', 'gender', 'group'), variable.name="sit", value.name="anx", measured=c('neutr', 'exam1', 'exam2'))
# fit model using ezANOVA
require(ez)
ezANOVA(data=dd.l, dv=.(anx), wid = .(subj), between=.(gender, group), within = .(sit), detailed=TRUE, type=2)## Warning: Converting "subj" to factor for ANOVA.## Warning: Converting "gender" to factor for ANOVA.## Warning: Converting "group" to factor for ANOVA.## $ANOVA
## Effect DFn DFd SSn SSd F p
## 1 (Intercept) 1 174 1.252237e+06 66454.4444 3278.7761815 7.929949e-115
## 2 gender 1 174 1.131627e+04 66454.4444 29.6297775 1.756852e-07
## 3 group 2 174 3.970208e+04 66454.4444 51.9766706 2.005531e-18
## 5 sit 2 348 5.023511e+03 149.0889 5862.8844835 9.865761e-269
## 4 gender:group 2 174 4.636144e+03 66454.4444 6.0694897 2.830981e-03
## 6 gender:sit 2 348 1.532844e+03 149.0889 1788.9658667 7.740180e-184
## 7 group:sit 4 348 1.777778e-01 149.0889 0.1037412 9.811539e-01
## 8 gender:group:sit 4 348 3.777778e-01 149.0889 0.2204501 9.269569e-01
## p<.05 ges
## 1 * 9.494984e-01
## 2 * 1.452297e-01
## 3 * 3.734711e-01
## 5 * 7.013428e-02
## 4 * 6.507812e-02
## 6 * 2.249671e-02
## 7 2.669187e-06
## 8 5.672006e-06
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 5 sit 0.558658 1.34317e-22 *
## 6 gender:sit 0.558658 1.34317e-22 *
## 7 group:sit 0.558658 1.34317e-22 *
## 8 gender:group:sit 0.558658 1.34317e-22 *
##
## $`Sphericity Corrections`
## Effect GGe p[GG] p[GG]<.05 HFe p[HF]
## 5 sit 0.6937979 2.030140e-187 * 0.6976012 1.983320e-188
## 6 gender:sit 0.6937979 1.642905e-128 * 0.6976012 3.374903e-129
## 7 group:sit 0.6937979 9.492048e-01 0.6976012 9.498437e-01
## 8 gender:group:sit 0.6937979 8.681917e-01 0.6976012 8.691952e-01
## p[HF]<.05
## 5 *
## 6 *
## 7
## 8# descriptives
ezStats(data=dd.l, dv=.(anx), wid = .(subj), between=.(gender, group), within = .(sit), type=2)## Warning: Converting "subj" to factor for ANOVA.## Warning: Converting "gender" to factor for ANOVA.## Warning: Converting "group" to factor for ANOVA.## gender group sit N Mean SD FLSD
## 1 f av neutr 30 48.10000 12.644612 0.3323904
## 2 f av exam1 30 57.76667 12.053568 0.3323904
## 3 f av exam2 30 54.23333 12.229651 0.3323904
## 4 f hi neutr 30 61.16667 12.605646 0.3323904
## 5 f hi exam1 30 70.80000 12.058478 0.3323904
## 6 f hi exam2 30 67.23333 12.316749 0.3323904
## 7 f lo neutr 30 33.16667 12.605646 0.3323904
## 8 f lo exam1 30 42.83333 12.026073 0.3323904
## 9 f lo exam2 30 39.30000 12.304050 0.3323904
## 10 m av neutr 30 38.86667 10.294737 0.3323904
## 11 m av exam1 30 41.96667 9.308814 0.3323904
## 12 m av exam2 30 45.90000 10.892421 0.3323904
## 13 m hi neutr 30 47.86667 10.173438 0.3323904
## 14 m hi exam1 30 50.86667 9.346706 0.3323904
## 15 m hi exam2 30 55.03333 10.857965 0.3323904
## 16 m lo neutr 30 33.90000 10.185690 0.3323904
## 17 m lo exam1 30 36.90000 9.393138 0.3323904
## 18 m lo exam2 30 40.90000 10.908238 0.3323904# sequence matters for output
ezStats(data=dd.l, dv=.(anx), wid = .(subj), between=.(group, gender), within = .(sit), type=2)## Warning: Converting "subj" to factor for ANOVA.
## Warning: Converting "group" to factor for ANOVA.## Warning: Converting "gender" to factor for ANOVA.## group gender sit N Mean SD FLSD
## 1 av f neutr 30 48.10000 12.644612 0.3323904
## 2 av f exam1 30 57.76667 12.053568 0.3323904
## 3 av f exam2 30 54.23333 12.229651 0.3323904
## 4 av m neutr 30 38.86667 10.294737 0.3323904
## 5 av m exam1 30 41.96667 9.308814 0.3323904
## 6 av m exam2 30 45.90000 10.892421 0.3323904
## 7 hi f neutr 30 61.16667 12.605646 0.3323904
## 8 hi f exam1 30 70.80000 12.058478 0.3323904
## 9 hi f exam2 30 67.23333 12.316749 0.3323904
## 10 hi m neutr 30 47.86667 10.173438 0.3323904
## 11 hi m exam1 30 50.86667 9.346706 0.3323904
## 12 hi m exam2 30 55.03333 10.857965 0.3323904
## 13 lo f neutr 30 33.16667 12.605646 0.3323904
## 14 lo f exam1 30 42.83333 12.026073 0.3323904
## 15 lo f exam2 30 39.30000 12.304050 0.3323904
## 16 lo m neutr 30 33.90000 10.185690 0.3323904
## 17 lo m exam1 30 36.90000 9.393138 0.3323904
## 18 lo m exam2 30 40.90000 10.908238 0.3323904# graphics
ezPlot(data=dd.l, dv=.(anx), wid = .(subj), between=.(gender, group), within = .(sit), x=group, split=.(sit))## Warning: Converting "subj" to factor for ANOVA.
## Warning: Converting "gender" to factor for ANOVA.## Warning: Converting "group" to factor for ANOVA.
# col gives us multiple plots
ezPlot(data=dd.l, dv=.(anx), wid = .(subj), between=.(gender, group), within = .(sit), x=group, col=.(gender), split=.(sit))## Warning: Converting "subj" to factor for ANOVA.## Warning: Converting "gender" to factor for ANOVA.## Warning: Converting "group" to factor for ANOVA.
ezPlot(data=dd.l, dv=.(anx), wid = .(subj), between=.(gender, group), within = .(sit), x=group, col=.(sit), split=.(gender))## Warning: Converting "subj" to factor for ANOVA.## Warning: Converting "gender" to factor for ANOVA.## Warning: Converting "group" to factor for ANOVA.
Ansatz ezANOVA: detaillierte Informationen über das aov-Objekt
Beispiel 2x3x3, STAI-State von verschiedenen Gruppen getrennt nach Geschlecht in drei unterschiedlichen Situationen
Im Beispiel:
- Ein dreistufiger Between-Subjects-Faktor (Gruppe).
- Ein zweistufiger Between-Subjects-Faktor (Geschlecht)
- Ein dreistufiger Within-Faktor (Situation).
Wir generieren ein aov-Objekt über den Parameter
return_aov = TRUE. Dies gibt uns die Möglichkeit, Details
des generierten Modells abzurufen.
dd <- read.delim("http://md.psych.bio.uni-goettingen.de/mv/data/virt/v_stait_exam_wide.txt")
# convert it to long format
require(reshape2)
dd.l <- melt(dd, id=c('subj', 'gender', 'group'), variable.name="sit", value.name="anx", measured=c('neutr', 'exam1', 'exam2'))
# fit model using ezANOVA
require(ez)
m.ez <- ezANOVA(data=dd.l, dv=.(anx), wid = .(subj), between=.(gender, group), within = .(sit), detailed=TRUE, type=2, return_aov = TRUE)## Warning: Converting "subj" to factor for ANOVA.## Warning: Converting "gender" to factor for ANOVA.## Warning: Converting "group" to factor for ANOVA.# summary shows list elements we can access
summary(m.ez)## Length Class Mode
## ANOVA 9 data.frame list
## Mauchly's Test for Sphericity 4 data.frame list
## Sphericity Corrections 7 data.frame list
## aov 3 aovlist list# we access elements by number
m.ez[2]## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 5 sit 0.558658 1.34317e-22 *
## 6 gender:sit 0.558658 1.34317e-22 *
## 7 group:sit 0.558658 1.34317e-22 *
## 8 gender:group:sit 0.558658 1.34317e-22 *# or by name
m.ez['ANOVA']## $ANOVA
## Effect DFn DFd SSn SSd F p
## 1 (Intercept) 1 174 1.252237e+06 66454.4444 3278.7761815 7.929949e-115
## 2 gender 1 174 1.131627e+04 66454.4444 29.6297775 1.756852e-07
## 3 group 2 174 3.970208e+04 66454.4444 51.9766706 2.005531e-18
## 5 sit 2 348 5.023511e+03 149.0889 5862.8844835 9.865761e-269
## 4 gender:group 2 174 4.636144e+03 66454.4444 6.0694897 2.830981e-03
## 6 gender:sit 2 348 1.532844e+03 149.0889 1788.9658667 7.740180e-184
## 7 group:sit 4 348 1.777778e-01 149.0889 0.1037412 9.811539e-01
## 8 gender:group:sit 4 348 3.777778e-01 149.0889 0.2204501 9.269569e-01
## p<.05 ges
## 1 * 9.494984e-01
## 2 * 1.452297e-01
## 3 * 3.734711e-01
## 5 * 7.013428e-02
## 4 * 6.507812e-02
## 6 * 2.249671e-02
## 7 2.669187e-06
## 8 5.672006e-06# we can get the model.matrix for our design.
m.mx <- model.matrix(m.ez['aov']$aov)
# As it has subj-1 dummy columns to code effect subject, it is big
# columns are:
m.mx[1,]## (Intercept) genderm grouphi
## 1 0 0
## grouplo sitexam1 sitexam2
## 0 0 0
## Error(subj/(sit))2 Error(subj/(sit))3 Error(subj/(sit))4
## 0 0 0
## Error(subj/(sit))5 Error(subj/(sit))6 Error(subj/(sit))7
## 0 0 0
## Error(subj/(sit))8 Error(subj/(sit))9 Error(subj/(sit))10
## 0 0 0
## Error(subj/(sit))11 Error(subj/(sit))12 Error(subj/(sit))13
## 0 0 0
## Error(subj/(sit))14 Error(subj/(sit))15 Error(subj/(sit))16
## 0 0 0
## Error(subj/(sit))17 Error(subj/(sit))18 Error(subj/(sit))19
## 0 0 0
## Error(subj/(sit))20 Error(subj/(sit))21 Error(subj/(sit))22
## 0 0 0
## Error(subj/(sit))23 Error(subj/(sit))24 Error(subj/(sit))25
## 0 0 0
## Error(subj/(sit))26 Error(subj/(sit))27 Error(subj/(sit))28
## 0 0 0
## Error(subj/(sit))29 Error(subj/(sit))30 Error(subj/(sit))31
## 0 0 0
## Error(subj/(sit))32 Error(subj/(sit))33 Error(subj/(sit))34
## 0 0 0
## Error(subj/(sit))35 Error(subj/(sit))36 Error(subj/(sit))37
## 0 0 0
## Error(subj/(sit))38 Error(subj/(sit))39 Error(subj/(sit))40
## 0 0 0
## Error(subj/(sit))41 Error(subj/(sit))42 Error(subj/(sit))43
## 0 0 0
## Error(subj/(sit))44 Error(subj/(sit))45 Error(subj/(sit))46
## 0 0 0
## Error(subj/(sit))47 Error(subj/(sit))48 Error(subj/(sit))49
## 0 0 0
## Error(subj/(sit))50 Error(subj/(sit))51 Error(subj/(sit))52
## 0 0 0
## Error(subj/(sit))53 Error(subj/(sit))54 Error(subj/(sit))55
## 0 0 0
## Error(subj/(sit))56 Error(subj/(sit))57 Error(subj/(sit))58
## 0 0 0
## Error(subj/(sit))59 Error(subj/(sit))60 Error(subj/(sit))61
## 0 0 0
## Error(subj/(sit))62 Error(subj/(sit))63 Error(subj/(sit))64
## 0 0 0
## Error(subj/(sit))65 Error(subj/(sit))66 Error(subj/(sit))67
## 0 0 0
## Error(subj/(sit))68 Error(subj/(sit))69 Error(subj/(sit))70
## 0 0 0
## Error(subj/(sit))71 Error(subj/(sit))72 Error(subj/(sit))73
## 0 0 0
## Error(subj/(sit))74 Error(subj/(sit))75 Error(subj/(sit))76
## 0 0 0
## Error(subj/(sit))77 Error(subj/(sit))78 Error(subj/(sit))79
## 0 0 0
## Error(subj/(sit))80 Error(subj/(sit))81 Error(subj/(sit))82
## 0 0 0
## Error(subj/(sit))83 Error(subj/(sit))84 Error(subj/(sit))85
## 0 0 0
## Error(subj/(sit))86 Error(subj/(sit))87 Error(subj/(sit))88
## 0 0 0
## Error(subj/(sit))89 Error(subj/(sit))90 Error(subj/(sit))91
## 0 0 0
## Error(subj/(sit))92 Error(subj/(sit))93 Error(subj/(sit))94
## 0 0 0
## Error(subj/(sit))95 Error(subj/(sit))96 Error(subj/(sit))97
## 0 0 0
## Error(subj/(sit))98 Error(subj/(sit))99 Error(subj/(sit))100
## 0 0 0
## Error(subj/(sit))101 Error(subj/(sit))102 Error(subj/(sit))103
## 0 0 0
## Error(subj/(sit))104 Error(subj/(sit))105 Error(subj/(sit))106
## 0 0 0
## Error(subj/(sit))107 Error(subj/(sit))108 Error(subj/(sit))109
## 0 0 0
## Error(subj/(sit))110 Error(subj/(sit))111 Error(subj/(sit))112
## 0 0 0
## Error(subj/(sit))113 Error(subj/(sit))114 Error(subj/(sit))115
## 0 0 0
## Error(subj/(sit))116 Error(subj/(sit))117 Error(subj/(sit))118
## 0 0 0
## Error(subj/(sit))119 Error(subj/(sit))120 Error(subj/(sit))121
## 0 0 0
## Error(subj/(sit))122 Error(subj/(sit))123 Error(subj/(sit))124
## 0 0 0
## Error(subj/(sit))125 Error(subj/(sit))126 Error(subj/(sit))127
## 0 0 0
## Error(subj/(sit))128 Error(subj/(sit))129 Error(subj/(sit))130
## 0 0 0
## Error(subj/(sit))131 Error(subj/(sit))132 Error(subj/(sit))133
## 0 0 0
## Error(subj/(sit))134 Error(subj/(sit))135 Error(subj/(sit))136
## 0 0 0
## Error(subj/(sit))137 Error(subj/(sit))138 Error(subj/(sit))139
## 0 0 0
## Error(subj/(sit))140 Error(subj/(sit))141 Error(subj/(sit))142
## 0 0 0
## Error(subj/(sit))143 Error(subj/(sit))144 Error(subj/(sit))145
## 0 0 0
## Error(subj/(sit))146 Error(subj/(sit))147 Error(subj/(sit))148
## 0 0 0
## Error(subj/(sit))149 Error(subj/(sit))150 Error(subj/(sit))151
## 0 0 0
## Error(subj/(sit))152 Error(subj/(sit))153 Error(subj/(sit))154
## 0 0 0
## Error(subj/(sit))155 Error(subj/(sit))156 Error(subj/(sit))157
## 0 0 0
## Error(subj/(sit))158 Error(subj/(sit))159 Error(subj/(sit))160
## 0 0 0
## Error(subj/(sit))161 Error(subj/(sit))162 Error(subj/(sit))163
## 0 0 0
## Error(subj/(sit))164 Error(subj/(sit))165 Error(subj/(sit))166
## 0 0 0
## Error(subj/(sit))167 Error(subj/(sit))168 Error(subj/(sit))169
## 0 0 0
## Error(subj/(sit))170 Error(subj/(sit))171 Error(subj/(sit))172
## 0 0 0
## Error(subj/(sit))173 Error(subj/(sit))174 Error(subj/(sit))175
## 0 0 0
## Error(subj/(sit))176 Error(subj/(sit))177 Error(subj/(sit))178
## 0 0 0
## Error(subj/(sit))179 Error(subj/(sit))180 genderm:grouphi
## 0 0 0
## genderm:grouplo genderm:sitexam1 genderm:sitexam2
## 0 0 0
## grouphi:sitexam1 grouplo:sitexam1 grouphi:sitexam2
## 0 0 0
## grouplo:sitexam2 genderm:grouphi:sitexam1 genderm:grouplo:sitexam1
## 0 0 0
## genderm:grouphi:sitexam2 genderm:grouplo:sitexam2
## 0 0# we can look at some significant parts to understand dummy coding
m.mx[c(2:5, 535:540),c(1:8, 185:187, 196:197)]## (Intercept) genderm grouphi grouplo sitexam1 sitexam2 Error(subj/(sit))2
## 2 1 0 0 0 1 0 0
## 3 1 0 0 0 0 1 0
## 4 1 0 0 0 0 0 1
## 5 1 0 0 0 1 0 1
## 535 1 1 1 0 0 0 0
## 536 1 1 1 0 1 0 0
## 537 1 1 1 0 0 1 0
## 538 1 1 1 0 0 0 0
## 539 1 1 1 0 1 0 0
## 540 1 1 1 0 0 1 0
## Error(subj/(sit))3 Error(subj/(sit))180 genderm:grouphi genderm:grouplo
## 2 0 0 0 0
## 3 0 0 0 0
## 4 0 0 0 0
## 5 0 0 0 0
## 535 0 0 1 0
## 536 0 0 1 0
## 537 0 0 1 0
## 538 0 1 1 0
## 539 0 1 1 0
## 540 0 1 1 0
## genderm:grouphi:sitexam2 genderm:grouplo:sitexam2
## 2 0 0
## 3 0 0
## 4 0 0
## 5 0 0
## 535 0 0
## 536 0 0
## 537 1 0
## 538 0 0
## 539 0 0
## 540 1 0# ... the major part of subject coding and some of the interaction codings are cut out
# we store aov element for convenience
o.aov <- m.ez[['aov']]
# and acccess some details of effect "subj"
o.aov[['subj']]$coefficients## genderm grouphi grouplo
## -11.122222 13.033333 -14.933333
## genderm:grouphi genderm:grouplo genderm:sitexam1
## -4.022222 9.922222 NA
## genderm:sitexam2 grouphi:sitexam1 grouplo:sitexam1
## NA NA NA
## grouphi:sitexam2 grouplo:sitexam2 genderm:grouphi:sitexam1
## NA NA NA
## genderm:grouplo:sitexam1 genderm:grouphi:sitexam2 genderm:grouplo:sitexam2
## NA NA NAhead(o.aov[['subj']]$residuals)## 2 3 4 5 6 7
## 20.412415 -1.885618 -43.833333 3.743884 -6.429965 7.394228head(o.aov[['subj']]$fitted.values)## 2 3 4 5 6
## 1.172396e-13 1.414424e-13 7.815970e-14 6.794565e-14 9.237056e-14
## 7
## -7.016610e-14head(o.aov[['subj']]$effects)## genderm grouphi grouplo genderm:grouphi genderm:grouplo
## -106.377943 175.362505 -94.604807 49.203743 -47.065233
##
## 7.394228# descriptives
ezStats(data=dd.l, dv=.(anx), wid = .(subj), between=.(gender, group), within = .(sit), type=2)## Warning: Converting "subj" to factor for ANOVA.## Warning: Converting "gender" to factor for ANOVA.## Warning: Converting "group" to factor for ANOVA.## gender group sit N Mean SD FLSD
## 1 f av neutr 30 48.10000 12.644612 0.3323904
## 2 f av exam1 30 57.76667 12.053568 0.3323904
## 3 f av exam2 30 54.23333 12.229651 0.3323904
## 4 f hi neutr 30 61.16667 12.605646 0.3323904
## 5 f hi exam1 30 70.80000 12.058478 0.3323904
## 6 f hi exam2 30 67.23333 12.316749 0.3323904
## 7 f lo neutr 30 33.16667 12.605646 0.3323904
## 8 f lo exam1 30 42.83333 12.026073 0.3323904
## 9 f lo exam2 30 39.30000 12.304050 0.3323904
## 10 m av neutr 30 38.86667 10.294737 0.3323904
## 11 m av exam1 30 41.96667 9.308814 0.3323904
## 12 m av exam2 30 45.90000 10.892421 0.3323904
## 13 m hi neutr 30 47.86667 10.173438 0.3323904
## 14 m hi exam1 30 50.86667 9.346706 0.3323904
## 15 m hi exam2 30 55.03333 10.857965 0.3323904
## 16 m lo neutr 30 33.90000 10.185690 0.3323904
## 17 m lo exam1 30 36.90000 9.393138 0.3323904
## 18 m lo exam2 30 40.90000 10.908238 0.3323904library(rstatix)
library(rstatix) ist Pipe-friendly
# we can do classical mixed factorial design using library(rstatix)
# we stay with omnibus tests for all effects
# library(rstatix)
library(tidyverse)## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──## ✔ ggplot2 3.3.5 ✔ dplyr 1.0.8
## ✔ tibble 3.1.6 ✔ stringr 1.4.0
## ✔ readr 2.1.2 ✔ forcats 0.5.1
## ✔ purrr 0.3.4## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()dd.w <- read.delim("http://md.psych.bio.uni-goettingen.de/mv/data/virt/v_stait_exam_wide.txt")
# convert it to long format
# require(reshape2)
dd <- reshape2::melt(dd.w, id=c('subj', 'gender', 'group'), variable.name="sit", value.name="anx", measured=c('neutr', 'exam1', 'exam2'))
dd %>% rstatix::factorial_design(dv = anx, wid = subj, within = sit, between = c(group, gender)) -> design
design # the defined design ## $dv
## [1] "anx"
##
## $between
## [1] "group" "gender"
##
## $wid
## [1] "subj"
##
## $within
## [1] "sit"
##
## $data
## # A tibble: 540 × 5
## subj gender group sit anx
## <fct> <fct> <fct> <fct> <int>
## 1 1 f av neutr 42
## 2 2 f av neutr 59
## 3 3 f av neutr 49
## 4 4 f av neutr 20
## 5 5 f av neutr 45
## 6 6 f av neutr 39
## 7 7 f av neutr 47
## 8 8 f av neutr 36
## 9 9 f av neutr 43
## 10 10 f av neutr 54
## # … with 530 more rows
##
## $idata
## sit
## 1 neutr
## 2 exam1
## 3 exam2
##
## $idesign
## ~sit
## <environment: 0x7fd21e668858>
##
## $lm_data
## # A tibble: 180 × 6
## subj group gender exam1 exam2 neutr
## <fct> <fct> <fct> <int> <int> <int>
## 1 1 av f 52 48 42
## 2 2 av f 68 65 59
## 3 3 av f 59 55 49
## 4 4 av f 31 27 20
## 5 5 av f 55 51 45
## 6 6 av f 49 45 39
## 7 7 av f 57 53 47
## 8 8 av f 47 43 36
## 9 9 av f 52 49 43
## 10 10 av f 63 60 54
## # … with 170 more rows
##
## $repeated
## [1] TRUE
##
## $lm_formula
## cbind(exam1, exam2, neutr) ~ group * gender
## <environment: 0x7fd22e7c63c8>
##
## $model
##
## Call:
## stats::lm(formula = lm_formula, data = data)
##
## Coefficients:
## exam1 exam2 neutr
## (Intercept) 50.1889 50.4333 43.8444
## group1 -0.3222 -0.3667 -0.3611
## group2 10.6444 10.7000 10.6722
## gender1 6.9444 3.1556 3.6333
## group1:gender1 0.9556 1.0111 0.9833
## group2:gender1 3.0222 2.9444 3.0167r.f <- car::Anova(design$model, idata = design$idata, idesign = design$idesign, type = 3)
summary(r.f, multivariate = FALSE)##
## Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
##
## Sum Sq num Df Error SS den Df F value Pr(>F)
## (Intercept) 1252237 1 66454 174 3278.7762 < 2.2e-16 ***
## group 39702 2 66454 174 51.9767 < 2.2e-16 ***
## gender 11316 1 66454 174 29.6298 1.757e-07 ***
## group:gender 4636 2 66454 174 6.0695 0.002831 **
## sit 5024 2 149 348 5862.8845 < 2.2e-16 ***
## group:sit 0 4 149 348 0.1037 0.981154
## gender:sit 1533 2 149 348 1788.9659 < 2.2e-16 ***
## group:gender:sit 0 4 149 348 0.2205 0.926957
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Mauchly Tests for Sphericity
##
## Test statistic p-value
## sit 0.55866 1.3432e-22
## group:sit 0.55866 1.3432e-22
## gender:sit 0.55866 1.3432e-22
## group:gender:sit 0.55866 1.3432e-22
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## sit 0.6938 <2e-16 ***
## group:sit 0.6938 0.9492
## gender:sit 0.6938 <2e-16 ***
## group:gender:sit 0.6938 0.8682
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## HF eps Pr(>F[HF])
## sit 0.6976012 1.983320e-188
## group:sit 0.6976012 9.498437e-01
## gender:sit 0.6976012 3.374903e-129
## group:gender:sit 0.6976012 8.691952e-01library(afex)
Wrapper Package, nutzt library(lme4) bei den Multilevel-Teilen
library(tidyverse)
dd.w <- read.delim("http://md.psych.bio.uni-goettingen.de/mv/data/virt/v_stait_exam_wide.txt")
# convert it to long format
# require(reshape2)
dd <- reshape2::melt(dd.w, id=c('subj', 'gender', 'group'), variable.name="sit", value.name="anx", measured=c('neutr', 'exam1', 'exam2'))
# pure repeated measure design
m.1 <- afex::mixed(anx ~ sit + (1|sit) + (1|subj), data=dd)## Contrasts set to contr.sum for the following variables: sit## 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: -8.8e-08m.1## Warning: lme4 reported (at least) the following warnings for 'full':
## * unable to evaluate scaled gradient
## * Model failed to converge: degenerate Hessian with 1 negative eigenvalues## Mixed Model Anova Table (Type 3 tests, S-method)
##
## Model: anx ~ sit + (1 | sit) + (1 | subj)
## Data: dd
## Effect df F p.value
## 1 sit 2, 358.00 0.22 .803
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1# mixed design
m.2 <- afex::mixed(anx ~ gender + group + sit + (1|sit) + (1|subj), data=dd)## Contrasts set to contr.sum for the following variables: gender, group, sit## 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: -2.5e-08m.2## Warning: lme4 reported (at least) the following warnings for 'full':
## * unable to evaluate scaled gradient
## * Model failed to converge: degenerate Hessian with 1 negative eigenvalues## Mixed Model Anova Table (Type 3 tests, S-method)
##
## Model: anx ~ gender + group + sit + (1 | sit) + (1 | subj)
## Data: dd
## Effect df F p.value
## 1 gender 1, 176.00 28.02 *** <.001
## 2 group 2, 176.00 49.15 *** <.001
## 3 sit 2, 358.00 0.22 .800
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1summary(m.2)$varcor## Groups Name Std.Dev.
## subj (Intercept) 11.5358
## sit (Intercept) 7.9034
## Residual 2.1679# we stay with the problem of having omnibus tests onlyMultilevel Ansatz mit library(nlme) und dem Befehl
nlme::lme().
Wir können Daten mit Multi-Level Ansätzen modellieren, wenn sie hierarchisch organisiert sind. Hier sind gruppenspezifische Parameterschätzungen möglich. Auch Gruppen-Variablen können in die Modelle mit einbezogen werden.
Wir können “klassische” varianzanalytische Messwiederholungsdesigns und Mixed-Designs als Multi-Level-Designs betrachten. Sie sind dann ein Spezialfall einer Multi-Level-Analyse, in der nur random Intercepts betrachtet werden.
Multi-Level-Ansatz wie in Unit ml
Benutzt wird lme() aus dem
package(nlme).
Ein dreistufiger Between-Subjects-Faktor (Gruppe). Ein zweistufiger Between-Subjects-Faktor (Geschlecht). Ein dreistufiger Within-Faktor (Situation).
dd <- read.delim("http://md.psych.bio.uni-goettingen.de/mv/data/virt/v_stait_exam_wide.txt")
# convert it to long format
# require(reshape2)
# dd.l <- melt(dd, id=c('subj', 'gender', 'group'), variable.name="sit", value.name="anx", measured=c('neutr', 'exam1', 'exam2'))
# convert it to long format using tidyr
require(tidyr)
dd.l <- dd %>% tidyr::gather(key=sit, value=anx, neutr, exam1, exam2)
dd.l$subj <- factor(dd.l$subj)
dd.l$sit <- factor(dd.l$sit)
head(dd.l)## subj gender group sit anx
## 1 1 f av neutr 42
## 2 2 f av neutr 59
## 3 3 f av neutr 49
## 4 4 f av neutr 20
## 5 5 f av neutr 45
## 6 6 f av neutr 39# fit model using nlme::lme
require(nlme)## Lade nötiges Paket: nlme##
## Attache Paket: 'nlme'## Das folgende Objekt ist maskiert 'package:dplyr':
##
## collapse# m.1 <- lme(fixed = anx ~ 1 + gender * group, random = ~1 + sit | subj, data=dd.l, method='ML')
# null model to compare with
m.0 <- lme(fixed= anx ~ 1, random = ~1 | subj, data=dd.l, method='ML')
# random intercept model, this is equivqlent to ezANOVA example above
m.1 <- lme(fixed = anx ~ 1 + gender * group * sit, random = ~1 | subj, data=dd.l, method='ML')
# we may also adapt a random intercept random slope model:
m.2 <- lme(fixed = anx ~ 1 + gender * group * sit, random = ~1 + sit | subj, data=dd.l, method='ML')
# compare models
anova(m.0, m.1)## Model df AIC BIC logLik Test L.Ratio p-value
## m.0 1 3 3764.877 3777.752 -1879.439
## m.1 2 20 2319.128 2404.960 -1139.564 1 vs 2 1479.749 <.0001# ... the ANOVA equvalent model differs from the zero model
anova(m.1, m.2)## Model df AIC BIC logLik Test L.Ratio p-value
## m.1 1 20 2319.128 2404.960 -1139.5641
## m.2 2 25 2014.035 2121.324 -982.0176 1 vs 2 315.0931 <.0001# ... integrating random slope in m.2 is even better and significant different from m.1
# see details
summary(m.1)## Linear mixed-effects model fit by maximum likelihood
## Data: dd.l
## AIC BIC logLik
## 2319.128 2404.96 -1139.564
##
## Random effects:
## Formula: ~1 | subj
## (Intercept) Residual
## StdDev: 11.08719 0.6435338
##
## Fixed effects: anx ~ 1 + gender * group * sit
## Value Std.Error DF t-value p-value
## (Intercept) 57.76667 2.062304 348 28.01074 0.0000
## genderm -15.80000 2.916539 174 -5.41738 0.0000
## grouphi 13.03333 2.916539 174 4.46877 0.0000
## grouplo -14.93333 2.916539 174 -5.12022 0.0000
## sitexam2 -3.53333 0.169000 348 -20.90726 0.0000
## sitneutr -9.66667 0.169000 348 -57.19912 0.0000
## genderm:grouphi -4.13333 4.124609 174 -1.00212 0.3177
## genderm:grouplo 9.86667 4.124609 174 2.39215 0.0178
## genderm:sitexam2 7.46667 0.239002 348 31.24096 0.0000
## genderm:sitneutr 6.56667 0.239002 348 27.47531 0.0000
## grouphi:sitexam2 -0.03333 0.239002 348 -0.13947 0.8892
## grouplo:sitexam2 0.00000 0.239002 348 0.00000 1.0000
## grouphi:sitneutr 0.03333 0.239002 348 0.13947 0.8892
## grouplo:sitneutr 0.00000 0.239002 348 0.00000 1.0000
## genderm:grouphi:sitexam2 0.26667 0.338001 348 0.78895 0.4307
## genderm:grouplo:sitexam2 0.06667 0.338001 348 0.19724 0.8438
## genderm:grouphi:sitneutr 0.06667 0.338001 348 0.19724 0.8438
## genderm:grouplo:sitneutr 0.10000 0.338001 348 0.29586 0.7675
## Correlation:
## (Intr) gendrm grouph groupl sitxm2 sitntr gndrm:grph
## genderm -0.707
## grouphi -0.707 0.500
## grouplo -0.707 0.500 0.500
## sitexam2 -0.041 0.029 0.029 0.029
## sitneutr -0.041 0.029 0.029 0.029 0.500
## genderm:grouphi 0.500 -0.707 -0.707 -0.354 -0.020 -0.020
## genderm:grouplo 0.500 -0.707 -0.354 -0.707 -0.020 -0.020 0.500
## genderm:sitexam2 0.029 -0.041 -0.020 -0.020 -0.707 -0.354 0.029
## genderm:sitneutr 0.029 -0.041 -0.020 -0.020 -0.354 -0.707 0.029
## grouphi:sitexam2 0.029 -0.020 -0.041 -0.020 -0.707 -0.354 0.029
## grouplo:sitexam2 0.029 -0.020 -0.020 -0.041 -0.707 -0.354 0.014
## grouphi:sitneutr 0.029 -0.020 -0.041 -0.020 -0.354 -0.707 0.029
## grouplo:sitneutr 0.029 -0.020 -0.020 -0.041 -0.354 -0.707 0.014
## genderm:grouphi:sitexam2 -0.020 0.029 0.029 0.014 0.500 0.250 -0.041
## genderm:grouplo:sitexam2 -0.020 0.029 0.014 0.029 0.500 0.250 -0.020
## genderm:grouphi:sitneutr -0.020 0.029 0.029 0.014 0.250 0.500 -0.041
## genderm:grouplo:sitneutr -0.020 0.029 0.014 0.029 0.250 0.500 -0.020
## gndrm:grpl gndr:2 gndrm:s grph:2 grpl:2 grph:s grpl:s
## genderm
## grouphi
## grouplo
## sitexam2
## sitneutr
## genderm:grouphi
## genderm:grouplo
## genderm:sitexam2 0.029
## genderm:sitneutr 0.029 0.500
## grouphi:sitexam2 0.014 0.500 0.250
## grouplo:sitexam2 0.029 0.500 0.250 0.500
## grouphi:sitneutr 0.014 0.250 0.500 0.500 0.250
## grouplo:sitneutr 0.029 0.250 0.500 0.250 0.500 0.500
## genderm:grouphi:sitexam2 -0.020 -0.707 -0.354 -0.707 -0.354 -0.354 -0.177
## genderm:grouplo:sitexam2 -0.041 -0.707 -0.354 -0.354 -0.707 -0.177 -0.354
## genderm:grouphi:sitneutr -0.020 -0.354 -0.707 -0.354 -0.177 -0.707 -0.354
## genderm:grouplo:sitneutr -0.041 -0.354 -0.707 -0.177 -0.354 -0.354 -0.707
## gndrm:grph:2 gndrm:grpl:2 gndrm:grph:
## genderm
## grouphi
## grouplo
## sitexam2
## sitneutr
## genderm:grouphi
## genderm:grouplo
## genderm:sitexam2
## genderm:sitneutr
## grouphi:sitexam2
## grouplo:sitexam2
## grouphi:sitneutr
## grouplo:sitneutr
## genderm:grouphi:sitexam2
## genderm:grouplo:sitexam2 0.500
## genderm:grouphi:sitneutr 0.500 0.250
## genderm:grouplo:sitneutr 0.250 0.500 0.500
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.97954106 -0.46734820 0.00540358 0.41253385 2.97715166
##
## Number of Observations: 540
## Number of Groups: 180summary(m.2)## Linear mixed-effects model fit by maximum likelihood
## Data: dd.l
## AIC BIC logLik
## 2014.035 2121.324 -982.0176
##
## Random effects:
## Formula: ~1 + sit | subj
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## (Intercept) 10.5977159 (Intr) sitxm2
## sitexam2 1.0972113 0.715
## sitneutr 0.7354644 0.930 0.908
## Residual 0.2770713
##
## Fixed effects: anx ~ 1 + gender * group * sit
## Value Std.Error DF t-value p-value
## (Intercept) 57.76667 1.968619 348 29.34375 0.0000
## genderm -15.80000 2.784048 174 -5.67519 0.0000
## grouphi 13.03333 2.784048 174 4.68143 0.0000
## grouplo -14.93333 2.784048 174 -5.36389 0.0000
## sitexam2 -3.53333 0.216350 348 -16.33158 0.0000
## sitneutr -9.66667 0.154746 348 -62.46790 0.0000
## genderm:grouphi -4.13333 3.937238 174 -1.04981 0.2953
## genderm:grouplo 9.86667 3.937238 174 2.50599 0.0131
## genderm:sitexam2 7.46667 0.305965 348 24.40368 0.0000
## genderm:sitneutr 6.56667 0.218844 348 30.00614 0.0000
## grouphi:sitexam2 -0.03333 0.305965 348 -0.10894 0.9133
## grouplo:sitexam2 0.00000 0.305965 348 0.00000 1.0000
## grouphi:sitneutr 0.03333 0.218844 348 0.15232 0.8790
## grouplo:sitneutr 0.00000 0.218844 348 0.00000 1.0000
## genderm:grouphi:sitexam2 0.26667 0.432700 348 0.61629 0.5381
## genderm:grouplo:sitexam2 0.06667 0.432700 348 0.15407 0.8776
## genderm:grouphi:sitneutr 0.06667 0.309492 348 0.21541 0.8296
## genderm:grouplo:sitneutr 0.10000 0.309492 348 0.32311 0.7468
## Correlation:
## (Intr) gendrm grouph groupl sitxm2 sitntr gndrm:grph
## genderm -0.707
## grouphi -0.707 0.500
## grouplo -0.707 0.500 0.500
## sitexam2 0.667 -0.471 -0.471 -0.471
## sitneutr 0.812 -0.574 -0.574 -0.574 0.834
## genderm:grouphi 0.500 -0.707 -0.707 -0.354 0.333 0.406
## genderm:grouplo 0.500 -0.707 -0.354 -0.707 0.333 0.406 0.500
## genderm:sitexam2 -0.471 0.667 0.333 0.333 -0.707 -0.590 -0.471
## genderm:sitneutr -0.574 0.812 0.406 0.406 -0.590 -0.707 -0.574
## grouphi:sitexam2 -0.471 0.333 0.667 0.333 -0.707 -0.590 -0.471
## grouplo:sitexam2 -0.471 0.333 0.333 0.667 -0.707 -0.590 -0.236
## grouphi:sitneutr -0.574 0.406 0.812 0.406 -0.590 -0.707 -0.574
## grouplo:sitneutr -0.574 0.406 0.406 0.812 -0.590 -0.707 -0.287
## genderm:grouphi:sitexam2 0.333 -0.471 -0.471 -0.236 0.500 0.417 0.667
## genderm:grouplo:sitexam2 0.333 -0.471 -0.236 -0.471 0.500 0.417 0.333
## genderm:grouphi:sitneutr 0.406 -0.574 -0.574 -0.287 0.417 0.500 0.812
## genderm:grouplo:sitneutr 0.406 -0.574 -0.287 -0.574 0.417 0.500 0.406
## gndrm:grpl gndr:2 gndrm:s grph:2 grpl:2 grph:s grpl:s
## genderm
## grouphi
## grouplo
## sitexam2
## sitneutr
## genderm:grouphi
## genderm:grouplo
## genderm:sitexam2 -0.471
## genderm:sitneutr -0.574 0.834
## grouphi:sitexam2 -0.236 0.500 0.417
## grouplo:sitexam2 -0.471 0.500 0.417 0.500
## grouphi:sitneutr -0.287 0.417 0.500 0.834 0.417
## grouplo:sitneutr -0.574 0.417 0.500 0.417 0.834 0.500
## genderm:grouphi:sitexam2 0.333 -0.707 -0.590 -0.707 -0.354 -0.590 -0.295
## genderm:grouplo:sitexam2 0.667 -0.707 -0.590 -0.354 -0.707 -0.295 -0.590
## genderm:grouphi:sitneutr 0.406 -0.590 -0.707 -0.590 -0.295 -0.707 -0.354
## genderm:grouplo:sitneutr 0.812 -0.590 -0.707 -0.295 -0.590 -0.354 -0.707
## gndrm:grph:2 gndrm:grpl:2 gndrm:grph:
## genderm
## grouphi
## grouplo
## sitexam2
## sitneutr
## genderm:grouphi
## genderm:grouplo
## genderm:sitexam2
## genderm:sitneutr
## grouphi:sitexam2
## grouplo:sitexam2
## grouphi:sitneutr
## grouplo:sitneutr
## genderm:grouphi:sitexam2
## genderm:grouplo:sitexam2 0.500
## genderm:grouphi:sitneutr 0.834 0.417
## genderm:grouplo:sitneutr 0.417 0.834 0.500
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -1.63923426 -0.39555906 -0.03025745 0.40201063 1.71650745
##
## Number of Observations: 540
## Number of Groups: 180# we might want to compare this with ezANOVA results (commented out)
# m.ez <- ezANOVA(data=dd.l, dv=.(anx), wid = .(subj), between=.(gender, group), within = .(sit), detailed=TRUE, type=2, return_aov = TRUE)
# summary(m.ez)
# m.ez
# we look at group descriptives
require(Rmisc)## Lade nötiges Paket: Rmisc## Lade nötiges Paket: lattice## Lade nötiges Paket: plyr## ------------------------------------------------------------------------------## You have loaded plyr after dplyr - this is likely to cause problems.
## If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
## library(plyr); library(dplyr)## ------------------------------------------------------------------------------##
## Attache Paket: 'plyr'## Die folgenden Objekte sind maskiert von 'package:dplyr':
##
## arrange, count, desc, failwith, id, mutate, rename, summarise,
## summarize## Das folgende Objekt ist maskiert 'package:purrr':
##
## compactsummarySE(data=dd.l, "anx", c("group", "gender", "sit"))## group gender sit N anx sd se ci
## 1 av f exam1 30 57.76667 12.053568 2.200670 4.500876
## 2 av f exam2 30 54.23333 12.229651 2.232819 4.566627
## 3 av f neutr 30 48.10000 12.644612 2.308580 4.721576
## 4 av m exam1 30 41.96667 9.308814 1.699549 3.475968
## 5 av m exam2 30 45.90000 10.892421 1.988675 4.067297
## 6 av m neutr 30 38.86667 10.294737 1.879553 3.844118
## 7 hi f exam1 30 70.80000 12.058478 2.201567 4.502710
## 8 hi f exam2 30 67.23333 12.316749 2.248720 4.599150
## 9 hi f neutr 30 61.16667 12.605646 2.301465 4.707025
## 10 hi m exam1 30 50.86667 9.346706 1.706467 3.490118
## 11 hi m exam2 30 55.03333 10.857965 1.982384 4.054431
## 12 hi m neutr 30 47.86667 10.173438 1.857407 3.798824
## 13 lo f exam1 30 42.83333 12.026073 2.195651 4.490610
## 14 lo f exam2 30 39.30000 12.304050 2.246402 4.594408
## 15 lo f neutr 30 33.16667 12.605646 2.301465 4.707025
## 16 lo m exam1 30 36.90000 9.393138 1.714944 3.507455
## 17 lo m exam2 30 40.90000 10.908238 1.991563 4.073203
## 18 lo m neutr 30 33.90000 10.185690 1.859644 3.803399detach("package:Rmisc", unload=TRUE)
# we can get at fixed effects
fixef(m.1)## (Intercept) genderm grouphi
## 5.776667e+01 -1.580000e+01 1.303333e+01
## grouplo sitexam2 sitneutr
## -1.493333e+01 -3.533333e+00 -9.666667e+00
## genderm:grouphi genderm:grouplo genderm:sitexam2
## -4.133333e+00 9.866667e+00 7.466667e+00
## genderm:sitneutr grouphi:sitexam2 grouplo:sitexam2
## 6.566667e+00 -3.333333e-02 -4.145480e-16
## grouphi:sitneutr grouplo:sitneutr genderm:grouphi:sitexam2
## 3.333333e-02 -1.570386e-14 2.666667e-01
## genderm:grouplo:sitexam2 genderm:grouphi:sitneutr genderm:grouplo:sitneutr
## 6.666667e-02 6.666667e-02 1.000000e-01# and at random effects, only the first few here because we have many subjects
head(ranef(m.1))## (Intercept)
## 1 -6.0265655
## 2 10.6214055
## 3 0.9655823
## 4 -27.3359684
## 5 -3.0299307
## 6 -9.0232003# we can predict individual values of level 1
head(predict(m.1, level=1))## 1 2 3 4 5 6
## 42.07343 58.72141 49.06558 20.76403 45.07007 39.07680dd.l$pred.vals <- predict(m.1, level=1) # Vorhersage auf Level 1
library(ICC)
ICCest(x = subj, y = anx, data = dd.l)## $ICC
## [1] 0.9223226
##
## $LowerCI
## [1] 0.9018286
##
## $UpperCI
## [1] 0.9393676
##
## $N
## [1] 180
##
## $k
## [1] 3
##
## $varw
## [1] 18.62778
##
## $vara
## [1] 221.1817Multilevel Ansatz mit lmer4::lmer() bzw.
library(lmerTest)
Wie vorgestellt in Field (2012, p. 573) Repeated-Measures Designs(GLM 4)
Benutzt wird lmer() aus dem package(lmer4).
Wegen der Diskussion um die korrekten P-Werte bei lmer4 Designs haben
die Autoren von library(lmer4) beschlossen, keine
p-Werte mehr auszugeben. Unter bestimmten Umständen (classical design,
balanced, nested, etc.) ist die Berechnung von p-Werten weniger
strittig. Wir können sie erhalten, wenn wir
library(lmerTest) benutzen, die ihrerseits auf
library(lmer4) aufbaut.
- Ein dreistufiger Between-Subjects-Faktor (Gruppe).
- Ein zweistufiger Between-Subjects-Faktor (Geschlecht).
- Ein dreiftufiger Within-Subjects-Faktor (Situation).
# # get data
# #require(reshape2)
# #require(psych)
# # require(lme4) # is loaded and used by lmerTest
# # require(lmerTest)
# dd <- read.delim("http://md.psych.bio.uni-goettingen.de/mv/data/virt/v_stait_exam_wide.txt")
# # we rename situation neutr to become reference group. This might be closer to our research design.
# # releveling has no impact in lme summaries.
# dd$aneutr <- dd$neutr
# dd$neutr <- NULL
#
# # convert it to long format
# # dd.l <- melt(dd, id=c('subj', 'gender', 'group'), variable.name="sit", value.name="anx", measured=c('neutr', 'exam1', 'exam2'))
#
# # convert it to long format using tidyr
# require(tidyr)
# dd.l <- dd %>% tidyr::gather(key=sit, value=anx, aneutr, exam1, exam2)
# dd.l$subj <- factor(dd.l$subj)
# dd.l$sit <- factor(dd.l$sit)
# head(dd.l)
dd <- read.delim("http://md.psych.bio.uni-goettingen.de/mv/data/virt/v_stait_exam_wide.txt")
# convert it to long format using tidyr
require(tidyr)
dd.l <- dd %>% tidyr::gather(key=sit, value=anx, neutr, exam1, exam2)
dd.l$subj <- factor(dd.l$subj)
dd.l$sit <- factor(dd.l$sit)
head(dd.l)## subj gender group sit anx
## 1 1 f av neutr 42
## 2 2 f av neutr 59
## 3 3 f av neutr 49
## 4 4 f av neutr 20
## 5 5 f av neutr 45
## 6 6 f av neutr 39#require(lme4)
require(lmerTest)## Lade nötiges Paket: lmerTest## Lade nötiges Paket: lme4## Lade nötiges Paket: Matrix##
## Attache Paket: 'Matrix'## Die folgenden Objekte sind maskiert von 'package:tidyr':
##
## expand, pack, unpack##
## Attache Paket: 'lme4'## Das folgende Objekt ist maskiert 'package:nlme':
##
## lmList##
## Attache Paket: 'lmerTest'## Das folgende Objekt ist maskiert 'package:lme4':
##
## lmer## Das folgende Objekt ist maskiert 'package:stats':
##
## step# m.f <- lme4::lmer(anx ~ group + gender + (1 | subj) + (1 | sit), data=dd.l)
# m.f
# summary(m.f)
m.t <- lmerTest::lmer(anx ~ group + gender + (1 | subj) + (1 | sit), data=dd.l)
m.t## Linear mixed model fit by REML ['lmerModLmerTest']
## Formula: anx ~ group + gender + (1 | subj) + (1 | sit)
## Data: dd.l
## REML criterion at convergence: 3167.463
## Random effects:
## Groups Name Std.Dev.
## subj (Intercept) 11.536
## sit (Intercept) 3.732
## Residual 2.168
## Number of obs: 540, groups: subj, 180; sit, 3
## Fixed Effects:
## (Intercept) grouphi grouplo genderm
## 52.383 11.022 -9.972 -9.156summary(m.t)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: anx ~ group + gender + (1 | subj) + (1 | sit)
## Data: dd.l
##
## REML criterion at convergence: 3167.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.89935 -0.68793 -0.08429 0.65393 1.67835
##
## Random effects:
## Groups Name Variance Std.Dev.
## subj (Intercept) 133.07 11.536
## sit (Intercept) 13.93 3.732
## Residual 4.70 2.168
## Number of obs: 540, groups: subj, 180; sit, 3
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 52.383 2.763 5.363 18.958 4.03e-06 ***
## grouphi 11.022 2.118 176.000 5.203 5.42e-07 ***
## grouplo -9.972 2.118 176.000 -4.707 5.07e-06 ***
## genderm -9.156 1.730 176.000 -5.293 3.55e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) grouph groupl
## grouphi -0.383
## grouplo -0.383 0.500
## genderm -0.313 0.000 0.000Grafik mit ggplot2
Hier ein expliziter Weg.
library(ez) erzeugt mit ezPlot() ebenfalls
ggplot-Grafiken.
require(reshape2)
dd <- read.delim("http://md.psych.bio.uni-goettingen.de/mv/data/virt/v_stait_exam_wide.txt")
# convert it to long format
dd.l <- melt(dd, id=c('subj', 'gender', 'group'), variable.name="sit", value.name="anx", measured=c('neutr', 'exam1', 'exam2'))
require(ggplot2)
ggplot(dd.l, aes(x=sit, y=anx)) +
# define colors
scale_fill_manual(values=c("#FFFFFF", "#CCCCCC", '#AAAAAA')) +
# create white background etc
#theme_bw() +
geom_point() +
# add verbose axis labels
labs(x = "Anxiety-Group", y = "STAI-Trait")
ggplot(dd.l, aes(x=sit, y=anx, fill=sit)) +
# define colors
scale_fill_manual(values=c("#FFFFFF", "#CCCCCC", '#AAAAAA')) +
# create white background etc
#theme_bw() +
stat_summary(fun.y=mean, geom="bar", position=position_dodge(width=0.7), colour='black') +
# add error information
stat_summary(fun.data = mean_cl_normal, geom = "errorbar", position=position_dodge(width=0.7), width=0.2) +
# put big points on top of bars
stat_summary(fun.data = mean_cl_normal, geom = "point", shape=21, size=5, position=position_dodge(width=0.7)) +
# add verbose axis labels
labs(x = "Anxiety-Group", y = "STAI-Trait")## Warning: `fun.y` is deprecated. Use `fun` instead.
lmer formulae und Erklärungen
What's the difference between (~1 +....) and (1 | ...) and (0 | ...) etc.?Say you have variable V1 predicted by categorical variable V2, which is treated as a random effect, and continuous variable V3, which is treated as a linear fixed effect. Using lmer syntax, simplest model (M1) is:
V1 ~ (1|V2) + V3
This model will estimate:
P1: A global intercept
P2: Random effect intercepts for V2 (i.e. for each level of V2, that level’s intercept’s deviation from the global intercept)
P3: A single global estimate for the effect (slope) of V3
The next most complex model (M2) is:
V1 ~ (1|V2) + V3 + (0+V3|V2)
This model estimates all the parameters from M1, but will additionally estimate:
P4: The effect of V3 within each level of V2 (more specifically, the degree to which the V3 effect within a given level deviates from the global effect of V3), while enforcing a zero correlation between the intercept deviations and V3 effect deviations across levels of V2.
This latter restriction is relaxed in a final most complex model (M3):
V1 ~ (1+V3|V2) + V3
In which all parameters from M2 are estimated while allowing correlation between the intercept deviations and V3 effect deviations within levels of V2. Thus, in M3, an additional parameter is estimated:
P5: The correlation between intercept deviations and V3 deviations across levels of V2
Usually model pairs like M2 and M3 are computed then compared to evaluate the evidence for correlations between fixed effects (including the global intercept).
Now consider adding another fixed effect predictor, V4. The model:
V1 ~ (1+V3V4|V2) + V3V4
would estimate:
P1: A global intercept
P2: A single global estimate for the effect of V3
P3: A single global estimate for the effect of V4
P4: A single global estimate for the interaction between V3 and V4
P5: Deviations of the intercept from P1 in each level of V2
P6: Deviations of the V3 effect from P2 in each level of V2
P7: Deviations of the V4 effect from P3 in each level of V2
P8: Deviations of the V3-by-V4 interaction from P4 in each level of V2
P9 Correlation between P5 and P6 across levels of V2
P10 Correlation between P5 and P7 across levels of V2
P11 Correlation between P5 and P8 across levels of V2
P12 Correlation between P6 and P7 across levels of V2
P13 Correlation between P6 and P8 across levels of V2
P14 Correlation between P7 and P8 across levels of V2
Phew, That’s a lot of parameters! And I didn’t even bother to list the variance parameters estimated by the model. What’s more, if you have a categorical variable with more than 2 levels that you want to model as a fixed effect, instead of a single effect for that variable you will always be estimating k-1 effects (where k is the number of levels), thereby exploding the number of parameters to be estimated by the model even further.
Ein paar unfertige Versuche
pitch and attractiveness
[http://seriousstats.wordpress.com/tag/repeated-measures-anova/]
# get data
dd <- read.csv('http://www2.ntupsychology.net/seriousstats/pitch.csv')
head(dd)## Constant Participant Face Context attract pitch base
## 1 1 1 1 1 2 99.88 97.79
## 2 1 1 2 1 1 108.50 100.40
## 3 1 1 3 1 3 99.06 101.00
## 4 1 1 4 1 4 97.02 100.28
## 5 1 1 5 1 2 102.04 97.79
## 6 1 1 6 1 3 98.38 101.00library(lme4)
pitch.me <- lmer(pitch ~ base + attract + (1|Face) + (1|Participant), data=dd)
pitch.me## Linear mixed model fit by REML ['lmerModLmerTest']
## Formula: pitch ~ base + attract + (1 | Face) + (1 | Participant)
## Data: dd
## REML criterion at convergence: 14121.79
## Random effects:
## Groups Name Std.Dev.
## Face (Intercept) 0.6665
## Participant (Intercept) 13.1046
## Residual 9.1909
## Number of obs: 1920, groups: Face, 32; Participant, 30
## Fixed Effects:
## (Intercept) base attract
## 90.0333 0.2054 0.4591Example: Politeness and pitch
http://www.bodowinter.com/tutorial/bw_LME_tutorial.pdf
The difference in politeness level is represented in the column called “attitude”. In that column, “pol” stands for polite and “inf” for informal. Sex is represented as “F” and “M” in the column “gender”. The dependent measure is “frequency”, which is the voice pitch measured in Hertz (Hz). To remind you, higher values mean higher pitch.
The interesting random effects for us are in the column “subject” and “scenario”, the latter being the name of the item column (remember the different scenarios like “asking for a favor”?).
# get data
# dd = read.csv("http://www.bodowinter.com/tutorial/politeness_data.csv")
dd <- readr::read_csv("http://md.psych.bio.uni-goettingen.de/mv/data/div/politeness_data.csv")
dd$subject <- factor(dd$subject)
dd$gender <- factor(dd$gender)
dd$attitude <- factor(dd$attitude)
require(psych)
psych::describe(dd)
# missing data?
which(is.na(dd)==T)
# a first boxplot
boxplot(frequency ~ attitude*gender,
col=c("white","lightgray"),dd)
# pitch.dat <- read.csv('http://www2.ntupsychology.net/seriousstats/pitch.csv')
# head(pitch.dat)
# library(lme4)
# pitch.me <- lmer(pitch ~ base + attract + (1|Face) + (1|Participant), data=pitch.dat)
# pitch.melmer() muss random effect haben, sonst läuft es nicht.
Modell mit Vp-Faktor und Faktor scenario:
# dd = read.csv("http://www.bodowinter.com/tutorial/politeness_data.csv")
dd <- readr::read_csv("http://md.psych.bio.uni-goettingen.de/mv/data/div/politeness_data.csv")
dd$subject <- factor(dd$subject)
dd$gender <- factor(dd$gender)
dd$attitude <- factor(dd$attitude)
# fit model
politeness.model = lmer(frequency ~ attitude + (1|subject) + (1|scenario), data=dd)
# look at the model
politeness.modelNote that the lmer() function (just like the lm() function in tutorial 1) took whatever comes first in the alphabet to be the reference level. “inf” comes before “pol”, so the slope represents the change from “inf” to “pol”.
Aufnahme von Geschlecht ins Modell. Geschlecht ist ein ‘fixed effekt’. Es gibt zwei Stufen (Ausprägungen) und wir erwarten im Mittel eine höhere Stimme bei Frauen.
# dd = read.csv("http://www.bodowinter.com/tutorial/politeness_data.csv")
dd <- readr::read_csv("http://md.psych.bio.uni-goettingen.de/mv/data/div/politeness_data.csv")
dd$subject <- factor(dd$subject)
dd$gender <- factor(dd$gender)
dd$attitude <- factor(dd$attitude)
# fit model
politeness.model = lmer(frequency ~ attitude + gender + (1|subject) + (1|scenario), data=dd)
# look at the model
politeness.modelGeschlecht bindet einen Teil der Variation, die vorher im subject-Faktor enthalten war.
Signifikanztests via Modellvergleich. Wegen der Diskussion um p-Werte in Mixed-Models hier mit Maximum-Likelihood Schätzung (REML-Flag).
# dd = read.csv("http://www.bodowinter.com/tutorial/politeness_data.csv")
dd <- readr::read_csv("http://md.psych.bio.uni-goettingen.de/mv/data/div/politeness_data.csv")
dd$subject <- factor(dd$subject)
dd$gender <- factor(dd$gender)
dd$attitude <- factor(dd$attitude)
# fit model
politeness.null = lmer(frequency ~ gender + (1|subject) + (1|scenario), data=dd, REML=FALSE)
politeness.model = lmer(frequency ~ attitude + gender + (1|subject) + (1|scenario), data=dd, REML=FALSE)
# perform likelihood ratio test of the difference
anova(politeness.null,politeness.model)
# random intercept model, so coefficients include subject and situation parameters
coef(politeness.model)random slope model (interaction)
# dd = read.csv("http://www.bodowinter.com/tutorial/politeness_data.csv")
# dd = read.csv("http://www.bodowinter.com/uploads/1/2/9/3/129362560/politeness_data.csv")
dd <- readr::read_csv("http://md.psych.bio.uni-goettingen.de/mv/data/div/politeness_data.csv")
dd$subject <- factor(dd$subject)
dd$gender <- factor(dd$gender)
dd$attitude <- factor(dd$attitude)
# fit model
politeness.model = lmer(frequency ~ attitude + gender + (1+attitude|subject) + (1+attitude|scenario), data=dd, REML=FALSE)
coef(politeness.model)
# test the effect of attitude, interaction included
politeness.null = lmer(frequency ~ gender + (1+attitude|subject) + (1+attitude|scenario), data=dd, REML=FALSE)
anova(politeness.null,politeness.model)Note that the null model needs to have the same random effects structure. So, if your full model is a random slope model, your null model also needs to be a random slope model.
influencial cases dfbeta() funktioniert nicht bei mixed models library(influence.ME) ansehen
require(lmerTest)
# politeness = read.csv("http://www.bodowinter.com/tutorial/politeness_data.csv")
dd <- readr::read_csv("http://md.psych.bio.uni-goettingen.de/mv/data/div/politeness_data.csv")
dd$subject <- factor(dd$subject)
dd$gender <- factor(dd$gender)
dd$attitude <- factor(dd$attitude)
# self written
#{} check somenumber
somenumber <- 5
all.res=numeric(nrow(politeness))
for(i in 1:nrow(politeness)){
#myfullmodel=lmer(frequency~attitude+(1+attitude|subject), data=politeness,POP[-i,])
myfullmodel=lmer(frequency~attitude+(1+attitude|subject), data=dd,REML=FALSE)
all.res[i]=fixef(myfullmodel)[somenumber]
}aov und formulae in reinen Messwiederholungsdesigns
aov()
Two factor design with repeated measures on both factors:
DV ~ IV1 * IV2 + Error(subject/(IV1*IV2))
Three factor mixed design with repeated measures on IV2 and IV3:
DV ~ IV1 * IV2 * IV3 + Error(subject/(IV2*IV3))
Links und Referenzen
http://dwoll.de/rexrepos/posts/anovaMixed.html
Wollschläger http://www.uni-kiel.de/psychologie/rexrepos/posts/anovaRBFpq.html
Bodo Winter: Mixed Models Tutorial and data seem to have moved. Some
links: [http://www.bodowinter.com/tutorial/bw_LME_tutorial.pdf]
[http://www.bodowinter.com/uploads/1/2/9/3/129362560/bw_lme_tutorial1.pdf]
Data
dd = read.csv("http://www.bodowinter.com/tutorial/politeness_data.csv")
dd = read.csv("http://www.bodowinter.com/uploads/1/2/9/3/129362560/politeness_data.csv")
Hao Zhang [http://www.stat.purdue.edu/~zhanghao/MAS/handout/anova%20with%20mixed%20effects.pdf]
Motor Control Lab http://gribblelab.wordpress.com/2009/03/09/repeated-measures-anova-using-r/
Missing values Mixed-Models-for-Repeated-Measures
Ein Video zum Unterschied fixed effects models vs. random effects models here
Zur ‘Correlation of Fixed Effects’ in lmer Erklärungsansatz Formel
Zur Diskussion p-Values und Vorgehen (statexchange Beitrag)[http://stats.stackexchange.com/questions/63464/is-this-an-acceptable-way-to-analyse-mixed-effect-models-with-lme4-in-r] {} Field (2012, p. 245) Kapitel 7: Regression
gar
ggplot(dd.r, aes(x=sit, y=anx)) +
# define colors
scale_fill_manual(values=c("#FFFFFF", "#CCCCCC", '#AAAAAA')) +
# create white background etc
#theme_bw() +
geom_point() +
# add verbose axis labels
labs(x = "Anxiety-Group", y = "STAI-Trait")
ggplot(dd.r, aes(x=sit, y=anx, fill=sit)) +
# define colors
scale_fill_manual(values=c("#FFFFFF", "#CCCCCC", '#AAAAAA')) +
# create white background etc
#theme_bw() +
stat_summary(fun.y=mean, geom="bar", position=position_dodge(width=0.7), colour='black') +
# add error information
stat_summary(fun.data = mean_cl_normal, geom = "errorbar", position=position_dodge(width=0.7), width=0.2) +
# put big points on top of bars
stat_summary(fun.data = mean_cl_normal, geom = "point", shape=21, size=5, position=position_dodge(width=0.7)) +
# add verbose axis labels
labs(x = "Anxiety-Group", y = "STAI-Trait")Multi-Level-Ansatz: Beispiel STAI-State verschiedene Gruppen in verschiedenen Situationen
Ein dreistufiger Between-Subjects-Faktor (Gruppe). Ein dreiftufiger Within-Faktor (Situation).
dd <- read.delim("http://md.psych.bio.uni-goettingen.de/mv/data/virt/v_stait_exam_wide.txt")
require("reshape2")
require(psych)## Lade nötiges Paket: psych##
## Attache Paket: 'psych'## Die folgenden Objekte sind maskiert von 'package:ggplot2':
##
## %+%, alpha# require(lme4) # is loaded and used by lmerTest
require(lmerTest)
head(dd)## subj gender group neutr exam1 exam2
## 1 1 f av 42 52 48
## 2 2 f av 59 68 65
## 3 3 f av 49 59 55
## 4 4 f av 20 31 27
## 5 5 f av 45 55 51
## 6 6 f av 39 49 45psych::describe(dd)## vars n mean sd median trimmed mad min max range skew kurtosis
## subj 1 180 90.50 52.11 90.5 90.50 66.72 1 180 179 0.00 -1.22
## gender* 2 180 1.50 0.50 1.5 1.50 0.74 1 2 1 0.00 -2.01
## group* 3 180 2.00 0.82 2.0 2.00 1.48 1 3 2 0.00 -1.52
## neutr 4 180 43.84 14.97 42.0 43.44 15.57 5 89 84 0.25 -0.06
## exam1 5 180 50.19 15.62 48.5 49.47 15.57 16 98 82 0.44 -0.13
## exam2 6 180 50.43 14.96 49.0 50.08 16.31 12 95 83 0.21 -0.13
## se
## subj 3.88
## gender* 0.04
## group* 0.06
## neutr 1.12
## exam1 1.16
## exam2 1.11# convert it to long format
dd.r <- reshape2::melt(dd, id=c('subj', 'gender', 'group'), variable.name='sit', value.name='anx', measured=c('neutr', 'exam1', 'exam2'))
# random intercept model (global intercept for subjects, not specific for fixed factor level)
m.1 <- lmer(anx ~ group + sit + (1|subj), data=dd.r)
m.1## Linear mixed model fit by REML ['lmerModLmerTest']
## Formula: anx ~ group + sit + (1 | subj)
## Data: dd.r
## REML criterion at convergence: 3184.494
## Random effects:
## Groups Name Std.Dev.
## subj (Intercept) 12.395
## Residual 2.168
## Number of obs: 540, groups: subj, 180
## Fixed Effects:
## (Intercept) grouphi grouplo sitexam1 sitexam2
## 43.494 11.022 -9.972 6.344 6.589summary(m.1)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: anx ~ group + sit + (1 | subj)
## Data: dd.r
##
## REML criterion at convergence: 3184.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.9374 -0.6756 -0.0809 0.6292 1.6401
##
## Random effects:
## Groups Name Variance Std.Dev.
## subj (Intercept) 153.6 12.395
## Residual 4.7 2.168
## Number of obs: 540, groups: subj, 180
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 43.4944 1.6137 179.3863 26.954 < 2e-16 ***
## grouphi 11.0222 2.2744 177.0000 4.846 2.74e-06 ***
## grouplo -9.9722 2.2744 177.0000 -4.384 1.99e-05 ***
## sitexam1 6.3444 0.2285 358.0000 27.764 < 2e-16 ***
## sitexam2 6.5889 0.2285 358.0000 28.834 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) grouph groupl sitxm1
## grouphi -0.705
## grouplo -0.705 0.500
## sitexam1 -0.071 0.000 0.000
## sitexam2 -0.071 0.000 0.000 0.500anova(m.1)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## group 400.8 200.38 2 177 42.638 7.684e-16 ***
## sit 5023.5 2511.76 2 358 534.451 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# descriptives - the Rcmdr way
tapply(dd.r$anx, list(situation=dd.r$sit), mean, na.rm=TRUE) # means## situation
## neutr exam1 exam2
## 43.84444 50.18889 50.43333tapply(dd.r$anx, list(situation=dd.r$sit), function(x) sum(!is.na(x))) # counts## situation
## neutr exam1 exam2
## 180 180 180# the psych package way
by(dd.r$anx, list(dd.r$sit), psych::describe)## : neutr
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 180 43.84 14.97 42 43.44 15.57 5 89 84 0.25 -0.06 1.12
## ------------------------------------------------------------
## : exam1
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 180 50.19 15.62 48.5 49.47 15.57 16 98 82 0.44 -0.13 1.16
## ------------------------------------------------------------
## : exam2
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 180 50.43 14.96 49 50.08 16.31 12 95 83 0.21 -0.13 1.11# visualization
require(ggplot2)
ggplot(dd, x=group, y=neutr, aes(group, neutr, fill = group)) +
stat_summary(fun.y=mean, geom="bar", position=position_dodge()) +
stat_summary(fun.data = mean_cl_normal, geom = "errorbar", position=position_dodge(width=0.9), width=0.2) ## Warning: `fun.y` is deprecated. Use `fun` instead.
ggplot(dd.r, aes(x=sit, y=anx)) +
# define colors
scale_fill_manual(values=c("#FFFFFF", "#CCCCCC", '#AAAAAA')) +
# create white background etc
#theme_bw() +
geom_point() +
# add verbose axis labels
labs(x = "Anxiety-Group", y = "STAI-Trait")
Check
Check, whether you have understood the following statements
In between factors we have different observations in every sub-category.
In within factors we have the same observations over the sub-cateogories of the factor. This implies dependencies. Several measures of the same observation are usually more similar than several measures of different observations.
In Anova designs the terms between-subjects-design and within-subjects-design can be found. Subjects are often our observations, but not always. Therefore we prefer the term within or between factor.
Mixed models mix between factors and within factors.
When we have continuous explanatory variables in our design, this is often called ANCOVA (Analysis of Covariance)
We recommend using multi level techniques (packages) to analyse mixed designs usually using the
package(nlme)but you might also consider usingpackage(lme4), which uses a different syntaxWhen we apply MLM to repeated measures designs we consider the observations to be a level above level 1. Measuring them repeatedly implies a hierarchy in our data. We allow random effects for our observations, f. e. random intercepts due to individual characteristics of them.
Random effects and their coefficients apply only to the adequate higher level subgroup whereas fixed effects are fix for all observations. So random coefficients are specific additional effects (corrections) over fixed effects in members of sub-levels.
It’s always a good idea to take a look at the coefficients of the adapted MLM to understand, what is happening. This also helps for visualizations of the effects.
All the aspects of categorial explanatory variables apply also (dummy-coding, contrasts, unbalanced designs, type of SS, etc)
Prediction includes random effects also. If we do not have the same observations again, this makes no sense, as the characteristics of a specific observation are used. But we can limit the prediction to levels used.
level=0uses fixed effects only.Example
# convert it to long format using tidyr
library(tidyverse)
dd.w <- read_tsv("http://md.psych.bio.uni-goettingen.de/mv/data/virt/v_stait_exam_wide.txt")## Rows: 180 Columns: 6
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (2): gender, group
## dbl (4): subj, neutr, exam1, exam2
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.dd <- dd.w %>% tidyr::gather(key=sit, value=anx, neutr, exam1, exam2)
dd$subj.f <- factor(dd$subj)
dd$group.f <- factor(dd$group, levels=c("av", "lo", "hi")) # tendency to show anxiety, actual anxiety in variable dd$anx
dd$sit.f <- factor(dd$sit, levels=c("neutr", "exam1", "exam2")) # situation
head(dd)## # A tibble: 6 × 8
## subj gender group sit anx subj.f group.f sit.f
## <dbl> <chr> <chr> <chr> <dbl> <fct> <fct> <fct>
## 1 1 f av neutr 42 1 av neutr
## 2 2 f av neutr 59 2 av neutr
## 3 3 f av neutr 49 3 av neutr
## 4 4 f av neutr 20 4 av neutr
## 5 5 f av neutr 45 5 av neutr
## 6 6 f av neutr 39 6 av neutrlibrary(nlme)
m.00 <- gls(anx ~ 1, data=dd, method="ML")
m.0 <- lme(anx ~ 1, random = ~1|subj, data=dd, method="ML")
m.g <- lme(anx ~ group.f, random = ~1|subj, data=dd, method="ML")
m.s <- lme(anx ~ sit.f, random = ~1|subj, data=dd, method="ML")
m.g.s <- lme(anx ~ group.f + sit.f, random = ~1|subj, data=dd, method="ML")
m.gs <- lme(anx ~ group.f * sit.f, random = ~1|subj, data=dd, method="ML")
summary(m.00)## Generalized least squares fit by maximum likelihood
## Model: anx ~ 1
## Data: dd
## AIC BIC logLik
## 4492.718 4501.301 -2244.359
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 48.15556 0.6652607 72.386 0
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -2.79415494 -0.72227898 -0.07481774 0.63738962 3.22723457
##
## Residual standard error: 15.44494
## Degrees of freedom: 540 total; 539 residualggplot(data=dd) + geom_point(aes(x=sit.f, y=anx, color=group.f), size=2) +
geom_point(aes(x=sit.f, y=predict(m.00)), color="grey") +
geom_point(aes(x=sit.f, y=m.00$coefficients[[1]]), size = 4, color = "grey") +
geom_line(aes(x=sit.f, y=m.00$coefficients[[1]]), color = "grey", linetype = "dashed") +
facet_grid(.~ group.f) 
summary(m.0)## Linear mixed-effects model fit by maximum likelihood
## Data: dd
## AIC BIC logLik
## 3764.877 3777.752 -1879.439
##
## Random effects:
## Formula: ~1 | subj
## (Intercept) Residual
## StdDev: 14.82965 4.315991
##
## Fixed effects: anx ~ 1
## Value Std.Error DF t-value p-value
## (Intercept) 48.15556 1.121871 360 42.92431 0
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -1.62656914 -0.80843768 0.06334385 0.79978749 1.48585661
##
## Number of Observations: 540
## Number of Groups: 180fixef(m.0)## (Intercept)
## 48.15556head(ranef(m.0))## (Intercept)
## 1 -0.7996447
## 2 15.4093704
## 3 6.0081416
## 4 -21.5471842
## 5 2.1179780
## 6 -3.7172675ggplot(data=dd) + geom_point(aes(x=sit.f, y=anx, color=group.f), size=2) +
geom_point(aes(x=sit.f, y=predict(m.0)), color="grey") +
geom_point(aes(x=sit.f, y=m.0$coefficients[[1]]), size = 4, color = "grey") +
geom_line(aes(x=sit.f, y=m.0$coefficients[[1]]), color = "grey", linetype = "dashed") +
facet_grid(.~ group.f) 
fixef(m.g)## (Intercept) group.flo group.fhi
## 47.805556 -9.972222 11.022222head(ranef(m.g))## (Intercept)
## 1 -0.4530083
## 2 15.5355190
## 3 6.2621732
## 4 -20.9183232
## 5 2.4249266
## 6 -3.3309432summary(m.g)## Linear mixed-effects model fit by maximum likelihood
## Data: dd
## AIC BIC logLik
## 3698.093 3719.551 -1844.047
##
## Random effects:
## Formula: ~1 | subj
## (Intercept) Residual
## StdDev: 12.09942 4.315991
##
## Fixed effects: anx ~ group.f
## Value Std.Error DF t-value p-value
## (Intercept) 47.80556 1.599259 360 29.892309 0
## group.flo -9.97222 2.261694 177 -4.409182 0
## group.fhi 11.02222 2.261694 177 4.873436 0
## Correlation:
## (Intr) grp.fl
## group.flo -0.707
## group.fhi -0.707 0.500
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -1.64314627 -0.80505553 0.05494354 0.80974947 1.46855051
##
## Number of Observations: 540
## Number of Groups: 180fixef(m.s)## (Intercept) sit.fexam1 sit.fexam2
## 43.844444 6.344444 6.588889head(ranef(m.s))## (Intercept)
## 1 -0.8165577
## 2 15.7352874
## 3 6.1352172
## 4 -22.0029194
## 5 2.1627744
## 6 -3.7958898(pl <- ggplot(data=dd) + geom_point(aes(x=sit.f, y=anx, color=group.f), size=2) +
geom_point(aes(x=sit.f, y=predict(m.s)), color="grey") +
stat_summary(fun = mean, aes(x=sit.f, y=predict(m.s)), geom="point", size = 4, color = "grey") +
stat_summary(fun = mean, aes(x=sit.f, y=predict(m.s), group=subj), geom="line",color = "grey", linetype = "dashed") +
facet_grid(.~ group.f)
)
summary(m.s)## Linear mixed-effects model fit by maximum likelihood
## Data: dd
## AIC BIC logLik
## 3271.095 3292.553 -1630.547
##
## Random effects:
## Formula: ~1 | subj
## (Intercept) Residual
## StdDev: 14.98565 2.161846
##
## Fixed effects: anx ~ sit.f
## Value Std.Error DF t-value p-value
## (Intercept) 43.84444 1.1316754 358 38.74295 0
## sit.fexam1 6.34444 0.2285143 358 27.76389 0
## sit.fexam2 6.58889 0.2285143 358 28.83360 0
## Correlation:
## (Intr) st.fx1
## sit.fexam1 -0.101
## sit.fexam2 -0.101 0.500
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -1.93326906 -0.66793218 -0.06934972 0.63027439 1.65419914
##
## Number of Observations: 540
## Number of Groups: 180fixef(m.gs)## (Intercept) group.flo group.fhi
## 43.48333333 -9.95000000 11.03333333
## sit.fexam1 sit.fexam2 group.flo:sit.fexam1
## 6.38333333 6.58333333 -0.05000000
## group.fhi:sit.fexam1 group.flo:sit.fexam2 group.fhi:sit.fexam2
## -0.06666667 -0.01666667 0.03333333head(ranef(m.gs))## (Intercept)
## 1 -0.4674021
## 2 16.0291421
## 3 6.4611465
## 4 -21.5829786
## 5 2.5019759
## 6 -3.4367800ggplot(data=dd) + geom_point(aes(x=sit.f, y=anx, group=group.f, color=group.f), size=2) +
geom_point(aes(x=sit.f, y=predict(m.gs))) +
stat_summary(fun = mean, aes(x=sit.f, y=predict(m.gs)), geom="point", size = 4, color = "grey") +
# stat_summary(fun = mean, aes(x=sit.f, y=predict(m.gs)), geom="line", color = "grey", linetype = "dashed") +
stat_summary(fun = mean, aes(x=sit.f, y=predict(m.gs), group=subj), geom="line",color = "grey", linetype = "dashed") +
facet_grid(.~ group.f)
summary(m.gs)## Linear mixed-effects model fit by maximum likelihood
## Data: dd
## AIC BIC logLik
## 3212.273 3259.48 -1595.136
##
## Random effects:
## Formula: ~1 | subj
## (Intercept) Residual
## StdDev: 12.29014 2.161732
##
## Fixed effects: anx ~ group.f * sit.f
## Value Std.Error DF t-value p-value
## (Intercept) 43.48333 1.6246026 354 26.765520 0.0000
## group.flo -9.95000 2.2975350 177 -4.330728 0.0000
## group.fhi 11.03333 2.2975350 177 4.802248 0.0000
## sit.fexam1 6.38333 0.3980072 354 16.038237 0.0000
## sit.fexam2 6.58333 0.3980072 354 16.540740 0.0000
## group.flo:sit.fexam1 -0.05000 0.5628671 354 -0.088831 0.9293
## group.fhi:sit.fexam1 -0.06667 0.5628671 354 -0.118441 0.9058
## group.flo:sit.fexam2 -0.01667 0.5628671 354 -0.029610 0.9764
## group.fhi:sit.fexam2 0.03333 0.5628671 354 0.059221 0.9528
## Correlation:
## (Intr) grp.fl grp.fh st.fx1 st.fx2 grp.fl:.1 grp.fh:.1
## group.flo -0.707
## group.fhi -0.707 0.500
## sit.fexam1 -0.122 0.087 0.087
## sit.fexam2 -0.122 0.087 0.087 0.500
## group.flo:sit.fexam1 0.087 -0.122 -0.061 -0.707 -0.354
## group.fhi:sit.fexam1 0.087 -0.061 -0.122 -0.707 -0.354 0.500
## group.flo:sit.fexam2 0.087 -0.122 -0.061 -0.354 -0.707 0.500 0.250
## group.fhi:sit.fexam2 0.087 -0.061 -0.122 -0.354 -0.707 0.250 0.500
## grp.fl:.2
## group.flo
## group.fhi
## sit.fexam1
## sit.fexam2
## group.flo:sit.fexam1
## group.fhi:sit.fexam1
## group.flo:sit.fexam2
## group.fhi:sit.fexam2 0.500
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -1.92919862 -0.67498189 -0.07762239 0.63649945 1.63275958
##
## Number of Observations: 540
## Number of Groups: 180# our prediction is a combination of fixed and random effects
head(model.matrix(m.gs, data=dd))## (Intercept) group.flo group.fhi sit.fexam1 sit.fexam2 group.flo:sit.fexam1
## 1 1 0 0 0 0 0
## 2 1 0 0 0 0 0
## 3 1 0 0 0 0 0
## 4 1 0 0 0 0 0
## 5 1 0 0 0 0 0
## 6 1 0 0 0 0 0
## group.fhi:sit.fexam1 group.flo:sit.fexam2 group.fhi:sit.fexam2
## 1 0 0 0
## 2 0 0 0
## 3 0 0 0
## 4 0 0 0
## 5 0 0 0
## 6 0 0 0head(unique(model.matrix(m.gs, data=dd)))## (Intercept) group.flo group.fhi sit.fexam1 sit.fexam2 group.flo:sit.fexam1
## 1 1 0 0 0 0 0
## 61 1 1 0 0 0 0
## 121 1 0 1 0 0 0
## 181 1 0 0 1 0 0
## 241 1 1 0 1 0 1
## 301 1 0 1 1 0 0
## group.fhi:sit.fexam1 group.flo:sit.fexam2 group.fhi:sit.fexam2
## 1 0 0 0
## 61 0 0 0
## 121 0 0 0
## 181 0 0 0
## 241 0 0 0
## 301 1 0 0# as.numeric(rownames(head(unique(model.matrix(m.gs, data=dd)))))
head(m.gs$fitted)## fixed subj
## 1 43.48333 43.01593
## 2 43.48333 59.51248
## 3 43.48333 49.94448
## 4 43.48333 21.90035
## 5 43.48333 45.98531
## 6 43.48333 40.04655m.gs$fitted[as.numeric(rownames(head(unique(model.matrix(m.gs, data=dd))))),]## fixed subj
## 1 43.48333 43.01593
## 61 33.53333 28.41940
## 121 54.51667 56.00685
## 181 49.86667 49.39926
## 241 39.86667 34.75274
## 301 60.83333 62.32352head(predict(m.gs))## 1 2 3 4 5 6
## 43.01593 59.51248 49.94448 21.90035 45.98531 40.04655dd.n <- tibble(group.f = factor(c("av","av","av","lo","lo","lo","hi","hi","hi")),
sit.f = factor(c("neutr","exam1", "exam2","neutr","exam1", "exam2","neutr","exam1", "exam2")),
subj = c(180,179,177,181,182,183,184,185,186)
)
# subj is 1:180, predict needs all variables in newdata that occur in the model behind, random variables also
predict(m.gs, dd.n)## 180 179 177 181 182 183 184 185
## 37.71504 35.19024 51.55685 NA NA NA NA NA
## 186
## NA
## attr(,"label")
## [1] "Predicted values"# prediction only on base of fixed effects or on base of a specific max. level can be done using the level parameter
predict(m.gs, dd.n, level=0) # here only the fixed effects are used to do prediction## [1] 43.48333 49.86667 50.06667 33.53333 39.86667 40.10000 54.51667 60.83333
## [9] 61.13333
## attr(,"label")
## [1] "Predicted values"