---
title: "Thermal traits across an elevational gradient"
author: "Esme Ashe-Jepson"
format: html
editor: visual
editor_options: 
  chunk_output_type: console
execute:
  message: false
  warning: false
  options:
    root.dir: "C:/Users/esa38zy/Documents/Research/Alps/Data/R"
bibliography: references.bib
---

# Preparation

First we can load our datasets that we have prepared in the separate 'Preparing dataset.qmd' file.

```{r}
#| label: load datasets

setwd("C:/Users/esa38zy/Documents/Research/Alps/Data/R")

load("Buff_clean.RData") # pre-formatted for analysis
load("CTmax_clean.RData") # CTmax data
```

Next we can load all the necessary packages.

```{r}
#| label: load packages

library(lme4) # for LME
library(ggplot2) # for plotting
library(dplyr) # for formatting data 
library(performance) # for checking models
library(car) # for term significance
library(lmerTest) # for checking LME
library(dplyr) # for data management
library(ggeffects) # for plotting models
library(survival) # for survival curve analysis
library(survminer) # for plotting survival curves
library(coxme) # for random effect Cox models (survival curves)
library(MuMIn) # for conditional R^2
library(sjPlot) # for setting a plotting theme

set_theme(base=theme_light(),
          axis.linecolor.x = "black",
          axis.linecolor.y = "black",
          title.size=2.2,
          axis.textsize=1.2,
          axis.title.size = 1.4)
```

We can check our data to make sure it has loaded correctly.

```{r}
#| label: check data

glimpse(buff.clean)

glimpse(ct.clean)
```

# Data analysis

## 1.1 Thermoregulation across an elevational gradient

To answer our first research question, we want to test whether thermoregulation (body temperature or the relationship between body temperature and air temperature) changes across an elevational gradient and therefore shows evidence of adaptation to local climatic conditions. We also want to know whether this is impacted by functional traits.

Because we have a lot of variables we want to include in this model that use many different scales (e.g. wing length in cm, elevation in 1000s of metres), we can scale all our numerical variables first to improve model performance.

```{r}
#| label: thermoregulation ~ elevation scale variables

buff.clean$scale.B.temp <- scale(buff.clean$B.temp)
buff.clean$scale.A.temp <- scale(buff.clean$A.temp)
buff.clean$scale.Elevation <- scale(buff.clean$Elevation)
buff.clean$scale.Length.centered <- scale(buff.clean$Length.centered)
buff.clean$scale.Nir_value <- scale(buff.clean$Nir_value)
buff.clean$scale.Species_WL <- scale(buff.clean$Species_WL)
```

Now we can build our model. We start with all fixed effects, all interactions with air temperature, all interactions with elevation, and finally all three-way interactions with elevation and air temperature. We also want to account for non-independence in our data (in repeated samples of species and within sites), so we include a random intercept and slope for each.

```{r}
#| label: thermoregulation ~ elevation full model

lmer.1.1 <- lmer(scale.B.temp ~ scale.A.temp + scale.Elevation + scale.Length.centered + Condition + Sex + scale.Nir_value + scale.Species_WL +
                         scale.A.temp:scale.Elevation + scale.A.temp:scale.Length.centered + scale.A.temp:Condition + scale.A.temp:Sex + scale.A.temp:scale.Nir_value + scale.A.temp:scale.Species_WL +
                         scale.Elevation:scale.Length.centered + scale.Elevation:Condition + scale.Elevation:Sex + scale.Elevation:scale.Nir_value + scale.Elevation:scale.Species_WL +
                         scale.A.temp:scale.Elevation:scale.Length.centered + scale.A.temp:scale.Elevation:Condition + scale.A.temp:scale.Elevation:Sex + scale.A.temp:scale.Elevation:scale.Nir_value + scale.A.temp:scale.Elevation:scale.Species_WL +
                         (1 + scale.A.temp | Species) + (1 + scale.A.temp | Site),
                       data = buff.clean)

```

We next check our model performance, starting with our random effects.

```{r}
#| label: thermoregulation ~ elevation check full model

ranova(lmer.1.1)
```

Our random effects are not currently improving our model. We can start by removing the random slope for species.

