1 Clay caterpillar analyses

Attack rate is measured as the proportion of 20 clay caterpillars on 10 trees that were attacked at a site on a unique date. There are 5 drainages, 4 sites within each, and either 3 to 4 sampling dates, resulting in 76 attack rates per site.

Bird diversity and abundance acquired from eBird along the same elevational gradient during the same time period of clay caterpillar data collection (Aug 2020), resulting in a total of 80 checklists, and 44 data points when averaged to the location due to uneven amounts of data per location. Diversity and abundance data account for sampling effort by using the residuals of each with time (duration length of each checklist).

1.1 Model selection using AIC (The smaller the AIC number the better the model)

## refitting model(s) with ML (instead of REML)
## Data: claycat
## Models:
## attack_AnovaR1: Propattack ~ elevation_m + pc1 + (1 | Site)
## attack_AnovaR2: Propattack ~ elevation_m + pc1 + (1 | Site) + (1 | DaysOut)
##                npar     AIC     BIC logLik deviance  Chisq Df Pr(>Chisq)   
## attack_AnovaR1    5 -49.093 -37.439 29.547  -59.093                        
## attack_AnovaR2    6 -55.203 -41.219 33.602  -67.203 8.1101  1   0.004402 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## refitting model(s) with ML (instead of REML)
## Data: claycat
## Models:
## attack_Anova0: Propattack ~ 1 + (1 | Site) + (1 | DaysOut)
## attack_Anova1: Propattack ~ elevation_m + (1 | Site) + (1 | DaysOut)
## attack_Anova2: Propattack ~ pc1 + (1 | Site) + (1 | DaysOut)
## attack_Anova3: Propattack ~ elevation_m + pc1 + (1 | Site) + (1 | DaysOut)
##               npar     AIC     BIC logLik deviance  Chisq Df Pr(>Chisq)   
## attack_Anova0    4 -52.226 -42.903 30.113  -60.226                        
## attack_Anova1    5 -52.172 -40.519 31.086  -62.172 1.9460  1   0.163016   
## attack_Anova2    5 -50.455 -38.801 30.228  -60.455 0.0000  0              
## attack_Anova3    6 -55.203 -41.219 33.602  -67.203 6.7481  1   0.009385 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

1.2 Models of proportion attack with elevation, aridity, and both

attack_Anova1 <- lmer(Propattack ~ elevation_m + (1 | Site) + (1 | DaysOut),
    data = claycat)  #Proportion attack rate at a unique site & date with elevation and site as a random factor
summary(attack_Anova1)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Propattack ~ elevation_m + (1 | Site) + (1 | DaysOut)
##    Data: claycat
## 
## REML criterion at convergence: -41.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.8139 -0.7004 -0.1366  0.5144  2.7178 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Site     (Intercept) 0.010478 0.10236 
##  DaysOut  (Intercept) 0.005257 0.07251 
##  Residual             0.018341 0.13543 
## Number of obs: 76, groups:  Site, 20; DaysOut, 3
## 
## Fixed effects:
##               Estimate Std. Error         df t value Pr(>|t|)
## (Intercept) -0.3271885  0.3374653 19.8293275  -0.970    0.344
## elevation_m  0.0001684  0.0001209 19.0771924   1.393    0.180
## 
## Correlation of Fixed Effects:
##             (Intr)
## elevation_m -0.989
attack_Anova2 <- lmer(Propattack ~ pc1 + (1 | Site) + (1 | DaysOut), data = claycat)  #Proportion attack rate at a unique site & date with aridity and site as a random factor
summary(attack_Anova2)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Propattack ~ pc1 + (1 | Site) + (1 | DaysOut)
##    Data: claycat
## 
## REML criterion at convergence: -49.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.8525 -0.6706 -0.1243  0.4969  2.7732 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Site     (Intercept) 0.011879 0.10899 
##  DaysOut  (Intercept) 0.005052 0.07107 
##  Residual             0.018379 0.13557 
## Number of obs: 76, groups:  Site, 20; DaysOut, 3
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)  
## (Intercept)  0.137398   0.050589  3.387301   2.716   0.0637 .
## pc1         -0.007655   0.016221 18.332822  -0.472   0.6426  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##     (Intr)
## pc1 0.013
attack_Anova3 <- lmer(Propattack ~ elevation_m + pc1 + (1 | Site) + (1 |
    DaysOut), data = claycat)  #Proportion attack rate at a unique site & date with elevation, aridity, and site as a random factor
