Snowshoe Hare
Non-stationary models with physical barriers
Set Directory
setwd("~/Michigan_State/Sean")Get Observation Data
hare.df = read.csv("./snowshoe.gwr.dataframe8-13.csv",
header = TRUE, sep=",")
dim(hare.df)## [1] 269 7head(hare.df)## Site Hare Lat Lon Forest Snow5yr MaxTemp
## 1 1 1 44.36486 -90.72816 61671 101.000 19.98825
## 2 2 0 44.27223 -90.84052 44749 100.096 19.96475
## 3 3 1 44.45853 -90.82058 62663 109.536 19.82050
## 4 4 0 45.44736 -92.11831 33176 113.450 18.59150
## 5 5 0 45.02406 -91.22783 28451 106.016 19.07525
## 6 6 0 44.21396 -90.01462 43213 105.176 19.84025Convert forest cover from pixel# to km^2
hare.df$PerForest = (hare.df$Forest/87266.46)*100Scale and center covariates
Snow5yr.sobj = scale(hare.df$Snow5yr)
MaxTem.sobj = scale(hare.df$MaxTemp)
PerForest.sobj = scale(hare.df$PerForest)
hare.df$sSnow5yr = as.numeric(Snow5yr.sobj)
hare.df$sMaxTemp = as.numeric(MaxTem.sobj)
hare.df$sPerForest = as.numeric(PerForest.sobj)Quick Peek
plot(density(hare.df$sSnow5yr))
plot(density(hare.df$sMaxTemp)) #interesting, bimodal
plot(density(hare.df$sPerForest))
Convert to Points
Hare.pnt = SpatialPointsDataFrame(hare.df[, c("Lon","Lat")], hare.df)
proj4string(Hare.pnt) = "+proj=longlat +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0"Domain Boundaries
World = map("world",
fill = TRUE,
plot = FALSE)
IDs = sapply(strsplit(World$names, ":"), function(x) x[1])
LL84 = "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
WorldP = map2SpatialPolygons(World, IDs = IDs,
proj4string = CRS(projection(LL84)))
#Add a dataframe
pid = sapply(slot(WorldP, "polygons"),
function(x) slot(x, "ID"))
p.df = data.frame( ID=1:length(WorldP),
row.names = pid)
World = SpatialPolygonsDataFrame(WorldP, p.df)
World = spTransform(World, proj4string(Hare.pnt ))
World = gBuffer(World, width = 0, byid = F)## Warning in gBuffer(World, width = 0, byid = F): Spatial object is not
## projected; GEOS expects planar coordinates#States
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"
StatesP = map2SpatialPolygons(States, IDs = IDs,
proj4string = CRS(projection(LL84)))
#Add a dataframe
pid = sapply(slot(StatesP, "polygons"),
function(x) slot(x, "ID"))
p.df = data.frame( ID=1:length(StatesP),
row.names = pid)
States = SpatialPolygonsDataFrame(StatesP, p.df)
States = spTransform(States, proj4string(Hare.pnt ))
Lakes = readOGR(dsn = "C:/Users/humph173/Documents/Michigan_State/Marten/ArcWork/Lakes",
layer = "Lake_2ks",
stringsAsFactors = FALSE)## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\humph173\Documents\Michigan_State\Marten\ArcWork\Lakes", layer: "Lake_2ks"
## with 1 features
## It has 15 fields
## Integer64 fields read as strings: OBJECTID Id area InPoly_FID SimPgnFlagLakesLL = spTransform(Lakes, proj4string(States))
Ext = c(-93.173415, -81.932841, 41.217132, 47.740649)
Domain = crop(States, Ext)
Domain$Name = str_cap_words(rownames(Domain@data))
Water0 = as(extent(Domain), "SpatialPolygons")
p.df = data.frame(ID=1:length(Water0))
Water0 = SpatialPolygonsDataFrame(Water0, p.df, match.ID = F)
proj4string(Water0) = proj4string(Domain)
Water1 = gDifference(Water0, spTransform(Lakes, proj4string(Water0)))
#Rasterized version
Ras = raster(res = 0.02, ext = extent(Water0),
crs = proj4string(Water0))
Domain.r = rasterize(Water0, Ras,
field = 0,
background = NA)
#Point grid version
Grd.pnt = rasterToPoints(Domain.r, spatial = TRUE)
Grd.pnt@data = Grd.pnt@data %>%
mutate(Long = Grd.pnt@coords[,1],
Lat = Grd.pnt@coords[,2]) %>%
select(-layer)Identify States
Samp.State = over(Hare.pnt, States)
Hare.pnt$State = Samp.State$IDQuick Plot
rng = seq(0, 255, 1)
mCols = brewer.pal(11, "RdYlBu")[-6]
cr0 = rev(colorRampPalette((mCols))(n = 256))
cr = colorRampPalette(c("tan", cr0),
bias = 1, space = "rgb")
MyMatrix = matrix(nrow=7, ncol=2)
rownames(MyMatrix) = rownames(coordinates(Domain))
MyMatrix[,1] = c(-89.07686, -86.11260, -91.91021, -83.82014, -91.31258, -83.46125, -89.51171)
MyMatrix[,2] = c(41.80000, 41.36514, 42.28710, 43.0, 47.5041, 41.61285, 43.63285)
levelplot(Domain.r,
margin = FALSE,
xlab = NULL,
ylab = NULL,
maxpixels = 1e5,
col.regions = cr, at = rng,
colorkey = FALSE, par.strip.text = list(fontface='bold', cex=1.5),
par.settings = list(axis.line = list(col = "black"),
strip.background = list(col = 'transparent'),
strip.border = list(col = 'transparent')),
scales = list(cex = 1.25)) +
latticeExtra::layer(sp.polygons(Domain, col = "black", lwd = 0.5)) +
latticeExtra::layer(sp.polygons(LakesLL, fill = "lightblue", col = "transparent", lwd = 0.5)) +
latticeExtra::layer(sp.polygons(Hare.pnt , col = "red", pch=factor(Hare.pnt$Hare), cex = 1)) +
latticeExtra::layer(sp.text(MyMatrix, txt = Domain$Name,
pos =c(1,3,3,2,2,1,1),
col="black",font=list(face="bold"), cex=1)) +
latticeExtra::layer({SpatialPolygonsRescale(layout.north.arrow(),
offset = c(-83, 45.5),
scale = 2)})
Project to KM
nProj = "+proj=utm +zone=16 +datum=NAD83 +units=km +no_defs +ellps=GRS80 +towgs84=0,0,0"
LakesKM = spTransform(Lakes, nProj)
FocalKM = spTransform(Water1, nProj)
Hare.pntP = spTransform(Hare.pnt, nProj)
Grd.pntP = spTransform(Grd.pnt, nProj)Construct Mesh
max.edge = 15
bound.outer = 100
bdry = inla.sp2segment(FocalKM)
mesh = inla.mesh.2d(boundary = bdry,
loc = Hare.pntP,
max.edge = c(1, 5)*max.edge,
cutoff = 15,
min.angle = 23,
offset = c(max.edge, bound.outer))mesh$n## [1] 1926plot(mesh, lwd=0.5)
Define Barriers
tl = length(mesh$graph$tv[,1])
posTri = matrix(0, tl, 2)
for (t in 1:tl){
temp = mesh$loc[mesh$graph$tv[t, ], ]
posTri[t,] = colMeans(temp)[c(1,2)]
}
posTri = SpatialPoints(posTri)
proj4string(posTri) = proj4string(FocalKM)
normal = over(FocalKM, posTri, returnList=T)
normal = unlist(normal)
barrier.triangles = setdiff(1:tl, normal)
poly.barrier = inla.barrier.polygon(mesh, barrier.triangles)Plot Barrier
plot(mesh, main="Mesh and Omega")
plot(poly.barrier, add=T, col='lightblue')
plot(mesh, add=T)
Proximity to Nearest Hare Detection
suppressMessages(library(spatstat))
By = as.factor(Hare.pntP$Hare)
P.pp = as.ppp(Hare.pntP)
NN = (nndist(P.pp, by=By, k = 1))
Hare.pntP$nHare = round(NN[,"1"]/10, 2)
range(Hare.pntP$nHare)## [1] 0.17 6.42Spatial Priors and Indicies
locs = cbind(Hare.pntP@coords[,1],Hare.pntP@coords[,2])
A.obs = inla.spde.make.A(mesh,
alpha = 2,
loc=locs)
A.forest = A.snow = A.temp = A.obs
barrier.model = inla.barrier.pcmatern(mesh,
barrier.triangles = barrier.triangles,
prior.range = c(1000, 0.5),
prior.sigma = c(1, 0.01))
spde = inla.spde2.pcmatern(mesh,
prior.range = c(1000, 0.5),
prior.sigma = c(1, 0.01))
field.obs = inla.spde.make.index("field.obs", spde$n.spde)
field.forest = inla.spde.make.index("field.forest", spde$n.spde)
field.forest.c = inla.spde.make.index("field.forest.c", spde$n.spde)
field.snow = inla.spde.make.index("field.snow", spde$n.spde)
field.snow.c = inla.spde.make.index("field.snow.c", spde$n.spde)
field.temp = inla.spde.make.index("field.temp", spde$n.spde)
field.temp.c = inla.spde.make.index("field.temp.c", spde$n.spde) Organize Data for Fitting
FE.df = Hare.pntP@data
FE.lst = list(c(field.forest,
list(intercept1 = 1)),
list(LocX = 1:nrow(FE.df),
StateX = FE.df[,"State"]))
forest.stk = inla.stack(data = list(Y = cbind(FE.df$sPerForest, NA, NA, NA)),
A = list(A.forest, 1),
effects = FE.lst,
tag = "for.0")
FE.lst = list(c(field.snow,
list(intercept2 = 1)),
list(LocX = 1:nrow(FE.df),
StateX = FE.df[,"State"]))
snow.stk = inla.stack(data = list(Y = cbind(NA, FE.df$sSnow5yr, NA, NA)),
A = list(A.snow, 1),
effects = FE.lst,
tag = "snw.0")
FE.lst = list(c(field.temp,
list(intercept3 = 1)),
list(LocX = 1:nrow(FE.df),
StateX = FE.df[,"State"]))
temp.stk = inla.stack(data = list(Y = cbind(NA, NA, FE.df$sMaxTemp, NA)),
A = list(A.temp, 1),
effects = FE.lst,
tag = "tmp.0")
FE.lst = list(c(field.obs,
field.forest.c,
field.snow.c,
field.temp.c,
list(intercept4 = 1)),
list(Loc = 1:nrow(FE.df),
State = FE.df[,"State"],
nHare = FE.df[,"nHare"]))
obs.stk = inla.stack(data = list(Y = cbind(NA, NA, NA, FE.df$Hare)),
A = list(A.obs, 1),
effects = FE.lst,
tag = "obs.0")
#For individual runs
Base.stk = inla.stack(data = list(PA = FE.df$Hare),
A = list(A.obs, 1),
effects = FE.lst,
tag = "base.0")
#For linear covariate comparison
FE.lst = list(c(field.obs,
list(intercept4 = 1)),
list(Snow = FE.df[,"sSnow5yr"],
Temp = FE.df[,"sMaxTemp"],
Forest = FE.df[,"sPerForest"]))
NVC.stk = inla.stack(data = list(PA = FE.df$Hare),
A = list(A.obs, 1),
effects = FE.lst,
tag = "lc.0")Combine Data Stacks
Joint.stk = inla.stack(forest.stk, snow.stk, temp.stk, obs.stk)Stationary (No Barriers Model)
Frm0 = PA ~ -1 + intercept4 +
f(field.obs, model=spde)
#theta0 = Model0$internal.summary.hyperpar$mean
theta0 = c(5.18945875, 0.06783803)
Model0 = inla(Frm0,
data = inla.stack.data(Base.stk),
family = "binomial",
verbose = TRUE,
control.fixed = list(prec = 0.001,
prec.intercept = 0.0001),
control.predictor = list(
A = inla.stack.A(Base.stk),
compute = TRUE,
link = 1),
control.mode = list(restart = TRUE, theta = theta0),
control.inla = list(strategy="gaussian",
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(Base.stk), ", " verbose = TRUE, control.compute = list(dic = TRUE, cpo = TRUE, ", " waic = TRUE), control.predictor = list(A = inla.stack.A(Base.stk), ", " compute = TRUE, link = 1), control.inla = list(strategy = \"gaussian\", ", " int.strategy = \"eb\"), control.results = list(return.marginals.random = TRUE, ", " return.marginals.predictor = TRUE), control.fixed = list(prec = 0.001, ", " prec.intercept = 1e-04), control.mode = list(restart = TRUE, ", " theta = theta0))")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 4.7277 4.3860 0.1117 9.2254
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept4 -0.9168 0.7303 -2.3506 -0.9168 0.5157 -0.9168 0
##
## Random effects:
## Name Model
## field.obs SPDE2 model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant 0.975quant
## Range for field.obs 395.987 208.1382 142.2926 347.411 933.359
## Stdev for field.obs 1.252 0.4236 0.5947 1.197 2.242
## mode
## Range for field.obs 272.599
## Stdev for field.obs 1.088
##
## Expected number of effective parameters(std dev): 11.90(0.00)
## Number of equivalent replicates : 22.60
##
## Deviance Information Criterion (DIC) ...............: 309.45
## Deviance Information Criterion (DIC, saturated) ....: 309.45
## Effective number of parameters .....................: 11.73
##
## Watanabe-Akaike information criterion (WAIC) ...: 308.56
## Effective number of parameters .................: 10.35
##
## Marginal log-Likelihood: -170.43
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computedBarrier Model
Frm1 = PA ~ -1 + intercept4 +
f(field.obs, model=barrier.model)
#theta1 = Model1$internal.summary.hyperpar$mean
theta1 = c(0.0158923, 4.9743292)
Model1 = inla(Frm1,
data = inla.stack.data(Base.stk),
family = "binomial",
verbose = TRUE,
control.fixed = list(prec = 0.001,
prec.intercept = 0.0001),
control.predictor = list(
A = inla.stack.A(Base.stk),
compute = TRUE,
link = 1),
control.mode = list(restart = TRUE, theta = theta1),
control.inla = list(strategy="gaussian",
int.strategy = "eb"),
control.results = list(return.marginals.random = TRUE,
return.marginals.predictor = TRUE),
control.compute=list(dic = TRUE, cpo = TRUE, waic = TRUE))summary(Model1)##
## Call:
## c("inla(formula = Frm1, family = \"binomial\", data = inla.stack.data(Base.stk), ", " verbose = TRUE, control.compute = list(dic = TRUE, cpo = TRUE, ", " waic = TRUE), control.predictor = list(A = inla.stack.A(Base.stk), ", " compute = TRUE, link = 1), control.inla = list(strategy = \"gaussian\", ", " int.strategy = \"eb\"), control.results = list(return.marginals.random = TRUE, ", " return.marginals.predictor = TRUE), control.fixed = list(prec = 0.001, ", " prec.intercept = 1e-04), control.mode = list(restart = TRUE, ", " theta = theta1))")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 0.9025 56.1582 0.1007 57.1615
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept4 -0.8921 0.6428 -2.1542 -0.8921 0.3689 -0.8921 0
##
## Random effects:
## Name Model
## field.obs RGeneric2
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant 0.975quant mode
## Theta1 for field.obs 0.0349 0.3581 -0.7023 0.0504 0.7018 0.1058
## Theta2 for field.obs 5.6634 0.4473 4.8319 5.6432 6.5888 5.5703
##
## Expected number of effective parameters(std dev): 12.11(0.00)
## Number of equivalent replicates : 22.22
##
## Deviance Information Criterion (DIC) ...............: 311.49
## Deviance Information Criterion (DIC, saturated) ....: 311.49
## Effective number of parameters .....................: 11.88
##
## Watanabe-Akaike information criterion (WAIC) ...: 310.85
## Effective number of parameters .................: 10.66
##
## Marginal log-Likelihood: -171.05
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computedCompare Base Models
Before adding covariates.
Pred.pnts = Grd.pntP
pLoc = cbind(Pred.pnts@coords[,1], Pred.pnts@coords[,2])
Ap = inla.spde.make.A(mesh, loc=pLoc)
Pred.pnts$M0rf = drop(Ap %*% Model0$summary.random$field.obs$mean)
Pred.pnts$Int0 = Model0$summary.fix[1,1]
Pred.pnts$Pred0 = plogis(Pred.pnts$M0rf + Pred.pnts$Int0)
Pred.pnts$M1rf = drop(Ap %*% Model1$summary.random$field.obs$mean)
Pred.pnts$Int1 = Model1$summary.fix[1,1]
Pred.pnts$Pred1 = plogis(Pred.pnts$M1rf + Pred.pnts$Int1)
M0rf.r = rasterize(spTransform(
Pred.pnts,
CRS(proj4string(Domain.r))),
Domain.r,
"M0rf",
background = NA)
M0est.r = rasterize(spTransform(
Pred.pnts,
CRS(proj4string(Domain.r))),
Domain.r,
"Pred0",
background = NA)
M1rf.r = rasterize(spTransform(
Pred.pnts,
CRS(proj4string(Domain.r))),
Domain.r,
"M1rf",
background = NA)
M1est.r = rasterize(spTransform(
Pred.pnts,
CRS(proj4string(Domain.r))),
Domain.r,
"Pred1",
background = NA)Compare Spatial Fields
RF.stk = stack(M0rf.r, M1rf.r)
names(RF.stk) = c("No.Barrier", "Barrier")
rng = seq(-1.42, 1.79, 0.01)
mCols = c("#A50026", "#D73027", "#F46D43", "#FDAE61", "#E0F3F8", "#ABD9E9", "#74ADD1", "#4575B4","#313695")
cr = rev(colorRampPalette(mCols)(n = 500))
cr = colorRampPalette(cr,
bias = 1, space = "rgb")
levelplot(RF.stk,
layout = c(2,1),
margin = FALSE,
xlab = "Random Field Density",
ylab = NULL,
maxpixels = 1e5,
col.regions = cr, at = rng,
colorkey = list(labels=list(cex=1.5),
space = "bottom"),
par.strip.text = list(fontface='bold', cex=1.5),
par.settings = list(axis.line = list(col = "black"),
strip.background = list(col = 'transparent'),
strip.border = list(col = 'transparent')),
scales = list(cex = 1)) +
latticeExtra::layer(sp.polygons(Domain, col = "black", lwd = 0.5)) +
latticeExtra::layer(sp.polygons(Hare.pnt , col = "black", pch=2, cex = 0.25)) +
latticeExtra::layer(sp.text(MyMatrix, txt = Domain$Name,
pos =c(1,3,3,2,2,1,1),
col="black",font=list(face="bold"), cex=1)) +
latticeExtra::layer({SpatialPolygonsRescale(layout.north.arrow(),
offset = c(-83, 45.5),
scale = 2)})
Estimated Occurrence
Mest.stk = stack(M0est.r,M1est.r)
names(Mest.stk) = c("No.Barrier", "Barrier")
rng = seq(0, 1, 0.001)
mCols = c("#A50026", "#D73027", "#F46D43", "#FDAE61", "#E0F3F8", "#ABD9E9", "#74ADD1", "#4575B4","#313695")
cr = rev(colorRampPalette(mCols)(n = 500))
cr = colorRampPalette(cr,
bias = 1, space = "rgb")
levelplot(Mest.stk,
layout = c(2,1),
margin = FALSE,
xlab = "Occurrence Probabilty",
ylab = NULL,
maxpixels = 1e5,
col.regions = cr, at = rng,
colorkey = list(labels=list(cex=1.5),
space = "bottom"),
par.strip.text = list(fontface='bold', cex=1.5),
par.settings = list(axis.line = list(col = "black"),
strip.background = list(col = 'transparent'),
strip.border = list(col = 'transparent')),
scales = list(cex = 1)) +
latticeExtra::layer(sp.polygons(Domain, col = "black", lwd = 0.5)) +
latticeExtra::layer(sp.polygons(Hare.pnt , col = "black", pch=2, cex = 0.25)) +
latticeExtra::layer(sp.text(MyMatrix, txt = Domain$Name,
pos =c(1,3,3,2,2,1,1),
col="black",font=list(face="bold"), cex=1)) +
latticeExtra::layer({SpatialPolygonsRescale(layout.north.arrow(),
offset = c(-83, 45.5),
scale = 2)})
Full Joint Models
Formula (no barriers)
Joint.frm = Y ~ -1 + intercept1 +
intercept2 +
intercept3 +
intercept4 +
f(field.forest,
model=spde) +
f(field.forest.c,
copy = "field.forest",
fixed = FALSE) +
f(field.snow,
model=spde) +
f(field.snow.c,
copy = "field.snow",
fixed = FALSE) +
f(field.temp,
model=spde) +
f(field.temp.c,
copy = "field.temp",
fixed = FALSE) +
f(field.obs,
model=spde)Execute Model
thetaJ = c(1.23558424, 2.82425785, 2.70605843, 4.21214628, -0.09113797, 5.63936215, 0.43178833,
5.12434219, 0.41778903, 4.97759354, -0.36802809, 1.00575593, 0.78103767, -0.07746242)
thetaJ = c(1.23558424, 2.82425785, 2.70605843, 4.21214628, -0.09113797, 5.63936215, 0.43178833,
5.12434219, 0.41778903, 4.97759354, -0.36802809, 1, 1.00575593, 0.78103767, -0.07746242)
ModelJ = inla(Joint.frm,
data = inla.stack.data(Joint.stk,
spde=spde,
barrier.model=barrier.model),
family = c(rep("gaussian",3), "binomial"),
verbose = TRUE,
control.fixed = list(prec = 0.001,
prec.intercept = 0.0001),
control.predictor = list(
A = inla.stack.A(Joint.stk),
compute = TRUE,
link = 1),
control.mode = list(restart = TRUE, theta = thetaJ),
control.inla = list(strategy = "gaussian",
int.strategy = "eb"),
control.results = list(return.marginals.random = TRUE,
return.marginals.predictor = TRUE),
control.compute=list(dic = TRUE, cpo = TRUE, waic = TRUE)) summary(ModelJ)##
## Call:
## c("inla(formula = Joint.frm, family = c(rep(\"gaussian\", 3), \"binomial\"), ", " data = inla.stack.data(Joint.stk, spde = spde, barrier.model = barrier.model), ", " verbose = TRUE, control.compute = list(dic = TRUE, cpo = TRUE, ", " waic = TRUE), control.predictor = list(A = inla.stack.A(Joint.stk), ", " compute = TRUE, link = 1), control.inla = list(strategy = \"gaussian\", ", " int.strategy = \"eb\"), control.results = list(return.marginals.random = TRUE, ", " return.marginals.predictor = TRUE), control.fixed = list(prec = 0.001, ", " prec.intercept = 1e-04), control.mode = list(restart = TRUE, ", " theta = thetaJ))")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 16.7221 213.1033 0.9339 230.7593
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept1 0.0618 0.3135 -0.5538 0.0618 0.6769 0.0618 0
## intercept2 -0.5382 0.8782 -2.2624 -0.5383 1.1845 -0.5382 0
## intercept3 -0.3564 0.7045 -1.7396 -0.3564 1.0257 -0.3564 0
## intercept4 -1.4095 0.8948 -3.1662 -1.4095 0.3458 -1.4095 0
##
## Random effects:
## Name Model
## field.forest SPDE2 model
## field.snow SPDE2 model
## field.temp SPDE2 model
## field.obs SPDE2 model
## field.forest.c Copy
## field.snow.c Copy
## field.temp.c Copy
##
## Model hyperparameters:
## mean sd 0.025quant
## Precision for the Gaussian observations 3.2951 0.4010 2.5731
## Precision for the Gaussian observations[2] 16.8480 2.0302 13.1752
## Precision for the Gaussian observations[3] 14.9655 1.8952 11.5373
## Range for field.forest 106.3074 25.4186 66.4175
## Stdev for field.forest 1.0986 0.1703 0.8100
## Range for field.snow 325.9284 79.6747 198.2785
## Stdev for field.snow 1.7199 0.3456 1.1451
## Range for field.temp 204.8080 46.1702 129.4749
## Stdev for field.temp 1.7475 0.3219 1.2045
## Range for field.obs 371.5285 194.1785 132.9072
## Stdev for field.obs 0.8278 0.3598 0.3125
## Beta for field.forest.c 0.9713 0.2160 0.5504
## Beta for field.snow.c 0.7604 0.2184 0.3211
## Beta for field.temp.c -0.0834 0.1888 -0.4609
## 0.5quant 0.975quant mode
## Precision for the Gaussian observations 3.2724 4.1498 3.2293
## Precision for the Gaussian observations[2] 16.7396 21.1592 16.5391
## Precision for the Gaussian observations[3] 14.8657 18.9859 14.6876
## Range for field.forest 102.8422 165.4906 96.0840
## Stdev for field.forest 1.0824 1.4760 1.0482
## Range for field.snow 316.0507 510.1317 297.1633
## Stdev for field.snow 1.6837 2.5007 1.6129
## Range for field.temp 199.5428 310.3377 189.3738
## Stdev for field.temp 1.7164 2.4674 1.6550
## Range for field.obs 326.8000 870.2437 257.0620
## Stdev for field.obs 0.7680 1.6997 0.6494
## Beta for field.forest.c 0.9701 1.3984 0.9658
## Beta for field.snow.c 0.7649 1.1783 0.7811
## Beta for field.temp.c -0.0808 0.2816 -0.0716
##
## Expected number of effective parameters(std dev): 288.95(0.00)
## Number of equivalent replicates : 3.724
##
## Deviance Information Criterion (DIC) ...............: 1054.06
## Deviance Information Criterion (DIC, saturated) ....: 1376.40
## Effective number of parameters .....................: 288.82
##
## Watanabe-Akaike information criterion (WAIC) ...: 1050.33
## Effective number of parameters .................: 229.18
##
## Marginal log-Likelihood: -836.80
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computedFormula (with barriers)
JointB.frm = Y ~ -1 + intercept1 +
intercept2 +
intercept3 +
intercept4 +
f(field.forest,
model=barrier.model) +
f(field.forest.c,
copy = "field.forest",
fixed = FALSE) +
f(field.snow,
model=barrier.model) +
f(field.snow.c,
copy = "field.snow",
fixed = FALSE) +
f(field.temp,
model=barrier.model) +
f(field.temp.c,
copy = "field.temp",
fixed = FALSE) +
f(field.obs,
model=barrier.model) Execute Model
thetaJB = c(1.23051826, 2.82031096, 2.72260259, -0.09904693, 4.11815429, 0.37785087, 5.45020805, 0.29230560,
4.81297640, -0.39232276, 4.84263158, 0.42749296, 1.06051386, 0.85966753, -0.0562059)
ModelJB = inla(JointB.frm,
data = inla.stack.data(Joint.stk,
spde=spde,
barrier.model=barrier.model),
family = c(rep("gaussian",3), "binomial"),
verbose = TRUE,
control.fixed = list(prec = 0.001,
prec.intercept = 0.0001),
control.predictor = list(
A = inla.stack.A(Joint.stk),
compute = TRUE,
link = 1),
control.mode = list(restart = TRUE, theta = thetaJB),
control.inla = list(strategy = "gaussian",
int.strategy = "eb"),
control.results = list(return.marginals.random = TRUE,
return.marginals.predictor = TRUE),
control.compute=list(dic = TRUE, cpo = TRUE, waic = TRUE)) summary(ModelJB)##
## Call:
## c("inla(formula = JointB.frm, family = c(rep(\"gaussian\", 3), \"binomial\"), ", " data = inla.stack.data(Joint.stk, spde = spde, barrier.model = barrier.model), ", " verbose = TRUE, control.compute = list(dic = TRUE, cpo = TRUE, ", " waic = TRUE), control.predictor = list(A = inla.stack.A(Joint.stk), ", " compute = TRUE, link = 1), control.inla = list(strategy = \"gaussian\", ", " int.strategy = \"eb\"), control.results = list(return.marginals.random = TRUE, ", " return.marginals.predictor = TRUE), control.fixed = list(prec = 0.001, ", " prec.intercept = 1e-04), control.mode = list(restart = TRUE, ", " theta = thetaJB), num.threads = 4)")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 2.8251 305.3424 1.3818 309.5493
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept1 0.0326 0.2949 -0.5465 0.0326 0.6112 0.0326 0
## intercept2 -0.7647 0.8134 -2.3616 -0.7647 0.8309 -0.7647 0
## intercept3 -0.1372 0.5524 -1.2218 -0.1372 0.9464 -0.1372 0
## intercept4 -1.4860 0.8645 -3.1833 -1.4861 0.2098 -1.4860 0
##
## Random effects:
## Name Model
## field.forest RGeneric2
## field.snow RGeneric2
## field.temp RGeneric2
## field.obs RGeneric2
## field.forest.c Copy
## field.snow.c Copy
## field.temp.c Copy
##
## Model hyperparameters:
## mean sd 0.025quant
## Precision for the Gaussian observations 3.2979 0.4003 2.5714
## Precision for the Gaussian observations[2] 16.7002 2.0264 13.0456
## Precision for the Gaussian observations[3] 15.2272 1.9372 11.7322
## Theta1 for field.forest 0.0710 0.1521 -0.2145
## Theta2 for field.forest 4.5148 0.2266 4.0961
## Theta1 for field.snow 0.4733 0.1883 0.1125
## Theta2 for field.snow 5.5846 0.2297 5.1453
## Theta1 for field.temp 0.4285 0.1677 0.1145
## Theta2 for field.temp 5.0045 0.2048 4.6217
## Theta1 for field.obs -0.2731 0.4293 -1.1469
## Theta2 for field.obs 5.6278 0.4518 4.7766
## Beta for field.forest.c 0.9754 0.2174 0.5505
## Beta for field.snow.c 0.7881 0.2248 0.3376
## Beta for field.temp.c -0.0764 0.1911 -0.4582
## 0.5quant 0.975quant mode
## Precision for the Gaussian observations 3.2775 4.1453 3.2409
## Precision for the Gaussian observations[2] 16.5878 21.0141 16.3761
## Precision for the Gaussian observations[3] 15.1212 19.3479 14.9284
## Theta1 for field.forest 0.0654 0.3825 0.0451
## Theta2 for field.forest 4.5037 4.9859 4.4637
## Theta1 for field.snow 0.4693 0.8523 0.4552
## Theta2 for field.snow 5.5796 6.0478 5.5616
## Theta1 for field.temp 0.4218 0.7725 0.3979
## Theta2 for field.temp 4.9962 5.4249 4.9661
## Theta1 for field.obs -0.2603 0.5411 -0.2149
## Theta2 for field.obs 5.6130 6.5472 5.5599
## Beta for field.forest.c 0.9746 1.4043 0.9720
## Beta for field.snow.c 0.7918 1.2208 0.8049
## Beta for field.temp.c -0.0741 0.2940 -0.0657
##
## Expected number of effective parameters(std dev): 292.43(0.00)
## Number of equivalent replicates : 3.68
##
## Deviance Information Criterion (DIC) ...............: 1054.36
## Deviance Information Criterion (DIC, saturated) ....: NULL
## Effective number of parameters .....................: 292.28
##
## Watanabe-Akaike information criterion (WAIC) ...: 1049.86
## Effective number of parameters .................: 231.15
##
## Marginal log-Likelihood: -841.92
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computedLinear Covariates and No Barriers
Simple regression for comparison.
Frm.lc = PA ~ -1 + intercept4 + Forest + Snow + Temp
#theta0 = Model0$internal.summary.hyperpar$mean
theta0 = c(5.18945875, 0.06783803)
Model.lc = inla(Frm.lc,
data = inla.stack.data(NVC.stk),
family = "binomial",
verbose = TRUE,
control.fixed = list(prec = 0.001,
prec.intercept = 0.0001),
control.predictor = list(
A = inla.stack.A(NVC.stk),
compute = TRUE,
link = 1),
#control.mode = list(restart = TRUE, theta = theta0),
control.inla = list(strategy="gaussian",
int.strategy = "eb"),
control.results = list(return.marginals.random = TRUE,
return.marginals.predictor = TRUE),
control.compute=list(dic = TRUE, cpo = TRUE, waic = TRUE)) summary(Model.lc)##
## Call:
## c("inla(formula = Frm.lc, family = \"binomial\", data = inla.stack.data(NVC.stk), ", " verbose = TRUE, control.compute = list(dic = TRUE, cpo = TRUE, ", " waic = TRUE), control.predictor = list(A = inla.stack.A(NVC.stk), ", " compute = TRUE, link = 1), control.inla = list(strategy = \"gaussian\", ", " int.strategy = \"eb\"), control.results = list(return.marginals.random = TRUE, ", " return.marginals.predictor = TRUE), control.fixed = list(prec = 0.001, ", " prec.intercept = 1e-04))")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 1.0034 0.7123 0.2464 1.9622
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept4 -0.9063 0.1600 -1.2206 -0.9063 -0.5924 -0.9063 0
## Forest 0.9657 0.1905 0.5916 0.9657 1.3395 0.9657 0
## Snow 0.1076 0.1875 -0.2605 0.1076 0.4755 0.1076 0
## Temp -0.6602 0.1881 -1.0294 -0.6602 -0.2912 -0.6602 0
##
## The model has no random effects
##
## The model has no hyperparameters
##
## Expected number of effective parameters(std dev): 4.00(0.00)
## Number of equivalent replicates : 67.24
##
## Deviance Information Criterion (DIC) ...............: 285.22
## Deviance Information Criterion (DIC, saturated) ....: 285.22
## Effective number of parameters .....................: 4.005
##
## Watanabe-Akaike information criterion (WAIC) ...: 285.34
## Effective number of parameters .................: 4.032
##
## Marginal log-Likelihood: -159.56
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computedLinear Covariates and Spatial Model
No Barriers.
Frm.spl = PA ~ -1 + intercept4 +
f(field.obs, model=spde) +
Forest + Snow + Temp
theta.spl = c(0.0158923, 4.9743292)
Model.spl = inla(Frm.spl,
data = inla.stack.data(NVC.stk),
family = "binomial",
verbose = TRUE,
control.fixed = list(prec = 0.001,
prec.intercept = 0.0001),
control.predictor = list(
A = inla.stack.A(NVC.stk),
compute = TRUE,
link = 1),
control.mode = list(restart = TRUE, theta = theta.spl),
control.inla = list(strategy="gaussian",
int.strategy = "eb"),
control.results = list(return.marginals.random = TRUE,
return.marginals.predictor = TRUE),
control.compute=list(dic = TRUE, cpo = TRUE, waic = TRUE)) summary(Model.spl)##
## Call:
## c("inla(formula = Frm.spl, family = \"binomial\", data = inla.stack.data(NVC.stk), ", " verbose = TRUE, control.compute = list(dic = TRUE, cpo = TRUE, ", " waic = TRUE), control.predictor = list(A = inla.stack.A(NVC.stk), ", " compute = TRUE, link = 1), control.inla = list(strategy = \"gaussian\", ", " int.strategy = \"eb\"), control.results = list(return.marginals.random = TRUE, ", " return.marginals.predictor = TRUE), control.fixed = list(prec = 0.001, ", " prec.intercept = 1e-04), control.mode = list(restart = TRUE, ", " theta = theta.spl))")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 4.8453 10.4640 0.1536 15.4629
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept4 -1.0406 0.3968 -1.8198 -1.0407 -0.2622 -1.0406 0
## Forest 0.9414 0.2011 0.5466 0.9414 1.3359 0.9414 0
## Snow 0.2716 0.2179 -0.1563 0.2715 0.6991 0.2716 0
## Temp -0.6166 0.2022 -1.0136 -0.6166 -0.2200 -0.6166 0
##
## Random effects:
## Name Model
## field.obs SPDE2 model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant 0.975quant
## Range for field.obs 494.8125 365.6389 143.2176 389.672 1460.748
## Stdev for field.obs 0.5592 0.3265 0.1253 0.497 1.361
## mode
## Range for field.obs 266.7724
## Stdev for field.obs 0.3406
##
## Expected number of effective parameters(std dev): 8.305(0.00)
## Number of equivalent replicates : 32.39
##
## Deviance Information Criterion (DIC) ...............: 278.14
## Deviance Information Criterion (DIC, saturated) ....: 278.14
## Effective number of parameters .....................: 8.277
##
## Watanabe-Akaike information criterion (WAIC) ...: 278.17
## Effective number of parameters .................: 7.998
##
## Marginal log-Likelihood: -160.13
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computedLinear Covariates Spatial Model (with barriers)
Frm.bl = PA ~ -1 + intercept4 +
f(field.obs, model=barrier.model) +
Forest + Snow + Temp
theta.bl = c(0.0158923, 4.9743292)
Model.bl = inla(Frm.bl,
data = inla.stack.data(NVC.stk),
family = "binomial",
verbose = TRUE,
control.fixed = list(prec = 0.001,
prec.intercept = 0.0001),
control.predictor = list(
A = inla.stack.A(NVC.stk),
compute = TRUE,
link = 1),
control.mode = list(restart = TRUE, theta = theta.bl),
control.inla = list(strategy="gaussian",
int.strategy = "eb"),
control.results = list(return.marginals.random = TRUE,
return.marginals.predictor = TRUE),
control.compute=list(dic = TRUE, cpo = TRUE, waic = TRUE)) summary(Model.bl)##
## Call:
## c("inla(formula = Frm.bl, family = \"binomial\", data = inla.stack.data(NVC.stk), ", " verbose = TRUE, control.compute = list(dic = TRUE, cpo = TRUE, ", " waic = TRUE), control.predictor = list(A = inla.stack.A(NVC.stk), ", " compute = TRUE, link = 1), control.inla = list(strategy = \"gaussian\", ", " int.strategy = \"eb\"), control.results = list(return.marginals.random = TRUE, ", " return.marginals.predictor = TRUE), control.fixed = list(prec = 0.001, ", " prec.intercept = 1e-04), control.mode = list(restart = TRUE, ", " theta = theta.spl))")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 1.0312 78.8265 0.1027 79.9604
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept4 -0.9439 0.3600 -1.6506 -0.9439 -0.2378 -0.9439 0
## Forest 0.9388 0.2008 0.5445 0.9388 1.3328 0.9388 0
## Snow 0.2566 0.2181 -0.1716 0.2566 0.6844 0.2566 0
## Temp -0.6222 0.2007 -1.0163 -0.6222 -0.2284 -0.6222 0
##
## Random effects:
## Name Model
## field.obs RGeneric2
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant 0.975quant mode
## Theta1 for field.obs -1.259 0.8209 -3.130 -1.139 0.0013 -0.6588
## Theta2 for field.obs 5.995 0.6375 4.947 5.921 7.4081 5.6424
##
## Expected number of effective parameters(std dev): 8.021(0.00)
## Number of equivalent replicates : 33.54
##
## Deviance Information Criterion (DIC) ...............: 279.54
## Deviance Information Criterion (DIC, saturated) ....: 279.54
## Effective number of parameters .....................: 7.991
##
## Watanabe-Akaike information criterion (WAIC) ...: 279.77
## Effective number of parameters .................: 7.894
##
## Marginal log-Likelihood: -160.42
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computedCompare Models
Models = c("Model0", "Model1",
"ModelJ", "ModelJB",
"Model.lc", "Model.spl",
"Model.bl")
Type = c("No Barriers, No Covariates",
"With Barriers, No Covariates",
"No Barriers, With Covariates",
"With Barriers, With Covariates",
"Non-spatial, Linear Covariates",
"No Barriers, Linear Covariates",
"With Barriers, Linear Covariates")
#Deviance information criteria
DICs = c(Model0$dic$dic,
Model1$dic$dic,
tapply(ModelJ$dic$local.dic, ModelJ$dic$family, sum)[4],
tapply(ModelJB$dic$local.dic, ModelJB$dic$family, sum)[4],
Model.lc$dic$dic, Model.spl$dic$dic, Model.bl$dic$dic)
#Watanabe-Akaike information criteria
WAICs = c(Model0$waic$waic,
Model1$waic$waic,
tapply(ModelJ$waic$local.waic, ModelJ$dic$family, sum)[4],
tapply(ModelJB$waic$local.waic, ModelJB$dic$family, sum)[4],
Model.lc$waic$waic, Model.spl$waic$waic, Model.bl$waic$waic)
#log Conditional Predictive Ordnance
CPO_func = function(x, nF){
tmp.frm = as.data.frame(
cbind(x$cpo$cpo, x$dic$family)) %>%
filter(V2 == nF)
cpo.out = -mean(log(tmp.frm$V1))
return(cpo.out)}
LCPOs = c(-mean(log(Model0$cpo$cpo)),
-mean(log(Model1$cpo$cpo)),
CPO_func(ModelJ, 4),
CPO_func(ModelJB, 4),
-mean(log(Model.lc$cpo$cpo)),
-mean(log(Model.spl$cpo$cpo)),
-mean(log(Model.bl$cpo$cpo)))
Model_mets = as.data.frame(cbind(Model = Models,
DIC = round(DICs, 2),
WAIC = round(WAICs, 2),
lCPO = round(LCPOs, 2),
Type = Type))
row.names(Model_mets) = 1:nrow(Model_mets)kable(Model_mets) %>%
kable_styling("striped", full_width = F) %>%
row_spec(0, font_size = 20) | Model | DIC | WAIC | lCPO | Type |
|---|---|---|---|---|
| Model0 | 309.45 | 308.56 | 0.57 | No Barriers, No Covariates |
| Model1 | 311.49 | 310.85 | 0.58 | With Barriers, No Covariates |
| ModelJ | 285.7 | 285.7 | 0.53 | No Barriers, With Covariates |
| ModelJB | 286.01 | 286.57 | 0.53 | With Barriers, With Covariates |
| Model.lc | 285.22 | 285.34 | 0.53 | Non-spatial, Linear Covariates |
| Model.spl | 278.14 | 278.17 | 0.52 | No Barriers, Linear Covariates |
| Model.bl | 279.54 | 279.77 | 0.52 | With Barriers, Linear Covariates |
Estimation (Fitted Values)
Mod.Coeffs = rbind(
ModelJ$summary.hyperpar[12:14, c(1:3,5)] %>% mutate(Model = "ModelJ"),
ModelJB$summary.hyperpar[12:14, c(1:3,5)] %>% mutate(Model = "ModelJB"),
Model.lc$summary.fixed[2:4, c(1:3,5)] %>% mutate(Model = "Model.lc"),
Model.spl$summary.fixed[2:4, c(1:3,5)] %>% mutate(Model = "Model.spl"),
Model.bl$summary.fixed[2:4, c(1:3,5)] %>% mutate(Model = "Model.bl"))
names(Mod.Coeffs) = c("Mean", "SD", "Q.025", "Q.975", "Model")
Mod.Coeffs$Coeff = rep(c("Forest", "Snow", "Temp"), 5)
Mod.Coeffs = Mod.Coeffs %>%
mutate(Mean = round(Mean, 2),
sd = round(SD, 2),
Q.025 = round(Q.025, 2),
Q.975 = round(Q.975, 2),
CI = paste("(",Q.025,", ",Q.975, ")", sep="")) %>%
select(Model, Coeff, Mean, CI)
Mod.Coeffs = arrange(Mod.Coeffs, Coeff)kable(Mod.Coeffs) %>%
kable_styling("striped", full_width = F) %>%
row_spec(0, font_size = 20) | Model | Coeff | Mean | CI |
|---|---|---|---|
| ModelJ | Forest | 0.97 | (0.55, 1.4) |
| ModelJB | Forest | 0.98 | (0.55, 1.4) |
| Model.lc | Forest | 0.97 | (0.59, 1.34) |
| Model.spl | Forest | 0.94 | (0.55, 1.34) |
| Model.bl | Forest | 0.94 | (0.54, 1.33) |
| ModelJ | Snow | 0.76 | (0.32, 1.18) |
| ModelJB | Snow | 0.79 | (0.34, 1.22) |
| Model.lc | Snow | 0.11 | (-0.26, 0.48) |
| Model.spl | Snow | 0.27 | (-0.16, 0.7) |
| Model.bl | Snow | 0.26 | (-0.17, 0.68) |
| ModelJ | Temp | -0.08 | (-0.46, 0.28) |
| ModelJB | Temp | -0.08 | (-0.46, 0.29) |
| Model.lc | Temp | -0.66 | (-1.03, -0.29) |
| Model.spl | Temp | -0.62 | (-1.01, -0.22) |
| Model.bl | Temp | -0.62 | (-1.02, -0.23) |
Check Fitted Values
Mj = inla.stack.index(Joint.stk, "obs.0")$data
MJ.comp = as.data.frame(
cbind(
plogis(ModelJ$summary.fitted.values$mean[Mj]),
FE.df$Hare))
names(MJ.comp) = c("Fit","Observation")
MJ.mse = mean((MJ.comp$Fit - MJ.comp$Observation)^2)
Mjb = inla.stack.index(Joint.stk, "obs.0")$data
MJB.comp = as.data.frame(
cbind(
plogis(ModelJB$summary.fitted.values$mean[Mjb]),
FE.df$Hare))
names(MJB.comp) = c("Fit","Observation")
MJB.mse = mean((MJB.comp$Fit - MJB.comp$Observation)^2)
Mlc = inla.stack.index(NVC.stk, "lc.0")$data
MLC.comp = as.data.frame(
cbind(
plogis(Model.lc$summary.fitted.values$mean[Mlc]),
FE.df$Hare))
names(MLC.comp) = c("Fit","Observation")
MLC.mse = mean((MLC.comp$Fit - MLC.comp$Observation)^2)
splc = inla.stack.index(NVC.stk, "lc.0")$data
SPLC.comp = as.data.frame(
cbind(
plogis(Model.spl$summary.fitted.values$mean[splc]),
FE.df$Hare))
names(SPLC.comp) = c("Fit","Observation")
SPLC.mse = mean((SPLC.comp$Fit - SPLC.comp$Observation)^2)
blc = inla.stack.index(NVC.stk, "lc.0")$data
BLC.comp = as.data.frame(
cbind(
plogis(Model.bl$summary.fitted.values$mean[blc]),
FE.df$Hare))
names(BLC.comp) = c("Fit","Observation")
BLC.mse = mean((BLC.comp$Fit - BLC.comp$Observation)^2)
Eval.df = as.data.frame(
cbind(
1:nrow(FE.df),
FE.df$Hare,
MJ.comp$Fit,
MJB.comp$Fit,
MLC.comp$Fit,
SPLC.comp$Fit,
BLC.comp$Fit))
names(Eval.df) = c("IDx", "Obs", "JMod", "JBMod", "LCMod", "SPLCMod", "BLCMod")
Mod.out = presence.absence.accuracy(Eval.df, threshold = 0.50)
Mod.out = Mod.out %>% select(model, PCC, sensitivity, specificity, AUC)
Mod.out$MSE = c(MJ.mse, MJB.mse, MLC.mse, SPLC.mse, BLC.mse)kable(Mod.out) %>%
kable_styling("striped", full_width = F) %>%
row_spec(0, font_size = 20) | model | PCC | sensitivity | specificity | AUC | MSE |
|---|---|---|---|---|---|
| JMod | 0.33829 | 1 | 0 | 0.8228794 | 0.2614891 |
| JBMod | 0.33829 | 1 | 0 | 0.8265218 | 0.2613945 |
| LCMod | 0.33829 | 1 | 0 | 0.7920731 | 0.2623348 |
| SPLCMod | 0.33829 | 1 | 0 | 0.8208421 | 0.2595588 |
| BLCMod | 0.33829 | 1 | 0 | 0.8173231 | 0.2599400 |
Get Spatial Random Fields
Pred.pnts = Grd.pntP
ModResult = ModelJB
ModResult2 = ModelJ
pLoc = cbind(Pred.pnts@coords[,1], Pred.pnts@coords[,2])
Ap = inla.spde.make.A(mesh, loc=pLoc)Forest
Pred.pnts$F.rf = drop(Ap %*% ModResult$summary.random$field.forest$mean)
Pred.pnts$F.int = ModResult$summary.fix[1,1]
Pred.pnts$estFor = Pred.pnts$F.rf + Pred.pnts$F.int
Pred.pnts$BestForNat = Pred.pnts$estFor*attr(PerForest.sobj, 'scaled:scale') +
attr(PerForest.sobj, 'scaled:center')
For.Best.r = rasterize(spTransform(
Pred.pnts,
CRS(proj4string(Domain.r))),
Domain.r,
"BestForNat",
background = NA)
Pred.pnts$F.rf = drop(Ap %*% ModResult2$summary.random$field.forest$mean)
Pred.pnts$F.int = ModResult2$summary.fix[1,1]
Pred.pnts$estFor = Pred.pnts$F.rf + Pred.pnts$F.int
Pred.pnts$nBestForNat = Pred.pnts$estFor*attr(PerForest.sobj, 'scaled:scale') +
attr(PerForest.sobj, 'scaled:center')
For.nBest.r = rasterize(spTransform(
Pred.pnts,
CRS(proj4string(Domain.r))),
Domain.r,
"nBestForNat",
background = NA)
Diff.For = For.Best.r-For.nBest.r
CFor.stk = stack(For.nBest.r, For.Best.r)
names(CFor.stk) = c("No.Barrier","Barrier")
rng = seq(0, 100, 0.1)
mCols = brewer.pal(11, "RdYlGn")
cr = colorRampPalette(mCols)(n = 500)
cr = colorRampPalette(cr,
bias = 1, space = "rgb")
levelplot(CFor.stk,
layout = c(2,1),
margin = FALSE,
xlab = NULL,
ylab = NULL,
main = "Forest (%)",
maxpixels = 1e5,
col.regions = cr, at = rng,
colorkey = list(labels=list(cex=1.5),
space = "bottom"),
par.strip.text = list(fontface='bold', cex=1.5),
par.settings = list(axis.line = list(col = "black"),
strip.background = list(col = 'transparent'),
strip.border = list(col = 'transparent')),
scales = list(cex = 1.25)) +
latticeExtra::layer(sp.polygons(Domain, col = "black", lwd = 0.5)) +
latticeExtra::layer(sp.polygons(Hare.pnt , col = "black", pch=2, cex = 0.25)) +
latticeExtra::layer(sp.text(MyMatrix, txt = Domain$Name,
pos =c(1,3,3,2,2,1,1),
col="black",font=list(face="bold"), cex=1)) +
latticeExtra::layer({SpatialPolygonsRescale(layout.north.arrow(),
offset = c(-83, 45.5),
scale = 2)})
rng = seq(-11, 7, 0.1)
mCols = brewer.pal(11, "RdYlGn")
cr = colorRampPalette(mCols)(n = 500)
cr = colorRampPalette(cr,
bias = 0.75, space = "rgb")
levelplot(Diff.For,
margin = FALSE,
xlab = NULL,
ylab = NULL,
main = "Forest Difference \n(Barrier minus NoBarrier)",
maxpixels = 1e5,
col.regions = cr, at = rng,
colorkey = list(labels=list(cex=1.5),
space = "bottom"),
par.strip.text = list(fontface='bold', cex=1.5),
par.settings = list(axis.line = list(col = "black"),
strip.background = list(col = 'transparent'),
strip.border = list(col = 'transparent')),
scales = list(cex = 1.25)) +
latticeExtra::layer(sp.polygons(Domain, col = "black", lwd = 0.5)) +
latticeExtra::layer(sp.polygons(Hare.pnt , col = "black", pch=2, cex = 0.25)) +
latticeExtra::layer(sp.text(MyMatrix, txt = Domain$Name,
pos =c(1,3,3,2,2,1,1),
col="black",font=list(face="bold"), cex=1)) +
latticeExtra::layer({SpatialPolygonsRescale(layout.north.arrow(),
offset = c(-83, 45.5),
scale = 2)})
Snow
Pred.pnts$S.rf = drop(Ap %*% ModResult$summary.random$field.snow$mean)
Pred.pnts$S.int = ModResult$summary.fix[2,1]
Pred.pnts$estSnw = Pred.pnts$S.rf + Pred.pnts$S.int
Pred.pnts$BestSnwNat = Pred.pnts$estSnw*attr(Snow5yr.sobj, 'scaled:scale') +
attr(Snow5yr.sobj, 'scaled:center')
Snw.Best.r = rasterize(spTransform(
Pred.pnts,
CRS(proj4string(Domain.r))),
Domain.r,
"BestSnwNat",
background = NA)
Pred.pnts$S.rf = drop(Ap %*% ModResult2$summary.random$field.snow$mean)
Pred.pnts$S.int = ModResult2$summary.fix[2,1]
Pred.pnts$estSnw = Pred.pnts$S.rf + Pred.pnts$S.int
Pred.pnts$nBestSnwNat = Pred.pnts$estSnw*attr(Snow5yr.sobj, 'scaled:scale') +
attr(Snow5yr.sobj, 'scaled:center')
Snw.nBest.r = rasterize(spTransform(
Pred.pnts,
CRS(proj4string(Domain.r))),
Domain.r,
"nBestSnwNat",
background = NA)
Diff.Snow = Snw.Best.r-Snw.nBest.r
CSnw.stk = stack(Snw.nBest.r, Snw.Best.r)
names(CSnw.stk) = c("No.Barrier","Barrier")
rng = seq(89, 130, 0.1)
mCols = brewer.pal(11, "BrBG")
cr = colorRampPalette(mCols)(n = 500)
cr = colorRampPalette(cr,
bias = 1, space = "rgb")
levelplot(CSnw.stk,
layout = c(2,1),
margin = FALSE,
xlab = NULL,
ylab = NULL,
main = "Snow (Days)",
maxpixels = 1e5,
col.regions = cr, at = rng,
colorkey = list(labels=list(cex=1.5),
space = "bottom"),
par.strip.text = list(fontface='bold', cex=1.5),
par.settings = list(axis.line = list(col = "black"),
strip.background = list(col = 'transparent'),
strip.border = list(col = 'transparent')),
scales = list(cex = 1.25)) +
latticeExtra::layer(sp.polygons(Domain, col = "black", lwd = 0.5)) +
latticeExtra::layer(sp.polygons(Hare.pnt , col = "black", pch=2, cex = 0.25)) +
latticeExtra::layer(sp.text(MyMatrix, txt = Domain$Name,
pos =c(1,3,3,2,2,1,1),
col="black",font=list(face="bold"), cex=1)) +
latticeExtra::layer({SpatialPolygonsRescale(layout.north.arrow(),
offset = c(-83, 45.5),
scale = 2)})
rng = seq(-11.5, 7, 0.1)
mCols = brewer.pal(11, "BrBG")
cr = colorRampPalette(mCols)(n = 500)
cr = colorRampPalette(cr,
bias = 0.75, space = "rgb")
levelplot(Diff.Snow,
margin = FALSE,
xlab = NULL,
ylab = NULL,
main = "Snow Difference \n(Barrier minus NoBarrier)",
maxpixels = 1e5,
col.regions = cr, at = rng,
colorkey = list(labels=list(cex=1.5),
space = "bottom"),
par.strip.text = list(fontface='bold', cex=1.5),
par.settings = list(axis.line = list(col = "black"),
strip.background = list(col = 'transparent'),
strip.border = list(col = 'transparent')),
scales = list(cex = 1.25)) +
latticeExtra::layer(sp.polygons(Domain, col = "black", lwd = 0.5)) +
latticeExtra::layer(sp.polygons(Hare.pnt , col = "black", pch=2, cex = 0.25)) +
latticeExtra::layer(sp.text(MyMatrix, txt = Domain$Name,
pos =c(1,3,3,2,2,1,1),
col="black",font=list(face="bold"), cex=1)) +
latticeExtra::layer({SpatialPolygonsRescale(layout.north.arrow(),
offset = c(-83, 45.5),
scale = 2)})
Max Temperature
Pred.pnts$T.rf = drop(Ap %*% ModResult$summary.random$field.temp$mean)
Pred.pnts$T.int = ModResult$summary.fix[3,1]
Pred.pnts$estTmp = Pred.pnts$T.rf + Pred.pnts$T.int
Pred.pnts$BestTmpNat = Pred.pnts$estTmp*attr(MaxTem.sobj, 'scaled:scale') +
attr(MaxTem.sobj, 'scaled:center')
Tmp.Best.r = rasterize(spTransform(
Pred.pnts,
CRS(proj4string(Domain.r))),
Domain.r,
"BestTmpNat",
background = NA)
Pred.pnts$T.rf = drop(Ap %*% ModResult2$summary.random$field.temp$mean)
Pred.pnts$T.int = ModResult2$summary.fix[3,1]
Pred.pnts$estTmp = Pred.pnts$T.rf + Pred.pnts$T.int
Pred.pnts$nBestTmpNat = Pred.pnts$estTmp*attr(MaxTem.sobj, 'scaled:scale') +
attr(MaxTem.sobj, 'scaled:center')
Tmp.nBest.r = rasterize(spTransform(
Pred.pnts,
CRS(proj4string(Domain.r))),
Domain.r,
"nBestTmpNat",
background = NA)
Diff.temp = Tmp.Best.r-Tmp.nBest.r
CTmp.stk = stack(Tmp.nBest.r, Tmp.Best.r)
names(CTmp.stk) = c("No.Barrier","Barrier")
rng = seq(16, 21, 0.01)
mCols = brewer.pal(11, "RdGy")
cr = rev(colorRampPalette(mCols)(n = 500))
cr = colorRampPalette(cr,
bias = 1, space = "rgb")
levelplot(CTmp.stk,
layout = c(2,1),
margin = FALSE,
xlab = NULL,
ylab = NULL,
main = "Average Max Temperature (Celsius)",
maxpixels = 1e5,
col.regions = cr, at = rng,
colorkey = list(labels=list(cex=1.5),
space = "bottom"),
par.strip.text = list(fontface='bold', cex=1.5),
par.settings = list(axis.line = list(col = "black"),
strip.background = list(col = 'transparent'),
strip.border = list(col = 'transparent')),
scales = list(cex = 1.25)) +
latticeExtra::layer(sp.polygons(Domain, col = "black", lwd = 0.5)) +
latticeExtra::layer(sp.polygons(Hare.pnt , col = "black", pch=2, cex = 0.25)) +
latticeExtra::layer(sp.text(MyMatrix, txt = Domain$Name,
pos =c(1,3,3,2,2,1,1),
col="black",font=list(face="bold"), cex=1)) +
latticeExtra::layer({SpatialPolygonsRescale(layout.north.arrow(),
offset = c(-83, 45.5),
scale = 2)})
rng = seq(-1, 2, 0.01)
mCols = brewer.pal(11, "RdGy")
cr = rev(colorRampPalette(mCols)(n = 500))
cr = colorRampPalette(cr,
bias = 1.8, space = "rgb")
levelplot(Diff.temp,
margin = FALSE,
xlab = NULL,
ylab = NULL,
main = "Temperature Difference \n(Barrier minus NoBarrier)",
maxpixels = 1e5,
col.regions = cr, at = rng,
colorkey = list(labels=list(cex=1.5),
space = "bottom"),
par.strip.text = list(fontface='bold', cex=1.5),
par.settings = list(axis.line = list(col = "black"),
strip.background = list(col = 'transparent'),
strip.border = list(col = 'transparent')),
scales = list(cex = 1.25)) +
latticeExtra::layer(sp.polygons(Domain, col = "black", lwd = 0.5)) +
latticeExtra::layer(sp.polygons(Hare.pnt , col = "black", pch=2, cex = 0.25)) +
latticeExtra::layer(sp.text(MyMatrix, txt = Domain$Name,
pos =c(1,3,3,2,2,1,1),
col="black",font=list(face="bold"), cex=1)) +
latticeExtra::layer({SpatialPolygonsRescale(layout.north.arrow(),
offset = c(-83, 45.5),
scale = 2)})
Hare
Pred.pnts$h.rf = drop(Ap %*% ModResult$summary.random$field.obs$mean)
Pred.pnts$h.Frf = drop(Ap %*% ModResult$summary.random$field.forest.c$mean)
Pred.pnts$h.Srf = drop(Ap %*% ModResult$summary.random$field.snow.c$mean)
Pred.pnts$h.Trf = drop(Ap %*% ModResult$summary.random$field.temp.c$mean)
Pred.pnts$obs.int = ModResult$summary.fix[4,1]
Pred.pnts$F.int = ModResult$summary.fix[1,1]
Pred.pnts$S.int = ModResult$summary.fix[2,1]
Pred.pnts$T.int = ModResult$summary.fix[3,1]
Pred.pnts@data = Pred.pnts@data %>%
mutate(Hare.prob = plogis(
h.rf + h.Frf + h.Srf + h.Trf +
obs.int + F.int + S.int + T.int))
Hare.prob.r = rasterize(spTransform(
Pred.pnts,
CRS(proj4string(Domain.r))),
Domain.r,
"Hare.prob",
background = NA)
##No Barrier estimate
Pred.pnts$h.rf = drop(Ap %*% ModResult2$summary.random$field.obs$mean)
Pred.pnts$h.Frf = drop(Ap %*% ModResult2$summary.random$field.forest.c$mean)
Pred.pnts$h.Srf = drop(Ap %*% ModResult2$summary.random$field.snow.c$mean)
Pred.pnts$h.Trf = drop(Ap %*% ModResult2$summary.random$field.temp.c$mean)
Pred.pnts$obs.int = ModResult2$summary.fix[4,1]
Pred.pnts$F.int = ModResult2$summary.fix[1,1]
Pred.pnts$S.int = ModResult2$summary.fix[2,1]
Pred.pnts$T.int = ModResult2$summary.fix[3,1]
Pred.pnts@data = Pred.pnts@data %>%
mutate(Hare.prob = plogis(
h.rf + h.Frf + h.Srf + h.Trf +
obs.int + F.int + S.int + T.int))
Hare.nBprob.r = rasterize(spTransform(
Pred.pnts,
CRS(proj4string(Domain.r))),
Domain.r,
"Hare.prob",
background = NA)
Diff.est = Hare.prob.r-Hare.nBprob.r
CEst.stk = stack(Hare.nBprob.r, Hare.prob.r)
names(CEst.stk) = c("No.Barrier","Barrier")Compare Estimate
rng = seq(0, 1, 0.01)
mCols = c("#A50026", "#D73027", "#F46D43", "#FDAE61", "#E0F3F8", "#ABD9E9", "#74ADD1", "#4575B4","#313695")
cr = rev(colorRampPalette(mCols)(n = 500))
cr = colorRampPalette(cr,
bias = 1, space = "rgb")
levelplot(CEst.stk,
layout = c(2,1),
margin = FALSE,
xlab = NULL,
ylab = NULL,
main = "Occurrence Probability",
maxpixels = 1e5,
col.regions = cr, at = rng,
colorkey = list(labels=list(cex=1.5),
space = "bottom"),
par.strip.text = list(fontface='bold', cex=1.5),
par.settings = list(axis.line = list(col = "black"),
strip.background = list(col = 'transparent'),
strip.border = list(col = 'transparent')),
scales = list(cex = 1.25)) +
latticeExtra::layer(sp.polygons(Domain, col = "black", lwd = 0.5)) +
latticeExtra::layer(sp.polygons(Hare.pnt , col = "red", pch=factor(Hare.pnt$Hare), cex = 1)) +
latticeExtra::layer(sp.text(MyMatrix, txt = Domain$Name,
pos =c(1,3,3,2,2,1,1),
col="white",font=list(face="bold"), cex=1)) +
latticeExtra::layer({SpatialPolygonsRescale(layout.north.arrow(),
offset = c(-83, 45.5),
scale = 2)})
Difference (Barrier minus No Barrier)
rng = seq(-0.3, 0.1, 0.01)
mCols = c("#A50026", "#D73027", "#F46D43", "#FDAE61", "#E0F3F8", "#ABD9E9", "#74ADD1", "#4575B4","#313695")
cr = rev(colorRampPalette(mCols)(n = 500))
cr = colorRampPalette(cr,
bias = 0.48, space = "rgb")
levelplot(Diff.est,
margin = FALSE,
xlab = NULL,
ylab = NULL,
main = "Occurrence Probability Difference \n(Barrier minus NoBarrier)",
maxpixels = 1e5,
col.regions = cr, at = rng,
colorkey = list(labels=list(cex=1.5),
space = "bottom"),
par.strip.text = list(fontface='bold', cex=1.5),
par.settings = list(axis.line = list(col = "black"),
strip.background = list(col = 'transparent'),
strip.border = list(col = 'transparent')),
scales = list(cex = 1.25)) +
latticeExtra::layer(sp.polygons(Domain, col = "black", lwd = 0.5)) +
latticeExtra::layer(sp.polygons(Hare.pnt , col = "black", pch=factor(Hare.pnt$Hare), cex = 1)) +
latticeExtra::layer(sp.text(MyMatrix, txt = Domain$Name,
pos =c(1,3,3,2,2,1,1),
col="black",font=list(face="bold"), cex=1)) +
latticeExtra::layer({SpatialPolygonsRescale(layout.north.arrow(),
offset = c(-83, 45.5),
scale = 2)})