Exp. 1 Elevation survey

Data wrangling

Data uploading & wrangling - averaged values per the site and survey date. (all main effects in analyses are continuous variables)

Model - hopper num v elevation & date

model_hopper4 <- lm(sqrt(num.gh)~elev*survey.date + veg.height, data = hopper_elev_survey2)


par(mfrow=c(1,2))
shapiro.test(resid(model_hopper4))

## 
##  Shapiro-Wilk normality test
## 
## data:  resid(model_hopper4)
## W = 0.96913, p-value = 0.3229

hist(resid(model_hopper4))
qqnorm(resid(model_hopper4))
qqline(resid(model_hopper4))#non-normal, so transforming

Model - lack of fit testing - abundance

plot(fitted(model_hopper4), resid(model_hopper4))
abline(h = 0, lty = 2)

m_lin  <- lm(sqrt(num.gh)~elev*survey.date + veg.height, data = hopper_elev_survey2)
m_quad <- lm(sqrt(num.gh) ~ elev * survey.date +
                               I(elev^2) +
                               veg.height,
             data = hopper_elev_survey2)


anova(m_lin, m_quad)

## Analysis of Variance Table
## 
## Model 1: sqrt(num.gh) ~ elev * survey.date + veg.height
## Model 2: sqrt(num.gh) ~ elev * survey.date + I(elev^2) + veg.height
##   Res.Df    RSS Df Sum of Sq    F    Pr(>F)    
## 1     36 38.524                                
## 2     35 23.992  1    14.532 21.2 5.268e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

AIC(m_lin, m_quad)

##        df      AIC
## m_lin   6 125.7987
## m_quad  7 108.3825

par(mfrow=c(1,2))
shapiro.test(resid(m_quad))

## 
##  Shapiro-Wilk normality test
## 
## data:  resid(m_quad)
## W = 0.97554, p-value = 0.5123

hist(resid(m_quad))
qqnorm(resid(m_quad))
qqline(resid(m_quad))#non-normal, so transforming

summary(m_quad)

## 
## Call:
## lm(formula = sqrt(num.gh) ~ elev * survey.date + I(elev^2) + 
##     veg.height, data = hopper_elev_survey2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.26827 -0.65945 -0.02858  0.45162  1.92151 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       4.264e+02  7.563e+01   5.638 2.31e-06 ***
## elev             -2.260e-01  4.167e-02  -5.423 4.45e-06 ***
## survey.date      -6.917e-01  2.425e-01  -2.852  0.00724 ** 
## I(elev^2)         2.908e-05  6.316e-06   4.604 5.27e-05 ***
## veg.height       -3.769e-02  2.719e-02  -1.386  0.17446    
## elev:survey.date  2.168e-04  7.829e-05   2.769  0.00893 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8279 on 35 degrees of freedom
## Multiple R-squared:  0.5744, Adjusted R-squared:  0.5136 
## F-statistic: 9.448 on 5 and 35 DF,  p-value: 9.072e-06

a1<-Anova(m_quad, type=3)
anova_df1 <- as.data.frame(a1)
kable(anova_df1, digits = 4, caption = "ANOVA Table (type 3)")# Print nicely formatted table

ANOVA Table (type 3)
	Sum Sq	Df	F value	Pr(>F)
(Intercept)	21.7922	1	31.7914	0.0000
elev	20.1585	1	29.4080	0.0000
survey.date	5.5763	1	8.1349	0.0072
I(elev^2)	14.5320	1	21.1999	0.0001
veg.height	1.3171	1	1.9215	0.1745
elev:survey.date	5.2570	1	7.6691	0.0089
Residuals	23.9917	35	NA	NA

avg_data <- hopper_elev_survey2 %>%
  group_by(elev) %>%
  summarise(mean_num_gh = mean(num.gh, na.rm = TRUE))

newdata <- data.frame(
  elev = seq(min(hopper_elev_survey2$elev),
             max(hopper_elev_survey2$elev),
             length.out = 200),
  survey.date = mean(hopper_elev_survey2$survey.date),  # fix at mean
  veg.height = mean(hopper_elev_survey2$veg.height)     # fix at mean
)

# Predict and back-transform
newdata$pred <- predict(m_quad, newdata = newdata)
newdata$pred_orig <- newdata$pred^2

ggplot() +
  geom_point(data = avg_data, aes(x = elev, y = mean_num_gh),
             color = "black", size = 2) +       # averaged points
  geom_line(data = newdata, aes(x = elev, y = pred_orig),
            color = "blue", size = 1.5) +      # smooth curve
  labs(x = "Elevation", y = "Number of GH",
       title = "Quadratic Regression with Mean Observed Points") +
  theme_minimal()

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

Model - hopper biomass v elevation & date

model_hopper6 <- lm(sqrt(length.tot)~elev*survey.date +veg.height, data = hopper_elev_survey2)

par(mfrow=c(1,2))
shapiro.test(resid(model_hopper6))

## 
##  Shapiro-Wilk normality test
## 
## data:  resid(model_hopper6)
## W = 0.97637, p-value = 0.5413

hist(resid(model_hopper6))
qqnorm(resid(model_hopper6))
qqline(resid(model_hopper6))#non-normal, so transforming

Model - lack of fit testing - biomass

plot(fitted(model_hopper6), resid(model_hopper6))
abline(h = 0, lty = 2)

m_lin  <- lm(sqrt(length.tot)~elev*survey.date + veg.height, data = hopper_elev_survey2)
m_quad <- lm(sqrt(length.tot) ~ elev * survey.date +
                               I(elev^2) +
                               veg.height,
             data = hopper_elev_survey2)


anova(m_lin, m_quad)

## Analysis of Variance Table
## 
## Model 1: sqrt(length.tot) ~ elev * survey.date + veg.height
## Model 2: sqrt(length.tot) ~ elev * survey.date + I(elev^2) + veg.height
##   Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
## 1     36 374.52                                  
## 2     35 219.42  1     155.1 24.739 1.736e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

AIC(m_lin, m_quad)

##        df      AIC
## m_lin   6 219.0480
## m_quad  7 199.1276

par(mfrow=c(1,2))
shapiro.test(resid(m_quad))

## 
##  Shapiro-Wilk normality test
## 
## data:  resid(m_quad)
## W = 0.97196, p-value = 0.3985

hist(resid(m_quad))
qqnorm(resid(m_quad))
qqline(resid(m_quad))#non-normal, so transforming

summary(m_quad)

## 
## Call:
## lm(formula = sqrt(length.tot) ~ elev * survey.date + I(elev^2) + 
##     veg.height, data = hopper_elev_survey2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.9512 -2.1276  0.2696  1.5811  6.1175 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       1.306e+03  2.287e+02   5.710 1.86e-06 ***
## elev             -7.117e-01  1.260e-01  -5.647 2.26e-06 ***
## survey.date      -1.837e+00  7.335e-01  -2.505    0.017 *  
## I(elev^2)         9.501e-05  1.910e-05   4.974 1.74e-05 ***
## veg.height       -6.802e-02  8.223e-02  -0.827    0.414    
## elev:survey.date  5.773e-04  2.368e-04   2.438    0.020 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.504 on 35 degrees of freedom
## Multiple R-squared:  0.5844, Adjusted R-squared:  0.5251 
## F-statistic: 9.844 on 5 and 35 DF,  p-value: 6.124e-06

