analyses.Rmd

---
title: "Detecting trait differences between a reintroduced and a dispersing eastern quoll population"
author: "Christina Gee"
date: "6 June 2024"
output:
  html_document:
    toc: true
    number_sections: true
    toc_depth: 3
    toc_float:
      collapsed: true
    theme: cerulean
    highlight: pygments
editor_options: 
  chunk_output_type: console
knit: (function(inputFile, encoding) { rmarkdown::render(inputFile, encoding = encoding, output_file = file.path(dirname(inputFile), 'tutorial.html')) })
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo=TRUE, eval=TRUE)
```

# * Background

Earth is suffering a biodiversity crisis due to anthropogenic threats including invasive species, habitat destruction, climate change, and human persecution. Important tools for addressing biodiversity loss are creating havens that exclude key threats and restoring the ecosystems within those havens using wildlife reintroductions. However, 25–50% of reintroductions have failed to create self-sustaining populations, due to a number of factors such as poor planning, under-funding, selecting inappropriate founders, failure to address initial causes of decline, and a lack of ongoing post-release monitoring.

Best-practice reintroduction guidelines published by the International Union for the Conservation of Nature state that founders should be selected on the basis of genetic, morphological, physiological, and behavioural representativeness. Each category is important, but several studies have shown the vital influence of individual behavioural syndromes (i.e., personality and plasticity) on reintroduction outcomes. Despite this, there are no known documented cases of selecting founders on the basis of personality and plasticity. This may be due to a lack of rapid field-based behavioural assessment (assay) methods available to practitioners. Most personality studies are conducted in a laboratory environment which is costly and time-intensive and can be stressful and invasive for the animals being assessed.

To address this problem, I designed and tested behavioural assay arenas for use in the field and conducted assays on two eastern quoll (*Dasyurus viverrinus*) populations at Mulligans Flat Woodland Sanctuary (MFWS) and Goorooyarroo Sanctuary (GS), to investigate how differences in personality and plasticity influence population dynamics for this endangered marsupial mesopredator, with a view to understand how knowledge of these differences can inform future reintroductions. I found significant demographic, morphological, and behavioural differences between the eastern quolls of MFWS and GS. The individuals at MFWS were more likely to be behaviourally proactive (lower latencies to emerge, greater activity levels) and plastic (more able to change behaviour over time). Individuals at GS were more likely to be reactive (greater latencies to emerge, and lower activity levels) and more rigid. These quolls were also more likely to be heavier, have greater body condition scores, and aged 1–2 years old (as opposed to MFWS where their ages covered the full life expectancy of up to 3–4 years old).

Understanding how personality and plasticity differences in individual animals affects reintroductions and ecosystem restoration is still in its infancy, but I have demonstrated a method for conducting rapid behavioural assessments in the field. This provides practitioners with a tool to explore behavioural syndromes and how they influence reintroduction outcomes and population dynamics.

# * Setup

First, I installed the [pacman Package Management Tool](https://cran.r-project.org/web/packages/pacman/index.html), which allows us to install and load subsequent packages in a condensed and efficient way. 

```{r, eval=FALSE}
#install.packages("pacman")
```

# Install packages
```{r, results='hide', warning=FALSE, message=FALSE}
# Install and load required packages
pacman::p_load(ggplot2, ggpubr, gtools, janitor, lme4,lmerTest, MuMIn,  
               readr, sf, readxl, rstudioapi, tidyverse, viridis, corrplot,
               broom, sjPlot, sjmisc, sjlabelled, knitr, kableExtra, broom.mixed,
               AICcmodavg, scales)

```


# Data preparation
```{r, eval=FALSE}
# Assign raw data filename to an object
raw_data <- "data.xlsx"

# Read in behavioural assay data
raw <- read_excel(raw_data, sheet="Assays", na="N/A") %>% 
  clean_names() %>%
  select("obs_id"="observation_id", 
         "assay_id", 
         "quoll_code",
         behaviour="behavior", 
         occurrences="total_number_of_occurences")

#Read in the data for all data sheets in excel 
meta <- read_excel(raw_data, sheet="Metadata", na="N/A") %>% 
  clean_names()

quoll <- read_excel(raw_data, sheet="Quoll", na="N/A") %>%
  clean_names()

capture <- read_excel(raw_data, sheet="Captures") %>%
  clean_names()

session <- read_excel(raw_data, sheet="Sessions") %>%
    clean_names()

#Clean up behaviour data from long to wide format 
occ <- pivot_wider(data = raw, 
                   names_from = behaviour, 
                   values_from = occurrences)  %>%
  clean_names() %>%
  mutate(across(all_of(c("section_1", "section_2", "section_3", "section_4", 
                         "moving", "jump_climb", "food", "legs", 
                         "grooming", "resting", "perch")), as.numeric)) %>%
  rename(section_1_o = "section_1", 
         section_2_o = "section_2", 
         section_3_o = "section_3", 
         section_4_o = "section_4", 
         moving_o = "moving", 
         jump_climb_o = "jump_climb", 
         food_o = "food", 
         legs_o = "legs", 
         grooming_o = "grooming", 
         resting_o = "resting", 
         perch_o = "perch")

raw <- read_excel(raw_data, sheet="Assays", , na="N/A") %>% 
  clean_names() %>%
  select("obs_id"="observation_id", 
         "assay_id", 
         quoll_code, 
         behaviour="behavior", 
         duration="total_duration_s")

dur <- pivot_wider(data = raw,
                   names_from = behaviour, 
                   values_from = duration) %>%
  clean_names() %>%
  mutate(across(all_of(c("section_1", "section_2", "section_3", "section_4", 
                         "moving", "jump_climb", "food", "legs", 
                         "grooming", "resting", "perch")), as.numeric)) %>%
  rename(section_1_d = "section_1", 
         section_2_d = "section_2", 
         section_3_d = "section_3", 
         section_4_d = "section_4", 
         moving_d = "moving", 
         jump_climb_d = "jump_climb", 
         food_d = "food", 
         legs_d = "legs", 
         grooming_d = "grooming", 
         resting_d = "resting", 
         perch_d = "perch")

assays <- left_join(occ, dur, by=c("obs_id", "assay_id", "quoll_code"))

#Adding all sections together to create a column with total section visitations
assays$sumsect_o <- rowSums(assays[, c("section_1_o", "section_2_o", "section_3_o", "section_4_o")])

#Adding together all activity occurrences without sections
assays$totalact_o <- rowSums(assays[, c(4:10)])
#view(assays)

# Append metadata for each assay to dataframe
behavdata <- left_join(meta, assays, by=c("assay_id", "quoll_code"))%>%
  filter(!is.na(assay_number))

behavdata <- left_join(behavdata, session, by=c("session"))%>%
  filter(!is.na(assay_number))

quolldata <- left_join(capture, quoll, by=c("quoll_code"))

data <- left_join(behavdata, quolldata, by=c("session", "quoll_code"))
  
data <- data %>%
  select(-min_temp, -session_desc, -capture_date, -new_animal) %>%
  mutate(assay_number = as.factor(assay_number),
         cond_n = as.numeric(ifelse(condition=="Excellent", 4,
                            ifelse(condition=="Good", 3,
                                   ifelse(condition=="Fair", 2, NA)))),
         age_n = as.numeric(ifelse(age_estimate=="<1yr", 0,
                               ifelse(age_estimate=="1-2yrs", 1,
                                      ifelse(age_estimate=="2-3yrs", 2, NA)))))

```

# Calculating time differences 
```{r}
#Change time columns to correct POSIXct format
data$snout <- as.POSIXct(data$snout)
data$process_end <- as.POSIXct(data$process_end)
data$first_light <- as.POSIXct(data$first_light)

#Calculate time differences
data$time_processtoassay <- data$snout - data$process_end
data$time_assaytolight <- data$first_light - data$snout

#Filtering out the two individuals who were assayed in Sept 2023 
data <- data %>%
  filter(!(quoll_code %in% c("53328", "87DE0")))
```

I conducted a correlation analysis for morphological traits, so I could reduce the number of variables. 

# Correlation analysis
```{r}
cor <- data  %>% 
  dplyr::select(age_n, body_weight_kg, cond_n, pes_mean_mm, head_length_mm)

#Running the correlation test
cor <- cor(cor)
view(cor)

#Saving them to a .csv file
write.csv(cor, "data - morphometric correlation matrix.csv")

#Creating and printing correlation matrix .jpgs which colour code the correlations
col_names <- c("Age estimate", "Body weight (kg)", "Condition", "Pes length (mm)", "Head length (mm)")
row_names <- c("Age estimate", "Body weight (kg)", "Condition", "Pes length (mm)", "Head length (mm)")
colnames(cor) <- col_names
rownames(cor) <- row_names

#Correlation of occurrences and duration with sections
jpeg(height=2000, width=2000, file="data corplot.jpeg", type="cairo")
corrplot(cor, method="color",
                  col= colorRampPalette(c("#B8DE29", "#FFFFFF","#2D718E"))(10),
         type="upper", #creates a correlation matrix
         cl.cex=3, number.cex=4, tl.cex=4, diag=FALSE,
         tl.col="black", tl.srt=45, addCoef.col="black")
dev.off()
```

#** SINGLE ASSAYS

I created a data subset to test the demographic, morphometric and behavioural variables for the first assay for all 36 individual eastern quolls. 

# Subsetting data to include first assays only 
```{r}
data1 <- data %>%
      filter(assay_number !=2)

data1$site01 <- ifelse(data1$study_site == "Goorooyarroo", 1, 0)
write.csv(data1, "data1 - first assay data raw.csv")

```

I ensured that the data frame is in the format needed. 
```{r}
data1 <- data1 %>%
  mutate(age_n = as.numeric(as.character(age_n)), 
         condition = factor(condition, 
                            levels = c("Excellent", "Good", "Fair", "Poor")))

data1 <- data1 %>%
  # assign as factor
  mutate(age_estimate = as.factor(age_estimate)) %>%
  # Assign new labels
  mutate(age_estimate = recode(age_estimate, 
                               "1-2yrs" = "1–2 years", 
                               "<1yr" = "<1 year", 
                               "2-3yrs" = "2–3 years"))  %>%
  # reorder new labels
  mutate(age_estimate = factor(age_estimate, 
                               levels = c("2–3 years", "1–2 years", "<1 year")))

```

#Demo/morpho analysis
```{r}
#GLM
mod <- glm(site01 ~ sex + body_weight_kg + age_n + cond_n, 
           data=data1, family=binomial)
summary(mod)
```

#Extract demo/morpho table with kableExtra
```{r}
tidy_results <- tidy(mod) %>%
  mutate_if(is.numeric, ~round(., digits = 3))

table <- tidy_results %>%
  mutate(p.value = ifelse(p.value <= 0.05, cell_spec(p.value, bold = TRUE), p.value)) %>%
  kbl(caption = "Study site", escape = FALSE) %>%
  kable_classic(full_width = F, html_font = "Times New Roman") %>%
  kable_styling(bootstrap_options = "striped", font_size = 15)
table

```

#Demo/morpho plots
```{r}
str(data1)
mutate(as.numeric(data1$site01))

#Age plot
age_plot  <- ggplot() + 
  geom_smooth(data=data1, aes(x=age_n, y=site01), colour = "#1caa85", fill = "#68C6AE" ) + 
  theme_minimal() +
  theme(panel.grid = element_blank(),
        axis.line.x = element_line(color = "black"),
        axis.line.y = element_line(color = "black"),
        text = element_text(size=16)) + 
  scale_y_continuous(breaks=seq(-0.4,1,0.2),labels=seq(-0.4,1,0.2), limits=c(-0.4,1)) +
  ylab("Probability of capture in MF (0) or GS (1)") +
  xlab("Age estimate") +
  annotate("text", x = 1.9, y =0.95, label = "A", size = 8, hjust = 1)

#Condition plot
cond_plot  <- ggplot() + 
 geom_smooth(data=data1, aes(x=cond_n, y=site01), colour = "#2D718E", fill = "#5089A1") + 
  theme_minimal() +
  theme(panel.grid = element_blank(),
        axis.line.x = element_line(color = "black"),
        axis.line.y = element_line(color = "black"),
        text = element_text(size=16)) + 
   scale_y_continuous(breaks=seq(-0.4,1,0.2),labels=seq(-0.4,1,0.2), limits=c(-0.4,1)) +
  ylab(" ") +
  xlab("Body condition score")+
  annotate("text", x = 3.9, y =0.95, label = "C", size = 8, hjust = 1)

weight_plot  <- ggplot() + 
 geom_smooth(data=data1, aes(x=body_weight_kg, y=site01), colour = "#482677", fill = "#674A8E") + 
  theme_minimal() +
  theme(panel.grid = element_blank(),
        axis.line.x = element_line(color = "black"),
        axis.line.y = element_line(color = "black"),
        text = element_text(size=16)) + 
   scale_y_continuous(breaks=seq(-0.4,1,0.2),labels=seq(-0.4,1,0.2), limits=c(-0.4,1)) +
  #scale_x_continuous(breaks=seq(0, 2, 0.5), labels=seq(0,2,0.5), limits=c(0,2))+
  ylab(" ") +
  xlab("Body weight (kg)")+
  annotate("text", x = 1.6, y =0.95, label = "B", size = 8, hjust = -0.5)

#Combine grooming plots
plot <- ggarrange(age_plot, weight_plot, cond_plot, nrow = 1, ncol = 3)
plot <- annotate_figure(plot, bottom = text_grob("Study site", size = 20))
print(plot)
ggsave(filename = "morpho plots.jpeg", plot, width=350, height= 150, units = "mm")

```

#Environmental vars influence GLM analysis
```{r}
mod <- glm(food_d ~ temp + moon_phase + time_assaytolight + time_processtoassay, data=data1)
mod <- glm(food_o ~ temp + moon_phase + time_assaytolight + time_processtoassay, data=data1)

mod <- glm(grooming_d ~ temp + moon_phase + time_assaytolight + time_processtoassay, data=data1)
mod <- glm(grooming_o ~ temp + moon_phase + time_assaytolight + time_processtoassay, data=data1)

mod <- glm(jump_climb_d ~ temp + moon_phase + time_assaytolight + time_processtoassay, data=data1)
mod <- glm(jump_climb_o ~ temp + moon_phase + time_assaytolight + time_processtoassay, data=data1)

mod <- glm(legs_d ~ temp + moon_phase + time_assaytolight + time_processtoassay, data=data1)
mod <- glm(legs_o  ~ temp + moon_phase + time_assaytolight + time_processtoassay, data=data1)

mod <- glm(moving_d ~ temp + moon_phase + time_assaytolight + time_processtoassay, data=data1)
mod <- glm(moving_o ~ temp + moon_phase + time_assaytolight + time_processtoassay, data=data1)

mod <- glm(perch_d ~ temp + moon_phase + time_assaytolight + time_processtoassay, data=data1)
mod <- glm(perch_o ~ temp + moon_phase + time_assaytolight + time_processtoassay, data=data1)

mod <- glm(bait_eaten_g ~ temp + moon_phase + time_assaytolight + time_processtoassay, data=data1)
mod <- glm(sumsect_o ~ temp + moon_phase + time_assaytolight + time_processtoassay, data=data1)
mod <- glm(totalact_o ~ temp + moon_phase + time_assaytolight + time_processtoassay, data=data1)
mod <- glm(emerge_max ~ temp + moon_phase + time_assaytolight + time_processtoassay, data=data1)

mod <- glm(stress_in_trap ~ temp + moon_phase, data=data1)
mod <- glm(stress_in_bag ~ temp + moon_phase, data=data1)
mod <- glm(stress_on_table ~ temp + moon_phase, data=data1)
summary(mod)
```

#Environmental vars tables  with kabelExtra
```{r}

#Tables with kableExtra
response_var <- all.vars(formula(mod))[1]
tidy_results <- tidy(mod) %>%
  mutate_if(is.numeric, ~round(., digits = 3))
table <- tidy_results %>%
  mutate(p.value = ifelse(p.value <= 0.05, cell_spec(p.value, bold = TRUE), p.value)) %>%
  kbl(caption = response_var, escape = FALSE) %>%
  kable_classic(full_width = F, html_font = "Times New Roman") %>%
  kable_styling(bootstrap_options = "striped", font_size = 12)
table

```


Analysis of behaviour indicators
  1. Test the influence of fixed (e.g., site) and random effects (e.g., arena) on behavioural responses (e.g., food interaction duration) using the `lmer()` function.
  2. If multiple fixed effects produce significance, perform model selection using the dredge() function to determine which fixed effects are within delta 2 AICc.

#Capture behaviour GLMMS
Testing behaviour in GLMMs including environmental variables as competing factors in fixed effects model.

```{r}

mod <- lmer(bait_eaten_g ~ study_site + sex + age_n + cond_n + body_weight_kg + (1|arena), 
            data=data1, na.action = "na.fail")

mod <- lmer(emerge_max ~ study_site + sex + age_n + cond_n + body_weight_kg + temp + (1|arena),
            data=data1, na.action = "na.fail")

mod <- lmer(food_d ~ study_site + sex + age_n + cond_n + body_weight_kg + (1|arena), 
            data=data1, na.action = "na.fail")

mod <- lmer(food_o ~ study_site + sex + age_n + cond_n + body_weight_kg + (1|arena), 
            data=data1, na.action = "na.fail")

mod <- lmer(grooming_d ~ study_site + sex + age_n + cond_n + body_weight_kg + temp + (1|arena),
            data=data1, na.action = "na.fail")

mod <- lmer(grooming_o ~ study_site + sex + age_n + cond_n + body_weight_kg + temp + (1|arena), 
            data=data1, na.action = "na.fail")

mod <- lmer(jump_climb_d ~ study_site + sex + age_n + cond_n + body_weight_kg + (1|arena), 
            data=data1, na.action = "na.fail")

