Code for ANOVAs, PLS, and generalized linear modeling
Gunnar R. Kramer
2023-01-12
Load libraries:
library(devtools)
library(tidyverse)
library(caret)
library(pls)
library(usdm)
library(raster)
library(maptools)
library(rgeos)
library(rgdal)
library(maps)
library(AICcmodavg)Maps of migration risk-factors, boxplots, and regional ANOVAs
Maps of migration risk-factors (Figure S4)
#Load risk-factor rasters (standardized extent, resolution, etc.)
setwd("/Users/gunnarkramer/Desktop/FINAL DATA FOR DRUM LAND ECOL/Rasters/")
TurbineRast<-raster("/Users/gunnarkramer/Desktop/FINAL DATA FOR DRUM LAND ECOL/Rasters/TurbineRastR.tif") #wind turbine raster
HurrPtsRast<-raster("/Users/gunnarkramer/Desktop/FINAL DATA FOR DRUM LAND ECOL/Rasters/HurrPtsRastR.tif") #hurricane raster
TotalTorn<-raster("/Users/gunnarkramer/Desktop/FINAL DATA FOR DRUM LAND ECOL/Rasters/TotalTornRastR.tif") #tornado raster
CropCov<-raster("/Users/gunnarkramer/Desktop/FINAL DATA FOR DRUM LAND ECOL/Rasters/CropCovresfix.tif") #agricultural cover raster
ForestShrubCov<-raster("/Users/gunnarkramer/Desktop/FINAL DATA FOR DRUM LAND ECOL/Rasters/ForestShrubCovresfix.tif") #forest and shrub cover raster
HFootP<-raster("/Users/gunnarkramer/Desktop/FINAL DATA FOR DRUM LAND ECOL/Rasters/HFootPresfix.tif") #human footprint raster
TallTowers<-raster("/Users/gunnarkramer/Desktop/FINAL DATA FOR DRUM LAND ECOL/Rasters/TallTowerRastR.tif") #communications towers raster
Forest2010<-raster("/Users/gunnarkramer/Desktop/FINAL DATA FOR DRUM LAND ECOL/Rasters/Forest_2010Hyde.tif") #forest cover 2010
Forest2000<-raster("/Users/gunnarkramer/Desktop/FINAL DATA FOR DRUM LAND ECOL/Rasters/Forest_2000Hyde.tif") #forest cover 2000
#Standardize, set NA's to zero
TurbineRast<-reclassify(TurbineRast, cbind(NA, NA, 0), right=TRUE)
HurrPtsRast<-reclassify(HurrPtsRast, cbind(NA, NA, 0), right=TRUE)
TotalTorn<-reclassify(TotalTorn, cbind(NA, NA, 0), right=TRUE)
CropCov<-reclassify(CropCov, cbind(NA, NA, 0), right=TRUE)
ForestShrubCov<-reclassify(ForestShrubCov, cbind(NA, NA, 0), right=TRUE)
HFootP<-reclassify(HFootP, cbind(NA, NA, 0), right=TRUE)
TallTowers<-reclassify(TallTowers, cbind(NA, NA, 0), right=TRUE)Some housekeeping and loading of boundary shapefiles, ocean mask layer
data(wrld_simpl)
#extent
x <- list(x=-110:-60, y = 0:60, z = outer(1:15, 1:15, "+"))
## use raster to quickly generate the polymask
r <- ForestShrubCov #get extent
p <- as(extent(r), "SpatialPolygons")
#get mask of ocean for plots
wmap <- gIntersection(wrld_simpl, p)
oceanmap <- gDifference(p, wmap)
#test
image(r)
plot(oceanmap, add = TRUE, col = "light blue")
#works!Load boundary shapefiles (not necessary for anlaysis)
#Canadian provinces
canada<-readOGR("/Users/gunnarkramer/Desktop/FINAL DATA FOR DRUM LAND ECOL/Rasters/Canada.shp")## OGR data source with driver: ESRI Shapefile
## Source: "/Users/gunnarkramer/Desktop/FINAL DATA FOR DRUM LAND ECOL/Rasters/Canada.shp", layer: "Canada"
## with 13 features
## It has 2 fieldscanada<-spTransform(canada, CRS("+proj=longlat +datum=WGS84"))
plot(canada)
#Great Lakes shapefile
GLakes.shp<-readOGR("/Users/gunnarkramer/Desktop/FINAL DATA FOR DRUM LAND ECOL/Rasters/Great_Lakes.shp")## OGR data source with driver: ESRI Shapefile
## Source: "/Users/gunnarkramer/Desktop/FINAL DATA FOR DRUM LAND ECOL/Rasters/Great_Lakes.shp", layer: "Great_Lakes"
## with 1 features
## It has 4 fields
## Integer64 fields read as strings: OBJECTID dissolveGLakes.shp<-spTransform(GLakes.shp, CRS("+proj=longlat +datum=WGS84"))
plot(GLakes.shp)
#make color ramp
pinkgreen<-colorRampPalette(c("#005AB5", "white", "#DC3220"), bias=2)
revpinkgreen<-colorRampPalette(c("#DC3220", "white", "#005AB5"), bias=2)
map.cols<-pinkgreen(8)
rev.cols<-revpinkgreen(8)
plot(1:8, pch=16, col=map.cols)
plot(1:8, pch=16, col=rev.cols)
col3<-c("#DC3220", "white", "#005AB5")
#calculate change in forest cover
deltaFor<-Forest2010-Forest2000
cellStats(deltaFor, range)## [1] -1 1deltaFor<-deltaFor*-1 #correct for direction (more forest is better for warblers)
plot(deltaFor)
Plot spatial distribution of migration risk-factors, breeding factors, and nonbreeding factors
##RASTERS ONLY -Horizontal
######
#png("VermMigration_RASTERPANEL_Horizontal_SCALED_29Nov2021_ONESCALE.png", width = (6.85039*2), height = (6.85039+1), units = 'in', res = 600)
par(mfrow=c(2,4), mar=c(1,2,4,0), pty="s")
plot(1,1, type='n', ylim=c(0,60), xlim=c(-110,-60), ylab="", xlab="", yaxt='n', xaxt='n')
plot(scale(ForestShrubCov), add=TRUE, legend=FALSE, col=rev.cols)
plot(oceanmap, add=TRUE, col="gray", border=FALSE)
plot(GLakes.shp, add=TRUE, col="gray", border=FALSE)
plot(wrld_simpl, add=T, border = "grey40", lwd=.5)
maps::map('state', add=TRUE, col = "grey40", lwd=.5)
plot(canada, add=TRUE, border="grey40", lwd=.5)
mtext("Forest and shrub cover", side=3, line=.5, cex=1.25, font=2)
box(cex=1.25)
plot(1,1, type='n', ylim=c(0,60), xlim=c(-110,-60), ylab="", xlab="", yaxt='n', xaxt='n')
plot(scale(deltaFor), add=TRUE, legend=FALSE, col=rev(col3))
plot(oceanmap, add=TRUE, col="gray", border=FALSE)
plot(GLakes.shp, add=TRUE, col="gray", border=FALSE)
plot(wrld_simpl, add=T, border = "grey40", lwd=.5)
maps::map('state', add=TRUE, col = "grey40", lwd=.5)
plot(canada, add=TRUE, border="grey40", lwd=.5)
mtext("Net change in forest cover 2000-2015", side=3, line=.5, cex=1.25, font=2)
box(cex=1.25)
plot(1,1, type='n', ylim=c(0,60), xlim=c(-110,-60), ylab="", xlab="", yaxt='n', xaxt='n')
plot(scale(CropCov), add=TRUE, legend=FALSE, col=map.cols)
plot(oceanmap, add=TRUE, col="gray", border=FALSE)
plot(GLakes.shp, add=TRUE, col="gray", border=FALSE)
plot(wrld_simpl, add=T, border = "grey40", lwd=.5)
maps::map('state', add=TRUE, col = "grey40", lwd=.5)
plot(canada, add=TRUE, border="grey40", lwd=.5)
mtext("Agriculture cover", side=3, line=.5, cex=1.25, font=2)
box(cex=1.25)
plot(1,1, type='n', ylim=c(0,60), xlim=c(-110,-60), ylab="", xlab="", yaxt='n', xaxt='n')
plot(scale(HFootP), add=TRUE, legend=FALSE, col=map.cols)
plot(GLakes.shp, add=TRUE, col="gray", border=FALSE)
plot(oceanmap, add=TRUE, col="gray", border=FALSE)
plot(wrld_simpl, add=T, border = "grey40", lwd=.5)
maps::map('state', add=TRUE, col = "grey40", lwd=.5)
plot(canada, add=TRUE, border="grey40", lwd=.5)
mtext("Human footprint", side=3, line=.5, cex=1.25, font=2)
box(cex=1.25)
plot(1,1, type='n', ylim=c(0,60), xlim=c(-110,-60), ylab="", xlab="", yaxt='n', xaxt='n')
plot(scale(TurbineRast), add=TRUE, legend=FALSE, col=map.cols)
plot(GLakes.shp, add=TRUE, col="gray", border=FALSE)
plot(oceanmap, add=TRUE, col="gray", border=FALSE)
plot(wrld_simpl, add=T, border = "grey40", lwd=.5)
maps::map('state', add=TRUE, col = "grey40", lwd=.5)
plot(canada, add=TRUE, border="grey40", lwd=.5)
mtext("Wind energy", side=3, line=.5, cex=1.25, font=2)
box(cex=1.25)
plot(1,1, type='n', ylim=c(0,60), xlim=c(-110,-60), ylab="", xlab="", yaxt='n', xaxt='n')
plot(scale(TallTowers), add=TRUE, legend=FALSE, col=map.cols)
plot(oceanmap, add=TRUE, col="gray", border=FALSE)
plot(GLakes.shp, add=TRUE, col="gray", border=FALSE)
plot(wrld_simpl, add=T, border = "grey40", lwd=.5)
maps::map('state', add=TRUE, col = "grey40", lwd=.5)
plot(canada, add=TRUE, border="grey40", lwd=.5)
mtext("Communications towers", side=3, line=.5, cex=1.25, font=2)
box(cex=1.25)
plot(1,1, type='n', ylim=c(0,60), xlim=c(-110,-60), ylab="", xlab="", yaxt='n', xaxt='n')
plot(scale(TotalTorn), add=TRUE, legend=FALSE, col=map.cols)
plot(GLakes.shp, add=TRUE, col="gray", border=FALSE)
plot(oceanmap, add=TRUE, col="gray", border=FALSE)
plot(wrld_simpl, add=T, border = "grey40", lwd=.5)
maps::map('state', add=TRUE, col = "grey40", lwd=.5)
plot(canada, add=TRUE, border="grey40", lwd=.5)
mtext("Tornados", side=3, line=.5, cex=1.25, font=2)
box(cex=1.25)
plot(1,1, type='n', ylim=c(0,60), xlim=c(-110,-60), ylab="", xlab="", yaxt='n', xaxt='n')
plot(scale(HurrPtsRast), add=TRUE, legend=FALSE, col=map.cols)
plot(oceanmap, add=TRUE, col="gray", border=FALSE)
plot(GLakes.shp, add=TRUE, col="gray", border=FALSE)
plot(wrld_simpl, add=T, border = "grey40", lwd=.5)
maps::map('state', add=TRUE, col = "grey40", lwd=.5)
plot(canada, add=TRUE, border="grey40", lwd=.5)
mtext("Hurricanes", side=3, line=.5, cex=1.25, font=2)
box(cex=1.25)
#dev.off()
#####- Test for differences in exposure among regional populations (BCR); ANOVAs
#BOXPLOTS
#ANNUAL DATA
lassodata25<-read.csv("/Users/gunnarkramer/Desktop/FINAL DATA FOR DRUM LAND ECOL/Data4Lasso25_data.csv", header=TRUE)
lassodata25$Individual<-as.character(lassodata25$Individual)
lassodata25$Species<-as.character(lassodata25$Species)
lassodata25$Pop<-as.character(lassodata25$Pop)
lassodata25$Trend00<-as.numeric(as.character(lassodata25$Trend00))
lassodata25$TrendGrp<-ifelse(lassodata25$Trend00 <=-1.5, "declining", ifelse(lassodata25$Trend00>-1.5 & lassodata25$Trend00<=1.5, "stable", ifelse(lassodata25$Trend00>1.5, "increasing", "NA")))
#set up colors for plots
PHTcol<-rgb(0/255, 146/255, 146/255, 1)
CHcol<-rgb(255/255, 9/255, 231/255, 1)
bwAMcol<-rgb(255/255, 157/255, 9/255, .5)
BHTcol<-rgb(146/255,0/255,0/255, 1)
AMcol<-rgb(255/255, 157/255, 9/255, 1)
bplotcols<-c(PHTcol,CHcol,bwAMcol,BHTcol,AMcol)
#Omit hybrids
lassodata25<-subset(lassodata25, Species!="HYBRID")
#create species x population column
lassodata25$spPop<-paste(lassodata25$Species, lassodata25$Pop, sep="")
#create stationary and or increasing vs. declining variable
lassodata25$IorS<-ifelse(lassodata25$TrendGrp!="declining", 1, 0)
#check
colnames(lassodata25)## [1] "Individual" "Species" "Pop"
## [4] "State" "BrLat" "BrLon"
## [7] "NonBrLat" "NonBrLon" "BCR_Trend"
## [10] "Trend00" "scale" "Turbines"
## [13] "Towers" "Hurricanes" "Tornados"
## [16] "Agriculture" "Forest.cover" "Human.footprint"
## [19] "Net.forest.change.2000" "OVERALL2" "Br_Forest"
## [22] "Br_HFootP" "Br_forestchange" "Nb_Forest"
## [25] "Nb_HFootP" "Nb_forestchange" "TrendGrp"
## [28] "spPop" "IorS"#set order of population variable (factor)
lassodata25$factPop<-factor(lassodata25$Pop, levels=c("PHT", "CH", "bwAM", "BHT", "AM"))Run ANOVAs
#ANOVAs for 25 scale regional population-level comparisons
#Forest and shrub cover
aovForest<-aov(scale(Forest.cover)~Pop, data=lassodata25)
summary(aovForest)## Df Sum Sq Mean Sq F value Pr(>F)
## Pop 4 13.87 3.468 3.986 0.00546 **
## Residuals 76 66.13 0.870
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# P = 0.005
TukeyHSD(aovForest)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = scale(Forest.cover) ~ Pop, data = lassodata25)
##
## $Pop
## diff lwr upr p adj
## BHT-AM -1.1259378 -1.9466631 -0.3052126 0.0023414
## bwAM-AM -0.5569012 -1.6003459 0.4865434 0.5712061
## CH-AM -0.6915535 -2.1818906 0.7987836 0.6940586
## PHT-AM -0.8785158 -2.1004730 0.3434415 0.2716921
## bwAM-BHT 0.5690366 -0.2778052 1.4158784 0.3381129
## CH-BHT 0.4343843 -0.9255675 1.7943361 0.8988143
## PHT-BHT 0.2474221 -0.8116076 1.3064517 0.9656016
## CH-bwAM -0.1346523 -1.6395296 1.3702250 0.9991187
## PHT-bwAM -0.3216145 -1.5612639 0.9180349 0.9500801
## PHT-CH -0.1869622 -1.8206876 1.4467631 0.9976870#BHT < AM
#Net change in forest and shrub cover 2000-2010
aovNetFor<-aov(scale(Net.forest.change.2000)~Pop, data=lassodata25)
summary(aovNetFor)## Df Sum Sq Mean Sq F value Pr(>F)
## Pop 4 6.39 1.5979 1.65 0.171
## Residuals 76 73.61 0.9685# P = 0.171
#Agricultural cover
aovCrop<-aov(scale(Agriculture)~Pop, data=lassodata25)
summary(aovCrop)## Df Sum Sq Mean Sq F value Pr(>F)
## Pop 4 14 3.499 4.029 0.00512 **
## Residuals 76 66 0.868
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# P = 0.00512
TukeyHSD(aovCrop)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = scale(Agriculture) ~ Pop, data = lassodata25)
##
## $Pop
## diff lwr upr p adj
## BHT-AM 0.84792412 0.02796754 1.6678807 0.0391053
## bwAM-AM 0.24353423 -0.79893320 1.2860017 0.9656109
## CH-AM 0.91905360 -0.56988772 2.4079949 0.4251668
## PHT-AM 1.45484828 0.23403542 2.6756611 0.0113971
## bwAM-BHT -0.60438990 -1.45043855 0.2416588 0.2777796
## CH-BHT 0.07112947 -1.28754864 1.4298076 0.9998946
## PHT-BHT 0.60692415 -0.45111365 1.6649620 0.5001994
## CH-bwAM 0.67551937 -0.82794855 2.1789873 0.7188874
## PHT-bwAM 1.21131405 -0.02717437 2.4498025 0.0583156
## PHT-CH 0.53579468 -1.09640060 2.1679900 0.8893976#PHT > AM; BHT > AM
#Human footprint
aovHFootP<-aov(scale(Human.footprint)~Pop, data=lassodata25)
summary(aovHFootP)## Df Sum Sq Mean Sq F value Pr(>F)
## Pop 4 3.74 0.9349 0.932 0.45
## Residuals 76 76.26 1.0034# P = 0.45
#Wind turbines
aovTurbine<-aov(scale(Turbines)~Pop, data=lassodata25)
summary(aovTurbine)## Df Sum Sq Mean Sq F value Pr(>F)
## Pop 4 11.24 2.8105 3.107 0.0201 *
## Residuals 76 68.76 0.9047
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# P = 0.0201
TukeyHSD(aovTurbine)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = scale(Turbines) ~ Pop, data = lassodata25)
##
## $Pop
## diff lwr upr p adj
## BHT-AM 0.8424962 0.005601311 1.6793911 0.0476759
## bwAM-AM 0.4049250 -0.659077255 1.4689273 0.8245232
## CH-AM 1.4126122 -0.107086999 2.9323115 0.0809546
## PHT-AM 1.0569459 -0.189085993 2.3029778 0.1348014
## bwAM-BHT -0.4375712 -1.301097149 0.4259548 0.6194106
## CH-BHT 0.5701160 -0.816629084 1.9568611 0.7800103
## PHT-BHT 0.2144497 -0.865444616 1.2943440 0.9810071
## CH-bwAM 1.0076872 -0.526838717 2.5422131 0.3614603
## PHT-bwAM 0.6520209 -0.612051723 1.9160935 0.6032232
## PHT-CH -0.3556663 -2.021578825 1.3102461 0.9752003#BHT>AM
#Communications towers
aovTowers<-aov(scale(Towers)~Pop, data=lassodata25)
summary(aovTowers)## Df Sum Sq Mean Sq F value Pr(>F)
## Pop 4 1.14 0.2854 0.275 0.893
## Residuals 76 78.86 1.0376# P = 0.893
#Tornados
aovTorn<-aov(scale(Tornados)~Pop, data=lassodata25)
summary(aovTorn)## Df Sum Sq Mean Sq F value Pr(>F)
## Pop 4 15.82 3.955 4.684 0.00197 **
## Residuals 76 64.18 0.844
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# P = 0.00197
TukeyHSD(aovTorn)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = scale(Tornados) ~ Pop, data = lassodata25)
##
## $Pop
## diff lwr upr p adj
## BHT-AM 0.80021538 -0.008335409 1.6087662 0.0537706
## bwAM-AM 0.35740652 -0.670559942 1.3853730 0.8671190
## CH-AM 1.54093299 0.072703207 3.0091628 0.0349059
## PHT-AM 1.52496675 0.321135697 2.7287978 0.0060179
## bwAM-BHT -0.44280886 -1.277088772 0.3914711 0.5764787
## CH-BHT 0.74071761 -0.599060960 2.0804962 0.5369908
## PHT-BHT 0.72475137 -0.318568873 1.7680716 0.3049591
## CH-bwAM 1.18352647 -0.299027842 2.6660808 0.1796810
## PHT-bwAM 1.16756023 -0.053700512 2.3888210 0.0678232
## PHT-CH -0.01596624 -1.625457289 1.5935248 0.9999999#PHT > AM; CH> AM
#Hurricanes
aovHurr<-aov(scale(Hurricanes)~Pop, data=lassodata25)
summary(aovHurr)## Df Sum Sq Mean Sq F value Pr(>F)
## Pop 4 10.66 2.6658 2.922 0.0264 *
## Residuals 76 69.34 0.9123
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# P = 0.026
TukeyHSD(aovHurr)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = scale(Hurricanes) ~ Pop, data = lassodata25)
##
## $Pop
## diff lwr upr p adj
## BHT-AM -0.6985064 -1.5389167 0.1419039 0.1490795
## bwAM-AM -0.2491952 -1.3176669 0.8192765 0.9658157
## CH-AM 0.3374667 -1.1886161 1.8635495 0.9718063
## PHT-AM 0.1512353 -1.1000306 1.4025012 0.9971385
## bwAM-BHT 0.4493112 -0.4178421 1.3164644 0.5990529
## CH-BHT 1.0359731 -0.3565971 2.4285433 0.2399309
## PHT-BHT 0.8497417 -0.2346887 1.9341722 0.1947935
## CH-bwAM 0.5866620 -0.9543098 2.1276337 0.8243278
## PHT-bwAM 0.4004306 -0.8689518 1.6698130 0.9028711
## PHT-CH -0.1862314 -1.8591416 1.4866788 0.9979234#nothing
#Overall (combined risk factors)
aovOVER<-aov(scale(OVERALL2)~Pop, data=lassodata25)
summary(aovOVER)## Df Sum Sq Mean Sq F value Pr(>F)
## Pop 4 2.19 0.5472 0.534 0.711
## Residuals 76 77.81 1.0238# P = 0.71
#####Boxplots of exposure to migration risk-factors grouped by Bird Conservation Region population
#png("VermMigrationBoxplots_HORIZONTAL.png", width = (6.85039*2), height = (6.85039+1), units = 'in', res = 300)
par(mfrow=c(2,4), mar=c(1,5,4,0), pty="s")
boxplot(scale(Forest.cover)~factPop, data=lassodata25, col=bplotcols, ylim=c(-3,4), main="", ylab="", xlab="", yaxt='n', xaxt='n', cex.lab=1.5, cex.axis=1.5)
axis(2, at=c(-3,-2,-1,0,1,2,3,4), labels=c("-3","","","0","","","3",""), cex.axis=1.25)
mtext("Scaled exposure", side=2, line=3, cex=1.25)
mtext("Forest and shrub cover", side=3, line=.5, cex=1.25, font=2)
box(cex=1.25)
legend("bottomleft", legend=c("BW PHT", "BW CH", "BW AM"), bty='n', col=c(bplotcols[1:3]), pch=15, pt.cex=2, cex=1.25)
legend("bottomright", legend=c("GW BHT", "GW AM"), bty='n', col=c(bplotcols[4:5]), pch=15, pt.cex=2, cex=1.25)
text(4,3.5, labels="a", cex=1.25, font=2)
text(5,3.5, labels="b", cex=1.25, font=2)
boxplot(scale(Net.forest.change.2000)~factPop, data=lassodata25, col=bplotcols, ylim=c(-3,4), main="", ylab="", xlab="", yaxt='n', xaxt='n', cex.axis=1.5)
axis(2, at=c(-3,-2,-1,0,1,2,3,4), labels=c("-3","","","0","","","3",""), cex.axis=1.25)
mtext("Net change in forest cover 2000-2010", side=3, line=.5, cex=1.25, font=2)
box(cex=1.25)
boxplot(scale(Agriculture)~factPop, data=lassodata25, col=bplotcols, ylim=c(-3,4), main="", ylab="", xlab="", yaxt='n', xaxt='n', cex.axis=1.5)
mtext("Agricultural cover", side=3, line=.5, cex=1.25, font=2)
axis(2, at=c(-3,-2,-1,0,1,2,3,4), labels=c("-3","","","0","","","3",""), cex.axis=1.25)
text(1,3.5, labels="a", cex=1.25, font=2)
text(4,3.5, labels="a", cex=1.25, font=2)
text(5,3.5, labels="b", cex=1.25, font=2)
box(cex=1.25)
boxplot(scale(Human.footprint)~factPop, data=lassodata25, col=bplotcols, ylim=c(-3,4), main="", ylab="", xlab="", yaxt='n', xaxt='n', cex.axis=1.5)
axis(2, at=c(-3,-2,-1,0,1,2,3,4), labels=c("-3","","","0","","","3",""), cex.axis=1.25)
mtext("Human footprint", side=3, line=.5, cex=1.25, font=2)
box(cex=1.25)
boxplot(scale(Turbines)~factPop, data=lassodata25, col=bplotcols, ylim=c(-3,4), main="", ylab="", xlab="", yaxt='n', xaxt='n', cex.axis=1.5)
axis(2, at=c(-3,-2,-1,0,1,2,3,4), labels=c("-3","","","0","","","3",""), cex.axis=1.25)
mtext("Wind energy", side=3, line=.5, cex=1.25, font=2)
text(4,3.75, labels="a", cex=1.25, font=2)
text(5,3.75, labels="b", cex=1.25, font=2)
box(cex=1.25)
boxplot(scale(Towers)~factPop, data=lassodata25, col=bplotcols, ylim=c(-3,4), main="", ylab="", xlab="", yaxt='n', xaxt='n', cex.axis=1.5)
axis(2, at=c(-3,-2,-1,0,1,2,3,4), labels=c("-3","","","0","","","3",""), cex.axis=1.25)
mtext("Communications towers", side=3, line=.5, cex=1.25, font=2)
box(cex=1.25)
boxplot(scale(Tornados)~factPop, data=lassodata25, col=bplotcols, ylim=c(-3,4), main="", ylab="", xlab="", yaxt='n', xaxt='n', cex.axis=1.5)
axis(2, at=c(-3,-2,-1,0,1,2,3,4), labels=c("-3","","","0","","","3",""), cex.axis=1.25)
mtext("Tornados", side=3, line=.5, cex=1.25, font=2)
text(1,3.5, labels="a", cex=1.25, font=2)
text(2,3.5, labels="a", cex=1.25, font=2)
text(5,3.5, labels="b", cex=1.25, font=2)
box(cex=1.25)
boxplot(scale(Hurricanes)~factPop, data=lassodata25, col=bplotcols, ylim=c(-3,4), main="", ylab="", xlab="", yaxt='n', xaxt='n', cex.axis=1.5)
axis(2, at=c(-3,-2,-1,0,1,2,3,4), labels=c("-3","","","0","","","3",""), cex.axis=1.25)
mtext("Hurricanes", side=3, line=.5, cex=1.25, font=2)
box(cex=1.25)
#dev.off()
#####Power analysis: ANOVAs
#ANOVA
aggregate(scale(lassodata25$Hurricanes)~lassodata25$Pop, FUN=mean) #average overall exposure by pop## lassodata25$Pop V1
## 1 AM 0.3952422
## 2 BHT -0.3032642
## 3 bwAM 0.1460470
## 4 CH 0.7327089
## 5 PHT 0.5464775aggregate(scale(lassodata25$Hurricanes)~lassodata25$Pop, FUN=var) #variance in overall exposure by pop## lassodata25$Pop V1
## 1 AM 0.7755378
## 2 BHT 0.8274563
## 3 bwAM 0.7974873
## 4 CH 3.1021605
## 5 PHT 0.9238888aggregate(scale(lassodata25$Hurricanes)~lassodata25$Pop, FUN=range) #range in overall exposure by pop## lassodata25$Pop V1.1 V1.2
## 1 AM -1.5669795 2.0680758
## 2 BHT -1.8447676 2.3220535
## 3 bwAM -1.0034665 2.1077599
## 4 CH -0.8209201 3.0760497
## 5 PHT -0.9002881 1.7902878aggregate(scale(lassodata25$Hurricanes)~lassodata25$Pop, FUN=length) #length in overall exposure by pop (n)## lassodata25$Pop V1
## 1 AM 13
## 2 BHT 45
## 3 bwAM 12
## 4 CH 4
## 5 PHT 7range(subset(scale(lassodata25$Hurricanes), lassodata25$Pop=="CH"))## [1] -0.8209201 3.0760497mean(c(0.7755378 , 0.8274563, 0.7974873, 0.9238888)) #average w/in group variance w/o CH outlier## [1] 0.8310926a<-aggregate(scale(lassodata25$Hurricanes)~lassodata25$Pop, FUN=mean) #average overall exposure by pop
groupmeans <- a$V1
p <- power.anova.test(groups = length(groupmeans[-5]),
between.var = var(groupmeans[-5]), within.var = 0.8310926,
power=0.9,sig.level=0.05, n=NULL)
p##
## Balanced one-way analysis of variance power calculation
##
## groups = 4
## n = 21.63869
## between.var = 0.1902655
## within.var = 0.8310926
## sig.level = 0.05
## power = 0.9
##
## NOTE: n is number in each groupn <- c(seq(2,10,by=1),seq(12,20,by=2),seq(25,50,by=5))
p <- power.anova.test(groups = length(groupmeans[-5]),
between.var = var(groupmeans[-5]), within.var = 0.8310926,
power=NULL, sig.level=0.05, n=n)
p##
## Balanced one-way analysis of variance power calculation
##
## groups = 4
## n = 2, 3, 4, 5, 6, 7, 8, 9, 10, 12, 14, 16, 18, 20, 25, 30, 35, 40, 45, 50
## between.var = 0.1902655
## within.var = 0.8310926
## sig.level = 0.05
## power = 0.08872084, 0.13906445, 0.19359421, 0.25058764, 0.30864246, 0.36655229, 0.42331245, 0.47811881, 0.53035668, 0.62550695, 0.70689953, 0.77438708, 0.82890990, 0.87199529, 0.94112547, 0.97454708, 0.98954390, 0.99588548, 0.99843952, 0.99942684
##
## NOTE: n is number in each group#png("VermMigrationANOVA_power.png", width = (6.85039), height = (6.85039), units = 'in', res = 300)
plot(p$power~n, ylab="Power", xlab="n")
#dev.off()Linear modeling and hierarchical GLMs
#Generalized linear modeling and hierarchical GLMs
#25th percentile data; univariate models (individual migration risk-factors) vs. intercept only
lassodata25<-read.csv("/Users/gunnarkramer/Desktop/FINAL DATA FOR DRUM LAND ECOL/Data4Lasso25_data.csv", header=TRUE)
#housekeeping
lassodata25$Individual<-as.character(lassodata25$Individual)
lassodata25$Species<-as.character(lassodata25$Species)
lassodata25$Pop<-as.character(lassodata25$Pop)
lassodata25$Trend00<-as.numeric(as.character(lassodata25$Trend00))
lassodata25$TrendGrp<-ifelse(lassodata25$Trend00 <=-1.5, "Declining", ifelse(lassodata25$Trend00>-1.5 & lassodata25$Trend00<=1.5, "Stable", ifelse(lassodata25$Trend00>1.5, "Increasing", "NA")))
#omit hybrids
lassodata25<-subset(lassodata25, Species!="HYBRID")
lassodata25$Species<-ifelse(lassodata25$Species=="BW","Blue-winged warbler", "Golden-winged warbler")
lassodata25$spPop<-paste(lassodata25$Species, lassodata25$Pop, sep="")
#start building models
#Intercept only
lm25a<-lm(Trend00~1, data=lassodata25)
summary(lm25a)##
## Call:
## lm(formula = Trend00 ~ 1, data = lassodata25)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.6717 -0.9217 0.5183 1.7783 20.1383
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.03827 0.61528 -0.062 0.951
##
## Residual standard error: 5.538 on 80 degrees of freedom#Turbines
lm25b<-lm(Trend00~Turbines, data=lassodata25)
summary(lm25b)##
## Call:
## lm(formula = Trend00 ~ Turbines, data = lassodata25)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.8033 -1.6432 0.4175 1.7929 21.2336
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.1474526 0.9450133 -1.214 0.228
## Turbines 0.0003860 0.0002511 1.537 0.128
##
## Residual standard error: 5.491 on 79 degrees of freedom
## Multiple R-squared: 0.02903, Adjusted R-squared: 0.01674
## F-statistic: 2.362 on 1 and 79 DF, p-value: 0.1283#Towers
lm25c<-lm(Trend00~Towers, data=lassodata25)
summary(lm25c)##
## Call:
## lm(formula = Trend00 ~ Towers, data = lassodata25)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.7941 -1.4838 0.5847 2.1370 19.9095
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.214512 1.366483 -0.889 0.377
## Towers 0.004446 0.004611 0.964 0.338
##
## Residual standard error: 5.54 on 79 degrees of freedom
## Multiple R-squared: 0.01163, Adjusted R-squared: -0.0008811
## F-statistic: 0.9296 on 1 and 79 DF, p-value: 0.3379#Hurricanes
lm25d<-lm(Trend00~Hurricanes, data=lassodata25)
summary(lm25d)##
## Call:
## lm(formula = Trend00 ~ Hurricanes, data = lassodata25)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.0163 -1.2123 0.5571 1.7703 20.4111
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.936108 1.595903 0.587 0.559
## Hurricanes -0.003265 0.004931 -0.662 0.510
##
## Residual standard error: 5.557 on 79 degrees of freedom
## Multiple R-squared: 0.005519, Adjusted R-squared: -0.00707
## F-statistic: 0.4384 on 1 and 79 DF, p-value: 0.5098#Tornados
lm25e<-lm(Trend00~Tornados, data=lassodata25)
summary(lm25e)##
## Call:
## lm(formula = Trend00 ~ Tornados, data = lassodata25)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.2388 -1.5545 0.4613 1.9316 20.7110
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.3791654 1.4715731 -1.617 0.1099
## Tornados 0.0006135 0.0003513 1.747 0.0846 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.468 on 79 degrees of freedom
## Multiple R-squared: 0.03718, Adjusted R-squared: 0.02499
## F-statistic: 3.05 on 1 and 79 DF, p-value: 0.0846#Agriculture
lm25f<-lm(Trend00~scale(Agriculture), data=lassodata25)
summary(lm25f)##
## Call:
## lm(formula = Trend00 ~ scale(Agriculture), data = lassodata25)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.1271 -1.6204 0.2306 1.8762 21.2834
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.03827 0.60203 -0.064 0.9495
## scale(Agriculture) 1.29383 0.60578 2.136 0.0358 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.418 on 79 degrees of freedom
## Multiple R-squared: 0.05459, Adjusted R-squared: 0.04262
## F-statistic: 4.562 on 1 and 79 DF, p-value: 0.03579#Forest and shrub cover
lm25g<-lm(Trend00~Forest.cover, data=lassodata25)
summary(lm25g)##
## Call:
## lm(formula = Trend00 ~ Forest.cover, data = lassodata25)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.3509 -1.4553 0.1497 1.6063 20.6162
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.943899 1.867603 1.576 0.1190
## Forest.cover -0.014712 0.008711 -1.689 0.0952 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.475 on 79 degrees of freedom
## Multiple R-squared: 0.03485, Adjusted R-squared: 0.02263
## F-statistic: 2.852 on 1 and 79 DF, p-value: 0.09519#Human footprint
lm25h<-lm(Trend00~Human.footprint, data=lassodata25)
summary(lm25h)##
## Call:
## lm(formula = Trend00 ~ Human.footprint, data = lassodata25)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.4580 -1.0292 0.5777 1.6469 20.3999
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.3616907 2.0117999 0.677 0.500
## Human.footprint -0.0003185 0.0004356 -0.731 0.467
##
## Residual standard error: 5.554 on 79 degrees of freedom
## Multiple R-squared: 0.006721, Adjusted R-squared: -0.005852
## F-statistic: 0.5345 on 1 and 79 DF, p-value: 0.4669#Net forest change 2000-2015
lm25i<-lm(Trend00~Net.forest.change.2000, data=lassodata25)
summary(lm25i)##
## Call:
## lm(formula = Trend00 ~ Net.forest.change.2000, data = lassodata25)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.4584 -1.9874 0.6235 2.0845 20.5826
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.394 1.067 -1.307 0.195
## Net.forest.change.2000 26.578 17.152 1.550 0.125
##
## Residual standard error: 5.49 on 79 degrees of freedom
## Multiple R-squared: 0.0295, Adjusted R-squared: 0.01721
## F-statistic: 2.401 on 1 and 79 DF, p-value: 0.1252#make model table
cand.set<-list()
cand.set[[1]]<-lm25a
cand.set[[2]]<-lm25b
cand.set[[3]]<-lm25c
cand.set[[4]]<-lm25d
cand.set[[5]]<-lm25e
cand.set[[6]]<-lm25f
cand.set[[7]]<-lm25g
cand.set[[8]]<-lm25h
cand.set[[9]]<-lm25i
AICcmodavg::aictab(cand.set, second.ord=TRUE, sort=TRUE)##
## Model selection based on AICc:
##
## K AICc Delta_AICc AICcWt Cum.Wt LL
## Mod6 3 507.90 0.00 0.30 0.30 -250.79
## Mod5 3 509.38 1.48 0.14 0.44 -251.53
## Mod7 3 509.57 1.67 0.13 0.57 -251.63
## Mod9 3 510.02 2.12 0.10 0.68 -251.85
## Mod2 3 510.06 2.16 0.10 0.78 -251.87
## Mod1 2 510.29 2.39 0.09 0.87 -253.07
## Mod3 3 511.50 3.60 0.05 0.92 -252.59
## Mod8 3 511.90 4.00 0.04 0.96 -252.79
## Mod4 3 512.00 4.10 0.04 1.00 -252.84#calculate -2*logLik and look at table
migmodtable<-AICcmodavg::aictab(cand.set, second.ord=TRUE, sort=TRUE)
migmodtable##
## Model selection based on AICc:
##
## K AICc Delta_AICc AICcWt Cum.Wt LL
## Mod6 3 507.90 0.00 0.30 0.30 -250.79
## Mod5 3 509.38 1.48 0.14 0.44 -251.53
## Mod7 3 509.57 1.67 0.13 0.57 -251.63
## Mod9 3 510.02 2.12 0.10 0.68 -251.85
## Mod2 3 510.06 2.16 0.10 0.78 -251.87
## Mod1 2 510.29 2.39 0.09 0.87 -253.07
## Mod3 3 511.50 3.60 0.05 0.92 -252.59
## Mod8 3 511.90 4.00 0.04 0.96 -252.79
## Mod4 3 512.00 4.10 0.04 1.00 -252.84migmodtable$LL*-2## [1] 501.5865 503.0649 503.2607 503.7084 503.7471 506.1336 505.1861 505.5874
## [9] 505.685410th percentile linear modeling
#Generalized linear modeling and hierarchical GLMs
#10th percentile data; univariate models (individual migration risk-factors) vs. intercept only
lassodata10<-read.csv("/Users/gunnarkramer/Desktop/FINAL DATA FOR DRUM LAND ECOL/Data4Lasso10_data.csv", header=TRUE)
#housekeeping
lassodata10$Individual<-as.character(lassodata10$Individual)
lassodata10$Species<-as.character(lassodata10$Species)
lassodata10$Pop<-as.character(lassodata10$Pop)
lassodata10$Trend00<-as.numeric(as.character(lassodata10$Trend00))
lassodata10$TrendGrp<-ifelse(lassodata10$Trend00 <=-1.5, "Declining", ifelse(lassodata10$Trend00>-1.5 & lassodata10$Trend00<=1.5, "Stable", ifelse(lassodata10$Trend00>1.5, "Increasing", "NA")))
#omit hybrids
lassodata10<-subset(lassodata10, Species!="HYBRID")
lassodata10$Species<-ifelse(lassodata10$Species=="BW","Blue-winged warbler", "Golden-winged warbler")
lassodata10$spPop<-paste(lassodata10$Species, lassodata10$Pop, sep="")
colnames(lassodata10)## [1] "Individual" "Species" "Pop"
## [4] "State" "BrLat" "BrLon"
## [7] "NBLat" "NBLon" "BCR_Trend"
## [10] "Trend00" "scale" "Turbines"
## [13] "Towers" "Hurricanes" "Tornados"
## [16] "Agriculture" "Forest.cover" "Human.footprint"
## [19] "Net.forest.change.2000" "OVERALL2" "Br_Forest"
## [22] "Br_HFootP" "Br_forestchange" "Nb_Forest"
## [25] "Nb_HFootP" "Nb_forestchange" "TrendGrp"
## [28] "spPop"#start building models
#Intercept only
lm10a<-lm(Trend00~1, data=lassodata10)
summary(lm10a)##
## Call:
## lm(formula = Trend00 ~ 1, data = lassodata10)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.6717 -0.9217 0.5183 1.7783 20.1383
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.03827 0.61528 -0.062 0.951
##
## Residual standard error: 5.538 on 80 degrees of freedom#Turbines
lm10b<-lm(Trend00~Turbines, data=lassodata10)
summary(lm10b)##
## Call:
## lm(formula = Trend00 ~ Turbines, data = lassodata10)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.549 -1.346 0.213 1.872 21.043
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.9430199 0.8029674 -1.174 0.2438
## Turbines 0.0010506 0.0006093 1.724 0.0886 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.471 on 79 degrees of freedom
## Multiple R-squared: 0.03627, Adjusted R-squared: 0.02407
## F-statistic: 2.973 on 1 and 79 DF, p-value: 0.08855#Towers
lm10c<-lm(Trend00~Towers, data=lassodata10)
summary(lm10c)##
## Call:
## lm(formula = Trend00 ~ Towers, data = lassodata10)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.0708 -1.0483 0.3318 1.6608 20.2175
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.360785 0.900732 0.401 0.690
## Towers -0.005314 0.008729 -0.609 0.544
##
## Residual standard error: 5.559 on 79 degrees of freedom
## Multiple R-squared: 0.004669, Adjusted R-squared: -0.00793
## F-statistic: 0.3706 on 1 and 79 DF, p-value: 0.5444#Hurricanes
lm10d<-lm(Trend00~Hurricanes, data=lassodata10)
summary(lm10d)##
## Call:
## lm(formula = Trend00 ~ Hurricanes, data = lassodata10)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.7740 -0.9555 0.4307 1.7225 20.2341
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.249962 1.411844 0.177 0.860
## Hurricanes -0.002447 0.010771 -0.227 0.821
##
## Residual standard error: 5.571 on 79 degrees of freedom
## Multiple R-squared: 0.0006527, Adjusted R-squared: -0.012
## F-statistic: 0.0516 on 1 and 79 DF, p-value: 0.8209#Tornados
lm10e<-lm(Trend00~Tornados, data=lassodata10)
summary(lm10e)##
## Call:
## lm(formula = Trend00 ~ Tornados, data = lassodata10)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.3096 -1.3806 0.5707 1.7983 20.8479
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.9008776 1.1221665 -0.803 0.424
## Tornados 0.0006624 0.0007204 0.920 0.361
##
## Residual standard error: 5.543 on 79 degrees of freedom
## Multiple R-squared: 0.01059, Adjusted R-squared: -0.001934
## F-statistic: 0.8456 on 1 and 79 DF, p-value: 0.3606#Agriculture
lm10f<-lm(Trend00~scale(Agriculture), data=lassodata10)
summary(lm10f)##
## Call:
## lm(formula = Trend00 ~ scale(Agriculture), data = lassodata10)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.996 -1.674 0.436 2.099 21.118
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.03827 0.60467 -0.063 0.9497
## scale(Agriculture) 1.19107 0.60844 1.958 0.0538 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.442 on 79 degrees of freedom
## Multiple R-squared: 0.04626, Adjusted R-squared: 0.03419
## F-statistic: 3.832 on 1 and 79 DF, p-value: 0.05381#Forest and shrub cover
lm10g<-lm(Trend00~Forest.cover, data=lassodata10)
summary(lm10g)##
## Call:
## lm(formula = Trend00 ~ Forest.cover, data = lassodata10)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.3082 -1.0453 0.8157 1.6404 20.5811
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.99445 1.79169 0.555 0.580
## Forest.cover -0.01435 0.02337 -0.614 0.541
##
## Residual standard error: 5.559 on 79 degrees of freedom
## Multiple R-squared: 0.00475, Adjusted R-squared: -0.007848
## F-statistic: 0.377 on 1 and 79 DF, p-value: 0.541#Human footprint
lm10h<-lm(Trend00~Human.footprint, data=lassodata10)
summary(lm10h)##
## Call:
## lm(formula = Trend00 ~ Human.footprint, data = lassodata10)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.209 -1.048 0.566 1.779 20.265
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.8653415 1.8197010 0.476 0.636
## Human.footprint -0.0006018 0.0011399 -0.528 0.599
##
## Residual standard error: 5.563 on 79 degrees of freedom
## Multiple R-squared: 0.003516, Adjusted R-squared: -0.009098
## F-statistic: 0.2787 on 1 and 79 DF, p-value: 0.599#Net forest change 2000-2015
lm10i<-lm(Trend00~Net.forest.change.2000, data=lassodata10)
summary(lm10i)##
## Call:
## lm(formula = Trend00 ~ Net.forest.change.2000, data = lassodata10)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.8980 -1.0588 0.6265 1.7918 20.4334
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.2612 0.7607 -0.343 0.732
## Net.forest.change.2000 4.1875 8.3275 0.503 0.616
##
## Residual standard error: 5.564 on 79 degrees of freedom
## Multiple R-squared: 0.003191, Adjusted R-squared: -0.009427
## F-statistic: 0.2529 on 1 and 79 DF, p-value: 0.6165#make model table
cand.set<-list()
cand.set[[1]]<-lm10a
cand.set[[2]]<-lm10b
cand.set[[3]]<-lm10c
cand.set[[4]]<-lm10d
cand.set[[5]]<-lm10e
cand.set[[6]]<-lm10f
cand.set[[7]]<-lm10g
cand.set[[8]]<-lm10h
cand.set[[9]]<-lm10i
#look at model table
AICcmodavg::aictab(cand.set, second.ord=TRUE, sort=TRUE)##
## Model selection based on AICc:
##
## K AICc Delta_AICc AICcWt Cum.Wt LL
## Mod6 3 508.61 0.00 0.32 0.32 -251.15
## Mod2 3 509.45 0.84 0.21 0.52 -251.57
## Mod1 2 510.29 1.68 0.14 0.66 -253.07
## Mod5 3 511.58 2.97 0.07 0.73 -252.64
## Mod7 3 512.06 3.45 0.06 0.79 -252.87
## Mod3 3 512.07 3.46 0.06 0.85 -252.88
## Mod8 3 512.16 3.55 0.05 0.90 -252.92
## Mod9 3 512.19 3.58 0.05 0.95 -252.94
## Mod4 3 512.39 3.78 0.05 1.00 -253.04#calculate -2*logLik
migmodtable<-AICcmodavg::aictab(cand.set, second.ord=TRUE, sort=TRUE)
migmodtable##
## Model selection based on AICc:
##
## K AICc Delta_AICc AICcWt Cum.Wt LL
## Mod6 3 508.61 0.00 0.32 0.32 -251.15
## Mod2 3 509.45 0.84 0.21 0.52 -251.57
## Mod1 2 510.29 1.68 0.14 0.66 -253.07
## Mod5 3 511.58 2.97 0.07 0.73 -252.64
## Mod7 3 512.06 3.45 0.06 0.79 -252.87
## Mod3 3 512.07 3.46 0.06 0.85 -252.88
## Mod8 3 512.16 3.55 0.05 0.90 -252.92
## Mod9 3 512.19 3.58 0.05 0.95 -252.94
## Mod4 3 512.39 3.78 0.05 1.00 -253.04migmodtable$LL*-2## [1] 502.2968 503.1409 506.1336 505.2712 505.7480 505.7546 505.8483 505.8748
## [9] 506.080750th percentile linear modeling
#Generalized linear modeling and hierarchical GLMs
#50th percentile data; univariate models (individual migration risk-factors) vs. intercept only
lassodata50<-read.csv("/Users/gunnarkramer/Desktop/FINAL DATA FOR DRUM LAND ECOL/Data4Lasso50_data.csv", header=TRUE)
#housekeeping
lassodata50$Individual<-as.character(lassodata50$Individual)
lassodata50$Species<-as.character(lassodata50$Species)
lassodata50$Pop<-as.character(lassodata50$Pop)
lassodata50$Trend00<-as.numeric(as.character(lassodata50$Trend00))
lassodata50$TrendGrp<-ifelse(lassodata50$Trend00 <=-1.5, "Declining", ifelse(lassodata50$Trend00>-1.5 & lassodata50$Trend00<=1.5, "Stable", ifelse(lassodata50$Trend00>1.5, "Increasing", "NA")))
#omit hybrids
lassodata50<-subset(lassodata50, Species!="HYBRID")
lassodata50$Species<-ifelse(lassodata50$Species=="BW","Blue-winged warbler", "Golden-winged warbler")
lassodata50$spPop<-paste(lassodata50$Species, lassodata50$Pop, sep="")
colnames(lassodata50)## [1] "Individual" "Species" "Pop"
## [4] "State" "BrLat" "BrLon"
## [7] "NBLat" "NBLon" "BCR_Trend"
## [10] "Trend00" "scale" "Turbines"
## [13] "Towers" "Hurricanes" "Tornados"
## [16] "Agriculture" "Forest.cover" "Human.footprint"
## [19] "Net.forest.change.2000" "OVERALL2" "Br_Forest"
## [22] "Br_HFootP" "Br_forestchange" "Nb_Forest"
## [25] "Nb_HFootP" "Nb_forestchange" "TrendGrp"
## [28] "spPop"#start building models
#Intercept only
lm50a<-lm(Trend00~1, data=lassodata50)
summary(lm50a)##
## Call:
## lm(formula = Trend00 ~ 1, data = lassodata50)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.6717 -0.9217 0.5183 1.7783 20.1383
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.03827 0.61528 -0.062 0.951
##
## Residual standard error: 5.538 on 80 degrees of freedom#Turbines
lm50b<-lm(Trend00~Turbines, data=lassodata50)
summary(lm50b)##
## Call:
## lm(formula = Trend00 ~ Turbines, data = lassodata50)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.3869 -1.8940 0.3188 1.8906 20.1971
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.0799063 1.2117246 -2.542 0.01299 *
## Turbines 0.0003567 0.0001242 2.873 0.00523 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.302 on 79 degrees of freedom
## Multiple R-squared: 0.09457, Adjusted R-squared: 0.08311
## F-statistic: 8.252 on 1 and 79 DF, p-value: 0.005226#Towers
lm50c<-lm(Trend00~Towers, data=lassodata50)
summary(lm50c)##
## Call:
## lm(formula = Trend00 ~ Towers, data = lassodata50)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.0214 -0.9983 0.5351 1.8483 20.2210
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.872946 2.024167 -0.431 0.667
## Towers 0.001268 0.002928 0.433 0.666
##
## Residual standard error: 5.566 on 79 degrees of freedom
## Multiple R-squared: 0.002368, Adjusted R-squared: -0.01026
## F-statistic: 0.1875 on 1 and 79 DF, p-value: 0.6662#Hurricanes
lm50d<-lm(Trend00~Hurricanes, data=lassodata50)
summary(lm50d)##
## Call:
## lm(formula = Trend00 ~ Hurricanes, data = lassodata50)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.9687 -1.7364 0.2195 1.8088 20.5215
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.735262 2.006541 1.363 0.177
## Hurricanes -0.004301 0.002964 -1.451 0.151
##
## Residual standard error: 5.5 on 79 degrees of freedom
## Multiple R-squared: 0.02597, Adjusted R-squared: 0.01364
## F-statistic: 2.106 on 1 and 79 DF, p-value: 0.1507#Tornados
lm50e<-lm(Trend00~Tornados, data=lassodata50)
summary(lm50e)##
## Call:
## lm(formula = Trend00 ~ Tornados, data = lassodata50)
##
## Residuals:
## Min 1Q Median 3Q Max
## -11.2229 -1.9711 0.4001 2.0041 20.5298
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -5.3209106 2.0072717 -2.651 0.0097 **
## Tornados 0.0005775 0.0002097 2.754 0.0073 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.323 on 79 degrees of freedom
## Multiple R-squared: 0.0876, Adjusted R-squared: 0.07605
## F-statistic: 7.585 on 1 and 79 DF, p-value: 0.007303#Agriculture
lm50f<-lm(Trend00~scale(Agriculture), data=lassodata50)
summary(lm50f)##
## Call:
## lm(formula = Trend00 ~ scale(Agriculture), data = lassodata50)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.8502 -1.5617 0.4393 1.7502 20.0542
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.03827 0.59954 -0.064 0.9493
## scale(Agriculture) 1.38311 0.60328 2.293 0.0245 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.396 on 79 degrees of freedom
## Multiple R-squared: 0.06238, Adjusted R-squared: 0.05052
## F-statistic: 5.256 on 1 and 79 DF, p-value: 0.02453#Forest and shrub cover
lm50g<-lm(Trend00~Forest.cover, data=lassodata50)
summary(lm50g)##
## Call:
## lm(formula = Trend00 ~ Forest.cover, data = lassodata50)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.7777 -2.4562 -0.0409 1.6669 20.0505
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.098164 1.980560 2.574 0.01192 *
## Forest.cover -0.010201 0.003754 -2.718 0.00808 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.329 on 79 degrees of freedom
## Multiple R-squared: 0.0855, Adjusted R-squared: 0.07392
## F-statistic: 7.386 on 1 and 79 DF, p-value: 0.008076#Human footprint
lm50h<-lm(Trend00~Human.footprint, data=lassodata50)
summary(lm50h)##
## Call:
## lm(formula = Trend00 ~ Human.footprint, data = lassodata50)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.4273 -1.7554 0.3204 1.9070 20.7445
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.8059732 2.3194538 1.641 0.1048
## Human.footprint -0.0003597 0.0002094 -1.717 0.0898 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.471 on 79 degrees of freedom
## Multiple R-squared: 0.03599, Adjusted R-squared: 0.02379
## F-statistic: 2.95 on 1 and 79 DF, p-value: 0.08982#Net forest change 2000-2015
lm50i<-lm(Trend00~Net.forest.change.2000, data=lassodata50)
summary(lm50i)##
## Call:
## lm(formula = Trend00 ~ Net.forest.change.2000, data = lassodata50)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.4049 -1.6035 0.4859 2.1161 20.1618
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.975 1.922 -1.548 0.126
## Net.forest.change.2000 58.824 36.505 1.611 0.111
##
## Residual standard error: 5.483 on 79 degrees of freedom
## Multiple R-squared: 0.03182, Adjusted R-squared: 0.01957
## F-statistic: 2.597 on 1 and 79 DF, p-value: 0.1111#make model table
cand.set<-list()
cand.set[[1]]<-lm50a
cand.set[[2]]<-lm50b
cand.set[[3]]<-lm50c
cand.set[[4]]<-lm50d
cand.set[[5]]<-lm50e
cand.set[[6]]<-lm50f
cand.set[[7]]<-lm50g
cand.set[[8]]<-lm50h
cand.set[[9]]<-lm50i
#look at model table
AICcmodavg::aictab(cand.set, second.ord=TRUE, sort=TRUE)##
## Model selection based on AICc:
##
## K AICc Delta_AICc AICcWt Cum.Wt LL
## Mod2 3 504.40 0.00 0.34 0.34 -249.04
## Mod5 3 505.02 0.62 0.25 0.59 -249.35
## Mod7 3 505.21 0.81 0.23 0.82 -249.45
## Mod6 3 507.23 2.83 0.08 0.91 -250.46
## Mod8 3 509.48 5.08 0.03 0.93 -251.58
## Mod9 3 509.83 5.43 0.02 0.96 -251.76
## Mod1 2 510.29 5.89 0.02 0.98 -253.07
## Mod4 3 510.31 5.92 0.02 0.99 -252.00
## Mod3 3 512.25 7.86 0.01 1.00 -252.97#calculate -2*logLik
migmodtable<-AICcmodavg::aictab(cand.set, second.ord=TRUE, sort=TRUE)
migmodtable##
## Model selection based on AICc:
##
## K AICc Delta_AICc AICcWt Cum.Wt LL
## Mod2 3 504.40 0.00 0.34 0.34 -249.04
## Mod5 3 505.02 0.62 0.25 0.59 -249.35
## Mod7 3 505.21 0.81 0.23 0.82 -249.45
## Mod6 3 507.23 2.83 0.08 0.91 -250.46
## Mod8 3 509.48 5.08 0.03 0.93 -251.58
## Mod9 3 509.83 5.43 0.02 0.96 -251.76
## Mod1 2 510.29 5.89 0.02 0.98 -253.07
## Mod4 3 510.31 5.92 0.02 0.99 -252.00
## Mod3 3 512.25 7.86 0.01 1.00 -252.97migmodtable$LL*-2## [1] 498.0863 498.7081 498.8940 500.9160 503.1645 503.5141 506.1336 504.0027
## [9] 505.9416Multiple linear regression (hierarchical GLMs); 25th percentile; annual exposure
########
#intercept only (not used)
hm1.null<-glm(Trend00 ~ 1, family = "gaussian", data=lassodata25)
summary(hm1.null)##
## Call:
## glm(formula = Trend00 ~ 1, family = "gaussian", data = lassodata25)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -8.6717 -0.9217 0.5183 1.7783 20.1383
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.03827 0.61528 -0.062 0.951
##
## (Dispersion parameter for gaussian family taken to be 30.6646)
##
## Null deviance: 2453.2 on 80 degrees of freedom
## Residual deviance: 2453.2 on 80 degrees of freedom
## AIC: 510.13
##
## Number of Fisher Scoring iterations: 2#Base model incorporating proxies for breeding/nonbreeding factors (breeding and nonbreeding latitude/longitude)
hm1.mc<-glm(Trend00 ~ BrLon + BrLat + NonBrLon + NonBrLat, family = "gaussian", data=lassodata25)
summary(hm1.mc)##
## Call:
## glm(formula = Trend00 ~ BrLon + BrLat + NonBrLon + NonBrLat,
## family = "gaussian", data = lassodata25)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -4.8747 -1.7047 -1.0403 0.6687 13.8656
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -53.74041 6.04236 -8.894 2.13e-13 ***
## BrLon -0.07836 0.07220 -1.085 0.281202
## BrLat 1.00650 0.14939 6.738 2.73e-09 ***
## NonBrLon 0.09334 0.12122 0.770 0.443676
## NonBrLat 0.68794 0.19232 3.577 0.000609 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 12.17733)
##
## Null deviance: 2453.17 on 80 degrees of freedom
## Residual deviance: 925.48 on 76 degrees of freedom
## AIC: 439.17
##
## Number of Fisher Scoring iterations: 2###Full model, breeding nonbreeding lat/lon + all migration risk-factors
hm1.mc.full<-glm(Trend00 ~ BrLon + BrLat + NonBrLon + NonBrLat +
scale(Turbines) + scale(Towers) +
scale(Hurricanes) + scale(Tornados) +
scale(Agriculture) + scale(Forest.cover) +
scale(Human.footprint) + scale(Net.forest.change.2000),
family = "gaussian", data=lassodata25)
summary(hm1.mc.full)##
## Call:
## glm(formula = Trend00 ~ BrLon + BrLat + NonBrLon + NonBrLat +
## scale(Turbines) + scale(Towers) + scale(Hurricanes) + scale(Tornados) +
## scale(Agriculture) + scale(Forest.cover) + scale(Human.footprint) +
## scale(Net.forest.change.2000), family = "gaussian", data = lassodata25)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -4.4355 -2.1425 -0.2716 0.8478 11.3765
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -57.94263 9.56922 -6.055 6.81e-08 ***
## BrLon -0.11329 0.08782 -1.290 0.201419
## BrLat 1.05863 0.15992 6.620 6.82e-09 ***
## NonBrLon 0.11795 0.14367 0.821 0.414511
## NonBrLat 0.74733 0.20879 3.579 0.000641 ***
## scale(Turbines) 0.82768 0.92166 0.898 0.372339
## scale(Towers) 1.65241 0.61727 2.677 0.009303 **
## scale(Hurricanes) 0.28351 0.68289 0.415 0.679330
## scale(Tornados) -0.89962 0.91122 -0.987 0.327006
## scale(Agriculture) -0.11805 1.19016 -0.099 0.921280
## scale(Forest.cover) 1.27261 1.14419 1.112 0.269954
## scale(Human.footprint) -2.17516 1.58053 -1.376 0.173268
## scale(Net.forest.change.2000) 0.07495 0.59457 0.126 0.900063
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 11.76653)
##
## Null deviance: 2453.17 on 80 degrees of freedom
## Residual deviance: 800.12 on 68 degrees of freedom
## AIC: 443.38
##
## Number of Fisher Scoring iterations: 2###Singly add individual migration risk-factor: TOWERS
hm1.mc.towers<-glm(Trend00 ~ BrLon + BrLat + NonBrLon + NonBrLat +
scale(Towers), family = "gaussian", data=lassodata25)
summary(hm1.mc.towers)##
## Call:
## glm(formula = Trend00 ~ BrLon + BrLat + NonBrLon + NonBrLat +
## scale(Towers), family = "gaussian", data = lassodata25)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -5.2233 -1.8105 -0.6896 0.8484 13.5802
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -54.14571 5.98184 -9.052 1.19e-13 ***
## BrLon -0.08641 0.07159 -1.207 0.231170
## BrLat 1.00067 0.14781 6.770 2.50e-09 ***
## NonBrLon 0.09231 0.11990 0.770 0.443772
## NonBrLat 0.67964 0.19029 3.572 0.000624 ***
## scale(Towers) 0.63393 0.38731 1.637 0.105876
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 11.91414)
##
## Null deviance: 2453.17 on 80 degrees of freedom
## Residual deviance: 893.56 on 75 degrees of freedom
## AIC: 438.33
##
## Number of Fisher Scoring iterations: 2###Singly add individual migration risk-factor: TURBINES
hm1.mc.Turbines<-glm(Trend00 ~ BrLon + BrLat + NonBrLon + NonBrLat +
scale(Turbines), family = "gaussian", data=lassodata25)
summary(hm1.mc.Turbines)##
## Call:
## glm(formula = Trend00 ~ BrLon + BrLat + NonBrLon + NonBrLat +
## scale(Turbines), family = "gaussian", data = lassodata25)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -5.3289 -1.7220 -0.7514 0.5981 13.2997
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -56.45831 6.63312 -8.512 1.27e-12 ***
## BrLon -0.09081 0.07329 -1.239 0.219158
## BrLat 0.99663 0.14973 6.656 4.07e-09 ***
## NonBrLon 0.06656 0.12419 0.536 0.593577
## NonBrLat 0.67669 0.19267 3.512 0.000756 ***
## scale(Turbines) -0.42834 0.43108 -0.994 0.323598
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 12.17936)
##
## Null deviance: 2453.17 on 80 degrees of freedom
## Residual deviance: 913.45 on 75 degrees of freedom
## AIC: 440.11
##
## Number of Fisher Scoring iterations: 2###Singly add individual migration risk-factor: HURRICANES
hm1.mc.Hurricanes<-glm(Trend00 ~ BrLon + BrLat + NonBrLon + NonBrLat +
scale(Hurricanes), family = "gaussian", data=lassodata25)
summary(hm1.mc.Hurricanes)##
## Call:
## glm(formula = Trend00 ~ BrLon + BrLat + NonBrLon + NonBrLat +
## scale(Hurricanes), family = "gaussian", data = lassodata25)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -5.1937 -1.8003 -0.8458 0.6620 13.5862
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -55.20956 6.32191 -8.733 4.79e-13 ***
## BrLon -0.07236 0.07274 -0.995 0.32309
## BrLat 1.01328 0.14996 6.757 2.64e-09 ***
## NonBrLon 0.06191 0.12754 0.485 0.62879
## NonBrLat 0.62920 0.20595 3.055 0.00311 **
## scale(Hurricanes) 0.35411 0.43722 0.810 0.42054
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 12.2327)
##
## Null deviance: 2453.17 on 80 degrees of freedom
## Residual deviance: 917.45 on 75 degrees of freedom
## AIC: 440.47
##
## Number of Fisher Scoring iterations: 2###Singly add individual migration risk-factor: TORNADOS
hm1.mc.Tornados<-glm(Trend00 ~ BrLon + BrLat + NonBrLon + NonBrLat +
scale(Tornados), family = "gaussian", data=lassodata25)
summary(hm1.mc.Tornados)##
## Call:
## glm(formula = Trend00 ~ BrLon + BrLat + NonBrLon + NonBrLat +
## scale(Tornados), family = "gaussian", data = lassodata25)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -4.9188 -1.6451 -0.8415 0.7002 13.6122
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -55.99301 6.62791 -8.448 1.67e-12 ***
## BrLon -0.09694 0.07569 -1.281 0.204207
## BrLat 0.99557 0.15026 6.626 4.63e-09 ***
## NonBrLon 0.08402 0.12197 0.689 0.493044
## NonBrLat 0.70862 0.19429 3.647 0.000487 ***
## scale(Tornados) -0.37083 0.44396 -0.835 0.406223
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 12.22597)
##
## Null deviance: 2453.17 on 80 degrees of freedom
## Residual deviance: 916.95 on 75 degrees of freedom
## AIC: 440.42
##
## Number of Fisher Scoring iterations: 2###Singly add individual migration risk-factor: AGRICULTURE
hm1.mc.Agriculture<-glm(Trend00 ~ BrLon + BrLat + NonBrLon + NonBrLat +
scale(Agriculture), family = "gaussian", data=lassodata25)
summary(hm1.mc.Agriculture)##
## Call:
## glm(formula = Trend00 ~ BrLon + BrLat + NonBrLon + NonBrLat +
## scale(Agriculture), family = "gaussian", data = lassodata25)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -5.1900 -1.7357 -0.6801 0.5726 13.2818
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -56.54429 6.69099 -8.451 1.65e-12 ***
## BrLon -0.09822 0.07503 -1.309 0.194482
## BrLat 1.01319 0.14959 6.773 2.46e-09 ***
## NonBrLon 0.08760 0.12139 0.722 0.472750
## NonBrLat 0.70681 0.19334 3.656 0.000473 ***
## scale(Agriculture) -0.42834 0.43846 -0.977 0.331739
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 12.18464)
##
## Null deviance: 2453.17 on 80 degrees of freedom
## Residual deviance: 913.85 on 75 degrees of freedom
## AIC: 440.15
##
## Number of Fisher Scoring iterations: 2###Singly add individual migration risk-factor: TOWERS
hm1.mc.Forest.cover<-glm(Trend00 ~ BrLon + BrLat + NonBrLon + NonBrLat +
scale(Forest.cover), family = "gaussian", data=lassodata25)
summary(hm1.mc.Forest.cover)##
## Call:
## glm(formula = Trend00 ~ BrLon + BrLat + NonBrLon + NonBrLat +
## scale(Forest.cover), family = "gaussian", data = lassodata25)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -4.8615 -1.7620 -0.9039 0.6538 13.6429
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -55.99937 6.75782 -8.287 3.40e-12 ***
## BrLon -0.08576 0.07307 -1.174 0.24422
## BrLat 1.00185 0.14994 6.682 3.64e-09 ***
## NonBrLon 0.06619 0.12677 0.522 0.60313
## NonBrLat 0.65919 0.19659 3.353 0.00125 **
## scale(Forest.cover) 0.33180 0.43943 0.755 0.45257
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 12.2466)
##
## Null deviance: 2453.17 on 80 degrees of freedom
## Residual deviance: 918.49 on 75 degrees of freedom
## AIC: 440.56
##
## Number of Fisher Scoring iterations: 2###Singly add individual migration risk-factor: HUMAN FOOTPRINT
hm1.mc.Human.footprint<-glm(Trend00 ~ BrLon + BrLat + NonBrLon + NonBrLat +
scale(Human.footprint), family = "gaussian", data=lassodata25)
summary(hm1.mc.Human.footprint)##
## Call:
## glm(formula = Trend00 ~ BrLon + BrLat + NonBrLon + NonBrLat +
## scale(Human.footprint), family = "gaussian", data = lassodata25)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -4.8982 -1.6512 -1.0360 0.6786 13.8402
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -53.55606 6.23746 -8.586 9.13e-13 ***
## BrLon -0.07741 0.07302 -1.060 0.292474
## BrLat 1.00868 0.15125 6.669 3.85e-09 ***
## NonBrLon 0.09643 0.12420 0.776 0.439941
## NonBrLat 0.69188 0.19582 3.533 0.000707 ***
## scale(Human.footprint) -0.05415 0.40675 -0.133 0.894453
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 12.33678)
##
## Null deviance: 2453.17 on 80 degrees of freedom
## Residual deviance: 925.26 on 75 degrees of freedom
## AIC: 441.15
##
## Number of Fisher Scoring iterations: 2###Singly add individual migration risk-factor: NET FOREST CHANGE
hm1.mc.Net.forest.change.2000<-glm(Trend00 ~ BrLon + BrLat + NonBrLon + NonBrLat +
scale(Net.forest.change.2000), family = "gaussian", data=lassodata25)
summary(hm1.mc.Net.forest.change.2000)##
## Call:
## glm(formula = Trend00 ~ BrLon + BrLat + NonBrLon + NonBrLat +
## scale(Net.forest.change.2000), family = "gaussian", data = lassodata25)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -4.8822 -1.5357 -1.0242 0.6527 13.8220
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -54.08140 6.27099 -8.624 7.73e-13 ***
## BrLon -0.07588 0.07351 -1.032 0.305293
## BrLat 1.00935 0.15087 6.690 3.52e-09 ***
## NonBrLon 0.08745 0.12482 0.701 0.485704
## NonBrLat 0.68431 0.19422 3.523 0.000729 ***
## scale(Net.forest.change.2000) -0.09293 0.41801 -0.222 0.824678
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 12.33157)
##
## Null deviance: 2453.17 on 80 degrees of freedom
## Residual deviance: 924.87 on 75 degrees of freedom
## AIC: 441.12
##
## Number of Fisher Scoring iterations: 2#AIC table
cand.set<-list()
cand.set[[1]]<-hm1.mc
cand.set[[2]]<-hm1.mc.towers
cand.set[[3]]<-hm1.mc.Turbines
cand.set[[4]]<-hm1.mc.Hurricanes
cand.set[[5]]<-hm1.mc.Tornados
cand.set[[6]]<-hm1.mc.Agriculture
cand.set[[7]]<-hm1.mc.Forest.cover
cand.set[[8]]<-hm1.mc.Human.footprint
cand.set[[9]]<-hm1.mc.Net.forest.change.2000
cand.set[[10]]<-hm1.mc.full
#build model table based on AICc
AICcmodavg::aictab(cand.set, second.ord=TRUE, sort=TRUE)##
## Model selection based on AICc:
##
## K AICc Delta_AICc AICcWt Cum.Wt LL
## Mod2 7 439.86 0.00 0.24 0.24 -212.16
## Mod1 6 440.31 0.44 0.19 0.44 -213.59
## Mod3 7 441.65 1.78 0.10 0.53 -213.06
## Mod6 7 441.68 1.82 0.10 0.63 -213.07
## Mod5 7 441.96 2.09 0.08 0.72 -213.21
## Mod4 7 442.00 2.14 0.08 0.80 -213.23
## Mod7 7 442.09 2.23 0.08 0.88 -213.28
## Mod9 7 442.65 2.79 0.06 0.94 -213.56
## Mod8 7 442.69 2.82 0.06 1.00 -213.58
## Mod10 14 449.75 9.88 0.00 1.00 -207.69#look at table and calculate -2*logLik
migmodtable<-AICcmodavg::aictab(cand.set, second.ord=TRUE, sort=TRUE)
migmodtable##
## Model selection based on AICc:
##
## K AICc Delta_AICc AICcWt Cum.Wt LL
## Mod2 7 439.86 0.00 0.24 0.24 -212.16
## Mod1 6 440.31 0.44 0.19 0.44 -213.59
## Mod3 7 441.65 1.78 0.10 0.53 -213.06
## Mod6 7 441.68 1.82 0.10 0.63 -213.07
## Mod5 7 441.96 2.09 0.08 0.72 -213.21
## Mod4 7 442.00 2.14 0.08 0.80 -213.23
## Mod7 7 442.09 2.23 0.08 0.88 -213.28
## Mod9 7 442.65 2.79 0.06 0.94 -213.56
## Mod8 7 442.69 2.82 0.06 1.00 -213.58
## Mod10 14 449.75 9.88 0.00 1.00 -207.69migmodtable$LL*-2## [1] 424.3300 427.1727 426.1134 426.1484 426.4227 426.4673 426.5593 427.1194
## [9] 427.1536 415.3838Multiple linear regression (hierarchical GLMs); 50th percentile; annual exposure
#50 percentile
hm1.null<-glm(Trend00 ~ 1, family = "gaussian", data=lassodata50)
summary(hm1.null)##
## Call:
## glm(formula = Trend00 ~ 1, family = "gaussian", data = lassodata50)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -8.6717 -0.9217 0.5183 1.7783 20.1383
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.03827 0.61528 -0.062 0.951
##
## (Dispersion parameter for gaussian family taken to be 30.6646)
##
## Null deviance: 2453.2 on 80 degrees of freedom
## Residual deviance: 2453.2 on 80 degrees of freedom
## AIC: 510.13
##
## Number of Fisher Scoring iterations: 2#Base model with migratory connectivity terms
hm1.mc<-glm(Trend00 ~ BrLon + BrLat + NBLon + NBLat, family = "gaussian", data=lassodata50)
summary(hm1.mc)##
## Call:
## glm(formula = Trend00 ~ BrLon + BrLat + NBLon + NBLat, family = "gaussian",
## data = lassodata50)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -4.8747 -1.7047 -1.0403 0.6687 13.8656
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -53.74041 6.04236 -8.894 2.13e-13 ***
## BrLon -0.07836 0.07220 -1.085 0.281202
## BrLat 1.00650 0.14939 6.738 2.73e-09 ***
## NBLon 0.09334 0.12122 0.770 0.443676
## NBLat 0.68794 0.19232 3.577 0.000609 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 12.17733)
##
## Null deviance: 2453.17 on 80 degrees of freedom
## Residual deviance: 925.48 on 76 degrees of freedom
## AIC: 439.17
##
## Number of Fisher Scoring iterations: 2###FULL MODEL MC + ALL MIG RISK TERMS
hm1.mc.full<-glm(Trend00 ~ BrLon + BrLat + NBLon + NBLat +
scale(Turbines) + scale(Towers) +
scale(Hurricanes) + scale(Tornados) +
scale(Agriculture) + scale(Forest.cover) +
scale(Human.footprint) + scale(Net.forest.change.2000),
family = "gaussian", data=lassodata50)
summary(hm1.mc.full)##
## Call:
## glm(formula = Trend00 ~ BrLon + BrLat + NBLon + NBLat + scale(Turbines) +
## scale(Towers) + scale(Hurricanes) + scale(Tornados) + scale(Agriculture) +
## scale(Forest.cover) + scale(Human.footprint) + scale(Net.forest.change.2000),
## family = "gaussian", data = lassodata50)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -5.2743 -2.0591 -0.5663 0.4407 13.7443
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -55.50443 12.23686 -4.536 2.40e-05 ***
## BrLon -0.10844 0.09436 -1.149 0.2545
## BrLat 0.99208 0.15934 6.226 3.41e-08 ***
## NBLon 0.10290 0.15570 0.661 0.5109
## NBLat 0.72348 0.22006 3.288 0.0016 **
## scale(Turbines) 0.57952 1.06759 0.543 0.5890
## scale(Towers) 0.73575 0.77733 0.947 0.3472
## scale(Hurricanes) 0.04336 0.86507 0.050 0.9602
## scale(Tornados) 0.41137 1.31607 0.313 0.7556
## scale(Agriculture) -1.43257 1.75873 -0.815 0.4182
## scale(Forest.cover) -0.50850 1.73040 -0.294 0.7698
## scale(Human.footprint) 0.26648 2.27085 0.117 0.9069
## scale(Net.forest.change.2000) -0.35897 0.65839 -0.545 0.5874
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 13.14995)
##
## Null deviance: 2453.2 on 80 degrees of freedom
## Residual deviance: 894.2 on 68 degrees of freedom
## AIC: 452.39
##
## Number of Fisher Scoring iterations: 2###Singly add individual migration risk-factor: TOWERS
hm1.mc.towers<-glm(Trend00 ~ BrLon + BrLat + NBLon + NBLat +
scale(Towers), family = "gaussian", data=lassodata50)
summary(hm1.mc.towers)##
## Call:
## glm(formula = Trend00 ~ BrLon + BrLat + NBLon + NBLat + scale(Towers),
## family = "gaussian", data = lassodata50)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -5.1781 -1.8102 -0.8763 0.6327 13.9655
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -54.09843 6.07673 -8.903 2.28e-13 ***
## BrLon -0.08045 0.07245 -1.110 0.270352
## BrLat 0.99985 0.15004 6.664 3.94e-09 ***
## NBLon 0.08483 0.12205 0.695 0.489175
## NBLat 0.67234 0.19391 3.467 0.000874 ***
## scale(Towers) 0.30266 0.39448 0.767 0.445352
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 12.2436)
##
## Null deviance: 2453.17 on 80 degrees of freedom
## Residual deviance: 918.27 on 75 degrees of freedom
## AIC: 440.54
##
## Number of Fisher Scoring iterations: 2###Singly add individual migration risk-factor: TURBINES
hm1.mc.Turbines<-glm(Trend00 ~ BrLon + BrLat + NBLon + NBLat +
scale(Turbines), family = "gaussian", data=lassodata50)
summary(hm1.mc.Turbines)##
## Call:
## glm(formula = Trend00 ~ BrLon + BrLat + NBLon + NBLat + scale(Turbines),
## family = "gaussian", data = lassodata50)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -4.9036 -1.6109 -1.0617 0.6708 13.8351
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -54.17004 7.19181 -7.532 9.22e-11 ***
## BrLon -0.08108 0.07661 -1.058 0.293322
## BrLat 1.00685 0.15040 6.695 3.45e-09 ***
## NBLon 0.09082 0.12407 0.732 0.466450
## NBLat 0.68583 0.19449 3.526 0.000723 ***
## scale(Turbines) -0.05242 0.46833 -0.112 0.911173
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 12.33763)
##
## Null deviance: 2453.17 on 80 degrees of freedom
## Residual deviance: 925.32 on 75 degrees of freedom
## AIC: 441.16
##
## Number of Fisher Scoring iterations: 2###Singly add individual migration risk-factor: HURRICANES
hm1.mc.Hurricanes<-glm(Trend00 ~ BrLon + BrLat + NBLon + NBLat +
scale(Hurricanes), family = "gaussian", data=lassodata50)
summary(hm1.mc.Hurricanes)##
## Call:
## glm(formula = Trend00 ~ BrLon + BrLat + NBLon + NBLat + scale(Hurricanes),
## family = "gaussian", data = lassodata50)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -5.0038 -1.7371 -0.9294 0.5931 13.7569
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -54.93136 6.54230 -8.396 2.1e-12 ***
## BrLon -0.07716 0.07261 -1.063 0.29134
## BrLat 1.00664 0.15014 6.705 3.3e-09 ***
## NBLon 0.07298 0.12873 0.567 0.57247
## NBLat 0.66193 0.20046 3.302 0.00147 **
## scale(Hurricanes) 0.21242 0.43407 0.489 0.62600
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 12.30042)
##
## Null deviance: 2453.17 on 80 degrees of freedom
## Residual deviance: 922.53 on 75 degrees of freedom
## AIC: 440.91
##
## Number of Fisher Scoring iterations: 2###Singly add individual migration risk-factor: TORNADOS
hm1.mc.Tornados<-glm(Trend00 ~ BrLon + BrLat + NBLon + NBLat +
scale(Tornados), family = "gaussian", data=lassodata50)
summary(hm1.mc.Tornados)##
## Call:
## glm(formula = Trend00 ~ BrLon + BrLat + NBLon + NBLat + scale(Tornados),
## family = "gaussian", data = lassodata50)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -4.8647 -1.6748 -1.0105 0.6444 13.8569
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -53.38903 6.81900 -7.829 2.52e-11 ***
## BrLon -0.07530 0.07748 -0.972 0.334194
## BrLat 1.00668 0.15037 6.694 3.45e-09 ***
## NBLon 0.09391 0.12211 0.769 0.444287
## NBLat 0.68516 0.19511 3.512 0.000758 ***
## scale(Tornados) 0.05156 0.45245 0.114 0.909579
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 12.33756)
##
## Null deviance: 2453.17 on 80 degrees of freedom
## Residual deviance: 925.32 on 75 degrees of freedom
## AIC: 441.16
##
## Number of Fisher Scoring iterations: 2###Singly add individual migration risk-factor: AGRICULTURE
hm1.mc.Agriculture<-glm(Trend00 ~ BrLon + BrLat + NBLon + NBLat +
scale(Agriculture), family = "gaussian", data=lassodata50)
summary(hm1.mc.Agriculture)##
## Call:
## glm(formula = Trend00 ~ BrLon + BrLat + NBLon + NBLat + scale(Agriculture),
## family = "gaussian", data = lassodata50)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -4.9394 -1.6094 -0.9540 0.6416 13.8121
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -54.78973 6.86750 -7.978 1.31e-11 ***
## BrLon -0.08780 0.07811 -1.124 0.264581
## BrLat 1.00453 0.15039 6.679 3.68e-09 ***
## NBLon 0.09009 0.12233 0.736 0.463772
## NBLat 0.69041 0.19360 3.566 0.000635 ***
## scale(Agriculture) -0.14873 0.45308 -0.328 0.743632
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 12.32199)
##
## Null deviance: 2453.17 on 80 degrees of freedom
## Residual deviance: 924.15 on 75 degrees of freedom
## AIC: 441.06
##
## Number of Fisher Scoring iterations: 2###Singly add individual migration risk-factor: Forest and shrub cover
hm1.mc.Forest.cover<-glm(Trend00 ~ BrLon + BrLat + NBLon + NBLat +
scale(Forest.cover), family = "gaussian", data=lassodata50)
summary(hm1.mc.Forest.cover)##
## Call:
## glm(formula = Trend00 ~ BrLon + BrLat + NBLon + NBLat + scale(Forest.cover),
## family = "gaussian", data = lassodata50)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -4.9150 -1.6320 -1.0314 0.7031 13.8560
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -53.27541 7.18160 -7.418 1.51e-10 ***
## BrLon -0.07730 0.07320 -1.056 0.294366
## BrLat 1.00765 0.15066 6.688 3.54e-09 ***
## NBLon 0.09946 0.13198 0.754 0.453420
## NBLat 0.69365 0.19918 3.482 0.000832 ***
## scale(Forest.cover) -0.05718 0.46965 -0.122 0.903419
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 12.33725)
##
## Null deviance: 2453.17 on 80 degrees of freedom
## Residual deviance: 925.29 on 75 degrees of freedom
## AIC: 441.16
##
## Number of Fisher Scoring iterations: 2###Singly add individual migration risk-factor: HUMAN FOOTPRINT
hm1.mc.Human.footprint<-glm(Trend00 ~ BrLon + BrLat + NBLon + NBLat +
scale(Human.footprint), family = "gaussian", data=lassodata50)
summary(hm1.mc.Human.footprint)##
## Call:
## glm(formula = Trend00 ~ BrLon + BrLat + NBLon + NBLat + scale(Human.footprint),
## family = "gaussian", data = lassodata50)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -4.8755 -1.7021 -1.0413 0.6705 13.8651
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -53.72827 6.52830 -8.230 4.35e-12 ***
## BrLon -0.07836 0.07269 -1.078 0.284494
## BrLat 1.00654 0.15055 6.686 3.58e-09 ***
## NBLon 0.09353 0.12784 0.732 0.466653
## NBLat 0.68813 0.19719 3.490 0.000813 ***
## scale(Human.footprint) -0.00219 0.42769 -0.005 0.995929
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 12.33969)
##
## Null deviance: 2453.17 on 80 degrees of freedom
## Residual deviance: 925.48 on 75 degrees of freedom
## AIC: 441.17
##
## Number of Fisher Scoring iterations: 2###Singly add individual migration risk-factor: NET FOREST CHANGE
hm1.mc.Net.forest.change.2000<-glm(Trend00 ~ BrLon + BrLat + NBLon + NBLat +
scale(Net.forest.change.2000), family = "gaussian", data=lassodata50)
summary(hm1.mc.Net.forest.change.2000)##
## Call:
## glm(formula = Trend00 ~ BrLon + BrLat + NBLon + NBLat + scale(Net.forest.change.2000),
## family = "gaussian", data = lassodata50)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -4.9103 -1.6784 -1.0515 0.6454 13.8716
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -53.67851 6.14218 -8.739 4.66e-13 ***
## BrLon -0.07856 0.07273 -1.080 0.283532
## BrLat 1.00507 0.15167 6.627 4.61e-09 ***
## NBLon 0.09313 0.12205 0.763 0.447838
## NBLat 0.68577 0.19589 3.501 0.000785 ***
## scale(Net.forest.change.2000) 0.02948 0.40741 0.072 0.942516
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 12.33883)
##
## Null deviance: 2453.17 on 80 degrees of freedom
## Residual deviance: 925.41 on 75 degrees of freedom
## AIC: 441.17
##
## Number of Fisher Scoring iterations: 2#AIC table
cand.set<-list()
cand.set[[1]]<-hm1.mc
cand.set[[2]]<-hm1.mc.towers
cand.set[[3]]<-hm1.mc.Turbines
cand.set[[4]]<-hm1.mc.Hurricanes
cand.set[[5]]<-hm1.mc.Tornados
cand.set[[6]]<-hm1.mc.Agriculture
cand.set[[7]]<-hm1.mc.Forest.cover
cand.set[[8]]<-hm1.mc.Human.footprint
cand.set[[9]]<-hm1.mc.Net.forest.change.2000
#cand.set[[10]]<-hm1.mc.cum
cand.set[[10]]<-hm1.mc.full
#look at model table
AICcmodavg::aictab(cand.set, second.ord=TRUE, sort=TRUE)##
## Model selection based on AICc:
##
## K AICc Delta_AICc AICcWt Cum.Wt LL
## Mod1 6 440.31 0.00 0.28 0.28 -213.59
## Mod2 7 442.07 1.77 0.12 0.39 -213.27
## Mod4 7 442.45 2.14 0.10 0.49 -213.46
## Mod6 7 442.59 2.28 0.09 0.58 -213.53
## Mod7 7 442.69 2.38 0.08 0.66 -213.58
## Mod5 7 442.69 2.39 0.08 0.75 -213.58
## Mod3 7 442.69 2.39 0.08 0.83 -213.58
## Mod9 7 442.70 2.39 0.08 0.92 -213.58
## Mod8 7 442.71 2.40 0.08 1.00 -213.59
## Mod10 14 458.75 18.44 0.00 1.00 -212.19#calculate -2*LogLik
migmodtable<-AICcmodavg::aictab(cand.set, second.ord=TRUE, sort=TRUE)
migmodtable$LL*-2## [1] 427.1727 426.5395 426.9145 427.0564 427.1567 427.1587 427.1592 427.1671
## [9] 427.1727 424.3876Multiple linear regression (hierarchical GLMs); 10th percentile; annual exposure
#10 percentile
#intercept only (not compared)
hm1.null<-glm(Trend00 ~ 1, family = "gaussian", data=lassodata10)
summary(hm1.null)##
## Call:
## glm(formula = Trend00 ~ 1, family = "gaussian", data = lassodata10)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -8.6717 -0.9217 0.5183 1.7783 20.1383
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.03827 0.61528 -0.062 0.951
##
## (Dispersion parameter for gaussian family taken to be 30.6646)
##
## Null deviance: 2453.2 on 80 degrees of freedom
## Residual deviance: 2453.2 on 80 degrees of freedom
## AIC: 510.13
##
## Number of Fisher Scoring iterations: 2#Base model with nonbreeding and breeding latitude and longitude.
hm1.mc<-glm(Trend00 ~ BrLon + BrLat + NBLon + NBLat, family = "gaussian", data=lassodata10)
summary(hm1.mc)##
## Call:
## glm(formula = Trend00 ~ BrLon + BrLat + NBLon + NBLat, family = "gaussian",
## data = lassodata10)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -4.9513 -2.4565 -0.7162 0.3071 17.7234
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -91.81258 28.64812 -3.205 0.00198 **
## BrLon -0.22773 0.08553 -2.663 0.00946 **
## BrLat 0.76662 0.52368 1.464 0.14734
## NBLon -0.41368 0.13339 -3.101 0.00270 **
## NBLat 0.00661 0.22703 0.029 0.97685
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 18.91735)
##
## Null deviance: 2453.2 on 80 degrees of freedom
## Residual deviance: 1437.7 on 76 degrees of freedom
## AIC: 474.85
##
## Number of Fisher Scoring iterations: 2###FULL MODEL breeding/nonbreeding lat/lot + ALL MIG RISK TERMS
hm1.mc.full<-glm(Trend00 ~ BrLon + BrLat + NBLon + NBLat +
scale(Turbines) + scale(Towers) +
scale(Hurricanes) + scale(Tornados) +
scale(Agriculture) + scale(Forest.cover) +
scale(Human.footprint) + scale(Net.forest.change.2000),
family = "gaussian", data=lassodata10)
summary(hm1.mc.full)##
## Call:
## glm(formula = Trend00 ~ BrLon + BrLat + NBLon + NBLat + scale(Turbines) +
## scale(Towers) + scale(Hurricanes) + scale(Tornados) + scale(Agriculture) +
## scale(Forest.cover) + scale(Human.footprint) + scale(Net.forest.change.2000),
## family = "gaussian", data = lassodata10)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -5.9357 -2.1428 -0.9905 0.6397 16.2797
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -103.55000 32.92219 -3.145 0.00246 **
## BrLon -0.24037 0.10312 -2.331 0.02273 *
## BrLat 0.83658 0.59706 1.401 0.16571
## NBLon -0.50947 0.15215 -3.348 0.00133 **
## NBLat -0.04184 0.26214 -0.160 0.87367
## scale(Turbines) -0.29329 0.87104 -0.337 0.73737
## scale(Towers) -0.63533 0.67195 -0.946 0.34775
## scale(Hurricanes) 0.07626 0.84194 0.091 0.92810
## scale(Tornados) -1.08079 0.88057 -1.227 0.22392
## scale(Agriculture) 0.76338 1.23320 0.619 0.53797
## scale(Forest.cover) 1.09708 0.99218 1.106 0.27274
## scale(Human.footprint) 0.17655 1.36806 0.129 0.89770
## scale(Net.forest.change.2000) 0.14287 0.63194 0.226 0.82182
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 19.17985)
##
## Null deviance: 2453.2 on 80 degrees of freedom
## Residual deviance: 1304.2 on 68 degrees of freedom
## AIC: 482.96
##
## Number of Fisher Scoring iterations: 2###Singly add individual migration risk-factor: TOWERS
hm1.mc.towers<-glm(Trend00 ~ BrLon + BrLat + NBLon + NBLat +
scale(Towers), family = "gaussian", data=lassodata10)
summary(hm1.mc.towers)##
## Call:
## glm(formula = Trend00 ~ BrLon + BrLat + NBLon + NBLat + scale(Towers),
## family = "gaussian", data = lassodata10)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -4.7898 -2.3055 -0.8417 0.4149 17.6735
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -89.14935 28.92557 -3.082 0.00288 **
## BrLon -0.21713 0.08682 -2.501 0.01457 *
## BrLat 0.71457 0.52928 1.350 0.18105
## NBLon -0.42396 0.13438 -3.155 0.00231 **
## NBLat 0.00225 0.22769 0.010 0.99214
## scale(Towers) -0.38904 0.49966 -0.779 0.43866
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 19.01588)
##
## Null deviance: 2453.2 on 80 degrees of freedom
## Residual deviance: 1426.2 on 75 degrees of freedom
## AIC: 476.2
##
## Number of Fisher Scoring iterations: 2###Singly add individual migration risk-factor: TURBINES
hm1.mc.Turbines<-glm(Trend00 ~ BrLon + BrLat + NBLon + NBLat +
scale(Turbines), family = "gaussian", data=lassodata10)
summary(hm1.mc.Turbines)##
## Call:
## glm(formula = Trend00 ~ BrLon + BrLat + NBLon + NBLat + scale(Turbines),
## family = "gaussian", data = lassodata10)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -5.2680 -2.4824 -0.7739 0.4184 17.3160
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -96.11971 29.56776 -3.251 0.00172 **
## BrLon -0.23738 0.08723 -2.721 0.00808 **
## BrLat 0.81098 0.53049 1.529 0.13054
## NBLon -0.43260 0.13726 -3.152 0.00233 **
## NBLat -0.01051 0.22956 -0.046 0.96360
## scale(Turbines) -0.33495 0.53287 -0.629 0.53154
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 19.06913)
##
## Null deviance: 2453.2 on 80 degrees of freedom
## Residual deviance: 1430.2 on 75 degrees of freedom
## AIC: 476.43
##
## Number of Fisher Scoring iterations: 2###Singly add individual migration risk-factor: HURRICANES
hm1.mc.Hurricanes<-glm(Trend00 ~ BrLon + BrLat + NBLon + NBLat +
scale(Hurricanes), family = "gaussian", data=lassodata10)
summary(hm1.mc.Hurricanes)##
## Call:
## glm(formula = Trend00 ~ BrLon + BrLat + NBLon + NBLat + scale(Hurricanes),
## family = "gaussian", data = lassodata10)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -5.5434 -2.5115 -0.9535 0.2633 17.4476
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -92.51143 28.71300 -3.222 0.00188 **
## BrLon -0.22498 0.08575 -2.624 0.01053 *
## BrLat 0.74659 0.52518 1.422 0.15929
## NBLon -0.44904 0.13998 -3.208 0.00197 **
## NBLat -0.05948 0.24043 -0.247 0.80528
## scale(Hurricanes) 0.45047 0.53113 0.848 0.39907
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 18.98748)
##
## Null deviance: 2453.2 on 80 degrees of freedom
## Residual deviance: 1424.1 on 75 degrees of freedom
## AIC: 476.08
##
## Number of Fisher Scoring iterations: 2###Singly add individual migration risk-factor: TORNADOS
hm1.mc.Tornados<-glm(Trend00 ~ BrLon + BrLat + NBLon + NBLat +
scale(Tornados), family = "gaussian", data=lassodata10)
summary(hm1.mc.Tornados)##
## Call:
## glm(formula = Trend00 ~ BrLon + BrLat + NBLon + NBLat + scale(Tornados),
## family = "gaussian", data = lassodata10)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -4.9974 -2.4502 -0.6325 0.5680 17.3423
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -88.56931 28.76555 -3.079 0.00290 **
## BrLon -0.24586 0.08702 -2.825 0.00605 **
## BrLat 0.65889 0.53223 1.238 0.21959
## NBLon -0.40606 0.13340 -3.044 0.00322 **
## NBLat 0.07241 0.23461 0.309 0.75844
## scale(Tornados) -0.58325 0.53404 -1.092 0.27827
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 18.86949)
##
## Null deviance: 2453.2 on 80 degrees of freedom
## Residual deviance: 1415.2 on 75 degrees of freedom
## AIC: 475.58
##
## Number of Fisher Scoring iterations: 2###Singly add individual migration risk-factor: AGRICULTURE
hm1.mc.Agriculture<-glm(Trend00 ~ BrLon + BrLat + NBLon + NBLat +
scale(Agriculture), family = "gaussian", data=lassodata10)
summary(hm1.mc.Agriculture)##
## Call:
## glm(formula = Trend00 ~ BrLon + BrLat + NBLon + NBLat + scale(Agriculture),
## family = "gaussian", data = lassodata10)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -4.9441 -2.4587 -0.6793 0.2979 17.6630
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -91.66859 28.87874 -3.174 0.00218 **
## BrLon -0.23045 0.09097 -2.533 0.01339 *
## BrLat 0.75841 0.53454 1.419 0.16010
## NBLon -0.41323 0.13435 -3.076 0.00293 **
## NBLat 0.01004 0.23151 0.043 0.96553
## scale(Agriculture) -0.05064 0.54703 -0.093 0.92649
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 19.1674)
##
## Null deviance: 2453.2 on 80 degrees of freedom
## Residual deviance: 1437.6 on 75 degrees of freedom
## AIC: 476.84
##
## Number of Fisher Scoring iterations: 2###Singly add individual migration risk-factor: Forest and shrub cover
hm1.mc.Forest.cover<-glm(Trend00 ~ BrLon + BrLat + NBLon + NBLat +
scale(Forest.cover), family = "gaussian", data=lassodata10)
summary(hm1.mc.Forest.cover)##
## Call:
## glm(formula = Trend00 ~ BrLon + BrLat + NBLon + NBLat + scale(Forest.cover),
## family = "gaussian", data = lassodata10)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -5.6879 -2.3312 -0.8685 0.4046 16.7684
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -101.96458 28.94798 -3.522 0.000732 ***
## BrLon -0.24943 0.08550 -2.917 0.004657 **
## BrLat 0.87330 0.52138 1.675 0.098106 .
## NBLon -0.46136 0.13483 -3.422 0.001011 **
## NBLat -0.05087 0.22695 -0.224 0.823266
## scale(Forest.cover) 0.86680 0.51585 1.680 0.097052 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 18.47409)
##
## Null deviance: 2453.2 on 80 degrees of freedom
## Residual deviance: 1385.6 on 75 degrees of freedom
## AIC: 473.86
##
## Number of Fisher Scoring iterations: 2###Singly add individual migration risk-factor: HUMAN FOOTPRINT
hm1.mc.Human.footprint<-glm(Trend00 ~ BrLon + BrLat + NBLon + NBLat +
scale(Human.footprint), family = "gaussian", data=lassodata10)
summary(hm1.mc.Human.footprint)##
## Call:
## glm(formula = Trend00 ~ BrLon + BrLat + NBLon + NBLat + scale(Human.footprint),
## family = "gaussian", data = lassodata10)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -5.0366 -2.4331 -0.7973 0.2238 17.7753
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -94.854209 29.436580 -3.222 0.00188 **
## BrLon -0.229265 0.086013 -2.665 0.00941 **
## BrLat 0.813120 0.534570 1.521 0.13245
## NBLon -0.424595 0.135844 -3.126 0.00252 **
## NBLat -0.009965 0.230599 -0.043 0.96565
## scale(Human.footprint) 0.249765 0.503360 0.496 0.62121
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 19.10686)
##
## Null deviance: 2453.2 on 80 degrees of freedom
## Residual deviance: 1433.0 on 75 degrees of freedom
## AIC: 476.59
##
## Number of Fisher Scoring iterations: 2###Singly add individual migration risk-factor: NET FOREST CHANGE
hm1.mc.Net.forest.change.2000<-glm(Trend00 ~ BrLon + BrLat + NBLon + NBLat +
scale(Net.forest.change.2000), family = "gaussian", data=lassodata10)
summary(hm1.mc.Net.forest.change.2000)##
## Call:
## glm(formula = Trend00 ~ BrLon + BrLat + NBLon + NBLat + scale(Net.forest.change.2000),
## family = "gaussian", data = lassodata10)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -5.2829 -2.2872 -0.7859 0.2034 17.7103
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -89.854553 29.022795 -3.096 0.00276 **
## BrLon -0.226199 0.085987 -2.631 0.01034 *
## BrLat 0.714773 0.535259 1.335 0.18579
## NBLon -0.421634 0.134865 -3.126 0.00252 **
## NBLat 0.008393 0.228142 0.037 0.97075
## scale(Net.forest.change.2000) -0.269680 0.510648 -0.528 0.59898
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 19.09856)
##
## Null deviance: 2453.2 on 80 degrees of freedom
## Residual deviance: 1432.4 on 75 degrees of freedom
## AIC: 476.55
##
## Number of Fisher Scoring iterations: 2#AIC table
cand.set<-list()
cand.set[[1]]<-hm1.mc
cand.set[[2]]<-hm1.mc.towers
cand.set[[3]]<-hm1.mc.Turbines
cand.set[[4]]<-hm1.mc.Hurricanes
cand.set[[5]]<-hm1.mc.Tornados
cand.set[[6]]<-hm1.mc.Agriculture
cand.set[[7]]<-hm1.mc.Forest.cover
cand.set[[8]]<-hm1.mc.Human.footprint
cand.set[[9]]<-hm1.mc.Net.forest.change.2000
cand.set[[10]]<-hm1.mc.full
#view model table
AICcmodavg::aictab(cand.set, second.ord=TRUE, sort=TRUE)##
## Model selection based on AICc:
##
## K AICc Delta_AICc AICcWt Cum.Wt LL
## Mod7 7 475.39 0.00 0.26 0.26 -229.93
## Mod1 6 475.99 0.59 0.19 0.46 -231.43
## Mod5 7 477.11 1.72 0.11 0.57 -230.79
## Mod4 7 477.61 2.22 0.09 0.65 -231.04
## Mod2 7 477.74 2.34 0.08 0.73 -231.10
## Mod3 7 477.96 2.57 0.07 0.81 -231.21
## Mod9 7 478.09 2.69 0.07 0.87 -231.28
## Mod8 7 478.12 2.73 0.07 0.94 -231.29
## Mod6 7 478.38 2.98 0.06 1.00 -231.42
## Mod10 14 489.32 13.93 0.00 1.00 -227.48#calculate -2*logLik
migmodtable<-AICcmodavg::aictab(cand.set, second.ord=TRUE, sort=TRUE)
migmodtable$LL*-2## [1] 459.8601 462.8535 461.5755 462.0804 462.2014 462.4279 462.5529 462.5880
## [9] 462.8443 454.9605Power analysis of linear generalized regression models
#Linear regression
#for base model
library(pwr)
pwr.f2.test(u = 1, v = 79, f2 = 0.13/(1-0.13), sig.level = 0.05, power = ) # u = numerator df, v = denominator df##
## Multiple regression power calculation
##
## u = 1
## v = 79
## f2 = 0.1494253
## sig.level = 0.05
## power = 0.9300872Simple scatterplots of raw data versus population trends (Fig. S6)
Analysis of future threats
#Future threat analysis
threat25dat<-read.csv("/Users/gunnarkramer/Desktop/FINAL DATA FOR DRUM LAND ECOL/VermThreats4R_scale25_data.csv", header=TRUE)
colnames(threat25dat)## [1] "Name" "Species" "Pop" "TrendStat" "Trend00"
## [6] "BioFuels" "AgThreat" "SolarThreat" "UrbanThreat" "WindThreat"
## [11] "TotalMining"threat25dat$Name<-as.character(threat25dat$Name)
threat25dat$Species<-as.character(threat25dat$Species)
threat25dat$Pop<-as.character(threat25dat$Pop)
threat25dat$Trend00<-as.numeric(as.character(threat25dat$Trend00))
threat25dat$TrendGrp<-ifelse(threat25dat$Trend00 <=-1.5, "declining", ifelse(threat25dat$Trend00>-1.5 & threat25dat$Trend00<=1.5, "stable", ifelse(threat25dat$Trend00>1.5, "increasing", "NA")))
threat25dat<-subset(threat25dat, Species!="HYBRID")
threat25dat$spPop<-paste(threat25dat$Species, threat25dat$Pop, sep="")
colnames(threat25dat)## [1] "Name" "Species" "Pop" "TrendStat" "Trend00"
## [6] "BioFuels" "AgThreat" "SolarThreat" "UrbanThreat" "WindThreat"
## [11] "TotalMining" "TrendGrp" "spPop"threat25dat$IorS<-ifelse(threat25dat$TrendGrp!="declining", "increasing or stable", "declining")
futureMining<-lm(scale(TotalMining)~IorS, data=threat25dat)
summary(futureMining)##
## Call:
## lm(formula = scale(TotalMining) ~ IorS, data = threat25dat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.20650 -0.66255 0.02969 0.53691 2.94824
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.1553 0.2591 -0.599 0.551
## IorSincreasing or stable 0.1906 0.2870 0.664 0.509
##
## Residual standard error: 1.004 on 79 degrees of freedom
## Multiple R-squared: 0.005549, Adjusted R-squared: -0.007039
## F-statistic: 0.4408 on 1 and 79 DF, p-value: 0.5087futureBioFuels<-lm(scale(BioFuels)~IorS, data=threat25dat)
summary(futureBioFuels)##
## Call:
## lm(formula = scale(BioFuels) ~ IorS, data = threat25dat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.3910 -0.7986 -0.1145 0.5643 3.6726
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.07107 0.25968 0.274 0.785
## IorSincreasing or stable -0.08722 0.28768 -0.303 0.763
##
## Residual standard error: 1.006 on 79 degrees of freedom
## Multiple R-squared: 0.001162, Adjusted R-squared: -0.01148
## F-statistic: 0.09193 on 1 and 79 DF, p-value: 0.7625futureAg<-lm(scale(AgThreat)~IorS, data=threat25dat)
summary(futureAg)##
## Call:
## lm(formula = scale(AgThreat) ~ IorS, data = threat25dat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.14081 -0.66226 -0.02669 0.68870 2.70069
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.3721 0.2557 1.456 0.149
## IorSincreasing or stable -0.4567 0.2832 -1.613 0.111
##
## Residual standard error: 0.9901 on 79 degrees of freedom
## Multiple R-squared: 0.03187, Adjusted R-squared: 0.01961
## F-statistic: 2.6 on 1 and 79 DF, p-value: 0.1108futureSolar<-glm(scale(SolarThreat)~IorS, data=threat25dat)
summary(futureSolar)##
## Call:
## glm(formula = scale(SolarThreat) ~ IorS, data = threat25dat)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.7265 -0.6262 -0.0274 0.4537 3.7110
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.5427 0.2509 2.163 0.0335 *
## IorSincreasing or stable -0.6661 0.2779 -2.397 0.0189 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 0.9440155)
##
## Null deviance: 80.000 on 80 degrees of freedom
## Residual deviance: 74.577 on 79 degrees of freedom
## AIC: 229.18
##
## Number of Fisher Scoring iterations: 2futureUrban<-glm(scale(UrbanThreat)~IorS, data=threat25dat)
summary(futureUrban)##
## Call:
## glm(formula = scale(UrbanThreat) ~ IorS, data = threat25dat)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.9853 -0.5607 -0.0977 0.4496 3.5250
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.3576 0.2560 1.397 0.166
## IorSincreasing or stable -0.4389 0.2836 -1.548 0.126
##
## (Dispersion parameter for gaussian family taken to be 0.9828527)
##
## Null deviance: 80.000 on 80 degrees of freedom
## Residual deviance: 77.645 on 79 degrees of freedom
## AIC: 232.44
##
## Number of Fisher Scoring iterations: 2futureWind<-glm(scale(WindThreat)~IorS, data=threat25dat)
summary(futureWind)##
## Call:
## glm(formula = scale(WindThreat) ~ IorS, data = threat25dat)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.6673 -0.6824 -0.1691 0.4946 2.5129
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.4253 0.2544 1.672 0.0985 .
## IorSincreasing or stable -0.5219 0.2818 -1.852 0.0677 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 0.9705142)
##
## Null deviance: 80.000 on 80 degrees of freedom
## Residual deviance: 76.671 on 79 degrees of freedom
## AIC: 231.42
##
## Number of Fisher Scoring iterations: 2threat25dat$combined<-scale(threat25dat$SolarThreat)+scale(threat25dat$UrbanThreat)+scale(threat25dat$AgThreat)+scale(threat25dat$WindThreat)+scale(threat25dat$BioFuels)+scale(threat25dat$TotalMining)
futureComb<-glm(scale(combined)~IorS, data=threat25dat)
summary(futureComb)##
## Call:
## glm(formula = scale(combined) ~ IorS, data = threat25dat)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.11226 -0.75276 -0.00529 0.59682 2.77399
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.3584 0.2560 1.400 0.165
## IorSincreasing or stable -0.4399 0.2836 -1.551 0.125
##
## (Dispersion parameter for gaussian family taken to be 0.9827193)
##
## Null deviance: 80.000 on 80 degrees of freedom
## Residual deviance: 77.635 on 79 degrees of freedom
## AIC: 232.43
##
## Number of Fisher Scoring iterations: 2threat25dat$IorScol<-ifelse(threat25dat$IorS == "declining", "#B6322B", "#C0BCB5")
#PLOT
#png("FutureThreats_Boxplot.png", width = (6.85039*2), height = (6.85039+1), units = 'in', res = 300)
par(mfrow=c(2,3), mar=c(1,5,4,0), pty="s")
boxplot(scale(AgThreat)~IorS, data=threat25dat, col=rev(unique(threat25dat$IorScol)), ylim=c(-3,4), main="", ylab="", xlab="", yaxt='n', xaxt='n', cex.lab=1.5, cex.axis=1.5)
axis(2, at=c(-3,-2,-1,0,1,2,3,4), labels=c("-3","","","0","","","3",""), cex.axis=1.5)
mtext("Scaled exposure", side=2, line=3, cex=1.5)
mtext("Agriculture", side=3, line=.5, cex=2, font=2)
box(cex=1.25)
#legend("bottomleft", legend=c("Declining", "Increasing or stationary"), bty='n', col=rev(unique(threat25dat$IorScol)), pch=15, pt.cex=2, cex=1.25)
#text(4,3.5, labels="a", cex=1.25, font=2)
#text(5,3.5, labels="b", cex=1.25, font=2)
boxplot(scale(BioFuels)~IorS, data=threat25dat, col=rev(unique(threat25dat$IorScol)), ylim=c(-3,4), main="", ylab="", xlab="", yaxt='n', xaxt='n', cex.lab=1.5, cex.axis=1.5)
axis(2, at=c(-3,-2,-1,0,1,2,3,4), labels=c("-3","","","0","","","3",""), cex.axis=1.5)
#mtext("Scaled exposure", side=2, line=3, cex=1.25)
mtext("Biofuels", side=3, line=.5, cex=2, font=2)
box(cex=1.25)
#legend("bottomleft", legend=c("Declining", "Increasing or stationary"), bty='n', col=rev(unique(threat25dat$IorScol)), pch=15, pt.cex=2, cex=1.25)
#text(4,3.5, labels="a", cex=1.25, font=2)
#text(5,3.5, labels="b", cex=1.25, font=2)
boxplot(scale(TotalMining)~IorS, data=threat25dat, col=rev(unique(threat25dat$IorScol)), ylim=c(-3,4), main="", ylab="", xlab="", yaxt='n', xaxt='n', cex.lab=1.5, cex.axis=1.5)
axis(2, at=c(-3,-2,-1,0,1,2,3,4), labels=c("-3","","","0","","","3",""), cex.axis=1.5)
#mtext("Scaled exposure", side=2, line=3, cex=1.25)
mtext("Mining", side=3, line=.5, cex=2, font=2)
box(cex=1.25)
#legend("bottomleft", legend=c("Declining", "Increasing or stationary"), bty='n', col=rev(unique(threat25dat$IorScol)), pch=15, pt.cex=2, cex=1.25)
#text(4,3.5, labels="a", cex=1.25, font=2)
#text(5,3.5, labels="b", cex=1.25, font=2)
boxplot(scale(SolarThreat)~IorS, data=threat25dat, col=rev(unique(threat25dat$IorScol)), ylim=c(-3,4), main="", ylab="", xlab="", yaxt='n', xaxt='n', cex.lab=1.5, cex.axis=1.5)
axis(2, at=c(-3,-2,-1,0,1,2,3,4), labels=c("-3","","","0","","","3",""), cex.axis=1.5)
mtext("Scaled exposure", side=2, line=3, cex=1.5)
mtext("Solar", side=3, line=.5, cex=2, font=2)
box(cex=1.25)
legend("bottomleft", legend=c("Declining", "Increasing or stationary"), bty='n', col=rev(unique(threat25dat$IorScol)), pch=15, pt.cex=2, cex=2)
#legend("bottomleft", legend=c("Declining", "Increasing or stationary"), bty='n', col=rev(unique(threat25dat$IorScol)), pch=15, pt.cex=2, cex=1.25)
#text(4,3.5, labels="a", cex=1.25, font=2)
text(1.5,3, labels="*", cex=3, font=2)
boxplot(scale(UrbanThreat)~IorS, data=threat25dat, col=rev(unique(threat25dat$IorScol)), ylim=c(-3,4), main="", ylab="", xlab="", yaxt='n', xaxt='n', cex.lab=1.5, cex.axis=1.5)
axis(2, at=c(-3,-2,-1,0,1,2,3,4), labels=c("-3","","","0","","","3",""), cex.axis=1.5)
#mtext("Scaled exposure", side=2, line=3, cex=1.25)
mtext("Urban", side=3, line=.5, cex=2, font=2)
box(cex=1.25)
#legend("bottomleft", legend=c("Declining", "Increasing or stationary"), bty='n', col=rev(unique(threat25dat$IorScol)), pch=15, pt.cex=2, cex=1.25)
#text(4,3.5, labels="a", cex=1.25, font=2)
boxplot(scale(WindThreat)~IorS, data=threat25dat, col=rev(unique(threat25dat$IorScol)), ylim=c(-3,4), main="", ylab="", xlab="", yaxt='n', xaxt='n', cex.lab=1.5, cex.axis=1.5)
axis(2, at=c(-3,-2,-1,0,1,2,3,4), labels=c("-3","","","0","","","3",""), cex.axis=1.5)
#mtext("Scaled exposure", side=2, line=3, cex=1.25)
mtext("Wind", side=3, line=.5, cex=2, font=2)
box(cex=1.25)
#legend("bottomleft", legend=c("Declining", "Increasing or stationary"), bty='n', col=rev(unique(threat25dat$IorScol)), pch=15, pt.cex=2, cex=1.25)
#text(4,3.5, labels="a", cex=1.25, font=2)
#dev.off()PLS regression (annual exposure)
- 10% core-use area PLS regression
#PLS 10th Percentile Core-Use Area Data
#read data
lassodata10<-read.csv("/Users/gunnarkramer/Desktop/FINAL DATA FOR DRUM LAND ECOL/Data4Lasso10_data.csv", header=TRUE)
#housekeeping
lassodata10$Individual<-as.character(lassodata10$Individual)
lassodata10$Species<-as.character(lassodata10$Species)
lassodata10$Pop<-as.character(lassodata10$Pop)
lassodata10$Trend00<-as.numeric(as.character(lassodata10$Trend00))
lassodata10$TrendGrp<-ifelse(lassodata10$Trend00 <=-1.5, "Declining", ifelse(lassodata10$Trend00>-1.5 & lassodata10$Trend00<=1.5, "Stable", ifelse(lassodata10$Trend00>1.5, "Increasing", "NA")))
#omit hybrids
lassodata10<-subset(lassodata10, Species!="HYBRID")
lassodata10$Species<-ifelse(lassodata10$Species=="BW","Blue-winged warbler", "Golden-winged warbler")
#create new column with combined info for species and population
lassodata10$spPop<-paste(lassodata10$Species, lassodata10$Pop, sep="")
###now Partial Least Squares analysis
#####PLS
pls.lassodata10<-lassodata10
colnames(pls.lassodata10)## [1] "Individual" "Species" "Pop"
## [4] "State" "BrLat" "BrLon"
## [7] "NBLat" "NBLon" "BCR_Trend"
## [10] "Trend00" "scale" "Turbines"
## [13] "Towers" "Hurricanes" "Tornados"
## [16] "Agriculture" "Forest.cover" "Human.footprint"
## [19] "Net.forest.change.2000" "OVERALL2" "Br_Forest"
## [22] "Br_HFootP" "Br_forestchange" "Nb_Forest"
## [25] "Nb_HFootP" "Nb_forestchange" "TrendGrp"
## [28] "spPop"#include migration risk-factors, breeding, nonbreeding and population trend (Trend00)
pls.lassodata10<-pls.lassodata10[, c(10,12:19, 21:26)]
colnames(pls.lassodata10) #check selection - Trend00 + 14 other factors## [1] "Trend00" "Turbines" "Towers"
## [4] "Hurricanes" "Tornados" "Agriculture"
## [7] "Forest.cover" "Human.footprint" "Net.forest.change.2000"
## [10] "Br_Forest" "Br_HFootP" "Br_forestchange"
## [13] "Nb_Forest" "Nb_HFootP" "Nb_forestchange"set.seed(314)
#subset training dataset
training.samples<-pls.lassodata10$Trend00 %>%
createDataPartition(p=0.8, list=FALSE)
train.data<-pls.lassodata10[training.samples, ]
test.data<-pls.lassodata10[-training.samples, ]
set.seed(314)
#train PLS model
model<-train(
Trend00~., data=train.data, method="pls",
scale=TRUE, center=TRUE,
trControl = trainControl("cv", number=5),
tuneLength = 14
)Look at PLS model results.
#look at model results
plot(model)
model$bestTune ## ncomp
## 3 3summary(model$finalModel) ## Data: X dimension: 67 14
## Y dimension: 67 1
## Fit method: oscorespls
## Number of components considered: 3
## TRAINING: % variance explained
## 1 comps 2 comps 3 comps
## X 18.45 32.91 41.07
## .outcome 27.73 38.00 41.99#look at model performance
predictions<- model %>% predict(test.data)
data.frame(
RMSE = caret::RMSE(predictions, test.data$Trend00),
Rsquare = caret::R2(predictions, test.data$Trend00)
)## RMSE Rsquare
## 1 3.295406 0.4818282#PLS model performance vs. test data
lm.out<-lm(test.data$Trend00~predictions)
newx<-seq(-23, 23, by=0.05)
conf_int<-stats::predict(lm.out, newdata=data.frame(predictions=newx), interval="confidence", level=0.95)
summary(lm.out)##
## Call:
## lm(formula = test.data$Trend00 ~ predictions)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.5670 -2.3051 -0.1647 2.3893 5.2037
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.2906 0.8402 -1.536 0.15046
## predictions 0.7601 0.2275 3.340 0.00588 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.143 on 12 degrees of freedom
## Multiple R-squared: 0.4818, Adjusted R-squared: 0.4386
## F-statistic: 11.16 on 1 and 12 DF, p-value: 0.005884Plot PLS model performance.
#Plot performance of PLS model (predicted vs. observed)
#png("PLS_perf_10_plot.png", width = (6.85039), height = (6.85039), units = 'in', res = 600)
plot(test.data$Trend00~predictions, type='n', xlim=c(-20,20), ylim=c(-20,20), ylab="Observed population trajectory (percent annual change)", xlab="Predicted population trajectory (percent annual change)", pch=19, cex=1.5, col=rgb(158/255, 52/255, 235/255), main="10th percentile core-use area")
polygon(c(rev(newx), newx), c(rev(conf_int[ ,3]), conf_int[ ,2]), col=rgb(158/255, 52/255, 235/255, 0.2), border = NA)
points(test.data$Trend00~predictions, xlim=c(-20,20), ylim=c(-20,20), ylab="Observed population trajectory (percent annual change)", xlab="Predicted population trajectory (percent annual change)", pch=19, cex=1.5, col=rgb(158/255, 52/255, 235/255))
abline(1,1, col="gray", lty=1, lwd=3)
abline(lm(test.data$Trend00~predictions), col=rgb(158/255, 52/255, 235/255), lwd=3, lty=3)
#matlines(newx, conf_int[,2:3], col=rgb(158/255, 52/255, 235/255, 0.5), lty=1, lwd=1)
legend(-20, 20, legend=c("1:1 line", "PLS regression model"), lty=c(1,3), lwd=c(3,3), col=c("gray", rgb(158/255, 52/255, 235/255)), box.lty=0)
text(10,-10, expression(italic('R')^'2'* ' = 0.48'), cex=1) #estimated above
#dev.off()VIP vs. coefficient plot.
#VIP vs. coefficient plot
coefficients<-coef(model$finalModel)
coefficients<-data.frame(coefficients)
coefficients$var<-row.names(coefficients)
coefficients<-coefficients[order(coefficients$var),]
pls.var.imp<-varImp(model$finalModel)
pls.var.imp$var<-rownames(pls.var.imp)
pls.var.imp<-pls.var.imp[order(pls.var.imp$var),]
vip.coef<-data.frame(cbind(coefficients$var, coefficients$.outcome.3.comps, pls.var.imp$Overall))
colnames(vip.coef)<-c("var","coef","vip")
vip.coef$coef<-abs(as.numeric(as.character(vip.coef$coef)))
vip.coef$vip<-as.numeric(as.character(vip.coef$vip))
vip.coef$type<-c(rep("br", 3), rep("mrf", 4), rep("nb", 3), rep("mrf", 4))
vip.coef$color<-c(rep("#35B779", 3), rep("#FDE725",4), rep("#443A83", 3), rep("#FDE725",4))
vip.coef$var<-c( "Forest", "Net change forest", "Human footprint",
"Agriculture", "Forest", "Human footprint", "Hurricanes",
"Forest", "Net change forest", "Human footprint",
"Net forest change", "Towers","Tornados", "Turbines")
#VIP vs. COEF plot
#png("VIP_vs_COEF_PLS_plot_10percentile_BrNBMig_NoAgriculture.png", width = (6.85039), height = (6.85039), units = 'in', res = 600)
plot(vip~coef, data=vip.coef, main="10th percentile core-use area", col=color, pch=19, ylab="Variable importance for the projection (VIP)", xlab="Regression coefficient (absolute value)", cex=1.5, ylim=c(0,1.2), xlim=c(0,2.2))
legend(x=0.05, y=1.2, legend=c("Breeding factors", "Migration risk-factors", "Nonbreeding factors"), box.lty=0, pt.bg=c("#35B779", "#FDE725", "#443A83"), pch=21, pt.cex=2, cex=1)
abline(v=1, col="gray70", lty=2)
abline(h=0.8, col="gray70", lty=2)
text(vip~coef, data=vip.coef[c(1,4,7),], labels=vip.coef$var[c(1,4,7)], cex=0.6, pos=2) # left
text(vip~coef, data=vip.coef[c(2, 5, 8,9, 12,14),], labels=vip.coef$var[c(2,5, 8,9, 12,14)], cex=0.6, pos=4) #right
text(vip~coef, data=vip.coef[c(6,10,13),], labels=vip.coef$var[c(6,10,13)], cex=0.6, pos=3) #up
text(vip~coef, data=vip.coef[c(3,11),], labels=vip.coef$var[c(3,11)], cex=0.6, pos=1) #bottom
#dev.off()- 25% core-use area PLS regression
#25 percentile data
#read data
lassodata25<-read.csv("/Users/gunnarkramer/Desktop/FINAL DATA FOR DRUM LAND ECOL/Data4Lasso25_data.csv", header=TRUE)
colnames(lassodata25)## [1] "Individual" "Species" "Pop"
## [4] "State" "BrLat" "BrLon"
## [7] "NonBrLat" "NonBrLon" "BCR_Trend"
## [10] "Trend00" "scale" "Turbines"
## [13] "Towers" "Hurricanes" "Tornados"
## [16] "Agriculture" "Forest.cover" "Human.footprint"
## [19] "Net.forest.change.2000" "OVERALL2" "Br_Forest"
## [22] "Br_HFootP" "Br_forestchange" "Nb_Forest"
## [25] "Nb_HFootP" "Nb_forestchange"#housekeeping
lassodata25$Individual<-as.character(lassodata25$Individual)
lassodata25$Species<-as.character(lassodata25$Species)
lassodata25$Pop<-as.character(lassodata25$Pop)
lassodata25$State<-as.character(lassodata25$State)
lassodata25$Trend00<-as.numeric(as.character(lassodata25$Trend00))
lassodata25$TrendGrp<-ifelse(lassodata25$Trend00 <=-1.5, "Declining", ifelse(lassodata25$Trend00>-1.5 & lassodata25$Trend00<=1.5, "Stable", ifelse(lassodata25$Trend00>1.5, "Increasing", "NA")))
#omit hybrids
lassodata25<-subset(lassodata25, Species!="HYBRID")
lassodata25$Species<-ifelse(lassodata25$Species=="BW","Blue-winged warbler", "Golden-winged warbler")
lassodata25$spPop<-paste(lassodata25$Species, lassodata25$Pop, sep="")Select data
###now pls analysis
#####PLS
pls.lassodata25<-lassodata25
colnames(pls.lassodata25)## [1] "Individual" "Species" "Pop"
## [4] "State" "BrLat" "BrLon"
## [7] "NonBrLat" "NonBrLon" "BCR_Trend"
## [10] "Trend00" "scale" "Turbines"
## [13] "Towers" "Hurricanes" "Tornados"
## [16] "Agriculture" "Forest.cover" "Human.footprint"
## [19] "Net.forest.change.2000" "OVERALL2" "Br_Forest"
## [22] "Br_HFootP" "Br_forestchange" "Nb_Forest"
## [25] "Nb_HFootP" "Nb_forestchange" "TrendGrp"
## [28] "spPop"#select relevant migration risk-factors, breeding factors, nonbreeding factors, and Trend00
pls.lassodata25<-pls.lassodata25[, c(10, 12:19, 21:26)]
colnames(pls.lassodata25)## [1] "Trend00" "Turbines" "Towers"
## [4] "Hurricanes" "Tornados" "Agriculture"
## [7] "Forest.cover" "Human.footprint" "Net.forest.change.2000"
## [10] "Br_Forest" "Br_HFootP" "Br_forestchange"
## [13] "Nb_Forest" "Nb_HFootP" "Nb_forestchange"Check VIF of 25th percentile factors
#Variable inflation factor (VIF) analysis
colnames(pls.lassodata25)## [1] "Trend00" "Turbines" "Towers"
## [4] "Hurricanes" "Tornados" "Agriculture"
## [7] "Forest.cover" "Human.footprint" "Net.forest.change.2000"
## [10] "Br_Forest" "Br_HFootP" "Br_forestchange"
## [13] "Nb_Forest" "Nb_HFootP" "Nb_forestchange"mc.lassodata25<-pls.lassodata25[,-c(1)] #only looking at correlation of explanatory variables, don't need Trend00
vif(mc.lassodata25) #need to remove "Trend00" from dataset before running## Variables VIF
## 1 Turbines 5.657621
## 2 Towers 2.688498
## 3 Hurricanes 3.267692
## 4 Tornados 5.670960
## 5 Agriculture 9.272558
## 6 Forest.cover 9.386344
## 7 Human.footprint 15.150604
## 8 Net.forest.change.2000 2.253778
## 9 Br_Forest 1.523311
## 10 Br_HFootP 2.110755
## 11 Br_forestchange 1.802065
## 12 Nb_Forest 1.922401
## 13 Nb_HFootP 1.338296
## 14 Nb_forestchange 1.397629Run PLS model
set.seed(314)
#subset training data
training.samples<-pls.lassodata25$Trend00 %>%
createDataPartition(p=0.8, list=FALSE)
train.data<-pls.lassodata25[training.samples, ]
test.data<-pls.lassodata25[-training.samples, ]
set.seed(314)
#train model
model<-train(
Trend00 ~., data=train.data, method="pls",
scale=TRUE, center=TRUE,
trControl = trainControl("cv", number=5),
tuneLength = 14
)
plot(model)
model$bestTune ## ncomp
## 2 2summary(model$finalModel) ## Data: X dimension: 67 14
## Y dimension: 67 1
## Fit method: oscorespls
## Number of components considered: 2
## TRAINING: % variance explained
## 1 comps 2 comps
## X 20.69 32.22
## .outcome 28.15 39.83predictions<- model %>% predict(test.data)
data.frame(
RMSE = caret::RMSE(predictions, test.data$Trend00),
Rsquare = caret::R2(predictions, test.data$Trend00)
)## RMSE Rsquare
## 1 3.961513 0.2986111Check and plot model performance
#PLS model performance vs. test data
lm.out<-lm(test.data$Trend00~predictions)
newx<-seq(-23, 23, by=0.05)
conf_int<-stats::predict(lm.out, newdata=data.frame(predictions=newx), interval="confidence", level=0.95)
summary(lm.out)##
## Call:
## lm(formula = test.data$Trend00 ~ predictions)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.3265 -2.7203 0.6776 2.2740 5.6985
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.9382 1.0101 -1.919 0.0791 .
## predictions 0.8648 0.3826 2.260 0.0432 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.656 on 12 degrees of freedom
## Multiple R-squared: 0.2986, Adjusted R-squared: 0.2402
## F-statistic: 5.109 on 1 and 12 DF, p-value: 0.04319#plot model performance (predicted vs. observed)
#png("PLS_perf_plot25_revised_9jan23_noAgriculture.png", width = (6.85039), height = (6.85039), units = 'in', res = 600)
plot(test.data$Trend00~predictions, type='n', xlim=c(-20,20), ylim=c(-20,20), ylab="Observed population trajectory (percent annual change)", xlab="Predicted population trajectory (percent annual change)", pch=19, cex=1.5, col=rgb(158/255, 52/255, 235/255), main="25th percentile core-use area")
polygon(c(rev(newx), newx), c(rev(conf_int[ ,3]), conf_int[ ,2]), col=rgb(158/255, 52/255, 235/255, 0.2), border = NA)
points(test.data$Trend00~predictions, xlim=c(-20,20), ylim=c(-20,20), ylab="Observed population trajectory (percent annual change)", xlab="Predicted population trajectory (percent annual change)", pch=19, cex=1.5, col=rgb(158/255, 52/255, 235/255))
abline(1,1, col="gray", lty=1, lwd=3)
abline(lm(test.data$Trend00~predictions), col=rgb(158/255, 52/255, 235/255), lwd=3, lty=3)
#matlines(newx, conf_int[,2:3], col=rgb(158/255, 52/255, 235/255, 0.5), lty=1, lwd=1)
legend(-20, 20, legend=c("1:1 line", "PLS regression model"), lty=c(1,3), lwd=c(3,3), col=c("gray", rgb(158/255, 52/255, 235/255)), box.lty=0)
text(10,-10, expression(italic('R')^'2'* ' = 0.30'), cex=1)
#dev.off()Calculate and plot VIP vs. coefficient estimates
coefficients<-coef(model$finalModel)
sum.coef<-sum(sapply(coefficients, abs))
coefficients<-coefficients*100/sum.coef
coefficients<-sort(coefficients[,1,1])
abs.coefficients<-rev(sort(abs(coefficients)))
# coefficients
sum(abs(coefficients))## [1] 100sort(abs(coefficients))## Tornados Hurricanes Turbines
## 0.1682573 2.5522789 2.9556321
## Agriculture Net.forest.change.2000 Human.footprint
## 3.2920171 3.8629477 5.5956061
## Br_HFootP Forest.cover Br_forestchange
## 6.9675081 7.1856068 7.6043783
## Towers Nb_forestchange Nb_HFootP
## 8.0141879 11.2927140 12.8160627
## Nb_Forest Br_Forest
## 13.1235824 14.5692207# VIP
pls.var.imp<-varImp(model$finalModel)
pls.var.imp<-as.data.frame(pls.var.imp)
pls.var.imp<-rev(sort(pls.var.imp[,1,1]))
#VIP vs. coefficient plot
coefficients<-coef(model$finalModel)
coefficients<-data.frame(coefficients)
coefficients$var<-row.names(coefficients)
coefficients<-coefficients[order(coefficients$var),]
pls.var.imp<-varImp(model$finalModel)
pls.var.imp$var<-rownames(pls.var.imp)
pls.var.imp<-pls.var.imp[order(pls.var.imp$var),]
vip.coef<-data.frame(cbind(coefficients$var, coefficients$.outcome.2.comps, pls.var.imp$Overall))
colnames(vip.coef)<-c("var","coef","vip")
vip.coef$coef<-abs(as.numeric(as.character(vip.coef$coef)))
vip.coef$vip<-as.numeric(as.character(vip.coef$vip))
vip.coef$type<-c("mrf", rep("br", 3), rep("mrf", 3), rep("nb", 3), rep("mrf", 4)) #br = breeding factor, mfr = migration risk factor, nb = nonbreeding factor
vip.coef$color<-c("#FDE725", rep("#35B779",3), rep("#FDE725", 3), rep("#443A83",3), rep("#FDE725", 4))
vip.coef$var<-c("Agriculture", "Forest", "Net change forest",
"Human footprint", "Forest", "Human footprint", "Hurricanes",
"Forest", "Net change forest", "Human footprint",
"Net forest change", "Tornados","Towers", "Turbines")
#VIP vs. COEF plot
#png("VIP_vs_COEF_PLS_plot_25percentile_BrNbMig_NoAg.png", width = (6.85039), height = (6.85039), units = 'in', res = 600)
plot(vip~coef, data=vip.coef, main="", col=color, pch=19, ylab="Variable importance for the projection (VIP)", xlab="Regression coefficient (absolute value)", cex=1.5, ylim=c(0,1), xlim=c(0,1.8))
legend(x=0, y=1.03, legend=c("Breeding factors", "Migration-risk factors", "Nonbreeding factors"), box.lty=0, pt.bg=c("#35B779", "#FDE725", "#443A83"), pch=21, pt.cex=2, cex=1)
abline(h=0.8, col="gray70", lty=2)
abline(v=1, col="gray70", lty=2)
text(vip~coef, data=vip.coef[c(4,7,13,14,15,16),], labels=vip.coef$var[c(4,7,13,14,15,16)], cex=0.6, pos=2) # left
text(vip~coef, data=vip.coef[c(1,2,3, 5, 6,9, 12),], labels=vip.coef$var[c(1,2,3,5,6,9, 12)], cex=0.6, pos=4) #right
text(vip~coef, data=vip.coef[c(8,10,11),], labels=vip.coef$var[c(8,10,11)], cex=0.6, pos=3) #up
#text(vip~coef, data=vip.coef[c(1),], labels=vip.coef$var[c(1)], cex=0.6, pos=1) #bottom
#dev.off()- 50% core-use area PLS regression
#PLS 50th Percentile Data
#read data
lassodata50<-read.csv("/Users/gunnarkramer/Desktop/FINAL DATA FOR DRUM LAND ECOL/Data4Lasso50_data.csv", header=TRUE)
#housekeeping
colnames(lassodata50)## [1] "Individual" "Species" "Pop"
## [4] "State" "BrLat" "BrLon"
## [7] "NBLat" "NBLon" "BCR_Trend"
## [10] "Trend00" "scale" "Turbines"
## [13] "Towers" "Hurricanes" "Tornados"
## [16] "Agriculture" "Forest.cover" "Human.footprint"
## [19] "Net.forest.change.2000" "OVERALL2" "Br_Forest"
## [22] "Br_HFootP" "Br_forestchange" "Nb_Forest"
## [25] "Nb_HFootP" "Nb_forestchange"lassodata50$Individual<-as.character(lassodata50$Individual)
lassodata50$Species<-as.character(lassodata50$Species)
lassodata50$Pop<-as.character(lassodata50$Pop)
lassodata50$Trend00<-as.numeric(as.character(lassodata50$Trend00))
lassodata50$TrendGrp<-ifelse(lassodata50$Trend00 <=-1.5, "Declining", ifelse(lassodata50$Trend00>-1.5 & lassodata50$Trend00<=1.5, "Stable", ifelse(lassodata50$Trend00>1.5, "Increasing", "NA")))
#omit hybrids
lassodata50<-subset(lassodata50, Species!="HYBRID")
lassodata50$Species<-ifelse(lassodata50$Species=="BW","Blue-winged warbler", "Golden-winged warbler")
lassodata50$spPop<-paste(lassodata50$Species, lassodata50$Pop, sep="")Select variables and run PLS regression.
#PLS
pls.lassodata50<-lassodata50
colnames(pls.lassodata50)## [1] "Individual" "Species" "Pop"
## [4] "State" "BrLat" "BrLon"
## [7] "NBLat" "NBLon" "BCR_Trend"
## [10] "Trend00" "scale" "Turbines"
## [13] "Towers" "Hurricanes" "Tornados"
## [16] "Agriculture" "Forest.cover" "Human.footprint"
## [19] "Net.forest.change.2000" "OVERALL2" "Br_Forest"
## [22] "Br_HFootP" "Br_forestchange" "Nb_Forest"
## [25] "Nb_HFootP" "Nb_forestchange" "TrendGrp"
## [28] "spPop"#include only migration risk-factors, breeding factors, nonbreeding factors and Trend00
pls.lassodata50<-pls.lassodata50[, c(10, 12:19, 21:26)]
colnames(pls.lassodata50) # Trend00 + 14 factors## [1] "Trend00" "Turbines" "Towers"
## [4] "Hurricanes" "Tornados" "Agriculture"
## [7] "Forest.cover" "Human.footprint" "Net.forest.change.2000"
## [10] "Br_Forest" "Br_HFootP" "Br_forestchange"
## [13] "Nb_Forest" "Nb_HFootP" "Nb_forestchange"set.seed(314)
#create training subset of data
training.samples<-pls.lassodata50$Trend00 %>%
createDataPartition(p=0.8, list=FALSE)
train.data<-pls.lassodata50[training.samples, ]
test.data<-pls.lassodata50[-training.samples, ]
set.seed(314)
#train PLS model
model<-train(
Trend00~., data=train.data, method="pls",
scale=TRUE, center=TRUE,
trControl = trainControl("cv", number=5),
tuneLength = 114
)
#view model summary
plot(model)
model$bestTune ## ncomp
## 3 3summary(model$finalModel) ## Data: X dimension: 67 14
## Y dimension: 67 1
## Fit method: oscorespls
## Number of components considered: 3
## TRAINING: % variance explained
## 1 comps 2 comps 3 comps
## X 21.05 32.72 56.67
## .outcome 31.70 39.47 40.62Look at performance and plot.
predictions<- model %>% predict(test.data)
data.frame(
RMSE = caret::RMSE(predictions, test.data$Trend00),
Rsquare = caret::R2(predictions, test.data$Trend00)
)## RMSE Rsquare
## 1 3.358473 0.5086996#PLS model performance vs. test data
lm.out<-lm(test.data$Trend00~predictions)
newx<-seq(-23, 23, by=0.05)
conf_int<-stats::predict(lm.out, newdata=data.frame(predictions=newx), interval="confidence", level=0.95)
summary(lm.out)##
## Call:
## lm(formula = test.data$Trend00 ~ predictions)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.309 -2.069 -0.642 2.804 4.175
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.7609 0.8257 -2.133 0.05431 .
## predictions 0.9309 0.2641 3.525 0.00419 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.06 on 12 degrees of freedom
## Multiple R-squared: 0.5087, Adjusted R-squared: 0.4678
## F-statistic: 12.42 on 1 and 12 DF, p-value: 0.004185#png("PLS_perf50.png", width = (6.85039), height = (6.85039), units = 'in', res = 600)
plot(test.data$Trend00~predictions, type='n', xlim=c(-20,20), ylim=c(-20,20), ylab="Observed population trajectory (percent annual change)", xlab="Predicted population trajectory (percent annual change)", pch=19, cex=1.5, col=rgb(158/255, 52/255, 235/255), main="50th percentile core-use area")
polygon(c(rev(newx), newx), c(rev(conf_int[ ,3]), conf_int[ ,2]), col=rgb(158/255, 52/255, 235/255, 0.2), border = NA)
points(test.data$Trend00~predictions, xlim=c(-20,20), ylim=c(-20,20), ylab="Observed population trajectory (percent annual change)", xlab="Predicted population trajectory (percent annual change)", pch=19, cex=1.5, col=rgb(158/255, 52/255, 235/255))
abline(1,1, col="gray", lty=1, lwd=3)
abline(lm(test.data$Trend00~predictions), col=rgb(158/255, 52/255, 235/255), lwd=3, lty=3)
#matlines(newx, conf_int[,2:3], col=rgb(158/255, 52/255, 235/255, 0.5), lty=1, lwd=1)
legend(-20, 20, legend=c("1:1 line", "PLS regression model"), lty=c(1,3), lwd=c(3,3), col=c("gray", rgb(158/255, 52/255, 235/255)), box.lty=0)
text(10,-10, expression(italic('R')^'2'* ' = 0.51'), cex=1)
#dev.off()VIP vs. coefficient plot
#VIP vs. coefficient plot
coefficients<-coef(model$finalModel)
coefficients<-data.frame(coefficients)
coefficients$var<-row.names(coefficients)
coefficients<-coefficients[order(coefficients$var),]
pls.var.imp<-varImp(model$finalModel)
pls.var.imp$var<-rownames(pls.var.imp)
pls.var.imp<-pls.var.imp[order(pls.var.imp$var),]
vip.coef<-data.frame(cbind(coefficients$var, coefficients$.outcome.3.comps, pls.var.imp$Overall))
colnames(vip.coef)<-c("var","coef","vip")
vip.coef$coef<-abs(as.numeric(as.character(vip.coef$coef)))
vip.coef$vip<-as.numeric(as.character(vip.coef$vip))
vip.coef$type<-c(rep("br", 3), rep("mrf", 4), rep("nb", 3), rep("mrf", 4))
vip.coef$color<-c(rep("#35B779",3), rep("#FDE725", 4), rep("#443A83",3), rep("#FDE725", 4))
vip.coef$var<-c("Forest", "Net change forest", "Human footprint",
"Agriculture", "Forest", "Human footprint", "Hurricanes",
"Forest", "Net change forest", "Human footprint",
"Net forest change", "Towers","Tornados", "Turbines")
#VIP vs. COEF plot
#png("VIP_vs_COEF_PLS_plot_50percentile_BrNBMig_updatedNOAgriculture.png", width = (6.85039), height = (6.85039), units = 'in', res = 600)
plot(vip~coef, data=vip.coef, main="50th percentile core-use area", col=color, pch=19, ylab="Variable importance for the projection (VIP)", xlab="Regression coefficient (absolute value)", cex=1.5, xlim=c(0,1.5), ylim=c(0,1.0))
legend(x=0, y=1.04, legend=c("Breeding factors", "Migration risk-factors", "Nonbreeding factors"), box.lty=0, pt.bg=c("#35B779", "#FDE725", "#443A83"), pch=21, pt.cex=2, cex=1)
abline(h=0.8, col="gray70", lty=2)
abline(v=1, col="gray70", lty=2)
text(vip~coef, data=vip.coef[c(1,9,14),], labels=vip.coef$var[c(1,9,14)], cex=0.6, pos=2) # left
text(vip~coef, data=vip.coef[c(3, 4, 5, 6, 7, 8,12,13),], labels=vip.coef$var[c(3,4,5,6,7, 8,12,13)], cex=0.6, pos=4) #right
text(vip~coef, data=vip.coef[c(2,10),], labels=vip.coef$var[c(2,10)], cex=0.6, pos=3) #up
text(vip~coef, data=vip.coef[c(11),], labels=vip.coef$var[c(11)], cex=0.6, pos=1) #bottom
#dev.off()PLS regression (Gulf of Mexico)
- Gulf of Mexico - autumn migration
Read data
Fall.Gulf.Risk<-read.csv("/Users/gunnarkramer/Desktop/FINAL DATA FOR DRUM LAND ECOL/FallMigration_PREGULF_Risk4R_25_data.csv", header=TRUE)
colnames(Fall.Gulf.Risk)## [1] "Individual" "Pop" "Species" "BrLat"
## [5] "BrLon" "NonBrLat" "NonBrLon" "Trend00"
## [9] "Season" "scale" "PreGulfArea" "Turbine"
## [13] "TallTowers" "Hurr" "TotalTorn" "CropCov"
## [17] "ForestShrubCov" "HFootP" "NetForest" "Br_Forest"
## [21] "Br_HFootP" "Br_forestchange" "Nb_Forest" "Nb_HFootP"
## [25] "Nb_forestchange"#housekeeping
Fall.Gulf.Risk$Individual<-as.character(Fall.Gulf.Risk$Individual)
Fall.Gulf.Risk$Species<-as.character(Fall.Gulf.Risk$Species)
Fall.Gulf.Risk$Pop<-as.character(Fall.Gulf.Risk$Pop)
Fall.Gulf.Risk$Trend00<-as.numeric(as.character(Fall.Gulf.Risk$Trend00))
#omit hybrids, individuals that did not have core-use areas within area of interest near Gulf of Mexico
Fall.Gulf.Risk<-subset(Fall.Gulf.Risk, Species!="HYBRID")
Fall.Gulf.Risk<-subset(Fall.Gulf.Risk, PreGulfArea=="YES")Subset data and train model
#####PLS
pls.fall.25<-Fall.Gulf.Risk
colnames(pls.fall.25)## [1] "Individual" "Pop" "Species" "BrLat"
## [5] "BrLon" "NonBrLat" "NonBrLon" "Trend00"
## [9] "Season" "scale" "PreGulfArea" "Turbine"
## [13] "TallTowers" "Hurr" "TotalTorn" "CropCov"
## [17] "ForestShrubCov" "HFootP" "NetForest" "Br_Forest"
## [21] "Br_HFootP" "Br_forestchange" "Nb_Forest" "Nb_HFootP"
## [25] "Nb_forestchange"#use only relevant migration risk-factors, breeding factors, nonbreeding factors and Trend00
pls.fall.25<-pls.fall.25[, c(8, 12:25)]
colnames(pls.fall.25) # Trend00 + 14 explanatory variables## [1] "Trend00" "Turbine" "TallTowers" "Hurr"
## [5] "TotalTorn" "CropCov" "ForestShrubCov" "HFootP"
## [9] "NetForest" "Br_Forest" "Br_HFootP" "Br_forestchange"
## [13] "Nb_Forest" "Nb_HFootP" "Nb_forestchange"set.seed(314)
#create training subset of data
training.samples<-pls.fall.25$Trend00 %>%
createDataPartition(p=0.8, list=FALSE)
train.data<-pls.fall.25[training.samples, ]
test.data<-pls.fall.25[-training.samples, ]
set.seed(314)
#train PLS model
model<-train(
Trend00~., data=train.data, method="pls",
scale=TRUE, center=TRUE,
trControl = trainControl("cv", number=5),
tuneLength = 14
)
#view model summary
plot(model)
model$bestTune ## ncomp
## 2 2summary(model$finalModel) ## Data: X dimension: 64 14
## Y dimension: 64 1
## Fit method: oscorespls
## Number of components considered: 2
## TRAINING: % variance explained
## 1 comps 2 comps
## X 25.98 47.57
## .outcome 25.67 36.79predictions<- model %>% predict(test.data)
data.frame(
RMSE = caret::RMSE(predictions, test.data$Trend00),
Rsquare = caret::R2(predictions, test.data$Trend00)
)## RMSE Rsquare
## 1 3.189483 0.563569#PLS model performance vs. test data
lm.out<-lm(test.data$Trend00~predictions)
newx<-seq(-23, 23, by=0.05)
conf_int<-stats::predict(lm.out, newdata=data.frame(predictions=newx), interval="confidence", level=0.95)
summary(lm.out)##
## Call:
## lm(formula = test.data$Trend00 ~ predictions)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.304 -2.237 -0.208 1.633 4.888
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.7577 0.7774 -2.261 0.04314 *
## predictions 1.0556 0.2682 3.936 0.00198 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.884 on 12 degrees of freedom
## Multiple R-squared: 0.5636, Adjusted R-squared: 0.5272
## F-statistic: 15.5 on 1 and 12 DF, p-value: 0.001975Plot model performance
#png("PLS_perf_plot25FallGulf.png", width = (6.85039), height = (6.85039), units = 'in', res = 600)
plot(test.data$Trend00~predictions, type='n', xlim=c(-20,20), ylim=c(-20,20), ylab="Observed population trajectory (percent annual change)", xlab="Predicted population trajectory (percent annual change)", pch=19, cex=1.5, col=rgb(158/255, 52/255, 235/255), main="Gulf of Mexico, autumn")
polygon(c(rev(newx), newx), c(rev(conf_int[ ,3]), conf_int[ ,2]), col=rgb(158/255, 52/255, 235/255, 0.2), border = NA)
points(test.data$Trend00~predictions, xlim=c(-20,20), ylim=c(-20,20), ylab="Observed population trajectory (percent annual change)", xlab="Predicted population trajectory (percent annual change)", pch=19, cex=1.5, col=rgb(158/255, 52/255, 235/255))
abline(1,1, col="gray", lty=1, lwd=3)
abline(lm(test.data$Trend00~predictions), col=rgb(158/255, 52/255, 235/255), lwd=3, lty=3)
#matlines(newx, conf_int[,2:3], col=rgb(158/255, 52/255, 235/255, 0.5), lty=1, lwd=1)
legend(-20, 20, legend=c("1:1 line", "PLS regression model"), lty=c(1,3), lwd=c(3,3), col=c("gray", rgb(158/255, 52/255, 235/255)), box.lty=0)
text(10,-10, expression(italic('R')^'2'* ' = 0.56'), cex=1)
#dev.off()Plot VIP vs. coefficient estimates
#VIP vs. coefficient plot
coefficients<-coef(model$finalModel)
coefficients<-data.frame(coefficients)
coefficients$var<-row.names(coefficients)
coefficients<-coefficients[order(coefficients$var),]
pls.var.imp<-varImp(model$finalModel)
pls.var.imp$var<-rownames(pls.var.imp)
pls.var.imp<-pls.var.imp[order(pls.var.imp$var),]
vip.coef<-data.frame(cbind(coefficients$var, coefficients$.outcome.2.comps, pls.var.imp$Overall))
colnames(vip.coef)<-c("var","coef","vip")
vip.coef$coef<-abs(as.numeric(as.character(vip.coef$coef)))
vip.coef$vip<-as.numeric(as.character(vip.coef$vip))
vip.coef$type<-c(rep("br", 3), rep("mrf", 4), rep("nb", 3), rep("mrf", 4)) #br = breeding factor, mrf = migration risk factor, nb = nonbreeding factor
vip.coef$color<-c(rep("#35B779",3), rep("#FDE725", 4), rep("#443A83",3), rep("#FDE725", 4))
vip.coef$var<-c("Forest", "Net change forest", "Human footprint",
"Agriculture", "Forest", "Human footprint", "Hurricanes",
"Forest", "Net change forest", "Human footprint",
"Net forest change", "Towers","Tornados", "Turbines")
#VIP vs. COEF plot
#png("VIP_vs_COEF_PLS_plot_25percentile_BrNBMig_FallPreGulf.png", width = (6.85039), height = (6.85039), units = 'in', res = 600)
plot(vip~coef, data=vip.coef, main="Gulf of Mexico, autumn", col=color, pch=19, ylab="Variable importance for the projection (VIP)", xlab="Regression coefficient (absolute value)", cex=1.5, xlim=c(0,1.4), ylim=c(0,1.0))
legend(x=0, y=1.04, legend=c("Breeding factors", "Migration risk-factors", "Nonbreeding factors"), box.lty=0, pt.bg=c("#35B779", "#FDE725", "#443A83"), pch=21, pt.cex=2, cex=1)
abline(h=0.8, col="gray70", lty=2)
abline(v=1, col="gray70", lty=2)
text(vip~coef, data=vip.coef[c(8,10,11),], labels=vip.coef$var[c(8,10,11)], cex=0.6, pos=2) # left
text(vip~coef, data=vip.coef[c(3, 6, 7),], labels=vip.coef$var[c(3,6,7)], cex=0.6, pos=4) #right
text(vip~coef, data=vip.coef[c(2,4,5,9,12,13,14),], labels=vip.coef$var[c(2,4,5,9,12,13,14)], cex=0.6, pos=3) #up
text(vip~coef, data=vip.coef[c(1),], labels=vip.coef$var[c(1)], cex=0.6, pos=1) #bottom
#explains 48% X variables; 38% outcome; R2 = 0.56
#dev.off()- Gulf of Mexico - spring migration
################SPRING MIGRATION
#read data
Spring.Gulf.Risk<-read.csv("/Users/gunnarkramer/Desktop/FINAL DATA FOR DRUM LAND ECOL/SpringMigration_PREGULF_Risk4R_data.csv", header=TRUE)
colnames(Spring.Gulf.Risk)## [1] "Individual" "Pop" "Species" "BrLat"
## [5] "BrLon" "NonBrLat" "NonBrLon" "Trend00"
## [9] "Season" "PreGulfArea" "CropCov" "ForestShrubCov"
## [13] "HFootP" "NetForest" "Br_Forest" "Br_HFootP"
## [17] "Br_forestchange" "Nb_Forest" "Nb_HFootP" "Nb_forestchange"#housekeeping
Spring.Gulf.Risk$Individual<-as.character(Spring.Gulf.Risk$Individual)
Spring.Gulf.Risk$Species<-as.character(Spring.Gulf.Risk$Species)
Spring.Gulf.Risk$Pop<-as.character(Spring.Gulf.Risk$Pop)
Spring.Gulf.Risk$Trend00<-as.numeric(as.character(Spring.Gulf.Risk$Trend00))
Spring.Gulf.Risk<-subset(Spring.Gulf.Risk, Species!="HYBRID")
Spring.Gulf.Risk<-subset(Spring.Gulf.Risk, PreGulfArea=="YES")
colnames(Spring.Gulf.Risk)## [1] "Individual" "Pop" "Species" "BrLat"
## [5] "BrLon" "NonBrLat" "NonBrLon" "Trend00"
## [9] "Season" "PreGulfArea" "CropCov" "ForestShrubCov"
## [13] "HFootP" "NetForest" "Br_Forest" "Br_HFootP"
## [17] "Br_forestchange" "Nb_Forest" "Nb_HFootP" "Nb_forestchange"Select relevant variables and train and run PLS regression
#####PLS
pls.spr.25<-Spring.Gulf.Risk
#use only relevant migration risk-factors, breeding, and nonbreeding factors. Some data sources not available for relevant regions near Gulf of Mexico in spring. See SI Methods and Table 2 in Kramer et al. 2023 for details.
pls.spr.25<-pls.spr.25[, c(8, 11:20)]
colnames(pls.spr.25) #Trend00 + 10 explanatory variables## [1] "Trend00" "CropCov" "ForestShrubCov" "HFootP"
## [5] "NetForest" "Br_Forest" "Br_HFootP" "Br_forestchange"
## [9] "Nb_Forest" "Nb_HFootP" "Nb_forestchange"set.seed(314)
#create training subset of data
training.samples<-pls.spr.25$Trend00 %>%
createDataPartition(p=0.8, list=FALSE)
train.data<-pls.spr.25[training.samples, ]
test.data<-pls.spr.25[-training.samples, ]
set.seed(314)
#train PLS model
model<-train(
Trend00~., data=train.data, method="pls",
scale=TRUE, center=TRUE,
trControl = trainControl("cv", number=5),
tuneLength = 10
)
#view model summary
plot(model)
model$bestTune ## ncomp
## 6 6summary(model$finalModel) ## Data: X dimension: 54 10
## Y dimension: 54 1
## Fit method: oscorespls
## Number of components considered: 6
## TRAINING: % variance explained
## 1 comps 2 comps 3 comps 4 comps 5 comps 6 comps
## X 23.93 41.76 54.29 65.79 75.68 83.41
## .outcome 24.86 35.15 39.30 40.22 40.42 40.49predictions<- model %>% predict(test.data)
data.frame(
RMSE = caret::RMSE(predictions, test.data$Trend00),
Rsquare = caret::R2(predictions, test.data$Trend00)
)## RMSE Rsquare
## 1 3.036188 0.6862633#PLS model performance vs. test data
lm.out<-lm(test.data$Trend00~predictions)
newx<-seq(-23, 23, by=0.05)
conf_int<-stats::predict(lm.out, newdata=data.frame(predictions=newx), interval="confidence", level=0.95)
summary(lm.out)##
## Call:
## lm(formula = test.data$Trend00 ~ predictions)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.224 -1.525 1.014 1.383 4.029
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.1290 0.8387 -1.346 0.21117
## predictions 0.7230 0.1630 4.437 0.00163 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.766 on 9 degrees of freedom
## Multiple R-squared: 0.6863, Adjusted R-squared: 0.6514
## F-statistic: 19.69 on 1 and 9 DF, p-value: 0.001631Plot model performance
#png("PLS_perf_plot25_springGulf.png", width = (6.85039), height = (6.85039), units = 'in', res = 600)
plot(test.data$Trend00~predictions, type='n', xlim=c(-20,20), ylim=c(-20,20), ylab="Observed population trajectory (percent annual change)", xlab="Predicted population trajectory (percent annual change)", pch=19, cex=1.5, col=rgb(158/255, 52/255, 235/255))
polygon(c(rev(newx), newx), c(rev(conf_int[ ,3]), conf_int[ ,2]), col=rgb(158/255, 52/255, 235/255, 0.2), border = NA)
points(test.data$Trend00~predictions, xlim=c(-20,20), ylim=c(-20,20), ylab="Observed population trajectory (percent annual change)", xlab="Predicted population trajectory (percent annual change)", pch=19, cex=1.5, col=rgb(158/255, 52/255, 235/255))
abline(1,1, col="gray", lty=1, lwd=3)
abline(lm(test.data$Trend00~predictions), col=rgb(158/255, 52/255, 235/255), lwd=3, lty=3)
#matlines(newx, conf_int[,2:3], col=rgb(158/255, 52/255, 235/255, 0.5), lty=1, lwd=1)
legend(-20, 20, legend=c("1:1 line", "PLS regression model"), lty=c(1,3), lwd=c(3,3), col=c("gray", rgb(158/255, 52/255, 235/255)), box.lty=0)
text(10,-10, expression(italic('R')^'2'* ' = 0.69'), cex=1)
#dev.off()Plot VIP versus coefficient estimates
#VIP vs. coefficient plot
coefficients<-coef(model$finalModel)
coefficients<-data.frame(coefficients)
coefficients$var<-row.names(coefficients)
coefficients<-coefficients[order(coefficients$var),]
pls.var.imp<-varImp(model$finalModel)
pls.var.imp$var<-rownames(pls.var.imp)
pls.var.imp<-pls.var.imp[order(pls.var.imp$var),]
vip.coef<-data.frame(cbind(coefficients$var, coefficients$.outcome.6.comps, pls.var.imp$Overall))
colnames(vip.coef)<-c("var","coef","vip")
vip.coef$coef<-abs(as.numeric(as.character(vip.coef$coef)))
vip.coef$vip<-as.numeric(as.character(vip.coef$vip))
vip.coef$type<-c(rep("br", 3), rep("mrf", 3), rep("nb", 3), rep("mrf", 1))
vip.coef$color<-c(rep("#35B779",3), rep("#FDE725", 3), rep("#443A83",3), rep("#FDE725", 1))
vip.coef$var<-c("Forest", "Net change forest", "Human footprint",
"Agriculture", "Forest", "Human footprint",
"Forest", "Net change forest", "Human footprint",
"Net forest change")
#VIP vs. COEF plot
#png("VIP_vs_COEF_PLS_plot_25percentile_BrNBMig_SpringPreGulf_NoAgriculture.png", width = (6.85039), height = (6.85039), units = 'in', res = 600)
plot(vip~coef, data=vip.coef, main="Gulf of Mexico, spring", col=color, pch=19, ylab="Variable importance for the projection (VIP)", xlab="Regression coefficient (absolute value)", cex=1.5, ylim=c(0,1.2), xlim=c(0,2.8))
legend(x=0.1, y=1.2, legend=c("Breeding factors", "Migration risk-factors", "Nonbreeding factors"), box.lty=0, pt.bg=c("#35B779", "#FDE725", "#443A83"), pch=21, pt.cex=2, cex=1)
abline(h=0.8, col="gray70", lty=2)
abline(v=1, col="gray70", lty=2)
text(vip~coef, data=vip.coef[c(1),], labels=vip.coef$var[c(1)], cex=0.6, pos=2) # left
text(vip~coef, data=vip.coef[c(3, 4, 7, 8, 10,11,12),], labels=vip.coef$var[c(3,4,7, 8,10, 11,12)], cex=0.6, pos=4) #right
text(vip~coef, data=vip.coef[c(2),], labels=vip.coef$var[c(2)], cex=0.6, pos=3) #up
text(vip~coef, data=vip.coef[c(5,6,9),], labels=vip.coef$var[c(5,6,9)], cex=0.6, pos=1) #bottom
#6 components ; explains 83% x; 40% outcome, R2=0.69
#dev.off()Generalized Linear Modeling
session_info()## ─ Session info ───────────────────────────────────────────────────────────────
## setting value
## version R version 4.2.2 (2022-10-31)
## os macOS Ventura 13.0.1
## system aarch64, darwin20
## ui X11
## language (EN)
## collate en_US.UTF-8
## ctype en_US.UTF-8
## tz America/New_York
## date 2023-01-16
## pandoc 2.19.2 @ /Applications/RStudio.app/Contents/MacOS/quarto/bin/tools/ (via rmarkdown)
##
## ─ Packages ───────────────────────────────────────────────────────────────────
## package * version date (UTC) lib source
## AICcmodavg * 2.3-1 2020-08-26 [1] CRAN (R 4.2.0)
## assertthat 0.2.1 2019-03-21 [1] CRAN (R 4.2.0)
## backports 1.4.1 2021-12-13 [1] CRAN (R 4.2.0)
## broom 1.0.2 2022-12-15 [1] CRAN (R 4.2.0)
## bslib 0.4.1 2022-11-02 [1] CRAN (R 4.2.0)
## cachem 1.0.6 2021-08-19 [1] CRAN (R 4.2.0)
## callr 3.7.3 2022-11-02 [1] CRAN (R 4.2.0)
## caret * 6.0-93 2022-08-09 [1] CRAN (R 4.2.0)
## cellranger 1.1.0 2016-07-27 [1] CRAN (R 4.2.0)
## class 7.3-20 2022-01-16 [1] CRAN (R 4.2.2)
## cli 3.4.1 2022-09-23 [1] CRAN (R 4.2.0)
## codetools 0.2-18 2020-11-04 [1] CRAN (R 4.2.2)
## colorspace 2.0-3 2022-02-21 [1] CRAN (R 4.2.0)
## crayon 1.5.2 2022-09-29 [1] CRAN (R 4.2.0)
## data.table 1.14.6 2022-11-16 [1] CRAN (R 4.2.0)
## DBI 1.1.3 2022-06-18 [1] CRAN (R 4.2.0)
## dbplyr 2.2.1 2022-06-27 [1] CRAN (R 4.2.0)
## devtools * 2.4.5 2022-10-11 [1] CRAN (R 4.2.0)
## digest 0.6.30 2022-10-18 [1] CRAN (R 4.2.0)
## dplyr * 1.0.10 2022-09-01 [1] CRAN (R 4.2.0)
## ellipsis 0.3.2 2021-04-29 [1] CRAN (R 4.2.0)
## evaluate 0.18 2022-11-07 [1] CRAN (R 4.2.0)
## fansi 1.0.3 2022-03-24 [1] CRAN (R 4.2.0)
## fastmap 1.1.0 2021-01-25 [1] CRAN (R 4.2.0)
## forcats * 0.5.2 2022-08-19 [1] CRAN (R 4.2.0)
## foreach 1.5.2 2022-02-02 [1] CRAN (R 4.2.0)
## foreign 0.8-83 2022-09-28 [1] CRAN (R 4.2.2)
## fs 1.5.2 2021-12-08 [1] CRAN (R 4.2.0)
## future 1.30.0 2022-12-16 [1] CRAN (R 4.2.0)
## future.apply 1.10.0 2022-11-05 [1] CRAN (R 4.2.0)
## gargle 1.2.1 2022-09-08 [1] CRAN (R 4.2.0)
## generics 0.1.3 2022-07-05 [1] CRAN (R 4.2.0)
## ggplot2 * 3.4.0 2022-11-04 [1] CRAN (R 4.2.0)
## globals 0.16.2 2022-11-21 [1] CRAN (R 4.2.0)
## glue 1.6.2 2022-02-24 [1] CRAN (R 4.2.0)
## googledrive 2.0.0 2021-07-08 [1] CRAN (R 4.2.0)
## googlesheets4 1.0.1 2022-08-13 [1] CRAN (R 4.2.0)
## gower 1.0.1 2022-12-22 [1] CRAN (R 4.2.0)
## gtable 0.3.1 2022-09-01 [1] CRAN (R 4.2.0)
## hardhat 1.2.0 2022-06-30 [1] CRAN (R 4.2.0)
## haven 2.5.1 2022-08-22 [1] CRAN (R 4.2.0)
## highr 0.9 2021-04-16 [1] CRAN (R 4.2.0)
## hms 1.1.2 2022-08-19 [1] CRAN (R 4.2.0)
## htmltools 0.5.4 2022-12-07 [1] CRAN (R 4.2.2)
## htmlwidgets 1.5.4 2021-09-08 [1] CRAN (R 4.2.0)
## httpuv 1.6.6 2022-09-08 [1] CRAN (R 4.2.0)
## httr 1.4.4 2022-08-17 [1] CRAN (R 4.2.0)
## ipred 0.9-13 2022-06-02 [1] CRAN (R 4.2.0)
## iterators 1.0.14 2022-02-05 [1] CRAN (R 4.2.0)
## jquerylib 0.1.4 2021-04-26 [1] CRAN (R 4.2.0)
## jsonlite 1.8.4 2022-12-06 [1] CRAN (R 4.2.0)
## knitr * 1.41 2022-11-18 [1] CRAN (R 4.2.0)
## later 1.3.0 2021-08-18 [1] CRAN (R 4.2.0)
## lattice * 0.20-45 2021-09-22 [1] CRAN (R 4.2.2)
## lava 1.7.1 2023-01-06 [1] CRAN (R 4.2.2)
## lifecycle 1.0.3 2022-10-07 [1] CRAN (R 4.2.0)
## listenv 0.9.0 2022-12-16 [1] CRAN (R 4.2.0)
## lubridate 1.9.0 2022-11-06 [1] CRAN (R 4.2.0)
## magrittr 2.0.3 2022-03-30 [1] CRAN (R 4.2.0)
## maps * 3.4.1 2022-10-30 [1] CRAN (R 4.2.0)
## maptools * 1.1-6 2022-12-14 [1] CRAN (R 4.2.0)
## MASS 7.3-58.1 2022-08-03 [1] CRAN (R 4.2.2)
## Matrix 1.5-1 2022-09-13 [1] CRAN (R 4.2.2)
## memoise 2.0.1 2021-11-26 [1] CRAN (R 4.2.0)
## mime 0.12 2021-09-28 [1] CRAN (R 4.2.0)
## miniUI 0.1.1.1 2018-05-18 [1] CRAN (R 4.2.0)
## ModelMetrics 1.2.2.2 2020-03-17 [1] CRAN (R 4.2.0)
## modelr 0.1.10 2022-11-11 [1] CRAN (R 4.2.0)
## munsell 0.5.0 2018-06-12 [1] CRAN (R 4.2.0)
## nlme 3.1-160 2022-10-10 [1] CRAN (R 4.2.2)
## nnet 7.3-18 2022-09-28 [1] CRAN (R 4.2.2)
## parallelly 1.33.0 2022-12-14 [1] CRAN (R 4.2.0)
## pbapply 1.7-0 2023-01-13 [1] CRAN (R 4.2.2)
## pillar 1.8.1 2022-08-19 [1] CRAN (R 4.2.0)
## pkgbuild 1.4.0 2022-11-27 [1] CRAN (R 4.2.0)
## pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 4.2.0)
## pkgload 1.3.2 2022-11-16 [1] CRAN (R 4.2.0)
## pls * 2.8-1 2022-07-16 [1] CRAN (R 4.2.0)
## plyr 1.8.8 2022-11-11 [1] CRAN (R 4.2.0)
## prettyunits 1.1.1 2020-01-24 [1] CRAN (R 4.2.0)
## pROC 1.18.0 2021-09-03 [1] CRAN (R 4.2.0)
## processx 3.8.0 2022-10-26 [1] CRAN (R 4.2.0)
## prodlim 2019.11.13 2019-11-17 [1] CRAN (R 4.2.0)
## profvis 0.3.7 2020-11-02 [1] CRAN (R 4.2.0)
## promises 1.2.0.1 2021-02-11 [1] CRAN (R 4.2.0)
## ps 1.7.2 2022-10-26 [1] CRAN (R 4.2.0)
## purrr * 0.3.5 2022-10-06 [1] CRAN (R 4.2.0)
## pwr * 1.3-0 2020-03-17 [1] CRAN (R 4.2.0)
## R6 2.5.1 2021-08-19 [1] CRAN (R 4.2.0)
## raster * 3.6-11 2022-11-28 [1] CRAN (R 4.2.0)
## Rcpp 1.0.9 2022-07-08 [1] CRAN (R 4.2.0)
## readr * 2.1.3 2022-10-01 [1] CRAN (R 4.2.0)
## readxl 1.4.1 2022-08-17 [1] CRAN (R 4.2.0)
## recipes 1.0.3 2022-11-09 [1] CRAN (R 4.2.0)
## remotes 2.4.2 2021-11-30 [1] CRAN (R 4.2.0)
## reprex 2.0.2 2022-08-17 [1] CRAN (R 4.2.0)
## reshape2 1.4.4 2020-04-09 [1] CRAN (R 4.2.0)
## rgdal * 1.6-3 2022-12-14 [1] CRAN (R 4.2.0)
## rgeos * 0.6-1 2022-12-14 [1] CRAN (R 4.2.0)
## rlang 1.0.6 2022-09-24 [1] CRAN (R 4.2.0)
## rmarkdown * 2.18 2022-11-09 [1] CRAN (R 4.2.0)
## rpart 4.1.19 2022-10-21 [1] CRAN (R 4.2.2)
## rstudioapi 0.14 2022-08-22 [1] CRAN (R 4.2.0)
## rvest 1.0.3 2022-08-19 [1] CRAN (R 4.2.0)
## sass 0.4.4 2022-11-24 [1] CRAN (R 4.2.0)
## scales 1.2.1 2022-08-20 [1] CRAN (R 4.2.0)
## sessioninfo 1.2.2 2021-12-06 [1] CRAN (R 4.2.0)
## shiny 1.7.3 2022-10-25 [1] CRAN (R 4.2.0)
## sp * 1.5-1 2022-11-07 [1] CRAN (R 4.2.0)
## stringi 1.7.8 2022-07-11 [1] CRAN (R 4.2.0)
## stringr * 1.5.0 2022-12-02 [1] CRAN (R 4.2.0)
## survival 3.4-0 2022-08-09 [1] CRAN (R 4.2.2)
## terra 1.6-47 2022-12-02 [1] CRAN (R 4.2.0)
## tibble * 3.1.8 2022-07-22 [1] CRAN (R 4.2.0)
## tidyr * 1.2.1 2022-09-08 [1] CRAN (R 4.2.0)
## tidyselect 1.2.0 2022-10-10 [1] CRAN (R 4.2.0)
## tidyverse * 1.3.2 2022-07-18 [1] CRAN (R 4.2.0)
## timechange 0.1.1 2022-11-04 [1] CRAN (R 4.2.0)
## timeDate 4022.108 2023-01-07 [1] CRAN (R 4.2.2)
## tzdb 0.3.0 2022-03-28 [1] CRAN (R 4.2.0)
## unmarked 1.2.5 2022-05-13 [1] CRAN (R 4.2.0)
## urlchecker 1.0.1 2021-11-30 [1] CRAN (R 4.2.0)
## usdm * 1.1-18 2017-06-25 [1] CRAN (R 4.2.0)
## usethis * 2.1.6 2022-05-25 [1] CRAN (R 4.2.0)
## utf8 1.2.2 2021-07-24 [1] CRAN (R 4.2.0)
## vctrs 0.5.1 2022-11-16 [1] CRAN (R 4.2.0)
## VGAM 1.1-7 2022-07-06 [1] CRAN (R 4.2.0)
## withr 2.5.0 2022-03-03 [1] CRAN (R 4.2.0)
## xfun 0.35 2022-11-16 [1] CRAN (R 4.2.0)
## xml2 1.3.3 2021-11-30 [1] CRAN (R 4.2.0)
## xtable 1.8-4 2019-04-21 [1] CRAN (R 4.2.0)
## yaml 2.3.6 2022-10-18 [1] CRAN (R 4.2.0)
##
## [1] /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library
##
## ──────────────────────────────────────────────────────────────────────────────