Light data
2025-06-24
hypothesis: hot group will spend more time in the light, avoiding the dark 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",])
hot<-(data[data$hot=="n",])Raw data plot:
library(ggplot2)
ggplot(data, aes(x = trial, y = light, 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 = "Light Across Trials by Treatment - Raw Data", y = "Light", x = "Trial", color="Hot?")
#Light variable not normal - does not meet model assumptions
describe(data$light)## vars n mean sd median trimmed mad min max range skew kurtosis
## X1 1 138 96.01 133.46 44 70.38 65.24 0 600.03 600.03 2.04 4.37
## se
## X1 11.36hist(data$light)
#Tried several transformations. This one is the best:
#Make new column with no 0 values
min_val <- min(data$light)
data$light_shift <- data$light + abs(min_val) + 1
#Transform
lambda <- 0.1262881
data$light_trans <- data$light_shift ^ lambda
describe(data$light_trans)## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 138 1.54 0.4 1.62 1.53 0.47 1 2.24 1.24 -0.14 -1.28 0.03hist(data$light_trans)
#much improved. #Model with transformed light data
modelhot<- lmer( light_trans ~ hot + (1 | subject), data = data)
model <- lmer( light_trans ~ hot * trial * exp + (1 | subject), data = data)
modeltrial <- lmer( light_trans ~ hot * trial + (1 | subject), data = data)
modelexp <- lmer( light_trans ~ hot * exp + (1 | subject), data = data)
anova(model, modeltrial, modelexp, modelhot)## refitting model(s) with ML (instead of REML)## Data: data
## Models:
## modelhot: light_trans ~ hot + (1 | subject)
## modelexp: light_trans ~ hot * exp + (1 | subject)
## modeltrial: light_trans ~ hot * trial + (1 | subject)
## model: light_trans ~ hot * trial * exp + (1 | subject)
## npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
## modelhot 4 114.32 126.03 -53.159 106.32
## modelexp 6 115.22 132.79 -51.612 103.22 3.0933 2 0.2130
## modeltrial 8 118.79 142.21 -51.397 102.79 0.4313 2 0.8060
## model 14 116.93 157.91 -44.465 88.93 13.8632 6 0.0312 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#including interaction terms for trial and experiment significantly improve model.
anova(model)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## hot 0.11141 0.11141 1 22.138 1.1559 0.29390
## trial 0.46891 0.23445 2 105.000 2.4324 0.09276 .
## exp 0.06713 0.06713 1 105.000 0.6964 0.40588
## hot:trial 0.22775 0.11388 2 105.000 1.1814 0.31088
## hot:exp 0.22474 0.22474 1 21.000 2.3317 0.14169
## trial:exp 0.68717 0.34358 2 105.000 3.5646 0.03178 *
## hot:trial:exp 0.30441 0.15220 2 105.000 1.5791 0.21102
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#significant interactions for trial and experiment, significant effect of 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 0.11585 0.127 105 0.914 0.6326
## trial1 - trial3 0.21812 0.127 105 1.721 0.2022
## trial2 - trial3 0.10227 0.127 105 0.807 0.6996
##
## hot = y, exp = 1:
## contrast estimate SE df t.ratio p.value
## trial1 - trial2 0.00695 0.132 105 0.052 0.9985
## trial1 - trial3 0.03067 0.132 105 0.232 0.9708
## trial2 - trial3 0.02373 0.132 105 0.179 0.9824
##
## hot = n, exp = 2:
## contrast estimate SE df t.ratio p.value
## trial1 - trial2 -0.21978 0.132 105 -1.660 0.2254
## trial1 - trial3 -0.26916 0.132 105 -2.033 0.1093
## trial2 - trial3 -0.04938 0.132 105 -0.373 0.9262
##
## hot = y, exp = 2:
## contrast estimate SE df t.ratio p.value
## trial1 - trial2 -0.29086 0.127 105 -2.295 0.0609
## trial1 - trial3 -0.04014 0.127 105 -0.317 0.9462
## trial2 - trial3 0.25072 0.127 105 1.978 0.1227
##
## 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 0.1219 0.16 79.2 0.762 0.4485
##
## trial = 2, exp = 1:
## contrast estimate SE df t.ratio p.value
## n - y 0.0130 0.16 79.2 0.081 0.9357
##
## trial = 3, exp = 1:
## contrast estimate SE df t.ratio p.value
## n - y -0.0656 0.16 79.2 -0.410 0.6829
##
## trial = 1, exp = 2:
## contrast estimate SE df t.ratio p.value
## n - y -0.3584 0.16 79.2 -2.241 0.0279
##
## trial = 2, exp = 2:
## contrast estimate SE df t.ratio p.value
## n - y -0.4295 0.16 79.2 -2.685 0.0088
##
## trial = 3, exp = 2:
## contrast estimate SE df t.ratio p.value
## n - y -0.1294 0.16 79.2 -0.809 0.4209
##
## Degrees-of-freedom method: kenward-roger#a lot of tests. Tukey applied to each set of comparisons.
#no significant differences between trials in either group.
#significant difference between treatment groups in exp2 trials 1 and 2
Looking across trials
#Trial 1 treatment model vs null
model <- lmer( light_trans ~ hot * exp +
(1 | subject),
data = T1)
nullmodel <- lmer( light_trans ~ 1 * exp +
(1 | subject),
data = T1)
anova(model, nullmodel)## refitting model(s) with ML (instead of REML)## Data: T1
## Models:
## nullmodel: light_trans ~ 1 * exp + (1 | subject)
## model: light_trans ~ hot * exp + (1 | subject)
## npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
## nullmodel 3 51.748 57.234 -22.874 45.748
## model 6 45.666 56.638 -16.833 33.666 12.082 3 0.007107 **
## ---
## 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 0.16135 0.16135 1 22.691 2.0089 0.169967
## exp 0.67634 0.67634 1 21.000 8.4209 0.008528 **
## hot:exp 0.24028 0.24028 1 21.000 2.9916 0.098370 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#Trial 2 treatment model vs null
model <- lmer( light_trans ~ hot * exp +
(1 | subject),
data = T2)## boundary (singular) fit: see help('isSingular')nullmodel <- lmer( light_trans ~ 1 * exp +
(1 | subject),
data = T2)
anova(model, nullmodel)## refitting model(s) with ML (instead of REML)## Data: T2
## Models:
## nullmodel: light_trans ~ 1 * exp + (1 | subject)
## model: light_trans ~ hot * exp + (1 | subject)
## npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
## nullmodel 3 52.692 58.178 -23.346 46.692
## model 6 51.039 62.011 -19.520 39.039 7.653 3 0.05376 .
## ---
## 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 0.23809 0.23809 1 42 1.5890 0.21442
## exp 0.06282 0.06282 1 42 0.4193 0.52083
## hot:exp 0.56184 0.56184 1 42 3.7498 0.05956 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#Treatment nearly significantly improves the model for T2
#near-significant result of hot*exp in trial 2
Trial 3 treatment model vs null
model <- lmer( light_trans ~ hot * exp +
(1 | subject),
data = T3)
nullmodel <- lmer(light_trans ~ 1 * exp +
(1 | subject),
data = T3)
anova(model, nullmodel)## refitting model(s) with ML (instead of REML)## Data: T3
## Models:
## nullmodel: light_trans ~ 1 * exp + (1 | subject)
## model: light_trans ~ hot * exp + (1 | subject)
## npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
## nullmodel 3 35.917 41.403 -14.958 29.917
## model 6 39.748 50.719 -13.874 27.747 2.1691 3 0.5381anova(model)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## hot 0.0000011 0.0000011 1 22.366 0.0000 0.9967
## exp 0.0151299 0.0151299 1 21.000 0.2384 0.6304
## hot:exp 0.0034275 0.0034275 1 21.000 0.0540 0.8185#Model not significantly better than null
#no effect of treatment, experiment, or interaction in trial 3
Separated by Experiment
# Raw data
# Experiment 1
p1<-ggplot(exp1, aes(x = factor(trial), y = light, 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 = "Light",
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 = light, 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 = "Light",
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")
#Transformed data
# Get the range of light_trans from both datasets
combined_range <- range(c(exp1$light_trans, exp2$light_trans), na.rm = TRUE)
# Updated plots with consistent y-axis limits
p1t <- ggplot(exp1, aes(x = trial, y = light_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 = light_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")
#is there any point to showing transformed data?
Experiment 1 treatment model vs null
model1 <- lmer(light_trans ~ hot * trial +
(1 | subject),
data = exp1)
nullmodel1 <- lmer(light_trans ~ 1 * trial +
(1 | subject),
data = exp1)
anova(model1, nullmodel1)## refitting model(s) with ML (instead of REML)## Data: exp1
## Models:
## nullmodel1: light_trans ~ 1 * trial + (1 | subject)
## model1: light_trans ~ hot * trial + (1 | subject)
## npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
## nullmodel1 3 69.890 76.592 -31.945 63.890
## model1 8 77.057 94.930 -30.528 61.057 2.8329 5 0.7257anova(model1)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## hot 0.003308 0.003308 1 21 0.0299 0.8643
## trial 0.177635 0.088817 2 42 0.8031 0.4547
## hot:trial 0.101708 0.050854 2 42 0.4598 0.6345#no significant effect of treatment in experiment 1. Trial is not significant either, but did not change the effect size of treatment, so I left it in the model for comparison with experiment 2.
#There is no interaction effect of trial*hot in experiment 1.
#Interaction does not improve the model - however, it may be biologically important to test this way bc we may assume learning to happen over trials? simpler model fits here much better but inorder to compare to experiment 2, I decided to move forward with the intraction model.
# Fit your model
model1 <- lmer(light_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.68 0.121 49.4 1.44 1.92
## 2 1.57 0.121 49.4 1.32 1.81
## 3 1.46 0.121 49.4 1.22 1.71
##
## hot = y:
## trial emmean SE df lower.CL upper.CL
## 1 1.56 0.126 49.4 1.31 1.81
## 2 1.55 0.126 49.4 1.30 1.81
## 3 1.53 0.126 49.4 1.27 1.78
##
## 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.11585 0.136 42 0.853 0.6723
## trial1 - trial3 0.21812 0.136 42 1.607 0.2541
## trial2 - trial3 0.10227 0.136 42 0.753 0.7333
##
## hot = y:
## contrast estimate SE df t.ratio p.value
## trial1 - trial2 0.00695 0.142 42 0.049 0.9987
## trial1 - trial3 0.03067 0.142 42 0.216 0.9745
## trial2 - trial3 0.02373 0.142 42 0.167 0.9847
##
## 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.1219 0.175 49.4 0.696 0.4897
##
## trial = 2:
## contrast estimate SE df t.ratio p.value
## n - y 0.0130 0.175 49.4 0.074 0.9413
##
## trial = 3:
## contrast estimate SE df t.ratio p.value
## n - y -0.0656 0.175 49.4 -0.375 0.7096
##
## Degrees-of-freedom method: kenward-roger#No significant differences in light time between treatment groups in exp 1
Estimated Marginal Means plot (Model-predicted group means, adjusting for covariates, trial interactions)
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="Light in Experiment 1", x = "Trial", y = "Estimated Mean", color = "Treatment (hot)")
print(p1)
#visually, a steady decrease in time in the light for control group only, not significant.
Experiment 2
# Experiment 2 treatment model vs null
model2 <- lmer( light_trans ~ hot + trial +
(1 | subject),
data = exp2)
nullmodel2 <- lmer( light_trans ~ 1 + trial +
(1 | subject),
data = exp2)
anova(model2, nullmodel2)## refitting model(s) with ML (instead of REML)## Data: exp2
## Models:
## nullmodel2: light_trans ~ 1 + trial + (1 | subject)
## model2: light_trans ~ hot + trial + (1 | subject)
## npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
## nullmodel2 5 50.632 61.802 -20.316 40.632
## model2 6 45.850 59.254 -16.925 33.850 6.7821 1 0.009208 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(model2)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## hot 0.49283 0.49283 1 21 7.2021 0.013904 *
## trial 0.76569 0.38284 2 44 5.5948 0.006841 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#significant improvement on the null model
#significant treatment effect is coming from experiment 2, trial 2 is the majority of the difference.
#Test interaction: Does hot group show a different trajectory across trials?
modelinteract <- lmer( light_trans ~ hot * trial +
(1 | subject),
data = exp2)
nullmodelinteract <- lmer( light_trans ~ 1 * trial +
(1 | subject),
data = exp2)
anova(modelinteract, nullmodelinteract)## refitting model(s) with ML (instead of REML)## Data: exp2
## Models:
## nullmodelinteract: light_trans ~ 1 * trial + (1 | subject)
## modelinteract: light_trans ~ hot * trial + (1 | subject)
## npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
## nullmodelinteract 3 57.055 63.757 -25.527 51.055
## modelinteract 8 45.321 63.194 -14.661 29.321 21.734 5 0.0005883 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(modelinteract)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: light_trans ~ hot * trial + (1 | subject)
## Data: exp2
##
## REML criterion at convergence: 47.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.0864 -0.6283 -0.0355 0.5096 1.8585
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 0.05286 0.2299
## Residual 0.06497 0.2549
## Number of obs: 69, groups: subject, 23
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.19806 0.10350 44.91782 11.576 4.49e-15 ***
## hoty 0.35845 0.14329 44.91782 2.502 0.0161 *
## trial2 0.21978 0.10868 42.00000 2.022 0.0496 *
## trial3 0.26916 0.10868 42.00000 2.477 0.0174 *
## hoty:trial2 0.07108 0.15046 42.00000 0.472 0.6391
## hoty:trial3 -0.22902 0.15046 42.00000 -1.522 0.1355
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) hoty trial2 trial3 hty:t2
## hoty -0.722
## trial2 -0.525 0.379
## trial3 -0.525 0.379 0.500
## hoty:trial2 0.379 -0.525 -0.722 -0.361
## hoty:trial3 0.379 -0.525 -0.361 -0.722 0.500#Model is significantly better than null
#interaction does not significantly improve model fit over the simpler model without the interaction. However, it may be biologically important to test this way bc we may assume learning to happen over trials? Significant hot effect persists with interaction term added. Control groups experiences significant increases in trials 2 and 3 from trial 1.
model2 <- lmer(light_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.20 0.1030 44.9 0.99 1.41
## 2 1.42 0.1030 44.9 1.21 1.63
## 3 1.47 0.1030 44.9 1.26 1.68
##
## hot = y:
## trial emmean SE df lower.CL upper.CL
## 1 1.56 0.0991 44.9 1.36 1.76
## 2 1.85 0.0991 44.9 1.65 2.05
## 3 1.60 0.0991 44.9 1.40 1.80
##
## 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.2198 0.109 42 -2.022 0.1193
## trial1 - trial3 -0.2692 0.109 42 -2.477 0.0449
## trial2 - trial3 -0.0494 0.109 42 -0.454 0.8927
##
## hot = y:
## contrast estimate SE df t.ratio p.value
## trial1 - trial2 -0.2909 0.104 42 -2.795 0.0208
## trial1 - trial3 -0.0401 0.104 42 -0.386 0.9214
## trial2 - trial3 0.2507 0.104 42 2.409 0.0523
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 3 estimates#Control group: Significantly more time in light in trial 3 vs 1
#Hot group: significantly more 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.358 0.143 44.9 -2.502 0.0161
##
## trial = 2:
## contrast estimate SE df t.ratio p.value
## n - y -0.430 0.143 44.9 -2.998 0.0044
##
## trial = 3:
## contrast estimate SE df t.ratio p.value
## n - y -0.129 0.143 44.9 -0.903 0.3712
##
## Degrees-of-freedom method: kenward-roger#Hot group in the light more often than control across trials, significantly so in trials 1 and 2.
Estimated Marginal Means plot
#(Model-predicted group means, adjusting for covariates, trial interactions)
# Exp1
emm1 <- emmeans(model1, ~ hot * trial)
emm_df1 <- as.data.frame(emm1)
# Convert trial to numeric (ensures spacing)
emm_df1 <- emm_df1 %>%
mutate(trial = as.numeric(as.factor(trial)))
# Exp2
emm2 <- emmeans(model2, ~ hot * trial)
emm_df2 <- as.data.frame(emm2)
emm_df2 <- emm_df2 %>%
mutate(trial = as.numeric(as.factor(trial)))
ymin <- min(emm_df1$lower.CL, emm_df2$lower.CL) - 0.3
ymax <- max(emm_df1$upper.CL, emm_df2$upper.CL) + 0.3
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 = "Light in Experiment 1", x = "Trial", y = "Estimated Mean", color = "Treatment (hot)") +
coord_cartesian(ylim = range(c(emm_df1$lower.CL, emm_df1$upper.CL,
emm_df2$lower.CL, emm_df2$upper.CL))) + # match y-axis
theme_minimal()
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 = "Light in Experiment 2", x = "Trial", y = "Estimated Mean", color = "Treatment (hot)") +
coord_cartesian(ylim = range(c(emm_df1$lower.CL, emm_df1$upper.CL,
emm_df2$lower.CL, emm_df2$upper.CL))) +
theme_minimal()
# Annotate p2
p2 <- p2 +
annotate("text", x = 1, y = max(emm_df2$emmean) + 0.01, label = "*", size = 6) +
annotate("text", x = 2, y = max(emm_df2$emmean) + 0.3, label = "*", size = 6)
# 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.25, label = "a", color = "red", size = 3) +
annotate("text", x = 2, y = hot_df2$emmean[2] + 0.25, label = "a,b", color = "red", size = 3) +
annotate("text", x = 3, y = hot_df2$emmean[3] + 0.25, label = "b", color = "red", size = 3)
# Add grey "a" below control group
control_df2 <- emm_df2 %>% filter(hot == "n")
p2 <- p2 +
annotate("text", x = 1, y = control_df2$emmean[1] - 0.25, label = "a", color = "grey50", size = 3) +
annotate("text", x = 3, y = control_df2$emmean[3] - 0.25, label = "a", color = "grey50", size = 3)
p1 <- p1 + coord_cartesian(ylim = c(ymin, ymax))## Coordinate system already present. Adding new coordinate system, which will
## replace the existing one.p2 <- p2 + coord_cartesian(ylim = c(ymin, ymax))## Coordinate system already present. Adding new coordinate system, which will
## replace the existing one.### --- 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 treatment groups (red for hot, grey for control).
Extras —————————————————————————————–
Side effect?
modelside<- lmer(light_trans~hot + side + (1|subject), data=data)
summary(modelside)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: light_trans ~ hot + side + (1 | subject)
## Data: data
##
## REML criterion at convergence: 115.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.15798 -0.74726 0.05628 0.58451 2.32453
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 0.05164 0.2272
## Residual 0.10112 0.3180
## Number of obs: 138, groups: subject, 23
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.40232 0.08025 26.65678 17.474 4.01e-16 ***
## hoty 0.13944 0.05414 114.00000 2.576 0.0113 *
## sider 0.14125 0.10925 21.00000 1.293 0.2101
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) hoty
## hoty -0.337
## sider -0.651 0.000#no significant effect of side
Time effect - does time of day the trial was conducted effect light?
#shrimp more active at night than during the day - looking at morning and evening vs midday.
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(light_trans ~ hot + hour2 + (1 | subject), data = data)#appears to be significant quadratic effect of time of day tested.
is there an issue with time data and this model?
ggplot(data, aes(x = hour_decimal, y = resid(model_time_quad))) +
geom_point() +
geom_smooth(method = "loess") +
labs(title = "Residuals vs. Time of Day")## `geom_smooth()` using formula = 'y ~ x'
#quadratic curve isn’t fitting the data evenly across time
ggplot(data, aes(x = hour_decimal, y = light_trans)) +
geom_point(alpha = 0.5) +
geom_smooth(method = "lm", formula = y ~ poly(x, 2), se = TRUE) +
labs(title = "Light vs. Time of Day with Quadratic Fit")
hist(data$hour_decimal, breaks = 12, main = "Distribution of Observations by Time of Day", xlab = "Hour")
#This model’s significant quadratic effect may be distorted by sparse data at the end of the day.The ~4 subjects I have at the end of the day (hour 16-18) may be contributing to noise rather than reflecting a true difference.
#Maybe I just state that I tried to keep observations toward the middle of the day and throw out these analyses. I had 2-3 trials which were performed after dusk, otherwise all were during the daytime.
#ALSO - looking at dataset - have 7 shrimp tested after 16:00, 5 were hot=y, 2 were hot=n. could be an effect of treatment. I think I might have to throw our time of day since I can’t separate it from treatment; I dont have even hot/control across time periods.
#Or is there a better way to do this?