mod <- lmer(jump_climb_o ~ study_site + sex + age_n + cond_n + body_weight_kg + (1|arena), 
            data=data1, na.action = "na.fail")

mod <- lmer(legs_d ~ study_site + sex + age_n + cond_n + body_weight_kg + (1|arena), 
            data=data1, na.action = "na.fail")

mod <- lmer(legs_o ~ study_site + sex + age_n + cond_n + body_weight_kg + (1|arena), 
            data=data1, na.action = "na.fail")

mod <- lmer(moving_d ~ study_site + sex + age_n + cond_n + body_weight_kg + time_assaytolight + (1|arena), 
            data=data1, na.action = "na.fail")

mod <- lmer(moving_o ~ study_site + sex + age_n + cond_n + body_weight_kg + time_assaytolight + (1|arena), 
            data=data1, na.action = "na.fail")

mod <- lmer(perch_d ~ study_site + sex + age_n + cond_n + body_weight_kg + temp + moon_phase + 
              time_assaytolight + time_processtoassay + (1|arena), 
            data=data1, na.action = "na.fail")

mod <- lmer(perch_o ~ study_site + sex + age_n + cond_n + body_weight_kg + (1|arena), 
            data=data1, na.action = "na.fail")

mod <- lmer(stress_in_bag ~ study_site + sex + age_n + cond_n + body_weight_kg + (1|arena), 
            data=data1, na.action = "na.fail")

mod <- lmer(stress_in_trap ~ study_site + sex + age_n + cond_n + body_weight_kg + (1|arena), 
            data=data1, na.action = "na.fail")

mod <- lmer(stress_on_table ~ study_site + sex + age_n + cond_n + body_weight_kg + (1|arena), 
            data=data1, na.action = "na.fail")

mod <- lmer(sumsect_o ~ study_site + sex + age_n + cond_n + body_weight_kg + (1|arena), 
            data=data1, na.action = "na.fail")

mod <- lmer(totalact_o ~ study_site + sex + age_n + cond_n + body_weight_kg + (1|arena), 
            data=data1, na.action = "na.fail")

```


```{r}
#Tables with kableExtra
response_var <- all.vars(formula(mod))[1]

tidy_results <- tidy(mod) %>%
  mutate_if(is.numeric, ~round(., digits = 3))
table <- tidy_results %>%
  mutate(p.value = ifelse(p.value <= 0.05, cell_spec(p.value, bold = TRUE), p.value)) %>%
  kbl(caption = response_var, escape = FALSE)
table
```

#Capture behaviour AICc model selection
```{r}

control_ml <- lmerControl(optimizer = "Nelder_Mead")

#Grooming duration 
mod_gd_site <- lmer(grooming_d ~ study_site + (1|arena), 
            data=data1, control = control_ml, na.action = "na.fail")
mod_gd_temp <- lmer(grooming_d ~ temp + (1|arena), 
            data=data1, control = control_ml, na.action = "na.fail")
models <- list(mod_gd_site, mod_gd_temp)
model.names <- c("Site", "Temp")
aictab(cand.set = models, modnames = model.names)

#Grooming occurrences
mod_go_site <- lmer(grooming_o ~ study_site + (1|arena), 
            data=data1, control = control_ml, na.action = "na.fail")
mod_go_temp <- lmer(grooming_o ~ temp + (1|arena), 
            data=data1, control = control_ml, na.action = "na.fail")
models <- list(mod_go_site, mod_go_temp)
model.names <- c("Site", "Temp")
aictab(cand.set = models, modnames = model.names)

#Latency to emerge
mod_emerge_temp <- lmer(emerge_max ~ temp + (1|arena), 
            data=data1, control = control_ml, na.action = "na.fail")
mod_emerge_weight <- lmer(emerge_max ~ body_weight_kg + (1|arena), 
            data=data1, control = control_ml, na.action = "na.fail")
mod_emerge_age <- lmer(emerge_max ~ age_n + (1|arena), 
            data=data1, control = control_ml, na.action = "na.fail")
models <- list(mod_emerge_temp, mod_emerge_weight, mod_emerge_age)
model.names <- c("Temp", "Weight", "Age")
aictab(cand.set = models, modnames = model.names)
```

#Capture behaviour plots
```{r}

#Grooming duration v study site plot
groomd_site <- ggplot(data=data1, aes(x = study_site, y = grooming_d, fill=study_site)) + 
 geom_violin() +
 geom_boxplot(col="grey20", fill="grey20", width=0.02, lwd=0.5) +
  stat_summary(fun="median", geom="point", col="black", fill="white", 
               shape=21, size=3, stroke=1) +
  scale_fill_manual (values = c("#7E67A0","#68C6A5")) +
  theme_minimal() +
  theme(legend.position = "none", 
        panel.grid = element_blank(),
        axis.line.x = element_line(color = "black"),
        axis.line.y = element_line(color = "black"),
        text = element_text(size=18),  
       axis.text.x = element_text(size=18)) +
  scale_y_continuous(breaks=seq(0,110,20),labels=seq(0,110,20), limits=c(0,110)) +
  xlab("") + ylab("Grooming duration (s)") + labs("")  +
  annotate("text", x = "Mulligans Flat", y =110, label = "A", size = 9, hjust = -2)
print(groomd_site)
ggsave(filename = "1st assay - study site v grooming dur.jpeg", groomd_site, width=150, height= 150, units = "mm")

#Grooming occurrences v study site plot
groomo_site <- ggplot(data=data1, aes(x = study_site, y = grooming_o, 
                        fill = study_site)) + 
  geom_violin() +
  geom_boxplot(col="grey20", fill="grey20", width=0.02, lwd=0.5) +
  stat_summary(fun="median", geom="point", col="black", fill="white", 
               shape=21, size=3, stroke=1) +
  scale_fill_manual (values = c("#7E67A0","#68C6A5")) +
  theme_minimal() +
  theme(legend.position = "none", 
        panel.grid = element_blank(),
        axis.line.x = element_line(color = "black"),
        axis.line.y = element_line(color = "black"),
        text = element_text(size=18),  
       axis.text.x = element_text(size=18)) +
  scale_y_continuous(breaks=seq(0,8,2),labels=seq(0,8,2), limits=c(0,8)) +
  xlab("") + ylab(expression(paste("Grooming occurrences ", italic(" (n)")))) + labs("") +
  annotate("text", x = "Mulligans Flat", y =8, label = "B", size = 9, hjust = -2)
print(groomo_site)
ggsave(filename = "1st assay - study site v grooming occ.jpeg", groomo_site, width=150, height= 150, units = "mm")

#Moving duration v study site plot
moved_site <- ggplot(data=data1, aes(x = study_site, y = moving_d, fill=study_site)) + 
 geom_violin() +
 geom_boxplot(col="grey20", fill="grey20", width=0.02, lwd=0.5) +
  stat_summary(fun="median", geom="point", col="black", fill="white", 
               shape=21, size=3, stroke=1) +
  scale_fill_manual (values = c("#7E67A0","#68C6A5")) +
  theme_minimal() +
  theme(legend.position = "none", 
        panel.grid = element_blank(),
        axis.line.x = element_line(color = "black"),
        axis.line.y = element_line(color = "black"),
        text = element_text(size=18),  
       axis.text.x = element_text(size=18)) +
  scale_y_continuous(breaks=seq(0,300,50),labels=seq(0,300,50), limits=c(0,300)) +
  xlab("") + ylab("Quadrupedal movement \nduration (s)") + labs("") +
    annotate("text", x = "Mulligans Flat", y =300, label = "C", size = 9, hjust = -2)
print(moved_site)
#ggsave(filename = "1st assay - study site v moving dur.jpeg", moved_site, width=100, height= 150, units = "mm")

#Combine grooming plots
plot <- ggarrange(groomd_site, groomo_site, moved_site, ncol = 3, nrow=1)
plot <- annotate_figure(plot, 
        #left = text_grob("Adjusted behaviours over two assays", rot=90, size = 16, face = "bold"),
        bottom = text_grob("Study site", size = 20))
print(plot)
ggsave(filename = "1st assay - groom and move violin plots.jpeg", plot, width=350, height= 150, units = "mm")

#Perching vs temp
perchtemp <- ggplot(data=data1, aes(x = temp, y = perch_d)) + 
 geom_point() +
  geom_smooth(method = "lm", se = TRUE, colour = "#1caa85", fill = "#68C6AE") +
  theme_minimal() +
    theme(legend.position = "right", 
        panel.grid = element_blank(),
        axis.line.x = element_line(color = "black"),
        axis.line.y = element_line(color = "black"),
        text = element_text(size=20)) + 
  scale_y_continuous(breaks=seq(0,1000,250),labels=seq(0,1000,250), limits=c(-30, 1000)) +
  xlab(expression(paste("Ambient temperature (", degree, "C)"))) + 
  ylab("Perching duration (s)") + labs("") +
  annotate("text", x =23, y =990, label = "A", size = 10, hjust = 0)
print(perchtemp)
ggsave(filename = "1st assay perch temp.jpeg", perchtemp, width=150, height= 150, units = "mm")

#Latency to emerge vs temp
emergetemp <- ggplot(data=data1, aes(x = temp, y = emerge_max)) + 
 geom_point() +
  geom_smooth(method = "lm", se = TRUE, colour = "#482677", fill = "#674A8E") +
  theme_minimal() +
    theme(legend.position = "right", 
        panel.grid = element_blank(),
        axis.line.x = element_line(color = "black"),
        axis.line.y = element_line(color = "black"),
        text = element_text(size=20)) + 
  scale_x_continuous(breaks=seq(10,25,5),labels=seq(10,25,5), limits=c(10,25)) +
  xlab(expression(paste("Ambient temperature (", degree, "C)"))) + 
  ylab("Latency to emerge (s)") + labs("") +
  annotate("text", x =23, y =59, label = "B", size = 10, hjust = -1)

print(emergetemp)
ggsave(filename = "1st assay emerge temp.jpeg", emergetemp, width=150, height= 150, units = "mm")

#Emerge weight
emergeweight <- ggplot(data=data1, aes(x = body_weight_kg, y = emerge_max)) + 
 geom_point() +
  geom_smooth(method = "lm", se = TRUE, colour = "#2D718E", fill = "#5089A1") +
  theme_minimal() +
    theme(legend.position = "right", 
        panel.grid = element_blank(),
        axis.line.x = element_line(color = "black"),
        axis.line.y = element_line(color = "black"),
        text = element_text(size=20)) + 
  scale_x_continuous(breaks=seq(0.25,1.75,0.5),labels=seq(0.25,1.75,0.5), limits=c(0.25,1.75)) +
  xlab("Body weight (kg)") + 
  ylab("Latency to emerge (s)") + 
  labs("") +
    annotate("text", x =1.6, y =59, label = "C", size = 10, hjust = 0)

print(emergeweight)
ggsave(filename = "1st assay emerge weight.jpeg", emergeweight, width=150, height= 150, units = "mm")


#Combine emerge plots
plot <- ggarrange(perchtemp, emergetemp, emergeweight, ncol = 3, nrow=1)
plot <- annotate_figure(plot, 
        bottom = text_grob("Study site", size = 20))
print(plot)
ggsave(filename = "1st assay emerge plots.jpeg", plot, width=350, height= 150, units = "mm")

```

## ** MULTIPLE ASSAYS

I created a data subset to test the behavioural variables over multiple assays for 17 individual eastern quolls. 

#Subsetting data to animals with repeated assays
```{r}
values_to_exclude <- c("524E3", "53221", "53511", "53512", "53515", "53526", "563A6", "56ABE", "56EB2", "5FC08", "6D409" , "6D53F" , "6E990" , "71A4D" , "71B69" , "86E19" , "86FFC" , "87F61", "BCA09")
data2 <- data %>%
  filter(!quoll_code %in% values_to_exclude)
write.csv(data2, "data2 - individuals with two assays raw data.csv")

```

#Capture and recapture environmental vars GLM analysis
```{r}

mod <- glm(bait_eaten_g ~ temp + moon_phase + time_assaytolight + time_processtoassay, data=data2)

mod <- glm(food_d ~ temp + moon_phase + time_assaytolight + time_processtoassay, data=data2)
mod <- glm(food_o ~ temp + moon_phase + time_assaytolight + time_processtoassay, data=data2)

mod <- glm(grooming_d ~ temp + moon_phase + time_assaytolight + time_processtoassay, data=data2)
mod <- glm(grooming_o ~ temp + moon_phase + time_assaytolight + time_processtoassay, data=data2)

mod <- glm(jump_climb_d ~ temp + moon_phase + time_assaytolight + time_processtoassay, data=data2)
mod <- glm(jump_climb_o ~ temp + moon_phase + time_assaytolight + time_processtoassay, data=data2)

mod <- glm(legs_d ~ temp + moon_phase + time_assaytolight + time_processtoassay, data=data2)
mod <- glm(legs_o ~ temp + moon_phase + time_assaytolight + time_processtoassay, data=data2)

mod <- glm(moving_d ~ temp + moon_phase + time_assaytolight + time_processtoassay, data=data2)
mod <- glm(moving_o ~ temp + moon_phase + time_assaytolight + time_processtoassay, data=data2)

mod <- glm(perch_d ~ temp + moon_phase + time_assaytolight + time_processtoassay, data=data2)
mod <- glm(perch_o ~ temp + moon_phase + time_assaytolight + time_processtoassay, data=data2)

mod <- glm(emerge_max ~ temp + moon_phase + time_assaytolight + time_processtoassay, data=data2)
mod <- glm(sumsect_o ~ temp + moon_phase + time_assaytolight + time_processtoassay, data=data2)
mod <- glm(totalact_o ~ temp + moon_phase + time_assaytolight + time_processtoassay, data=data2)

mod <- glm(stress_in_trap ~ temp + moon_phase, data=data2)
mod <- glm(stress_in_bag ~ temp + moon_phase, data=data2)
mod <- glm(stress_on_table ~ temp + moon_phase, data=data2)

```

#Capture and recapture environmental vars tables with kableExtra
```{r}
response_var <- all.vars(formula(mod))[1]
tidy_results <- tidy(mod) %>%
  mutate_if(is.numeric, ~round(., digits = 3))
table <- tidy_results %>%
  mutate(p.value = ifelse(p.value <= 0.05, cell_spec(p.value, bold = TRUE), p.value)) %>%
  kbl(caption = response_var, escape = FALSE)
table

```

#Capture and recapture behaviour GLMMs with environmental vars
```{r}
data2 <- data2 %>%
  mutate(age_n = as.numeric(as.character(age_n)))

mod <- lmer(bait_eaten_g ~ study_site + sex + age_n + cond_n + body_weight_kg + moon_phase + (1|arena),
            data=data2, na.action = "na.fail")

mod <- lmer(emerge_max ~ study_site + sex + age_n + cond_n + body_weight_kg + temp + (1|arena),
            data=data2, na.action = "na.fail")

mod <- lmer(food_d ~ study_site + sex + age_n + cond_n + body_weight_kg + moon_phase +(1|arena), 
            data=data2, na.action = "na.fail")

mod <- lmer(food_o ~ study_site + sex + age_n + cond_n + body_weight_kg + time_assaytolight + 
              time_processtoassay + (1|arena), data=data2, na.action = "na.fail")

mod <- lmer(grooming_d ~ study_site + sex + age_n + cond_n + body_weight_kg + (1|arena),
            data=data2, na.action = "na.fail")

mod <- lmer(grooming_o ~ study_site + sex + age_n + cond_n + body_weight_kg + (1|arena), 
            data=data2, na.action = "na.fail")

mod <- lmer(jump_climb_d ~ study_site + sex + age_n + cond_n + body_weight_kg + time_assaytolight+ (1|arena), 
            data=data2, na.action = "na.fail")

mod <- lmer(jump_climb_o ~ study_site + sex + age_n + cond_n + body_weight_kg + (1|arena), 
            data=data2, na.action = "na.fail")

mod <- lmer(legs_d ~ study_site + sex + age_n + cond_n + body_weight_kg + moon_phase + (1|arena), 
            data=data2, na.action = "na.fail")

mod <- lmer(legs_o ~ study_site + sex + age_n + cond_n + body_weight_kg + time_assaytolight + 
              time_processtoassay + (1|arena), 
            data=data2, na.action = "na.fail")

mod <- lmer(moving_d ~ study_site + sex + age_n + cond_n + body_weight_kg + moon_phase + 
              time_assaytolight + time_processtoassay + (1|arena), 
            data=data2, na.action = "na.fail")

mod <- lmer(moving_o ~ study_site + sex + age_n + cond_n + body_weight_kg + moon_phase + 
              time_assaytolight + time_processtoassay + (1|arena), 
            data=data2, na.action = "na.fail")

mod <- lmer(perch_d ~ study_site + sex + age_n + cond_n + body_weight_kg + (1|arena), 
            data=data2, na.action = "na.fail")

mod <- lmer(perch_o ~ study_site + sex + age_n + cond_n + body_weight_kg + (1|arena), 
            data=data2, na.action = "na.fail")

mod <- lmer(stress_in_bag ~ study_site + sex + age_n + cond_n + body_weight_kg + (1|arena), 
            data=data2, na.action = "na.fail")

mod <- lmer(stress_in_trap ~ study_site + sex + age_n + cond_n + body_weight_kg + (1|arena), 
            data=data2, na.action = "na.fail")

mod <- lmer(stress_on_table ~ study_site + sex + age_n + cond_n + body_weight_kg + 
            (1|arena), 
            data=data2, na.action = "na.fail")

mod <- lmer(sumsect_o ~ study_site + sex + age_n + cond_n + body_weight_kg + moon_phase + 
              time_processtoassay +  (1|arena), 
            data=data2, na.action = "na.fail")

mod <- lmer(totalact_o ~ study_site + sex + age_n + cond_n + body_weight_kg + moon_phase + 
              time_processtoassay + time_assaytolight + (1|arena), 
            data=data2, na.action = "na.fail")