```{r}
#| label: thermoregulation ~ elevation random effects 1

lmer.1.1.2 <- lmer(scale.B.temp ~ scale.A.temp + scale.Elevation + scale.Length.centered + Condition + Sex + scale.Nir_value + scale.Species_WL +
                           scale.A.temp:scale.Elevation + scale.A.temp:scale.Length.centered + scale.A.temp:Condition + scale.A.temp:Sex + scale.A.temp:scale.Nir_value + scale.A.temp:scale.Species_WL +
                           scale.Elevation:scale.Length.centered + scale.Elevation:Condition + scale.Elevation:Sex + scale.Elevation:scale.Nir_value + scale.Elevation:scale.Species_WL +
                           scale.A.temp:scale.Elevation:scale.Length.centered + scale.A.temp:scale.Elevation:Condition + scale.A.temp:scale.Elevation:Sex + scale.A.temp:scale.Elevation:scale.Nir_value + scale.A.temp:scale.Elevation:scale.Species_WL +
                           (1 | Species) + (1 + scale.A.temp | Site),
                         data = buff.clean)

ranova(lmer.1.1.2)
```

We still have a random effect that is not improving our model. Next let's remove the random slope for site.

```{r}
#| label: thermoregulation ~ elevation random effects 2

lmer.1.1.3 <- lmer(scale.B.temp ~ scale.A.temp + scale.Elevation + scale.Length.centered + Condition + Sex + scale.Nir_value + scale.Species_WL +
                           scale.A.temp:scale.Elevation + scale.A.temp:scale.Length.centered + scale.A.temp:Condition + scale.A.temp:Sex + scale.A.temp:scale.Nir_value + scale.A.temp:scale.Species_WL +
                           scale.Elevation:scale.Length.centered + scale.Elevation:Condition + scale.Elevation:Sex + scale.Elevation:scale.Nir_value + scale.Elevation:scale.Species_WL +
                           scale.A.temp:scale.Elevation:scale.Length.centered + scale.A.temp:scale.Elevation:Condition + scale.A.temp:scale.Elevation:Sex + scale.A.temp:scale.Elevation:scale.Nir_value + scale.A.temp:scale.Elevation:scale.Species_WL +
                           (1 | Species) + (1 | Site),
                         data = buff.clean)

ranova(lmer.1.1.3)

AIC(lmer.1.1, lmer.1.1.2, lmer.1.1.3) # all models perform similarly but lmer.1.1.3 is the simplest
```

Now we can finish checking our model now that our random effects structure is improved.

```{r}
#| label: thermoregulation ~ elevation check model

check_model(lmer.1.1.3)

qqnorm(resid(lmer.1.1.3))
qqline(resid(lmer.1.1.3))

plot(fitted(lmer.1.1.3), resid(lmer.1.1.3))

# Test for multicollinearity
lmer.1.1.3.vif <- lmer(scale.B.temp ~ scale.A.temp + scale.Elevation + scale.Length.centered + Condition + Sex + scale.Nir_value + scale.Species_WL +
                           (1 | Species) + (1 | Site),
                         data = buff.clean)
vif(lmer.1.1.3.vif)

r.squaredGLMM(lmer.1.1.3)
```

Our model performs well, and explains 74% of the variation in body temperature of the butterflies. Now we can test for term significance.

```{r}
#| label: thermoregulation ~ elevation term significance

Anova(lmer.1.1.3)
```

We can translate these values into real world numbers to get a sense of how butterfly body temperatures change across the elevational gradient.

```{r}
#| label: thermoregulation ~ elevation real world numbers

sd(buff.clean$Elevation)
sd(buff.clean$B.temp)
0.06996/340.3949 #  per metre
0.000205526*3.260723 # increase in body temperature per 1 metre increase in elevation
0.0006701634*1000 # increase across the 1000m elevation
```

We can see that though significant, the actual change in body temperature is very small.

Now we want to plot our results. @fig-thermoregulation-elevation shows how butterfly body temperatures change across elevations. @fig-thermoregulation-elevation-condition shows how butterfly body temperatures change across different wing conditions.

```{r}
#| label: fig-thermoregulation-elevation
#| fig-cap: | 
#|  The relationship between body temperature and air temperature across elevations.

preds <- ggpredict(lmer.1.1.3, 
                   terms = c("scale.A.temp [-1, 0, 1]", "scale.Elevation [-1, 0, 1]"))

ggplot(preds, aes(x = x, y = predicted, colour = group)) +
  geom_line(size = 1.2) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = group), alpha = 0.2, colour = NA) +
  labs(x = "Air temperature (scaled)",
       y = "Predicted body temperature (scaled)",
       colour = "Elevation \n(scaled)",
       fill = "Elevation \n(scaled)")
```

```{r}
#| label: fig-thermoregulation-elevation-condition
#| fig-cap: | 
#|  The relationship between body temperature and air temperature between butterflies of different conditions.

preds <- ggpredict(lmer.1.1.3, 
                   terms = c("scale.A.temp [-1, 0, 1]", "Condition"))

ggplot(preds, aes(x = x, y = predicted, colour = group)) +
  geom_line(size = 1.2) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = group), alpha = 0.2, colour = NA) +
  labs(x = "Air temperature (scaled)",
       y = "Predicted body temperature (scaled)",
       colour = "Condition",
       fill = "Condition")
```