summary(attack_Anova3)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Propattack ~ elevation_m + pc1 + (1 | Site) + (1 | DaysOut)
##    Data: claycat
## 
## REML criterion at convergence: -41.5
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.72129 -0.67314 -0.09072  0.52511  2.68508 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Site     (Intercept) 0.007601 0.08718 
##  DaysOut  (Intercept) 0.005606 0.07487 
##  Residual             0.018318 0.13534 
## Number of obs: 76, groups:  Site, 20; DaysOut, 3
## 
## Fixed effects:
##               Estimate Std. Error         df t value Pr(>|t|)  
## (Intercept) -1.9420893  0.7747507 18.3462858  -2.507   0.0218 *
## elevation_m  0.0007545  0.0002805 18.1336799   2.690   0.0149 *
## pc1          0.0815079  0.0359479 17.3930863   2.267   0.0364 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) elvtn_
## elevation_m -0.998       
## pc1         -0.919  0.921

This model removes the variation from elevation in a model with aridity and attack rates.

attack_model_A <- lmer(Propattack ~ elevation_m + (1 | Site), data = claycat)  #attack rate with elevation 
claycat$ele_residuals <- resid(attack_model_A)  #Getting the residuals that account for elevation

attack_model_A <- lmer(ele_residuals ~ pc1 + (1 | Site), data = claycat)  #model the attack rate/elevation residuals with aridity
## boundary (singular) fit: see help('isSingular')
summary(attack_model_A)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: ele_residuals ~ pc1 + (1 | Site)
##    Data: claycat
## 
## REML criterion at convergence: -71.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.0847 -0.7774 -0.2240  0.5773  3.1212 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Site     (Intercept) 0.00000  0.0000  
##  Residual             0.01947  0.1395  
## Number of obs: 76, groups:  Site, 20
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)
## (Intercept) 3.405e-04  1.602e-02 7.400e+01   0.021    0.983
## pc1         4.423e-03  8.918e-03 7.400e+01   0.496    0.621
## 
## Correlation of Fixed Effects:
##     (Intr)
## pc1 0.043 
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')

This model removes the variation from aridity in a model with elevation and attack rates.

attack_model_B <- lmer(Propattack ~ pc1 + (1 | Site), data = claycat)  #attack rate with aridity 
claycat$pc1_residuals <- resid(attack_model_B)  #Getting the residuals that account for aridity


attack_model_B <- lmer(pc1_residuals ~ elevation_m + (1 | Site), data = claycat)
## boundary (singular) fit: see help('isSingular')
# model the attack rate/aridity residuals with elevation
summary(attack_model_B)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: pc1_residuals ~ elevation_m + (1 | Site)
##    Data: claycat
## 
## REML criterion at convergence: -62.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.1776 -0.7704 -0.2233  0.6118  3.0949 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Site     (Intercept) 0.00000  0.000   
##  Residual             0.01932  0.139   
## Number of obs: 76, groups:  Site, 20
## 
## Fixed effects:
##               Estimate Std. Error         df t value Pr(>|t|)
## (Intercept) -1.063e-01  1.946e-01  7.400e+01  -0.546    0.587
## elevation_m  3.842e-05  7.009e-05  7.400e+01   0.548    0.585
## 
## Correlation of Fixed Effects:
##             (Intr)
## elevation_m -0.997
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')

1.3 New approach: including days exposed within the response variable

With this, for each sampling date, there is a CatDays number, which is the number of days total that caterpillars were exposed for between the sampling dates (where if 20 caterpillars were exposed for 3 days, there are 60 caterpillar days over that time). With this, you divide this number by the number of attacks over that time as the response variable.

catday_model1 = glmer(c(Totalattack/CatDays) ~ scale(pc1) + scale(elevation_m) +
    (1 | Site), data = claycat, family = "binomial", weights = CatDays)
