Middle and PauseMID
Jesi
2025-07-05
data<-Behavioral.Trials.All.Data#Separate trials
T1<- data[data$trial=="1",]
T2<- data[data$trial=="2",]
T3<- data[data$trial=="3",]
data23<- rbind(T2,T3)
#separate experiments
exp1<-(data[data$exp=="1",])
exp2<-(data[data$exp=="2",])
#separate treatment groups
hot<-(data[data$hot=="y",])
control<-(data[data$hot=="n",])# Reorder subjects based on the 'hot' group
data <- data %>%
mutate(subject = factor(subject, levels = unique(subject[order(hot)])))Middle, pauseMiddle
#Plot raw values
p1<- ggplot(data, aes(x = trial, y = middle, color = hot)) +
stat_summary(fun = mean, geom = "line", aes(group = hot)) +
stat_summary(fun.data = mean_se, geom = "errorbar", width = 0.2) +
stat_summary(fun = mean, geom = "point", size = 3) +
scale_color_manual(values = c("n" = "grey50", "y" = "red")) +
labs(title = "Middle Across Trials by Treatment - Raw Data", y = "Light", x = "Trial", color="Hot?")
p2<- ggplot(data, aes(x = trial, y = pauseMID, color = hot)) +
stat_summary(fun = mean, geom = "line", aes(group = hot)) +
stat_summary(fun.data = mean_se, geom = "errorbar", width = 0.2) +
stat_summary(fun = mean, geom = "point", size = 3) +
scale_color_manual(values = c("n" = "grey50", "y" = "red")) +
labs(title = "Middle Pause Across Trials by Treatment - Raw Data", y = "Light", x = "Trial", color="Hot?")
library(patchwork)## Warning: package 'patchwork' was built under R version 4.3.3# Combine plots and collect shared legend
(p1 + p2) +
plot_layout(ncol = 2, guides = "collect") &
theme(legend.position = "bottom")
#plots nearly identical
cor.test(data$middle, data$pauseMID, method = "spearman")## Warning in cor.test.default(data$middle, data$pauseMID, method = "spearman"):
## Cannot compute exact p-value with ties##
## Spearman's rank correlation rho
##
## data: data$middle and data$pauseMID
## S = 62884, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.8564251#strong correlation between time paused in the middle and pause middle - Not a large area, if not paused, only spent a couple of seconds in middle.
#Just going to focus on middle because the transformed histogram looks more normal.
#Determine Transformation
#make new column with no 0 values
min_val <- min(data$middle)
data$mid_shift <- data$middle + abs(min_val) + 1
powerTransform(data$middle + abs(min(data$middle)) + 1)## Estimated transformation parameter
## data$middle + abs(min(data$middle)) + 1
## 0.07227417#Perform transformation
lambda <- 0.07227417
data$mid_trans <- data$mid_shift ^ lambdahist(data$middle)
hist(data$mid_trans)
#Great improvement via transformation
Modelling
model <- lmer( dark ~ hot * trial * exp + (1 | subject), data = data)
modeltrial <- lmer( dark ~ hot * trial + (1 | subject), data = data)
modelexp <- lmer( dark ~ hot * exp + (1 | subject), data = data)
anova(model, modeltrial, modelexp)## Data: data
## Models:
## modelexp: dark ~ hot * exp + (1 | subject)
## modeltrial: dark ~ hot * trial + (1 | subject)
## model: dark ~ hot * trial * exp + (1 | subject)
## npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
## modelexp 6 1839.8 1857.3 -913.87 1827.8
## modeltrial 8 1851.0 1874.4 -917.47 1835.0 0.000 2 1.00000
## model 14 1847.6 1888.6 -909.80 1819.6 15.355 6 0.01767 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#model including trial and experiment interactions significantly improves fit
model <- lmer( mid_trans ~ hot * trial * exp +
(1 | subject),
data = data)
nullmodel <- lmer( mid_trans ~ 1 * trial * exp +
(1 | subject),
data = data)
anova(model, nullmodel)## Data: data
## Models:
## nullmodel: mid_trans ~ 1 * trial * exp + (1 | subject)
## model: mid_trans ~ hot * trial * exp + (1 | subject)
## npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
## nullmodel 3 -81.168 -72.386 43.584 -87.168
## model 14 -89.812 -48.830 58.906 -117.812 30.644 11 0.001254 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(model)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: mid_trans ~ hot * trial * exp + (1 | subject)
## Data: data
##
## REML criterion at convergence: -66.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.2258 -0.6550 -0.0494 0.7233 2.0608
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 0.001572 0.03965
## Residual 0.025931 0.16103
## Number of obs: 138, groups: subject, 23
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.58717 0.10802 123.97428 14.693 < 2e-16 ***
## hoty -0.45762 0.15829 97.34183 -2.891 0.004737 **
## trial2 -0.22380 0.14833 105.00000 -1.509 0.134362
## trial3 -0.23407 0.14833 105.00000 -1.578 0.117564
## exp -0.23504 0.06923 123.97428 -3.395 0.000922 ***
## hoty:trial2 0.26136 0.21256 105.00000 1.230 0.221613
## hoty:trial3 0.36523 0.21256 105.00000 1.718 0.088700 .
## hoty:exp 0.36274 0.10066 89.33206 3.604 0.000516 ***
## trial2:exp 0.18581 0.09506 105.00000 1.955 0.053280 .
## trial3:exp 0.19823 0.09506 105.00000 2.085 0.039463 *
## hoty:trial2:exp -0.31080 0.13444 105.00000 -2.312 0.022734 *
## hoty:trial3:exp -0.27726 0.13444 105.00000 -2.062 0.041636 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) hoty trial2 trial3 exp hty:t2 hty:t3 hty:xp trl2:x
## hoty -0.714
## trial2 -0.687 0.469
## trial3 -0.687 0.469 0.500
## exp -0.947 0.685 0.650 0.650
## hoty:trial2 0.479 -0.671 -0.698 -0.349 -0.454
## hoty:trial3 0.479 -0.671 -0.349 -0.698 -0.454 0.500
## hoty:exp 0.689 -0.954 -0.447 -0.447 -0.727 0.634 0.634
## trial2:exp 0.650 -0.444 -0.947 -0.474 -0.687 0.661 0.331 0.472
## trial3:exp 0.650 -0.444 -0.474 -0.947 -0.687 0.331 0.661 0.472 0.500
## hty:trl2:xp -0.460 0.637 0.670 0.335 0.485 -0.949 -0.474 -0.668 -0.707
## hty:trl3:xp -0.460 0.637 0.335 0.670 0.485 -0.474 -0.949 -0.668 -0.354
## trl3:x hty:2:
## hoty
## trial2
## trial3
## exp
## hoty:trial2
## hoty:trial3
## hoty:exp
## trial2:exp
## trial3:exp
## hty:trl2:xp -0.354
## hty:trl3:xp -0.707 0.500#significant interactions for trial and experiment
#treatment significantly improves model
#Look within groups to compare trials to each other - trial 2 looks different from 1 and 3 for hot group
# Relevel to set trial 2 as the reference
hot$trial <- relevel(hot$trial, ref = "2")
#run model
model_hot <- lmer(mid_trans ~ trial + (1 | subject), data = hot)
# Summary gives direct pairwise tests vs. trial 2
summary(model_hot)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: mid_trans ~ trial + (1 | subject)
## Data: hot
##
## REML criterion at convergence: -46.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.06487 -0.64857 -0.03543 0.80229 1.66743
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 0.002388 0.04886
## Residual 0.022804 0.15101
## Number of obs: 69, groups: subject, 23
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.17123 0.03310 64.83520 35.390 < 2e-16 ***
## trial1 0.15264 0.04453 44.00000 3.428 0.001331 **
## trial3 0.16354 0.04453 44.00000 3.673 0.000647 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) trial1
## trial1 -0.673
## trial3 -0.673 0.500#In the hot group, time in the middle in trial 2 significantly differs from trials 1 and 3 (as well as from control (see above))
#Trial 1 > Trial 2 (p = 0.0013) #Trial 3 > Trial 2 (p = 0.00065)
#Look within groups to compare trials to each other
#run model
model_control <- lmer(mid_trans ~ trial + (1 | subject), data = control)
# Summary gives direct pairwise tests vs. trial 2
summary(model_control)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: mid_trans ~ trial + (1 | subject)
## Data: control
##
## REML criterion at convergence: -27.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.79891 -0.72875 0.03641 0.53402 1.85557
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 0.004975 0.07053
## Residual 0.029096 0.17058
## Number of obs: 69, groups: subject, 23
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.23972 0.03849 63.30097 32.211 <2e-16 ***
## trial2 0.05088 0.05030 44.00000 1.012 0.317
## trial3 0.05897 0.05030 44.00000 1.172 0.247
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) trial2
## trial2 -0.653
## trial3 -0.653 0.500#no differences between trials for the control group. Non-significant increase over trials.
By experiment
# Raw data
# Experiment 1
p1<-ggplot(exp1, aes(x = factor(trial), y = middle, color = hot, group = hot)) +
geom_jitter(width = 0.2, alpha = 0.5, size = 1.8) + # raw data
stat_summary(fun = mean, geom = "line", aes(group = hot), size = 1) + # mean line
stat_summary(fun.data = mean_se, geom = "errorbar", width = 0.2) + # error bars
stat_summary(fun = mean, geom = "point", size = 3, shape = 18) + # mean points
scale_color_manual(values = c("n" = "grey50", "y" = "red")) +
labs(title = "Exp 1 - Raw Data + Means and SE" ,
y = "Middle",
x = "Trial",
color = "Hot?") +
theme_minimal()## 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.# Experiment 2
p2<-ggplot(exp2, aes(x = factor(trial), y = middle, color = hot, group = hot)) +
geom_jitter(width = 0.2, alpha = 0.5, size = 1.8) +
stat_summary(fun = mean, geom = "line", aes(group = hot), size = 1) +
stat_summary(fun.data = mean_se, geom = "errorbar", width = 0.2) +
stat_summary(fun = mean, geom = "point", size = 3, shape = 18) +
scale_color_manual(values = c("n" = "grey50", "y" = "red")) +
labs(title = "Exp 2 - Raw Data + Means and SE",
y = "",
x = "Trial",
color = "Hot?") +
theme_minimal()
library(patchwork)
# Combine plots and collect shared legend
(p1 + p2) +
plot_layout(ncol = 2, guides = "collect") &
theme(legend.position = "bottom")
#Transformed data
# Get the range of mid_trans from both datasets
combined_range <- range(c(exp1$mid_trans, exp2$mid_trans), na.rm = TRUE)
# Updated plots with consistent y-axis limits
p1t <- ggplot(exp1, aes(x = trial, y = mid_trans, color = hot)) +
stat_summary(fun = mean, geom = "line", aes(group = hot)) +
stat_summary(fun.data = mean_se, geom = "errorbar", width = 0.2) +
stat_summary(fun = mean, geom = "point", size = 3) +
scale_color_manual(values = c("n" = "grey50", "y" = "red")) +
coord_cartesian(ylim = combined_range) +
labs(title = "Exp 1 - Transformed Data + Means, SE", y = "Light", x = "Trial", color = "Hot?")
p2t <- ggplot(exp2, aes(x = trial, y = mid_trans, color = hot)) +
stat_summary(fun = mean, geom = "line", aes(group = hot)) +
stat_summary(fun.data = mean_se, geom = "errorbar", width = 0.2) +
stat_summary(fun = mean, geom = "point", size = 3) +
scale_color_manual(values = c("n" = "grey50", "y" = "red")) +
coord_cartesian(ylim = combined_range) +
labs(title = "Exp 2 - Transformed Data + Means, SE", y = "Light", x = "Trial", color = "Hot?")
# Combine plots and collect shared legend
(p1t + p2t) +
plot_layout(ncol = 2, guides = "collect") &
theme(legend.position = "bottom")
Experiment 1
# Fit your model
model1 <- lmer(mid_trans ~ hot * trial + (1 | subject), data = exp1)
# Get estimated marginal means for each trial and hot combination
emm1 <- emmeans(model1, ~ trial | hot)
summary(emm1)## hot = n:
## trial emmean SE df lower.CL upper.CL
## 1 1.35 0.0496 57.7 1.25 1.45
## 2 1.31 0.0496 57.7 1.21 1.41
## 3 1.32 0.0496 57.7 1.22 1.42
##
## hot = y:
## trial emmean SE df lower.CL upper.CL
## 1 1.26 0.0518 57.7 1.15 1.36
## 2 1.17 0.0518 57.7 1.07 1.27
## 3 1.31 0.0518 57.7 1.21 1.41
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95#Compare Trials within groups
pairs(emm1, by = "hot")## hot = n:
## contrast estimate SE df t.ratio p.value
## trial1 - trial2 0.03798 0.0622 42 0.611 0.8151
## trial1 - trial3 0.03584 0.0622 42 0.576 0.8335
## trial2 - trial3 -0.00214 0.0622 42 -0.034 0.9993
##
## hot = y:
## contrast estimate SE df t.ratio p.value
## trial1 - trial2 0.08743 0.0650 42 1.346 0.3783
## trial1 - trial3 -0.05213 0.0650 42 -0.802 0.7036
## trial2 - trial3 -0.13956 0.0650 42 -2.148 0.0923
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 3 estimates#No significant differences between trials within groups
#Hot vs Control at each trial
emm1.1 <- emmeans(model1, ~ hot | trial)
pairs(emm1.1, by = "trial")## trial = 1:
## contrast estimate SE df t.ratio p.value
## n - y 0.09488 0.0718 57.7 1.322 0.1913
##
## trial = 2:
## contrast estimate SE df t.ratio p.value
## n - y 0.14433 0.0718 57.7 2.011 0.0490
##
## trial = 3:
## contrast estimate SE df t.ratio p.value
## n - y 0.00691 0.0718 57.7 0.096 0.9236
##
## Degrees-of-freedom method: kenward-roger#significant difference between treatment groups in trial 2 in experiment 1
Experiment 2
model2 <- lmer(mid_trans ~ hot * trial + (1 | subject), data = exp2)
# Get estimated marginal means for each trial and hot combination
emm2 <- emmeans(model2, ~ trial | hot)
summary(emm2)## hot = n:
## trial emmean SE df lower.CL upper.CL
## 1 1.12 0.0481 63 1.02 1.21
## 2 1.26 0.0481 63 1.17 1.36
## 3 1.28 0.0481 63 1.18 1.38
##
## hot = y:
## trial emmean SE df lower.CL upper.CL
## 1 1.38 0.0461 63 1.29 1.48
## 2 1.17 0.0461 63 1.08 1.26
## 3 1.36 0.0461 63 1.27 1.45
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95#Compare Trials within groups
pairs(emm2, by = "hot")## hot = n:
## contrast estimate SE df t.ratio p.value
## trial1 - trial2 -0.1478 0.0680 42 -2.173 0.0877
## trial1 - trial3 -0.1624 0.0680 42 -2.387 0.0550
## trial2 - trial3 -0.0146 0.0680 42 -0.214 0.9751
##
## hot = y:
## contrast estimate SE df t.ratio p.value
## trial1 - trial2 0.2124 0.0651 42 3.261 0.0061
## trial1 - trial3 0.0269 0.0651 42 0.413 0.9105
## trial2 - trial3 -0.1855 0.0651 42 -2.848 0.0182
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 3 estimates#Control group: no significant differences between trials. Trend toward more time in the middle between 1 and 3.
#Hot group: significantly less time in trial 2 than in trials 1 and 3. No difference between 1 and 3.
#Hot vs Control at each trial
emm2.2 <- emmeans(model2, ~ hot | trial)
pairs(emm2.2, by = "trial")## trial = 1:
## contrast estimate SE df t.ratio p.value
## n - y -0.2679 0.0666 63 -4.022 0.0002
##
## trial = 2:
## contrast estimate SE df t.ratio p.value
## n - y 0.0924 0.0666 63 1.387 0.1702
##
## trial = 3:
## contrast estimate SE df t.ratio p.value
## n - y -0.0786 0.0666 63 -1.180 0.2426
##
## Degrees-of-freedom method: kenward-roger#significant difference between treatments in trial 1 exp 2
# --- Estimated Marginal Means ---
# Experiment 1
emm1 <- emmeans(model1, ~ hot * trial)
emm_df1 <- as.data.frame(emm1)
emm_df1 <- emm_df1 %>%
mutate(trial = as.numeric(as.factor(trial)))
# Experiment 2
emm2 <- emmeans(model2, ~ hot * trial)
emm_df2 <- as.data.frame(emm2)
emm_df2 <- emm_df2 %>%
mutate(trial = as.numeric(as.factor(trial)))
# --- Compute shared y-axis limit that includes asterisks ---
asterisk_height <- 1.65 # maximum asterisk height
y_max <- max(c(emm_df1$upper.CL, emm_df2$upper.CL, asterisk_height)) + 0.05 * asterisk_height
y_min <- min(c(emm_df1$lower.CL, emm_df2$lower.CL))
# --- Plot for Experiment 1 ---
p1 <- ggplot(emm_df1, aes(x = factor(trial), y = emmean, group = hot)) +
geom_line(aes(color = hot), size = 1.2) +
geom_point(aes(color = hot), size = 3) +
geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL, color = hot), width = 0.2) +
scale_color_manual(values = c("n" = "grey50", "y" = "red")) +
labs(title = "Dark in Experiment 1", x = "Trial", y = "Estimated Mean", color = "Treatment (hot)") +
coord_cartesian(ylim = c(y_min, y_max)) +
theme_minimal()
# --- Annotate Experiment 1 with manual asterisks ---
p1 <- p1 +
annotate("text", x = 2, y = 1.45, label = "*", size = 6) # trial 2
# --- Plot for Experiment 2 ---
p2 <- ggplot(emm_df2, aes(x = factor(trial), y = emmean, group = hot)) +
geom_line(aes(color = hot), size = 1.2) +
geom_point(aes(color = hot), size = 3) +
geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL, color = hot), width = 0.2) +
scale_color_manual(values = c("n" = "grey50", "y" = "red")) +
labs(title = "Dark in Experiment 2", x = "Trial", y = "", color = "Treatment (hot)") +
coord_cartesian(ylim = c(y_min, y_max)) +
theme_minimal()
# --- Annotate Experiment 2 with manual asterisks ---
p2 <- p2 +
annotate("text", x = 1, y = 1.5, label = "*", size = 6) # trial 1
# Add red "a", "a,b", "b" for hot group
hot_df2 <- emm_df2 %>% filter(hot == "y")
p2 <- p2 +
annotate("text", x = 1, y = hot_df2$emmean[1] + 0.153, label = "a", color = "red", size = 3) +
annotate("text", x = 2, y = hot_df2$emmean[2] + -0.15, label = "a,b", color = "red", size = 3) +
annotate("text", x = 3, y = hot_df2$emmean[3] + 0.15, label = "b", color = "red", size = 3)
# --- Combine and Display ---
(p1 + p2) +
plot_layout(ncol = 2, guides = "collect") &
theme(legend.position = "bottom")
##asterisks show significant differences between treatment groups, letters show significant differences between trials within the hot treatment group