```

```{r}
#Tables with kableExtra
response_var <- all.vars(formula(mod))[1]
tidy_results <- tidy(mod) %>%
  mutate_if(is.numeric, ~round(., digits = 3))
table <- tidy_results %>%
  mutate(p.value = ifelse(p.value <= 0.05, cell_spec(p.value, bold = TRUE), p.value)) %>%
  kbl(caption = response_var, escape = FALSE)
table

```

#Capture and recapture behaviour AICc model selection
```{r}

control_ml <- lmerControl(optimizer = "Nelder_Mead")

#Food occ   
mod_fo_process <- lmer(food_o ~ time_processtoassay + (1|arena), 
            data=data2, control = control_ml, na.action = "na.fail")
mod_fo_light <- lmer(food_o ~ time_assaytolight + (1|arena), 
            data=data2, control = control_ml, na.action = "na.fail")
models <- list(mod_fo_process, mod_fo_light)
model.names <- c("Process to assay", "Assay to light")
aictab(cand.set = models, modnames = model.names)

#Legs occ   
mod_lo_process <- lmer(legs_o ~ time_processtoassay + (1|arena), 
            data=data2, control = control_ml, na.action = "na.fail")
mod_lo_light <- lmer(legs_o ~ time_assaytolight + (1|arena), 
            data=data2, control = control_ml, na.action = "na.fail")
models <- list(mod_lo_process, mod_lo_light)
model.names <- c("Process to assay", "Assay to light")
aictab(cand.set = models, modnames = model.names)

#Move occ   
mod_mo_process <- lmer(moving_o ~ time_processtoassay + (1|arena), 
            data=data2, control = control_ml, na.action = "na.fail")
mod_mo_light <- lmer(moving_o ~ time_assaytolight + (1|arena), 
            data=data2, control = control_ml, na.action = "na.fail")
models <- list(mod_mo_process, mod_mo_light)
model.names <- c("Process to assay", "Assay to light")
aictab(cand.set = models, modnames = model.names)

#Move dur   
mod_md_process <- lmer(moving_d ~ time_processtoassay + (1|arena), 
            data=data2, control = control_ml, na.action = "na.fail")
mod_md_light <- lmer(moving_d ~ time_assaytolight + (1|arena), 
            data=data2, control = control_ml, na.action = "na.fail")
models <- list(mod_md_process, mod_md_light)
model.names <- c("Process to assay", "Assay to light")
aictab(cand.set = models, modnames = model.names)

#Table stress   
mod_table_process <- lmer(stress_on_table ~ time_processtoassay + (1|arena), 
            data=data2, control = control_ml, na.action = "na.fail")
mod_table_light <- lmer(stress_on_table ~ time_assaytolight + (1|arena), 
            data=data2, control = control_ml, na.action = "na.fail")
models <- list(mod_table_process, mod_table_light)
model.names <- c("Process to assay", "Assay to light")
aictab(cand.set = models, modnames = model.names)

```

#Plots for two assays
```{r}

#Grooming duration v sex 
groomdur <- ggplot(data=data2, aes(x = sex, y = grooming_d, fill = sex)) + 
  geom_boxplot() +
  geom_boxplot(col="grey20", fill="grey20", width=0.02, lwd=0.5) +
  stat_summary(fun="median", geom="point", col="black", fill="black", 
             shape=21, size=1.5, stroke=1) +
  scale_fill_manual (values = c("#7E67A0","#4DBC94")) +
  theme_minimal() +
  theme(legend.position = "none", 
        panel.grid = element_blank(),
        axis.line.x = element_line(color = "black"),
        axis.line.y = element_line(color = "black"),
        text = element_text(size=20)) + 
  scale_y_continuous(breaks=seq(0,80,20),labels=seq(0,80,20), limits=c(0,80)) +
  xlab("Sex") + 
  ylab("Grooming duration (s)") + labs("") +
  scale_x_discrete(labels = c("M" = "Male", "F" = "Female"))+
  annotate("text", x="M", y =80, label = "C", size = 10, hjust = -2)
print(groomdur)
ggsave(filename = "2 assays - sex v grooming.jpeg", groomdur, width=150, height= 150, units = "mm")

#Bait eaten v age
baitage <- ggplot(data=data2, aes(x = age_estimate, y = bait_eaten_g, fill = age_estimate)) + 
  geom_violin() +
  geom_boxplot(col="grey20", fill="grey20", width=0.02, lwd=0.5) +
  stat_summary(fun="median", geom="point", col="black", fill="black", 
             shape=21, size=1.5, stroke=1) +
  scale_fill_manual (values = c("#7E67A0","#68C6A5","#D0E970")) +
  theme_minimal() +
  theme(legend.position = "none", 
        panel.grid = element_blank(),
        axis.line.x = element_line(color = "black"),
        axis.line.y = element_line(color = "black"),
        text = element_text(size=19)) + 
  scale_y_continuous(breaks=seq(0,50,25),labels=seq(0,50,25), limits=c(0,50)) +
  xlab("Age estimate") + 
  ylab("Bait eaten (g)") + labs("") +
  scale_x_discrete(labels = c("<1yr" = "<1 year", "1-2yrs" = "1-2 years", "2-3yrs" ="2-3 years"))+
  annotate("text", x =2, y =50, label = "A", size = 10, hjust = -6)
print(baitage)

#Perhing v site
perchocc <- ggplot(data=data2, aes(x = study_site, y = perch_o, fill = study_site)) + 
  geom_violin() +
  geom_boxplot(col="grey20", fill="grey20", width=0.02, lwd=0.5) +
  stat_summary(fun="median", geom="point", col="black", fill="black", 
             shape=21, size=1.5, stroke=1) +
  scale_fill_manual (values = c("#7E67A0","#68C6A5")) +
  theme_minimal() +
  theme(legend.position = "none", 
        panel.grid = element_blank(),
        axis.line.x = element_line(color = "black"),
        axis.line.y = element_line(color = "black"),
        text = element_text(size=19)) + 
  scale_y_continuous(breaks=seq(0,60,20),labels=seq(0,60,20), limits=c(0,60)) +
  xlab("Study site") + 
  ylab(expression(paste("Perching occurrences ", italic(" (n)")))) +
  annotate("text", x ="Mulligans Flat", y =60, label = "B", size = 10, hjust = -2)
print(perchocc)
ggsave(filename = "2 assays raw - site v perching.jpeg", perchocc, width=150, height= 150, units = "mm")

#Combine  plots
plot <- ggarrange(baitage, perchocc, groomdur, ncol = 3, nrow=1)
print(plot)
ggsave(filename = "2 assays -bait perch emerge.jpeg", plot, width=350, height= 150, units = "mm")

#Total section vs moon 
sum_moon <- ggplot(data=data2, aes(x = moon_phase, y = sumsect_o)) + 
 geom_point() +
  geom_smooth(method = "lm", se = TRUE, colour = "#440154", fill = "#82568D") +
  theme_minimal() +
    theme(legend.position = "right", 
        panel.grid = element_blank(),
        axis.line.x = element_line(color = "black"),
        axis.line.y = element_line(color = "black"),
        text = element_text(size=18)) + 
  scale_y_continuous(breaks=seq(0,100,25),labels=seq(0,100,25), limits=c(0,100)) +
  xlab("Moon illumination (%)") + 
    ylab(expression(paste("Sum of arena section visits", italic(" (n)")))) + labs("") +
  annotate("text", x =100, y =100, label = "A", size = 10, hjust = 1)
print(sum_moon)
ggsave(filename = "2 assay raw - emerge temp.jpeg", emergetemp, width=150, height= 150, units = "mm")

#Food occ v time assay to light 
foodo_light <- ggplot(data=data2, aes(x = time_assaytolight, y = food_o)) + 
 geom_point() +
  geom_smooth(method = "lm", se = TRUE, colour = "#472D7B", fill = "#8473A7") +
  theme_minimal() +
    theme(legend.position = "right", 
        panel.grid = element_blank(),
        axis.line.x = element_line(color = "black"),
        axis.line.y = element_line(color = "black"),
        text = element_text(size=18)) + 
  scale_y_continuous(breaks=seq(0,50,10),labels=seq(0,50,10), limits=c(-5,55)) +
  xlab("Time: assay to first light (mins)") + 
  ylab(expression(paste("Food interaction occurrences ", italic(" (n)")))) + labs("") +
  annotate("text", x =400, y =55, label = "B", size = 10, hjust = 1)
print(foodo_light)
ggsave(filename = "2 assay raw - food light.jpeg", foodo_light, width=150, height= 150, units = "mm")

#Food occ v time process to assay
foodo_process <- ggplot(data=data2, aes(x = time_processtoassay, y = food_o)) + 
 geom_point() +
  geom_smooth(method = "lm", se = TRUE, colour = "#3B528B", fill = "#7C8CB2") +
  theme_minimal() +
    theme(legend.position = "right", 
        panel.grid = element_blank(),
        axis.line.x = element_line(color = "black"),
        axis.line.y = element_line(color = "black"),
        text = element_text(size=18)) + 
  scale_y_continuous(breaks=seq(0,50,10),labels=seq(0,50,10), limits=c(-10,55)) +
  scale_x_continuous(breaks=seq(50,250,50),labels=seq(50,250,50), limits=c(40,250)) +
  xlab("Time: processing to assay (mins)") + 
  ylab(expression(paste("Food interaction occurrences ", italic(" (n)")))) + labs("") +
  annotate("text", x =250, y =55, label = "C", size = 10, hjust = 1)
print(foodo_process)
ggsave(filename = "2 assay raw - food process.jpeg", foodo_process, width=150, height= 150, units = "mm")

#Move dur v time assay to light 
moved_light <- ggplot(data=data2, aes(x = time_assaytolight, y = moving_d)) + 
 geom_point() +
  geom_smooth(method = "lm", se = TRUE, colour = "#2C728E", fill = "#72A1B4") +
  theme_minimal() +
    theme(legend.position = "right", 
        panel.grid = element_blank(),
        axis.line.x = element_line(color = "black"),
        axis.line.y = element_line(color = "black"),
        text = element_text(size=18)) + 
  scale_y_continuous(breaks=seq(0,250,50),labels=seq(0,250,50), limits=c(0,250)) +
  xlab("Time: assay to first light (mins)") + 
  ylab("Quadrupedal movement duration (s)") + labs("") +
  annotate("text", x =400, y =245, label = "D", size = 10, hjust = 1)
print(moved_light)
ggsave(filename = "2 assay raw - legs light.jpeg", legso_light, width=150, height= 150, units = "mm")

#Move occ v time assay to light 
moveo_light <- ggplot(data=data2, aes(x = time_assaytolight, y = moving_o)) + 
 geom_point() +
  geom_smooth(method = "lm", se = TRUE, colour = "#21908C", fill = "#6BB5B2") +
  theme_minimal() +
    theme(legend.position = "right", 
        panel.grid = element_blank(),
        axis.line.x = element_line(color = "black"),
        axis.line.y = element_line(color = "black"),
        text = element_text(size=18)) + 
  scale_y_continuous(breaks=seq(0,200,50),labels=seq(0,200,50), limits=c(0,200)) +
  xlab("Time: assay to first light (mins)") + 
  ylab(expression(paste("Quadrupedal movement occurrences ", italic(" (n)")))) + labs("") +
  annotate("text", x =395, y =195, label = "E", size = 10, hjust = 1)
print(moveo_light)
ggsave(filename = "2 assay raw - legs light.jpeg", legso_light, width=150, height= 150, units = "mm")

#Move occ v time process to assay
moveo_process <- ggplot(data=data2, aes(x = time_processtoassay, y = moving_o)) + 
 geom_point() +
  geom_smooth(method = "lm", se = TRUE, colour = "#27AD81", fill = "#6FC8AB") +
  theme_minimal() +
    theme(legend.position = "right", 
        panel.grid = element_blank(),
        axis.line.x = element_line(color = "black"),
        axis.line.y = element_line(color = "black"),
        text = element_text(size=18)) + 
  scale_x_continuous(breaks=seq(50,300,50),labels=seq(50,300,50), limits=c(50,300)) +
  xlab("Time: processing to assay (mins)") + 
  ylab(expression(paste("Quadrupedal movement occurrences ", italic(" (n)")))) + labs("") +
  annotate("text", x =300, y =195, label = "F", size = 10, hjust = 1)
print(moveo_process)
ggsave(filename = "2 assay raw - legs process.jpeg", legso_process, width=150, height= 150, units = "mm")

#Legs occ v time assay to light 
legso_light <- ggplot(data=data2, aes(x = time_assaytolight, y = legs_o)) + 
 geom_point() +
  geom_smooth(method = "lm", se = TRUE, colour = "#5DC863", fill = "#93DA97") +
  theme_minimal() +
    theme(legend.position = "right", 
        panel.grid = element_blank(),
        axis.line.x = element_line(color = "black"),
        axis.line.y = element_line(color = "black"),
        text = element_text(size=18)) + 
  scale_y_continuous(breaks=seq(0,150,50),labels=seq(0,150,50), limits=c(0,150)) +
  xlab("Time: assay to first light (mins)") + 
  ylab(expression(paste("Hind leg stand occurrences", italic(" (n)")))) + labs("") +
  annotate("text", x =400, y =150, label = "G", size = 10, hjust = 1)
print(legso_light)
ggsave(filename = "2 assay raw - legs light.jpeg", legso_light, width=150, height= 150, units = "mm")

#Legs occ v time process to assay
legso_process <- ggplot(data=data2, aes(x = time_processtoassay, y = legs_o)) + 
 geom_point() +
  geom_smooth(method = "lm", se = TRUE, colour = "#AADC32", fill = "#C6E876") +
  theme_minimal() +
    theme(legend.position = "right", 
        panel.grid = element_blank(),
        axis.line.x = element_line(color = "black"),
        axis.line.y = element_line(color = "black"),
        text = element_text(size=18)) + 
  scale_y_continuous(breaks=seq(0,150,50),labels=seq(0,150,50), limits=c(-20,150)) +
  scale_x_continuous(breaks=seq(50,280,50),labels=seq(50,280,50), limits=c(30,290)) +
  xlab("Time: processing to assay (mins)") + 
  ylab(expression(paste("Hind leg stand occurrences", italic(" (n)")))) + labs("") +
  annotate("text", x =290, y =150, label = "H", size = 10, hjust = 1)
print(legso_process)
ggsave(filename = "2 assay raw - legs process.jpeg", legso_process, width=150, height= 150, units = "mm")

#Latency to emerge vs temp
emergetemp <- ggplot(data=data2, aes(x = temp, y = emerge_max)) + 
 geom_point() +
  geom_smooth(method = "lm", se = TRUE, colour = "#AADC32", fill = "#C6E876") +
  theme_minimal() +
    theme(legend.position = "right", 
        panel.grid = element_blank(),
        axis.line.x = element_line(color = "black"),
        axis.line.y = element_line(color = "black"),
        text = element_text(size=20)) + 
  scale_x_continuous(breaks=seq(10,25,5),labels=seq(10,25,5), limits=c(10,25)) +
  xlab(expression(paste("Ambient temperature (", degree, "C)"))) + 
  ylab("Latency to emerge (s)") + labs("") +
  annotate("text", x =25, y =60, label = "I", size = 10, hjust = 1)
print(emergetemp)
ggsave(filename = "2 assay raw - emerge temp.jpeg", emergetemp, width=150, height= 150, units = "mm")

#Combine  plots
plot <- ggarrange(sum_moon, foodo_light, foodo_process, 
                  moved_light, moveo_light, moveo_process, 
                  legso_light, legso_process, emergetemp, 
                  ncol = 3, nrow=3) 
print(plot)
ggsave(filename = "2 assays - fixed effects.jpeg", plot, width=350, height= 350, units = "mm")

```

#Environmental vars data mutate

I minimised environmental influence on the behaviour models by standardising their values: the means of the environmental variables were all scaled to zero, and the standard deviation to one, as in Wilson et al. (2022).

```{r}
#Mutating the environmental variables 
data2_mutate <- data2 %>%
  mutate_at(vars(5, 44, 65, 66), ~ as.vector(scale(.)))
data2$condition <- factor(data2$condition, levels = c("Excellent", "Good", "Fair", "Poor"))
str(data2_mutate$age_n)
write.csv(data2_mutate, "data2_mutate - both assays env vars adjusted.csv")

```

#Adjusted behaviour variables - extracting values

I fitted GLMs of the adjusted behaviours against assay number and the standardised environmental variables, and at the same time, I extracted intercept (‘personality’ measure) and slope (‘plasticity’ measure) coefficients for each behavioural variable. 

```{r}
mod <- glm(food_o ~ temp + moon_phase + time_assaytolight + time_processtoassay, 
           family=gaussian, data=
             data2_mutate)
             data2_mutate$foodo_p <- predict(mod, type = "response")
             
mod <- glm(food_d ~ temp + moon_phase + time_assaytolight + time_processtoassay, 
           family=gaussian, data=
             data2_mutate)
             data2_mutate$foodd_p <- predict(mod, type = "response")
             
mod <- glm(grooming_o ~ temp + moon_phase + time_assaytolight + time_processtoassay, 
           family=gaussian, data=
             data2_mutate)
             data2_mutate$groomo_p <- predict(mod, type = "response")

mod <- glm(grooming_d ~ temp + moon_phase + time_assaytolight + time_processtoassay, 
           family=gaussian, data=
             data2_mutate)
             data2_mutate$groomd_p <- predict(mod, type = "response")
             
mod <- glm(jump_climb_o ~ temp + moon_phase + time_assaytolight + time_processtoassay, 
           family=gaussian, data=
             data2_mutate)
             data2_mutate$jumpo_p <- predict(mod, type = "response")

