library(lme4)
## Loading required package: Matrix
str(sleepstudy)
## 'data.frame': 180 obs. of 3 variables:
## $ Reaction: num 250 259 251 321 357 ...
## $ Days : num 0 1 2 3 4 5 6 7 8 9 ...
## $ Subject : Factor w/ 18 levels "308","309","310",..: 1 1 1 1 1 1 1 1 1 1 ...
xtabs(~ Subject + Days, sleepstudy)
## Days
## Subject 0 1 2 3 4 5 6 7 8 9
## 308 1 1 1 1 1 1 1 1 1 1
## 309 1 1 1 1 1 1 1 1 1 1
## 310 1 1 1 1 1 1 1 1 1 1
## 330 1 1 1 1 1 1 1 1 1 1
## 331 1 1 1 1 1 1 1 1 1 1
## 332 1 1 1 1 1 1 1 1 1 1
## 333 1 1 1 1 1 1 1 1 1 1
## 334 1 1 1 1 1 1 1 1 1 1
## 335 1 1 1 1 1 1 1 1 1 1
## 337 1 1 1 1 1 1 1 1 1 1
## 349 1 1 1 1 1 1 1 1 1 1
## 350 1 1 1 1 1 1 1 1 1 1
## 351 1 1 1 1 1 1 1 1 1 1
## 352 1 1 1 1 1 1 1 1 1 1
## 369 1 1 1 1 1 1 1 1 1 1
## 370 1 1 1 1 1 1 1 1 1 1
## 371 1 1 1 1 1 1 1 1 1 1
## 372 1 1 1 1 1 1 1 1 1 1
(model6 <- lmer(Reaction ~ 1 + Days + (1 + Days|Subject), sleepstudy, REML=FALSE))
## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: Reaction ~ 1 + Days + (1 + Days | Subject)
## Data: sleepstudy
## AIC BIC logLik deviance df.resid
## 1763.9393 1783.0971 -875.9697 1751.9393 174
## Random effects:
## Groups Name Std.Dev. Corr
## Subject (Intercept) 23.781
## Days 5.717 0.08
## Residual 25.592
## Number of obs: 180, groups: Subject, 18
## Fixed Effects:
## (Intercept) Days
## 251.41 10.47
head(ranef(model6)[["Subject"]])
## (Intercept) Days
## 308 2.815682 9.0755339
## 309 -40.048485 -8.6440673
## 310 -38.433152 -5.5133788
## 330 22.832296 -4.6587504
## 331 21.549990 -2.9445201
## 332 8.815586 -0.2352092
(model7 <- lmer(Reaction ~ 1 + Days + (1|Subject) + (0+Days|Subject), sleepstudy, REML=FALSE))
## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: Reaction ~ 1 + Days + (1 | Subject) + (0 + Days | Subject)
## Data: sleepstudy
## AIC BIC logLik deviance df.resid
## 1762.0033 1777.9680 -876.0016 1752.0033 175
## Random effects:
## Groups Name Std.Dev.
## Subject (Intercept) 24.171
## Subject.1 Days 5.799
## Residual 25.556
## Number of obs: 180, groups: Subject, 18
## Fixed Effects:
## (Intercept) Days
## 251.41 10.47
head(ranef(model7)[["Subject"]])
## (Intercept) Days
## 308 1.854656 9.2364346
## 309 -40.022307 -8.6174730
## 310 -38.723163 -5.4343800
## 330 23.903320 -4.8581940
## 331 22.396322 -3.1048405
## 332 9.052001 -0.2821598
anova(model7, model6)
## Data: sleepstudy
## Models:
## model7: Reaction ~ 1 + Days + (1 | Subject) + (0 + Days | Subject)
## model6: Reaction ~ 1 + Days + (1 + Days | Subject)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## model7 5 1762.0 1778.0 -876.00 1752.0
## model6 6 1763.9 1783.1 -875.97 1751.9 0.0639 1 0.8004
pr07 <- profile(model7)
confint(pr07)
## 2.5 % 97.5 %
## .sig01 15.258647 37.786472
## .sig02 3.964074 8.769159
## .sigma 22.880555 28.787598
## (Intercept) 237.572148 265.238062
## Days 7.334067 13.600505
str(rr1 <- ranef(model7))
## List of 1
## $ Subject:'data.frame': 18 obs. of 2 variables:
## ..$ (Intercept): num [1:18] 1.85 -40.02 -38.72 23.9 22.4 ...
## ..$ Days : num [1:18] 9.24 -8.62 -5.43 -4.86 -3.1 ...
## - attr(*, "class")= chr "ranef.mer"
plot(ranef(model7))
## $Subject
library(MEMSS)
##
## Attaching package: 'MEMSS'
## The following objects are masked from 'package:datasets':
##
## CO2, Orange, Theoph
str(Machines)
## 'data.frame': 54 obs. of 3 variables:
## $ Worker : Factor w/ 6 levels "1","2","3","4",..: 1 1 1 2 2 2 3 3 3 4 ...
## $ Machine: Factor w/ 3 levels "A","B","C": 1 1 1 1 1 1 1 1 1 1 ...
## $ score : num 52 52.8 53.1 51.8 52.8 53.1 60 60.2 58.4 51.1 ...
xtabs(~ Machine + Worker, Machines)
## Worker
## Machine 1 2 3 4 5 6
## A 3 3 3 3 3 3
## B 3 3 3 3 3 3
## C 3 3 3 3 3 3
model10 <- lmer(score ~ Machine + (1|Worker), Machines, REML=FALSE)
model11 <- lmer(score ~ Machine + (Machine|Worker), Machines, REML=FALSE)
model12 <- lmer(score ~ Machine + (1|Worker) + (1|Machine:Worker), Machines, REML=FALSE)
anova(model10, model11, model12)
## Data: Machines
## Models:
## model10: score ~ Machine + (1 | Worker)
## model12: score ~ Machine + (1 | Worker) + (1 | Machine:Worker)
## model11: score ~ Machine + (Machine | Worker)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## model10 5 303.70 313.65 -146.85 293.70
## model12 6 237.27 249.20 -112.64 225.27 68.4338 1 < 2e-16 ***
## model11 10 236.42 256.31 -108.21 216.42 8.8516 4 0.06492 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
str(ergoStool)
## 'data.frame': 36 obs. of 3 variables:
## $ effort : num 12 15 12 10 10 14 13 12 7 14 ...
## $ Type : Factor w/ 4 levels "T1","T2","T3",..: 1 2 3 4 1 2 3 4 1 2 ...
## $ Subject: Factor w/ 9 levels "A","B","C","D",..: 1 1 1 1 2 2 2 2 3 3 ...
model13 <- lmer(effort ~ 1 + (1|Subject) + (1|Type), ergoStool, REML=FALSE)
model13
## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: effort ~ 1 + (1 | Subject) + (1 | Type)
## Data: ergoStool
## AIC BIC logLik deviance df.resid
## 144.0224 150.3564 -68.0112 136.0224 32
## Random effects:
## Groups Name Std.Dev.
## Subject (Intercept) 1.305
## Type (Intercept) 1.505
## Residual 1.101
## Number of obs: 36, groups: Subject, 9; Type, 4
## Fixed Effects:
## (Intercept)
## 10.25
model13a <- lmer(effort ~ 1 + (1|Subject), ergoStool, REML=FALSE)
anova(model13a, model13)
## Data: ergoStool
## Models:
## model13a: effort ~ 1 + (1 | Subject)
## model13: effort ~ 1 + (1 | Subject) + (1 | Type)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## model13a 3 164.15 168.90 -79.075 158.15
## model13 4 144.02 150.36 -68.011 136.02 22.128 1 2.551e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(model14 <- lmer(effort ~ 1 + Type + (1|Subject), ergoStool, REML = 0))
## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: effort ~ 1 + Type + (1 | Subject)
## Data: ergoStool
## AIC BIC logLik deviance df.resid
## 134.1444 143.6456 -61.0722 122.1444 30
## Random effects:
## Groups Name Std.Dev.
## Subject (Intercept) 1.256
## Residual 1.037
## Number of obs: 36, groups: Subject, 9
## Fixed Effects:
## (Intercept) TypeT2 TypeT3 TypeT4
## 8.5556 3.8889 2.2222 0.6667
model14a <- profile(lmer(effort ~ 1 + Type + (1|Subject), within(ergoStool, Type <- relevel(Type, "T2")), REML=FALSE))
pr14b <- profile(lmer(effort ~ 1 + Type + (1|Subject), within(ergoStool, Type <- relevel(Type, "T3")), REML=FALSE))
summary(model15 <- aov(effort ~ Subject + Type, ergoStool))
## Df Sum Sq Mean Sq F value Pr(>F)
## Subject 8 66.50 8.312 6.866 0.000106 ***
## Type 3 81.19 27.065 22.356 3.93e-07 ***
## Residuals 24 29.06 1.211
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
coef(summary.lm(model15))
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.055556e+01 0.6352554 1.661624e+01 1.146367e-14
## SubjectB 6.005604e-15 0.7780258 7.719029e-15 1.000000e+00
## SubjectC -1.500000e+00 0.7780258 -1.927957e+00 6.577394e-02
## SubjectD -3.000000e+00 0.7780258 -3.855913e+00 7.577323e-04
## SubjectE -3.750000e+00 0.7780258 -4.819892e+00 6.565314e-05
## SubjectF -2.000000e+00 0.7780258 -2.570609e+00 1.678081e-02
## SubjectG -1.500000e+00 0.7780258 -1.927957e+00 6.577394e-02
## SubjectH -4.000000e+00 0.7780258 -5.141218e+00 2.907651e-05
## SubjectI -2.250000e+00 0.7780258 -2.891935e+00 8.010522e-03
## TypeT2 3.888889e+00 0.5186838 7.497610e+00 9.753420e-08
## TypeT3 2.222222e+00 0.5186838 4.284348e+00 2.562927e-04
## TypeT4 6.666667e-01 0.5186838 1.285304e+00 2.109512e-01
TukeyHSD(model15, which = "Type")
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = effort ~ Subject + Type, data = ergoStool)
##
## $Type
## diff lwr upr p adj
## T2-T1 3.8888889 2.4580431 5.3197347 0.0000006
## T3-T1 2.2222222 0.7913764 3.6530680 0.0013709
## T4-T1 0.6666667 -0.7641791 2.0975125 0.5807508
## T3-T2 -1.6666667 -3.0975125 -0.2358209 0.0182037
## T4-T2 -3.2222222 -4.6530680 -1.7913764 0.0000115
## T4-T3 -1.5555556 -2.9864013 -0.1247098 0.0295822
kappa(model15, exact = TRUE)
## [1] 10.76924