options(contrasts=c("contr.sum","contr.poly"))
load(url("http://pnb.mcmaster.ca/bennett/psy710/datasets/smile-2023.rda") )
sapply(smile.df,class)
## smile leniency subjID
## "factor" "numeric" "factor"
sapply(smile.wide,class)
## neutral contempt miserable real
## "numeric" "numeric" "numeric" "numeric"
class(smile.mat)
## [1] "matrix" "array"The command loads the long-format data frame smile.df,
the wide-format data frame smile.wide, and the wide-format
matrix smile.mat.
smile. Adjust the \(p\) value with the Geisser-Greenhouse or
Huynh-Feldt adjustment to correct for deviations from sphericity.# with aov_car:
library(afex)
smile.aov.01 <- aov_car(leniency~smile+Error(subjID/smile),data=smile.df)
summary(smile.aov.01)
##
## Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
##
## Sum Sq num Df Error SS den Df F value Pr(>F)
## (Intercept) 200000 1 5653 19 672.25 2.7e-16
## smile 1320 3 6934 57 3.62 0.018
##
##
## Mauchly Tests for Sphericity
##
## Test statistic p-value
## smile 0.634 0.152
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## smile 0.768 0.03
##
## HF eps Pr(>F[HF])
## smile 0.8799 0.02352lmer command in the lmerTest
package to evaluate the effect of smile with a mixed model.
Evaluate the effect with an \(F\) test
that assumes sphericity and with a chi-square test that does
not assume sphereicity.# with lmer:
library(lmerTest)
smile.lmer.01 <- lmer(leniency~smile+(1|subjID),data=smile.df)
print(anova(smile.lmer.01)) # assumes sphericity
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## smile 1320 440 3 57 3.62 0.018
# following test does not assume sphericity:
# drop fixed effect
smile.lmer.00 <- lmer(leniency~1+(1|subjID),data=smile.df)
print(anova(smile.lmer.00,smile.lmer.01) ) # compare models
## Data: smile.df
## Models:
## smile.lmer.00: leniency ~ 1 + (1 | subjID)
## smile.lmer.01: leniency ~ smile + (1 | subjID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## smile.lmer.00 3 641 649 -318 635
## smile.lmer.01 6 637 651 -312 625 10.5 3 0.015
# model with fixed effect fits significantly betterlme command in the nlme package to
evaluate the effect of smile with two mixed models: one
that assumes the variance-covariance matrix of the residuals is compound
symmetric (i.e., the correlations between conditions is constant), and
another that assumes the variance-covariance matrix is symmetric (i.e.,
the correlations between conditions are allowed to differ). Do you
prefer one model over the other? Why or why not?# with lme:
library(nlme)
# compound symmetry
smile.lme.01 <- lme(leniency~smile,
data=smile.df,
random=~1|subjID,
correlation=corCompSymm(form=~1|subjID))
print(anova(smile.lme.01))
## numDF denDF F-value p-value
## (Intercept) 1 57 672.3 <.0001
## smile 3 57 3.6 0.0184
# general/symmetric var-covar matrix
smile.lme.02 <- lme(leniency~smile,
data=smile.df,
random=~1|subjID,
correlation=corSymm(form = ~1|subjID))
print( anova(smile.lme.02) )
## numDF denDF F-value p-value
## (Intercept) 1 57 598.4 <.0001
## smile 3 57 4.9 0.0041
# now compare the two models:
print(anova(smile.lme.01,smile.lme.02) )
## Model df AIC BIC logLik Test L.Ratio p-value
## smile.lme.01 1 7 626.3 642.6 -306.1
## smile.lme.02 2 12 632.2 660.2 -304.1 1 vs 2 4.077 0.5384
# The goodness-of-fit does not differ significantly between models
# Based on that result, I prefer the simpler model smile.lme.01
# Note: the L.Ratio is distributed approximately as a Chi-square variable.
# Although the chi-square test is slightly conservative, I would still prefer the simpler modelsubjID.library(VCA)
vca.aov.01 <- anovaMM(leniency~smile+(subjID),Data=smile.df)
print(vca.aov.01)
##
##
## ANOVA-Type Estimation of Mixed Model:
## --------------------------------------
##
## [Fixed Effects]
##
## int smilecontempt smilemiserable smileneutral smilereal
## 56 -11 -5 -8 0
##
##
## [Variance Components]
##
## Name DF SS MS VC %Total SD CV[%]
## 1 total 62.735217 165.609058 100 12.868918 25.737837
## 2 subjID 19 5652.639108 297.507321 43.966088 26.548118 6.630693 13.261386
## 3 error 57 6933.649276 121.64297 121.64297 73.451882 11.029187 22.058374
##
## Mean: 50 (N = 80)
##
## Experimental Design: balanced | Method: ANOVA
( subject.icc <- 43.96 / (43.96+121.64) )
## [1] 0.2655
# Note that this value is the same as one estimated from my lmer model:
library(performance)
smile.lmer.01 <- lmer(leniency~smile+(1|subjID),data=smile.df)
print(icc(smile.lmer.01,by_group = T))
## # ICC by Group
##
## Group | ICC
## --------------
## subjID | 0.265real smile condition did not differ from average
leniency in the other three conditions.with(smile.df,tapply(leniency,smile,mean)) # means
## neutral contempt miserable real
## 48 45 51 56
names(smile.wide)
## [1] "neutral" "contempt" "miserable" "real"
w <- c(-1/3, -1/3, -1/3, 1) # real smile vs other 3 conditions
# smile.mat <- as.matrix(smile.wide)
contrast.scores <- smile.mat %*% w
t.test(contrast.scores) # contrast is significant
##
## One Sample t-test
##
## data: contrast.scores
## t = 2.7, df = 19, p-value = 0.01
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## 1.768 14.232
## sample estimates:
## mean of x
## 8
# The contrast is significant.
# I reject null hypothesis that leniency in real smile condition
# is equal to average leniency in the other conditions
# with emmeans & aov_car: (local error term)
library(emmeans)
smile.aov.emm <- emmeans(smile.aov.01,specs="smile")
contrast(smile.aov.emm,method=list("real-vs-others"=w))
## contrast estimate SE df t.ratio p.value
## real-vs-others 8 2.98 19 2.687 0.0146
# with emmeans & lmer: (global error term)
smile.lmer.emm <- emmeans(smile.lmer.01,specs="smile")
contrast(smile.lmer.emm,method=list("real-vs-others"=w))
## contrast estimate SE df t.ratio p.value
## real-vs-others 8 2.85 57 2.809 0.0068
##
## Degrees-of-freedom method: kenward-roger