mod <- glm(jump_climb_d ~ temp + moon_phase + time_assaytolight + time_processtoassay, 
           family=gaussian, data=
             data2_mutate)
             data2_mutate$jumpd_p <- predict(mod, type = "response")

mod <- glm(legs_d ~ temp + moon_phase + time_assaytolight + time_processtoassay, 
           family=gaussian, data=
             data2_mutate)
             data2_mutate$legsd_p <- predict(mod, type = "response")

mod <- glm(legs_o ~ temp + moon_phase + time_assaytolight + time_processtoassay, 
           family=gaussian, data=
             data2_mutate)
             data2_mutate$legso_p <- predict(mod, type = "response")

mod <- glm(moving_d ~ temp + moon_phase + time_assaytolight + time_processtoassay, 
           family=gaussian, data=
             data2_mutate)
             data2_mutate$moved_p <- predict(mod, type = "response")
             
mod <- glm(moving_o ~ temp + moon_phase + time_assaytolight + time_processtoassay, 
           family=gaussian, data=
             data2_mutate)
             data2_mutate$moveo_p <- predict(mod, type = "response")

mod <- glm(perch_d ~ temp + moon_phase + time_assaytolight + time_processtoassay, 
           family=gaussian, data=
             data2_mutate)
             data2_mutate$perchd_p <- predict(mod, type = "response")
             
mod <- glm(perch_o ~ temp + moon_phase + time_assaytolight + time_processtoassay, 
           family=gaussian, data=
             data2_mutate)
             data2_mutate$percho_p <- predict(mod, type = "response")

mod <- glm(bait_eaten_g ~ temp + moon_phase + time_assaytolight + time_processtoassay, 
           family=gaussian, data=
             data2_mutate)
             data2_mutate$bait_p <- predict(mod, type = "response")
 
mod <- glm(sumsect_o ~ temp + moon_phase + time_assaytolight + time_processtoassay, 
           family=gaussian, data=
             data2_mutate)
             data2_mutate$sumsecto_p <- predict(mod, type = "response")            
             
mod <- glm(totalact_o ~ temp + moon_phase + time_assaytolight + time_processtoassay, 
           family=gaussian, data=
             data2_mutate)
             data2_mutate$totalacto_p <- predict(mod, type = "response")
             
mod <- glm(stress_in_trap ~ temp + moon_phase, family = gaussian, data=data2_mutate)
             data2_mutate$trapstress_p <- predict(mod, type = "response")

mod <- glm(stress_in_bag ~ temp + moon_phase,
           family=gaussian, data=data2_mutate)
             data2_mutate$bagstress_p <- predict(mod, type = "response")

mod <- glm(stress_on_table ~ temp + moon_phase,
           family=gaussian, data=data2_mutate)
             data2_mutate$tablestress_p <- predict(mod, type = "response")

mod <- glm(emerge_max ~ temp + moon_phase + time_assaytolight + time_processtoassay,
           family=gaussian, data=data2_mutate)
             data2_mutate$emerge_p <- predict(mod, type = "response")

```

#Adjusted behaviour variables - create df 

Once the adjusted variables were extracted, I tidied up the data frame. 

```{r}

data2_pred <- data2_mutate %>%
  select(-food_o, -food_d, -grooming_o, -grooming_d, 
                     -jump_climb_o, -jump_climb_d, -legs_o, -legs_d, 
                     -moving_o, -moving_d, -perch_d, -perch_o, -resting_o,
                     -resting_d, -bait_eaten_g, 
                     -sumsect_o, -totalact_o, -emerge_max, -time_to_emerge_s,
                     -arena_entry) %>%
mutate(assay_number = as.factor(assay_number))
write.csv(data2_pred,"data2_pred - both assays adjusted env variables.csv")

```

#Adjusted behaviour variables - test vs assay number
```{r}
mod <- glm(foodd_p ~ assay_number, family=gaussian, data=data2_pred)

mod <- glm(foodo_p ~ assay_number, family=gaussian, data=data2_pred)
mod <- glm(bait_p ~ assay_number, family=gaussian, data=data2_pred)

mod <- glm(groomd_p ~ assay_number, family=gaussian, data=data2_pred)
mod <- glm(groomo_p ~ assay_number, family=gaussian, data=data2_pred)

mod <- glm(jumpd_p ~ assay_number, family=gaussian, data=data2_pred)
mod <- glm(jumpo_p ~ assay_number, family=gaussian, data=data2_pred)

mod <- glm(legsd_p ~ assay_number, family=gaussian, data=data2_pred)
mod <- glm(legso_p ~ assay_number, family=gaussian, data=data2_pred)

mod <- glm(moved_p ~ assay_number, family=gaussian, data=data2_pred)
mod <- glm(moveo_p ~ assay_number, family=gaussian, data=data2_pred)
mod <- glm(emerge_p ~ assay_number, family=gaussian, data=data2_pred)

mod <- glm(perchd_p ~ assay_number, family=gaussian, data=data2_pred)
mod <- glm(percho_p ~ assay_number, family=gaussian, data=data2_pred)

mod <- glm(sumsecto_p ~ assay_number, family=gaussian, data=data2_pred)
mod <- glm(totalacto_p ~ assay_number, family=gaussian, data=data2_pred)

mod <- glm(bagstress_p ~ assay_number, family=gaussian, data=data2_pred)
mod <- glm(trapstress_p ~ assay_number, family=gaussian, data=data2_pred)
mod <- glm(tablestress_p ~ assay_number, family=gaussian, data=data2_pred)

```

```{r}
#Table with kableExtra
response_var <- all.vars(formula(mod))[1]
tidy_results <- tidy(mod) %>%
  mutate_if(is.numeric, ~round(., digits = 3))
table <- tidy_results %>%
  mutate(p.value = ifelse(p.value <= 0.05, cell_spec(p.value, bold = TRUE), p.value)) %>%
  kbl(caption = response_var, escape = FALSE)
table

```

#Adjusted behaviour variables vs assay number plots
```{r}
#Latency to emerge
ed_plot <- ggplot(data=data2_pred, aes(x = assay_number, y = emerge_p, 
                        fill = assay_number)) + 
 geom_violin() +
 geom_boxplot(col="grey20", fill="grey20", width=0.02, lwd=0.5) +
  stat_summary(fun="median", geom="point", col="black", fill="white", 
               shape=21, size=3, stroke=1) +
  scale_fill_manual (values = c("#7E67A0","#CDE868")) +
  theme_minimal() +
  theme(legend.position = "none", 
        panel.grid = element_blank(),
        axis.line.x = element_line(color = "black"),
        axis.line.y = element_line(color = "black"),
        text = element_text(size=19), 
       axis.text.x = element_text(size=19), 
        plot.margin = margin(0.5, 0.5, 0.5, 0.5, "cm")) +  
  scale_y_continuous(breaks=seq(0,50,10),labels=seq(0,50,10), limits=c(0,50)) +
  xlab("") + 
  ylab(expression(paste("Latency to emerge (s)")))+
  labs("")+
  annotate("text", x = "2", y =50, label = "D", size = 10, hjust = -2)
print(ed_plot)
ggsave(filename = "assay number - emerge.jpeg", ed_plot, width=150, height= 150, units = "mm")

#Food interaction occurrences
fo_plot <- ggplot(data=data2_pred, aes(x = assay_number, y = foodo_p, 
                        fill = assay_number)) + 
  geom_violin() +
 geom_boxplot(col="grey20", fill="grey20", width=0.02, lwd=0.5) +
  stat_summary(fun="median", geom="point", col="black", fill="white", 
               shape=21, size=3, stroke=1) +
  scale_fill_manual (values = c("#7E67A0","#CDE868")) +
  theme_minimal() +
  theme(legend.position = "none", 
        panel.grid = element_blank(),
        axis.line.x = element_line(color = "black"),
        axis.line.y = element_line(color = "black"),
        text = element_text(size=19), 
       axis.text.x = element_text(size=19),  
        plot.margin = margin(0.5, 0.5, 0.5, 1, "cm")) + 
  xlab("") + ylab(expression(paste("Food interaction \noccurrences", italic(" (n) ")))) +
  annotate("text", x = "2", y =30, label = "A", size = 10, hjust = -2) + 
  labs("") 
print(fo_plot)
ggsave(filename = "assay number - food occ - site.jpeg", fo_plot, width=200, height= 150, units = "mm")

#Jump climb duration
jd_plot <- ggplot(data=data2_pred, aes(x = assay_number, y = jumpd_p, 
                        fill = assay_number)) + 
  geom_violin() +
 geom_boxplot(col="grey20", fill="grey20", width=0.02, lwd=0.5) +
  stat_summary(fun="median", geom="point", col="black", fill="white", 
               shape=21, size=3, stroke=1) +
  scale_fill_manual (values = c("#7E67A0","#CDE868")) +
  theme_minimal() +
  theme(legend.position = "none", 
        panel.grid = element_blank(),
        axis.line.x = element_line(color = "black"),
        axis.line.y = element_line(color = "black"),
       text = element_text(size=19), 
       axis.text.x = element_text(size=19),  
        plot.margin = margin(0.5, 0.5, 0.5, 1, "cm")) + 
   xlab("") + ylab(expression(paste("Jumping + climbing \nduration (s)" )))  +
  annotate("text", x = "2", y =200, label = "B", size = 10, hjust = -2) +
  labs("") 
print(jd_plot)
ggsave(filename = "assay number - jump dur - site.jpeg", jd_plot, width=200, height= 150, units = "mm")

#Jump climb occurrences
jo_plot <- ggplot(data=data2_pred, aes(x = assay_number, y = jumpo_p, 
                        fill = assay_number)) + 
  geom_violin() +
 geom_boxplot(col="grey20", fill="grey20", width=0.02, lwd=0.5) +
  stat_summary(fun="median", geom="point", col="black", fill="white", 
               shape=21, size=3, stroke=1) +
  scale_fill_manual (values = c("#7E67A0","#CDE868")) +
  theme_minimal() +
  theme(legend.position = "none", 
        panel.grid = element_blank(),
        axis.line.x = element_line(color = "black"),
        axis.line.y = element_line(color = "black"),
        text = element_text(size=19), 
       axis.text.x = element_text(size=19),  
        plot.margin = margin(0.5, 0.5, 0.5, 1, "cm")) + 
  xlab("") + ylab(expression(paste("Jumping + \nclimbing occurrences ", italic(" (n) ")))) +
  annotate("text", x = "2", y =80, label = "C", size = 10, hjust = -2) + 
  labs("") 
print(jo_plot)
ggsave(filename = "assay number - jump occ - site.jpeg", jo_plot, width=200, height= 150, units = "mm")

#Perching occurrences
po_plot <- ggplot(data=data2_pred, aes(x = assay_number, y = percho_p, 
                        fill = assay_number)) + 
  geom_violin() +
 geom_boxplot(col="grey20", fill="grey20", width=0.02, lwd=0.5) +
  stat_summary(fun="median", geom="point", col="black", fill="white", 
               shape=21, size=3, stroke=1) +
  scale_fill_manual (values = c("#7E67A0","#CDE868")) +
  theme_minimal() +
  theme(legend.position = "none", 
        panel.grid = element_blank(),
        axis.line.x = element_line(color = "black"),
        axis.line.y = element_line(color = "black"),
        text = element_text(size=19), 
       axis.text.x = element_text(size=19),  
        plot.margin = margin(0.5, 0.5, 0.5, 1, "cm")) +   
  scale_y_continuous(breaks=seq(0,30,10),labels=seq(0,30,10), limits=c(0,30)) +
  xlab("") + ylab(expression(paste("Perching \noccurrences ", italic(" (n) ")))) + 
    annotate("text", x = "2", y =30, label = "E", size = 10, hjust = -2) + 
  labs("") 
print(po_plot)
ggsave(filename = "assay number - perch occ.jpeg", po_plot, width=150, height= 150, units = "mm")

#Moving duration
md_plot <- ggplot(data=data2_pred, aes(x = assay_number, y = moved_p, 
                        fill = assay_number)) + 
 geom_violin() +
 geom_boxplot(col="grey20", fill="grey20", width=0.02, lwd=0.5) +
  stat_summary(fun="median", geom="point", col="black", fill="white", 
               shape=21, size=3, stroke=1) +
  scale_fill_manual (values = c("#7E67A0","#CDE868")) +
  theme_minimal() +
  theme(legend.position = "none", 
        panel.grid = element_blank(),
        axis.line.x = element_line(color = "black"),
        axis.line.y = element_line(color = "black"),
       text = element_text(size=19), 
       axis.text.x = element_text(size=19),  
        plot.margin = margin(0.5, 0.5, 0.5, 1, "cm")) +   
  scale_y_continuous(breaks=seq(0,160,50),labels=seq(0,160,50), limits=c(0,160)) +
  xlab("") + ylab(expression(paste("Quadrupedal \nmovement duration (s)" ))) +
  annotate("text", x = "2", y =160, label = "F", size = 10, hjust = -2) + 
  labs("") 
print(md_plot)
ggsave(filename = "assay number - move dur - site.jpeg", md_plot, width=200, height= 150, units = "mm")

#Moving occurrences
mo_plot <- ggplot(data=data2_pred, aes(x = assay_number, y = moveo_p, 
                        fill = assay_number)) + 
  geom_violin() +
 geom_boxplot(col="grey20", fill="grey20", width=0.02, lwd=0.5) +
  stat_summary(fun="median", geom="point", col="black", fill="white", 
               shape=21, size=3, stroke=1) +
  scale_fill_manual (values = c("#7E67A0","#CDE868")) +
  theme_minimal() +
  theme(legend.position = "none", 
        panel.grid = element_blank(),
        axis.line.x = element_line(color = "black"),
        axis.line.y = element_line(color = "black"),
        text = element_text(size=19), 
       axis.text.x = element_text(size=19),  
        plot.margin = margin(0.5, 0.5, 0.5, 1, "cm")) +   
  scale_y_continuous(breaks=seq(0,150,50),labels=seq(0,150,50), limits=c(0,150)) +
  xlab("") + ylab(expression(paste("Quadrupedal\n movement occurrences", italic(" (n) ")))) +
  annotate("text", x = "2", y =150, label = "G", size = 10, hjust = -2) + 
  labs("") 
print(mo_plot)
ggsave(filename = "assay number - move occ - site.jpeg", mo_plot, width=200, height= 150, units = "mm")

#Legs occurrences
lo_plot <- ggplot(data=data2_pred, aes(x = assay_number, y = legso_p, 
                        fill = assay_number)) + 
 geom_violin() +
 geom_boxplot(col="grey20", fill="grey20", width=0.02, lwd=0.5) +
  stat_summary(fun="median", geom="point", col="black", fill="white", 
               shape=21, size=3, stroke=1) +
  scale_fill_manual (values = c("#7E67A0","#CDE868")) +
  theme_minimal() +
  theme(legend.position = "none", 
        panel.grid = element_blank(),
        axis.line.x = element_line(color = "black"),
        axis.line.y = element_line(color = "black"),
        text = element_text(size=19), 
       axis.text.x = element_text(size=19),  
        plot.margin = margin(0.5, 0.5, 0.5, 1, "cm")) +   
  scale_y_continuous(breaks=seq(0,110,20),labels=seq(0,110,20), limits=c(0,110)) +
  xlab("") + ylab(expression(paste("Hindleg stand \noccurrences ", italic("(n) ")))) +
  annotate("text", x = "2", y =110, label = "H", size = 10, hjust = -2) +
  labs("") 
print(lo_plot)
ggsave(filename = "assay number - legs occ - site.jpeg", lo_plot, width=200, height= 150, units = "mm")

#Sum total activity 
act_plot <- ggplot(data=data2_pred, aes(x = assay_number, y = totalacto_p, 
                        fill = assay_number)) + 
  geom_violin() +
 geom_boxplot(col="grey20", fill="grey20", width=0.02, lwd=0.5) +
  stat_summary(fun="median", geom="point", col="black", fill="white", 
               shape=21, size=3, stroke=1) +
  scale_fill_manual (values = c("#7E67A0","#CDE868")) +
  theme_minimal() +
  theme(legend.position = "none", 
        panel.grid = element_blank(),
        axis.line.x = element_line(color = "black"),
        axis.line.y = element_line(color = "black"),
        text = element_text(size=19), 
       axis.text.x = element_text(size=19),  
        plot.margin = margin(0.5, 0.5, 0.5, 1, "cm")) +   
  scale_y_continuous(breaks=seq(0,400,100),labels=seq(0,400,100), limits=c(0,400)) +
  xlab("") + ylab(expression(paste("Total activity \noccurrences", italic(" (n) ")))) +
  annotate("text", x = "2", y =400, label = "I", size = 10, hjust = -2) + 
  labs("") 
print(act_plot)
ggsave(filename = "assay number - tot sect - site.jpeg", act_plot, width=200, height= 150, units = "mm")

#Sum section visits
sum_plot <- ggplot(data=data2_pred, aes(x = assay_number, y = sumsecto_p, 
                        fill = assay_number)) + 
  geom_violin() +
 geom_boxplot(col="grey20", fill="grey20", width=0.02, lwd=0.5) +
  stat_summary(fun="median", geom="point", col="black", fill="white", 
               shape=21, size=3, stroke=1) +
  scale_fill_manual (values = c("#7E67A0","#CDE868")) +
  theme_minimal() +
  theme(legend.position = "none", 
        panel.grid = element_blank(),
        axis.line.x = element_line(color = "black"),
        axis.line.y = element_line(color = "black"),
        text = element_text(size=19), 
       axis.text.x = element_text(size=19),  
        plot.margin = margin(0.5, 0.5, 0.5, 1, "cm")) +  
  scale_y_continuous(breaks=seq(0,80,20),labels=seq(0,80,20), limits=c(0,80)) +
  xlab("") + ylab(expression(paste("Total arena \nsection visits", italic(" (n) ")))) +
  annotate("text", x = "2", y =80, label = "J", size = 10, hjust = -2) + 
  labs("") 
