TouchTask <- read.csv("C:/Users/gibbsj1/OneDrive - College of Charleston/Jesi CofC/Thesis/Behavioral Design/TouchTask/TouchTask.csv")
TouchTask$disturbance<- as.integer(TouchTask$disturbance)
df_summary <- TouchTask %>%
group_by(video, treat) %>%
summarise(mean_rating = mean(disturbance), .groups = "drop")
ggplot(df_summary, aes(x = treat, y = mean_rating)) +
geom_boxplot() +
theme_minimal() +
ylab("Mean Volunteer Rating per Video") +
xlab("Treatment")
# Get unique subjects per treatment
subjects_hot <- unique(TouchTask$subject[TouchTask$treat == "hot"])
subjects_control <- unique(TouchTask$subject[TouchTask$treat == "control"])
# Assign warm colors to hot, cool colors to control
warm_colors <- colorRampPalette(c("indianred", "red2", "firebrick","firebrick4"))(length(subjects_hot))
cool_colors <- colorRampPalette(c( "gray80", "grey", "gray60","gray40"))(length(subjects_control))
# Combine into one named vector of colors for the subject factor
subject_colors <- c(setNames(warm_colors, subjects_hot),
setNames(cool_colors, subjects_control))
df_summary <- TouchTask %>%
group_by(treat) %>%
summarise(
mean_rating = mean(disturbance),
se_rating = sd(disturbance) / sqrt(n()),
.groups = "drop"
)
# Plot raw points colored by subject + mean ± SE
p<- ggplot(TouchTask, aes(x = treat, y = disturbance, color = subject)) +
geom_jitter(width = 0.3, alpha = 0.7, size = 2) + #height=0 removes vertical jitter
scale_color_manual(values = subject_colors) +
geom_point(data = df_summary, aes(x = treat, y = mean_rating),
color = "black", size = 3, inherit.aes = FALSE) +
geom_errorbar(data = df_summary, aes(x = treat, ymin = mean_rating - se_rating, ymax = mean_rating + se_rating),
width = 0.2, color = "black", inherit.aes = FALSE) +
theme_minimal() +
ylab("Degree of Disturbance") +
xlab("Treatment") +
ggtitle("Volunteer Ratings by Subject with Mean ± SE") +
guides(color = guide_legend(title = "Subject"))
p + scale_y_continuous(breaks = 1:5) +
annotate("text", x = Inf, y = -Inf, label = "p = 1.24e-12",
hjust = 1.1, vjust = -0.5, size = 4) +
guides(color = "none")
# Fit a cumulative link mixed model
library(ordinal)
## Warning: package 'ordinal' was built under R version 4.3.3
##
## Attaching package: 'ordinal'
## The following object is masked from 'package:dplyr':
##
## slice
clmm_model <- clmm(
factor(disturbance) ~ treat +
(1 | subject) + (1 | volunteer) + (1 | video),
data = TouchTask
)
# Summary
summary(clmm_model)
## Cumulative Link Mixed Model fitted with the Laplace approximation
##
## formula: factor(disturbance) ~ treat + (1 | subject) + (1 | volunteer) +
## (1 | video)
## data: TouchTask
##
## link threshold nobs logLik AIC niter max.grad cond.H
## logit flexible 252 -260.62 537.24 382(2420) 1.79e-05 1.8e+02
##
## Random effects:
## Groups Name Variance Std.Dev.
## video (Intercept) 5.0507 2.247
## volunteer (Intercept) 0.2353 0.485
## subject (Intercept) 0.1340 0.366
## Number of groups: video 84, volunteer 21, subject 12
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## treathot 6.0610 0.8536 7.1 1.24e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Threshold coefficients:
## Estimate Std. Error z value
## 1|2 1.2454 0.5149 2.419
## 2|3 4.5583 0.6909 6.598
## 3|4 5.2453 0.7289 7.196
## 4|5 8.2735 0.9305 8.891
Was the iron turned on?
# Convert 'y'/'n' to 1/0
TouchTask$on_bin <- ifelse(TouchTask$on == "y", 1, 0)
library(dplyr)
# Barplot with 95% CI Error Bars
TouchTask %>%
group_by(treat) %>%
summarise(
prop = mean(on_bin),
n = n(),
se = sqrt(prop * (1 - prop) / n),
lower = prop - 1.96 * se,
upper = prop + 1.96 * se
) %>%
ggplot(aes(x = treat, y = prop, fill = treat)) +
geom_col() +
geom_errorbar(aes(ymin = lower, ymax = upper), width = 0.2) +
scale_fill_manual(values = c("control" = "gray", "hot" = "red3")) +
ylim(0, 1) +
ylab("Proportion Volunteers Responding 'Hot'") +
xlab("Treatment") +
theme_minimal() +
theme(legend.position = "none")
library(ggplot2)
library(dplyr)
df_summary <- TouchTask %>%
group_by(treat) %>%
summarise(
prop = mean(on_bin),
n = n(),
se = sqrt(prop * (1 - prop) / n),
lower = prop - 1.96 * se,
upper = prop + 1.96 * se,
.groups = "drop"
)
p<-ggplot(TouchTask, aes(x = treat, y = on_bin, color = treat)) +
geom_jitter(width = 0.4, height = 0, size = 2, alpha = 0.6) + # raw binary responses
geom_point(data = df_summary, aes(x = treat, y = prop, color = treat),
size = 4, shape = 18, inherit.aes = FALSE) + # proportion point
geom_errorbar(data = df_summary, aes(x = treat, ymin = lower, ymax = upper, color = treat),
width = 0.2, size = 1, inherit.aes = FALSE) + # 95% CI bars
scale_color_manual(values = c("control" = "gray60", "hot" = "red2")) +
scale_y_continuous(breaks = c(0, 1), labels = c("No", "Yes"), limits = c(-0.05, 1.05)) +
theme_minimal() +
theme(legend.position = "none") +
ylab("Volunteer Response") +
xlab("Treatment") +
ggtitle(" Was the Iron Hot? Responses, Mean ± 95% CI")
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
p+ annotate("text", x = Inf, y = -Inf, label = "p = 5.34e-07", hjust = 1.1, vjust = -0.5, size = 4)
library(Matrix)
##
## Attaching package: 'Matrix'
## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack
library(lme4)
# Fit the GLMM with binomial link
glmm_on <- glmer(
on_bin ~ treat +
(1 | subject) + (1 | volunteer) + (1 | video),
data = TouchTask,
family = binomial(link = "logit")
)
# View results
summary(glmm_on)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: on_bin ~ treat + (1 | subject) + (1 | volunteer) + (1 | video)
## Data: TouchTask
##
## AIC BIC logLik -2*log(L) df.resid
## 220.6 238.2 -105.3 210.6 247
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.1105 -0.2199 -0.1481 0.4280 3.0074
##
## Random effects:
## Groups Name Variance Std.Dev.
## video (Intercept) 1.37616 1.1731
## volunteer (Intercept) 0.67941 0.8243
## subject (Intercept) 0.05052 0.2248
## Number of obs: 252, groups: video, 84; volunteer, 21; subject, 12
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.3681 0.7129 -4.724 2.31e-06 ***
## treathot 4.6560 0.9287 5.014 5.34e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## treathot -0.900