Simulator Sickness - Mixed ANOVA
Nate hatfield
July 25th, 2024
1 Overview
This code will walkthrough how to run a mixed effects ANOVA. For this particular example, I will explore the effects of different visual patterns on subjective simulator sickness scores over time. The experimental setup follows a 4 (between subjects: control, horizon, dot, coil) x 5 (within subjects: baseline, 5 mins, 10 mins, 15 mins, 20 mins) mixed design.
2 Preparing Data for Analysis
Lets first start by ensuring the grouping factor is classified as categorical and that all N/As are removed from this dataset for the sake of demonstration. Please note that handling outliers, missing data, and non-normal distributions are all part of the data preparation process, though skipped in this exercise.
#Load the packages / libraries
library(rstudioapi)
library(tidyverse) #this package is used for data cleaning
#Import the dataset and preparing it for analysis.
Prompt_tbl <- read_csv("/Users/4792241/Prompt/Data/Prompt_Data.csv")%>%
mutate(ID = as.factor(ID))%>%
mutate(UI_Condition = as.factor(Condition))%>%
mutate(Duration = as.factor(Duration)) %>%
drop_na()
#Reorder Labels
Prompt_tbl$UI_Condition <- ordered(Prompt_tbl$UI_Condition, levels=c("Control", "Horizon", "Dot", "Coil"))
Prompt_tbl$Duration <- ordered(Prompt_tbl$Duration, levels=c("BL", "5", "10", "15", "20"))
#Look at the structure
str(Prompt_tbl)## tibble [160 × 5] (S3: tbl_df/tbl/data.frame)
## $ ID : Factor w/ 32 levels "1","2","3","4",..: 1 1 1 1 1 2 2 2 2 2 ...
## $ Condition : chr [1:160] "Control" "Control" "Control" "Control" ...
## $ Duration : Ord.factor w/ 5 levels "BL"<"5"<"10"<..: 1 2 3 4 5 1 2 3 4 5 ...
## $ SSQ : num [1:160] 2 23 25 32 42 3 17 31 47 53 ...
## $ UI_Condition: Ord.factor w/ 4 levels "Control"<"Horizon"<..: 1 1 1 1 1 1 1 1 1 1 ...#Examine first and last few cases in dataframe
psych::headtail(Prompt_tbl)## ID Condition Duration SSQ UI_Condition
## 1 1 Control BL 2 Control
## 2 1 Control 5 23 Control
## 3 1 Control 10 25 Control
## 4 1 Control 15 32 Control
## 5 <NA> <NA> <NA> ... <NA>
## 6 32 Coil 5 2 Coil
## 7 32 Coil 10 5 Coil
## 8 32 Coil 15 8 Coil
## 9 32 Coil 20 12 Coil3 Inspect Data with Visualizations
#Quickly inspect data
boxplot(SSQ~UI_Condition*Duration,data=Prompt_tbl)
#Quickly inspect data
boxplot(SSQ~UI_Condition*Duration,data=Prompt_tbl)
# Create the boxplot using ggplot2 with color-coded UI_Condition
ggplot(Prompt_tbl, aes(x = interaction(UI_Condition, Duration), y = SSQ, fill = UI_Condition)) +
geom_boxplot() +
labs(x = "UI Condition and Duration",
y = "SSQ",
title = "Boxplot of SSQ by UI Condition and Duration",
fill = "UI Condition") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) # Rotates the x-axis labels for better readability
#Create Interaction Plot
interaction.plot(Prompt_tbl$UI_Condition, Prompt_tbl$Duration, Prompt_tbl$SSQ)
interaction.plot(Prompt_tbl$Duration, Prompt_tbl$UI_Condition, Prompt_tbl$SSQ)
4 Create Summary Statistics
#Load packages
library(Rmisc)
#Create Summary Data
summary1 <- summarySE(data=Prompt_tbl,
measurevar="SSQ",
groupvars=c("Duration","UI_Condition"),
conf.interval=.95)
summary1## Duration UI_Condition N SSQ sd se ci
## 1 BL Control 8 2.250 0.7071068 0.2500000 0.5911561
## 2 BL Horizon 8 2.000 1.6035675 0.5669467 1.3406159
## 3 BL Dot 8 1.125 1.2464235 0.4406772 1.0420361
## 4 BL Coil 8 1.375 0.9161254 0.3238992 0.7659000
## 5 5 Control 8 16.125 4.3895167 1.5519285 3.6697278
## 6 5 Horizon 8 8.625 2.7742438 0.9808433 2.3193258
## 7 5 Dot 8 6.875 1.2464235 0.4406772 1.0420361
## 8 5 Coil 8 1.875 0.9910312 0.3503824 0.8285228
## 9 10 Control 8 23.250 5.5741751 1.9707685 4.6601270
## 10 10 Horizon 8 11.375 2.6692696 0.9437293 2.2315652
## 11 10 Dot 8 9.250 1.9820624 0.7007649 1.6570456
## 12 10 Coil 8 3.625 1.9226098 0.6797452 1.6073420
## 13 15 Control 8 33.250 7.1264096 2.5195663 5.9578276
## 14 15 Horizon 8 18.375 3.7772817 1.3354708 3.1578865
## 15 15 Dot 8 13.750 2.7124054 0.9589801 2.2676276
## 16 15 Coil 8 5.250 1.9820624 0.7007649 1.6570456
## 17 20 Control 8 42.500 6.8243262 2.4127637 5.7052795
## 18 20 Horizon 8 22.750 3.1509636 1.1140339 2.6342715
## 19 20 Dot 8 19.000 3.4226139 1.2100767 2.8613768
## 20 20 Coil 8 10.000 2.2038927 0.7791937 1.8425004#bOXPLOT TO SHOW MEANS OF EACH CONDITION
grey_palette <- c("#FFFFFF", "#D3D3D3", "#696969", "black") # Adjust the number of colors as needed
# Create the boxplot with mean points and mean text labels
ggplot(data = Prompt_tbl, aes(x = Duration, y = SSQ, fill = UI_Condition)) +
geom_boxplot() +
geom_point(
data = summary1,
aes(y = SSQ, x = Duration),
position = position_dodge(width = .75),
color = "red",
size = 1.5) +
geom_text(
data = summary1,
aes(y = SSQ, x = Duration, label = round(SSQ, 1)),
position = position_dodge(width = .75),
vjust = -0.5,
color = "black",
size = 2.5) +
scale_fill_manual(values = grey_palette) +
xlab("Duration") +
ylab("Average Total SSQ Score") +
ggtitle("Boxplot of SSQ Score") +
labs(subtitle = "Red Points represent the mean") + # Add subtitle here
theme_minimal() +
theme(text = element_text(size = 12))
#Plot effects over time
ggplot(summary1, aes(x = Duration, y = SSQ, group = UI_Condition)) +
geom_line(aes(linetype = UI_Condition), color = "black") + # Keep lines black
geom_point(aes(shape = UI_Condition), color = "black", size = 3) + # Set points color to black
geom_errorbar(aes(ymin = SSQ - ci, ymax = SSQ + ci), colour = "black", width = .1) +
xlab("Duration") +
ylab("Mean Total SSQ Score") +
ggtitle("Effects of UI Condition and Duration on Simulator Sickness Scores") +
labs(subtitle = "Error Bars Represent 95% Confidence Interval") +
ylim(0, 60) +
theme_classic() +
theme(
text = element_text(size = 12),
panel.background = element_blank(),
plot.background = element_blank(),
legend.background = element_blank(),
legend.title = element_text(size = 12, face = "bold"),
legend.text = element_text(size = 10)
) +
scale_color_brewer(palette = "Set1") +
scale_shape_manual(values = c(16, 17, 18, 19)) +
scale_linetype_manual(values = c("solid", "twodash", "longdash", "dotted"))
ggplot(summary1, aes(x = Duration, y = SSQ, group = UI_Condition, color = UI_Condition)) +
geom_line(aes(linetype = UI_Condition)) +
geom_point(size = 2) +
geom_errorbar(aes(ymin = SSQ - ci, ymax = SSQ + ci), colour = "black", width = .1) +
xlab("Duration") +
ylab("Mean Total SSQ Score") +
ggtitle("Effects of UI Condition and Duration on Simulator Sickness Scores") +
labs(subtitle = "Error Bars Represent 95% Confidence Interval") + # Add subtitle here
ylim(0, 60) + # Set y-axis limits from 0 to 60
theme_classic() +
theme(
text = element_text(size = 12),
panel.background = element_blank(), # Remove background
plot.background = element_blank(), # Remove plot background
legend.background = element_blank(), # Remove legend background
legend.title = element_text(size = 12, face = "bold"),
legend.text = element_text(size = 10)
) +
scale_color_brewer(palette = "Set2") +
scale_linetype_manual(values = c("solid", "twodash", "longdash", "dotted"))
ggplot(summary1, aes(x = Duration, y = SSQ, group = UI_Condition)) +
geom_line(aes(linetype = UI_Condition), color = "black") + # Keep lines black
geom_point(aes(shape = UI_Condition), color = "black", size = 3) + # Set points color to black
geom_errorbar(aes(ymin = SSQ - ci, ymax = SSQ + ci), colour = "black", width = .1) +
xlab("Duration") +
ylab("Mean Total SSQ Score") +
ggtitle("Effects of UI Condition and Duration on Simulator Sickness Scores") +
labs(subtitle = "Error Bars Represent 95% Confidence Interval") +
ylim(0, 60) +
theme_classic() +
theme(
text = element_text(size = 12),
panel.background = element_blank(),
plot.background = element_blank(),
legend.background = element_blank(),
legend.title = element_text(size = 12, face = "bold"),
legend.text = element_text(size = 10)
) +
scale_color_brewer(palette = "Set1") +
scale_shape_manual(values = c(16, 17, 18, 19)) +
scale_linetype_manual(values = c("solid", "twodash", "longdash", "dotted"))+
# Add extremely transparent shaded regions
geom_rect(aes(xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = 5), fill = "lightblue", alpha = 0.01) +
geom_rect(aes(xmin = -Inf, xmax = Inf, ymin = 5, ymax = 10), fill = "lightgreen", alpha = 0.01) +
geom_rect(aes(xmin = -Inf, xmax = Inf, ymin = 10, ymax = 15), fill = "yellow", alpha = 0.01) +
geom_rect(aes(xmin = -Inf, xmax = Inf, ymin = 15, ymax = 20), fill = "orange", alpha = 0.01) +
geom_rect(aes(xmin = -Inf, xmax = Inf, ymin = 20, ymax = Inf), fill = "red", alpha = 0.01)
ggplot(summary1, aes(x = Duration, y = SSQ, group = UI_Condition)) +
geom_line(aes(linetype = UI_Condition), color = "black") + # Keep lines black
geom_point(aes(shape = UI_Condition), color = "black", size = 2) + # Set points color to black and decrease size
geom_errorbar(aes(ymin = SSQ - ci, ymax = SSQ + ci), colour = "black", width = .1) +
geom_text(aes(label = round(SSQ, 1)), vjust = -1, size = 3.5, color = "black") + # Add mean text labels
xlab("Duration") +
ylab("Mean Total SSQ Score") +
ggtitle("Effects of UI Condition and Duration on Simulator Sickness Scores") +
labs(subtitle = "Error Bars Represent 95% Confidence Interval") +
ylim(0, 60) +
theme_classic() +
theme(
text = element_text(size = 12),
panel.background = element_blank(),
plot.background = element_blank(),
legend.background = element_blank(),
legend.title = element_text(size = 12, face = "bold"),
legend.text = element_text(size = 10)
) +
facet_wrap(~UI_Condition)+
scale_color_brewer(palette = "Set1") +
scale_shape_manual(values = c(16, 17, 18, 19)) +
scale_linetype_manual(values = c("solid", "twodash", "longdash", "dotted")) +
# Add extremely transparent shaded regions
geom_rect(aes(xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = 5), fill = "lightblue", alpha = 0.03) +
geom_rect(aes(xmin = -Inf, xmax = Inf, ymin = 5, ymax = 10), fill = "lightgreen", alpha = 0.03) +
geom_rect(aes(xmin = -Inf, xmax = Inf, ymin = 10, ymax = 15), fill = "yellow", alpha = 0.03) +
geom_rect(aes(xmin = -Inf, xmax = Inf, ymin = 15, ymax = 20), fill = "orange", alpha = 0.03) +
geom_rect(aes(xmin = -Inf, xmax = Inf, ymin = 20, ymax = Inf), fill = "red", alpha = 0.03)
summary_Duration <- summarySE(data=Prompt_tbl,
measurevar="SSQ",
groupvars=c("Duration"),
conf.interval=.95)
summary_Duration## Duration N SSQ sd se ci
## 1 BL 32 1.68750 1.203154 0.2126897 0.4337834
## 2 5 32 8.37500 5.801835 1.0256292 2.0917846
## 3 10 32 11.87500 7.946393 1.4047371 2.8649802
## 4 15 32 17.65625 11.125702 1.9667649 4.0112434
## 5 20 32 23.56250 12.730298 2.2504200 4.5897618ggplot(summary_Duration, aes(x = Duration, y = SSQ)) +
geom_line(aes(linetype = Duration), show.legend = FALSE) +
geom_point(size = 2) +
geom_errorbar(aes(ymin = SSQ - ci, ymax = SSQ + ci), colour = "black", width = 0.1) +
geom_text(aes(label = round(SSQ, 1)), hjust = 1.5, color = "black", size = 3.5) + # Add mean text labels
xlab("Duration") +
ylab("Mean Total SSQ Score") +
ggtitle("Average Simulator Sickness Over Time Collapsed Across all UI Conditions") +
labs(subtitle = "Error Bars Represent 95% Confidence Interval") + # Add subtitle here
ylim(0, 60) + # Set y-axis limits from 0 to 60
theme_classic() +
theme(
text = element_text(size = 12),
panel.background = element_blank(), # Remove background
plot.background = element_blank(), # Remove plot background
legend.background = element_blank() # Remove legend background
)
summary_UICondition <- summarySE(data=Prompt_tbl,
measurevar="SSQ",
groupvars=c("UI_Condition"),
conf.interval=.95)
ggplot(summary_UICondition, aes(x = UI_Condition, y = SSQ)) +
geom_line(aes(linetype = UI_Condition), show.legend = FALSE) +
geom_point(size = 2) +
geom_errorbar(aes(ymin = SSQ - ci, ymax = SSQ + ci), colour = "black", width = 0.1) +
geom_text(aes(label = round(SSQ, 1)), hjust = 1.5, color = "black", size = 3.5) + # Add mean text labels
xlab("Duration") +
ylab("Mean Total SSQ Score") +
ggtitle("Average Simulator Sickness Per UI Condition Collapsed Across all Time Duration") +
labs(subtitle = "Error Bars Represent 95% Confidence Interval") + # Add subtitle here
ylim(0, 60) + # Set y-axis limits from 0 to 60
theme_classic() +
theme(
text = element_text(size = 12),
panel.background = element_blank(), # Remove background
plot.background = element_blank(), # Remove plot background
legend.background = element_blank() # Remove legend background
)
# Perform Inferential Statistics
#Libraries
library(compute.es)
library(ez)
library(ggplot2)
library(multcomp)
library(nlme)
library(pastecs)
library(reshape)
library(ggpubr)
library(rstatix)
library(Rmisc)
library(viridis)
##Perform Base ANOVA
fit1.aov <- stats::aov(SSQ~Duration*UI_Condition + Error(ID/Duration),Prompt_tbl)
summary(fit1.aov)##
## Error: ID
## Df Sum Sq Mean Sq F value Pr(>F)
## UI_Condition 3 7674 2558.0 80.59 6.93e-14 ***
## Residuals 28 889 31.7
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: ID:Duration
## Df Sum Sq Mean Sq F value Pr(>F)
## Duration 4 9062 2265.6 351.13 <2e-16 ***
## Duration:UI_Condition 12 2621 218.5 33.86 <2e-16 ***
## Residuals 112 723 6.5
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#Report means from the fit object
print(model.tables(fit1.aov,"means"),digits=3) #report the means of each level## Tables of means
## Grand mean
##
## 12.63125
##
## Duration
## Duration
## BL 5 10 15 20
## 1.69 8.37 11.88 17.66 23.56
##
## UI_Condition
## UI_Condition
## Control Horizon Dot Coil
## 23.47 12.62 10.00 4.42
##
## Duration:UI_Condition
## UI_Condition
## Duration Control Horizon Dot Coil
## BL 2.2 2.0 1.1 1.4
## 5 16.1 8.6 6.9 1.9
## 10 23.3 11.4 9.3 3.6
## 15 33.2 18.4 13.7 5.2
## 20 42.5 22.7 19.0 10.0##Report deviation effects for each level
print(model.tables(fit1.aov,"effects", n=TRUE),digits=3) # report deviation effects for each level## Tables of effects
##
## Duration
## Duration
## BL 5 10 15 20
## -10.94 -4.26 -0.76 5.02 10.93
##
## UI_Condition
## UI_Condition
## Control Horizon Dot Coil
## 10.84 -0.01 -2.63 -8.21
##
## Duration:UI_Condition
## UI_Condition
## Duration Control Horizon Dot Coil
## BL -10.28 0.32 2.07 7.89
## 5 -3.09 0.26 1.13 1.71
## 10 0.53 -0.49 0.01 -0.04
## 15 4.75 0.73 -1.27 -4.20
## 20 8.09 -0.81 -1.93 -5.36##Use EXANOVA
fit2.ez <- ezANOVA(
data=Prompt_tbl,
dv=.(SSQ),
wid=.(ID),
within=.(Duration),
between=.(UI_Condition),
type=3,
detailed=TRUE,
return_aov=TRUE
)
print(fit2.ez)## $ANOVA
## Effect DFn DFd SSn SSd F p p<.05
## 1 (Intercept) 1 28 25527.756 888.725 804.27261 3.613323e-22 *
## 2 UI_Condition 3 28 7674.119 888.725 80.59311 6.931801e-14 *
## 3 Duration 4 112 9062.275 722.650 351.12945 2.248589e-62 *
## 4 UI_Condition:Duration 12 112 2621.475 722.650 33.85747 9.919584e-32 *
## ges
## 1 0.9406254
## 2 0.8264632
## 3 0.8490324
## 4 0.6193168
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 3 Duration 0.4070385 0.004801773 *
## 4 UI_Condition:Duration 0.4070385 0.004801773 *
##
## $`Sphericity Corrections`
## Effect GGe p[GG] p[GG]<.05 HFe p[HF]
## 3 Duration 0.7047335 1.085896e-44 * 0.7917565 6.630686e-50
## 4 UI_Condition:Duration 0.7047335 5.371898e-23 * 0.7917565 1.423881e-25
## p[HF]<.05
## 3 *
## 4 *
##
## $aov
##
## Call:
## aov(formula = formula(aov_formula), data = data)
##
## Grand Mean: 12.63125
##
## Stratum 1: ID
##
## Terms:
## UI_Condition Residuals
## Sum of Squares 7674.119 888.725
## Deg. of Freedom 3 28
##
## Residual standard error: 5.633842
## Estimated effects are balanced
##
## Stratum 2: ID:Duration
##
## Terms:
## Duration UI_Condition:Duration Residuals
## Sum of Squares 9062.275 2621.475 722.650
## Deg. of Freedom 4 12 112
##
## Residual standard error: 2.540124
## Estimated effects are balanced#Rstatic approach
fit6.aovtest <- anova_test(data=Prompt_tbl,
dv=SSQ,
wid=ID,
within=Duration,
between=UI_Condition)
fit6.aovtest## ANOVA Table (type II tests)
##
## $ANOVA
## Effect DFn DFd F p p<.05 ges
## 1 UI_Condition 3 28 80.593 6.93e-14 * 0.826
## 2 Duration 4 112 351.129 2.25e-62 * 0.849
## 3 UI_Condition:Duration 12 112 33.857 9.92e-32 * 0.619
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 1 Duration 0.407 0.005 *
## 2 UI_Condition:Duration 0.407 0.005 *
##
## $`Sphericity Corrections`
## Effect GGe DF[GG] p[GG] p[GG]<.05 HFe DF[HF]
## 1 Duration 0.705 2.82, 78.93 1.09e-44 * 0.792 3.17, 88.68
## 2 UI_Condition:Duration 0.705 8.46, 78.93 5.37e-23 * 0.792 9.5, 88.68
## p[HF] p[HF]<.05
## 1 6.63e-50 *
## 2 1.42e-25 *#Brief results summaries
get_anova_table(fit6.aovtest, correction = "GG")## ANOVA Table (type II tests)
##
## Effect DFn DFd F p p<.05 ges
## 1 UI_Condition 3.00 28.00 80.593 6.93e-14 * 0.826
## 2 Duration 2.82 78.93 351.129 1.09e-44 * 0.849
## 3 UI_Condition:Duration 8.46 78.93 33.857 5.37e-23 * 0.619#Contrasts
contrasts(Prompt_tbl$UI_Condition) <- contr.sum
contrasts(Prompt_tbl$UI_Condition)## [,1] [,2] [,3]
## Control 1 0 0
## Horizon 0 1 0
## Dot 0 0 1
## Coil -1 -1 -1contrasts(Prompt_tbl$Duration) <- contr.sum
contrasts(Prompt_tbl$Duration)## [,1] [,2] [,3] [,4]
## BL 1 0 0 0
## 5 0 1 0 0
## 10 0 0 1 0
## 15 0 0 0 1
## 20 -1 -1 -1 -1options(repos = c(CRAN = "https://cloud.r-project.org/"))
install.packages("afex")##
## The downloaded binary packages are in
## /var/folders/l3/0wmdqfhn3371l3q_w5kbk3tm0000gq/T//RtmpnZ7yf3/downloaded_packageslibrary(afex)
fit3.afex <- aov_car(SSQ~UI_Condition*Duration + Error(ID/Duration), data=Prompt_tbl)
summary(fit3.afex)##
## Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
##
## Sum Sq num Df Error SS den Df F value Pr(>F)
## (Intercept) 25527.8 1 888.73 28 804.273 < 2.2e-16 ***
## UI_Condition 7674.1 3 888.73 28 80.593 6.932e-14 ***
## Duration 9062.3 4 722.65 112 351.130 < 2.2e-16 ***
## UI_Condition:Duration 2621.5 12 722.65 112 33.858 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Mauchly Tests for Sphericity
##
## Test statistic p-value
## Duration 0.40704 0.0048018
## UI_Condition:Duration 0.40704 0.0048018
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## Duration 0.70473 < 2.2e-16 ***
## UI_Condition:Duration 0.70473 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## HF eps Pr(>F[HF])
## Duration 0.7917565 6.630686e-50
## UI_Condition:Duration 0.7917565 1.423881e-25anova(fit3.afex, correction="GG")## Anova Table (Type 3 tests)
##
## Response: SSQ
## num Df den Df MSE F ges Pr(>F)
## UI_Condition 3.0000 28.00 31.740 80.593 0.82646 6.932e-14 ***
## Duration 2.8189 78.93 9.156 351.130 0.84903 < 2.2e-16 ***
## UI_Condition:Duration 8.4568 78.93 9.156 33.858 0.61932 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1install.packages("gt")##
## There is a binary version available but the source version is later:
## binary source needs_compilation
## gt 0.10.1 0.11.0 FALSElibrary(gt)
gt(nice(anova(fit3.afex, correction="none")))| Effect | df | MSE | F | ges | p.value |
|---|---|---|---|---|---|
| UI_Condition | 3, 28 | 31.74 | 80.59 *** | .826 | <.001 |
| Duration | 4, 112 | 6.45 | 351.13 *** | .849 | <.001 |
| UI_Condition:Duration | 12, 112 | 6.45 | 33.86 *** | .619 | <.001 |
gt(nice(anova(fit3.afex, correction="GG")))| Effect | df | MSE | F | ges | p.value |
|---|---|---|---|---|---|
| UI_Condition | 3, 28 | 31.74 | 80.59 *** | .826 | <.001 |
| Duration | 2.82, 78.93 | 9.16 | 351.13 *** | .849 | <.001 |
| UI_Condition:Duration | 8.46, 78.93 | 9.16 | 33.86 *** | .619 | <.001 |
fit4.afex <- aov_car(SSQ~UI_Condition*Duration + Error(ID/Duration), data=Prompt_tbl,
anova_table = list(es = "pes",
correction="none"))
gt(nice(fit4.afex))| Effect | df | MSE | F | pes | p.value |
|---|---|---|---|---|---|
| UI_Condition | 3, 28 | 31.74 | 80.59 *** | .896 | <.001 |
| Duration | 4, 112 | 6.45 | 351.13 *** | .926 | <.001 |
| UI_Condition:Duration | 12, 112 | 6.45 | 33.86 *** | .784 | <.001 |
#Interaction Term Contrasts
install.packages("emmeans")##
## There is a binary version available but the source version is later:
## binary source needs_compilation
## emmeans 1.10.2 1.10.3 FALSElibrary(emmeans)
emm1.int <- emmeans(fit3.afex, specs=c("Duration","UI_Condition"))
contrast(emm1.int, interaction=c(Duration="poly", UI_Condition= c("trt.vs.ctrl")))## Duration_poly UI_Condition_trt.vs.ctrl estimate SE df t.ratio p.value
## linear Horizon - Control -46.38 5.59 28 -8.299 <.0001
## quadratic Horizon - Control 6.12 4.97 28 1.233 0.2279
## cubic Horizon - Control -4.75 3.07 28 -1.546 0.1333
## quartic Horizon - Control -1.75 6.60 28 -0.265 0.7929
## linear Dot - Control -55.00 5.59 28 -9.843 <.0001
## quadratic Dot - Control 7.50 4.97 28 1.510 0.1423
## cubic Dot - Control -1.88 3.07 28 -0.610 0.5466
## quartic Dot - Control 6.38 6.60 28 0.965 0.3426
## linear Coil - Control -77.00 5.59 28 -13.780 <.0001
## quadratic Coil - Control 14.75 4.97 28 2.969 0.0061
## cubic Coil - Control -4.12 3.07 28 -1.343 0.1902
## quartic Coil - Control 17.88 6.60 28 2.707 0.0114emm <- emmeans(fit3.afex, specs = c("Duration", "UI_Condition"))
# Define the contrasts for the interaction
# Here we specify polynomial contrasts for the within-subjects factor and treatment vs. control comparisons for the between-subjects factor
interaction_contrasts <- list(
Duration = "poly",
UI_Condition = list(
"Control - Horizon" = c(-1, 1, 0, 0),
"Dot - Control" = c(-1, 0, 1, 0),
"Dot - Horizon" = c(0, -1, 1, 0),
"Dot - Coil" = c(0, 0, -1, 1),
"All - Control" = c(1/3, 1/3, 1/3, -1) # Comparing all conditions to Control
)
)
# Perform the interaction contrasts
interaction_results <- contrast(emm, interaction = interaction_contrasts)
# View the results
summary(interaction_results)## Duration_poly UI_Condition_custom estimate SE df t.ratio p.value
## linear Control - Horizon -46.38 5.59 28 -8.299 <.0001
## quadratic Control - Horizon 6.12 4.97 28 1.233 0.2279
## cubic Control - Horizon -4.75 3.07 28 -1.546 0.1333
## quartic Control - Horizon -1.75 6.60 28 -0.265 0.7929
## linear Dot - Control -55.00 5.59 28 -9.843 <.0001
## quadratic Dot - Control 7.50 4.97 28 1.510 0.1423
## cubic Dot - Control -1.88 3.07 28 -0.610 0.5466
## quartic Dot - Control 6.38 6.60 28 0.965 0.3426
## linear Dot - Horizon -8.62 5.59 28 -1.544 0.1339
## quadratic Dot - Horizon 1.38 4.97 28 0.277 0.7840
## cubic Dot - Horizon 2.88 3.07 28 0.936 0.3574
## quartic Dot - Horizon 8.12 6.60 28 1.230 0.2288
## linear Dot - Coil -22.00 5.59 28 -3.937 0.0005
## quadratic Dot - Coil 7.25 4.97 28 1.459 0.1556
## cubic Dot - Coil -2.25 3.07 28 -0.732 0.4701
## quartic Dot - Coil 11.50 6.60 28 1.741 0.0926
## linear All - Control 43.21 4.56 28 9.470 <.0001
## quadratic All - Control -10.21 4.06 28 -2.517 0.0179
## cubic All - Control 1.92 2.51 28 0.764 0.4512
## quartic All - Control -16.33 5.39 28 -3.029 0.0052#Main Effect Contrast
emm1.duration <- emmeans(fit3.afex, specs="Duration")
emm1.duration## Duration emmean SE df lower.CL upper.CL
## BL 1.69 0.207 28 1.26 2.11
## X5 8.38 0.480 28 7.39 9.36
## X10 11.88 0.598 28 10.65 13.10
## X15 17.66 0.772 28 16.07 19.24
## X20 23.56 0.756 28 22.01 25.11
##
## Results are averaged over the levels of: UI_Condition
## Confidence level used: 0.95contrast(emm1.duration, "poly")## contrast estimate SE df t.ratio p.value
## linear 53.031 1.98 28 26.843 <.0001
## quadratic 0.719 1.76 28 0.409 0.6855
## cubic 3.312 1.09 28 3.049 0.0050
## quartic -7.625 2.33 28 -3.266 0.0029
##
## Results are averaged over the levels of: UI_Condition#Simple main effects
emm1.sme <- emmeans(fit3.afex, "Duration", by="UI_Condition")
emm1.sme## UI_Condition = Control:
## Duration emmean SE df lower.CL upper.CL
## BL 2.25 0.413 28 1.4035 3.10
## X5 16.12 0.960 28 14.1582 18.09
## X10 23.25 1.197 28 20.7988 25.70
## X15 33.25 1.545 28 30.0861 36.41
## X20 42.50 1.511 28 39.4046 45.60
##
## UI_Condition = Horizon:
## Duration emmean SE df lower.CL upper.CL
## BL 2.00 0.413 28 1.1535 2.85
## X5 8.62 0.960 28 6.6582 10.59
## X10 11.38 1.197 28 8.9238 13.83
## X15 18.38 1.545 28 15.2111 21.54
## X20 22.75 1.511 28 19.6546 25.85
##
## UI_Condition = Dot:
## Duration emmean SE df lower.CL upper.CL
## BL 1.12 0.413 28 0.2785 1.97
## X5 6.88 0.960 28 4.9082 8.84
## X10 9.25 1.197 28 6.7988 11.70
## X15 13.75 1.545 28 10.5861 16.91
## X20 19.00 1.511 28 15.9046 22.10
##
## UI_Condition = Coil:
## Duration emmean SE df lower.CL upper.CL
## BL 1.38 0.413 28 0.5285 2.22
## X5 1.88 0.960 28 -0.0918 3.84
## X10 3.62 1.197 28 1.1738 6.08
## X15 5.25 1.545 28 2.0861 8.41
## X20 10.00 1.511 28 6.9046 13.10
##
## Confidence level used: 0.95test(emm1.sme, joint=TRUE, adjust="none")## UI_Condition df1 df2 F.ratio p.value
## Control 5 28 176.937 <.0001
## Horizon 5 28 55.416 <.0001
## Dot 5 28 36.225 <.0001
## Coil 5 28 11.118 <.0001#paired t-tests
pairs(emm1.sme, adjust="tukey")## UI_Condition = Control:
## contrast estimate SE df t.ratio p.value
## BL - X5 -13.88 1.090 28 -12.727 <.0001
## BL - X10 -21.00 1.268 28 -16.559 <.0001
## BL - X15 -31.00 1.596 28 -19.423 <.0001
## BL - X20 -40.25 1.526 28 -26.373 <.0001
## X5 - X10 -7.12 0.911 28 -7.824 <.0001
## X5 - X15 -17.12 1.318 28 -12.993 <.0001
## X5 - X20 -26.38 1.501 28 -17.575 <.0001
## X10 - X15 -10.00 0.935 28 -10.694 <.0001
## X10 - X20 -19.25 1.197 28 -16.087 <.0001
## X15 - X20 -9.25 1.154 28 -8.015 <.0001
##
## UI_Condition = Horizon:
## contrast estimate SE df t.ratio p.value
## BL - X5 -6.62 1.090 28 -6.077 <.0001
## BL - X10 -9.38 1.268 28 -7.393 <.0001
## BL - X15 -16.38 1.596 28 -10.260 <.0001
## BL - X20 -20.75 1.526 28 -13.596 <.0001
## X5 - X10 -2.75 0.911 28 -3.020 0.0394
## X5 - X15 -9.75 1.318 28 -7.397 <.0001
## X5 - X20 -14.12 1.501 28 -9.412 <.0001
## X10 - X15 -7.00 0.935 28 -7.486 <.0001
## X10 - X20 -11.38 1.197 28 -9.506 <.0001
## X15 - X20 -4.38 1.154 28 -3.791 0.0061
##
## UI_Condition = Dot:
## contrast estimate SE df t.ratio p.value
## BL - X5 -5.75 1.090 28 -5.274 0.0001
## BL - X10 -8.12 1.268 28 -6.407 <.0001
## BL - X15 -12.62 1.596 28 -7.910 <.0001
## BL - X20 -17.88 1.526 28 -11.712 <.0001
## X5 - X10 -2.38 0.911 28 -2.608 0.0960
## X5 - X15 -6.88 1.318 28 -5.216 0.0001
## X5 - X20 -12.12 1.501 28 -8.079 <.0001
## X10 - X15 -4.50 0.935 28 -4.812 0.0004
## X10 - X20 -9.75 1.197 28 -8.148 <.0001
## X15 - X20 -5.25 1.154 28 -4.549 0.0008
##
## UI_Condition = Coil:
## contrast estimate SE df t.ratio p.value
## BL - X5 -0.50 1.090 28 -0.459 0.9904
## BL - X10 -2.25 1.268 28 -1.774 0.4078
## BL - X15 -3.88 1.596 28 -2.428 0.1374
## BL - X20 -8.62 1.526 28 -5.651 <.0001
## X5 - X10 -1.75 0.911 28 -1.922 0.3297
## X5 - X15 -3.38 1.318 28 -2.561 0.1058
## X5 - X20 -8.12 1.501 28 -5.414 0.0001
## X10 - X15 -1.62 0.935 28 -1.738 0.4284
## X10 - X20 -6.38 1.197 28 -5.327 0.0001
## X15 - X20 -4.75 1.154 28 -4.116 0.0026
##
## P value adjustment: tukey method for comparing a family of 5 estimates# Simple main effects of UI_Condition, the BG factor at levels of the repeated measure factor.
emm2.sme <- emmeans(fit3.afex, "UI_Condition", by="Duration")
emm2.sme## Duration = BL:
## UI_Condition emmean SE df lower.CL upper.CL
## Control 2.25 0.413 28 1.4035 3.10
## Horizon 2.00 0.413 28 1.1535 2.85
## Dot 1.12 0.413 28 0.2785 1.97
## Coil 1.38 0.413 28 0.5285 2.22
##
## Duration = X5:
## UI_Condition emmean SE df lower.CL upper.CL
## Control 16.12 0.960 28 14.1582 18.09
## Horizon 8.62 0.960 28 6.6582 10.59
## Dot 6.88 0.960 28 4.9082 8.84
## Coil 1.88 0.960 28 -0.0918 3.84
##
## Duration = X10:
## UI_Condition emmean SE df lower.CL upper.CL
## Control 23.25 1.197 28 20.7988 25.70
## Horizon 11.38 1.197 28 8.9238 13.83
## Dot 9.25 1.197 28 6.7988 11.70
## Coil 3.62 1.197 28 1.1738 6.08
##
## Duration = X15:
## UI_Condition emmean SE df lower.CL upper.CL
## Control 33.25 1.545 28 30.0861 36.41
## Horizon 18.38 1.545 28 15.2111 21.54
## Dot 13.75 1.545 28 10.5861 16.91
## Coil 5.25 1.545 28 2.0861 8.41
##
## Duration = X20:
## UI_Condition emmean SE df lower.CL upper.CL
## Control 42.50 1.511 28 39.4046 45.60
## Horizon 22.75 1.511 28 19.6546 25.85
## Dot 19.00 1.511 28 15.9046 22.10
## Coil 10.00 1.511 28 6.9046 13.10
##
## Confidence level used: 0.95test(contrast(emm2.sme), joint=TRUE, adjust="none")## Duration df1 df2 F.ratio p.value note
## BL 3 28 1.617 0.2078 d
## X5 3 28 37.831 <.0001 d
## X10 3 28 47.627 <.0001 d
## X15 3 28 57.687 <.0001 d
## X20 3 28 82.337 <.0001 d
##
## d: df1 reduced due to linear dependencepairs(emm2.sme)## Duration = BL:
## contrast estimate SE df t.ratio p.value
## Control - Horizon 0.250 0.584 28 0.428 0.9732
## Control - Dot 1.125 0.584 28 1.925 0.2408
## Control - Coil 0.875 0.584 28 1.497 0.4525
## Horizon - Dot 0.875 0.584 28 1.497 0.4525
## Horizon - Coil 0.625 0.584 28 1.069 0.7105
## Dot - Coil -0.250 0.584 28 -0.428 0.9732
##
## Duration = X5:
## contrast estimate SE df t.ratio p.value
## Control - Horizon 7.500 1.358 28 5.523 <.0001
## Control - Dot 9.250 1.358 28 6.812 <.0001
## Control - Coil 14.250 1.358 28 10.495 <.0001
## Horizon - Dot 1.750 1.358 28 1.289 0.5774
## Horizon - Coil 6.750 1.358 28 4.971 0.0002
## Dot - Coil 5.000 1.358 28 3.682 0.0051
##
## Duration = X10:
## contrast estimate SE df t.ratio p.value
## Control - Horizon 11.875 1.692 28 7.017 <.0001
## Control - Dot 14.000 1.692 28 8.273 <.0001
## Control - Coil 19.625 1.692 28 11.597 <.0001
## Horizon - Dot 2.125 1.692 28 1.256 0.5978
## Horizon - Coil 7.750 1.692 28 4.580 0.0005
## Dot - Coil 5.625 1.692 28 3.324 0.0125
##
## Duration = X15:
## contrast estimate SE df t.ratio p.value
## Control - Horizon 14.875 2.184 28 6.810 <.0001
## Control - Dot 19.500 2.184 28 8.927 <.0001
## Control - Coil 28.000 2.184 28 12.819 <.0001
## Horizon - Dot 4.625 2.184 28 2.117 0.1723
## Horizon - Coil 13.125 2.184 28 6.009 <.0001
## Dot - Coil 8.500 2.184 28 3.891 0.0030
##
## Duration = X20:
## contrast estimate SE df t.ratio p.value
## Control - Horizon 19.750 2.137 28 9.242 <.0001
## Control - Dot 23.500 2.137 28 10.996 <.0001
## Control - Coil 32.500 2.137 28 15.208 <.0001
## Horizon - Dot 3.750 2.137 28 1.755 0.3157
## Horizon - Coil 12.750 2.137 28 5.966 <.0001
## Dot - Coil 9.000 2.137 28 4.211 0.0013
##
## P value adjustment: tukey method for comparing a family of 4 estimates#Contrasts on simple main effects
test(contrast(emm1.sme, "poly"), adjust="none")## UI_Condition = Control:
## contrast estimate SE df t.ratio p.value
## linear 97.62 3.95 28 24.708 <.0001
## quadratic -6.38 3.51 28 -1.815 0.0803
## cubic 6.00 2.17 28 2.762 0.0100
## quartic -13.25 4.67 28 -2.838 0.0084
##
## UI_Condition = Horizon:
## contrast estimate SE df t.ratio p.value
## linear 51.25 3.95 28 12.971 <.0001
## quadratic -0.25 3.51 28 -0.071 0.9438
## cubic 1.25 2.17 28 0.575 0.5696
## quartic -15.00 4.67 28 -3.212 0.0033
##
## UI_Condition = Dot:
## contrast estimate SE df t.ratio p.value
## linear 42.62 3.95 28 10.788 <.0001
## quadratic 1.12 3.51 28 0.320 0.7512
## cubic 4.12 2.17 28 1.899 0.0680
## quartic -6.88 4.67 28 -1.472 0.1521
##
## UI_Condition = Coil:
## contrast estimate SE df t.ratio p.value
## linear 20.62 3.95 28 5.220 <.0001
## quadratic 8.38 3.51 28 2.384 0.0241
## cubic 1.88 2.17 28 0.863 0.3954
## quartic 4.62 4.67 28 0.990 0.3304#Simple Main Effects (Trend on duration separately for each UI Condition)
Control <- subset(Prompt_tbl, UI_Condition=="Control")
fitcontrol.afex <- aov_car(SSQ~Duration + Error(ID/Duration), data=Control,
anova_table = list(correction="none"))
anova(fitcontrol.afex, correction="none")## Anova Table (Type 3 tests)
##
## Response: SSQ
## num Df den Df MSE F ges Pr(>F)
## Duration 4 28 14.007 137.37 0.88123 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1control.emm <- emmeans(fitcontrol.afex, "Duration")
control.emm## Duration emmean SE df lower.CL upper.CL
## BL 2.25 0.25 7 1.66 2.84
## X5 16.12 1.55 7 12.46 19.79
## X10 23.25 1.97 7 18.59 27.91
## X15 33.25 2.52 7 27.29 39.21
## X20 42.50 2.41 7 36.79 48.21
##
## Confidence level used: 0.95test(contrast(control.emm, "poly"), adjust="none")## contrast estimate SE df t.ratio p.value
## linear 97.62 6.56 7 14.882 <.0001
## quadratic -6.38 4.02 7 -1.585 0.1570
## cubic 6.00 3.41 7 1.758 0.1221
## quartic -13.25 5.16 7 -2.567 0.0371###RSTATIX SIMPLE MAIN EFFECTS####
simple.UI <- Prompt_tbl %>%
group_by(Duration) %>%
anova_test(dv=SSQ, wid=ID, between=UI_Condition) %>%
get_anova_table() %>%
adjust_pvalue(method="holm")
simple.UI## # A tibble: 5 × 9
## Duration Effect DFn DFd F p `p<.05` ges p.adj
## <ord> <chr> <dbl> <dbl> <dbl> <dbl> <chr> <dbl> <dbl>
## 1 BL UI_Condition 3 28 1.62 2.08e- 1 "" 0.148 2.08e- 1
## 2 5 UI_Condition 3 28 37.8 5.53e-10 "*" 0.802 1.11e- 9
## 3 10 UI_Condition 3 28 47.6 4.01e-11 "*" 0.836 1.20e-10
## 4 15 UI_Condition 3 28 57.7 4.17e-12 "*" 0.861 1.67e-11
## 5 20 UI_Condition 3 28 82.3 5.30e-14 "*" 0.898 2.65e-13simple.Duration <- Prompt_tbl %>%
group_by(UI_Condition) %>%
anova_test(dv=SSQ, wid=ID, within=Duration) %>%
#get_anova_table(correction="GG") %>%
get_anova_table(correction="none") %>%
adjust_pvalue(method="holm")
simple.Duration## # A tibble: 4 × 9
## UI_Condition Effect DFn DFd F p `p<.05` ges p.adj
## <ord> <chr> <dbl> <dbl> <dbl> <dbl> <chr> <dbl> <dbl>
## 1 Control Duration 4 28 137. 5.68e-18 * 0.881 2.27e-17
## 2 Horizon Duration 4 28 88.4 1.84e-15 * 0.88 3.68e-15
## 3 Dot Duration 4 28 112. 8.27e-17 * 0.89 2.48e-16
## 4 Coil Duration 4 28 38.6 5.20e-11 * 0.794 5.20e-11