print(sum_plot)
ggsave(filename = "assay number - sum sect - site.jpeg", sum_plot, width=200, height= 150, units = "mm")

#Bag stress
bag_plot <- ggplot(data=data2_pred, aes(x = assay_number, y = bagstress_p, 
                        fill = assay_number)) + 
  geom_violin() +
 geom_boxplot(col="grey20", fill="grey20", width=0.02, lwd=0.5) +
  stat_summary(fun="median", geom="point", col="black", fill="white", 
               shape=21, size=3, stroke=1) +
  scale_fill_manual (values = c("#7E67A0","#CDE868")) +
  theme_minimal() +
  theme(legend.position = "none", 
        panel.grid = element_blank(),
        axis.line.x = element_line(color = "black"),
        axis.line.y = element_line(color = "black"),
        text = element_text(size=19), 
       axis.text.x = element_text(size=19), 
        plot.margin = margin(0.5, 0.5, 0.5, 0.5, "cm")) +   
  scale_y_continuous(breaks=seq(1,1.5,0.1),labels=seq(1,1.5,0.1), limits=c(0.9,1.5)) +
  xlab("") + ylab("Adjusted bag\n stress rating") +
  annotate("text", x = "2", y =1.5, label = "K", size = 10, hjust = -2) + 
  labs("") 
print(bag_plot)
ggsave(filename = "assay number - sum sect - site.jpeg", sum_plot, width=200, height= 150, units = "mm")

#Table stress
table_plot <- ggplot(data=data2_pred, aes(x = assay_number, y = tablestress_p, 
                        fill = assay_number)) + 
  geom_violin() +
 geom_boxplot(col="grey20", fill="grey20", width=0.02, lwd=0.5) +
  stat_summary(fun="median", geom="point", col="black", fill="white", 
               shape=21, size=3, stroke=1) +
  scale_fill_manual (values = c("#7E67A0","#CDE868")) +
  theme_minimal() +
  theme(legend.position = "none", 
        panel.grid = element_blank(),
        axis.line.x = element_line(color = "black"),
        axis.line.y = element_line(color = "black"),
        text = element_text(size=19), 
       axis.text.x = element_text(size=19), 
        plot.margin = margin(0.5, 0.5, 0.5, 0.5, "cm")) +  
  scale_y_continuous(breaks=seq(1,1.5,0.1),labels=seq(1,1.5,0.1), limits=c(0.9,1.5)) +
  xlab("") + ylab("Adjusted table \nstress rating") +
  annotate("text", x = "2", y =1.5, label = "L", size = 10, hjust = -2) + 
  labs("") 
print(table_plot)
ggsave(filename = "assay number - sum sect - site.jpeg", sum_plot, width=200, height= 150, units = "mm")

plot <- ggarrange(fo_plot, jd_plot, jo_plot, 
                  ed_plot, po_plot, md_plot, 
                  mo_plot, lo_plot, act_plot, 
                  sum_plot, bag_plot, table_plot,  nrow=4, ncol=3)
plot <- annotate_figure(plot, 
        bottom = text_grob("Assay / capture number", size = 25))
print(plot)
ggsave(filename = "2 assays - activity by assay number.jpeg", plot, width=400, height= 550, units = "mm")

```

#** PERSONALITY AND PLASTICITY 
I created data frame for personality and plasticity analysis with individual eastern quoll demographic and morphological traits.  

```{r}

quoll2 <- data2_pred %>%
  select(-assay_date, -initial_repeat_assay, -trap_location, -snout, 
                     -haunches, -trap_session_date, -trap_session_morning, -first_light,
                     -ssr, -process_start, -process_end, -pouch, -teeth_colour, 
                     -teeth_sharpness, -teeth_shine, -stress_in_trap, -stress_in_bag,
                      -stress_on_table)
unique(quoll2$age_estimate)
quoll2 <- quoll2 %>%
  # assign as factor
  mutate(condition = as.factor(condition), 
         age_estimate = as.factor(age_estimate)) %>%
  # Assign new labels
  mutate(age_estimate = recode(age_estimate, 
                               "1-2yrs" = "1–2 years", 
                               "<1yr" = "<1 year", 
                               "2-3yrs" = "2–3 years"))  %>%
  # reorder new labels
   mutate(condition = factor(condition, 
                               levels = c("Excellent", "Good", 
                                          "Fair", "Poor")), 
          age_estimate = factor(age_estimate, 
                               levels = c("2–3 years", 
                                          "1–2 years", "<1 year")))
str(quoll2)
write.csv(quoll2, "quoll2.csv")

```

#Intercept and slope 

Using the adjusted values I fitted models to extract the intercept (personality) and  slope (plasticity) over multiple assays. 
```{r}

# Create a list of dependent behavioural variables
dependent <- c("foodd_p","foodo_p","groomd_p","groomo_p",
               "jumpd_p","jumpo_p","legsd_p","legso_p", "moved_p",
               "moveo_p", "perchd_p", "percho_p", "bait_p", "sumsecto_p", "totalacto_p", 
               "trapstress_p", "bagstress_p", "tablestress_p", "emerge_p")

# Define the levels in an object
quolls <- levels(factor(quoll2$quoll_code))

# Create a blank dataframe with four columns
coefficients <- data.frame(quoll_code=NULL, intercept=NULL, 
                           slope=NULL, variable=NULL)

# Loop to extract the intercept and slope of each behavioural response
for(v in 1:NROW(dependent)) {
  for(i in 1:NROW(quolls)){
    quoll_code <- quolls[i]
    variable <- dependent[v]
    mod <- lm(eval(as.name(paste0(variable))) ~ 
                assay_number + temp + moon_phase + time_assaytolight + time_processtoassay, 
              data=subset(quoll2, quoll_code==quolls[i]))
    i <- coef(mod)[1]
    s <- coef(mod)[2]
    results <- data.frame(cbind(quoll_code, i, s, variable))
    coefficients <- rbind(coefficients, results)
  }
}

# Convert the df to wide format
coefficients <- reshape(coefficients, idvar="quoll_code", 
                        timevar="variable", direction="wide")

# Define names of the coefficients
names(coefficients) <- gsub("\\.", "_", names(coefficients))

# Remove duplicated coefficients
coefficients <- coefficients[ , !duplicated(colnames(coefficients))]

# Check structure of coefficients
str(coefficients)

#If they are being read as characters, change to numeric 
cols_to_convert <- setdiff(names(coefficients), "quoll_code")
coefficients[cols_to_convert] <- lapply(coefficients[cols_to_convert], as.numeric)

# Join behavioural df with post-release df
ppdata <- left_join(quoll2, coefficients, by="quoll_code", "session") %>%
  distinct(quoll_code, .keep_all = TRUE)

# Remove duplicated column names
ppdata <- ppdata[ , !duplicated(colnames(ppdata))]

# Write to csv file
write_csv(ppdata, "ppdata - intercept and slope data.csv")

```

# Analyse behaviour indicators for personality and plasticity 

  1. I tested the influence of fixed (e.g., site) and random effects (e.g., arena) on behavioural responses (e.g., food interaction duration) using the `lmer()` function.
  2. If multiple fixed effects produced significance, I performed model selection to determine which fixed effects are within delta 2 AICc.

#Intercept GLMMs - testing behaviour intercept (personality) against independent variables
```{r}
mod <- lmer(i_bait_p ~ sex + body_weight_kg + study_site + cond_n + age_n + (1|arena), 
            data=ppdata, na.action = "na.fail")

mod <- lmer(i_emerge_p ~ sex + body_weight_kg + study_site + cond_n + age_n + (1|arena), 
            data=ppdata, na.action = "na.fail")

mod <- lmer(i_foodd_p ~ sex + body_weight_kg + study_site + cond_n + age_n + (1|arena), 
            data=ppdata, na.action = "na.fail")

mod <- lmer(i_foodo_p ~ sex + body_weight_kg + study_site + cond_n + age_n + (1|arena), 
            data=ppdata, na.action = "na.fail")

mod <- lmer(i_groomd_p ~ sex + body_weight_kg + study_site + cond_n + age_n + (1|arena), 
            data=ppdata, na.action = "na.fail")

mod <- lmer(i_groomo_p ~ sex + body_weight_kg + study_site + cond_n + age_n + (1|arena), 
            data=ppdata, na.action = "na.fail")

mod <- lmer(i_jumpd_p ~ sex + body_weight_kg + study_site + cond_n + age_n + (1|arena), 
            data=ppdata, na.action = "na.fail")

mod <- lmer(i_jumpo_p ~ sex + body_weight_kg + study_site + cond_n + age_n + (1|arena), 
            data=ppdata, na.action = "na.fail")

mod <- lmer(i_legsd_p ~ sex + body_weight_kg + study_site + cond_n + age_n + (1|arena), 
            data=ppdata, na.action = "na.fail")

mod <- lmer(i_legso_p ~ sex + body_weight_kg + study_site + cond_n + age_n + (1|arena), 
            data=ppdata, na.action = "na.fail")

mod <- lmer(i_moved_p ~ sex + body_weight_kg + study_site + cond_n + age_n + (1|arena), 
            data=ppdata, na.action = "na.fail")

mod <- lmer(i_moveo_p ~ sex + body_weight_kg + study_site + cond_n + age_n + (1|arena), 
            data=ppdata, na.action = "na.fail")

mod <- lmer(i_perchd_p ~ sex + body_weight_kg + study_site + cond_n + age_n + (1|arena), 
            data=ppdata, na.action = "na.fail")

mod <- lmer(i_percho_p ~ sex + body_weight_kg + study_site + cond_n + age_n + (1|arena), 
            data=ppdata, na.action = "na.fail")

mod <- lmer(i_sumsecto_p ~ sex + body_weight_kg + study_site + cond_n + age_n + (1|arena), 
            data=ppdata, na.action = "na.fail")

mod <- lmer(i_totalacto_p ~ sex + body_weight_kg + study_site + cond_n + age_n + (1|arena), 
            data=ppdata, na.action = "na.fail")

mod <- lmer(i_trapstress_p ~ sex + body_weight_kg + study_site + cond_n + age_n + (1|arena), 
            data=ppdata, na.action = "na.fail")

mod <- lmer(i_bagstress_p ~ sex + body_weight_kg + study_site + cond_n + age_n + (1|arena), 
            data=ppdata, na.action = "na.fail")

mod <- lmer(i_tablestress_p ~ sex + body_weight_kg + study_site + cond_n + age_n + (1|arena), 
            data=ppdata, na.action = "na.fail")

```

#Intercept tables
```{r}
response_var <- all.vars(formula(mod))[1]
tidy_results <- tidy(mod) %>%
  mutate_if(is.numeric, ~round(., digits = 3))
table <- tidy_results %>%
  mutate(p.value = ifelse(p.value <= 0.05, cell_spec(p.value, bold = TRUE), p.value)) %>%
  kbl(caption = response_var, escape = FALSE) %>%
  kable_classic(full_width = F)
table
```

#Intercept AICc model selection
```{r}
#Perching duration
mod_pdi_site <- lmer(i_perchd_p ~ study_site + (1|arena), 
            data=ppdata, control = control_ml, na.action = "na.fail")
mod_pdi_weight <- lmer(i_perchd_p ~ sex*body_weight_kg  + (1|arena), 
            data=ppdata, control = control_ml, na.action = "na.fail")
models <- list(mod_pdi_site, mod_pdi_weight)
model.names <- c("Study site", "Interaction of Sex and body weight kg")
aictab(cand.set = models, modnames = model.names)

```

#Intercept plots
```{r}

#Bag stress
i_bs_site <- ggplot(data=ppdata, aes(x = study_site, y = i_bagstress_p, fill=study_site)) + 
geom_violin() +
 geom_boxplot(col="grey20", fill="grey20", width=0.02, lwd=0.5) +
  stat_summary(fun="median", geom="point", col="black", fill="white", 
               shape=21, size=3, stroke=1) +
  scale_fill_manual (values = c("#6B9BAF","#CDE868")) +
  theme_minimal() +
  theme(legend.position = "none", 
        panel.grid = element_blank(),
        axis.line.x = element_line(color = "black"),
        axis.line.y = element_line(color = "black"),
        text = element_text(size=19), 
       axis.text.x = element_text(size=19),
        plot.margin = margin(0.5, 0.5, 0.5, 0.5, "cm")) +   
  scale_y_continuous(breaks=seq(0.9,1.3,0.1),labels=seq(0.9,1.3,0.1), limits=c(0.9,1.3)) +
  xlab(" ") + 
  ylab("Adjusted bag stress \nrating intercept") + labs("") +
  annotate("text", x = "Mulligans Flat", y = 1.3, label = parse(text = "A"), size = 10, hjust = -2 )
print(i_bs_site)
ggsave(filename = "intercept - study site v bag stress.jpeg", i_bs_site, width=200, height= 150, units = "mm")

#Table stress
i_ts_site <- ggplot(data=ppdata, aes(x = study_site, y = i_tablestress_p, fill=study_site)) + 
 geom_violin() +
 geom_boxplot(col="grey20", fill="grey20", width=0.02, lwd=0.5) +
  stat_summary(fun="median", geom="point", col="black", fill="white", 
               shape=21, size=3, stroke=1) +
  scale_fill_manual (values = c("#6B9BAF","#CDE868")) +
  theme_minimal() +
  theme(legend.position = "none", 
        panel.grid = element_blank(),
        axis.line.x = element_line(color = "black"),
        axis.line.y = element_line(color = "black"),
       text = element_text(size=19), 
       axis.text.x = element_text(size=19), 
        plot.margin = margin(0.5, 0.5, 0.5, 0.5, "cm")) + 
  scale_y_continuous(breaks=seq(0.9,1.3,0.1),labels=seq(0.9,1.3,0.1), limits=c(0.9,1.3)) +
  xlab("") + 
  ylab("Adjusted table stress \nrating intercept") + labs("") +
  annotate("text", x = "Mulligans Flat", y = 1.3, label = parse(text = "B"), size = 10, hjust = -2)
print(i_ts_site)
ggsave(filename = "intercept - study site v table stress.jpeg", i_ts_site, width=200, height= 150, units = "mm")

#Latency to emerge
i_e_site <- ggplot(data=ppdata, aes(x = study_site, y = i_emerge_p, fill=study_site)) + 
 geom_violin() +
 geom_boxplot(col="grey20", fill="grey20", width=0.02, lwd=0.5) +
  stat_summary(fun="median", geom="point", col="black", fill="white", 
               shape=21, size=3, stroke=1) +
  scale_fill_manual (values = c("#6B9BAF","#CDE868")) +
  theme_minimal() +
  theme(legend.position = "none", 
        panel.grid = element_blank(),
        axis.line.x = element_line(color = "black"),
        axis.line.y = element_line(color = "black"),
        text = element_text(size=19), 
       axis.text.x = element_text(size=19), 
        plot.margin = margin(0.5, 0.5, 0.5, 0.5, "cm")) + 
  scale_y_continuous(breaks=seq(0,50,10),labels=seq(0,50,10), limits=c(0,50)) +
  xlab(" ") + 
  ylab("Adjusted latency \nto emerge (s) intercept") + labs("") +
  annotate("text", x = "Mulligans Flat", y = 50, label = parse(text = "C"), size = 10, hjust = -2)
print(i_e_site)
ggsave(filename = "intercept - study site v latency to emerge.jpeg", i_e_site, width=200, height= 150, units = "mm")


#Bait eaten
i_b_site <- ggplot(data=ppdata, aes(x = study_site, y = i_bait_p, fill=study_site)) + 
geom_violin() +
 geom_boxplot(col="grey20", fill="grey20", width=0.02, lwd=0.5) +
  stat_summary(fun="median", geom="point", col="black", fill="white", 
               shape=21, size=3, stroke=1) +
  scale_fill_manual (values = c("#6B9BAF","#CDE868")) +
  theme_minimal() +
  theme(legend.position = "none", 
        panel.grid = element_blank(),
        axis.line.x = element_line(color = "black"),
        axis.line.y = element_line(color = "black"),
       text = element_text(size=19), 
       axis.text.x = element_text(size=19),
        plot.margin = margin(0.5, 0.5, 0.5, 0.5, "cm")) + 
  xlab(" ") + 
  ylab("Adjusted \nbait eaten (g) intercept") + labs("") +
  annotate("text", x = "Mulligans Flat", y = 20, label = parse(text = "F"), size = 10, hjust = -2)
print(i_b_site)
ggsave(filename = "intercept - study site v bait.jpeg", i_b_site, width=200, height= 150, units = "mm")

#Groom duration
i_gd_site <- ggplot(data=ppdata, aes(x = study_site, y = i_groomd_p, fill=study_site)) + geom_violin() +
 geom_boxplot(col="grey20", fill="grey20", width=0.02, lwd=0.5) +
  stat_summary(fun="median", geom="point", col="black", fill="white", 
               shape=21, size=3, stroke=1) +
  scale_fill_manual (values = c("#6B9BAF","#CDE868")) +
  theme_minimal() +
  theme(legend.position = "none", 
        panel.grid = element_blank(),
        axis.line.x = element_line(color = "black"),
        axis.line.y = element_line(color = "black"),
       text = element_text(size=19), 
       axis.text.x = element_text(size=19), 
        plot.margin = margin(0.5, 0.5, 0.5, 0.5, "cm")) + 
  xlab("") + 
  ylab("Adjusted grooming\n duration (s) intercept") + labs("") +
  annotate("text", x = "Mulligans Flat", y = 25, label = parse(text = "E"), size = 10, hjust = -2)
print(i_gd_site)
ggsave(filename = "intercept - study site v grooming dur.jpeg", i_gd_site, width=150, height= 150, units = "mm")