a1<-Anova(m_quad, type=3)
anova_df1 <- as.data.frame(a1)
kable(anova_df1, digits = 4, caption = "ANOVA Table (type 3)")# Print nicely formatted table

ANOVA Table (type 3)
	Sum Sq	Df	F value	Pr(>F)
(Intercept)	204.4198	1	32.6068	0.0000
elev	199.8928	1	31.8847	0.0000
survey.date	39.3421	1	6.2754	0.0170
I(elev^2)	155.0971	1	24.7394	0.0000
veg.height	4.2892	1	0.6842	0.4138
elev:survey.date	37.2782	1	5.9462	0.0200
Residuals	219.4231	35	NA	NA

avg_data <- hopper_elev_survey2 %>%
  group_by(elev) %>%
  summarise(mean_len_gh = mean(length.tot, na.rm = TRUE))

newdata <- data.frame(
  elev = seq(min(hopper_elev_survey2$elev),
             max(hopper_elev_survey2$elev),
             length.out = 200),
  survey.date = mean(hopper_elev_survey2$survey.date),  # fix at mean
  veg.height = mean(hopper_elev_survey2$veg.height)     # fix at mean
)

# Predict and back-transform
newdata$pred <- predict(m_quad, newdata = newdata)
newdata$pred_orig <- newdata$pred^2

ggplot() +
  geom_point(data = avg_data, aes(x = elev, y = mean_len_gh),
             color = "black", size = 2) +       # averaged points
  geom_line(data = newdata, aes(x = elev, y = pred_orig),
            color = "blue", size = 1.5) +      # smooth curve
  labs(x = "Elevation", y = "biomass of GH",
       title = "Quadratic Regression with Mean Observed Points") +
  theme_minimal()

Table 1. Analysis of Grasshopper abundance and biomass by Elevation

Abundance

ANOVA Table (type 3)
	Sum Sq	Df	F value	Pr(>F)
(Intercept)	7.3147	1	6.8355	0.0130
elev	6.6779	1	6.2405	0.0172
survey.date	6.3703	1	5.9530	0.0197
veg.height	0.0004	1	0.0004	0.9845
elev:survey.date	5.9284	1	5.5401	0.0242
Residuals	38.5237	36	NA	NA

Biomass

ANOVA Table (type 3)
	Sum Sq	Df	F value	Pr(>F)
(Intercept)	53.2169	1	5.1154	0.0299
elev	48.7537	1	4.6863	0.0371
survey.date	46.2678	1	4.4474	0.0420
veg.height	2.8597	1	0.2749	0.6033
elev:survey.date	43.1468	1	4.1474	0.0491
Residuals	374.5202	36	NA	NA

Fig 1a. Hopper vs Elevation

## `geom_smooth()` using formula = 'y ~ x'

Partial r =

## [1] 0.2355848

Supplementary (exp 1)

Model - hopper length v elevation & date

model_hopper5 <- lm(length.ave~elev*survey.date +veg.height , data = hopper_elev_survey2)


par(mfrow=c(1,2))
shapiro.test(resid(model_hopper5))

## 
##  Shapiro-Wilk normality test
## 
## data:  resid(model_hopper5)
## W = 0.96896, p-value = 0.4153

hist(resid(model_hopper5))
qqnorm(resid(model_hopper5))
qqline(resid(model_hopper5))#non-normal, so transforming

Table S2. Analysis of Grasshopper length by Elevation

Length

ANOVA Table (type 3)
	Sum Sq	Df	F value	Pr(>F)
(Intercept)	6.6703	1	1.3627	0.2523
elev	6.6803	1	1.3647	0.2519
survey.date	7.1308	1	1.4567	0.2369
veg.height	18.6441	1	3.8088	0.0604
elev:survey.date	6.5123	1	1.3304	0.2578
Residuals	146.8520	30	NA	NA

Fig S1. Hopper biomass survey per date

## `geom_smooth()` using formula = 'y ~ x'

## Warning: Some predictor variables are on very different scales: consider
## rescaling
## Warning: Some predictor variables are on very different scales: consider
## rescaling

Partial r =

## [1] 0.1116878

Fig S2. Hopper biomass (across survey preiod)

## `geom_smooth()` using formula = 'y ~ x'

Exp. 2 Elevation Feeding Trials

Data wrangling

Data uploading & wrangling - models subset data to hour 17 of the bioassay - all main and random effects are factors (plant species, GH elev, plant elev, and the random effect of arena), percent herbivory is numeric

Model - Percent herbivory no interactions

# Subset the data
dat_sub <- subset(feed2.dat, hours == "17")

# Create bounded proportion variable (if not already created)
dat_sub$prop_herb <- (dat_sub$percent_herb + 0.01) / 100.02  # adjust if needed

# Fit the model
model <- glmmTMB(
  prop_herb ~ plant_species + plant_elev + gh_elev + (1 | gh_id),
  data = dat_sub,
  family = beta_family()
)


plot(residuals(model, type = "pearson") ~ fitted(model))

qqnorm(residuals(model, type = "pearson"))
qqline(residuals(model, type = "pearson"))

modelx <- glmmTMB(
  prop_herb ~ plant_species + plant_elev + gh_elev + (1 | gh_id),
  data = dat_sub,
  family = beta_family()
)

at<-Anova(modelx, type = 2)
anova_df <- as.data.frame(at)
kable(anova_df, digits = 3, caption = "ANOVA Table (type 2)")

ANOVA Table (type 2)
	Chisq	Df	Pr(>Chisq)
plant_species	3.673	4	0.452
plant_elev	3.234	1	0.072
gh_elev	4.591	1	0.032

modely <- glmmTMB(
  prop_herb ~ plant_species*plant_elev * gh_elev + (1 | gh_id),
  data = dat_sub,
  family = beta_family()
)

at<-Anova(modely, type = 2)
anova_df <- as.data.frame(at)
kable(anova_df, digits = 3, caption = "ANOVA Table (type 2)")

ANOVA Table (type 2)
	Chisq	Df	Pr(>Chisq)
plant_species	3.976	4	0.409
plant_elev	3.587	1	0.058
gh_elev	5.109	1	0.024
plant_species:plant_elev	1.279	4	0.865
plant_species:gh_elev	3.356	4	0.500
plant_elev:gh_elev	1.464	1	0.226
plant_species:plant_elev:gh_elev	1.791	4	0.774

modelz <- glmmTMB(
  prop_herb ~ plant_species*plant_elev+ plant_species*gh_elev+ plant_elev*gh_elev + (1 | gh_id),
  data = dat_sub,
  family = beta_family()
)

at<-Anova(modelz, type = 2)
anova_df <- as.data.frame(at)
kable(anova_df, digits = 3, caption = "ANOVA Table (type 2)")

ANOVA Table (type 2)
	Chisq	Df	Pr(>Chisq)
plant_species	3.957	4	0.412
plant_elev	3.578	1	0.059
gh_elev	5.062	1	0.024
plant_species:plant_elev	1.275	4	0.866
plant_species:gh_elev	3.345	4	0.502
plant_elev:gh_elev	1.453	1	0.228