Next we want to show our three-way interactions. These can be presented in several ways. For each case, I will show this in two combinations so that the patterns can be seen more clearly.

```{r}
#| label: fig-thermoregulation-elevation-spWL-1
#| fig-cap: | 
#|  The relationship between body temperature and air temperature across elevations between butterfly species of differing sizes.

preds_speciesWL <- ggpredict(
  lmer.1.1.3,
  terms = c("scale.A.temp [-1, 0, 1]",   
            "scale.Elevation [-1, 0, 1]", 
            "scale.Species_WL [-1, 0, 1]"))

ggplot(preds_speciesWL, aes(x = x, y = predicted, group = interaction(group, facet))) +
  geom_line(aes(colour = interaction(group, facet), linewidth = facet)) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = interaction(group, facet), group = interaction(group, facet)),
              alpha = 0.1, colour = NA) +
  facet_wrap(~ facet, labeller = labeller(facet = ~ paste("Species wing length:", .))) +
  scale_colour_manual(
    values = c(
      "-1.-1" = "#FF0000", "0.-1" = "#10E600", "1.-1" = "#1C1CFF",
      "-1.0"  = "#FF0000", "0.0"  = "#10E600", "1.0"  = "#1C1CFF",
      "-1.1"  = "#FF0000", "0.1"  = "#10E600", "1.1"  = "#1C1CFF"
    ),
    breaks = c("-1.0","0.0","1.0"),
    labels = c("-1","0","1"),
    name = "Elevation \n(scaled)"
  ) +
    scale_fill_manual(
    values = c(
      "-1.-1" = "#FF0000", "0.-1" = "#10E600", "1.-1" = "#1C1CFF",
      "-1.0"  = "#FF0000", "0.0"  = "#10E600", "1.0"  = "#1C1CFF",
      "-1.1"  = "#FF0000", "0.1"  = "#10E600", "1.1"  = "#1C1CFF"
    ),
    guide = "none") +
  scale_linewidth_manual(
    values = c(
      "-1" = 1,
      "0" = 2,
      "1" = 3
    ),
    guide = "none"
  ) +
  labs(x = "Air temperature (scaled)", y = "Predicted body temperature (scaled)")
```

Now we can present this in another way, with elevation as the facets and wing lengths as the grouping variable.

```{r}
#| label: fig-thermoregulation-elevation-spWL-2
#| fig-cap: | 
#|  The relationship between body temperature and air temperature across elevations between butterfly species of differing sizes.

preds_speciesWL.2 <- ggpredict(
  lmer.1.1.3,
  terms = c("scale.A.temp [-1, 0, 1]",   
            "scale.Species_WL [-1, 0, 1]", 
            "scale.Elevation [-1, 0, 1]"))

ggplot(preds_speciesWL.2, aes(x = x, y = predicted, group = interaction(group, facet))) +
  geom_line(aes(colour = facet, linewidth = interaction(group, facet))) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = facet, group = interaction(group, facet)),
              alpha = 0.1, colour = NA) +
  facet_wrap(~ facet, labeller = labeller(facet = ~ paste("Elevation:", .))) +
  scale_linewidth_manual(
    values = c(
      "-1.-1" = 1, "0.-1" = 2, "1.-1" = 3,
      "-1.0"  = 1, "0.0"  = 2, "1.0"  = 3,
      "-1.1"  = 1, "0.1"  = 2, "1.1"  = 3
    ),
    name = "Species mean \nwing length \n(scaled)",
        breaks = c("-1.0","0.0","1.0"),
    labels = c("-1","0","1")
  ) +
    scale_colour_manual(
    values = c(
      "-1" = "#FF0000",
      "0"  = "#10E600",
      "1"  = "#1C1CFF"
    ),
    guide = "none"
  ) +
  scale_fill_manual(
    values = c(
      "-1" = "#FF0000",
      "0"  = "#10E600",
      "1"  = "#1C1CFF"
    ),
    guide = "none"
  ) +
  labs(x = "Air temperature (scaled)", y = "Predicted body temperature (scaled)")
```

Now we can do the same for the NIR reflectance.