#Hind legs duration
i_ld_site <- ggplot(data=ppdata, aes(x = study_site, y = i_legsd_p, fill=study_site)) + 
 geom_violin() +
 geom_boxplot(col="grey20", fill="grey20", width=0.02, lwd=0.5) +
  stat_summary(fun="median", geom="point", col="black", fill="white", 
               shape=21, size=3, stroke=1) +
  scale_fill_manual (values = c("#6B9BAF","#CDE868")) +
  theme_minimal() +
  theme(legend.position = "none", 
        panel.grid = element_blank(),
        axis.line.x = element_line(color = "black"),
        axis.line.y = element_line(color = "black"),
       text = element_text(size=19), 
       axis.text.x = element_text(size=19), 
        plot.margin = margin(0.5, 0.5, 0.5, 0.5, "cm")) + 
  scale_y_continuous(breaks=seq(100,300,100),labels=seq(100,300,100), limits=c(100,300)) +
  xlab(" ") + 
  ylab("Adjusted hindleg stands \nduration (s) intercept") + labs("") +
  annotate("text", x = "Mulligans Flat", y = 300, label = parse(text = "D"), size = 10, hjust = -2)
print(i_ld_site)
ggsave(filename = "intercept - study site v legs dur.jpeg", i_ld_site, width=200, height= 150, units = "mm")

#Jump/climb duration
i_jd_site <- ggplot(data=ppdata, aes(x = study_site, y = i_jumpd_p, fill=study_site)) + 
geom_violin() +
 geom_boxplot(col="grey20", fill="grey20", width=0.02, lwd=0.5) +
  stat_summary(fun="median", geom="point", col="black", fill="white", 
               shape=21, size=3, stroke=1) +
  scale_fill_manual (values = c("#6B9BAF","#CDE868")) +
  theme_minimal() +
  theme(legend.position = "none", 
        panel.grid = element_blank(),
        axis.line.x = element_line(color = "black"),
        axis.line.y = element_line(color = "black"),
       text = element_text(size=19), 
       axis.text.x = element_text(size=19), 
        plot.margin = margin(0.5, 0.5, 0.5, 0.5, "cm")) + 
  scale_y_continuous(breaks=seq(0,200,50),labels=seq(0,200,50), limits=c(0,220)) +
  xlab(" ") + 
  ylab("Adjusted jumping + climbing \nduration (s) intercept") + labs("") +
  annotate("text", x = "Mulligans Flat", y = 220, label = parse(text = "G"), size = 10, hjust = -2)
print(i_jd_site)
ggsave(filename = "intercept - study site v jumpclimb dur.jpeg", i_jd_site, width=200, height= 150, units = "mm")

#Jump/climb occurrences
i_jo_site <- ggplot(data=ppdata, aes(x = study_site, y = i_jumpo_p, fill=study_site)) + 
 geom_violin() +
 geom_boxplot(col="grey20", fill="grey20", width=0.02, lwd=0.5) +
  stat_summary(fun="median", geom="point", col="black", fill="white", 
               shape=21, size=3, stroke=1) +
  scale_fill_manual (values = c("#6B9BAF","#CDE868")) +
  theme_minimal() +
  theme(legend.position = "none", 
        panel.grid = element_blank(),
        axis.line.x = element_line(color = "black"),
        axis.line.y = element_line(color = "black"),
     text = element_text(size=19), 
       axis.text.x = element_text(size=19), 
        plot.margin = margin(0.5, 0.5, 0.5, 1, "cm")) + 
  scale_y_continuous(breaks=seq(0,80,20),labels=seq(0,80,20), limits=c(0,90)) +
  xlab(" ") + 
  ylab(expression(paste("Adjusted jumping + climbing\n occurrences intercept",italic("(n)")))) +
  labs("") +
  annotate("text", x = "Mulligans Flat", y = 90, label = parse(text = "H"), size = 10, hjust = -2)
print(i_jo_site)
ggsave(filename = "intercept - study site v jumpclimb occ.jpeg", i_jo_site, width=200, height= 150, units = "mm")


#Perch duration
i_pd_site <- ggplot(data=ppdata, aes(x = study_site, y = i_perchd_p, fill=study_site)) + geom_violin() +
 geom_boxplot(col="grey20", fill="grey20", width=0.02, lwd=0.5) +
  stat_summary(fun="median", geom="point", col="black", fill="white", 
               shape=21, size=3, stroke=1) +
  scale_fill_manual (values = c("#6B9BAF","#CDE868")) +
  theme_minimal() +
  theme(legend.position = "none", 
        panel.grid = element_blank(),
        axis.line.x = element_line(color = "black"),
        axis.line.y = element_line(color = "black"),
       text = element_text(size=19), 
       axis.text.x = element_text(size=19),
        plot.margin = margin(0.5, 0.5, 0.5, 0.5, "cm")) + 
  xlab(" ") + 
  ylab("Adjusted perching\n duration (s) intercept") + labs("") +
  annotate("text", x = "Mulligans Flat", y = 300, label = parse(text = "J"), size = 10, hjust = -2)
print(i_pd_site)
ggsave(filename = "intercept - study site v perch dur.jpeg", i_pd_site, width=200, height= 150, units = "mm")

#Perch occurrences
i_po_site <- ggplot(data=ppdata, aes(x = study_site, y = i_percho_p, fill=study_site)) + geom_violin() +
 geom_boxplot(col="grey20", fill="grey20", width=0.02, lwd=0.5) +
  stat_summary(fun="median", geom="point", col="black", fill="white", 
               shape=21, size=3, stroke=1) +
  scale_fill_manual (values = c("#6B9BAF","#CDE868")) +
  theme_minimal() +
  theme(legend.position = "none", 
        panel.grid = element_blank(),
        axis.line.x = element_line(color = "black"),
        axis.line.y = element_line(color = "black"),
        text = element_text(size=19), 
       axis.text.x = element_text(size=19), 
        plot.margin = margin(0.5, 0.5, 0.5, 1, "cm")) + 
  xlab(" ") + 
  ylab(expression(paste("Adjusted perching \noccurrences intercept ", italic("(n)")))) +
  annotate("text", x = "Mulligans Flat", y = 30, label = parse(text = "K"), size = 10, hjust = -2)
print(i_po_site)
ggsave(filename = "intercept - study site v perch occ.jpeg", i_po_site, width=200, height= 150, units = "mm")

#Move duration
i_md_site <- ggplot(data=ppdata, aes(x = study_site, y = i_moved_p, fill=study_site)) + 
 geom_violin() +
 geom_boxplot(col="grey20", fill="grey20", width=0.02, lwd=0.5) +
  stat_summary(fun="median", geom="point", col="black", fill="white", 
               shape=21, size=3, stroke=1) +
  scale_fill_manual (values = c("#6B9BAF","#CDE868")) +
  theme_minimal() +
  theme(legend.position = "none", 
        panel.grid = element_blank(),
        axis.line.x = element_line(color = "black"),
        axis.line.y = element_line(color = "black"),
        text = element_text(size=19), 
       axis.text.x = element_text(size=19), 
        plot.margin = margin(0.5, 0.5, 0.5, 0.5, "cm")) + 
  xlab(" ") + 
  ylab("Adjusted quadrupedal movement\n duration intercept  (s)") + 
  labs("") +
  annotate("text", x = "Mulligans Flat", y = 160, label = parse(text = "I"), size = 10, hjust = -2)
  print(i_md_site)
ggsave(filename = "intercept - study site v move dur.jpeg", i_md_site, width=200, height= 150, units = "mm")

#Sum of section visits
i_sumsect_site <- ggplot(data=ppdata, aes(x = study_site, y = i_sumsecto_p, fill=study_site)) + 
 geom_violin() +
 geom_boxplot(col="grey20", fill="grey20", width=0.02, lwd=0.5) +
  stat_summary(fun="median", geom="point", col="black", fill="white", 
               shape=21, size=3, stroke=1) +
 scale_fill_manual (values = c("#6B9BAF","#CDE868")) +
  theme_minimal() +
  theme(legend.position = "none", 
        panel.grid = element_blank(),
        axis.line.x = element_line(color = "black"),
        axis.line.y = element_line(color = "black"),
        text = element_text(size=19), 
       axis.text.x = element_text(size=19),
        plot.margin = margin(0.5, 0.5, 0.5, 1, "cm")) + 
  scale_y_continuous(breaks=seq(30,70,10),labels=seq(30,70,10), limits=c(30,70)) +
  xlab("") + 
 ylab(expression(paste("Adjusted total arena \nsection visits intercept ", italic("(n)")))) +  annotate("text", x = "Mulligans Flat", y = 70, label = parse(text = "L"), size = 9, hjust = -2)
print(i_sumsect_site)
ggsave(filename = "intercept - study site v sum sections.jpeg", i_sumsect_site, width=200, height= 150, units = "mm")

#Total activity occurrences
i_act_site <- ggplot(data=ppdata, aes(x = study_site, y = i_totalacto_p, fill=study_site)) + 
 geom_violin() +
 geom_boxplot(col="grey20", fill="grey20", width=0.02, lwd=0.5) +
  stat_summary(fun="median", geom="point", col="black", fill="white", 
               shape=21, size=3, stroke=1) +
  scale_fill_manual (values = c("#6B9BAF","#CDE868")) +
  theme_minimal() +
  theme(legend.position = "none", 
        panel.grid = element_blank(),
        axis.line.x = element_line(color = "black"),
        axis.line.y = element_line(color = "black"),
        text = element_text(size=19), 
       axis.text.x = element_text(size=19),
        plot.margin = margin(0.5, 0.5, 0.5, 1, "cm")) + 
  scale_y_continuous(breaks=seq(150,350,100),labels=seq(150,350,100), limits=c(150,350)) +
  xlab("") + 
  ylab(expression(paste("Adjusted total activity \noccurrences intercept", 
                        italic("(n)")))) +
  annotate("text", x = "Mulligans Flat", y = 350, label = parse(text = "M"), size = 9, hjust = -2)
print(i_act_site)
ggsave(filename = "intercept - study site v total activity.jpeg", i_act_site, width=150, height= 150, units = "mm")

#Combine plots
plot <- ggarrange(i_bs_site, i_ts_site, i_e_site, 
                  i_ld_site, i_gd_site, i_b_site, 
                  i_jd_site, i_jo_site,i_md_site, 
                  i_pd_site,i_po_site, i_sumsect_site,
                  i_act_site, nrow=5,ncol=3)
plot <- annotate_figure(plot,
        bottom = text_grob("Study site", size = 25))
print(plot)
ggsave(filename = "intercept plots.jpeg", plot, width=450, height= 650, units = "mm")

```

#Slope GLMMs - testing behaviour slope (plasticity) against independent variables
```{r}
mod <- lmer(s_bait_p ~ sex + body_weight_kg + study_site + cond_n + age_n + (1|arena), 
            data=ppdata, na.action = "na.fail")

mod <- lmer(s_emerge_p ~ sex + body_weight_kg + study_site + cond_n + age_n + (1|arena), 
            data=ppdata, na.action = "na.fail")

mod <- lmer(s_foodd_p ~ sex + body_weight_kg + study_site + cond_n + age_n + (1|arena), 
            data=ppdata, na.action = "na.fail")

mod <- lmer(s_foodo_p ~ sex + body_weight_kg + study_site + cond_n + age_n + (1|arena), 
            data=ppdata, na.action = "na.fail")

mod <- lmer(s_groomd_p ~ sex + body_weight_kg + study_site + cond_n + age_n + (1|arena), 
            data=ppdata, na.action = "na.fail")

mod <- lmer(s_groomo_p ~ sex + body_weight_kg + study_site + cond_n + age_n + (1|arena), 
            data=ppdata, na.action = "na.fail")

mod <- lmer(s_jumpd_p ~ sex + body_weight_kg + study_site + cond_n + age_n + (1|arena), 
            data=ppdata, na.action = "na.fail")

mod <- lmer(s_jumpo_p ~ sex + body_weight_kg + study_site + cond_n + age_n + (1|arena), 
            data=ppdata, na.action = "na.fail")

mod <- lmer(s_legsd_p ~ sex + body_weight_kg + study_site + cond_n + age_n + (1|arena), 
            data=ppdata, na.action = "na.fail")

mod <- lmer(s_legso_p ~ sex + body_weight_kg + study_site + cond_n + age_n + (1|arena), 
            data=ppdata, na.action = "na.fail")

mod <- lmer(s_moved_p ~ sex + body_weight_kg + study_site + cond_n + age_n + (1|arena), 
            data=ppdata, na.action = "na.fail")

mod <- lmer(s_moveo_p ~ sex + body_weight_kg + study_site + cond_n + age_n + (1|arena), 
            data=ppdata, na.action = "na.fail")

mod <- lmer(s_perchd_p ~ sex + body_weight_kg + study_site + cond_n + age_n + (1|arena), 
            data=ppdata, na.action = "na.fail")

mod <- lmer(s_percho_p ~ sex + body_weight_kg + study_site + cond_n + age_n + (1|arena), 
            data=ppdata, na.action = "na.fail")

mod <- lmer(s_sumsecto_p ~ sex + body_weight_kg + study_site + cond_n + age_n + (1|arena), 
            data=ppdata, na.action = "na.fail")

mod <- lmer(s_totalacto_p ~ sex + body_weight_kg + study_site + cond_n + age_n + (1|arena), 
            data=ppdata, na.action = "na.fail")

mod <- lmer(s_bagstress_p ~ sex + body_weight_kg + study_site + cond_n + age_n + (1|arena), 
            data=ppdata, na.action = "na.fail")

mod <- lmer(s_trapstress_p ~ sex + body_weight_kg + study_site + cond_n + age_n + (1|arena), 
            data=ppdata, na.action = "na.fail")

mod <- lmer(s_tablestress_p ~ sex + body_weight_kg + study_site + cond_n + age_n + (1|arena), 
            data=ppdata, na.action = "na.fail")
```

#Tables with kableExtra
```{r}
response_var <- all.vars(formula(mod))[1]
tidy_results <- tidy(mod) %>%
  mutate_if(is.numeric, ~round(., digits = 3))
table <- tidy_results %>%
  mutate(p.value = ifelse(p.value <= 0.05, cell_spec(p.value, bold = TRUE), p.value)) %>%
  kbl(caption = response_var, escape = FALSE) %>%
  kable_classic(full_width = F)
table
```

#Slope AICc model selection
```{r}
#Food occ
mod_fos_cond <- lmer(s_foodo_p ~ cond_n + (1|arena), 
            data=ppdata, control = control_ml, na.action = "na.fail")

mod_fos_sex <- lmer(s_foodo_p ~ sex + (1|arena), 
            data=ppdata, control = control_ml, na.action = "na.fail")

mod_fos_weight <- lmer(s_foodo_p ~ body_weight_kg  + (1|arena), 
            data=ppdata, control = control_ml, na.action = "na.fail")

models <- list(mod_fos_cond, mod_fos_sex, mod_fos_weight)
model.names <- c("Condition", "Sex", "Body weight")
aictab(cand.set = models, modnames = model.names)

#Legs occ
mod_los_cond <- lmer(s_legso_p ~ cond_n + (1|arena), 
            data=ppdata, control = control_ml, na.action = "na.fail")

mod_los_sex <- lmer(s_legso_p ~ sex + (1|arena), 
            data=ppdata, control = control_ml, na.action = "na.fail")

mod_los_weight <- lmer(s_legso_p ~ body_weight_kg  + (1|arena), 
            data=ppdata, control = control_ml, na.action = "na.fail")

mod_los_site <- lmer(s_legso_p ~ study_site  + (1|arena), 
            data=ppdata, control = control_ml, na.action = "na.fail")

models <- list(mod_los_cond, mod_los_sex, mod_los_weight, mod_los_site)
model.names <- c("Condition", "Sex", "Body weight", "Study site")
aictab(cand.set = models, modnames = model.names)

#Jump occ
mod_jos_cond <- lmer(s_jumpo_p ~ cond_n + (1|arena), 
            data=ppdata, control = control_ml, na.action = "na.fail")

mod_jos_weight <- lmer(s_jumpo_p ~ body_weight_kg  + (1|arena), 
            data=ppdata, control = control_ml, na.action = "na.fail")

mod_jos_site <- lmer(s_jumpo_p ~ study_site  + (1|arena), 
            data=ppdata, control = control_ml, na.action = "na.fail")

models <- list(mod_jos_cond, mod_jos_weight, mod_jos_site)
model.names <- c("Condition", "Body weight", "Study site")
aictab(cand.set = models, modnames = model.names)

#Move dur
mod_mds_cond <- lmer(s_moved_p ~ cond_n + (1|arena), 
            data=ppdata, control = control_ml, na.action = "na.fail")

mod_mds_weight <- lmer(s_moved_p ~ body_weight_kg  + (1|arena), 
            data=ppdata, control = control_ml, na.action = "na.fail")

mod_mds_site <- lmer(s_moved_p ~ study_site  + (1|arena), 
            data=ppdata, control = control_ml, na.action = "na.fail")

models <- list(mod_mds_cond, mod_mds_weight, mod_mds_site)
model.names <- c("Condition", "Body weight", "Study site")
aictab(cand.set = models, modnames = model.names)

#Move occ
mod_mos_cond <- lmer(s_moveo_p ~ cond_n + (1|arena), 
            data=ppdata, control = control_ml, na.action = "na.fail")

mod_mos_weight <- lmer(s_moveo_p ~ body_weight_kg  + (1|arena), 
            data=ppdata, control = control_ml, na.action = "na.fail")

mod_mos_site <- lmer(s_moveo_p ~ study_site  + (1|arena), 
            data=ppdata, control = control_ml, na.action = "na.fail")

mod_mos_sex <- lmer(s_moveo_p ~ sex  + (1|arena), 
            data=ppdata, control = control_ml, na.action = "na.fail")

models <- list(mod_mos_cond, mod_mos_weight, mod_mos_site, mod_mos_sex)
model.names <- c("Condition", "Body weight", "Study site", "Sex")
aictab(cand.set = models, modnames = model.names)