modelx <- lmer(log(percent_herb+1) ~ plant_species+plant_elev + gh_elev + (1|gh_id), data = feed2.dat, subset = hours == "17")

modely <- lmer(log(percent_herb+1) ~ plant_species*plant_elev * gh_elev + (1|gh_id), data = feed2.dat, subset = hours == "17")

modelz <- lmer(log(percent_herb+1) ~ plant_species*plant_elev+ plant_species*gh_elev+ plant_elev*gh_elev + (1|gh_id), data = feed2.dat, subset = hours == "17")


anova(modelx,modely,modelz)

## refitting model(s) with ML (instead of REML)

## Data: feed2.dat
## Subset: hours == "17"
## Models:
## modelx: log(percent_herb + 1) ~ plant_species + plant_elev + gh_elev + (1 | gh_id)
## modelz: log(percent_herb + 1) ~ plant_species * plant_elev + plant_species * gh_elev + plant_elev * gh_elev + (1 | gh_id)
## modely: log(percent_herb + 1) ~ plant_species * plant_elev * gh_elev + (1 | gh_id)
##        npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## modelx    9 359.11 382.56 -170.56   341.11                     
## modelz   18 370.28 417.17 -167.14   334.28 6.8322  9     0.6546
## modely   22 374.47 431.78 -165.23   330.47 3.8120  4     0.4320

AIC(modelx,modely,modelz)

##        df      AIC
## modelx  9 361.8095
## modely 22 358.4124
## modelz 18 363.0879

model <- lmer(log(percent_herb+1) ~ plant_species+plant_elev + gh_elev + (1|gh_id), data = feed2.dat, subset = hours == "17")

par(mfrow=c(1,2))
shapiro.test(resid(model)) # p < 0.05 means non-normal residuals

## 
##  Shapiro-Wilk normality test
## 
## data:  resid(model)
## W = 0.98259, p-value = 0.2105

hist(resid(model)) # visualize the residuals
qqnorm(resid(model)) # if normal residuals, should fall along the 1:1 line
qqline(resid(model)) # add the 1:1 line

Table 2. Analysis of grasshopper herbivory

ANOVA Table (type 2)
	Chisq	Df	Pr(>Chisq)
plant_species	9.681	4	0.046
plant_elev	10.954	1	0.001
gh_elev	8.187	1	0.004

Fig 1b. Plant elevation main effect

#this is raw in the first v

#this should be in the main fig?

lsmeans1_plant_elev <- emmeans(model, ~ plant_elev, type = "response")
lsmeans1_plant_elev_df <- as.data.frame(lsmeans1_plant_elev)

ggplot(lsmeans1_plant_elev_df, aes(x = plant_elev, y = response)) +
  geom_col(fill = "#BFD8B8",position = position_dodge(), color = "black") +
  geom_errorbar(aes(ymin = response - SE, ymax = response + SE), 
                width = 0.2, position = position_dodge(0.9)) +  
  labs( x = "Source Plant Elevation (m)", y = "Herbivory (%)")+
  theme_minimal(base_size = 16) +  
  theme(
    axis.text.x = element_text(hjust = 0.5, vjust = 0.5, size = 14),  
    axis.text.y = element_text(size = 14),  
    axis.title.x = element_text(size = 16, face = "bold"),  
    axis.title.y = element_text(size = 16, face = "bold"),  
    panel.grid.major = element_blank(),  
    panel.grid.minor = element_blank(),  
    panel.border = element_blank(),  
    axis.line = element_line(color = "black", linewidth = 1),  # Add black x and y axis lines
    legend.position = "none"  # Remove legend
  )+
  scale_x_discrete(labels = c("L" = "Low (2872) ",
                              "H" = "High (3360) "), 
                   limits = c("L","H"))

Supplementary

Model - Percent herbivory with three way interaction

model <- lmer(log(percent_herb+1) ~ plant_species*plant_elev*gh_elev +(1|gh_id), data = feed2.dat, subset = hours == "17")


par(mfrow=c(1,2))
shapiro.test(resid(model)) # p < 0.05 means non-normal residuals

## 
##  Shapiro-Wilk normality test
## 
## data:  resid(model)
## W = 0.98856, p-value = 0.5504

hist(resid(model)) # visualize the residuals
qqnorm(resid(model)) # if normal residuals, should fall along the 1:1 line
qqline(resid(model)) # add the 1:1 line

Model - Percent herbivory with two way interaction

model1 <- lmer(log(percent_herb+1) ~ plant_species*plant_elev + plant_elev*gh_elev +  gh_elev*plant_species+ (1|gh_id), data = feed2.dat, subset = hours == "17")

par(mfrow=c(1,2))
shapiro.test(resid(model)) # p < 0.05 means non-normal residuals

## 
##  Shapiro-Wilk normality test
## 
## data:  resid(model)
## W = 0.98856, p-value = 0.5504

hist(resid(model)) # visualize the residuals
qqnorm(resid(model)) # if normal residuals, should fall along the 1:1 line
qqline(resid(model)) # add the 1:1 line

Table S3. Percent herbivory with three way interaction

ANOVA Table (type 2)
	Chisq	Df	Pr(>Chisq)
plant_species	9.427	4	0.051
plant_elev	10.329	1	0.001
gh_elev	7.972	1	0.005
plant_species:plant_elev	1.315	4	0.859
plant_species:gh_elev	2.846	4	0.584
plant_elev:gh_elev	1.718	1	0.190
plant_species:plant_elev:gh_elev	3.169	4	0.530

Table S4. Percent herbivory with two way interaction

ANOVA Table (type 2)
	Chisq	Df	Pr(>Chisq)
plant_species	9.427	4	0.051
plant_elev	10.528	1	0.001
gh_elev	7.972	1	0.005
plant_species:plant_elev	1.341	4	0.854
plant_elev:gh_elev	1.751	1	0.186
plant_species:gh_elev	2.846	4	0.584

Fig S3. Grasshopper elevation main effect

change to this : raw GH elevation

feed2.22.plant_elev.summary <-feed2.dat %>%  
  filter(hours=="17") 

ggplot(feed2.22.plant_elev.summary, aes(x = gh_elev, y = percent_herb)) +
    labs( x = "Source Grasshopper Elevation (m)", y = "Herbivory (%)")+
  geom_violin(fill = "#BFD8B8", color = "black", alpha = 0.5, trim = FALSE) +  # violin shape
  geom_boxplot(width = 0.1, fill = "white", color = "black", outlier.shape = NA) +  # boxplot inside violin
  stat_summary(fun = mean, geom = "point", shape = 18, size = 3, color = "red") + # mean point
  theme_minimal(base_size = 16) +  
  theme(
    axis.text.x = element_text(hjust = 0.5, vjust = 0.5, size = 14),  
    axis.text.y = element_text(size = 14),  
    axis.title.x = element_text(size = 16, face = "bold"),  
    axis.title.y = element_text(size = 16, face = "bold"),  
    panel.grid.major = element_blank(),  
    panel.grid.minor = element_blank(),  
    panel.border = element_blank(),  
    axis.line = element_line(color = "black", linewidth = 1),  
    legend.position = "none"  
  ) +
  scale_x_discrete(labels = c("L" = "Low (2872) ",
                              "H" = "High (3360) "), 
                   limits = c("L","H"))+
  coord_cartesian(ylim = c(-15, 110))   # bound y-axis from 0 to 100