summary(catday_model1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: c(Totalattack/CatDays) ~ scale(pc1) + scale(elevation_m) + (1 |  
##     Site)
##    Data: claycat
## Weights: CatDays
## 
##      AIC      BIC   logLik deviance df.resid 
##    541.4    550.7   -266.7    533.4       72 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.6891 -1.2736 -0.4005  2.3227  6.0255 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.2064   0.4543  
## Number of obs: 76, groups:  Site, 20
## 
## Fixed effects:
##                    Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         -3.1941     0.1260 -25.342  < 2e-16 ***
## scale(pc1)           0.7864     0.3256   2.415  0.01572 *  
## scale(elevation_m)   0.8806     0.3157   2.790  0.00528 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) scl(1)
## scale(pc1)  -0.040       
## scl(lvtn_m) -0.032  0.918
catday_model2 = glmer(c(Totalattack/CatDays) ~ scale(elevation_m) + (1 |
    Site), data = claycat, family = "binomial", weights = CatDays)
summary(catday_model2)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: c(Totalattack/CatDays) ~ scale(elevation_m) + (1 | Site)
##    Data: claycat
## Weights: CatDays
## 
##      AIC      BIC   logLik deviance df.resid 
##    544.4    551.4   -269.2    538.4       73 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.6177 -1.3165 -0.3686  2.2445  6.1669 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.3052   0.5525  
## Number of obs: 76, groups:  Site, 20
## 
## Fixed effects:
##                    Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         -3.2045     0.1447 -22.147   <2e-16 ***
## scale(elevation_m)   0.1737     0.1423   1.221    0.222    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## scl(lvtn_m) 0.024
catday_model3 = glmer(c(Totalattack/CatDays) ~ pc1 + (1 | Site), data = claycat,
    family = "binomial", weights = CatDays)
summary(catday_model3)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: c(Totalattack/CatDays) ~ pc1 + (1 | Site)
##    Data: claycat
## Weights: CatDays
## 
##      AIC      BIC   logLik deviance df.resid 
##    545.7    552.7   -269.9    539.7       73 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.5893 -1.3049 -0.3481  2.2031  6.2238 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.3379   0.5813  
## Number of obs: 76, groups:  Site, 20
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -3.21476    0.15032 -21.386   <2e-16 ***
## pc1         -0.02458    0.08429  -0.292    0.771    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##     (Intr)
## pc1 0.013
# plotSimulatedResiduals(simulateResiduals(m.total))
# Anova(catday_model1)

# predicted.total = predict(catday_model1,type='response') #this are
# predicted values already backtransformed in original scale (this
# can also be used to fix the scale after using the residuals in the
# graph ) claycat$pred.total = predicted.total


# ggplot(data = claycat, aes(x = elevation_m, y =pred.total)) +
# geom_point(shape = 21, fill = 'white', colour = 'black') +
# geom_smooth(method='lm', formula= y~x) + labs(x = 'Tree species
# richness\n(log-transformed)', y = 'Predation rate (%)', title = 'A)
# Total predation ')

# ggplot(data = claycat, aes(x = pc1, y =pred.total)) +
# geom_point(shape = 21, fill = 'white', colour = 'black') +
# geom_smooth(method='lm', formula= y~x) + labs(x = 'Tree species
# richness\n(log-transformed)', y = 'Predation rate (%)', title = 'A)
# Total predation ')

##Graphs of proportion attack with elevation and aridity

These are the proportion attack rates for each sampling date averaged to the site. (not using residuals)

## Joining, by = "Site"
## 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.
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'

These are the graphs for the attack rate for each date averaged to a site, using attack that either accounts for elevation or aridity (using the residuals of a model with each).

## Joining, by = "Site"
## Joining, by = "Site"
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'

1.4 Correlation of elevation and aridity scaled by average attack at a site (across sampling dates)

claycat_ave <- claycat %>%
    filter(!duplicated(Site))

attack_Anova <- lm(elevation_m ~ pc1, data = claycat)
summary(attack_Anova)
## 
## Call:
## lm(formula = elevation_m ~ pc1, data = claycat)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -173.292  -83.750    1.054   62.170  139.713 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 2758.311     10.271  268.56   <2e-16 ***
## pc1         -116.815      5.718  -20.43   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 89.46 on 74 degrees of freedom
## Multiple R-squared:  0.8494, Adjusted R-squared:  0.8474 
## F-statistic: 417.4 on 1 and 74 DF,  p-value: < 2.2e-16
## `geom_smooth()` using formula = 'y ~ x'

2 eBird analyses

2.1 Models of bird diversity with elevation, aridity, and both

eBird_div_model1 <- lm(div_residual_ave ~ elevation, data = ebird)
# ebird checklists: bird diversity accounting for sampling effort
# (time) averaged per location (elevation) with elevation
summary(eBird_div_model1)
## 
## Call:
## lm(formula = div_residual_ave ~ elevation, data = ebird)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.3076 -0.3296 -0.1126  0.5022  1.2143 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)
## (Intercept)  0.2204060  1.1545121   0.191    0.850
## elevation   -0.0001208  0.0004140  -0.292    0.772
## 
## Residual standard error: 0.5393 on 42 degrees of freedom
## Multiple R-squared:  0.002022,   Adjusted R-squared:  -0.02174 
## F-statistic: 0.08509 on 1 and 42 DF,  p-value: 0.7719
eBird_div_model2 <- lm(div_residual_ave ~ pc1, data = ebird)
# ebird checklists: bird diversity accounting for sampling effort
# (time) averaged per location (elevation) with aridity
summary(eBird_div_model2)
## 
## Call:
## lm(formula = div_residual_ave ~ pc1, data = ebird)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.2897 -0.3117 -0.0834  0.4790  1.0948 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.12236    0.08229  -1.487    0.145
## pc1         -0.01996    0.03995  -0.500    0.620
## 
## Residual standard error: 0.5383 on 42 degrees of freedom
## Multiple R-squared:  0.00591,    Adjusted R-squared:  -0.01776 
## F-statistic: 0.2497 on 1 and 42 DF,  p-value: 0.6199
eBird_div_model3 <- lm(div_residual_ave ~ elevation + pc1, data = ebird)
# ebird checklists: bird diversity accounting for sampling effort
# (time) averaged per location (elevation) with elevation and aridity
summary(eBird_div_model3)
## 
## Call:
## lm(formula = div_residual_ave ~ elevation + pc1, data = ebird)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.29885 -0.32448 -0.08318  0.51153  1.01352 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)
## (Intercept)  1.5133359  1.6839441   0.899    0.374
## elevation   -0.0005932  0.0006100  -0.973    0.337
## pc1         -0.0621350  0.0589807  -1.053    0.298
## 
## Residual standard error: 0.5386 on 41 degrees of freedom
## Multiple R-squared:  0.02832,    Adjusted R-squared:  -0.01907 
## F-statistic: 0.5976 on 2 and 41 DF,  p-value: 0.5549

2.2 Models of bird abundance with elevation, aridity, and both

eBird_abund_model1 <- lm(abun_residual_ave ~ elevation, data = ebird)
# ebird checklists: bird abundance accounting for sampling effort
# (time) averaged per location (elevation) with elevation
summary(eBird_abund_model1)
## 
## Call:
## lm(formula = abun_residual_ave ~ elevation, data = ebird)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -66.325 -13.287  -2.023  10.806 132.078 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)  
## (Intercept) 137.95822   65.89352   2.094   0.0424 *
## elevation    -0.05207    0.02363  -2.203   0.0331 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 30.78 on 42 degrees of freedom
## Multiple R-squared:  0.1036, Adjusted R-squared:  0.08228 
## F-statistic: 4.855 on 1 and 42 DF,  p-value: 0.0331
eBird_abund_model2 <- lm(abun_residual_ave ~ pc1, data = ebird)
# ebird checklists: bird abundance accounting for sampling effort
# (time) averaged per location (elevation) with aridity
summary(eBird_abund_model2)
## 
## Call:
## lm(formula = abun_residual_ave ~ pc1, data = ebird)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -60.216 -14.934  -4.885  10.531 133.753 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)   -6.066      4.913  -1.234    0.224
## pc1            2.361      2.386   0.989    0.328
## 
## Residual standard error: 32.14 on 42 degrees of freedom
## Multiple R-squared:  0.02278,    Adjusted R-squared:  -0.0004866 
## F-statistic: 0.9791 on 1 and 42 DF,  p-value: 0.3281
eBird_abund_model3 <- lm(abun_residual_ave ~ elevation + pc1, data = ebird)
# ebird checklists: bird abundance accounting for sampling effort
# (time) averaged per location (elevation) with elevation and aridity
summary(eBird_abund_model3)
## 
## Call:
## lm(formula = abun_residual_ave ~ elevation + pc1, data = ebird)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -64.997 -15.646  -2.435  10.194 131.165 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)  
## (Intercept) 198.67415   96.53002   2.058   0.0460 *
## elevation    -0.07425    0.03497  -2.124   0.0398 *
## pc1          -2.91786    3.38099  -0.863   0.3931  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 30.88 on 41 degrees of freedom
## Multiple R-squared:  0.1196, Adjusted R-squared:  0.07666 
## F-statistic: 2.785 on 2 and 41 DF,  p-value: 0.07342

2.3 Models accounting for elevation or ardity with bird diversity

eBird_Model_A <- lm(div_residual_ave ~ elevation, data = ebird)
# bird diversity (accounting for sampling effort) with elevation
ebird$ele_div_residuals <- resid(eBird_Model_A)  #Getting the residuals that account for elevation

eBird_Model_A <- lm(ele_div_residuals ~ pc1, data = ebird)  #model the diversity/elevation residuals with aridity
summary(eBird_Model_A)
## 
## Call:
## lm(formula = ele_div_residuals ~ pc1, data = ebird)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.29160 -0.32490 -0.08708  0.46757  1.07826 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.009747   0.081949  -0.119    0.906
## pc1         -0.028550   0.039790  -0.718    0.477
## 
## Residual standard error: 0.5361 on 42 degrees of freedom
## Multiple R-squared:  0.01211,    Adjusted R-squared:  -0.01141 
## F-statistic: 0.5149 on 1 and 42 DF,  p-value: 0.477
eBird_Model_C <- lm(div_residual_ave ~ pc1, data = ebird)
# bird diversity (accounting for sampling effort) with aridity
ebird$pc1_div_residuals <- resid(eBird_Model_C)  #Getting the residuals that account for aridity

eBird_Model_C <- lm(pc1_div_residuals ~ elevation, data = ebird)  #model the diversity/aridity residuals with elevation
summary(eBird_Model_C)
## 
## Call:
## lm(formula = pc1_div_residuals ~ elevation, data = ebird)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.3048 -0.3167 -0.1013  0.4867  1.1498 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)
## (Intercept)  0.7581991  1.1462765   0.661    0.512
## elevation   -0.0002726  0.0004111  -0.663    0.511
## 
## Residual standard error: 0.5355 on 42 degrees of freedom
## Multiple R-squared:  0.01036,    Adjusted R-squared:  -0.0132 
## F-statistic: 0.4397 on 1 and 42 DF,  p-value: 0.5109

2.4 Models accounting for elevation or ardity with bird abundance

eBird_Model_B <- lm(abun_residual_ave ~ elevation, data = ebird)
# bird abundance (accounting for sampling effort) with elevation
ebird$ele_abun_residuals <- resid(eBird_Model_B)  #Getting the residuals that account for elevation


eBird_Model_B <- lm(ele_abun_residuals ~ pc1, data = ebird)  #model the abundance/elevation residuals with aridity
summary(eBird_Model_B)
## 
## Call:
## lm(formula = ele_abun_residuals ~ pc1, data = ebird)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -63.568 -14.151  -3.502   9.853 131.938 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  -0.4577     4.6865  -0.098    0.923
## pc1          -1.3407     2.2755  -0.589    0.559
## 
## Residual standard error: 30.66 on 42 degrees of freedom
## Multiple R-squared:  0.008198,   Adjusted R-squared:  -0.01542 
## F-statistic: 0.3472 on 1 and 42 DF,  p-value: 0.5589
eBird_Model_D <- lm(abun_residual_ave ~ pc1, data = ebird)
# bird abundance (accounting for sampling effort) with aridity
ebird$pc1_abun_residuals <- resid(eBird_Model_D)  #Getting the residuals that account for aridity

eBird_Model_D <- lm(pc1_abun_residuals ~ elevation, data = ebird)  #model the abundance/aridity residuals with elevation
summary(eBird_Model_D)
## 
## Call:
## lm(formula = pc1_abun_residuals ~ elevation, data = ebird)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -67.400 -12.804  -2.388  10.937 132.816 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept) 94.90390   67.21610   1.412    0.165
## elevation   -0.03412    0.02410  -1.415    0.164
## 
## Residual standard error: 31.4 on 42 degrees of freedom
## Multiple R-squared:  0.04553,    Adjusted R-squared:  0.0228 
## F-statistic: 2.003 on 1 and 42 DF,  p-value: 0.1643

2.5 Graphs of eBirdf abundance or diversity with elevation and aridity

These do not use the residuals to account for elevation or aridity.

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

These do use the residuals to account for elevation or aridity.

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