data<-Behavioral.Trials.All.Data
head(data)
## exp observer day time tank sex molt length snaplat side hot shrimpID subject
## 1 1 Jesi 2 11:20 3 f n 30 2.906 r n R 1
## 2 1 Jesi 3 9:45 1 m n 31 2.386 l n DDD 3
## 3 1 Jesi 1 12:30 2 m n 26 2.978 r n ZZ 8
## 4 1 Jesi 3 12:15 4 f n 26 43.036 l n EE 10
## 5 1 Jesi 4 11:30 2 f y 27 3.742 r n WW 12
## 6 1 Jesi 4 13:00 3 f n 29 1.598 l n S 13
## trial groomPOST pausePOST groomALL snap tailflip drop dark avgdark middle
## 1 1 6 295.66 17 0 0 0 401.36 36.49 54.16
## 2 1 0 269.57 0 0 0 0 0.00 0.00 600.02
## 3 1 4 0.00 4 0 0 0 355.51 19.75 56.92
## 4 1 2 293.33 6 0 0 0 539.92 77.13 11.12
## 5 1 4 113.10 7 0 0 0 592.12 296.06 8.12
## 6 1 12 275.71 12 0 0 0 339.15 67.83 195.05
## avgmid light avglight pauseCPA groomCPA pauseLIGHT pauseDARK pauseMID
## 1 2.26 144.37 11.11 68.23 11 19.19 23.12 25.93
## 2 600.02 0.00 0.00 595.29 0 0.00 0.00 595.29
## 3 1.67 187.30 11.02 11.14 0 0.00 11.14 0.00
## 4 0.86 49.01 7.00 421.97 4 0.00 421.97 0.00
## 5 8.12 0.00 0.00 345.78 3 0.00 345.78 0.00
## 6 21.67 62.49 10.42 443.16 0 0.00 267.62 175.54
## escapeDARK escapeLIGHT escape
## 1 0 0 0
## 2 0 0 0
## 3 0 0 0
## 4 0 0 0
## 5 0 0 0
## 6 0 0 0
#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",])
# Reshape to long format - all data
long_all <- data %>%
pivot_longer(cols = c(dark, light, middle),
names_to = "zone",
values_to = "seconds")
# Reshape to long format - experiment 1
long_exp1 <- exp1 %>%
pivot_longer(cols = c(dark, light, middle),
names_to = "zone",
values_to = "seconds")
# Reshape to long format - experiment 2
long_exp2 <- exp2 %>%
pivot_longer(cols = c(dark, light, middle),
names_to = "zone",
values_to = "seconds")
ggplot(long_all, aes(x = trial, y = seconds, group = zone)) +
stat_summary(fun = mean, geom = "line", aes(color = zone), size = 1) +
stat_summary(fun.data = mean_se, geom = "errorbar", width = 0.2, color = "black") +
stat_summary(fun = mean, geom = "point", size = 3, color = "black") +
scale_color_manual(
values = c("dark" = "darkblue", "light" = "orange", "middle" = "gray50"),
name = "Zone"
) +
facet_wrap(~ hot, labeller = labeller(hot = c("n" = "Control", "y" = "Hot"))) +
labs(title = "Time Spent in Each Zone Across All Trials",
x = "Trial",
y = "Seconds",
color = "Zone") +
theme_minimal(base_size = 14)
## 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.
ggplot(long_exp1, aes(x = trial, y = seconds, group = zone)) +
stat_summary(fun = mean, geom = "line", aes(color = zone), size = 1) +
stat_summary(fun.data = mean_se, geom = "errorbar", width = 0.2, color = "black") +
stat_summary(fun = mean, geom = "point", size = 3, color = "black") +
scale_color_manual(
values = c("dark" = "darkblue", "light" = "orange", "middle" = "gray50"),
name = "Zone"
) +
coord_cartesian(ylim = c(0, 600)) +
facet_wrap(~ hot, labeller = labeller(hot = c("n" = "Control", "y" = "Hot"))) +
labs(title = "Time Spent in Each Zone in Experiment 1",
x = "Trial",
y = "Seconds",
color = "Zone") +
theme_minimal(base_size = 14)
ggplot(long_exp2, aes(x = trial, y = seconds, group = zone)) +
stat_summary(fun = mean, geom = "line", aes(color = zone), size = 1) +
stat_summary(fun.data = mean_se, geom = "errorbar", width = 0.2, color = "black") +
stat_summary(fun = mean, geom = "point", size = 3, color = "black") +
scale_color_manual(
values = c("dark" = "darkblue", "light" = "orange", "middle" = "gray50"),
name = "Zone"
) +
facet_wrap(~ hot, labeller = labeller(hot = c("n" = "Control", "y" = "Hot"))) +
labs(title = "Time Spent in Each Zone in Experiment 2",
x = "Trial",
y = "Seconds",
color = "Zone") +
theme_minimal(base_size = 14)
#Only graphing at dark and light (a bit easier to look at)
# Reshape to long format
long_all <- data %>%
pivot_longer(cols = c(dark, light),
names_to = "zone",
values_to = "seconds")
# Reshape to long format
long_exp1 <- exp1 %>%
pivot_longer(cols = c(dark, light),
names_to = "zone",
values_to = "seconds")
# Reshape to long format
long_exp2 <- exp2 %>%
pivot_longer(cols = c(dark, light),
names_to = "zone",
values_to = "seconds")
#Dark and light for control and hot overlaid
ggplot(long_all, aes(x = trial, y = seconds, group = interaction(zone, hot))) +
stat_summary(fun = mean, geom = "line",
aes(color = zone, linetype = hot),
size = 1) +
stat_summary(fun.data = mean_se,
geom = "errorbar",
width = 0.2,
color = "black") +
stat_summary(fun = mean,
geom = "point",
size = 3,
color = "black") +
scale_color_manual(
values = c("dark" = "darkblue", "light" = "orange"),
name = "Zone"
) +
scale_linetype_manual(
values = c("n" = "dashed", "y" = "solid"),
name = "Treatment",
labels = c("Control", "Hot")
) +
labs(
title = "Time Spent in Each Zone All Trials",
x = "Trial",
y = "Seconds"
) +
theme_minimal(base_size = 14)
ggplot(long_exp1, aes(x = trial, y = seconds, group = interaction(zone, hot))) +
stat_summary(fun = mean, geom = "line",
aes(color = zone, linetype = hot),
size = 1) +
stat_summary(fun.data = mean_se,
geom = "errorbar",
width = 0.2,
color = "black") +
stat_summary(fun = mean,
geom = "point",
size = 3,
color = "black") +
scale_color_manual(
values = c("dark" = "darkblue", "light" = "orange"),
name = "Zone"
) +
scale_linetype_manual(
values = c("n" = "dashed", "y" = "solid"),
name = "Treatment",
labels = c("Control", "Hot")
) +
coord_cartesian(ylim = c(0, 600)) +
labs(
title = "Time Spent in Each Zone - Exp 1",
x = "Trial",
y = "Seconds"
) +
theme_minimal(base_size = 14)
# Plot
ggplot(long_exp2, aes(x = trial, y = seconds, group = interaction(zone, hot))) +
stat_summary(fun = mean, geom = "line",
aes(color = zone, linetype = hot),
size = 1) +
stat_summary(fun.data = mean_se,
geom = "errorbar",
width = 0.2,
color = "black") +
stat_summary(fun = mean,
geom = "point",
size = 3,
color = "black") +
scale_color_manual(
values = c("dark" = "darkblue", "light" = "orange"),
name = "Zone"
) +
scale_linetype_manual(
values = c("n" = "dashed", "y" = "solid"),
name = "Treatment",
labels = c("Control", "Hot")
) +
labs(
title = "Time Spent in Each Zone Exp 2",
x = "Trial",
y = "Seconds"
) +
theme_minimal(base_size = 14)
#Plots comparing experiments within treatments
#Hot Only
hot_only <- data %>%
filter(hot == "y") %>%
mutate(exp = factor(exp, levels = c(1, 2), labels = c("exp1", "exp2"))) %>%
pivot_longer(cols = c(dark, light, middle),
names_to = "zone",
values_to = "seconds")
# Overlay plot
ggplot(hot_only, aes(x = trial, y = seconds, group = interaction(zone, exp))) +
stat_summary(fun = mean, geom = "line",
aes(color = zone, linetype = exp),
size = 1) +
stat_summary(fun.data = mean_se,
geom = "errorbar",
width = 0.2,
color = "black") +
stat_summary(fun = mean,
geom = "point",
size = 3,
color = "black") +
scale_color_manual(
values = c("dark" = "darkblue", "light" = "orange", "middle" = "gray50"),
name = "Zone"
) +
scale_linetype_manual(
values = c("exp1" = "solid", "exp2" = "dashed"),
name = "Experiment",
labels = c("Experiment 1", "Experiment 2")
) +
labs(
title = "Hot Only: Time in Each Zone by Experiment",
x = "Trial",
y = "Seconds"
) +
theme_minimal(base_size = 14)
#Control only
control_only <- data %>%
filter(hot == "n") %>%
mutate(exp = factor(exp, levels = c(1, 2), labels = c("exp1", "exp2"))) %>%
pivot_longer(cols = c(dark, light, middle),
names_to = "zone",
values_to = "seconds")
# Overlay plot
ggplot(control_only, aes(x = trial, y = seconds, group = interaction(zone, exp))) +
stat_summary(fun = mean, geom = "line",
aes(color = zone, linetype = exp),
size = 1) +
stat_summary(fun.data = mean_se,
geom = "errorbar",
width = 0.2,
color = "black") +
stat_summary(fun = mean,
geom = "point",
size = 3,
color = "black") +
scale_color_manual(
values = c("dark" = "darkblue", "light" = "orange", "middle" = "gray50"),
name = "Zone"
) +
scale_linetype_manual(
values = c("exp1" = "solid", "exp2" = "dashed"),
name = "Experiment",
labels = c("Experiment 1", "Experiment 2")
) +
labs(
title = "Control Only: Time in Each Zone by Experiment",
x = "Trial",
y = "Seconds"
) +
theme_minimal(base_size = 14)
#Same plots, easier to see dark and light
hot_only <- data %>%
filter(hot == "y") %>%
mutate(exp = factor(exp, levels = c(1, 2), labels = c("exp1", "exp2"))) %>%
pivot_longer(cols = c(dark, light),
names_to = "zone",
values_to = "seconds")
# Overlay plot
ggplot(hot_only, aes(x = trial, y = seconds, group = interaction(zone, exp))) +
stat_summary(fun = mean, geom = "line",
aes(color = zone, linetype = exp),
size = 1) +
stat_summary(fun.data = mean_se,
geom = "errorbar",
width = 0.2,
color = "black") +
stat_summary(fun = mean,
geom = "point",
size = 3,
color = "black") +
scale_color_manual(
values = c("dark" = "darkblue", "light" = "orange"),
name = "Zone"
) +
coord_cartesian(ylim = c(0, 600)) +
scale_linetype_manual(
values = c("exp1" = "solid", "exp2" = "dashed"),
name = "Experiment",
labels = c("Experiment 1", "Experiment 2")
) +
labs(
title = "Hot Only: Time in Each Zone by Experiment",
x = "Trial",
y = "Seconds"
) +
theme_minimal(base_size = 14)
control_only <- data %>%
filter(hot == "n") %>%
mutate(exp = factor(exp, levels = c(1, 2), labels = c("exp1", "exp2"))) %>%
pivot_longer(cols = c(dark, light),
names_to = "zone",
values_to = "seconds")
# Overlay plot
ggplot(control_only, aes(x = trial, y = seconds, group = interaction(zone, exp))) +
stat_summary(fun = mean, geom = "line",
aes(color = zone, linetype = exp),
size = 1) +
stat_summary(fun.data = mean_se,
geom = "errorbar",
width = 0.2,
color = "black") +
stat_summary(fun = mean,
geom = "point",
size = 3,
color = "black") +
scale_color_manual(
values = c("dark" = "darkblue", "light" = "orange"),
name = "Zone"
) +
scale_linetype_manual(
values = c("exp1" = "solid", "exp2" = "dashed"),
name = "Experiment",
labels = c("Experiment 1", "Experiment 2")
) +
labs(
title = "Control Only: Time in Each Zone by Experiment",
x = "Trial",
y = "Seconds"
) +
theme_minimal(base_size = 14)
#Learning Effects?
# First, summarize treatments by subject
treat_order <- data %>%
select(subject, exp, hot) %>%
distinct() %>%
pivot_wider(names_from = exp, values_from = hot, names_prefix = "exp") %>%
mutate(order_group = case_when(
exp1 == "y" & exp2 == "n" ~ "HotFirst",
exp1 == "n" & exp2 == "y" ~ "ControlFirst",
TRUE ~ NA_character_
))
# Join it back to the long-format data
data_long <- data %>%
pivot_longer(cols = c(dark, light, middle),
names_to = "zone",
values_to = "seconds") %>%
left_join(treat_order %>% select(subject, order_group), by = "subject") %>%
filter(!is.na(order_group)) # keep only HotFirst and ControlFirst
model <- lmer(seconds ~ zone * order_group + trial + (1 | subject), data = data_long)
## boundary (singular) fit: see help('isSingular')
anova(model)
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## zone 8463835 4231918 2 406 158.9740 < 2.2e-16 ***
## order_group 213 213 1 406 0.0080 0.9287
## trial 380 190 2 406 0.0071 0.9929
## zone:order_group 591256 295628 2 406 11.1054 2.015e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#The zone * treatment-order interaction indicates that shrimp who received hot treatment first behave differently in zones compared to those who received control first — even if total activity or trial progression doesn’t differ
emm <- emmeans(model, ~ order_group | zone)
summary(emm)
## zone = dark:
## order_group emmean SE df lower.CL upper.CL
## ControlFirst 348.9 19.2 155 311.0 387
## HotFirst 454.4 20.1 155 414.7 494
##
## zone = light:
## order_group emmean SE df lower.CL upper.CL
## ControlFirst 124.2 19.2 155 86.2 162
## HotFirst 65.3 20.1 155 25.6 105
##
## zone = middle:
## order_group emmean SE df lower.CL upper.CL
## ControlFirst 126.9 19.2 155 88.9 165
## HotFirst 76.1 20.1 155 36.4 116
##
## Results are averaged over the levels of: trial
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
pairs(emm)
## zone = dark:
## contrast estimate SE df t.ratio p.value
## ControlFirst - HotFirst -105.5 27.8 155 -3.793 0.0002
##
## zone = light:
## contrast estimate SE df t.ratio p.value
## ControlFirst - HotFirst 59.0 27.8 155 2.120 0.0356
##
## zone = middle:
## contrast estimate SE df t.ratio p.value
## ControlFirst - HotFirst 50.8 27.8 155 1.827 0.0696
##
## Results are averaged over the levels of: trial
## Degrees-of-freedom method: kenward-roger
#Dark - Significant difference between treatment-orders. Significantly less time in the dark for controlfirst/hotsecond. Shrimp with prior experience of the test perhaps learned to avoid the dark?
#Light - Significant difference between treatment-orders. Significantly more time in the light for controlfirst/hotsecond. Shrimp with prior experience of the test perhaps learned the light side was safer?
#Middle - not significant, trend to more time in middle for controlfirst/hotsecond. Part of dark avoidance? Different strategy/fear from prior experience?
#Plotting results:
emm_df <- as.data.frame(emm)
# First, store your p-values for annotation
p_vals <- tibble(
zone = c("dark", "light", "middle"),
label = c("p = .0002", "p = .036", "p = .070"),
y_pos = c(525, 200, 200) # vertical position above each bar group
)
ggplot(emm_df, aes(x = zone, y = emmean, fill = order_group)) +
geom_col(position = position_dodge(0.7), width = 0.6) +
geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL),
width = 0.2,
position = position_dodge(0.7),
color = "black") +
scale_fill_manual(
values = c("ControlFirst" = "gray60", "HotFirst" = "red"),
name = "Treatment Order",
labels = c("ControlFirst" = "Control First", "HotFirst" = "Hot First")
) +
labs(
title = "Estimated Time in Each Zone by Treatment Order",
x = "Zone",
y = "Estimated Time (s)"
) +
geom_text(data = p_vals, aes(x = zone, y = y_pos, label = label),
inherit.aes = FALSE, size = 4) +
theme_minimal(base_size = 14)
# Step 0: Recreate order_group from experiment-wise treatment info
order_lookup <- data %>%
select(subject, exp, hot) %>%
distinct() %>%
pivot_wider(names_from = exp, values_from = hot, names_prefix = "exp") %>%
mutate(order_group = case_when(
exp1 == "y" & exp2 == "n" ~ "HotFirst",
exp1 == "n" & exp2 == "y" ~ "ControlFirst",
TRUE ~ NA_character_
))
# Join order_group back into your main data
data <- left_join(data, order_lookup %>% select(subject, order_group), by = "subject")
# Step 1: Reshape to long
long_data <- data %>%
pivot_longer(cols = c(dark, light, middle), names_to = "zone", values_to = "seconds") %>%
group_by(subject, zone, hot, order_group) %>%
summarise(mean_sec = mean(seconds), .groups = "drop")
# Step 2: Pivot to get hot and control side by side
diff_data <- long_data %>%
pivot_wider(names_from = hot, values_from = mean_sec, names_prefix = "hot_") %>%
mutate(diff = hot_y - hot_n) # compute difference per subject and zone
# Step 3: Assign which group the difference belongs to
plot_data <- diff_data %>%
mutate(diff_group = case_when(
order_group == "HotFirst" ~ "Hot First Exp",
order_group == "ControlFirst" ~ "Hot Second Exp",
TRUE ~ NA_character_
))
p_labels <- tibble(
zone = c("dark", "light"),
diff_group = c("Hot First Exp", "Hot First Exp"),
label = c("p = .087", "p = .020"),
y_pos = c(10, -10), # Adjust based on bar height
x_pos = c(.82, 1.82)
)
ggplot(plot_data, aes(x = zone, y = diff, fill = diff_group)) +
stat_summary(fun = mean, geom = "col",
position = position_dodge(0.7), width = 0.6) +
stat_summary(fun.data = mean_se, geom = "errorbar",
position = position_dodge(0.7), width = 0.2, color = "black") +
geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
geom_text(data = p_labels, aes(x = x_pos, y = y_pos, label = label),
position = position_dodge(width = 0.7),
inherit.aes = FALSE, size = 3) +
scale_fill_manual(values = c("Hot First Exp" = "red",
"Hot Second Exp" = "gray50"),
name = "Treatment Comparison") +
labs(title = "Zone Time Differences: Hot − Control",
x = "Zone",
y = "Difference in Time (s)") +
theme_minimal(base_size = 14)
#bars below the 0 line show decreased time in that zone in the hot condition. Bars above the line show increased time when in the hot condition. Separted by treament order.
#Both treatment orders show an increase in time spent in the light when exposed to the noxious heat. When subjects experienced hot treatment first, they spent less time in the dark and those that recieved the hot treatment second spent less time in the middle.
#See below for statistics
# One-sample t-tests: Is diff ≠ 0 within each group?
within_results <- plot_data %>%
group_by(zone, diff_group) %>%
summarise(
n = n(),
mean_diff = mean(diff),
sd_diff = sd(diff),
t_result = list(t.test(diff, mu = 0)),
.groups = "drop"
) %>%
mutate(
t_stat = map_dbl(t_result, ~ .x$statistic),
p_value = map_dbl(t_result, ~ .x$p.value),
df = map_dbl(t_result, ~ .x$parameter),
ci_low = map_dbl(t_result, ~ .x$conf.int[1]),
ci_high = map_dbl(t_result, ~ .x$conf.int[2])
) %>%
select(zone, diff_group, n, mean_diff, sd_diff, t_stat, df, p_value, ci_low, ci_high)
#between group tests - is there a difference between treatment orders?
between_results <- plot_data %>%
group_by(zone) %>%
summarise(
t_result = list(t.test(diff ~ diff_group, var.equal = TRUE)),
.groups = "drop"
) %>%
mutate(
t_stat = map_dbl(t_result, ~ .x$statistic),
p_value = map_dbl(t_result, ~ .x$p.value),
df = map_dbl(t_result, ~ .x$parameter),
ci_low = map_dbl(t_result, ~ .x$conf.int[1]),
ci_high = map_dbl(t_result, ~ .x$conf.int[2])
) %>%
select(zone, t_stat, df, p_value, ci_low, ci_high)
within_results
## # A tibble: 6 × 10
## zone diff_group n mean_diff sd_diff t_stat df p_value ci_low ci_high
## <chr> <chr> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 dark Hot First… 11 -77.6 136. -1.90 10 0.0872 -169. 13.6
## 2 dark Hot Secon… 12 4.22 195. 0.0751 11 0.941 -119. 128.
## 3 light Hot First… 11 62.7 75.6 2.75 10 0.0205 11.9 113.
## 4 light Hot Secon… 12 31.9 97.4 1.13 11 0.281 -30.0 93.8
## 5 middle Hot First… 11 6.85 142. 0.160 10 0.876 -88.5 102.
## 6 middle Hot Secon… 12 -36.4 177. -0.714 11 0.490 -149. 75.8
between_results
## # A tibble: 3 × 6
## zone t_stat df p_value ci_low ci_high
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 dark -1.16 21 0.259 -229. 65.0
## 2 light 0.842 21 0.409 -45.3 107.
## 3 middle 0.644 21 0.527 -96.5 183.
#Significant difference between hot and control for light in Hotfirst (p=0.020). More time spent in light when recieving the hot treatment if you recieved it first. Near-significant decrease in dark time for hot treatment in hotfirst group (0.087).
#No significant differences in the effect of hot between Hotfirst and Hotsecond in dark, light, or middle. Didnt change their outcomes in relation to hot significantly. Biggest difference between treatment-order groups was for dark (clear from the figure).