Fig S4. Species main effect

#### change to this :

feed2.22.plant_elev.summary <-feed2.dat %>%  
  filter(hours=="17") 

ggplot(feed2.22.plant_elev.summary, aes(x = plant_species, y = percent_herb)) +
  labs( x = "", y = "Herbivory (%)")+
  geom_violin(fill = "#966e83", color = "black", alpha = 0.5, trim = FALSE) +  # violin shape
  geom_boxplot(width = 0.1, fill = "white", color = "black", outlier.shape = NA) +  # boxplot inside violin
  stat_summary(fun = mean, geom = "point", shape = 18, size = 3, color = "red") + # mean point
  theme_minimal(base_size = 16) +  
  theme(
    axis.text.x = element_text(angle = 45, hjust = 0.5, vjust = 0.5, size = 14),  
    axis.text.y = element_text(size = 14),  
    axis.title.x = element_text(size = 16, face = "bold"),  
    axis.title.y = element_text(size = 16, face = "bold"),  
    panel.grid.major = element_blank(),  
    panel.grid.minor = element_blank(),  
    panel.border = element_blank(),  
    axis.line = element_line(color = "black", linewidth = 1),  
    legend.position = "none"  
  ) +
  scale_x_discrete(labels = c(
    "D" = "D. nuttallianum",
    "E" = "E. speciosus",
    "F" = "F. virginiana",
    "H" = "H. quinquenervis",
    "P" = "P. pulcherrima"))

Fig S5. Interaction terms figure

#### change to this :

feed2.22.plant_elev.summary <-feed2.dat %>%  
  filter(hours=="17") 

pd <- position_dodge(width = 0.7)

ggplot(
  feed2.22.plant_elev.summary,
  aes(
    x = plant_species,
    y = percent_herb,
    fill = plant_elev,
    group = interaction(plant_species, plant_elev)
  )
) +
  facet_wrap(
    ~ gh_elev,
    ncol = 1,
    labeller = labeller(
      gh_elev = c(
        "H" = "Source Grasshopper Elevation – High (3360 m)",
        "L" = "Source Grasshopper Elevation – Low (2872 m)"
      )
    )
  ) +
geom_violin(
  position = pd,
  width = 1.1,
  alpha = 0.5,
  color = "black",
  trim = FALSE
) +
geom_boxplot(
  position = pd,
  width = 0.18,
  fill = "white",
  color = "black",
  outlier.shape = NA
) +
stat_summary(
  fun = mean,
  geom = "point",
  shape = 18,
  size = 3,
  color = "red",
  position = pd,
  show.legend = FALSE
) +
  labs(x = "", y = "Herbivory (%)") +
  scale_fill_manual(
    name = "Source Plant Elevation",
    values = c("H" = "#BFD8B8", "L" = "#4F7F5A")
  ) +
  scale_x_discrete(labels = c(
    "D" = "D. nuttallianum",
    "E" = "E. speciosus",
    "F" = "F. virginiana",
    "H" = "H. quinquenervis",
    "P" = "P. pulcherrima"
  )) +
    theme(
    axis.text.x = element_text(angle = 45, hjust = 1, size = 14),  
    strip.text = element_text(size = 14),
    axis.text.y = element_text(size = 14),  
    axis.title.x = element_text(size = 16, face = "bold"),  
    axis.title.y = element_text(size = 16, face = "bold"),  
    panel.grid.major = element_blank(),  
    panel.grid.minor = element_blank(),  
    panel.border = element_blank(),  
    axis.line = element_line(color = "black"),
    strip.background = element_rect(fill = "grey", color = NA)  
    )+
  coord_cartesian(ylim = c(-15, 110))

## Warning: `position_dodge()` requires non-overlapping x intervals.
## `position_dodge()` requires non-overlapping x intervals.

### Fig XX. for supp - raw elevation v herbivory

Exp. 3 Climate Manipulation Feeding Trials

Data wrangling

Data uploading & wrangling

this code allows us to analyze by family. we are not persuing this right now.

Model - Mass with main effects (no interactions)

GNA.model1 <- lmer(sqrt(Mass) ~ Snow+Temp+Species+ (1|Snow:Plot) +(1|Plant.ID), 
                   data = GNA.weights)

par(mfrow=c(1,2))
shapiro.test(resid(GNA.model1))

## 
##  Shapiro-Wilk normality test
## 
## data:  resid(GNA.model1)
## W = 0.98671, p-value = 2.541e-05

hist(resid(GNA.model1))
qqnorm(resid(GNA.model1))
qqline(resid(GNA.model1))#non-normal, so transforming

what about these models? (sqrt and glmm gamma log link)

GNA.model_log <- lmer(
  log(Mass + 1) ~ Snow + Temp + Species +
    (1 | Snow:Plot) + (1 | Plant.ID),
  data = GNA.weights
)


par(mfrow=c(1,2))
shapiro.test(resid(GNA.model_log))

## 
##  Shapiro-Wilk normality test
## 
## data:  resid(GNA.model_log)
## W = 0.9897, p-value = 0.0002985

hist(resid(GNA.model_log))
qqnorm(resid(GNA.model_log))
qqline(resid(GNA.model_log))#non-normal, so transforming

at_log <- Anova(GNA.model_log, test.statistic = "F", type = 2)
anova_log <- as.data.frame(at_log)

kable(anova_log, digits = 4,
      caption = "Type II ANOVA – log(Mass + 1)")

Type II ANOVA – log(Mass + 1)
	F	Df	Df.res	Pr(>F)
Snow	5.4206	1	3.9987	0.0804
Temp	7.5818	1	241.0003	0.0063
Species	12.6585	13	254.0184	0.0000

GNA.model_gamma <- glmmTMB(
  Mass ~ Snow + Temp + Species +
    (1 | Snow:Plot) + (1 | Plant.ID),
  data = GNA.weights,
  family = Gamma(link = "log")
)

res <- simulateResiduals(GNA.model_gamma)
plot(res)#this cheks the scaled residuals with the fitted, looking for a strong deviation

overdispersion_test <- sum(residuals(GNA.model_gamma, type = "pearson")^2) / df.residual(GNA.model_gamma) #this checkes for overdispersion (meaning that there's a high amount of variance for a binomial mode) Values >> 1 suggest overdispersion
overdispersion_test

## [1] 0.2007613

at_gamma <- Anova(GNA.model_gamma, type = 2)
anova_gamma <- as.data.frame(at_gamma)

kable(anova_gamma, digits = 4,
      caption = "Type II ANOVA – Gamma GLMM (log link)")

Type II ANOVA – Gamma GLMM (log link)
	Chisq	Df	Pr(>Chisq)
Snow	9.2433	1	0.0024
Temp	9.5811	1	0.0020
Species	203.1293	13	0.0000

Model - Survival (binary) with 3 way interaction

