Claycat_Analyses
Lydia Dean
2022-09-06
Predation rate of insectivorous birds across an elevational gradient
This file contains the statistical models, analyses, and visualizations.
Attack data
linear mixed models
Attack at a site not averaged across dates: climate alone
claycatModel4 <- lmer(Propattack ~ pc1+ (1| DaysOut), data = prop.attack)
Anova(claycatModel4)## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Propattack
## Chisq Df Pr(>Chisq)
## pc1 0.0293 1 0.864summary(claycatModel4)## Linear mixed model fit by REML ['lmerMod']
## Formula: Propattack ~ pc1 + (1 | DaysOut)
## Data: prop.attack
##
## REML criterion at convergence: -37.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.1976 -0.5850 -0.2581 0.2801 3.8726
##
## Random effects:
## Groups Name Variance Std.Dev.
## DaysOut (Intercept) 0.004089 0.06395
## Residual 0.030140 0.17361
## Number of obs: 76, groups: DaysOut, 3
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.142622 0.042541 3.353
## pc1 -0.003418 0.019959 -0.171
##
## Correlation of Fixed Effects:
## (Intr)
## pc1 0.007Attack at a site not averaged across dates: elevation alone
claycatModel5 <- lmer(Propattack ~ elevation_m+ (1| DaysOut), data = prop.attack)
Anova(claycatModel5)## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Propattack
## Chisq Df Pr(>Chisq)
## elevation_m 3.7273 1 0.05353 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(claycatModel5)## Linear mixed model fit by REML ['lmerMod']
## Formula: Propattack ~ elevation_m + (1 | DaysOut)
## Data: prop.attack
##
## REML criterion at convergence: -41.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.3802 -0.6276 -0.1898 0.4159 3.7949
##
## Random effects:
## Groups Name Variance Std.Dev.
## DaysOut (Intercept) 0.004831 0.06951
## Residual 0.028587 0.16908
## Number of obs: 76, groups: DaysOut, 3
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.14071 0.04513 3.118
## elevation_m 0.03770 0.01953 1.931
##
## Correlation of Fixed Effects:
## (Intr)
## elevation_m -0.013Attack at a site not averaged across dates: both climate and elevation
claycatModel6 <- lm(Propattack ~ elevation_m + pc1, data = prop.attack)
Anova(claycatModel6)## Anova Table (Type II tests)
##
## Response: Propattack
## Sum Sq Df F value Pr(>F)
## elevation_m 0.39057 1 13.927 0.0003735 ***
## pc1 0.30901 1 11.018 0.0014099 **
## Residuals 2.04731 73
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(claycatModel6)##
## Call:
## lm(formula = Propattack ~ elevation_m + pc1, data = prop.attack)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.24869 -0.13430 -0.01863 0.09804 0.49697
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.16118 0.01921 8.391 2.6e-12 ***
## elevation_m 0.16209 0.04343 3.732 0.000373 ***
## pc1 0.14417 0.04343 3.319 0.001410 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1675 on 73 degrees of freedom
## Multiple R-squared: 0.1602, Adjusted R-squared: 0.1372
## F-statistic: 6.965 on 2 and 73 DF, p-value: 0.001704Best model is with both elevation and climate
anova(claycatModel4,claycatModel5,claycatModel6)## refitting model(s) with ML (instead of REML)## Data: prop.attack
## Models:
## claycatModel4: Propattack ~ pc1 + (1 | DaysOut)
## claycatModel5: Propattack ~ elevation_m + (1 | DaysOut)
## claycatModel6: Propattack ~ elevation_m + pc1
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## claycatModel4 4 -40.355 -31.032 24.178 -48.355
## claycatModel5 4 -43.871 -34.548 25.936 -51.871 3.5159 0
## claycatModel6 4 -51.001 -41.678 29.500 -59.001 7.1296 0GLM
best glm with random and fixed effects
attack.model3 <- glmer(Attack.Binary ~ DaysOut + elevation_m + pc1 + (1|Drainage) + (1|Tree:Site:Drainage), data = CompleteData, family = "binomial") #elevation and climate
summary(attack.model3)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Attack.Binary ~ DaysOut + elevation_m + pc1 + (1 | Drainage) +
## (1 | Tree:Site:Drainage)
## Data: CompleteData
##
## AIC BIC logLik deviance df.resid
## 815.2 843.0 -401.6 803.2 754
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.3171 -0.5381 -0.4751 0.7592 2.8260
##
## Random effects:
## Groups Name Variance Std.Dev.
## Tree:Site:Drainage (Intercept) 0.1787 0.4228
## Drainage (Intercept) 0.2261 0.4755
## Number of obs: 760, groups: Tree:Site:Drainage, 200; Drainage, 5
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.65870 0.27629 -2.384 0.017122 *
## DaysOut -0.15039 0.04552 -3.304 0.000954 ***
## elevation_m 0.97897 0.41052 2.385 0.017093 *
## pc1 0.98275 0.38992 2.520 0.011722 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) DaysOt elvtn_
## DaysOut -0.517
## elevation_m 0.038 -0.027
## pc1 0.026 -0.032 0.850attack.model3 <- glmer(Attack.Binary ~ elevation_m + pc1 + (1|Drainage) + (1|Tree:Site:Drainage), data = CompleteData, family = "binomial") #elevation and climate
summary(attack.model3)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula:
## Attack.Binary ~ elevation_m + pc1 + (1 | Drainage) + (1 | Tree:Site:Drainage)
## Data: CompleteData
##
## AIC BIC logLik deviance df.resid
## 824.7 847.9 -407.4 814.7 755
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.1634 -0.5061 -0.4664 0.8196 2.1470
##
## Random effects:
## Groups Name Variance Std.Dev.
## Tree:Site:Drainage (Intercept) 0.1462 0.3824
## Drainage (Intercept) 0.2120 0.4605
## Number of obs: 760, groups: Tree:Site:Drainage, 200; Drainage, 5
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.1584 0.2298 -5.041 4.62e-07 ***
## elevation_m 0.9645 0.3995 2.414 0.0158 *
## pc1 0.9629 0.3792 2.539 0.0111 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) elvtn_
## elevation_m 0.026
## pc1 0.009 0.850logistic regressions
Logistic regressions of elevation, climate, and both. Plots each tree as attacked or not
#fit logistic regression model
attack.modela <- glm(Attack.Binary ~ elevation_m, data = CompleteData, family = binomial(link = "logit")) #elevation and climate
#define new data frame that contains predictor variable
newdata <- data.frame(elevation_m=seq(min(CompleteData$elevation_m), max(CompleteData$elevation_m),len=500))
#use fitted model to predict values of vs
newdata$Attack.Binary = predict(attack.modela, newdata, type="response")
#plot logistic regression curve
plot(Attack.Binary ~ elevation_m , data = CompleteData, col="steelblue")
lines(Attack.Binary ~ elevation_m , newdata, lwd=2)
eBird models
Using diversity and abundance across an elevational gradient
eBird data:
Accounting for variation- less than 5 km, last 10 years, fewer than 10 observers, less than 5 hours Filters- Gunnison, July & Aug 2020, Traveling/Stationary, Latitude, elevation >2400 <3300, birds with diet >40% insectivorous. https://www.google.com/maps/d/u/0/edit?mid=1Rs-jSKK1LY4ugp40EZVwcb3SUd4YIq8&ll=38.91495734634999%2C-106.96873191794234&z=9 (map of checklist locations and clay caterpillar sites)
4 models using residuals
The residuals of diversity in a model with time
ebirdModel4 <- lm(div_residual_ave ~ elevation, data = ebird_claycat_finalA) #model the residuals with elevation
summary(ebirdModel4)##
## Call:
## lm(formula = div_residual_ave ~ elevation, data = ebird_claycat_finalA)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.3578 -0.3038 0.1509 0.1725 1.2992
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.5951628 0.6809606 3.811 0.000275 ***
## elevation -0.0009539 0.0002496 -3.822 0.000265 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4523 on 78 degrees of freedom
## Multiple R-squared: 0.1577, Adjusted R-squared: 0.1469
## F-statistic: 14.6 on 1 and 78 DF, p-value: 0.000265## `geom_smooth()` using formula 'y ~ x'
Residuals of abundance with time
ebirdModel5 <- lm(abun_residual_ave ~ elevation, data = ebird_claycat_finalA) #model the residuals with elevation
summary(ebirdModel5)##
## Call:
## lm(formula = abun_residual_ave ~ elevation, data = ebird_claycat_finalA)
##
## Residuals:
## Min 1Q Median 3Q Max
## -71.574 -9.810 5.226 8.226 130.654
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 211.47327 36.40108 5.810 1.30e-07 ***
## elevation -0.07773 0.01334 -5.826 1.21e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 24.18 on 78 degrees of freedom
## Multiple R-squared: 0.3032, Adjusted R-squared: 0.2943
## F-statistic: 33.94 on 1 and 78 DF, p-value: 1.212e-07## `geom_smooth()` using formula 'y ~ x'
Diversity divided by abundance, in a model with time, and plotted the residuals.
ebirdModel6 <- lm(div_abund_residual_ave ~ elevation, data = ebird_claycat_finalA) #model the residuals with elevation
summary(ebirdModel6)##
## Call:
## lm(formula = div_abund_residual_ave ~ elevation, data = ebird_claycat_finalA)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.091120 -0.039162 -0.005346 0.026901 0.268234
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.617e-01 9.735e-02 -4.743 9.35e-06 ***
## elevation 1.697e-04 3.568e-05 4.756 8.89e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.06465 on 78 degrees of freedom
## Multiple R-squared: 0.2248, Adjusted R-squared: 0.2149
## F-statistic: 22.62 on 1 and 78 DF, p-value: 8.887e-06## `geom_smooth()` using formula 'y ~ x'
Residuals in a model with abundance and time
ebirdModel7 <- lm(residuals_ave ~ elevation, data = ebird_claycat_finalA) #model the residuals with elevation
summary(ebirdModel7)##
## Call:
## lm(formula = residuals_ave ~ elevation, data = ebird_claycat_finalA)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.2797 -0.2834 0.1056 0.1899 1.0256
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.8756064 0.5850995 1.497 0.139
## elevation -0.0003218 0.0002145 -1.501 0.137
##
## Residual standard error: 0.3886 on 78 degrees of freedom
## Multiple R-squared: 0.02806, Adjusted R-squared: 0.0156
## F-statistic: 2.252 on 1 and 78 DF, p-value: 0.1375## `geom_smooth()` using formula 'y ~ x'
Accounting for variation in sampling effort
Model that log transformed effort and offset in the model.
ebirdModela<- glm(Shannon ~ elevation + offset(duration_minutes), data = ebird_claycat_finalA)
summary(ebirdModela)##
## Call:
## glm(formula = Shannon ~ elevation + offset(duration_minutes),
## data = ebird_claycat_finalA)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -142.185 -35.200 6.359 34.610 78.613
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -148.24258 72.50470 -2.045 0.0443 *
## elevation 0.02568 0.02658 0.966 0.3370
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 2318.799)
##
## Null deviance: 183031 on 79 degrees of freedom
## Residual deviance: 180866 on 78 degrees of freedom
## AIC: 850.91
##
## Number of Fisher Scoring iterations: 2Anova(ebirdModela)## Analysis of Deviance Table (Type II tests)
##
## Response: Shannon
## LR Chisq Df Pr(>Chisq)
## elevation 0.93341 1 0.334ADD MOST RECENT
add the resulting climate variable for the ebird data show the pca and the 2 regressions
show the four plots of proportion models with reciprocal residuals
# averaged between dates
claycatmodelA <- lm(Ave.Prop.Date ~ pc1 , data = ave.prop.attack) #diversity/abundance with distance
ave.prop.attack$pc1_residuals <- resid(claycatmodelA)
summary(claycatmodelA$fitted.values)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.1576 0.1577 0.1585 0.1596 0.1616 0.1641claycatmodelA <- lm(pc1_residuals ~ elevation_m, data = ave.prop.attack) #model the residuals with elevation
summary(claycatmodelA)##
## Call:
## lm(formula = pc1_residuals ~ elevation_m, data = ave.prop.attack)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.14759 -0.06388 -0.01709 0.03531 0.41475
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.00196 0.02745 0.071 0.944
## elevation_m 0.03013 0.02700 1.116 0.279
##
## Residual standard error: 0.1225 on 18 degrees of freedom
## Multiple R-squared: 0.06472, Adjusted R-squared: 0.01276
## F-statistic: 1.245 on 1 and 18 DF, p-value: 0.2791ggplot(data = ave.prop.attack, aes(x = elevation_m, y =pc1_residuals ))+
geom_point(size=6) +
geom_smooth(colour="black",size =2,show.legend = F,method= "lm", se=FALSE)+
theme(text = element_text(size=25))+
labs(x= "Elevation", y = "residuals (div ~ distance)") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())## `geom_smooth()` using formula 'y ~ x'
# averaged between dates
claycatmodelA <- lm(Ave.Prop.Date ~ elevation_m , data = ave.prop.attack) #diversity/abundance with distance
ave.prop.attack$ele_residuals <- resid(claycatmodelA)
summary(claycatmodelA$fitted.values)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.1174 0.1339 0.1516 0.1596 0.1889 0.2163claycatmodelA <- lm(ele_residuals ~ pc1, data = ave.prop.attack) #model the residuals with elevation
summary(claycatmodelA)##
## Call:
## lm(formula = ele_residuals ~ pc1, data = ave.prop.attack)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.13985 -0.05216 -0.03010 0.02246 0.41130
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.0009368 0.0265377 -0.035 0.972
## pc1 0.0272494 0.0269111 1.013 0.325
##
## Residual standard error: 0.1186 on 18 degrees of freedom
## Multiple R-squared: 0.05389, Adjusted R-squared: 0.00133
## F-statistic: 1.025 on 1 and 18 DF, p-value: 0.3247ggplot(data = ave.prop.attack, aes(x = pc1, y =ele_residuals ))+
geom_point(size=6) +
geom_smooth(colour="black",size =2,show.legend = F,method= "lm", se=FALSE)+
theme(text = element_text(size=25))+
labs(x= "Elevation", y = "residuals (div ~ distance)") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())## `geom_smooth()` using formula 'y ~ x'
#not averaged between dates
claycatmodelA <- lm(Propattack ~ elevation_m , data = prop.attack) #diversity/abundance with distance
prop.attack$ele_residuals <- resid(claycatmodelA)
summary(claycatmodelA$fitted.values)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.1160 0.1331 0.1555 0.1612 0.1946 0.2170prop.attackx<- prop.attack %>%
group_by(elevation_m) %>%
summarise(ele_residual_ave = mean(ele_residuals)) #average the residuals across location (elevation)
prop.attack<-left_join(prop.attack, prop.attackx)## Joining, by = "elevation_m"claycatmodelA <- lm(ele_residual_ave ~ pc1, data = prop.attack) #model the residuals with elevation
summary(claycatmodelA)##
## Call:
## lm(formula = ele_residual_ave ~ pc1, data = prop.attack)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.14053 -0.05205 -0.03142 0.02172 0.41020
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.823e-18 1.336e-02 0.000 1.0000
## pc1 2.824e-02 1.337e-02 2.112 0.0381 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1165 on 74 degrees of freedom
## Multiple R-squared: 0.05685, Adjusted R-squared: 0.04411
## F-statistic: 4.461 on 1 and 74 DF, p-value: 0.03806ggplot(data = prop.attack, aes(x = pc1, y =ele_residual_ave ))+
geom_point(size=6) +
geom_smooth(colour="black",size =2,show.legend = F,method= "lm", se=FALSE)+
theme(text = element_text(size=25))+
labs(x= "Elevation", y = "residuals (div ~ distance)") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())## `geom_smooth()` using formula 'y ~ x'
#not averaged between dates
claycatmodelA <- lm(Propattack ~ pc1 , data = prop.attack) #diversity/abundance with distance
summary(claycatmodelA)##
## Call:
## lm(formula = Propattack ~ pc1, data = prop.attack)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.16286 -0.11228 -0.06149 0.08817 0.63893
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.161184 0.020820 7.742 3.98e-11 ***
## pc1 -0.001175 0.020834 -0.056 0.955
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1815 on 74 degrees of freedom
## Multiple R-squared: 4.3e-05, Adjusted R-squared: -0.01347
## F-statistic: 0.003182 on 1 and 74 DF, p-value: 0.9552Anova(claycatmodelA)## Anova Table (Type II tests)
##
## Response: Propattack
## Sum Sq Df F value Pr(>F)
## pc1 0.0001 1 0.0032 0.9552
## Residuals 2.4379 74prop.attack$pc1_residuals <- resid(claycatmodelA)
summary(claycatmodelA$fitted.values)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.1601 0.1602 0.1608 0.1612 0.1622 0.1635prop.attackx<- prop.attack %>%
group_by(elevation_m) %>%
summarise(pc1_residual_ave = mean(pc1_residuals)) #average the residuals across location (elevation)
prop.attack<-left_join(prop.attack, prop.attackx)## Joining, by = "elevation_m"claycatmodelA <- lm(pc1_residual_ave ~ elevation_m, data = prop.attack) #model the residuals with elevation
summary(claycatmodelA)##
## Call:
## lm(formula = pc1_residual_ave ~ elevation_m, data = prop.attack)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.14737 -0.07021 -0.02473 0.03567 0.41381
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.971e-18 1.379e-02 0.0 1.0000
## elevation_m 3.175e-02 1.380e-02 2.3 0.0242 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1202 on 74 degrees of freedom
## Multiple R-squared: 0.06674, Adjusted R-squared: 0.05413
## F-statistic: 5.292 on 1 and 74 DF, p-value: 0.02424ggplot(data = prop.attack, aes(x = elevation_m, y =pc1_residual_ave ))+
geom_point(size=6) +
geom_smooth(colour="black",size =2,show.legend = F,method= "lm", se=FALSE)+
theme(text = element_text(size=25))+
labs(x= "Elevation", y = "residuals (div ~ distance)") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())## `geom_smooth()` using formula 'y ~ x'
climate data w