#Table stress
mod_tabs_cond <- lmer(s_tablestress_p ~ cond_n + (1|arena), 
            data=ppdata, control = control_ml, na.action = "na.fail")

mod_tabs_weight <- lmer(s_tablestress_p ~ body_weight_kg  + (1|arena), 
            data=ppdata, control = control_ml, na.action = "na.fail")

mod_tabs_site <- lmer(s_tablestress_p ~ study_site  + (1|arena), 
            data=ppdata, control = control_ml, na.action = "na.fail")

models <- list(mod_tabs_cond, mod_tabs_weight, mod_tabs_site)
model.names <- c("Condition", "Body weight", "Study site")
aictab(cand.set = models, modnames = model.names)

#Sum section
mod_sums_weight <- lmer(s_sumsecto_p ~ body_weight_kg  + (1|arena), 
            data=ppdata, control = control_ml, na.action = "na.fail")

mod_sums_sex <- lmer(s_sumsecto_p ~ sex  + (1|arena), 
            data=ppdata, control = control_ml, na.action = "na.fail")

models <- list(mod_sums_weight, mod_sums_sex)
model.names <- c("Body weight", "Sex")
aictab(cand.set = models, modnames = model.names)

#Total act
mod_tots_cond <- lmer(s_totalacto_p ~ cond_n + (1|arena), 
            data=ppdata, control = control_ml, na.action = "na.fail")

mod_tots_weight <- lmer(s_totalacto_p ~ body_weight_kg  + (1|arena), 
            data=ppdata, control = control_ml, na.action = "na.fail")

mod_tots_sex <- lmer(s_totalacto_p ~ sex  + (1|arena), 
            data=ppdata, control = control_ml, na.action = "na.fail")

models <- list(mod_tots_cond, mod_tots_weight, mod_tots_sex)
model.names <- c("Condition", "Body weight", "Sex")
aictab(cand.set = models, modnames = model.names)

```

#Plasticity x weight plots
```{r}

#Food occ v weight
food <- ggplot(data=ppdata, aes(x = body_weight_kg, y = s_foodo_p)) + 
 geom_point() +
  geom_smooth(method = "lm", se = TRUE, colour = "#440154", fill = "#82568D") +
  theme_minimal() +
    theme(legend.position = "right", 
        panel.grid = element_blank(),
        axis.line.x = element_line(color = "black"),
        axis.line.y = element_line(color = "black"),
        text = element_text(size=19),  
       axis.text.x = element_text(size=19), 
        plot.margin = margin(0.5, 0.5, 0.5, 1, "cm")) + 
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey") +  
  scale_y_continuous(breaks=seq(-30,10,10),labels=seq(-30,10,10), limits=c(-30,10)) +
  xlab("") + 
    ylab(expression(paste("Adjusted food \ninteraction occurrences slope", italic(" (n)")))) + labs("") +
  annotate("text", x =1.6, y =10, label = "A", size = 10, hjust = 1)
print(food)

#Legs occ v weight
legso_weight <- ggplot(data=ppdata, aes(x = body_weight_kg, y = s_legso_p)) + 
 geom_point() +
  geom_smooth(method = "lm", se = TRUE, colour = "#472D7B", fill = "#8473A7") +
  theme_minimal() +
    theme(legend.position = "right", 
        panel.grid = element_blank(),
        axis.line.x = element_line(color = "black"),
        axis.line.y = element_line(color = "black"),
        text = element_text(size=19),  
       axis.text.x = element_text(size=19), 
        plot.margin = margin(0.5, 0.5, 0.5, 0.5, "cm")) + 
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey") +  
  scale_y_continuous(breaks=seq(-100,25,25),labels=seq(-100,25,25), limits=c(-100,25)) +
  xlab("") + 
    ylab(expression(paste("Adjusted hindleg stand \noccurrences slope", italic(" (n)"),))) +
  labs("") +
  annotate("text", x =1.6, y =25, label = "B", size = 10, hjust = 1)
print(legso_weight)

#Legs dur v weight
legsd_weight <- ggplot(data=ppdata, aes(x = body_weight_kg, y = s_legsd_p)) + 
 geom_point() +
  geom_smooth(method = "lm", se = TRUE, colour = "#3B528B", fill = "#7C8CB2") +
  theme_minimal() +
    theme(legend.position = "right", 
        panel.grid = element_blank(),
        axis.line.x = element_line(color = "black"),
        axis.line.y = element_line(color = "black"),
        text = element_text(size=19),  
       axis.text.x = element_text(size=19), 
        plot.margin = margin(0.5, 0.5, 0.5, 0.5, "cm")) + 
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey") +  
  scale_y_continuous(breaks=seq(-40,50,20),labels=seq(-40,50,20), limits=c(-40, 50)) +
  xlab("") + 
    ylab("Adjusted hindleg stand \nduration slope (s)") + labs("") +
  annotate("text", x =1.6, y =50, label = "C", size = 10, hjust = 1)
print(legsd_weight)

#Move dur v weight
moveweight <- ggplot(data=ppdata, aes(x = body_weight_kg, y = s_moved_p)) + 
 geom_point() +
  geom_smooth(method = "lm", se = TRUE, colour = "#2C728E", fill = "#72A1B4") +
  theme_minimal() +
    theme(legend.position = "right", 
        panel.grid = element_blank(),
        axis.line.x = element_line(color = "black"),
        axis.line.y = element_line(color = "black"),
        text = element_text(size=19),  
       axis.text.x = element_text(size=19), 
        plot.margin = margin(0.5, 0.5, 0.5, 0.5, "cm")) +    
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey") +  
  scale_y_continuous(breaks=seq(-150,50,50),labels=seq(-150,50,50), limits=c(-170, 70)) +
  xlab("") + 
    ylab("Adjusted quadrupedal \nmovement duration slope (s)") +
  annotate("text", x =1.6, y =70, label = "D", size = 10, hjust = 1)
print(moveweight)

#Move occ v weight
moveo_weight <- ggplot(data=ppdata, aes(x = body_weight_kg, y = s_moveo_p)) + 
 geom_point() +
  geom_smooth(method = "lm", se = TRUE, colour = "#21908C", fill = "#6BB5B2") +
  theme_minimal() +
    theme(legend.position = "right", 
        panel.grid = element_blank(),
        axis.line.x = element_line(color = "black"),
        axis.line.y = element_line(color = "black"),
        text = element_text(size=19),  
       axis.text.x = element_text(size=19), 
        plot.margin = margin(0.5, 0.5, 0.5, 0.5, "cm")) + 
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey") +  
  scale_y_continuous(breaks=seq(-150,50,50),labels=seq(-150,50,50), limits=c(-150, 50)) +
  xlab("") + 
    ylab(expression(paste("Adjusted quadrupedal \nmovement occurrences slope ", italic(" (n)"),))) +
  annotate("text", x =1.6, y =50, label = "E", size = 10, hjust = 1)
print(moveo_weight)

#Sums sect v weight
sum_weight <- ggplot(data=ppdata, aes(x = body_weight_kg, y = s_sumsecto_p)) + 
 geom_point() +
  geom_smooth(method = "lm", se = TRUE, colour = "#27AD81", fill = "#6FC8AB") +
  theme_minimal() +
    theme(legend.position = "right", 
        panel.grid = element_blank(),
        axis.line.x = element_line(color = "black"),
        axis.line.y = element_line(color = "black"),
        text = element_text(size=19),  
       axis.text.x = element_text(size=19), 
        plot.margin = margin(0.5, 0.5, 0.5, 0.5, "cm")) + 
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey") +  
  scale_y_continuous(breaks=seq(-40,0,20),labels=seq(-40,0,20), limits=c(-40, 0)) +
  xlab("") + 
    ylab(expression(paste("Adjusted total\n arena section visits slope", italic(" (n)"),))) +
  annotate("text", x =1.6, y =0, label = "F", size = 10, hjust = 1)
#print(sum_weight)

#Total act v weight
tot_weight <- ggplot(data=ppdata, aes(x = body_weight_kg, y = s_totalacto_p)) + 
 geom_point() +
  geom_smooth(method = "lm", se = TRUE, colour = "#5DC863", fill = "#93DA97") +
  theme_minimal() +
    theme(legend.position = "right", 
        panel.grid = element_blank(),
        axis.line.x = element_line(color = "black"),
        axis.line.y = element_line(color = "black"),
        text = element_text(size=19),  
       axis.text.x = element_text(size=19), 
        plot.margin = margin(0.5, 0.5, 0.5, 1, "cm")) +  
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey") +  
  scale_y_continuous(breaks=seq(-300,10,50),labels=seq(-300,10,50), limits=c(-300, 20)) +
  xlab("") + 
    ylab(expression(paste("Adjusted total\n activity slope", italic(" (n)"),))) +
  annotate("text", x =1.6, y =0, label = "G", size = 10, hjust = 1)
print(tot_weight)

#Combine plots
plot <- ggarrange(food, legso_weight, legsd_weight, moveweight, moveo_weight, sum_weight, tot_weight, nrow=3,ncol=3)
plot <- annotate_figure(plot,
        bottom = text_grob("Body weight (kg)", size = 25))
print(plot)
ggsave(filename = "slope - weight plots.jpeg", plot, width=350, height= 400, units = "mm")

```

#Plasticity x site plots
```{r}

custom_colors <- c("Goorooyarroo" = "#2D718E", "Mulligans Flat" = "#B8DE29")
#custom_colors <- c("Goorooyarroo" = "#674A8E", "Mulligans Flat" = "#2D718E")

#Food interaction dur violin
s_fd_site <- ggplot(data=ppdata, aes(x = study_site, y = s_foodd_p, fill=study_site)) + 
 geom_violin() +
  geom_boxplot(col="grey20", fill="grey20", width=0.02, lwd=0.5) +
  stat_summary(fun="median", geom="point", col="black", fill="white", 
               shape=21, size=3, stroke=1) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey") +  
  scale_fill_manual (values = c("#6B9BAF","#CDE868")) +
  theme_minimal() +
  theme(legend.position = "none", 
        panel.grid = element_blank(),
        axis.line.x = element_line(color = "black"),
        axis.line.y = element_line(color = "black"),
        text = element_text(size=19),  
       axis.text.x = element_text(size=18), 
        plot.margin = margin(0.5, 0.5, 0.5, 0.5, "cm")) +    
  scale_y_continuous(breaks=seq(-100,60,20),labels=seq(-100,60,20), limits=c(-100,60)) +
  xlab("") + 
  ylab("Adjusted food \ninteraction duration slope (s)") + labs("") +
  annotate("text", x ="Mulligans Flat", y =60, label = "Ai", size = 10, hjust = -2)
print(s_fd_site)

#Food dur vs site smooth
food <- ggplot(data=d2pred, aes(x = assay_number, y =foodd_p, color = study_site, fill = study_site)) + 
 geom_point() +
  geom_smooth(method = "lm", se = TRUE, aes(fill = study_site)) +
  theme_minimal() +
    theme(legend.position = "none", 
        panel.grid = element_blank(),
        axis.line.x = element_line(color = "black"),
        axis.line.y = element_line(color = "black"),
        text = element_text(size=19),  
       axis.text.x = element_text(size=18), 
        plot.margin = margin(0.5, 0.5, 0.5, 0.5, "cm")) + 
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey") + 
  scale_y_continuous(breaks=seq(-50,100,25),labels=seq(-50,100,25), limits=c(-50,110)) +
  scale_x_continuous(breaks = c(min(d2pred$assay_number), max(d2pred$assay_number))) +
  xlab("") + 
  ylab("Adjusted food \ninteraction duration (s)") +
  labs(color = "Study site", fill = "Study site") +
  scale_color_manual(values = custom_colors) +
  scale_fill_manual(values = custom_colors) +
  annotate("text", x =max(d2pred$assay_number), y =110, label = "Aii", size = 10, hjust = 1)
print(food)

#Jump + climb dur violin
s_jd_site <- ggplot(data=ppdata, aes(x = study_site, y = s_jumpd_p, fill=study_site)) + 
 geom_violin() +
  geom_boxplot(col="grey20", fill="grey20", width=0.02, lwd=0.5) +
  stat_summary(fun="median", geom="point", col="black", fill="white", 
               shape=21, size=3, stroke=1) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey") +
 scale_fill_manual (values = c("#6B9BAF","#CDE868")) +
  theme_minimal() +
  theme(legend.position = "none", 
        panel.grid = element_blank(),
        axis.line.x = element_line(color = "black"),
        axis.line.y = element_line(color = "black"),
        text = element_text(size=19),  
       axis.text.x = element_text(size=18), 
        plot.margin = margin(0.5, 0.5, 0.5, 0.5, "cm")) + 
  scale_y_continuous(breaks=seq(-150,50,50),labels=seq(-150,50,50), limits=c(-150,50)) +
  xlab("") + 
  ylab("Adjusted jumping + \nclimbing duration slope (s)") +
  annotate("text", x ="Mulligans Flat", y =50, label = "Bi", size = 10, hjust = -2)
print(s_jd_site)

#Jump climb dur smooth
jump <- ggplot(data=d2pred, aes(x = assay_number, y =jumpd_p, color = study_site, fill = study_site)) + 
 geom_point() +
  geom_smooth(method = "lm", se = TRUE, aes(fill = study_site)) +
  theme_minimal() +
    theme(legend.position = "none", 
        panel.grid = element_blank(),
        axis.line.x = element_line(color = "black"),
        axis.line.y = element_line(color = "black"),
        text = element_text(size=19),  
       axis.text.x = element_text(size=18), 
        plot.margin = margin(0.5, 0.5, 0.5, 1, "cm")) + 
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey") +  
  scale_x_continuous(breaks = c(min(d2pred$assay_number), max(d2pred$assay_number))) +
  xlab("") + 
  ylab("Adjusted jumping + \nclimbing duration (s)") +
  labs(color = "Study site", fill = "Study site") +
  scale_color_manual(values = custom_colors) +
  scale_fill_manual(values = custom_colors) +
  annotate("text", x =max(d2pred$assay_number), y =200, label = "Bii", size = 10, hjust = 1)
print(jump)

#Jump + climb occurrences violin
s_jo_site <- ggplot(data=ppdata, aes(x = study_site, y = s_jumpo_p, fill=study_site)) + 
 geom_violin() +
  geom_boxplot(col="grey20", fill="grey20", width=0.02, lwd=0.5) +
  stat_summary(fun="median", geom="point", col="black", fill="white", 
               shape=21, size=3, stroke=1) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey") +  # Add horizontal line at y = 0
  scale_fill_manual (values = c("#6B9BAF","#CDE868")) +
  theme_minimal() +
  theme(legend.position = "none", 
        panel.grid = element_blank(),
        axis.line.x = element_line(color = "black"),
        axis.line.y = element_line(color = "black"),
        text = element_text(size=19),  
       axis.text.x = element_text(size=18), 
        plot.margin = margin(0.5, 0.5, 0.5, 1, "cm")) + 
  scale_y_continuous(breaks=seq(-60,0,20),labels=seq(-60,0,20), limits=c(-60,0)) +
  xlab(" ") + 
  ylab(expression(paste("Adjusted jumping + \nclimbing occurrences slope  ", italic("(n)"))))+
  annotate("text", x ="Mulligans Flat", y =0, label = "Ci", size = 10, hjust = -2)
print(s_jo_site)

#Jump climb occ smooth
jump2 <- ggplot(data=d2pred, aes(x = assay_number, y =jumpo_p, color = study_site, fill = study_site)) + 
 geom_point() +
  geom_smooth(method = "lm", se = TRUE, aes(fill = study_site)) +
  theme_minimal() +
    theme(legend.position = "none", 
        panel.grid = element_blank(),
        axis.line.x = element_line(color = "black"),
        axis.line.y = element_line(color = "black"),
        text = element_text(size=19),  
       axis.text.x = element_text(size=18), 
        plot.margin = margin(0.5, 0.5, 0.5, 1, "cm")) + 
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey") +
  scale_x_continuous(breaks = c(min(d2pred$assay_number), max(d2pred$assay_number))) +
  xlab("") + 
  ylab(expression(paste("Adjusted jumping + \nclimbing occurrences ", italic("(n)"), ""))) +
  labs(color = "Study site", fill = "Study site") +
  scale_color_manual(values = custom_colors) +
  scale_fill_manual(values = custom_colors) +
  annotate("text", x =max(d2pred$assay_number), y =80, label = "Cii", size = 10, hjust = 1)
print(jump2)

#Latency to emerge violin
s_e_site <- ggplot(data=ppdata, aes(x = study_site, y = s_emerge_p, fill=study_site)) + 
 geom_violin() +
  geom_boxplot(col="grey20", fill="grey20", width=0.02, lwd=0.5) +
  stat_summary(fun="median", geom="point", col="black", fill="white", 
               shape=21, size=3, stroke=1) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey") +  # Add horizontal line at y = 0
  scale_y_continuous(breaks=seq(-10,50,10),labels=seq(-10,50,10), limits=c(-10,50)) +
  scale_fill_manual (values = c("#6B9BAF","#CDE868")) +
  theme_minimal() +
  theme(legend.position = "none", 
        panel.grid = element_blank(),
        axis.line.x = element_line(color = "black"),
        axis.line.y = element_line(color = "black"),
        text = element_text(size=19),  
       axis.text.x = element_text(size=18), 
        plot.margin = margin(0.5, 0.5, 0.5, 0.5, "cm")) +  
  xlab("") + 
  ylab("Adjusted latency to \nemerge slope (s)") + labs("") +
  annotate("text", x ="Mulligans Flat", y =50, label = "Di", size = 10, hjust = -2)
print(s_e_site)