GNA.model2 <- glmer(Survival_binary ~  Snow * Temp * Species +(1|Plant.ID) + (1|Plot) , 
                    data = GNA.survive.binary, 
                    family = binomial(link = "logit"))

## Warning in (function (fn, par, lower = rep.int(-Inf, n), upper = rep.int(Inf, :
## failure to converge in 10000 evaluations

## Warning in optwrap(optimizer, devfun, start, rho$lower, control = control, :
## convergence code 4 from Nelder_Mead: failure to converge in 10000 evaluations

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## unable to evaluate scaled gradient

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge: degenerate Hessian with 1 negative eigenvalues

summary(GNA.model2)

## Warning in vcov.merMod(object, use.hessian = use.hessian): variance-covariance matrix computed from finite-difference Hessian is
## not positive definite or contains NA values: falling back to var-cov estimated from RX

## Warning in vcov.merMod(object, correlation = correlation, sigm = sig): variance-covariance matrix computed from finite-difference Hessian is
## not positive definite or contains NA values: falling back to var-cov estimated from RX

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Survival_binary ~ Snow * Temp * Species + (1 | Plant.ID) + (1 |  
##     Plot)
##    Data: GNA.survive.binary
## 
##      AIC      BIC   logLik deviance df.resid 
##   2055.2   2366.6   -969.6   1939.2     1527 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.5668 -0.8129 -0.5224  1.0338  5.3041 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev. 
##  Plant.ID (Intercept) 8.052e-02 0.2837629
##  Plot     (Intercept) 1.654e-07 0.0004067
## Number of obs: 1585, groups:  Plant.ID, 317; Plot, 6
## 
## Fixed effects:
##                                              Estimate Std. Error z value
## (Intercept)                                  -1.02629    0.43014  -2.386
## SnowSnow_M                                    0.31683    0.59109   0.536
## TempTemp_W                                    0.61116    0.58111   1.052
## SpeciesDephinium                             -0.87347    0.70185  -1.245
## SpeciesDugaldia                               1.88839    0.59847   3.155
## SpeciesErigeron                              -0.38327    0.72236  -0.531
## SpeciesFrageria                               0.89059    0.57657   1.545
## SpeciesGaliium                                0.82122    0.63841   1.286
## SpeciesGrass                                  0.16086    0.59864   0.269
## SpeciesLathyrus                              -0.18592    0.62206  -0.299
## SpeciesPotentilla                             1.02661    0.57614   1.782
## SpeciesPseudocymopsos                         0.16360    0.59854   0.273
## SpeciesSenecio                                0.88983    0.57667   1.543
## SpeciesTeraxicum                              0.75419    0.57827   1.304
## SpeciesValerianaF                             1.58325    0.58526   2.705
## SpeciesValerianaM                             0.61303    0.58117   1.055
## SnowSnow_M:TempTemp_W                        -0.45796    0.81213  -0.564
## SnowSnow_M:SpeciesDephinium                  -1.09040    1.10529  -0.987
## SnowSnow_M:SpeciesDugaldia                   -1.45192    0.81974  -1.771
## SnowSnow_M:SpeciesErigeron                    0.81913    0.91415   0.896
## SnowSnow_M:SpeciesFrageria                   -1.21058    0.82585  -1.466
## SnowSnow_M:SpeciesGaliium                    -0.31623    0.89126  -0.355
## SnowSnow_M:SpeciesGrass                      -0.47946    0.84134  -0.570
## SnowSnow_M:SpeciesLathyrus                    1.30921    0.83909   1.560
## SnowSnow_M:SpeciesPotentilla                 -0.31717    0.80209  -0.395
## SnowSnow_M:SpeciesPseudocymopsos              0.13135    0.82182   0.160
## SnowSnow_M:SpeciesSenecio                    -0.94774    0.83538  -1.135
## SnowSnow_M:SpeciesTeraxicum                  -0.31915    0.80527  -0.396
## SnowSnow_M:SpeciesValerianaF                 -2.28401    0.85528  -2.670
## SnowSnow_M:SpeciesValerianaM                 -0.03999    0.80623  -0.050
## TempTemp_W:SpeciesDephinium                  -2.11504    1.31456  -1.609
## TempTemp_W:SpeciesDugaldia                   -0.91483    0.81762  -1.119
## TempTemp_W:SpeciesErigeron                   -0.60796    1.42102  -0.428
## TempTemp_W:SpeciesFrageria                   -0.61143    0.79541  -0.769
## TempTemp_W:SpeciesGaliium                    -1.27008    0.90564  -1.402
## TempTemp_W:SpeciesGrass                      -1.15521    0.85765  -1.347
## TempTemp_W:SpeciesLathyrus                    0.04469    0.83497   0.054
## TempTemp_W:SpeciesPotentilla                 -1.02541    0.79835  -1.284
## TempTemp_W:SpeciesPseudocymopsos             -0.77857    0.83438  -0.933
## TempTemp_W:SpeciesSenecio                    -0.33821    0.79553  -0.425
## TempTemp_W:SpeciesTeraxicum                  -0.75412    0.80005  -0.943
## TempTemp_W:SpeciesValerianaF                 -0.75362    0.80507  -0.936
## TempTemp_W:SpeciesValerianaM                  0.21577    0.80201   0.269
## SnowSnow_M:TempTemp_W:SpeciesDephinium      -10.88773  428.74628  -0.025
## SnowSnow_M:TempTemp_W:SpeciesDugaldia         0.62095    1.13693   0.546
## SnowSnow_M:TempTemp_W:SpeciesErigeron         0.14136    1.63742   0.086
## SnowSnow_M:TempTemp_W:SpeciesFrageria         0.45766    1.15113   0.398
## SnowSnow_M:TempTemp_W:SpeciesGaliium          1.11719    1.25975   0.887
## SnowSnow_M:TempTemp_W:SpeciesGrass            1.32152    1.18622   1.114
## SnowSnow_M:TempTemp_W:SpeciesLathyrus        -1.02560    1.15085  -0.891
## SnowSnow_M:TempTemp_W:SpeciesPotentilla       0.87236    1.11958   0.779
## SnowSnow_M:TempTemp_W:SpeciesPseudocymopsos   0.01063    1.16449   0.009
## SnowSnow_M:TempTemp_W:SpeciesSenecio          0.95186    1.17297   0.811
## SnowSnow_M:TempTemp_W:SpeciesTeraxicum        0.73924    1.12226   0.659
## SnowSnow_M:TempTemp_W:SpeciesValerianaF       1.87424    1.15858   1.618
## SnowSnow_M:TempTemp_W:SpeciesValerianaM      -0.50623    1.12353  -0.451
##                                             Pr(>|z|)   
## (Intercept)                                  0.01704 * 
## SnowSnow_M                                   0.59195   
## TempTemp_W                                   0.29294   
## SpeciesDephinium                             0.21330   
## SpeciesDugaldia                              0.00160 **
## SpeciesErigeron                              0.59571   
## SpeciesFrageria                              0.12243   
## SpeciesGaliium                               0.19832   
## SpeciesGrass                                 0.78815   
## SpeciesLathyrus                              0.76503   
## SpeciesPotentilla                            0.07477 . 
## SpeciesPseudocymopsos                        0.78459   
## SpeciesSenecio                               0.12282   
## SpeciesTeraxicum                             0.19216   
## SpeciesValerianaF                            0.00683 **
## SpeciesValerianaM                            0.29151   
## SnowSnow_M:TempTemp_W                        0.57282   
## SnowSnow_M:SpeciesDephinium                  0.32387   
## SnowSnow_M:SpeciesDugaldia                   0.07653 . 
## SnowSnow_M:SpeciesErigeron                   0.37022   
## SnowSnow_M:SpeciesFrageria                   0.14269   
## SnowSnow_M:SpeciesGaliium                    0.72273   
## SnowSnow_M:SpeciesGrass                      0.56877   
## SnowSnow_M:SpeciesLathyrus                   0.11869   
## SnowSnow_M:SpeciesPotentilla                 0.69252   
## SnowSnow_M:SpeciesPseudocymopsos             0.87301   
## SnowSnow_M:SpeciesSenecio                    0.25658   
## SnowSnow_M:SpeciesTeraxicum                  0.69186   
## SnowSnow_M:SpeciesValerianaF                 0.00757 **
## SnowSnow_M:SpeciesValerianaM                 0.96044   
## TempTemp_W:SpeciesDephinium                  0.10763   
## TempTemp_W:SpeciesDugaldia                   0.26319   
## TempTemp_W:SpeciesErigeron                   0.66877   
## TempTemp_W:SpeciesFrageria                   0.44208   
## TempTemp_W:SpeciesGaliium                    0.16079   
## TempTemp_W:SpeciesGrass                      0.17800   
## TempTemp_W:SpeciesLathyrus                   0.95732   
## TempTemp_W:SpeciesPotentilla                 0.19900   
## TempTemp_W:SpeciesPseudocymopsos             0.35077   
## TempTemp_W:SpeciesSenecio                    0.67073   
## TempTemp_W:SpeciesTeraxicum                  0.34589   
## TempTemp_W:SpeciesValerianaF                 0.34923   
## TempTemp_W:SpeciesValerianaM                 0.78790   
## SnowSnow_M:TempTemp_W:SpeciesDephinium       0.97974   
## SnowSnow_M:TempTemp_W:SpeciesDugaldia        0.58495   
## SnowSnow_M:TempTemp_W:SpeciesErigeron        0.93120   
## SnowSnow_M:TempTemp_W:SpeciesFrageria        0.69094   
## SnowSnow_M:TempTemp_W:SpeciesGaliium         0.37517   
## SnowSnow_M:TempTemp_W:SpeciesGrass           0.26525   
## SnowSnow_M:TempTemp_W:SpeciesLathyrus        0.37284   
## SnowSnow_M:TempTemp_W:SpeciesPotentilla      0.43587   
## SnowSnow_M:TempTemp_W:SpeciesPseudocymopsos  0.99272   
## SnowSnow_M:TempTemp_W:SpeciesSenecio         0.41708   
## SnowSnow_M:TempTemp_W:SpeciesTeraxicum       0.51008   
## SnowSnow_M:TempTemp_W:SpeciesValerianaF      0.10573   
## SnowSnow_M:TempTemp_W:SpeciesValerianaM      0.65230   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 
## Correlation matrix not shown by default, as p = 56 > 12.
## Use print(x, correlation=TRUE)  or
##     vcov(x)        if you need it

## optimizer (Nelder_Mead) convergence code: 4 (failure to converge in 10000 evaluations)
## unable to evaluate scaled gradient
## Model failed to converge: degenerate  Hessian with 1 negative eigenvalues
## failure to converge in 10000 evaluations

res <- simulateResiduals(GNA.model2)
plot(res)#this cheks the scaled residuals with the fitted, looking for a strong deviation

overdispersion_test <- sum(residuals(GNA.model2, type = "pearson")^2) / df.residual(GNA.model2) #this checkes for overdispersion (meaning that there's a high amount of variance for a binomial mode) Values >> 1 suggest overdispersion

Table 3. Effects of early snowmelt and warming on S. exigua mass

ANOVA Table (type 2)
	F	Df	Df.res	Pr(>F)
Snow	5.3113	1	3.9966	0.0826
Temp	5.9627	1	242.6669	0.0153
Species	11.9429	13	254.5648	0.0000

Table 4. Effects of early snowmelt and warming on S. exigua survival

## Warning in vcov.merMod(mod, complete = FALSE): variance-covariance matrix computed from finite-difference Hessian is
## not positive definite or contains NA values: falling back to var-cov estimated from RX

ANOVA Table (type 3)
	Chisq	Df	Pr(>Chisq)
(Intercept)	5.6927	1	0.0170
Snow	0.2873	1	0.5919
Temp	1.1061	1	0.2929
Species	35.5369	13	0.0007
Snow:Temp	0.3180	1	0.5728
Snow:Species	29.1932	13	0.0061
Temp:Species	9.1244	13	0.7635
Snow:Temp:Species	10.6865	13	0.6371

Fig 1d. warming Treatment vs Mass

Fig 1c. snow Treatment vs Mass

Descriptive information

(mass) Absolute differences in treatments

percent decrease from warming

## [1] 22.76977

percent decrease from snowmelt

## [1] 22.45524

(mass) Fold change among the species

## [1] 31.86979

Supplementary

Model - Mass with 3 way interactions

# Create design matrix of predictors (without interaction expansion)
GNA.weights$Species <- as.factor(GNA.weights$Species)
GNA.weights$Snow <- as.factor(GNA.weights$Snow)
GNA.weights$Temp <- as.factor(GNA.weights$Temp)

GNA.lasso <- glmmLasso(fix = Mass**.5 ~ as.factor(Species) +as.factor(Snow) + as.factor(Temp) ,
  rnd = list(Plot = ~1, Plant.ID = ~1), data = GNA.weights, family = gaussian(link = "identity"), lambda = 10)

summary(GNA.lasso)

## Call:
## glmmLasso(fix = Mass^0.5 ~ as.factor(Species) + as.factor(Snow) + 
##     as.factor(Temp), rnd = list(Plot = ~1, Plant.ID = ~1), data = GNA.weights, 
##     lambda = 10, family = gaussian(link = "identity"))
## 
## 
## Fixed Effects:
## 
## Coefficients:
##                                   Estimate StdErr z.value p.value
## (Intercept)                       3.688161     NA      NA      NA
## as.factor(Species)Dephinium      -2.477432     NA      NA      NA
## as.factor(Species)Dugaldia       -0.017403     NA      NA      NA
## as.factor(Species)Erigeron       -1.474856     NA      NA      NA
## as.factor(Species)Frageria       -0.983447     NA      NA      NA
## as.factor(Species)Galiium        -0.731830     NA      NA      NA
## as.factor(Species)Grass          -1.283360     NA      NA      NA
## as.factor(Species)Lathyrus       -0.593030     NA      NA      NA
## as.factor(Species)Potentilla     -0.540102     NA      NA      NA
## as.factor(Species)Pseudocymopsos  0.245578     NA      NA      NA
## as.factor(Species)Senecio        -0.403596     NA      NA      NA
## as.factor(Species)Teraxicum       1.153780     NA      NA      NA
## as.factor(Species)ValerianaF     -1.424272     NA      NA      NA
## as.factor(Species)ValerianaM     -1.445819     NA      NA      NA
## as.factor(Snow)Snow_M            -0.067440     NA      NA      NA
## as.factor(Temp)Temp_W            -0.265291     NA      NA      NA
## 
## Random Effects:
## 
## StdDev:
## [[1]]
##           Plot
## Plot 0.2681851
## 
## [[2]]
##           Plant.ID
## Plant.ID 0.4135637

GNA.model5 <- lmer(Mass**.5 ~ Snow*Temp*Species + (1|Plot) + (1|Plant.ID), 
                   data = GNA.weights)

## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient

par(mfrow=c(1,2))
shapiro.test(resid(GNA.model5))

## 
##  Shapiro-Wilk normality test
## 
## data:  resid(GNA.model5)
## W = 0.98464, p-value = 5.36e-06

hist(resid(GNA.model5))
qqnorm(resid(GNA.model5))
qqline(resid(GNA.model5))#non-normal, so transforming

Model - Mass with 2 way interactions

GNA.model6 <- lmer(Mass**.5 ~ Snow+Temp+Species + Species*Temp + Species*Snow+ Snow*Temp+ (1|Plot) + (1|Plant.ID), 
                   data = GNA.weights)

par(mfrow=c(1,2))
shapiro.test(resid(GNA.model5))

## 
##  Shapiro-Wilk normality test
## 
## data:  resid(GNA.model5)
## W = 0.98464, p-value = 5.36e-06

hist(resid(GNA.model5))
qqnorm(resid(GNA.model5))
qqline(resid(GNA.model5))#non-normal, so transforming

Table S5. Effects of early snowmelt, warming, species, and their three-way interactions on S. exigua mass

ANOVA Table (type 3)
	F	Df	Df.res	Pr(>F)
(Intercept)	55.875	1	233.680	0.000
Snow	2.689	1	234.143	0.102
Temp	0.723	1	231.105	0.396
Species	4.126	13	219.088	0.000
Snow:Temp	0.024	1	236.966	0.878
Snow:Species	1.224	13	224.694	0.263
Temp:Species	1.024	13	217.898	0.429
Snow:Temp:Species	1.120	12	213.366	0.345

Table S6. Effects of early snowmelt, warming, species, and their two-way interactions on S. exigua mass

ANOVA Table (type 3)
	F	Df	Df.res	Pr(>F)
(Intercept)	75.406	1	223.756	0.000
Snow	5.360	1	175.378	0.022
Temp	1.274	1	247.545	0.260
Species	5.987	13	232.131	0.000
Temp:Species	1.015	13	225.946	0.438
Snow:Species	1.102	13	233.351	0.358
Snow:Temp	0.161	1	221.275	0.689

Fig S6. gain across species

### Fig XX. raw gain across species

Fig S7. Species specific response to the snow treatment

## Warning in vcov.merMod(object, complete = FALSE, ...): variance-covariance matrix computed from finite-difference Hessian is
## not positive definite or contains NA values: falling back to var-cov estimated from RX

## NOTE: Results may be misleading due to involvement in interactions

## Species = Achillea:
##  contrast        odds.ratio       SE  df asymp.LCL asymp.UCL null z.ratio
##  Snow_C / Snow_M      0.916 3.70e-01 Inf    0.4132  2.00e+00    1  -0.216
##  p.value
##   0.8287
## 
## Species = Dephinium:
##  contrast        odds.ratio       SE  df asymp.LCL asymp.UCL null z.ratio
##  Snow_C / Snow_M    630.439 1.35e+05 Inf    0.0000 1.88e+185    1   0.030
##  p.value
##   0.9760
## 
## Species = Dugaldia:
##  contrast        odds.ratio       SE  df asymp.LCL asymp.UCL null z.ratio
##  Snow_C / Snow_M      2.868 1.14e+00 Inf    1.3150  6.00e+00    1   2.648
##  p.value
##   0.0081
## 
## Species = Erigeron:
##  contrast        odds.ratio       SE  df asymp.LCL asymp.UCL null z.ratio
##  Snow_C / Snow_M      0.376 2.70e-01 Inf    0.0934  2.00e+00    1  -1.375
##  p.value
##   0.1691
## 
## Species = Frageria:
##  contrast        odds.ratio       SE  df asymp.LCL asymp.UCL null z.ratio
##  Snow_C / Snow_M      2.445 1.00e+00 Inf    1.0990  5.00e+00    1   2.191
##  p.value
##   0.0284
## 
## Species = Galiium:
##  contrast        odds.ratio       SE  df asymp.LCL asymp.UCL null z.ratio
##  Snow_C / Snow_M      0.719 3.50e-01 Inf    0.2797  2.00e+00    1  -0.686
##  p.value
##   0.4929
## 
## Species = Grass:
##  contrast        odds.ratio       SE  df asymp.LCL asymp.UCL null z.ratio
##  Snow_C / Snow_M      0.764 3.30e-01 Inf    0.3274  2.00e+00    1  -0.623
##  p.value
##   0.5335
## 
## Species = Lathyrus:
##  contrast        odds.ratio       SE  df asymp.LCL asymp.UCL null z.ratio
##  Snow_C / Snow_M      0.413 1.70e-01 Inf    0.1857  1.00e+00    1  -2.169
##  p.value
##   0.0301
## 
## Species = Potentilla:
##  contrast        odds.ratio       SE  df asymp.LCL asymp.UCL null z.ratio
##  Snow_C / Snow_M      0.813 3.10e-01 Inf    0.3821  2.00e+00    1  -0.537
##  p.value
##   0.5914
## 
## Species = Pseudocymopsos:
##  contrast        odds.ratio       SE  df asymp.LCL asymp.UCL null z.ratio
##  Snow_C / Snow_M      0.799 3.30e-01 Inf    0.3526  2.00e+00    1  -0.538
##  p.value
##   0.5905
## 
## Species = Senecio:
##  contrast        odds.ratio       SE  df asymp.LCL asymp.UCL null z.ratio
##  Snow_C / Snow_M      1.468 6.20e-01 Inf    0.6405  3.00e+00    1   0.907
##  p.value
##   0.3642
## 
## Species = Teraxicum:
##  contrast        odds.ratio       SE  df asymp.LCL asymp.UCL null z.ratio
##  Snow_C / Snow_M      0.871 3.40e-01 Inf    0.4076  2.00e+00    1  -0.357
##  p.value
##   0.7210
## 
## Species = ValerianaF:
##  contrast        odds.ratio       SE  df asymp.LCL asymp.UCL null z.ratio
##  Snow_C / Snow_M      3.522 1.46e+00 Inf    1.5672  8.00e+00    1   3.047
##  p.value
##   0.0023
## 
## Species = ValerianaM:
##  contrast        odds.ratio       SE  df asymp.LCL asymp.UCL null z.ratio
##  Snow_C / Snow_M      1.228 4.80e-01 Inf    0.5737  3.00e+00    1   0.529
##  p.value
##   0.5970
## 
## Results are averaged over the levels of: Temp 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log odds ratio scale 
## Tests are performed on the log odds ratio scale

## Warning in vcov.merMod(object, complete = FALSE, ...): variance-covariance matrix computed from finite-difference Hessian is
## not positive definite or contains NA values: falling back to var-cov estimated from RX

## NOTE: Results may be misleading due to involvement in interactions

Fig XX. raw Species specific response to the snow treatment

GNA.model2 <- glmer(Survival_binary ~  Snow * Temp * Species +(1|Plant.ID) + (1|Plot) , 
                    data = GNA.survive.binary, 
                    family = binomial(link = "logit"))

## Warning in (function (fn, par, lower = rep.int(-Inf, n), upper = rep.int(Inf, :
## failure to converge in 10000 evaluations

## Warning in optwrap(optimizer, devfun, start, rho$lower, control = control, :
## convergence code 4 from Nelder_Mead: failure to converge in 10000 evaluations

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## unable to evaluate scaled gradient

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge: degenerate Hessian with 1 negative eigenvalues

ggplot(GNA.survive.binary, aes(x = Snow, y = Survival_binary, fill=Snow))  +
  geom_violin(fill = "#BFD8B8", color = "black", alpha = 0.5, trim = FALSE) +  # violin shape
  geom_boxplot(width = 0.1, fill = "white", color = "black", outlier.shape = NA) +  # boxplot inside violin
  stat_summary(fun = mean, geom = "point", shape = 18, size = 3, color = "red") +
  facet_wrap(~Species,labeller = labeller(Species = label_map), scales = "free_x") +  # Creates separate plots per species
  labs(
       x = "", 
       y = "Proportion of S. Exigua survived") + 
  geom_segment(data = brackets_df,
               aes(x = x_start, xend = x_end, y = y_bracket, yend = y_bracket),
               inherit.aes = FALSE,
               linewidth = 0.6) +
   geom_text(data = sig_labels,
            aes(x = 1.5, y = y_pos, label = label),  # x = 1.5 centers between 2 bars
            inherit.aes = FALSE,
            size = 6)+
  theme_minimal() +
  theme(axis.text.x = element_text(hjust = 0.5, vjust = 0.5, size = 14),  
    axis.text.y = element_text(size = 14),  
    axis.title.x = element_text(size = 16, face = "bold"),  
    axis.title.y = element_text(size = 16, face = "bold"),  
    strip.text = element_text(size = 14),
    strip.background = element_rect(fill = "gray90", color = NA),
    panel.grid.major = element_blank(),  
    panel.grid.minor = element_blank(),  
    panel.border = element_blank(),  
    axis.line = element_line(color = "black"),
    legend.position = "none")+  # Remove legend) +  # Rotate x-axis labels
  scale_x_discrete(labels = c("Snow_C" = "Control",
                              "Snow_M" = "Snowmelt"))+
  scale_fill_manual(values = c("Snow_C" = "#BFD8B8",  
                              "Snow_M" = "#BFD8B8"))

Fig XX. for supp - raw warming Treatment vs Mass

(survival)

## # A tibble: 1 × 2
##   prop_survived     n
##           <dbl> <int>
## 1         0.382  1585

##   meansurvival  n
## 1    0.3711912 14

Unused code

we could treat snowmelt as continuous, but qualitatively produced similar.

Model - Mass with main effects (no interactions) and the snow treatment as a continuous variable

## 
##  Shapiro-Wilk normality test
## 
## data:  resid(model)
## W = 0.98856, p-value = 0.5504

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: sqrt(percent_herb)
##                                   Chisq Df Pr(>Chisq)  
## plant_species                    9.2801  4    0.05447 .
## plant_elev                       6.2765  1    0.01223 *
## gh_elev                          4.4930  1    0.03403 *
## plant_species:plant_elev         3.5825  4    0.46545  
## plant_species:gh_elev            0.3763  4    0.98437  
## plant_elev:gh_elev               1.6460  1    0.19950  
## plant_species:plant_elev:gh_elev 3.2220  4    0.52139  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 
##  Shapiro-Wilk normality test
## 
## data:  resid(model)
## W = 0.97475, p-value = 0.05156

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: percent_herb
##                                   Chisq Df Pr(>Chisq)   
## plant_species                    7.6569  4   0.104986   
## plant_elev                       4.0478  1   0.044229 * 
## gh_elev                          7.1901  1   0.007331 **
## plant_species:plant_elev         3.9902  4   0.407337   
## plant_species:gh_elev            0.6596  4   0.956212   
## plant_elev:gh_elev               1.5246  1   0.216932   
## plant_species:plant_elev:gh_elev 3.3385  4   0.502848   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 
##  Shapiro-Wilk normality test
## 
## data:  resid(model)
## W = 0.98917, p-value = 0.5979

## `summarise()` has grouped output by 'plant_elev'. You can override using the
## `.groups` argument.

## `summarise()` has grouped output by 'gh_elev'. You can override using the
## `.groups` argument.

## `summarise()` has grouped output by 'plant_species'. You can override using the
## `.groups` argument.

## `summarise()` has grouped output by 'plant_species'. You can override using the
## `.groups` argument.

## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'

plasticity_resistance_analyses_MS

Lydia Dean

2025-02-02

Exp. 1 Elevation survey

Data wrangling

Model - hopper num v elevation & date

Model - lack of fit testing - abundance

Model - hopper biomass v elevation & date

Model - lack of fit testing - biomass

Table 1. Analysis of Grasshopper abundance and biomass by Elevation

Fig 1a. Hopper vs Elevation

Supplementary (exp 1)

Model - hopper length v elevation & date

Table S2. Analysis of Grasshopper length by Elevation

Fig S1. Hopper biomass survey per date

Fig S2. Hopper biomass (across survey preiod)

Exp. 2 Elevation Feeding Trials

Data wrangling

Model - Percent herbivory no interactions

Table 2. Analysis of grasshopper herbivory

Fig 1b. Plant elevation main effect

Supplementary

Model - Percent herbivory with three way interaction

Model - Percent herbivory with two way interaction

Table S3. Percent herbivory with three way interaction

Table S4. Percent herbivory with two way interaction

Fig S3. Grasshopper elevation main effect

change to this : raw GH elevation

Fig S4. Species main effect

Fig S5. Interaction terms figure

Exp. 3 Climate Manipulation Feeding Trials

Data wrangling

Model - Mass with main effects (no interactions)

what about these models? (sqrt and glmm gamma log link)

Model - Survival (binary) with 3 way interaction

Table 3. Effects of early snowmelt and warming on S. exigua mass

Table 4. Effects of early snowmelt and warming on S. exigua survival

Fig 1d. warming Treatment vs Mass

Fig 1c. snow Treatment vs Mass

Descriptive information

(mass) Absolute differences in treatments

(mass) Fold change among the species

Supplementary

Model - Mass with 3 way interactions

Model - Mass with 2 way interactions

Table S5. Effects of early snowmelt, warming, species, and their three-way interactions on S. exigua mass

Table S6. Effects of early snowmelt, warming, species, and their two-way interactions on S. exigua mass

Fig S6. gain across species

Fig S7. Species specific response to the snow treatment

Fig XX. raw Species specific response to the snow treatment

Fig XX. for supp - raw warming Treatment vs Mass

Fig XX. for supp - raw warming Treatment vs Mass

(survival)

Unused code