```{r}
#| label: fig-thermoregulation-elevation-nir-1
#| fig-cap: | 
#|  The relationship between body temperature and air temperature across elevations for butterflies of different colours.

preds_NIR <- ggpredict(
  lmer.1.1.3,
  terms = c("scale.A.temp [-1, 0, 1]",
            "scale.Elevation [-1, 0, 1]", 
            "scale.Nir_value [-1, 0, 1]"))

ggplot(preds_NIR, aes(x = x, y = predicted, group = interaction(group, facet))) +
  geom_line(aes(colour = interaction(group, facet)), linewidth = 1.2) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = interaction(group, facet)),
              alpha = 0.15, colour = NA) +
  facet_wrap(~ facet, labeller = labeller(facet = ~ paste("NIR reflection:", .))) +
  scale_colour_manual(
    values = c(
      "-1.-1" = "#800000", "0.-1" = "#006600", "1.-1" = "#0000A3",
      "-1.0"  = "#FF0000", "0.0"  = "#10E600", "1.0"  = "#1C1CFF",
      "-1.1"  = "#FFB8B8", "0.1"  = "#73FF73", "1.1"  = "#B8CBFF"
    ),
    breaks = c("-1.0","0.0","1.0"),
    labels = c("-1","0","1"),
    name = "Elevation \n(scaled)"
  ) +
  scale_fill_manual(
    values = c(
      "-1.-1" = "#800000", "0.-1" = "#006600", "1.-1" = "#0000A3",
      "-1.0"  = "#FF0000", "0.0"  = "#10E600", "1.0"  = "#1C1CFF",
      "-1.1"  = "#FFB8B8", "0.1"  = "#73FF73", "1.1"  = "#B8CBFF"
    ),
    guide = "none"
  ) +
  labs(x = "Air temperature (scaled)", y = "Predicted body temperature (scaled)")
```

And again, this can be presented in another way.

```{r}
#| label: fig-thermoregulation-elevation-nir-2
#| fig-cap: | 
#|  The relationship between body temperature and air temperature across elevations for butterflies of different colours.

preds_NIR.2 <- ggpredict(
  lmer.1.1.3,
  terms = c("scale.A.temp [-1, 0, 1]",
            "scale.Nir_value [-1, 0, 1]", 
            "scale.Elevation [-1, 0, 1]"))

ggplot(preds_NIR.2, aes(x = x, y = predicted, group = interaction(group, facet))) +
  geom_line(aes(colour = interaction(group, facet)), linewidth = 1.2) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = interaction(group, facet)),
              alpha = 0.15, colour = NA) +
  facet_wrap(~ facet, labeller = labeller(facet = ~ paste("Elevation:", .))) +
  scale_colour_manual(
    values = c(
      "-1.-1" = "#7F0000", "0.-1" = "#FF0000", "1.-1" = "#FFB8B8",
      "-1.0"  = "#003F00", "0.0"  = "#009100", "1.0"  = "#73FF73",
      "-1.1"  = "#00007F", "0.1"  = "#1915ED", "1.1"  = "#96B3FF"
    ),
    breaks = c("-1.0","0.0","1.0"),
    labels = c("-1","0","1"),
    name = "NIR \nreflectance \n(scaled)"
  ) +
  scale_fill_manual(
    values = c(
      "-1.-1" = "#7F0000", "0.-1" = "#FF0000", "1.-1" = "#FFB8B8",
      "-1.0"  = "#003F00", "0.0"  = "#009100", "1.0"  = "#73FF73",
      "-1.1"  = "#00007F", "0.1"  = "#1915ED", "1.1"  = "#96B3FF"
    ),
    guide = "none"
  ) +
  labs(x = "Air temperature (scaled)", y = "Predicted body temperature (scaled)")

```

## 1.2 Thermal tolerance across an elevational gradient

Next we want to know whether the upper thermal limit of butterflies changes across an elevational gradient, and therefore can adapt to local climatic conditions. We can used a mixed effects Cox model to test for changes in survival with temperature between groups while also controlling for non-independence in our data (species and sites again).

```{r}
#| label: ctmax ~ elevation full model

cox.1 <- coxme(Surv(Temp.end, A1) ~ Elevation + Sex + Condition + Nir_value + Length.centered + Species_WL +
                 Elevation:Sex + Elevation:Condition + Elevation:Nir_value + Elevation:Length.centered + Elevation:Species_WL +
                 (1 | Species) + (1 | Site), 
               data = ct.clean)
```

We cannot use the 'ranova' function in the 'lmerTest' package to test how our random effects are performing. Instead we must do this manually and compare AIC.

```{r}
#| label: ctmax ~ elevation random effects

# No random effects
cox.2 <- coxph(Surv(Temp.end, A1) ~ Elevation + Sex + Condition + Nir_value + Length.centered + Species_WL +
                 Elevation:Sex + Elevation:Condition + Elevation:Nir_value + Elevation:Length.centered + Elevation:Species_WL,
               data = ct.clean)

# Only species
cox.3 <- coxme(Surv(Temp.end, A1) ~ Elevation + Sex + Condition + Nir_value + Length.centered + Species_WL +
                 Elevation:Sex + Elevation:Condition + Elevation:Nir_value + Elevation:Length.centered + Elevation:Species_WL +
                 (1 | Species), 
               data = ct.clean)


# Only site
cox.4 <- coxme(Surv(Temp.end, A1) ~ Elevation + Sex + Condition + Nir_value + Length.centered + Species_WL +
                 Elevation:Sex + Elevation:Condition + Elevation:Nir_value + Elevation:Length.centered + Elevation:Species_WL +
                 (1 | Site), 
               data = ct.clean)

AIC(cox.1, cox.2, cox.3, cox.4)
```

