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",])

Dark, Light, Middle

# 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)

Replotting John’s Graph - look at differences between zone time when recieving hot treatment and zone time when recieving control treatment, within treatment-order groups.

# 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).