Load libraries
library(readxl)
library(reshape2)
library(lme4)
library(stargazer)
library(effects)
library(MASS)
library(lme4)
library(car)
library(MuMIn)
library(plotrix)
library(ggplot2)Read in all data (active and passive) from DataAll.xlsx and bind together to one dataframe (df)
df.active <- read_xlsx(path = 'DataAll.xlsx',sheet = 1,n_max = 34)
# build long-format df
df.active.re <- data.frame(PPnum = rep(df.active$PPnum,3*2),
Intelligibility = rep(rep(c(20,50,80),each=34),2),
State = rep(c('Alone','Together'),each=3*34),
PPD = c(df.active$`A - 20`,
df.active$`A - 50`,
df.active$`A 80`,
df.active$`T - 20`,
df.active$`T - 50`,
df.active$`T - 80`))
# include task identifier
df.active.re$Task <- 'Active'
df.passive <- read_xlsx(path = 'DataAll.xlsx',sheet = 2,n_max = 34)
# build long-format df
df.passive.re <- data.frame(PPnum = rep(df.passive$PPnum,3*2),
Intelligibility = rep(rep(c(20,50,80),each=34),2),
State = rep(c('Alone','Together'),each=3*34),
PPD = c(df.passive$`A - 20`,
df.passive$`A - 50`,
df.passive$`A 80`,
df.passive$`T - 20`,
df.passive$`T - 50`,
df.passive$`T - 80`))
# include task identifier
df.passive.re$Task <- 'Passive'
# full dataframe
df <- rbind(df.active.re,df.passive.re)Summary of data (means + SD)
df.agg <- aggregate(PPD~Intelligibility+State+Task,mean,data=df)
df.agg$CI <- qnorm(0.975)*aggregate(PPD~Intelligibility+State+Task,std.error,data=df)$PPD
print(ggplot(df.agg,aes(x=Intelligibility,y=PPD,fill=State))+
geom_col(position = position_dodge2())+
geom_errorbar(
aes(ymin = PPD - CI, ymax = PPD + CI),
position = position_dodge2()) +
scale_x_continuous(breaks = c(20,50,80))+
theme_bw()+
labs(title='Mean pupil dilation ±95% CI',y='PPD (mm)') +
facet_wrap(~Task))
Subsetting and scaling of continuous predictor
# subsetting: no subsetting applied
df.sub <- df
# scale and center intelligibility (continuous variable)
df.sub$Intelligibility <- scale(df.sub$Intelligibility)Define hypotheses and fit associated LMM’s
fit.0 <- lmer(PPD~1+(1|PPnum),data = df.sub)
cat(paste('AIC null model: ', AIC(fit.0)))
AIC null model: -429.543872160715
fit.1 <- lmer(PPD~Intelligibility*State*Task+(1|PPnum),data = df.sub,REML = T)
cat(paste('AIC model 1: ', AIC(fit.1)))
AIC model 1: -632.638051924798
fit.2 <- lmer(PPD~Intelligibility*State*Task+(Task|PPnum),data = df.sub,REML = T)
cat(paste('AIC model 2: ', AIC(fit.2)))
AIC model 2: -725.207622133983
fit.3 <- lmer(PPD~Intelligibility*State*Task+(State|PPnum),data = df.sub,REML = T)
cat(paste('AIC model 3: ', AIC(fit.3)))
AIC model 3: -631.299678395452Evaluate best fitting model (model 2)
# test if selected model is better than the null (i.e. are the added parameters necessary?)
anova(fit.0,fit.2)
Data: df.sub
Models:
fit.0: PPD ~ 1 + (1 | PPnum)
fit.2: PPD ~ Intelligibility * State * Task + (Task | PPnum)
Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
fit.0 3 -436.29 -424.96 221.14 -442.29
fit.2 12 -787.19 -741.90 405.60 -811.19 368.9 9 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Print fitted coefficients together with p-values
stargazer(fit.2,type = 'text',align = T)
=====================================================================
Dependent variable:
---------------------------
PPD
---------------------------------------------------------------------
Intelligibility -0.034***
(0.007)
StateTogether 0.028***
(0.009)
TaskPassive -0.152***
(0.018)
Intelligibility:StateTogether 0.009
(0.009)
Intelligibility:TaskPassive 0.036***
(0.009)
StateTogether:TaskPassive -0.023*
(0.013)
Intelligibility:StateTogether:TaskPassive -0.015
(0.013)
Constant 0.220***
(0.022)
---------------------------------------------------------------------
Observations 322
Log Likelihood 374.604
Akaike Inf. Crit. -725.208
Bayesian Inf. Crit. -679.913
=====================================================================
Note: *p<0.1; **p<0.05; ***p<0.01The main effects are:
Intelligibility: increase with one SD (30) intelligibility results in 0.034 decreased pupil dilation
Social state: performing the task together results in 0.028 increased pupil dilation
Task: when performing the task passively, pupil dilation is reduced with 0.152
In addition, the interaction between the effect of intelligibility and task (i.e. active or passive) is positive and significant.
Residual QQ-plot
print(qqPlot(resid(fit.2),main='Linear'))
291 143
229 111 Estimate proportion of explained variance
cat(paste('Marginal R squared value associated with fixed effects: ',round(r.squaredGLMM(fit.2)[1],3)))
Marginal R squared value associated with fixed effects: 0.413
cat(paste('Conditional R2 value associated with fixed + random effects: ',round(r.squaredGLMM(fit.2)[2],3)))
Conditional R2 value associated with fixed + random effects: 0.8Variance inflation factor (<5 acceptable)
round(vif(fit.2),3)
Intelligibility State
3.989 2.003
Task Intelligibility:State
1.152 3.990
Intelligibility:Task State:Task
3.995 2.154
Intelligibility:State:Task
3.996