Our first model with both species and site as random effects performs the best.

Next we must check the performance of our model.

```{r}
#| label: ctmax ~ elevation check model

cox.zph(cox.1)

plot(cox.zph(cox.1), var = "Condition")

plot(cox.zph(cox.1), var = "Sex")

```

Though some variables are significant, we can see that their effect is very weak (the lines are flat and horizontal). Therefore we can conclude that our model performs within acceptable ranges and we can test term significance.

```{r}
#| label: ctmax ~ elevation term significance

Anova(cox.1)
```

Now we can plot our results. First we need to recreate our model without random effects for plotting.

```{r}
#| label: ctmax ~ elevation plotting model

cox.plot <- coxph(Surv(Temp.end, A1) ~ Elevation + Sex + Condition + Nir_value + Length.centered + Species_WL +
                     Elevation:Sex + Elevation:Condition + Elevation:Nir_value + Elevation:Length.centered + Elevation:Species_WL,
                   data = ct.clean)
```

@fig-ctmax-elevation-sex shows how male butterflies have significantly lower upper thermal limits than females.

```{r}
#| label: fig-ctmax-elevation-sex
#| fig-cap: | 
#|  The difference in upper thermal limit between male and female butterflies.

newdat_sex <- data.frame(
  Elevation = mean(ct.clean$Elevation, na.rm = TRUE),
  Length.centered = mean(ct.clean$Length.centered, na.rm = TRUE),
  Nir_value = mean(ct.clean$Nir_value, na.rm = TRUE),
  Species_WL = mean(ct.clean$Species_WL, na.rm = TRUE),
  Sex = c("F", "M"),
  Condition = "2")

# Interpolate survival onto a regular grid to get smooth CIs
grid_time <- seq(40, 50, by = 0.1)

pred_list <- lapply(1:nrow(newdat_sex), function(i) {
  sf <- survfit(cox.plot, newdata = newdat_sex[i, ])
  
  # Linear interpolation of surv and CIs
  surv_interp <- approx(sf$time, sf$surv, xout = grid_time, method = "constant", yleft = 1, yright = 0)$y
  lower_interp <- approx(sf$time, sf$lower, xout = grid_time, method = "constant", yleft = 1, yright = 0)$y
  upper_interp <- approx(sf$time, sf$upper, xout = grid_time, method = "constant", yleft = 1, yright = 0)$y
  
  data.frame(
    time = grid_time,
    surv = surv_interp,
    surv_lower = lower_interp,
    surv_upper = upper_interp,
    Sex = newdat_sex$Sex[i]
  )
})

surv_df <- bind_rows(pred_list)

# For dashed lines need to extract LT50
lt50_df <- surv_df %>%
  group_by(Sex) %>%
  summarize(LT50 = min(time[surv <= 0.5]), .groups = "drop")

# Plot survival curves with confidence intervals and LT50 lines
ggplot(surv_df, aes(x = time, y = surv, color = Sex, fill = Sex)) +
  geom_line(linewidth = 1) +
  geom_ribbon(aes(ymin = surv_lower, ymax = surv_upper), alpha = 0.2, color = NA) +
  geom_segment(data = lt50_df,
               aes(x = 40, xend = LT50, y = 0.5, yend = 0.5),
               linetype = "dashed", color = "black") +
  geom_segment(data = lt50_df,
               aes(x = LT50, xend = LT50, y = 0, yend = 0.5, color = Sex),
               linetype = "dashed") +
  labs(x = "Temperature (°C)",
       y = "Survival probability",
       color = "Sex",
       fill = "Sex") +
  scale_color_manual(values = c("red", "blue")) +
  scale_fill_manual(values = c("red", "blue")) +
  xlim(40, 50) +
  ylim(0, 1)
```

Because we have a significant interaction between condition and elevation, we may want to plot separately for different elevations. @fig-ctmax-condition-elevation shows how the upper thermal limits of butterflies changes across elevations and conditions.

