Blue-winged Teal (Spatula discors)
Spatio-temporal “toy model” of seasonal, annual, and spatial variabilty across North America.
Contact:
John: jmh09r@my.fsu.edu
date() ## [1] "Wed Feb 14 22:14:46 2018"Options
library(knitr)
library(rgl)
opts_knit$set(verbose = TRUE)
knit_hooks$set(webgl = hook_webgl)Load Packages:
suppressMessages(library(INLA))
suppressMessages(library(rgdal))
suppressMessages(library(reshape2))
suppressMessages(library(auk))
suppressMessages(library(raster))
suppressMessages(library(fields))
suppressMessages(library(FactoMineR))
suppressMessages(library(adegenet))
suppressMessages(library(factoextra))
suppressMessages(library(corrplot))
suppressMessages(library(rasterVis))
suppressMessages(library(rgeos))
suppressMessages(library(maptools))
suppressMessages(library(mapproj))
suppressMessages(library(spdep))
suppressMessages(library(ggplot2))
suppressMessages(library(ggthemes))
suppressMessages(library(GISTools))
suppressMessages(library(lattice))
suppressMessages(library(gridExtra))
suppressMessages(library(xtable))
suppressMessages(library(PresenceAbsence))
suppressMessages(library(dismo))
suppressMessages(library(missMDA))
suppressMessages(library(dplyr))
suppressMessages(library(data.table))
suppressMessages(library(readr))
suppressMessages(library(stringr))
source("RScripts.R") #script plots longitude and latitude on ggplotSet Directory
setwd("D:/Patuxent/Patux_level1")load("D:/Patuxent/Patux_level1/Wrk_ST_Feb14_2018.RData")Function to Clean eBird Data
Using eBird data.
Citation:
eBird Basic Dataset. Version: EBD_relNov-2017. Cornell Lab of Ornithology, Ithaca, New York. Nov 2017.
CleanBird = function(x, y) {
df1 = read_delim(x, col_names = TRUE, delim = "\t", quoted_na = FALSE, trim_ws = TRUE)
colnames(df1) = make.names(colnames(df1))
df1 = df1 %>%
filter(is.na(LONGITUDE) == FALSE,
is.na(LATITUDE) == FALSE,
APPROVED == 1) %>%
mutate(Spp = as.factor(SCIENTIFIC.NAME),
Lat = as.numeric(LATITUDE),
Long = as.numeric(LONGITUDE),
Count = as.integer(OBSERVATION.COUNT),
Date = OBSERVATION.DATE) %>%
select(Spp, Long, Lat, Count, Date)
df1 = df1[complete.cases(df1),]
df1 = df1[!is.na(as.numeric(as.character(df1$Long))),]
df1 = df1[!is.na(as.numeric(as.character(df1$Lat))),]
assign(y, df1, envir = .GlobalEnv)
}Clean the Data
Selecting years and designating seasons. Using 2012-2016 for this example.
CleanBird("D:/Patuxent/eBird/ebd_buwtea_relNov_2017/ebd_buwtea_relNov-2017.txt", "BWTE.cln")
BWTE.cln = as.data.frame(BWTE.cln)
BWTE.cln$Year = as.integer(substr(BWTE.cln$Date,0,4))
BWTE.cln$Month = as.integer(format(as.Date(BWTE.cln$Date), "%m"))
BWTE.cln$Day = as.integer(format(as.Date(BWTE.cln$Date), "%d"))
dim(BWTE.cln)[1] #No of Records
BWTE.cln = BWTE.cln %>%
filter(Year >= 2012 & Year <= 2016)
BWTE.cln$Season = ifelse(BWTE.cln$Month == 12 | BWTE.cln$Month == 1, "Winter",
ifelse(BWTE.cln$Month == 3, "Spring",
ifelse(BWTE.cln$Month >= 4 & BWTE.cln$Month <= 8, "Breeding",
ifelse(BWTE.cln$Month >= 9 & BWTE.cln$Month <= 11, "Fall",
ifelse(BWTE.cln$Month == 2 & BWTE.cln$Day <= 15, "Winter",
ifelse(BWTE.cln$Month == 2 & BWTE.cln$Day >= 16, "Spring",
NA))))))
BWTE.cln %>%
group_by(Season) %>%
summarise(no_rows = length(Season))
str(BWTE.cln)Get Boundaries of North America
Used to define model boundaries and to visualize results.
NAmer = readOGR(dsn = "./Arc/NAmer",
layer = "NAmer",
stringsAsFactors = FALSE)
tmp.ras = raster(res=0.25)
extent(tmp.ras) = extent(NAmer)
Domain.r = rasterize(NAmer,
tmp.ras,
field = 0,
background = NA)
Grd.pnt = rasterToPoints(Domain.r, spatial = TRUE)
Grd.pnt@data = Grd.pnt@data %>%
mutate(Long = Grd.pnt@coords[,1],
Lat = Grd.pnt@coords[,2])Convert Observations to Points
Removing duplicates for each unique year and season combination (just to speed thing up).
Removing observations over oceans.
BWTE.cln$uStep = paste(BWTE.cln$Year, BWTE.cln$Season, sep = "")
uStep.lvs = levels(factor(BWTE.cln$uStep))
for(i in 1:length(uStep.lvs)){
Tmp.fr = BWTE.cln %>%
filter(uStep == uStep.lvs[i])
Tmp.fr$dups = duplicated(Tmp.fr[, c("Long","Lat")])
Tmp.fr = Tmp.fr %>%
filter(dups == FALSE)
if(i == 1){BWTE.cln2 = Tmp.fr}
else{BWTE.cln2 = rbind(BWTE.cln2, Tmp.fr)}
}
Obs.pnt = SpatialPointsDataFrame(BWTE.cln2[,2:3], BWTE.cln2)
proj4string(Obs.pnt) = proj4string(NAmer)
Obs.pnt$NAmer = is.na(over(Obs.pnt, NAmer))[,1]
Obs.pnt = subset(Obs.pnt, NAmer == FALSE)Quick Plot
Visual inspection of points.
wmap_df = fortify(NAmer)## Regions defined for each Polygonstmp.df = Obs.pnt@data %>%
select(Long, Lat, Season) %>%
mutate(Season = as.factor(Season))
myCol = c("darkblue", "red", "black", "green")
names(myCol) = levels(tmp.df$Season)
colScale = scale_colour_manual(name = "Season", values = myCol,
labels = c("Breeding", "Fall Migration", "Spring Migration", "Winter"))
ggplot(wmap_df, aes(long,lat, group=group)) +
geom_polygon(fill = "tan", col="black") +
geom_point(data=tmp.df,
aes(Long, Lat,
group=NULL, fill=NULL,
col = tmp.df$Season), size=0.75) +
colScale +
xlab("Longitude") +
ylab("Latitude") +
ggtitle("Blue-winged Teal (2012-2016)") +
scale_x_longitude(xmin=-170, xmax=-50, step=10) +
scale_y_latitude(ymin=15, ymax=85, step=10) +
coord_equal() +
theme(panel.grid.minor = element_blank(),
panel.grid.major = element_blank(),
panel.background = element_blank(),
plot.background = element_blank(),
panel.border = element_blank(),
legend.title = element_text(size=12, face="bold"),
axis.title.x = element_text(size=16, face="bold"),
axis.title.y = element_text(size=16, face="bold"),
plot.title = element_text(size=20, face="bold", hjust = 0.5))
State Boundaries for Figures
Getting individual state boundaries for plotting later results.
States = map("state",
fill = TRUE,
plot = FALSE)
IDs = sapply(strsplit(States$names, ":"),
function(x) x[1])
LL84 = "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
States = map2SpatialPolygons(States, IDs = IDs,
proj4string = CRS(LL84))
#Add a dataframe
pid = sapply(slot(States, "polygons"),
function(x) slot(x, "ID"))
p.df = data.frame(ID=1:length(States),
row.names = pid)
States = SpatialPolygonsDataFrame(States, p.df)Stochastic Partial Differential Equations (SPDE)
Point-process model incorporating INLA.
Construct Mesh
Defining study area using a “nested” mesh and projecting to 3D Cartesian coordinates.
SA.dom2 = gConvexHull(NAmer)
bdry = inla.sp2segment(SA.dom2)
MaxEdge = 3.0
mesh0 = inla.mesh.2d(#loc = Obs.pnt,
boundary = bdry,
cutoff = 0.5,
max.edge = MaxEdge,
min.angle = 25)
mesh1 = inla.mesh.2d(loc=mesh0$loc,
max.edge = MaxEdge*2,
offset=MaxEdge*2,
min.angle = 25)
mesh2 = inla.mesh.2d(loc=mesh1$loc,
max.edge = MaxEdge*4,
offset=MaxEdge*4,
min.angle = 25)
MeshLocs = cbind(mesh2$loc[,1], mesh2$loc[,2]) #Locations of mesh nodes
xyz = as.data.frame(
inla.mesh.map(MeshLocs,
projection = "longlat",
inverse = TRUE))
true.radius.of.earth = 6371
radius.of.earth = 1
mesh.dom = inla.mesh.create(loc = xyz,
#boundary = bdry,
cutoff = 50/true.radius.of.earth, #50
refine=list(max.edge = 2000/true.radius.of.earth, min.angle = 25)) Inspecting Mesh
This is an interactive plot; you can click, drag, and zoom on the mesh.
mesh.dom$n #Number of mesh verticies## [1] 2496plot(mesh.dom, rgl = TRUE, col = "white", main = " ")You must enable Javascript to view this page properly.
Mesh Verticies to Quadrature
Duplicating mesh verticies for each year and season.
dd = as.data.frame(cbind(mesh.dom$loc[,1],
mesh.dom$loc[,2],
mesh.dom$loc[,3]))
dd = as.data.frame(
inla.mesh.map(dd,
projection = "longlat",
inverse = FALSE))
names(dd) = c("Long", "Lat")
Years.v = levels(factor(Obs.pnt$Year))
Season.v = levels(factor(Obs.pnt$Season))
for (i in 1:length(Years.v)){
tmp.dd = dd %>%
mutate(Spp = "Qpnt",
OBS = 0,
Year = Years.v[i])
for(j in 1:length(Season.v)){
tmp.dd = tmp.dd %>%
mutate(Season = Season.v[j])
if(j == 1){Ys.dd = tmp.dd}
else(Ys.dd = rbind(Ys.dd, tmp.dd))}
if (i == 1){Q.df = Ys.dd}
else{Q.df = rbind(Q.df, Ys.dd)}
}Combine Observations and Quadrature
Obs.df = Obs.pnt@data %>%
select(Long, Lat, Spp, Year, Season) %>%
mutate(OBS = 1)
Mod.df = rbind(Obs.df, Q.df)
Mod.df$uStep = paste(Mod.df$Year,
Mod.df$Season,
sep = "")
Mod.pnts = SpatialPointsDataFrame(Mod.df[,1:2], Mod.df)
proj4string(Mod.pnts) = proj4string(Obs.pnt)Environmental Covariates
Climate Data
Bio = getData('worldclim',
var='bio',
res=2.5)ENVIREM Data
#Elev
TWet = raster("./ENVIREM/elev/current_2-5arcmin_topoWet.bil")
TRI = raster("./ENVIREM/elev/current_2-5arcmin_tri.bil")
#Set1
aPET = raster("./ENVIREM/Curr/current_2-5arcmin_annualPET.bil")
AIT = raster("./ENVIREM/Curr/current_2-5arcmin_aridityIndexThornthwaite.bil")
CMI = raster("./ENVIREM/Curr/current_2-5arcmin_climaticMoistureIndex.bil")
Cont = raster("./ENVIREM/Curr/current_2-5arcmin_continentality.bil")
#Set2
EmbQ = raster("./ENVIREM/Curr/current_2-5arcmin_embergerQ.bil")
GD0 = raster("./ENVIREM/Curr/current_2-5arcmin_growingDegDays0.bil")
GD5 = raster("./ENVIREM/Curr/current_2-5arcmin_growingDegDays5.bil")
mTc = raster("./ENVIREM/Curr/current_2-5arcmin_maxTempColdest.bil")
#Set3
mTw = raster("./ENVIREM/Curr/current_2-5arcmin_minTempWarmest.bil")
mC10 = raster("./ENVIREM/Curr/current_2-5arcmin_monthCountByTemp10.bil")
PETcq = raster("./ENVIREM/Curr/current_2-5arcmin_PETColdestQuarter.bil")
PETdq = raster("./ENVIREM/Curr/current_2-5arcmin_PETDriestQuarter.bil")
#Set4
sPET = raster("./ENVIREM/Curr/current_2-5arcmin_PETseasonality.bil")
PETwq = raster("./ENVIREM/Curr/current_2-5arcmin_PETWarmestQuarter.bil")
PETwetq = raster("./ENVIREM/Curr/current_2-5arcmin_PETWettestQuarter.bil")
ThrmI = raster("./ENVIREM/Curr/current_2-5arcmin_thermicityIndex.bil")
Envirem.stk = stack(TWet,TRI,aPET,AIT,CMI,Cont,EmbQ,GD0,GD5,mTc,mTw,mC10,PETcq,PETdq,sPET,PETwq,PETwetq,ThrmI)
names(Envirem.stk) = c("TWet","TRI","aPET","AIT","CMI","Cont","EmbQ","GD0","GD5",
"mTc","mTw","mC10","PETcq","PETdq","sPET","PETwq","PETwetq","ThrmI")Extract Climate and Edaphic Variables
Extrcating for model fitting data (Mod.pnts) and the grid used for prediction and visualization (Grd.pnt).
BioEx = as.data.frame(
extract(Bio,
spTransform(Mod.pnts,
CRS(proj4string(Bio))), method = "simple"))
EnviR = as.data.frame(
extract(Envirem.stk,
spTransform(Mod.pnts,
CRS(proj4string(Envirem.stk))), method = "simple"))
BioEx2 = as.data.frame(
extract(Bio,
spTransform(Grd.pnt,
CRS(proj4string(Bio))), method = "simple"))
EnviR2 = as.data.frame(
extract(Envirem.stk,
spTransform(Grd.pnt,
CRS(proj4string(Envirem.stk))), method = "simple"))Impute, Scale, and Center
BioEx = imputePCA(BioEx, ncp=18)
BioEx = as.data.frame(BioEx$completeObs)
EnviR = imputePCA(EnviR, ncp=17)
EnviR = as.data.frame(EnviR$completeObs)
Cov.df = cbind(BioEx, EnviR)
for(i in 1:ncol(Cov.df)){
Cov.df[,i] = as.numeric(scale(Cov.df[,i], scale = T, center = T))
}Covariate Correlation
CorTab = cor(Cov.df)
corrplot(CorTab)
As above, but for prediction grid
BioEx2 = imputePCA(BioEx2, ncp=18)
BioEx2 = as.data.frame(BioEx2$completeObs)
EnviR2 = imputePCA(EnviR2, ncp=17)
EnviR2 = as.data.frame(EnviR2$completeObs)
Cov2.df = cbind(BioEx2, EnviR2)
for(i in 1:ncol(Cov2.df)){
Cov2.df[,i] = as.numeric(scale(Cov2.df[,i], scale = T, center = T))
}
Grd.pnt@data = cbind(Grd.pnt@data, Cov2.df)Add Covariates to Observations
Mod.pnts@data = cbind(Mod.pnts@data, Cov.df)GlobCover
Getting a land use raster and creating a distance to water/wetlands raster. Also, removing oceans so that the absence of observations over the oceans doesn’t swamp observations over continental surface waters.
#GLC = raster("D:/BulkData_dl/Global_LC/Globcover2009_V2.3_Global_/GLOBCOVER_L4_200901_200912_V2.3.tif")
#GLCc = crop(GLC, extent(NAmer))
#tmp.r = rasterize(NAmer,
# tmp.ras,
# field = 0,
# background = NA)
#GLCc = GLCc + tmp.r #Removing oceans
#writeRaster(GLCc, "./GlobCover.tif")
#160, 170, 180, 210 = Water/wetlands
#GLCc[GLCc<160] = NA
#GLCc[GLCc==190] = NA
#GLCc[GLCc==200] = NA
#GLCc[GLCc>210] = NA
#GLCcdm = distance(GLCc)/1000 #convert to km
#GLCcdm = GLCcdm + tmp.r
#writeRaster(GLCcdm, "./Arc/LC/WetDist/GLCcdm.tif")
GlobCover = raster("./Arc/LC/LClass/GlobCover.tif") #pre-processed above
GLCcdm = raster("./Arc/LC/WetDist/GLCcdm.tif") #pre-processed above
Mod.pnts$LC = as.factor(
as.character(
extract(GlobCover,
spTransform(Mod.pnts,
CRS(proj4string(GlobCover))),
method = "simple")))
Mod.pnts$LC[is.na(Mod.pnts$LC)] = "230"
levels(Mod.pnts$LC)
Mod.pnts$wDs = as.numeric(
extract(GLCcdm,
spTransform(Mod.pnts,
CRS(proj4string(GLCcdm))),
method = "simple"))
range(Mod.pnts$wDs, na.rm=T)
Mod.pnts$wDs[is.na(Mod.pnts$wDs)] = mean(Mod.pnts$wDs, na.rm=T)
Mod.pnts$wDsR = round(Mod.pnts$wDs, 0) #rounding to whole km
range(Mod.pnts$wDsR)
Mod.pnts$wDsR[Mod.pnts$wDsR>400] = 400
range(Mod.pnts$wDsR)
#Prediction Grid
Grd.pnt$LC = as.factor(
as.character(
extract(GlobCover,
spTransform(Grd.pnt,
CRS(proj4string(GlobCover))),
method = "simple")))
Grd.pnt$LC[is.na(Grd.pnt$LC)] = "230"
levels(Grd.pnt$LC)
Grd.pnt$wDs = as.numeric(
extract(GLCcdm,
spTransform(Grd.pnt,
CRS(proj4string(GLCcdm))),
method = "simple"))
range(Grd.pnt$wDs, na.rm=T)
Grd.pnt$wDs[is.na(Grd.pnt$wDs)] = mean(Grd.pnt$wDs, na.rm=T)
Grd.pnt$wDsR = round(Grd.pnt$wDs, 0) #rounding to whole km
Grd.pnt$wDsR[Grd.pnt$wDsR>400] = 400
range(Grd.pnt$wDsR)Human Population Density
To consider sampling bias.
Source:
http://sedac.ciesin.columbia.edu/data/set/gpw-v3-population-density
Pop = raster("./PopDensity.tif")
Mod.pnts$Pop = extract(Pop, spTransform(Mod.pnts, CRS(proj4string(Pop))), method = "simple")/1000 #scaling
Mod.pnts$Pop[is.na(Mod.pnts$Pop)] = mean(Mod.pnts$Pop, na.rm=T)
Grd.pnt$Pop = extract(Pop, spTransform(Grd.pnt, CRS(proj4string(Pop))), method = "simple")/1000
Grd.pnt$Pop[is.na(Grd.pnt$Pop)] = mean(Grd.pnt$Pop, na.rm=T)Distance to Conspecific
Distance to nearest neighbor.
suppressMessages(library(spatstat))
nProj = "+proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0
+x_0=0 +y_0=0 +k=1 +units=km +nadgrids=@null +no_defs"
P.pp = spTransform(Mod.pnts, nProj)
By = as.factor(Mod.pnts$OBS)
P.pp = as.ppp(P.pp)
NN = (nndist(P.pp, by=By, k = 1))
Mod.pnts$nConSp = round(NN[,2]/100,0)
range(Mod.pnts$nConSp) Select Time Steps
This model using an autoregressive term for seasonal correlation.
k = length(levels(factor(Mod.pnts$Season)))
Mod.pnts$Step = as.integer(as.factor(Mod.pnts$Season))Convert to 3D Cartesian Coordinates
Also indexing the A matrix by season.
locs = cbind(Mod.pnts$Long, Mod.pnts$Lat)
locs = inla.mesh.map(locs,
projection = "longlat",
inverse = TRUE)
A.bwte = inla.spde.make.A(mesh.dom,
alpha = 2,
loc=locs,
group = Mod.pnts$Step)Prior for SPDE Model0
PC Prior and creating spatial index for each location+season.
FE.df = Mod.pnts@data
spde0 = inla.spde2.pcmatern(mesh.dom, alpha = 2,
prior.range=c(0.01, 0.5),
prior.sigma=c(1, 0.5),
constr = TRUE)
field1 = inla.spde.make.index("field1",
spde0$n.spde,
n.group=k) Construct Data Stack
Organizing data for model fitting.
FE.lst = list(c(field1,
list(intercept1 = 1)),
list(bio1 = FE.df[,"bio1"],
bio2 = FE.df[,"bio2"],
bio7 = FE.df[,"bio7"],
bio13 = FE.df[,"bio13"],
CMI = FE.df[,"CMI"],
ThrmI = FE.df[,"ThrmI"],
PETdq = FE.df[,"PETdq"],
mC10 = FE.df[,"mC10"],
Cont = FE.df[,"Cont"],
TRI = FE.df[,"TRI"],
wDs = FE.df[,"wDsR"],
nConSp = FE.df[,"nConSp"],
Pop = FE.df[,"Pop"],
Year = FE.df[,"Year"],
Season = FE.df[,"Season"],
Step2 = FE.df[,"Step"],
Step3 = FE.df[,"Step"],
LC = FE.df[,"LC"]))
Stack1 = inla.stack(data = list(PA = FE.df$OBS),
A = list(A.bwte, 1),
effects = FE.lst,
tag = "BWTE.0")Setting priors, constructing formula, and executing model
h.spec = list(theta=list(prior='pccor1', param=c(0, 0.9)))
YearPrior = list(theta=list(prior = "normal", param=c(0, 10)))
LCPrior = list(theta=list(prior = "normal", param=c(0, 10)))
pcprior = list(prec = list(prior="pc.prec", param = c(0.5, 0.01)))
Frm0 = PA ~ -1 + intercept1 +
f(field1, model=spde,
group = field1.group,
control.group=list(model="ar1"),
hyper=h.spec) +
f(Year, model = "iid",
hyper = YearPrior) +
f(LC, model = "iid",
replicate = Step2,
hyper = LCPrior) +
f(wDs, model = "rw1",
hyper = pcprior,
replicate = Step3,
scale.model = TRUE) +
mC10 + bio13 + Cont + TRI + Pop
#theta0 = Model0$internal.summary.hyperpar$mean
theta0 = c(-0.3745091, 3.1833081, 6.5489076, 0.1851724, -1.1987900, -3.4403777)
Model0 = inla(Frm0,
data = inla.stack.data(Stack1, spde=spde0),
family = "binomial",
verbose = TRUE,
control.fixed = list(prec = 0.01,
prec.intercept = 0.001),
control.predictor = list(
A = inla.stack.A(Stack1),
compute = TRUE,
link = 1),
control.mode = list(restart = TRUE, theta = theta0),
control.inla = list(int.strategy = "eb"),
control.results = list(return.marginals.random = TRUE,
return.marginals.predictor = TRUE),
control.compute=list(dic = TRUE, cpo = TRUE, waic = TRUE)) summary(Model0)##
## Call:
## c("inla(formula = Frm0, family = \"binomial\", data = inla.stack.data(Stack1, ", " spde = spde0), verbose = TRUE, control.compute = list(dic = TRUE, ", " cpo = TRUE, waic = TRUE), control.predictor = list(A = inla.stack.A(Stack1), ", " compute = TRUE, link = 1), control.inla = list(int.strategy = \"eb\"), ", " control.results = list(return.marginals.random = TRUE, return.marginals.predictor = TRUE), ", " control.fixed = list(prec = 0.01, prec.intercept = 0.001))" )
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 2.7019 61439.0234 22.1233 61463.8487
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept1 -21.8866 3.6141 -29.0037 -21.8795 -14.8164 -21.8648 0
## mC10 1.8968 0.1166 1.6686 1.8966 2.1263 1.8961 0
## bio13 -0.0627 0.0538 -0.1684 -0.0627 0.0426 -0.0626 0
## Cont -0.2372 0.1571 -0.5482 -0.2363 0.0687 -0.2346 0
## TRI -0.5383 0.0232 -0.5841 -0.5383 -0.4930 -0.5381 0
## Pop 0.2860 0.0799 0.1355 0.2837 0.4494 0.2792 0
##
## Random effects:
## Name Model
## field1 SPDE2 model
## Year IID model
## LC IID model
## wDs RW1 model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant 0.975quant mode
## Range for field1 0.7032 0.1523 0.4642 0.6822 1.0588 0.6397
## Stdev for field1 24.6824 5.4002 16.2461 23.9261 37.3232 22.3968
## GroupRho for field1 0.9971 0.0006 0.9957 0.9971 0.9982 0.9972
## Precision for Year 1.2646 0.4020 0.6310 1.2152 2.1948 1.1185
## Precision for LC 0.3050 0.0461 0.2217 0.3028 0.4022 0.2992
## Precision for wDs 0.0324 0.0049 0.0241 0.0320 0.0432 0.0310
##
## Expected number of effective parameters(std dev): 855.04(0.00)
## Number of equivalent replicates : 167.42
##
## Watanabe-Akaike information criterion (WAIC) ...: 24818.88
## Effective number of parameters .................: 789.16
##
## Marginal log-Likelihood: -14978.67
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computedQuick Check
Correlation of predicted to actual observations.
idat = inla.stack.index(Stack1, "BWTE.0")$data
cor(FE.df$OBS, Model0$summary.linear.predictor$mean[idat])## [1] 0.840957Results
Land Use
Notice that there appears to be seasonal variabilty in land use.
LU.tab = as.data.frame(Model0$summary.random$LC[, c(1:2,4,6)])
LU.tab$Indx = as.factor(rep(c(1:4), each = 23))
names(LU.tab) = c("ID", "Mean", "Q025", "Q975", "Indx")
levels(LU.tab$Indx) = c("Breeding", "Fall", "Spring", "Winter")
SLabs = paste(LU.tab $ID[1:23])
LU.plt = ggplot(LU.tab, aes(x=ID, y=Mean)) +
geom_point(size=2, pch=19, col = "red") +
geom_linerange(aes(ymin=Q025, ymax=Q975), colour="black") +
geom_hline(yintercept = 0,
linetype = "dotted",
colour = "red",
size = 1) +
scale_x_discrete(labels=SLabs) +
facet_wrap(~Indx, scales = "free", ncol=5) +
theme_classic() +
xlab("Land Use Category") +
ylab("Effect (logit)") +
theme(axis.title.y = element_text(face="bold", size=16),
axis.title.x = element_text(face="bold", size=14),
axis.text.y = element_text(face="bold", size=16),
axis.text.x = element_text(face="bold", size=10,
vjust=1, hjust=1, angle=45),
strip.text = element_text(face="bold", size = 20))
LU.plt
Distance to Water
mic.df = as.data.frame(Model0$summary.random$wDs)
names(mic.df) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975")
mic.df$POS = as.factor(rep(1:4, each = 288))
levels(mic.df$POS) = c("Breeding", "Fall", "Spring", "Winter")
ggplot(mic.df, aes(ID, Mean)) +
geom_smooth(col = "black",
linetype= "solid",
method = "loess",
se = FALSE,
lwd = 1) +
geom_smooth(data = mic.df, aes(ID, Q025),
col = "grey",
method = "loess",
se = FALSE,
linetype= "dashed") +
geom_smooth(data = mic.df, aes(ID, Q975),
col = "grey",
method = "loess",
se = FALSE,
linetype= "dashed") +
geom_hline(yintercept = 0,
linetype = "dotted",
col = "red",
size = 1) +
geom_vline(xintercept = 0,
linetype = "dotted",
col = "red",
size = 1) +
xlab("Distance from Water/Wetland (km)") +
ylab("Effect (logit)") +
facet_wrap(~POS, scales = "free", ncol=5) +
theme_classic() +
theme(axis.text=element_text(size=16),
strip.text = element_text(face="bold", size = 20),
axis.title.y = element_text(face="bold", size = 20),
axis.title.x = element_text(face="bold", size = 20, vjust=-2))
Zero Threshold
Distance at which the probabilty of occurence goes to zero. These are pretty close across all seasons.
s0LC = list()
s0LC = split(mic.df[,c(1:2)], mic.df[,9])
b1.df = data.frame(s0LC[1])
names(b1.df) = c("ID", "Mean")
b2.df = data.frame(s0LC[2])
names(b2.df) = c("ID", "Mean")
b3.df = data.frame(s0LC[3])
names(b3.df) = c("ID", "Mean")
b4.df = data.frame(s0LC[4])
names(b4.df) = c("ID", "Mean")
FindZero = function(X.df) {
Cross = which(diff(sign(X.df$Mean))!=0)
mDis = (X.df$ID[Cross] + X.df$ID[Cross+1])/2
return(mDis)
}
###Breeding
FindZero(b1.df) ## [1] 152###Fall
FindZero(b2.df) ## [1] 135###Spring
FindZero(b3.df) ## [1] 133###Winter
FindZero(b4.df) ## [1] 127.5Fixed Effects
Int.df = as.data.frame(Model0$marginals.fix[[1]])
mC10.df = as.data.frame(Model0$marginals.fix[[2]])
bio13.df = as.data.frame(Model0$marginals.fix[[3]])
Cont.df = as.data.frame(Model0$marginals.fix[[4]])
Tri.df = as.data.frame(Model0$marginals.fix[[5]])
Pop.df = as.data.frame(Model0$marginals.fix[[6]])
CI.df = as.data.frame(Model0$summary.fixed)[,c(3,5)]
In.plt = ggplot(Int.df , aes(x, y)) +
geom_line(size = 1, col = "black") +
geom_vline(xintercept = CI.df[1,1], col = "gray") +
geom_vline(xintercept = CI.df[1,2], col = "gray") +
geom_vline(xintercept = 0, col = "red") +
xlab(NULL) +
ylab("Density") +
ggtitle("Intercept") +
theme_classic() +
theme(axis.text=element_text(size=16),
axis.title.y = element_text(size = 20),
axis.title.x = element_text(size = 20, vjust=-2),
plot.title = element_text(size = 20, hjust = 0.5))
mC10.plt = ggplot(mC10.df , aes(x, y)) +
geom_line(size = 1, col = "black") +
geom_vline(xintercept = CI.df[2,1], col = "gray") +
geom_vline(xintercept = CI.df[2,2], col = "gray") +
geom_vline(xintercept = 0, col = "red") +
xlab(NULL) +
ylab(NULL) +
ggtitle(expression("Months Count (>10"*~degree*C*")")) +
theme_classic() +
theme(axis.text=element_text(size=16),
axis.title.y = element_text(size = 20),
axis.title.x = element_text(size = 20, vjust=-2),
plot.title = element_text(size = 20, hjust = 0.5))
bio13.plt = ggplot(bio13.df, aes(x, y)) +
geom_line(size = 1, col = "black") +
geom_vline(xintercept = CI.df[3,1], col = "gray") +
geom_vline(xintercept = CI.df[3,2], col = "gray") +
geom_vline(xintercept = 0, col = "red") +
xlab(NULL) +
ylab("Density") +
ggtitle("Precipitation Seasonality (CoV)") +
theme_classic() +
theme(axis.text=element_text(size=16),
axis.title.y = element_text(size = 20),
axis.title.x = element_text(size = 20, vjust=-2),
plot.title = element_text(size = 20, hjust = 0.5))
Cont.plt = ggplot(Cont.df, aes(x, y)) +
geom_line(size = 1, col = "black") +
geom_vline(xintercept = CI.df[4,1], col = "gray") +
geom_vline(xintercept = CI.df[4,2], col = "gray") +
geom_vline(xintercept = 0, col = "red") +
xlab(NULL) +
ylab(NULL) +
ggtitle(expression("Continentality ("*~degree*C*" )")) +
theme_classic() +
theme(axis.text=element_text(size=16),
axis.title.y = element_text(size = 20),
axis.title.x = element_text(size = 20, vjust=-2),
plot.title = element_text(size = 20, hjust = 0.5))
TRI.plt = ggplot(Tri.df , aes(x, y)) +
geom_line(size = 1, col = "black") +
geom_vline(xintercept = CI.df[5,1], col = "gray") +
geom_vline(xintercept = CI.df[5,2], col = "gray") +
geom_vline(xintercept = 0, col = "red") +
xlab(NULL) +
ylab("Density") +
ggtitle("Terrain Rougness Index") +
theme_classic() +
theme(axis.text=element_text(size=16),
axis.title.y = element_text(size = 20),
axis.title.x = element_text(size = 20, vjust=-2),
plot.title = element_text(size = 20, hjust = 0.5))
Pop.plt = ggplot(Pop.df , aes(x, y)) +
geom_line(size = 1, col = "black") +
geom_vline(xintercept = CI.df[6,1], col = "gray") +
geom_vline(xintercept = CI.df[6,2], col = "gray") +
geom_vline(xintercept = 0, col = "red") +
xlab(NULL) +
ylab("Density") +
ggtitle("Population Density (persons/km^2)") +
theme_classic() +
theme(axis.text=element_text(size=16),
axis.title.y = element_text(size = 20),
axis.title.x = element_text(size = 20, vjust=-2),
plot.title = element_text(size = 20, hjust = 0.5))
grid.arrange(In.plt, mC10.plt, bio13.plt, Cont.plt, TRI.plt, Pop.plt, ncol = 2)
Year
This pattern is interesting… Increasing trend with year? I’m curious to see what a longer period would like.
PhySig.df2 = as.data.frame(Model0$summary.random$Year)
names(PhySig.df2) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975", "mode", "kld")
SE.m3 = ggplot(PhySig.df2, aes(x=ID, y=Mean)) +
geom_point(size=2, pch=19, col = "red") +
geom_linerange(aes(ymin=Q025, ymax=Q975), colour="black") +
geom_hline(yintercept = 0,
linetype = "dotted",
colour = "red",
size = 1) +
theme_classic() +
xlab("Year") +
ylab("Effect (logit)") +
theme(axis.title.y = element_text(face="bold", size=14),
axis.title.x = element_text(face="bold", size=14),
axis.text.y = element_text(face="bold", size=12),
axis.text.x = element_text(face="bold", size=14))
SE.m3
Gaussian Random Fields
Organizing Data for Plots
F.s1 = cbind(Model0$summary.random$field1$mean,
field1$field1.group)
s1xmean = list()
s1xmean = split(F.s1[,1], F.s1[,2])
Pred.pnts = Grd.pnt
GLoc = inla.mesh.map(Pred.pnts@data[,2:3],
projection = "longlat",
inverse = TRUE)
Grd_A = inla.spde.make.A(mesh.dom, loc=GLoc)
for(i in 1:k){
Pred.pnts$m1RE = drop(Grd_A %*% s1xmean[[i]])
tmp.r = rasterize(Pred.pnts,
Domain.r,
"m1RE",
background = NA)
if(i == 1){M0.stk = tmp.r}
else{M0.stk = stack(M0.stk, tmp.r)}
}
names(M0.stk) = Season.v[1:k]Gaussian Random Fields (Seasonal)
Essentially, the spatial density of observations.
rng = seq(-23.83, 18.01, 0.01)
mCols = c("#A50026", "#D73027", "#F46D43", "#FDAE61", "#E0F3F8", "#ABD9E9", "#74ADD1", "#4575B4","#313695")
#brewer.pal(11, "RdBu")
cr = colorRampPalette(c(rev(mCols)),
bias = 0.8, space = "rgb")
RF0osi = levelplot(M0.stk[[1:4]],
margin = FALSE,
layout=c(2, 2),
xlab = " ", ylab = NULL,
col.regions = cr, at = rng,
colorkey = list(labels=list(cex=1.5),
space = "right"),
par.settings = list(axis.line = list(col = "black"),
strip.background = list(col = 'transparent'),
strip.border = list(col = 'transparent')),
scales = list(cex = 1.5)) +
latticeExtra::layer(sp.polygons(NAmer, col = "black", lwd = 1)) +
latticeExtra::layer(sp.polygons(States, col = "black", lwd = 1))
RF0osi
Summing all Components (Prediction)
Getting results from the model object (Model0) is a bit tedious… Here, just plotting the average year effect across all years; but, all years can be plotted by season.
Pred.pnts = Grd.pnt
ModResult = Model0
ModMesh = mesh.dom
pLoc = cbind(Pred.pnts$Long, Pred.pnts$Lat)
pLoc = inla.mesh.map(pLoc,
projection = "longlat",
inverse = TRUE)
Ap = inla.spde.make.A(ModMesh, loc=pLoc)
for(i in 1:k){
Pred.pnts$RF = drop(Ap %*% s1xmean[[i]])
Pred.pnts$Int = ModResult$summary.fix[1,1]
#Get fixed effects
mydf = as.data.frame(ModResult$summary.fixed) #Write results to data frame
#Find closest match for the Water RW1 results
LuTable = ModResult$summary.random$wDs
LuTable$POS = rep(1:4, each = 288)
s0LC = list()
s0LC = split(LuTable[,1:2], LuTable[,9])
Rw.lu = as.data.frame(s0LC[i])[,1]
Pred.pnts$RW = sapply(Pred.pnts$wDs, function(x)which.min(abs(x - Rw.lu)))
Rw.lu = as.data.frame(as.data.frame(s0LC[i])[,2])
names(Rw.lu) = "Mean"
tmp.l = dim(Rw.lu)[1]
Rw.lu$POS = as.integer(seq(1, tmp.l, 1))
Pred.pnts$Ds_lu = as.numeric(with(Rw.lu,
Mean[match(Pred.pnts$RW, POS)]))
#Land Cover
LuTable = as.data.frame(ModResult$summary.random$LC[,1:2])
LuTable$ID = rep(LuTable$ID[1:23], 4)
LuTable$ID = as.factor(LuTable$ID)
LuTable$Splt = rep(c(1:4), each = 23)
s1LC = list()
s1LC = split(LuTable[,1:2], LuTable[,3])
Tmp.LU = as.data.frame(s1LC[i])
names(Tmp.LU) = c("ID", "mean")
Pred.pnts$LCEff = with(Tmp.LU,
mean[match(Pred.pnts$LC,
ID)])
#Year
Pred.pnts$YReff = mean(ModResult$summary.random$Year$mean)
Pred.pnts@data = Pred.pnts@data %>%
mutate(LP = Int + RF + LCEff + Ds_lu + YReff +
mC10*ModResult$summary.fix[2,1] +
bio13*ModResult$summary.fix[3,1] +
Cont*ModResult$summary.fix[4,1] +
TRI*ModResult$summary.fix[5,1] +
Pop*ModResult$summary.fix[6,1],
LP = ifelse(is.na(LP), 0, LP),
LPe = plogis(LP))
tmp.r = rasterize(Pred.pnts,
Domain.r,
"LPe",
background = NA)
names(tmp.r) = Season.v[i]
if(i == 1){M0pF.stk = tmp.r}
else{M0pF.stk = stack(M0pF.stk, tmp.r)}
}Predicted Distributions
rng = seq(0, 1, 0.001)
mCols = brewer.pal(11, "RdYlBu")[-6]
cr0 = rev(colorRampPalette((mCols))(n = 500))
cr = colorRampPalette(c("tan", cr0),
bias = 1, space = "rgb")
M1Jop = levelplot(M0pF.stk,
margin = FALSE,
xlab = " ", ylab = NULL,
col.regions = cr, at = rng,
colorkey = list(labels=list(at=c(0, 0.25, 0.5, 0.75, 1),
labels=c("0%", "25%", "50%", "75%", "100%"), cex=1.5),
labels=list(cex=12),
space = "right"),
par.strip.text = list(fontface='bold', cex=1),
par.settings = list(axis.line = list(col = "black"),
strip.background = list(col = 'transparent'),
strip.border = list(col = 'transparent')),
scales = list(cex = 1.5)) +
latticeExtra::layer(sp.polygons(NAmer, col = "black", lwd = 1)) +
latticeExtra::layer(sp.polygons(States, col = "black", lwd = 1))
M1Jop
Session Info
sessionInfo()## R version 3.4.3 (2017-11-30)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 8.1 x64 (build 9600)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=English_United States.1252
## [2] LC_CTYPE=English_United States.1252
## [3] LC_MONETARY=English_United States.1252
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United States.1252
##
## attached base packages:
## [1] grid stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] bindrcpp_0.2 stringr_1.2.0 readr_1.1.1
## [4] data.table_1.10.4-3 dplyr_0.7.4 missMDA_1.11
## [7] dismo_1.1-4 PresenceAbsence_1.1.9 xtable_1.8-2
## [10] gridExtra_2.3 GISTools_0.7-4 MASS_7.3-47
## [13] ggthemes_3.4.0 spdep_0.7-4 spData_0.2.6.7
## [16] mapproj_1.2-5 maptools_0.9-2 rgeos_0.3-26
## [19] rasterVis_0.41 latticeExtra_0.6-28 RColorBrewer_1.1-2
## [22] lattice_0.20-35 corrplot_0.84 factoextra_1.0.5
## [25] ggplot2_2.2.1 adegenet_2.1.0 ade4_1.7-8
## [28] FactoMineR_1.39 fields_9.0 maps_3.2.0
## [31] spam_2.1-1 dotCall64_0.9-04 raster_2.6-7
## [34] auk_0.1.0 reshape2_1.4.2 rgdal_1.2-16
## [37] INLA_17.08.19 Matrix_1.2-12 sp_1.2-5
## [40] rgl_0.98.1 knitr_1.17
##
## loaded via a namespace (and not attached):
## [1] nlme_3.1-131 gmodels_2.16.2 rprojroot_1.2
## [4] tools_3.4.3 backports_1.1.1 R6_2.2.2
## [7] vegan_2.4-5 rpart_4.1-11 lazyeval_0.2.1
## [10] mgcv_1.8-22 colorspace_1.3-2 nnet_7.3-12
## [13] permute_0.9-4 compiler_3.4.3 mice_2.46.0
## [16] flashClust_1.01-2 expm_0.999-2 labeling_0.3
## [19] scales_0.5.0 mvtnorm_1.0-6 hexbin_1.27.1
## [22] digest_0.6.12 foreign_0.8-69 rmarkdown_1.8
## [25] pkgconfig_2.0.1 htmltools_0.3.6 htmlwidgets_0.9
## [28] rlang_0.1.4 shiny_1.0.5 bindr_0.1
## [31] zoo_1.8-0 jsonlite_1.5 gtools_3.5.0
## [34] magrittr_1.5 leaps_3.0 Rcpp_0.12.14
## [37] munsell_0.4.3 ape_5.0 scatterplot3d_0.3-40
## [40] stringi_1.1.6 plyr_1.8.4 parallel_3.4.3
## [43] gdata_2.18.0 ggrepel_0.7.0 deldir_0.1-14
## [46] splines_3.4.3 hms_0.4.0 igraph_1.1.2
## [49] boot_1.3-20 markdown_0.8 seqinr_3.4-5
## [52] LearnBayes_2.15 glue_1.2.0 evaluate_0.10.1
## [55] httpuv_1.3.5 gtable_0.2.0 assertthat_0.2.0
## [58] mime_0.5 coda_0.19-1 survival_2.41-3
## [61] viridisLite_0.2.0 tibble_1.3.4 cluster_2.0.6