#Latency to emerge smooth
emerge <- ggplot(data=d2pred, aes(x = assay_number, y = emerge_p, color = study_site, fill = study_site)) + 
 geom_point() +
  geom_smooth(method = "lm", se = TRUE, aes(fill = study_site)) +
  theme_minimal() +
    theme(legend.position = "none", 
        panel.grid = element_blank(),
        axis.line.x = element_line(color = "black"),
        axis.line.y = element_line(color = "black"),
        text = element_text(size=19),  
       axis.text.x = element_text(size=18), 
        plot.margin = margin(0.5, 0.5, 0.5, 0.5, "cm")) +  
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey") +
  scale_x_continuous(breaks = c(min(d2pred$assay_number), max(d2pred$assay_number))) +
  scale_y_continuous(breaks=seq(-10,50,10),labels=seq(-10,50,10), limits=c(-10,50)) +
  xlab("") + 
  ylab("Adjusted latency \nto emerge (s)") +
  labs(color = "Study site", fill = "Study site") +
  scale_color_manual(values = custom_colors) +
  scale_fill_manual(values = custom_colors) +
  annotate("text", x =max(d2pred$assay_number), y =50, label = "Dii", size = 10, hjust = 1)
print(emerge)

#Perching duration violin
s_pd_site <- ggplot(data=ppdata, aes(x = study_site, y = s_perchd_p, fill=study_site)) + 
 geom_violin() +
  geom_boxplot(col="grey20", fill="grey20", width=0.02, lwd=0.5) +
  stat_summary(fun="median", geom="point", col="black", fill="white",  
               shape=21, size=3, stroke=1) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey") +  # Add horizontal line at y = 0
  scale_fill_manual (values = c("#6B9BAF","#CDE868")) +
  theme_minimal() +
  theme(legend.position = "none", 
        panel.grid = element_blank(),
        axis.line.x = element_line(color = "black"),
        axis.line.y = element_line(color = "black"),
        text = element_text(size=19), 
       axis.text.x = element_text(size=18),
        plot.margin = margin(0.5, 0.5, 0.5, 0.5, "cm")) +
  scale_y_continuous(breaks=seq(-200,200,100),labels=seq(-200,200,100), limits=c(-200,200)) +
  xlab(" ") + 
  ylab("Adjusted perching \nduration slope (s)") + labs("") +
  annotate("text", x ="Mulligans Flat", y =200, label = "Ei", size = 10, hjust = -2)
print(s_pd_site)

#Perch dur smooth
perch <- ggplot(data=d2pred, aes(x = assay_number, y =perchd_p, color = study_site, fill = study_site)) + 
 geom_point() +
  geom_smooth(method = "lm", se = TRUE, aes(fill = study_site)) +
  theme_minimal() +
    theme(legend.position = "none", 
        panel.grid = element_blank(),
        axis.line.x = element_line(color = "black"),
        axis.line.y = element_line(color = "black"),
        text = element_text(size=19), 
       axis.text.x = element_text(size=18),
        plot.margin = margin(0.5, 0.5, 0.5, 0.5, "cm")) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey") +
  scale_x_continuous(breaks = c(min(d2pred$assay_number), max(d2pred$assay_number))) +
  xlab("") + 
  ylab("Adjusted perching \nduration (s)") +
  labs(color = "Study site", fill = "Study site") +
  scale_color_manual(values = custom_colors) +
  scale_fill_manual(values = custom_colors) +
  annotate("text", x =max(d2pred$assay_number), y =300, label = "Eii", size = 10, hjust = 1)
 print(perch)
 
#Perching occurrences violin
s_po_site <- ggplot(data=ppdata, aes(x = study_site, y = s_percho_p, fill=study_site)) + 
 geom_violin() +
  geom_boxplot(col="grey20", fill="grey20", width=0.02, lwd=0.5) +
  stat_summary(fun="median", geom="point", col="black", fill="white", 
               shape=21, size=3, stroke=1) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey") +  # Add horizontal line at y = 0
  scale_fill_manual (values = c("#6B9BAF","#CDE868")) +
  theme_minimal() +
  theme(legend.position = "none", 
        panel.grid = element_blank(),
        axis.line.x = element_line(color = "black"),
        axis.line.y = element_line(color = "black"),
        text = element_text(size=19), 
       axis.text.x = element_text(size=18),
        plot.margin = margin(0.5, 0.5, 0.5, 1, "cm")) +
  xlab("") + 
  ylab(expression(paste("Adjusted perching \noccurrences slope", italic(" (n)")))) +
  annotate("text", x ="Mulligans Flat", y =5, label = "Fi", size = 10, hjust = -2)
print(s_po_site)

#Perch occ smooth 
perch2 <- ggplot(data=d2pred, aes(x = assay_number, y =percho_p, color = study_site, fill = study_site)) + 
 geom_point() +
  geom_smooth(method = "lm", se = TRUE, aes(fill = study_site)) +
  theme_minimal() +
    theme(legend.position = "none", 
        panel.grid = element_blank(),
        axis.line.x = element_line(color = "black"),
        axis.line.y = element_line(color = "black"),
        text = element_text(size=19), 
       axis.text.x = element_text(size=18),
        plot.margin = margin(0.5, 0.5, 0.5, 1, "cm")) +
  scale_x_continuous(breaks = c(min(d2pred$assay_number), max(d2pred$assay_number))) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey") +
  xlab("") + 
  ylab(expression(paste("Adjusted perching\n occurrences ", italic("(n)"), ""))) +
  labs(color = "Study site", fill = "Study site") +
  scale_color_manual(values = custom_colors) +
  scale_fill_manual(values = custom_colors) +
  annotate("text", x =max(d2pred$assay_number), y =30, label = "Fii", size = 10, hjust = 1)
print(perch2)
 
#Movement duration v site
s_md_site <- ggplot(data=ppdata, aes(x = study_site, y = s_moved_p, fill=study_site)) + 
 geom_violin() +
  geom_boxplot(col="grey20", fill="grey20", width=0.02, lwd=0.5) +
  stat_summary(fun="median", geom="point", col="black", fill="white", 
               shape=21, size=3, stroke=1) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey") +  # Add horizontal line at y = 0
  scale_fill_manual (values = c("#6B9BAF","#CDE868")) +
  theme_minimal() +
  theme(legend.position = "none", 
        panel.grid = element_blank(),
        axis.line.x = element_line(color = "black"),
        axis.line.y = element_line(color = "black"),
        text = element_text(size=19), 
       axis.text.x = element_text(size=18),
        plot.margin = margin(0.5, 0.5, 0.5, 0.5, "cm")) +
  xlab("") + 
  ylab("Adjusted quadrupedal \nmovement duration slope (s)") + labs("") +
    annotate("text", x ="Mulligans Flat", y =50, label = "Gi", size = 10, hjust = -2)
print(s_md_site)

#Move dur vs site 
move <- ggplot(data=d2pred, aes(x = assay_number, y =moved_p, color = study_site, fill = study_site)) + 
 geom_point() +
  geom_smooth(method = "lm", se = TRUE, aes(fill = study_site)) +
  theme_minimal() +
    theme(legend.position = "none", 
        panel.grid = element_blank(),
        axis.line.x = element_line(color = "black"),
        axis.line.y = element_line(color = "black"),
        text = element_text(size=19), 
       axis.text.x = element_text(size=18),
        plot.margin = margin(0.5, 0.5, 0.5, 0.5, "cm")) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey") +
  scale_x_continuous(breaks = c(min(d2pred$assay_number), max(d2pred$assay_number))) +
  xlab("") + 
  ylab("Adjusted quadrupedal\n movement duration (s)") +
  labs(color = "Study site", fill = "Study site") +
  scale_color_manual(values = custom_colors) +
  scale_fill_manual(values = custom_colors) +
  annotate("text", x =max(d2pred$assay_number), y =160, label = "Gii", size = 10, hjust = 1)
print(move)

#Bag stress rating violin
s_bs_site <- ggplot(data=ppdata, aes(x = study_site, y = s_bagstress_p, fill=study_site)) + 
 geom_violin() +
  geom_boxplot(col="grey20", fill="grey20", width=0.02, lwd=0.5) +
  stat_summary(fun="median", geom="point", col="black", fill="white", 
               shape=21, size=3, stroke=1) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey") +  # Add horizontal line at y = 0
  scale_fill_manual (values = c("#6B9BAF","#CDE868")) +
  theme_minimal() +
  theme(legend.position = "none", 
        panel.grid = element_blank(),
        axis.line.x = element_line(color = "black"),
        axis.line.y = element_line(color = "black"),
        text = element_text(size=19), 
       axis.text.x = element_text(size=18),
        plot.margin = margin(0.5, 0.5, 0.5, 0.5, "cm")) +
  xlab("") + 
  ylab("Adjusted bag stress\n rating slope") + labs("") +
  annotate("text", x ="Mulligans Flat", y =0.2, label = "Hi", size = 10, hjust = -2)
print(s_bs_site)

#Bag stress vs site smooth
bag <- ggplot(data=d2pred, aes(x = assay_number, y =bagstress_p, color = study_site, fill = study_site)) + 
 geom_point() +
  geom_smooth(method = "lm", se = TRUE, aes(fill = study_site)) +
  theme_minimal() +
    theme(legend.position = "none", 
        panel.grid = element_blank(),
        axis.line.x = element_line(color = "black"),
        axis.line.y = element_line(color = "black"),
        text = element_text(size=19), 
       axis.text.x = element_text(size=18),
        plot.margin = margin(0.5, 0.5, 0.5, 0.5, "cm")) +
  #geom_hline(yintercept = 0, linetype = "dashed", color = "grey") +
  scale_y_continuous(breaks=seq(0.9,1.3,0.1),labels=seq(0.9,1.3,0.1), limits=c(0.9,1.3))+
  scale_x_continuous(breaks = c(min(d2pred$assay_number), max(d2pred$assay_number))) +
  xlab("") + 
  ylab("Adjusted bag \nstress rating") +
  labs(color = "Study site", fill = "Study site") +
  scale_color_manual(values = custom_colors) +
  scale_fill_manual(values = custom_colors) +
  annotate("text", x =max(d2pred$assay_number), y =1.3, label = "Hii", size = 10, hjust = 1)
print(bag)

#Table stress site violin
s_ts_site <- ggplot(data=ppdata, aes(x = study_site, y = s_tablestress_p, fill=study_site)) + 
 geom_violin() +
  geom_boxplot(col="grey20", fill="grey20", width=0.02, lwd=0.5) +
  stat_summary(fun="median", geom="point", col="black", fill="white",  
               shape=21, size=3, stroke=1) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey") +  # Add horizontal line at y = 0
 scale_fill_manual (values = c("#6B9BAF","#CDE868")) +
  theme_minimal() +
  theme(legend.position = "none", 
        panel.grid = element_blank(),
        axis.line.x = element_line(color = "black"),
        axis.line.y = element_line(color = "black"),
        text = element_text(size=19), 
       axis.text.x = element_text(size=18),
        plot.margin = margin(0.5, 0.5, 0.5, 0.5, "cm")) + 
  scale_y_continuous(breaks=seq(-0.1,0.2,0.1),labels=seq(-0.1,0.2,0.1), limits=c(-0.1,0.2)) +
  xlab("") + 
  ylab("Adjusted table \nstress rating slope") + labs("") +
  annotate("text", x ="Mulligans Flat", y =0.2, label = "Ii", size = 10, hjust = -3)
print(s_ts_site)

#Table stress vs site 
table <- ggplot(data=d2pred, aes(x = assay_number, y =tablestress_p, color = study_site, fill = study_site)) + 
 geom_point() +
  geom_smooth(method = "lm", se = TRUE, aes(fill = study_site)) +
  theme_minimal() +
    theme(legend.position = "none", 
        panel.grid = element_blank(),
        axis.line.x = element_line(color = "black"),
        axis.line.y = element_line(color = "black"),
        text = element_text(size=19), 
       axis.text.x = element_text(size=18),
        plot.margin = margin(0.5, 0.5, 0.5, 0.5, "cm")) +
   scale_x_continuous(breaks = c(min(d2pred$assay_number), max(d2pred$assay_number))) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey") +
  xlab("") + 
  ylab("Adjusted table\n stress rating") +
  labs(color = "Study site", fill = "Study site") +
  scale_color_manual(values = custom_colors) +
  scale_fill_manual(values = custom_colors) +
  annotate("text", x =max(d2pred$assay_number), y =2, label = "Iii", size = 10, hjust = 1)
print(table)

#Combine plots
plot <- ggarrange(s_fd_site,	s_jd_site,	s_jo_site,
                  food,	jump,	jump2,
                  s_e_site,	s_pd_site, s_po_site,
                  emerge, perch, perch2,
                  s_md_site, s_bs_site, s_ts_site,
                  move, bag, table, 
                  nrow=6,ncol=3)
ggsave(filename = "plasticity plots.jpeg", plot, width=450, height= 700, units = "mm")

```

# **Session information**
```{r}
# Display version information for R, OS, and packages
sessionInfo()
```

R version 4.3.3 (2024-02-29 ucrt)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19045)

Matrix products: default

locale:
[1] LC_COLLATE=English_Australia.utf8 
[2] LC_CTYPE=English_Australia.utf8   
[3] LC_MONETARY=English_Australia.utf8
[4] LC_NUMERIC=C                      
[5] LC_TIME=English_Australia.utf8    

time zone: Australia/Sydney
tzcode source: internal

attached base packages:
[1] stats     graphics 
[3] grDevices utils    
[5] datasets  methods  
[7] base     

other attached packages:
 [1] scales_1.3.0       
 [2] glmmTMB_1.1.9      
 [3] AICcmodavg_2.3-3   
 [4] broom.mixed_0.2.9.5
 [5] kableExtra_1.4.0   
 [6] knitr_1.45         
 [7] sjlabelled_1.2.0   
 [8] sjmisc_2.8.9       
 [9] sjPlot_2.8.15      
[10] broom_1.0.5        
[11] corrplot_0.92      
[12] viridis_0.6.5      
[13] viridisLite_0.4.2  
[14] lubridate_1.9.3    
[15] forcats_1.0.0      
[16] stringr_1.5.1      
[17] dplyr_1.1.4        
[18] purrr_1.0.2        
[19] tidyr_1.3.1        
[20] tibble_3.2.1       
[21] tidyverse_2.0.0    
[22] rstudioapi_0.16.0  
[23] readxl_1.4.3       
[24] sf_1.0-16          
[25] readr_2.1.5        
[26] MuMIn_1.47.5       
[27] janitor_2.2.0      
[28] gtools_3.9.5       
[29] ggpubr_0.6.0       
[30] ggplot2_3.5.0      
[31] lmerTest_3.1-3     
[32] lme4_1.1-35.1      
[33] Matrix_1.6-5       

loaded via a namespace (and not attached):
 [1] DBI_1.2.2          
 [2] gridExtra_2.3      
 [3] rlang_1.1.3        
 [4] magrittr_2.0.3     
 [5] furrr_0.3.1        
 [6] snakecase_0.11.1   
 [7] e1071_1.7-14       
 [8] compiler_4.3.3     
 [9] mgcv_1.9-1         
[10] systemfonts_1.0.6  
[11] vctrs_0.6.5        
[12] pkgconfig_2.0.3    
[13] fastmap_1.1.1      
[14] backports_1.4.1    
[15] labeling_0.4.3     
[16] utf8_1.2.4         
[17] rmarkdown_2.26     
[18] tzdb_0.4.0         
[19] nloptr_2.0.3       
[20] ragg_1.3.0         
[21] xfun_0.43          
[22] ggeffects_1.5.1    
[23] VGAM_1.1-10        
[24] parallel_4.3.3     
[25] R6_2.5.1           
[26] stringi_1.8.3      
[27] parallelly_1.37.1  
[28] car_3.1-2          
[29] boot_1.3-29        
[30] cellranger_1.1.0   
[31] numDeriv_2016.8-1.1
[32] estimability_1.5   
[33] Rcpp_1.0.12        
[34] modelr_0.1.11      
[35] pacman_0.5.1       
[36] splines_4.3.3      
[37] timechange_0.3.0   
[38] tidyselect_1.2.1   
[39] abind_1.4-5        
[40] yaml_2.3.8         
[41] TMB_1.9.11         
[42] codetools_0.2-19   
[43] listenv_0.9.1      
[44] lattice_0.22-5     
[45] withr_3.0.0        
[46] bayestestR_0.13.2  
[47] coda_0.19-4.1      
[48] evaluate_0.23      
[49] survival_3.5-8     
[50] future_1.33.2      
[51] units_0.8-5        
[52] proxy_0.4-27       
[53] xml2_1.3.6         
[54] pillar_1.9.0       
[55] carData_3.0-5      
[56] KernSmooth_2.23-22 
[57] stats4_4.3.3       
[58] insight_0.19.10    
[59] generics_0.1.3     
[60] hms_1.1.3          
[61] munsell_0.5.0      
[62] minqa_1.2.6        
[63] globals_0.16.3     
[64] xtable_1.8-4       
[65] unmarked_1.4.1     
[66] class_7.3-22       
[67] glue_1.7.0         
[68] emmeans_1.10.0     
[69] tools_4.3.3        
[70] ggsignif_0.6.4     
[71] mvtnorm_1.2-4      
[72] cowplot_1.1.3      
[73] grid_4.3.3         
[74] colorspace_2.1-0   
[75] nlme_3.1-164       
[76] performance_0.11.0 
[77] cli_3.6.2          
[78] textshaping_0.3.7  
[79] fansi_1.0.6        
[80] svglite_2.1.3      
[81] sjstats_0.18.2     
[82] gtable_0.3.4       
[83] rstatix_0.7.2      
[84] digest_0.6.35      
[85] classInt_0.4-10    
[86] farver_2.1.1       
[87] htmltools_0.5.8    
[88] lifecycle_1.0.4    
[89] MASS_7.3-60.0.1    

```