```{r}
#| label: fig-ctmax-condition-elevation
#| fig-cap: The relationship between wing condition and upper thermal limit across elevations.

elev_mean <- mean(ct.clean$Elevation, na.rm = TRUE)
elev_sd   <- sd(ct.clean$Elevation, na.rm = TRUE)
elev_levels <- c("-1" = elev_mean - elev_sd,
                 "0" = elev_mean,
                 "1" = elev_mean + elev_sd)

newdat <- expand.grid(
  Elevation = elev_levels,
  Condition = levels(ct.clean$Condition),
  Sex = "M",
  Length.centered = mean(ct.clean$Length.centered, na.rm = TRUE),
  Nir_value = mean(ct.clean$Nir_value, na.rm = TRUE),
  Species_WL = mean(ct.clean$Species_WL, na.rm = TRUE)
)
newdat$Elevation_Fac <- factor(newdat$Elevation,
                               levels = elev_levels,
                               labels = names(elev_levels))

# Interpolate survival onto a regular grid to get smooth CIs
grid_time <- seq(40, 50, by = 0.1)

pred_list <- lapply(1:nrow(newdat), function(i) {
  sf <- survfit(cox.plot, newdata = newdat[i, ])
  
  # Linear interpolation of surv and CIs
  surv_interp <- approx(sf$time, sf$surv, xout = grid_time, method = "constant", yleft = 1, yright = 0)$y
  lower_interp <- approx(sf$time, sf$lower, xout = grid_time, method = "constant", yleft = 1, yright = 0)$y
  upper_interp <- approx(sf$time, sf$upper, xout = grid_time, method = "constant", yleft = 1, yright = 0)$y
  
  data.frame(
    time = grid_time,
    surv = surv_interp,
    surv_lower = lower_interp,
    surv_upper = upper_interp,
    Condition = newdat$Condition[i],
    Elevation_Fac = newdat$Elevation_Fac[i]
  )
})

surv_df <- bind_rows(pred_list)

# For dashed lines need to extract LT50
lt50_df <- surv_df %>%
  group_by(Condition, Elevation_Fac) %>%
  summarize(LT50 = min(time[surv <= 0.5]), .groups = "drop")

# Plot survival curves with confidence intervals and LT50 lines
ggplot(surv_df, aes(x = time, y = surv, color = Condition, fill = Condition)) +
  geom_line(linewidth = 1) +
  geom_ribbon(aes(ymin = surv_lower, ymax = surv_upper), alpha = 0.2, color = NA) +
  facet_wrap(~Elevation_Fac, labeller = labeller(Elevation_Fac = ~ paste("Elevation:", .))) +
  geom_segment(data = lt50_df,
               aes(x = 40, xend = LT50, y = 0.5, yend = 0.5),
               linetype = "dashed", color = "black") +
  geom_segment(data = lt50_df,
               aes(x = LT50, xend = LT50, y = 0, yend = 0.5, color = Condition),
               linetype = "dashed") +
  labs(x = "Temperature (°C)",
       y = "Survival probability",
       color = "Condition",
       fill = "Condition") +
  scale_color_manual(values = c("blue", "purple", "red")) +
  scale_fill_manual(values = c("blue", "purple", "red")) +
  xlim(40, 50)
```

## 2. Functional traits across an elevational gradient

Now that we have tested whether thermal traits change across an elevational gradient, we may want to consider how functional traits may also respond. Here we an run a series of simple models to test for change. Because we have exactly the same sampling effort for every site (6 hours exactly), we can assume that every butterfly caught represents a random individual within that community.

First we can test whether butterfly species change in size across the elevational gradient.

```{r}
#| label: sp wing length ~ elevation model

m.1 <- lmer(Wing.length ~ Elevation + (1 | Site), data = buff.clean)

ranova(m.1)

check_model(m.1)

Anova(m.1)
```

Now we can plot our model. @fig-sp-wing-length-elevation-plot shows the raw values. @fig-sp-wing-length-elevation-mean-plot shows the same data but as mean values.

```{r}
#| label: fig-sp-wing-length-elevation-plot
#| fig-cap: The relationship between species wing length and elevation.

preds <- ggpredict(
  m.1, terms = c("Elevation"))

ggplot(preds, aes(x = x, y = predicted)) +
  geom_line(size = 1.2, color = "red") +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = group), 
              alpha = 0.2, colour = NA) +
  labs(
    x = "Elevation (m)",
    y = "Predicted wing length (cm)") +
    geom_jitter(data = buff.clean, aes(x = Elevation, y = Wing.length),
                alpha = 0.6, width = 10, height = 0.1) +
  theme(legend.position = "none")
```

