Dark
Jesi
2025-07-05
hypothesis: hot group will spend less time in the dark, avoiding the the side of the tank where they experinced a noxious stimulus. there may be learning effects over trials and/or experiments
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 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)])))Raw data plot:
library(ggplot2)
ggplot(data, aes(x = trial, y = dark, 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 = "Dark Across Trials by Treatment - Raw Data", y = "dark", x = "Trial", color="Hot?")
#Transformation not necessary - boxcox transformation did not make it significantly better and actually increased skewness. Use raw variable. #Model with dark data
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)## refitting model(s) with ML (instead of REML)## 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
anova(model)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## hot 148089 148089 1 23.499 4.8298 0.03808 *
## trial 129499 64749 2 105.000 2.1118 0.12614
## exp 57699 57699 1 105.000 1.8818 0.17305
## hot:trial 79712 39856 2 105.000 1.2999 0.27691
## hot:exp 200555 200555 1 21.000 6.5410 0.01834 *
## trial:exp 140767 70383 2 105.000 2.2955 0.10574
## hot:trial:exp 67206 33603 2 105.000 1.0959 0.33802
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#significant interactions for trial and experiment
# Get EMMs for all combinations
emm <- emmeans(model, ~ hot * trial * exp)
# Compare trials within hot and exp
trial_contrasts <- contrast(emm, method = "pairwise", by = c("hot", "exp"), adjust = "tukey")
summary(trial_contrasts)## hot = n, exp = 1:
## contrast estimate SE df t.ratio p.value
## trial1 - trial2 -81.7 71.5 105 -1.142 0.4903
## trial1 - trial3 -105.2 71.5 105 -1.472 0.3086
## trial2 - trial3 -23.5 71.5 105 -0.329 0.9420
##
## hot = y, exp = 1:
## contrast estimate SE df t.ratio p.value
## trial1 - trial2 -85.8 74.7 105 -1.149 0.4866
## trial1 - trial3 28.2 74.7 105 0.378 0.9244
## trial2 - trial3 114.0 74.7 105 1.526 0.2827
##
## hot = n, exp = 2:
## contrast estimate SE df t.ratio p.value
## trial1 - trial2 73.7 74.7 105 0.987 0.5865
## trial1 - trial3 103.6 74.7 105 1.388 0.3509
## trial2 - trial3 29.9 74.7 105 0.401 0.9154
##
## hot = y, exp = 2:
## contrast estimate SE df t.ratio p.value
## trial1 - trial2 59.9 71.5 105 0.837 0.6807
## trial1 - trial3 44.9 71.5 105 0.629 0.8048
## trial2 - trial3 -14.9 71.5 105 -0.209 0.9763
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 3 estimates# Compare hot vs. control within each trial × exp:
hot_contrasts <- contrast(emm, method = "pairwise", by = c("trial", "exp"), adjust = "tukey")
summary(hot_contrasts)## trial = 1, exp = 1:
## contrast estimate SE df t.ratio p.value
## n - y -111.8 78.4 116 -1.426 0.1566
##
## trial = 2, exp = 1:
## contrast estimate SE df t.ratio p.value
## n - y -115.9 78.4 116 -1.478 0.1421
##
## trial = 3, exp = 1:
## contrast estimate SE df t.ratio p.value
## n - y 21.6 78.4 116 0.275 0.7837
##
## trial = 1, exp = 2:
## contrast estimate SE df t.ratio p.value
## n - y 166.3 78.4 116 2.121 0.0361
##
## trial = 2, exp = 2:
## contrast estimate SE df t.ratio p.value
## n - y 152.5 78.4 116 1.944 0.0543
##
## trial = 3, exp = 2:
## contrast estimate SE df t.ratio p.value
## n - y 107.7 78.4 116 1.372 0.1726
##
## Degrees-of-freedom method: kenward-roger#no significant differences between trials
#Significant difference between treatment groups in experiment 2, trial 1 (trend in trial 2)- control spent more time in the dark
Looking across trials
#Trial 1 treatment model vs null
model <- lmer( dark ~ hot * exp +
(1 | subject),
data = T1)
nullmodel <- lmer( dark ~ 1 * exp +
(1 | subject),
data = T1)
anova(model, nullmodel)## refitting model(s) with ML (instead of REML)## Data: T1
## Models:
## nullmodel: dark ~ 1 * exp + (1 | subject)
## model: dark ~ hot * exp + (1 | subject)
## npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
## nullmodel 3 618.98 624.47 -306.49 612.98
## model 6 611.60 622.57 -299.80 599.60 13.387 3 0.00387 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(model)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## hot 156672 156672 1 25.1 5.6632 0.02524 *
## exp 190100 190100 1 21.0 6.8716 0.01595 *
## hot:exp 197001 197001 1 21.0 7.1210 0.01438 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#great improvement on the null model.
Significant effect of treatment, exp, and interaction in trial 1
#Trial 2 treatment model vs null
model <- lmer( dark ~ hot * exp +
(1 | subject),
data = T2)
nullmodel <- lmer(dark ~ 1 * exp +
(1 | subject),
data = T2)
anova(model, nullmodel)## refitting model(s) with ML (instead of REML)## Data: T2
## Models:
## nullmodel: dark ~ 1 * exp + (1 | subject)
## model: dark ~ hot * exp + (1 | subject)
## npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
## nullmodel 3 613.15 618.64 -303.57 607.15
## model 6 614.62 625.59 -301.31 602.62 4.5306 3 0.2096anova(model)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## hot 65377 65377 1 22.681 3.4695 0.07551 .
## exp 5458 5458 1 21.000 0.2897 0.59610
## hot:exp 74617 74617 1 21.000 3.9599 0.05977 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#does not imrpove upon null model #near significant effect of treatment and interaction in trial 2
#Trial 3 treatment model vs null
#interaction not
model <- lmer( dark ~ hot * exp +
(1 | subject),
data = T3)## boundary (singular) fit: see help('isSingular')nullmodel <- lmer( dark ~ 1 * exp +
(1 | subject),
data = T3)## boundary (singular) fit: see help('isSingular')anova(model, nullmodel)## refitting model(s) with ML (instead of REML)## Data: T3
## Models:
## nullmodel: dark ~ 1 * exp + (1 | subject)
## model: dark ~ hot * exp + (1 | subject)
## npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
## nullmodel 3 622.77 628.25 -308.38 616.77
## model 6 626.90 637.87 -307.45 614.90 1.8641 3 0.6011anova(model)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## hot 4774.5 4774.5 1 42 0.1165 0.7345
## exp 2907.5 2907.5 1 42 0.0710 0.7912
## hot:exp 21259.6 21259.6 1 42 0.5189 0.4753#does not improve on null model
#no significant differences in trial 3
Separated by Experiment
# Experiment 1
p1<-ggplot(exp1, aes(x = factor(trial), y = dark, 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 - Data + Means and SE" ,
y = "Dark",
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 = dark, 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 - Data + Means and SE",
y = "",
x = "Trial",
color = "Hot?") +
theme_minimal()
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")
Experiment 1 treatment model vs null
#interaction does not improve model in exp 1
model <- lmer(dark ~ hot * trial +
(1 | subject),
data = exp1)
nullmodel <- lmer(dark ~ 1 * trial +
(1 | subject),
data = exp1)
anova(model, nullmodel)## refitting model(s) with ML (instead of REML)## Data: exp1
## Models:
## nullmodel: dark ~ 1 * trial + (1 | subject)
## model: dark ~ hot * trial + (1 | subject)
## npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
## nullmodel 3 927.36 934.06 -460.68 921.36
## model 8 930.41 948.28 -457.21 914.41 6.9492 5 0.2244anova(model)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## hot 36205 36205 1 21 1.2987 0.2673
## trial 80606 40303 2 42 1.4457 0.2471
## hot:trial 70270 35135 2 42 1.2603 0.2941#no significant effect of treatment or trial in exp 1
# Fit your model
model1 <- lmer(dark ~ 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 285 57.3 53.7 170 400
## 2 366 57.3 53.7 251 481
## 3 390 57.3 53.7 275 505
##
## hot = y:
## trial emmean SE df lower.CL upper.CL
## 1 396 59.9 53.7 276 516
## 2 482 59.9 53.7 362 602
## 3 368 59.9 53.7 248 488
##
## 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 -81.7 68.2 42 -1.198 0.4609
## trial1 - trial3 -105.2 68.2 42 -1.544 0.2814
## trial2 - trial3 -23.5 68.2 42 -0.345 0.9364
##
## hot = y:
## contrast estimate SE df t.ratio p.value
## trial1 - trial2 -85.8 71.2 42 -1.205 0.4571
## trial1 - trial3 28.2 71.2 42 0.396 0.9172
## trial2 - trial3 114.0 71.2 42 1.601 0.2565
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 3 estimates#no significant differences in dark between trials within treatments in exp1
#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 -111.8 82.9 53.7 -1.349 0.1831
##
## trial = 2:
## contrast estimate SE df t.ratio p.value
## n - y -115.9 82.9 53.7 -1.398 0.1678
##
## trial = 3:
## contrast estimate SE df t.ratio p.value
## n - y 21.6 82.9 53.7 0.260 0.7957
##
## Degrees-of-freedom method: kenward-roger#no significant differences in dark between treatment groups within trials in exp 1
emm1 <- emmeans(model1, ~ hot * trial)
emm_df1 <- as.data.frame(emm1)
# Convert trial to numeric (if factor)
emm_df1 <- emm_df1 %>%
mutate(trial = as.numeric(as.factor(trial))) # ensures proper spacing on x-axis
# Base plot
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")) +
theme_minimal() +
labs(title="Dark in Experiment 1", x = "Trial", y = "Estimated Mean", color = "Treatment (hot)")
print(p1)
Experiment 2
# Experiment 2 treatment model vs null
model <- lmer( dark ~ hot * trial +
(1 | subject),
data = exp2)
nullmodel <- lmer( dark ~ 1 * trial +
(1 | subject),
data = exp2)
anova(model, nullmodel)## refitting model(s) with ML (instead of REML)## Data: exp2
## Models:
## nullmodel: dark ~ 1 * trial + (1 | subject)
## model: dark ~ hot * trial + (1 | subject)
## npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
## nullmodel 3 919.96 926.66 -456.98 913.96
## model 8 918.46 936.34 -451.23 902.46 11.493 5 0.04244 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(model)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## hot 244533 244533 1 21 8.9567 0.006935 **
## trial 76784 38392 2 42 1.4062 0.256369
## hot:trial 10803 5402 2 42 0.1978 0.821256
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#Model significantly better than null.
#significant treatment effect in experiment 2. Hot significantly more likely to spend less time in the dark. Non-significant but consistent decreases in time spent in the dark over trials 2 and 3.
model2 <- lmer(dark ~ 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 552 53.2 61.1 446 659
## 2 479 53.2 61.1 372 585
## 3 449 53.2 61.1 342 555
##
## hot = y:
## trial emmean SE df lower.CL upper.CL
## 1 386 51.0 61.1 284 488
## 2 326 51.0 61.1 224 428
## 3 341 51.0 61.1 239 443
##
## 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 73.7 70.5 42 1.046 0.5523
## trial1 - trial3 103.6 70.5 42 1.471 0.3150
## trial2 - trial3 29.9 70.5 42 0.425 0.9057
##
## hot = y:
## contrast estimate SE df t.ratio p.value
## trial1 - trial2 59.9 67.5 42 0.887 0.6511
## trial1 - trial3 44.9 67.5 42 0.666 0.7843
## trial2 - trial3 -14.9 67.5 42 -0.221 0.9734
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 3 estimates#no significant differences between trials within treatment groups
#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 166 73.7 61.1 2.258 0.0275
##
## trial = 2:
## contrast estimate SE df t.ratio p.value
## n - y 152 73.7 61.1 2.070 0.0427
##
## trial = 3:
## contrast estimate SE df t.ratio p.value
## n - y 108 73.7 61.1 1.461 0.1491
##
## Degrees-of-freedom method: kenward-roger#in experiment 2, significant difference between treatments in trials 1 and 2, and a weak trend in trial 3. Hot group spends less time in the dark across trials, significantly so in trials 1 and 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 <- 705 # 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()
# --- 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 = 675, label = "*", size = 6) + # trial 1
annotate("text", x = 2, y = 605, label = "*", size = 6) # trial 2
# --- Combine and Display ---
(p1 + p2) +
plot_layout(ncol = 2, guides = "collect") &
theme(legend.position = "bottom")
#asterisks show significant differences between treatment groups in trials 1 and 2
#Extras——————————————————————————————
Side effect?
modelside<- lmer(dark~ side + (1|subject), data=data)
summary(modelside)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: dark ~ side + (1 | subject)
## Data: data
##
## REML criterion at convergence: 1819.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.5471 -0.4395 0.2173 0.6470 1.7054
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 7627 87.33
## Residual 30920 175.84
## Number of obs: 138, groups: subject, 23
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 395.921 32.635 21.001 12.132 5.95e-11 ***
## sider 7.208 47.190 21.001 0.153 0.88
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## sider -0.692#no effect from which side of experimental tank was nearest the window
Time effect
library(hms)## Warning: package 'hms' was built under R version 4.3.3library(lubridate)## Warning: package 'lubridate' was built under R version 4.3.3##
## Attaching package: 'lubridate'## The following object is masked from 'package:hms':
##
## hms## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, union# Step 1: Add ":00" to make it HH:MM:SS
data$time_clean <- paste0(data$time, ":00")
# Step 2: Parse as hms time
data$time_parsed <- hms::as_hms(data$time_clean)
# Step 3: Extract decimal hour
data$hour_decimal <- hour(data$time_parsed) + minute(data$time_parsed)/60data$hour2 <- data$hour_decimal^2
model_time_quad <- lmer(dark ~ hour_decimal + hour2 + (1 | subject), data = data)
summary(model_time_quad)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: dark ~ hour_decimal + hour2 + (1 | subject)
## Data: data
##
## REML criterion at convergence: 1817.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.6292 -0.4295 0.2460 0.6712 1.6494
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 6732 82.05
## Residual 31098 176.35
## Number of obs: 138, groups: subject, 23
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -250.559 476.152 134.881 -0.526 0.600
## hour_decimal 101.110 75.156 134.606 1.345 0.181
## hour2 -3.781 2.881 133.932 -1.313 0.192
##
## Correlation of Fixed Effects:
## (Intr) hr_dcm
## hour_deciml -0.995
## hour2 0.982 -0.996#no linear or quadratic trend associated with time of day in regard to dark