```{r}
#| label: fig-sp-wing-length-elevation-mean-plot
#| fig-cap: The relationship between mean species wing length and elevation.

buff.means <- buff.clean %>%
  group_by(Site, Elevation) %>%
  summarise(Com_WL = mean(Wing.length, na.rm = TRUE),
            Com_Nir = mean(Nir_value, na.rm = TRUE),
            Com_SpWL = mean(Length.centered, na.rm = TRUE)) %>%
  ungroup()

preds <- ggpredict(
  m.1, terms = c("Elevation"))

ggplot(preds, aes(x = x, y = predicted)) +
  geom_line(size = 1.2, color = "red") +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = group), 
              alpha = 0.2, colour = NA) +
  labs(
    x = "Elevation (m)",
    y = "Predicted mean wing length (cm)") +
  geom_point(data = buff.means, aes(x = Elevation, y = Com_WL),
             size = 2) +
  theme(legend.position = "none")
```

Next we can test whether the individuals within species change in size across the elevational gradient.

```{r}
#| label: centred wing length ~ elevation model

ind.wl.1 <- lmer(Length.centered ~ Elevation + (1 | Species), data = buff.clean) 

ranova(ind.wl.1)

ind.wl.2 <- lm(Length.centered ~ Elevation, data = buff.clean)

AIC(ind.wl.1, ind.wl.2)

check_model(ind.wl.2)

Anova(ind.wl.2)
```

We can now plot our result. @fig-centred-wing-length-elevation shows this relationship on the raw data, @fig-centred-wing-length-elevation-mean shows the mean change.

```{r}
#| label: fig-centred-wing-length-elevation
#| fig-cap: The relationship between individual wing length (within species) and elevation.

preds <- ggpredict(
  ind.wl.2, terms = c("Elevation"))

ggplot(preds, aes(x = x, y = predicted)) +
  geom_line(size = 1.2, color = "red") +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = group), 
              alpha = 0.2, colour = NA) +
  labs(
    x = "Elevation (m)",
    y = "Predicted species-centred wing length (cm)") +
      geom_jitter(data = buff.clean, aes(x = Elevation, y = Length.centered),
                alpha = 0.6, width = 10, height = 0.01) +
  theme(legend.position = "none")
```

```{r}
#| label: fig-centred-wing-length-elevation-mean
#| fig-cap: The relationship between mean individual wing length (within species) and elevation.

buff.means <- buff.clean %>%
  group_by(Site, Elevation) %>%
  summarise(Com_WL = mean(Wing.length, na.rm = TRUE),
            Com_Nir = mean(Nir_value, na.rm = TRUE),
            Com_SpWL = mean(Length.centered, na.rm = TRUE)) %>%
  ungroup()

preds <- ggpredict(
  ind.wl.2, terms = c("Elevation"))

ggplot(preds, aes(x = x, y = predicted)) +
  geom_line(size = 1.2, color = "red") +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = group), 
              alpha = 0.2, colour = NA) +
  labs(
    x = "Elevation (m)",
    y = "Predicted mean species-centred wing length (cm)") +
  geom_point(data = buff.means, aes(x = Elevation, y = Com_SpWL),
             size = 2) +
  theme(legend.position = "none")
```

Finally we may want to know if species differ in colouration across the elevational gradient.

```{r}
#| label: nir ~ elevation model

m.2 <- lmer(Nir_value ~ Elevation + (1 | Site), data = buff.clean)

ranova(m.2)

m.2.2 <- lm(Nir_value ~ Elevation, data = buff.clean)

AIC(m.2, m.2.2)

check_model(m.2.2)

Anova(m.2.2)
```

Now we can plot our result. @fig-nir-elevation shows the relationship on the raw data, @fig-nir-elevation-mean shows the mean change.

```{r}
#| label: fig-nir-elevation
#| fig-cap: The relationship between colouration (NIR reflectance) and elevation.

preds <- ggpredict(
  m.2.2, terms = c("Elevation"))

ggplot(preds, aes(x = x, y = predicted)) +
  geom_line(size = 1.2, color = "red") +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = group), 
              alpha = 0.2, colour = NA) +
  labs(
    x = "Elevation (m)",
    y = "Predicted NIR reflectance (%)") +
  geom_jitter(data = buff.clean, aes(x = Elevation, y = Nir_value),
             size = 1.5, alpha = 0.5, width = 10, height = 1) +
  theme(legend.position = "none")
```

```{r}
#| label: fig-nir-elevation-mean
#| fig-cap: The relationship between mean colouration (NIR reflectance) and elevation.

buff.means <- buff.clean %>%
  group_by(Site, Elevation) %>%
  summarise(Com_WL = mean(Wing.length, na.rm = TRUE),
            Com_Nir = mean(Nir_value, na.rm = TRUE),
            Com_SpWL = mean(Length.centered, na.rm = TRUE)) %>%
  ungroup()

preds <- ggpredict(
  m.2.2, terms = c("Elevation"))

ggplot(preds, aes(x = x, y = predicted)) +
  geom_line(size = 1.2, color = "red") +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = group), 
              alpha = 0.2, colour = NA) +
  labs(
    x = "Elevation (m)",
    y = "Predicted NIR reflectance (%)") +
  geom_jitter(data = buff.means, aes(x = Elevation, y = Com_Nir),
             size = 1.5, alpha = 0.5, width = 10, height = 1) +
  theme(legend.position = "none")
```

## 3. Elevational range shifts and thermal traits

Finally, we may want to consider whether the thermal or functional traits we have recorded can explain recent responses to climate change. We can test our data against another dataset from [@kerner2023] that looked at the same species in the same mountain range.

First we need to load and merge the datasets.

```{r}
#| label: shift dataset

shift <- read.csv("Janika_shift.csv", header = T)

# Subset to the same species in the thermal traits dataset
shift.2 <- shift %>%
  filter(Species %in% buff.clean$Species)

buff.clean.3 <- left_join(buff.clean, shift.2, by = "Species")
buff.clean.3$scale.ElevShift <- scale(buff.clean.3$ElevShift)

# We need a new dataset with a slope value per species
sp <- unique(buff.clean.3$Species)

# Empty vectors to store results
slopes <- numeric(length(sp))

# Loop over species
for (i in seq_along(sp)) {
  dat.sp <- subset(buff.clean.3, Species == sp[i])
  fit <- lm(B.temp ~ A.temp, data = dat.sp)
  slopes[i] <- coef(fit)["A.temp"]
}

# Combine into dataframe
species_slopes <- data.frame(Species = sp, slope = slopes)

species_traits <- aggregate(cbind(ElevShift, Wing.length, Nir_value) ~ Species,
                            data = buff.clean.3, FUN = mean) 

# include ctmax? what's the overlap in species
load("CTmax_clean_final.RData")

species_ct <- aggregate(cbind(Temp.end) ~ Species,
                        data = ct.clean, FUN = mean) 
# only 8 species overlap, probably not enough

species_data <- merge(species_slopes, species_traits, by = "Species")
species_data <- left_join(species_data, species_ct, by = "Species")

# inverse slope so that high values = better and low values = worse
species_data$slope.inv <- 1-species_data$slope

# V. cardui slope looks very weird, exclude
species_data <- subset(species_data, Species != "Vanessa cardui")
```

Now we can test for a relationship with a linear model.

```{r}
#| label: shift ~ traits model

shift.lm.1 <- lm(ElevShift ~ slope.inv + Wing.length + Nir_value,
                 data = species_data)

# Does including CTmax improve model, even if we lose species?

dat_sub <- subset(species_data, !is.na(Temp.end)) # subset to smaller dataset to fairly test models

shift.lm.1.2 <- lm(ElevShift ~ slope.inv + Wing.length + Nir_value,
                   data = dat_sub)

shift.lm.2 <- lm(ElevShift ~ slope.inv + Wing.length + Nir_value + Temp.end,
                 data = dat_sub)

AIC(shift.lm.1.2, shift.lm.2)
# Including CTmax does not improve the model enough to justify keeping it

check_model(shift.lm.1)

r.squaredGLMM(shift.lm.1)

Anova(shift.lm.1)
```

We can now plot our results. @fig-slope-shift shows that species with a lower thermoregulatory capacity have shown greater upwards elevational range shifts. @fig-winglength-shift shows that larger species have shown greater upwards elevational range shifts.

```{r}
#| label: fig-slope-shift
#| fig-cap: The relationship between thermoregulatory capacity and elevational range shifts.

preds <- ggpredict(shift.lm.1, 
                   terms = c("slope.inv"))

ggplot(preds, aes(x = x, y = predicted)) +
  geom_line(size = 1.2, color = "red") +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = group), alpha = 0.2, colour = NA) +
  labs(x = "Thermoregulatory capacity",
       y = "Predicted elevational range shift (m)") +
  theme(legend.position = "none") +
  geom_abline(slope = 0, intercept = 0, linetype = "dashed", color = "black") +
  geom_point(data = species_data, aes(x = slope.inv, y = ElevShift), color = "black")
```

```{r}
#| label: fig-winglength-shift
#| fig-cap: The relationship between wing length and elevational range shifts.

preds <- ggpredict(shift.lm.1, 
                   terms = c("Wing.length"))

ggplot(preds, aes(x = x, y = predicted)) +
  geom_line(size = 1.2, color = "red") +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = group), alpha = 0.2, colour = NA) +
  labs(x = "Wing length (cm)",
       y = "Predicted elevational range shift (m)") +
  theme(legend.position = "none") +
  geom_abline(slope = 0, intercept = 0, linetype = "dashed", color = "black") +
  geom_point(data = species_data, aes(x = Wing.length, y = ElevShift), color = "black")
```