A Bayesian geostatistical approach to modelling global distributions of Lygodium microphyllum under projected climate warming
Authors:
John M. Humphreys (jmh09r@my.fsu.edu)
James B. Elsner
Thomas H. Jagger
Stephanie Pau
Portions of this study have been submitted for publication.
Please contact authors at email address above to cite.
Pre-Processing
Pre-processing script available here: http://rpubs.com/JMHumphreys/Lyg020217pp
Data Availabilty
Downloadable script and data to fit the full model (Model5): https://github.com/JMHumphreys/SDM_Lygodium
date()## [1] "Mon Mar 06 23:11:27 2017"Options
library(knitr)
library(rgl)
opts_knit$set(verbose = FALSE)
knit_hooks$set(webgl = hook_webgl)Load Packages:
suppressMessages(library(INLA))
suppressMessages(library(rgdal))
suppressMessages(library(reshape2))
suppressMessages(library(raster))
suppressMessages(library(fields))
suppressMessages(library(rasterVis))
suppressMessages(library(rgeos))
suppressMessages(library(Thermimage))
suppressMessages(library(maptools))
suppressMessages(library(mapproj))
suppressMessages(library(spdep))
suppressMessages(library(ggplot2))
suppressMessages(library(GISTools))
suppressMessages(library(lattice))
suppressMessages(library(gridExtra))
suppressMessages(library(dismo))
suppressMessages(library(dplyr))
suppressMessages(library(threejs))
suppressMessages(library(xtable))Set Directory
setwd("D:/Lygodium/LygProject/LygProject")Get Data
load("Results022517.RData")Locations of Observations in 2D
Cp = levelplot(Domain.r,
margin = FALSE,
#sub = "Density of Random Effects",
xlab = NULL, ylab = NULL,
col.regions = cr, at = rng,
colorkey = NULL,
panel = panel.levelplot, space = "right", par.strip.text = list(fontface='bold'),
col = cr, par.settings = list(axis.line = list(col = "black"),
strip.background = list(col = 'transparent'),
strip.border = list(col = 'transparent')),
scales = list(col = "black", cex = 2, relation="free")) +
latticeExtra::layer(sp.polygons(Countries, col = "black", lwd = 1)) +
latticeExtra::layer(sp.polygons(LygC, col = "red", pch=19, cex=0.5))
Cp
View Occurrence Locations in 3D
Interactive plot of mesh (click, zoom, drag, etc..)
earth = "D:/DRAFT_Paper/Images020417/world.topo.bathy.200412.3x5400x2700.jpg"
globejs(img=earth,
lat=LygC$Lat,
long=LygC$Long,
value=25, pointsize=2,
color="red", bg="white",
atmosphere=TRUE)Count of Occurrences by Latitude
myLL = cbind(Lyg.df$Long, Lyg.df$Lat)
LygLat = SpatialPointsDataFrame(myLL, Lyg.df)
proj4string(LygLat) = LL84
LygLat$Interval = as.factor(cut(LygLat$Lat,
seq(-90, 90, 10),
labels = FALSE))
Lat.dat = LygPA@data
Lat.dat$Interval = as.factor(cut(Lat.dat$Lat,
seq(-90, 90, 10),
labels = FALSE))
Lat.dat2 = as.data.frame(
Lat.dat %>%
group_by(Interval) %>%
summarise(Count = sum(OBS)))
FreqBar = ggplot(Lat.dat2, aes(x = factor(Interval), y = Count)) +
geom_bar(stat = "identity") +
scale_x_discrete(labels = seq(-90, 90, 10)) +
ylim(0, 250) +
xlab(expression(""*~degree*"Latitude")) +
ylab("Count of Observations") +
theme_classic() +
theme(axis.title.y = element_text(size = 18),
axis.title.x = element_text(size = 18))
FreqBar
View Mesh
Interactive plot of mesh (click, zoom, drag, etc..)
Mesh is designed to have larger triangulations over the oceans and smaller ones over terrestrial areas.
O.pnts$ID = !is.na(over(O.pnts, Countries))[,1]
cp = colorRampPalette(c("lightblue", "tan"))
plot(mesh1,
rgl=TRUE,
col=O.pnts$ID,
color.palette=cp,
draw.vertices=FALSE,
draw.edges=TRUE)You must enable Javascript to view this page properly.
Convert to Cartesian Coordinates
Converting the 2D long-lat coordinates from the observation and quadrature locations to 3D Cartesian coordinates (xyz) and then joining them the mesh.
locs = cbind(LygPA$Long, LygPA$Lat)
locs = inla.mesh.map(locs,
projection = "longlat",
inverse = TRUE)
A = inla.spde.make.A(mesh1, loc=locs) #project point locations to meshPrior for SPDE Model
An informative “PC” Prior is specified based on the mesh parameters. Also creating a spatial index (“field0”).
FE.df = LygPA@data
spde = inla.spde2.pcmatern(mesh1,
prior.range=c(0.9, 0.9), #Specifing priori probability of 90%
prior.sigma=c(1, 0.01)) #that the range will be less than 0.9.
field0 = inla.spde.make.index("field0",
spde$n.spde) #spatial indexConstruct Data Stack
Data stacks are used to put all fixed and random variables into a more easily managed object.
FE.lst = list(c(field0, #Spatial Index
list(intercept0 = 1)), #Intercept
list(mTcm = FE.df[,"mTcm"], #Min Temp of Coldest Month
Pwq = FE.df[,"Pwq"], #Precipitation of Wettest Quarter
Pop = FE.df[,"Pop"], #Human Population density
CTI = FE.df[,"CTI"], #Compound Topographic Index
NN = FE.df[,"NN"])) #Distance to nearest plant occurrence
#Stack0 for model fitting
Stack0 = inla.stack(data = list(OBS = FE.df$OBS), #Response. Presence (1) or Quadrature (0)
A = list(A, 1), #Intersection of points with mesh
effects = FE.lst, #Index, Intercept, and Covariates
tag = "obs0") #Just a label SPDE Spatial Effect Only (Model0)
Course spatial effect with no other covariates.
Frm0 = OBS ~ -1 + intercept0 +
f(field0, model=spde)
#theta0 = Model0$internal.summary.hyperpar$mean
theta0 = c(-0.5992131, 1.6128842) #previous run
Model0 = inla(Frm0,
data = inla.stack.data(Stack0, spde=spde),
family = "binomial",
verbose = TRUE,
control.fixed = list(prec = 0.1, prec.intercept = 0.1),
control.predictor = list(
A = inla.stack.A(Stack0),
compute = TRUE,
link = 1),
control.mode = list(restart = TRUE, theta = theta0),
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(Stack0, ", " spde = spde), verbose = TRUE, control.compute = list(dic = TRUE, ", " cpo = TRUE, waic = TRUE), control.predictor = list(A = inla.stack.A(Stack0), ", " compute = TRUE, link = 1), control.results = list(return.marginals.random = TRUE, ", " return.marginals.predictor = TRUE), control.fixed = list(prec = 0.1, ", " prec.intercept = 0.1), control.mode = list(restart = TRUE, ", " theta = theta0))")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 1.6056 3444.7690 2.8695 3449.2441
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept0 -9.1654 0.9675 -11.201 -9.1189 -7.3893 -9.0235 0
##
## Random effects:
## Name Model
## field0 SPDE2 model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant 0.975quant mode
## Range for field0 0.5297 0.0738 0.403 0.5232 0.6924 0.5094
## Stdev for field0 3.9055 0.4029 3.191 3.8777 4.7723 3.8166
##
## Expected number of effective parameters(std dev): 150.92(7.872)
## Number of equivalent replicates : 172.66
##
## Deviance Information Criterion (DIC) ...: 2288.54
## Effective number of parameters .........: 124.18
##
## Watanabe-Akaike information criterion (WAIC) ...: 2272.87
## Effective number of parameters .................: 99.31
##
## Marginal log-Likelihood: -1271.19
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computedModify SPDE Model
Before adding covariates to the model, the “spde” specification is modified. The “extraconstr” control directs the SPDE spatial effect to a “design matrix” that holds the covariates and intercept; this helps reduce confounding by ensuring that the spatial effect is orthogonal to the covariates.
n.data = dim(LygPA@data)[1]
n.covariates = 6
X = cbind(rep(1,n.data),
FE.df$mTcm,
FE.df$Pwq,
FE.df$Pop,
FE.df$CTI,
FE.df$NN)
Q = qr.Q(qr(X))
spde = inla.spde2.pcmatern(mesh1,
prior.range=c(0.9, 0.9), #Specifing priori probability of 90%
prior.sigma=c(1, 0.01), #that the range will be less than 0.9.
extraconstr = list(A = as.matrix(t(Q)%*%A),
e = rep(0, n.covariates))) Adding Covariates (Model1)
Including all fixed and random effects except Invasive Cohorts (CoH).
A “shrinking” PC Prior is specified for the mTcm random walk model.
pcprior1 = list(prec = list(prior="pc.prec", param = c(3, 0.01)))
pcprior2 = list(prec = list(prior="pc.prec", param = c(0.1, 0.01)))
Frm1 = OBS ~ -1 + intercept0 +
f(field0, model=spde) +
f(mTcm, model = "rw1", #Min Temp of Coldest Month
hyper = pcprior1,
scale.model = TRUE) +
f(NN, model = "rw1", #Distance to nearest presence
hyper = pcprior2,
scale.model = TRUE) +
Pop + CTI + Pwq #Population Density, "Wetness index", Precipitation of Wettest Quarter
#theta1 = Model1$internal.summary.hyperpar$mean
theta1 = c(-0.9202037, 0.2753616, -0.3887835, -0.7783262)
Model1 = inla(Frm1,
data = inla.stack.data(Stack0, spde=spde),
family = "binomial",
verbose = TRUE,
control.fixed = list(prec = 0.1, prec.intercept = 0.1),
control.predictor = list(
A = inla.stack.A(Stack0),
compute = TRUE,
link = 1),
control.mode = list(restart = TRUE, theta = theta1),
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(Stack0, ", " spde = spde), verbose = TRUE, control.compute = list(dic = TRUE, ", " cpo = TRUE, waic = TRUE), control.predictor = list(A = inla.stack.A(Stack0), ", " compute = TRUE, link = 1), control.results = list(return.marginals.random = TRUE, ", " return.marginals.predictor = TRUE), control.fixed = list(prec = 0.1, ", " prec.intercept = 0.1))")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 1.3144 19144.0973 2.9158 19148.3274
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept0 -9.1316 1.1489 -11.6128 -9.0485 -7.0990 -8.8748 0
## Pop 0.4393 0.1705 0.1176 0.4345 0.7880 0.4247 0
## CTI -0.9185 0.3328 -1.5759 -0.9171 -0.2692 -0.9144 0
## Pwq 1.3236 0.2312 0.8671 1.3243 1.7758 1.3257 0
##
## Random effects:
## Name Model
## field0 SPDE2 model
## mTcm RW1 model
## NN RW1 model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant 0.975quant mode
## Range for field0 0.4113 0.1065 0.2467 0.3961 0.6610 0.3666
## Stdev for field0 1.3317 0.1995 0.9866 1.3155 1.7673 1.2830
## Precision for mTcm 0.7592 0.3826 0.2658 0.6781 1.7250 0.5409
## Precision for NN 0.4686 0.0960 0.3121 0.4575 0.6874 0.4351
##
## Expected number of effective parameters(std dev): 82.22(9.695)
## Number of equivalent replicates : 316.94
##
## Deviance Information Criterion (DIC) ...: 1952.85
## Effective number of parameters .........: 78.56
##
## Watanabe-Akaike information criterion (WAIC) ...: 1946.80
## Effective number of parameters .................: 68.74
##
## Marginal log-Likelihood: -9421.89
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computedFull Model (Model2)
mTcm as a linear effect.
Frm2 = OBS ~ -1 + intercept0 +
f(field0, model=spde) +
f(NN, model = "rw1",
hyper = pcprior2,
scale.model = TRUE) +
Pop + CTI + mTcm + Pwq
#theta2 = Model2$internal.summary.hyperpar$mean
theta2 = c(-1.0737464, 0.2696745, -0.9042602)
Model2 = inla(Frm2,
data = inla.stack.data(Stack0, spde=spde),
family = "binomial",
verbose = TRUE,
control.fixed = list(prec = 0.1, prec.intercept = 0.1),
control.predictor = list(
A = inla.stack.A(Stack0),
compute = TRUE,
link = 1),
control.mode = list(restart = TRUE, theta = theta2),
control.results = list(return.marginals.random = TRUE,
return.marginals.predictor = TRUE),
control.compute=list(dic = TRUE, cpo = TRUE, waic = TRUE)) summary(Model2) ##
## Call:
## c("inla(formula = Frm2, family = \"binomial\", data = inla.stack.data(Stack0, ", " spde = spde), verbose = TRUE, control.compute = list(dic = TRUE, ", " cpo = TRUE, waic = TRUE), control.predictor = list(A = inla.stack.A(Stack0), ", " compute = TRUE, link = 1), control.results = list(return.marginals.random = TRUE, ", " return.marginals.predictor = TRUE), control.fixed = list(prec = 0.1, ", " prec.intercept = 0.1))")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 1.7852 5705.8224 2.1268 5709.7344
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept0 -8.4046 1.0373 -10.6554 -8.3248 -6.5829 -8.1562 0
## Pop 0.4147 0.1675 0.0993 0.4097 0.7580 0.3995 0
## CTI -1.0545 0.3313 -1.7090 -1.0532 -0.4082 -1.0504 0
## mTcm 1.1709 0.2033 0.7857 1.1658 1.5848 1.1556 0
## Pwq 1.2423 0.2351 0.7784 1.2429 1.7021 1.2440 0
##
## Random effects:
## Name Model
## field0 SPDE2 model
## NN RW1 model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant 0.975quant mode
## Range for field0 0.3680 0.1357 0.1405 0.3585 0.6553 0.3252
## Stdev for field0 1.3317 0.2393 0.8758 1.3322 1.8027 1.3583
## Precision for NN 0.4126 0.0806 0.2720 0.4071 0.5874 0.3976
##
## Expected number of effective parameters(std dev): 70.19(10.67)
## Number of equivalent replicates : 371.26
##
## Deviance Information Criterion (DIC) ...: 1968.62
## Effective number of parameters .........: 67.86
##
## Watanabe-Akaike information criterion (WAIC) ...: 1963.33
## Effective number of parameters .................: 59.45
##
## Marginal log-Likelihood: -7601.84
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computedRemoving Non-linear Effects for NN and mTcm (Model3)
both NN and mTcm as linear effects.
Frm3 = OBS ~ -1 + intercept0 +
f(field0, model=spde) +
Pop + CTI + mTcm + Pwq + NN
#theta3 = Model3$internal.summary.hyperpar$mean
theta3 = c(-0.9339132, 0.3567992)
Model3 = inla(Frm3,
data = inla.stack.data(Stack0, spde=spde),
family = "binomial",
verbose = TRUE,
control.fixed = list(prec = 0.1, prec.intercept = 0.1),
control.predictor = list(
A = inla.stack.A(Stack0),
compute = TRUE,
link = 1),
control.mode = list(restart = TRUE, theta = theta3),
control.results = list(return.marginals.random = TRUE,
return.marginals.predictor = TRUE),
control.compute=list(dic = TRUE, cpo = TRUE, waic = TRUE)) summary(Model3) ##
## Call:
## c("inla(formula = Frm3, family = \"binomial\", data = inla.stack.data(Stack0, ", " spde = spde), verbose = TRUE, control.compute = list(dic = TRUE, ", " cpo = TRUE, waic = TRUE), control.predictor = list(A = inla.stack.A(Stack0), ", " compute = TRUE, link = 1), control.results = list(return.marginals.random = TRUE, ", " return.marginals.predictor = TRUE), control.fixed = list(prec = 0.1, ", " prec.intercept = 0.1))")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 1.6665 9143.8050 3.8676 9149.3391
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept0 -2.0686 0.5841 -3.2248 -2.0655 -0.9306 -2.0593 0
## Pop 0.4438 0.1652 0.1367 0.4373 0.7869 0.4239 0
## CTI -1.0337 0.3274 -1.6804 -1.0324 -0.3948 -1.0298 0
## mTcm 1.2662 0.2156 0.8572 1.2611 1.7045 1.2508 0
## Pwq 1.3346 0.2435 0.8534 1.3356 1.8099 1.3377 0
## NN -7.3302 0.7259 -8.7786 -7.3224 -5.9266 -7.3066 0
##
## Random effects:
## Name Model
## field0 SPDE2 model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant 0.975quant mode
## Range for field0 0.4079 0.1139 0.2313 0.392 0.6757 0.362
## Stdev for field0 1.4417 0.1942 1.0926 1.431 1.8542 1.413
##
## Expected number of effective parameters(std dev): 62.93(7.956)
## Number of equivalent replicates : 414.05
##
## Deviance Information Criterion (DIC) ...: 2061.03
## Effective number of parameters .........: 61.54
##
## Watanabe-Akaike information criterion (WAIC) ...: 2064.13
## Effective number of parameters .................: 61.47
##
## Marginal log-Likelihood: -1080.24
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computedDropping Clustering Term (Model4)
Excluding “NN”
Frm4 = OBS ~ -1 + intercept0 +
f(field0, model=spde) +
f(mTcm, model = "rw1",
hyper = pcprior1,
scale.model = TRUE) +
Pop + CTI + Pwq
#theta4 = Model4$internal.summary.hyperpar$mean
theta4 = c(-1.033992, 1.076775, -1.173984)
Model4 = inla(Frm4,
data = inla.stack.data(Stack0, spde=spde),
family = "binomial",
verbose = TRUE,
control.fixed = list(prec = 0.1, prec.intercept = 0.1),
control.predictor = list(
A = inla.stack.A(Stack0),
compute = TRUE,
link = 1),
control.mode = list(restart = TRUE, theta = theta4),
control.results = list(return.marginals.random = TRUE,
return.marginals.predictor = TRUE),
control.compute=list(dic = TRUE, cpo = TRUE, waic = TRUE)) summary(Model4) ##
## Call:
## c("inla(formula = Frm4, family = \"binomial\", data = inla.stack.data(Stack0, ", " spde = spde), verbose = TRUE, control.compute = list(dic = TRUE, ", " cpo = TRUE, waic = TRUE), control.predictor = list(A = inla.stack.A(Stack0), ", " compute = TRUE, link = 1), control.results = list(return.marginals.random = TRUE, ", " return.marginals.predictor = TRUE), control.fixed = list(prec = 0.1, ", " prec.intercept = 0.1))")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 2.2412 5279.4870 3.9241 5285.6522
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept0 -9.4234 0.9978 -11.5655 -9.3583 -7.6429 -9.2271 0
## Pop 0.6724 0.1777 0.3282 0.6708 1.0255 0.6677 0
## CTI -1.3001 0.3440 -1.9796 -1.2987 -0.6289 -1.2959 0
## Pwq 1.8449 0.2617 1.3298 1.8453 2.3576 1.8461 0
##
## Random effects:
## Name Model
## field0 SPDE2 model
## mTcm RW1 model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant 0.975quant mode
## Range for field0 0.3599 0.0561 0.2635 0.3551 0.4833 0.3454
## Stdev for field0 2.9498 0.2951 2.4196 2.9327 3.5766 2.8968
## Precision for mTcm 0.3321 0.1313 0.1491 0.3077 0.6553 0.2650
##
## Expected number of effective parameters(std dev): 152.96(8.474)
## Number of equivalent replicates : 170.36
##
## Deviance Information Criterion (DIC) ...: 2064.75
## Effective number of parameters .........: 133.76
##
## Watanabe-Akaike information criterion (WAIC) ...: 2052.97
## Effective number of parameters .................: 112.19
##
## Marginal log-Likelihood: -2975.33
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computedConstructing Biotic Variable (CoH)
Here, species richness for 31 invasive species associated with L. microphyllum (in Florida) is estimated. First raw counts are taken of observations and then the samples are corrected for spatial correlation and sampling bias.
Data Pre-processing
r = raster(res = 1,
crs = proj4string(Countries))
extent(r) = extent(Countries)
World.b = gBuffer(Countries,
byid = FALSE,
width = 0.5)
World.r = rasterize(World.b, r,
field = 0,
background = NA)
CohSpp.df = read.csv("./InvCohortSpp.csv",
header = TRUE, sep=",")
SppLst = as.data.frame(
CohSpp.df %>%
group_by(species) %>%
summarise(no_rows = length(species)))
names(SppLst) = c("Species", "Count")
SppLst## Species Count
## 1 Abrus precatorius 4549
## 2 Alternanthera philoxeroides 1464
## 3 Ardisia elliptica 588
## 4 Bischofia javanica 1514
## 5 Callistemon viminalis 1333
## 6 Epipremnum pinnatum 618
## 7 Eragrostis atrovirens 879
## 8 Ficus microcarpa 787
## 9 Hymenachne amplexicaulis 2920
## 10 Ixora coccinea 448
## 11 Limnophila sessiliflora 251
## 12 Mangifera indica 2729
## 13 Melaleuca quinquenervia 19490
## 14 Neyraudia reynaudiana 1876
## 15 Oryza rufipogon 5838
## 16 Panicum maximum 7926
## 17 Panicum repens 8168
## 18 Pouzolzia zeylanica 518
## 19 Psidium guajava 5465
## 20 Pteris vittata 2606
## 21 Sacciolepis indica 2193
## 22 Schefflera actinophylla 1790
## 23 Schinus terebinthifolia 2222
## 24 Senna pendula 2289
## 25 Solanum diphyllum 843
## 26 Spathodea campanulata 741
## 27 Spirodela punctata 780
## 28 Syngonium podophyllum 2005
## 29 Syzygium cumini 1287
## 30 Terminalia catappa 1546
## 31 Youngia japonica 2344#SpTable = as.data.frame(SppLst[,1])
#print(xtable(SppLst), include.rownames = TRUE)
###Loop to rasterize and then stack species individually
SppNames = levels(as.factor(as.character(CohSpp.df$species)))
Spp.stk = stack()
for (i in 1:length(SppNames)) {
Spp1 = CohSpp.df %>% filter(species == SppNames[i])
coords1 = cbind(Spp1$decimallongitude,
Spp1$decimallatitude)
Ras1 = rasterize(coords1, World.r,
field = 1,
background = 0)
Ras1 = Ras1 + World.r
names(Ras1) = paste0(SppNames[i])
if (length(SppNames) == 1){
Spp.stk = Ras1}
else {Spp.stk = stack(Spp.stk, Ras1)}
}
Richness = sum(Spp.stk)
rng = seq(0, 27, 1)
cr = colorRampPalette(c("tan", "blue", "lightblue","green", "lightgreen", "yellow",
"orange", "red", "darkred"),
space = "rgb")
Cp = levelplot(Richness,
margin = FALSE,
sub = "Invasive Cohort Richness\n(raw)",
xlab = NULL, ylab = NULL,
col.regions = cr, at = rng,
colorkey = list(at = rng, col = cr),
par.settings = list(fontsize = list(text = 13)))
Cp
#Convert to points and extract Population Denisty (Pop)
Rich.pnts = rasterToPoints(Richness, sp=TRUE)
names(Rich.pnts) = "SppR"
Rich.pnts$Pop = extract(Pop, Rich.pnts, method = "simple")/1000
Rich.pnts$Pop = ifelse(is.na(Rich.pnts$Pop) == TRUE,
mean(Rich.pnts$Pop, na.rm = TRUE),
Rich.pnts$Pop)
Rich.pnts@data = Rich.pnts@data %>%
mutate(Long = Rich.pnts@coords[,1],
Lat = Rich.pnts@coords[,2])
#Convert to Cartesian and project to mesh
Rlocs = cbind(Rich.pnts$Long, Rich.pnts$Lat)
Rlocs = inla.mesh.map(Rlocs,
projection = "longlat",
inverse = TRUE)
AR = inla.spde.make.A(mesh1, loc=Rlocs)
#Set prior and create spatial index
FEr.df = Rich.pnts@data
spdeR = inla.spde2.pcmatern(mesh1,
prior.range=c(0.9, 0.9),
prior.sigma=c(1, 0.01))
fieldR = inla.spde.make.index("fieldR",
spdeR$n.spde)
FER.lst = list(c(fieldR,
list(interceptR = 1)),
list(Pop = FEr.df[,"Pop"]))
#Stack for model fitting
StackR = inla.stack(data = list(Rich = FEr.df$SppR),
A = list(AR, 1),
effects = FER.lst,
tag = "R0")
#Stack1 for predicting global distributions
StackR2 = inla.stack(data = list(Rich = NA),
A = list(AR, 1),
effects = FER.lst,
tag = "pR0")
#Combined stack for fitting and prediction
RStack.global = inla.stack(StackR, StackR2)Model for Richness (ModelR)
Including a spde spatial term and Pop covariate.
FrmR = Rich ~ -1 + interceptR +
f(fieldR, model=spdeR) + Pop
#thetaR = ModelR$internal.summary.hyperpar$mean
thetaR = c(-0.5478317, 1.2696178)
ModelR = inla(FrmR,
data = inla.stack.data(RStack.global, spde=spdeR),
family = "poisson",
verbose = TRUE,
control.fixed = list(prec = 0.1, prec.intercept = 0.1),
control.predictor = list(
A = inla.stack.A(RStack.global),
compute = TRUE,
link = 1),
control.mode = list(restart = TRUE, theta = thetaR),
control.results = list(return.marginals.random = TRUE,
return.marginals.predictor = TRUE),
control.compute=list(dic = TRUE, cpo = TRUE, waic = TRUE)) summary(ModelR)##
## Call:
## c("inla(formula = FrmR, family = \"poisson\", data = inla.stack.data(RStack.global, ", " spde = spdeR), verbose = TRUE, control.compute = list(dic = TRUE, ", " cpo = TRUE, waic = TRUE), control.predictor = list(A = inla.stack.A(RStack.global), ", " compute = TRUE, link = 1), control.results = list(return.marginals.random = TRUE, ", " return.marginals.predictor = TRUE), control.fixed = list(prec = 0.1, ", " prec.intercept = 0.1))")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 1.429 4118.430 2.500 4122.359
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## interceptR -3.3413 0.7653 -4.8757 -3.3346 -1.8438 -3.3204 0
## Pop 0.3262 0.0315 0.2633 0.3265 0.3871 0.3272 0
##
## Random effects:
## Name Model
## fieldR SPDE2 model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant 0.975quant mode
## Range for fieldR 0.5811 0.0589 0.4785 0.5764 0.7094 0.5654
## Stdev for fieldR 3.5720 0.3012 3.0387 3.5509 4.2176 3.5016
##
## Expected number of effective parameters(std dev): 463.87(7.631)
## Number of equivalent replicates : 50.92
##
## Deviance Information Criterion (DIC) ...: 18318.03
## Effective number of parameters .........: 450.14
##
## Watanabe-Akaike information criterion (WAIC) ...: 18662.24
## Effective number of parameters .................: 706.80
##
## Marginal log-Likelihood: -9654.06
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computedCompare Raw and Corrected Invasive Cohort Richness
I0 = inla.stack.index(RStack.global, "pR0") #Index to pull results
E0 = ModelR$summary.linear.predictor[I0$data,] #pull results
Pred.pnts = Rich.pnts
Pred.pnts$Cp = exp(E0[,"mean"])
CoHp = rasterize(Pred.pnts, World.r, "Cp", #Rasterize predicted values
background = NA)
CoHp = sum(CoHp, World.r)
CoH.stk = stack(Richness, CoHp)
names(CoH.stk) = c("Raw", "Corrected")
rng = seq(0, 26, 0.001)
cr = colorRampPalette(c("tan", "blue", "lightblue", "yellow", "red", "darkred"),
bias = 1.5, space = "rgb")
Cp = levelplot(CoH.stk,
margin = FALSE,
sub = "Invasive Cohort Richness\n(31 Species)",
xlab = NULL, ylab = NULL,
col.regions = cr, at = rng,
colorkey = list(at=seq(0, 33, 0.001),
labels=list(at=seq(0, 30, 5),
labels=c("0", "5", "10", "15", "20", "25", "30")),
panel = panel.levelplot, space = "right",
col = cr, par.settings = list(fontsize = list(text = 18))))
Cp +
latticeExtra::layer(sp.polygons(Countries, col = "dimgray", lwd = 1))
#Add field to dataset
LygPA$CoH = drop(A %*% ModelR$summary.random$fieldR$mean) Update Data Stacks to include CoH
FE.lst = list(c(field0,
list(intercept0 = 1)),
list(mTcm = FE.df[,"mTcm"],
Pwq = FE.df[,"Pwq"],
Pop = FE.df[,"Pop"],
CTI = FE.df[,"CTI"],
NN = FE.df[,"NN"],
CoH = FE.df[,"CoH"])) #Adding Invasive Cohort Richness
#Stack0 for model fitting
Stack0 = inla.stack(data = list(OBS = FE.df$OBS),
A = list(A, 1),
effects = FE.lst,
tag = "obs0")
#Stack1 for predicting global distributions
Stack1 = inla.stack(data = list(OBS = NA), #Setting the response to "NA" for prediction
A = list(A, 1),
effects = FE.lst,
tag = "pred0")
#Combined stack for global prediction
CStack.global = inla.stack(Stack0, Stack1)Modify SPDE Model
Before adding the CoH covariate to the model, the “spde” specification is updated.
n.data = dim(LygPA@data)[1]
n.covariates = 7
X = cbind(rep(1,n.data),
FE.df$mTcm,
FE.df$Pwq,
FE.df$Pop,
FE.df$CTI,
FE.df$NN,
FE.df$CoH) #Adding CoH
Q = qr.Q(qr(X))
spde = inla.spde2.pcmatern(mesh1,
prior.range=c(0.9, 0.9),
prior.sigma=c(1, 0.01),
extraconstr = list(A = as.matrix(t(Q)%*%A),
e = rep(0, n.covariates))) Full Model (Model5)
Adding the CoH covariate to best performing model (Model1).
#theta5 = Model5$internal.summary.hyperpar$mean
theta5 = c(-0.89178665, -0.06363225, -0.32293059, -0.64617359)
Frm5 = OBS ~ -1 + intercept0 +
f(field0, model=spde) +
f(mTcm, model = "rw1",
hyper = pcprior1,
scale.model = TRUE) +
f(NN, model = "rw1",
hyper = pcprior2,
scale.model = TRUE) +
Pop + CTI + Pwq + CoH
Model5 = inla(Frm5,
data = inla.stack.data(Stack0, spde=spde),
family = "binomial",
verbose = TRUE,
control.fixed = list(prec = 0.1, prec.intercept = 0.1),
control.predictor = list(
A = inla.stack.A(Stack0),
compute = TRUE,
link = 1),
control.mode = list(restart = TRUE, theta = theta5),
control.results = list(return.marginals.random = TRUE,
return.marginals.predictor = TRUE),
control.compute=list(dic = TRUE, cpo = TRUE, waic = TRUE)) summary(Model5) ##
## Call:
## c("inla(formula = Frm5, family = \"binomial\", data = inla.stack.data(Stack0, ", " spde = spde), verbose = TRUE, control.compute = list(dic = TRUE, ", " cpo = TRUE, waic = TRUE), control.predictor = list(A = inla.stack.A(Stack0), ", " compute = TRUE, link = 1), control.results = list(return.marginals.random = TRUE, ", " return.marginals.predictor = TRUE), control.fixed = list(prec = 0.1, ", " prec.intercept = 0.1), control.mode = list(restart = TRUE, ", " theta = theta1))")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 1.3204 18901.7118 2.7386 18905.7708
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept0 -9.2261 1.5352 -12.2967 -9.2083 -6.2573 -9.1739 0
## Pop 0.3328 0.1556 0.0441 0.3266 0.6555 0.3138 0
## CTI -0.8211 0.3281 -1.4692 -0.8198 -0.1807 -0.8172 0
## Pwq 0.9616 0.2332 0.5004 0.9626 1.4169 0.9645 0
## CoH 0.7397 0.0952 0.5569 0.7382 0.9309 0.7354 0
##
## Random effects:
## Name Model
## field0 SPDE2 model
## mTcm RW1 model
## NN RW1 model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant 0.975quant mode
## Range for field0 0.4343 0.1543 0.2153 0.4064 0.8123 0.3565
## Stdev for field0 0.9570 0.1922 0.6394 0.9366 1.3910 0.8962
## Precision for mTcm 0.8250 0.4544 0.2694 0.7203 1.9876 0.5540
## Precision for NN 0.5366 0.1189 0.3452 0.5221 0.8098 0.4933
##
## Expected number of effective parameters(std dev): 63.03(8.888)
## Number of equivalent replicates : 413.40
##
## Deviance Information Criterion (DIC) ...: 1931.42
## Effective number of parameters .........: 62.68
##
## Watanabe-Akaike information criterion (WAIC) ...: 1927.74
## Effective number of parameters .................: 56.41
##
## Marginal log-Likelihood: -9393.84
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computedNon-Spatial (Model6)
To assess the significance of the spatial effect.
Frm6 = OBS ~ -1 + intercept0 + NN + mTcm + Pop + CTI + Pwq + CoH
Model6 = inla(Frm6,
data = inla.stack.data(Stack0),
family = "binomial",
verbose = TRUE,
control.fixed = list(prec = 0.1, prec.intercept = 0.1),
control.predictor = list(
A = inla.stack.A(Stack0),
compute = TRUE,
link = 1),
control.results = list(return.marginals.random = TRUE,
return.marginals.predictor = TRUE),
control.compute=list(dic = TRUE, cpo = TRUE, waic = TRUE)) summary(Model6) ##
## Call:
## c("inla(formula = Frm6, family = \"binomial\", data = inla.stack.data(Stack0), ", " verbose = TRUE, control.compute = list(dic = TRUE, cpo = TRUE, ", " waic = TRUE), control.predictor = list(A = inla.stack.A(Stack0), ", " compute = TRUE, link = 1), control.results = list(return.marginals.random = TRUE, ", " return.marginals.predictor = TRUE), control.fixed = list(prec = 0.1, ", " prec.intercept = 0.1))")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 0.3671 87.5464 4.1186 92.0321
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept0 -4.5114 0.5285 -5.5553 -4.5094 -3.4799 -4.5052 0
## mTcm 0.8413 0.1208 0.6066 0.8404 1.0808 0.8386 0
## NN -9.6922 0.6949 -11.0868 -9.6819 -8.3556 -9.6610 0
## Pop 0.2025 0.1193 -0.0113 0.1950 0.4573 0.1794 0
## CTI -0.1561 0.2792 -0.7053 -0.1557 0.3907 -0.1550 0
## Pwq 0.3157 0.1943 -0.0689 0.3168 0.6943 0.3188 0
## CoH 0.8574 0.0597 0.7418 0.8568 0.9759 0.8557 0
##
## The model has no random effects
##
## The model has no hyperparameters
##
## Expected number of effective parameters(std dev): 6.912(0.00)
## Number of equivalent replicates : 3769.75
##
## Deviance Information Criterion (DIC) ...: 2201.34
## Effective number of parameters .........: 6.886
##
## Watanabe-Akaike information criterion (WAIC) ...: 2205.32
## Effective number of parameters .................: 10.33
##
## Marginal log-Likelihood: -1121.40
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computedModel Results
Compare Model Metrics
Models = c("Model0", "Model1",
"Model2", "Model3",
"Model4", "Model5",
"Model6")
Fixed = c("W (Course spatial structure only) ",
"f(NN) + f(mTcm) + Pop + CTI + Pwq + W",
"f(NN) + mTcm + Pop + CTI + Pwq + W",
"NN + mTcm + Pop + CTI + Pwq + W",
"f(mTcm) + Pop + CTI + Pwq + W",
"f(NN) + f(mTcm) + Pop + CTI + Pwq + CoH + W",
"NN + mTcm + Pop + CTI + Pwq + CoH")
#Deviance information criteria
DICs = c(Model0$dic$dic, Model1$dic$dic,
Model2$dic$dic, Model3$dic$dic,
Model4$dic$dic, Model5$dic$dic,
Model6$dic$dic)
#Watanabe-Akaike information criteria
WAICs = c(Model0$waic$waic, Model1$waic$waic,
Model2$waic$waic, Model3$waic$waic,
Model4$waic$waic, Model5$waic$waic,
Model6$waic$waic)
#log Conditional Predictive Ordnance
LCPOs = c(-mean(log(Model0$cpo$cpo)), -mean(log(Model1$cpo$cpo)),
-mean(log(Model2$cpo$cpo)), -mean(log(Model3$cpo$cpo)),
-mean(log(Model4$cpo$cpo)), -mean(log(Model5$cpo$cpo)),
-mean(log(Model6$cpo$cpo)))
Model_mets = as.data.frame(cbind(Model = Models,
DIC = round(DICs, 2),
WAIC = round(WAICs, 2),
LCPO = round(LCPOs, 3),
COVARIATES = Fixed))Model Comparison
print.data.frame(Model_mets, right = FALSE)## Model DIC WAIC LCPO COVARIATES
## 1 Model0 2288.54 2272.87 0.044 W (Course spatial structure only)
## 2 Model1 1952.85 1946.8 0.037 f(NN) + f(mTcm) + Pop + CTI + Pwq + W
## 3 Model2 1968.62 1963.33 0.038 f(NN) + mTcm + Pop + CTI + Pwq + W
## 4 Model3 2061.03 2064.13 0.04 NN + mTcm + Pop + CTI + Pwq + W
## 5 Model4 2064.75 2052.97 0.039 f(mTcm) + Pop + CTI + Pwq + W
## 6 Model5 1931.42 1927.74 0.037 f(NN) + f(mTcm) + Pop + CTI + Pwq + CoH + W
## 7 Model6 2201.34 2205.32 0.042 NN + mTcm + Pop + CTI + Pwq + CoH#ModComp = xtable(Model_mets)
#print(ModComp, include.rownames = FALSE,
# floating.environment = "sidewaystable")Selected Model (Fixed Effects)
SelMod.tab = as.data.frame(Model5$summary.fixed)
SelMod.tab2 = cbind(SelMod.tab[,1:3], SelMod.tab[,5])
names(SelMod.tab2) = c("Mean", "sd", "2.5% Q", "97.5% Q")
print.data.frame(SelMod.tab2, right = FALSE)## Mean sd 2.5% Q 97.5% Q
## intercept0 -9.2260725 1.53523445 -12.29672012 -6.2572756
## Pop 0.3327789 0.15558055 0.04405673 0.6555362
## CTI -0.8210867 0.32813052 -1.46920598 -0.1807039
## Pwq 0.9616265 0.23321119 0.50039125 1.4169286
## CoH 0.7396538 0.09516221 0.55686433 0.9308740#SelMod.tab2 = xtable(SelMod.tab2)
#print(SelMod.tab2, include.rownames = TRUE)Compare Random Effects
proj = inla.mesh.projector(mesh1,
dims=c(400, 800))
plotdata = inla.mesh.project(proj, Model0$summary.random$field0[,"mean"])
plotdata2 = inla.mesh.project(proj, Model5$summary.random$field0[,"mean"])
R600 = raster(nrows = 400, ncols = 800)
R601 = R600
values(R600) = plotdata
values(R601) = plotdata2
M600 = as.matrix(R600)
M601 = as.matrix(R601)
M0.spm = rotate90.matrix(M600) #rotation before rasterization
M0.spmr = raster(M0.spm)
extent(M0.spmr) = c(-180, 180, -90, 90)
proj4string(M0.spmr) = CRS(proj4string(Domain.r))
M1.spm = rotate90.matrix(M601) #rotation before rasterization
M1.spmr = raster(M1.spm)
extent(M1.spmr) = c(-180, 180, -90, 90)
proj4string(M1.spmr) = CRS(proj4string(Domain.r))
RandStack = stack(M0.spmr, M1.spmr)
names(RandStack) = c("A", "B")
rng = seq(-4, 13, 0.01)
cr = colorRampPalette(c("darkblue", "blue", "lightblue",
"yellow", "orange", "red", "darkred"),
bias = 1.5, space = "rgb")
Cp = levelplot(RandStack,
margin = FALSE,
#sub = "Density of Random Effects",
xlab = NULL, ylab = NULL,
col.regions = cr, at = rng,
colorkey = list(at=seq(-5, 13, 0.01),
labels=list(at=c(-4, 0, 4, 8, 12)),
labels=c("-4", "0", "4", "8", "12")),
panel = panel.levelplot, space = "right", par.strip.text = list(fontface='bold'),
col = cr, par.settings = list(axis.line = list(col = "black"),
strip.background = list(col = 'transparent'),
strip.border = list(col = 'transparent')),
scales = list(col = "black", cex = 2))
Cp +
latticeExtra::layer(sp.polygons(Countries, col = "black", lwd = 1))
Compare Standard Deviation
plotdata = inla.mesh.project(proj, Model0$summary.random$field0[,"sd"])
plotdata2 = inla.mesh.project(proj, Model5$summary.random$field0[,"sd"])
R600 = raster(nrows = 400, ncols = 800)
R601 = R600
values(R600) = plotdata
values(R601) = plotdata2
M600 = as.matrix(R600)
M601 = as.matrix(R601)
M0.spm = rotate90.matrix(M600) #rotation before rasterization
M0.spmr = raster(M0.spm)
extent(M0.spmr) = c(-180, 180, -90, 90)
proj4string(M0.spmr) = CRS(proj4string(Domain.r))
M1.spm = rotate90.matrix(M601) #rotation before rasterization
M1.spmr = raster(M1.spm)
extent(M1.spmr) = c(-180, 180, -90, 90)
proj4string(M1.spmr) = CRS(proj4string(Domain.r))
RandStack = stack(M0.spmr, M1.spmr)
names(RandStack) = c("A", "B")
rng = seq(0, 4.4, 0.01)
cr = colorRampPalette(c("lightblue", "yellow", "darkred"),
space = "rgb")
Cp = levelplot(RandStack,
margin = FALSE,
#sub = "Standard Deviation",
xlab = NULL, ylab = NULL,
col.regions = cr, at = rng,
colorkey = list(at=seq(0, 4.4, 0.01),
labels=list(at=c(0, 1, 2, 3, 4)),
labels=c("0", "1", "2", "3", "4")),
panel = panel.levelplot, space = "right", par.strip.text = list(fontface='bold'),
col = cr, par.settings = list(axis.line = list(col = "black"),
strip.background = list(col = 'transparent'),
strip.border = list(col = 'transparent')),
scales = list(col = "black", cex = 2))
Cp +
latticeExtra::layer(sp.polygons(Countries, col = "black", lwd = 1))
Temperature Limitation
Reduced survival in cold temperatures.
mic.df = as.data.frame(Model5$summary.random$mTcm[,1:6])
names(mic.df) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975")
ggplot(mic.df, aes(ID*10, Mean)) +
geom_smooth(col = "black",
linetype= "solid",
method = "loess",
span = 0.25,
se = FALSE,
lwd = 1) +
geom_smooth(data = mic.df, aes(ID*10, Q025),
col = "grey",
method = "loess",
span = 0.25,
se = FALSE,
linetype= "dashed") +
geom_smooth(data = mic.df, aes(ID*10, Q975),
col = "grey",
method = "loess",
span = 0.25,
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(expression("Limiting Temperature ("*~degree*C*" )")) +
ylab("Logit") +
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))
Find Limiting Temperture
Point at which logit becomes positive.
FindZero = function(X.df) {
Cross = which(diff(sign(X.df$Mean))!=0)
mTemp = (X.df$ID[Cross] + X.df$ID[Cross+1])/2
return(mTemp)
}
###Temperature Limit (degrees C) #3.85
FindZero(mic.df)[2]*10## [1] 3.85Fine-scale Spatial Structure
Accounting for the observed clustering pattern (possible related to dispersal ability).
mic.df = as.data.frame(Model5$summary.random$NN[,1:6])
names(mic.df) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975")
Full.p = ggplot(mic.df, aes(ID*1000, Mean)) +
geom_smooth(method = "loess",
se = FALSE, col = "black",
linetype= "solid") +
geom_smooth(data = mic.df, aes(ID*1000, Q025),
method = "loess",
se = FALSE, col = "grey",
linetype= "dashed") +
geom_smooth(data = mic.df, aes(ID*1000, Q975),
method = "loess",
se = FALSE, col = "grey",
linetype= "dashed") +
geom_hline(yintercept = 0,
linetype = "dotted",
col = "red",
size = 1) +
geom_vline(xintercept = 0,
linetype = "dotted",
col = "red",
size = 1) +
xlim(c(0,28000)) +
xlab(NULL) +
ylab("Logit") +
ggtitle("A") +
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(hjust = 0.5))
###Zoom in some
mic.df = as.data.frame(Model5$summary.random$NN[,1:6])
names(mic.df) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975")
mic.df = mic.df %>%
filter(ID <= 2)
Zoom.p = ggplot(mic.df, aes(ID*1000, Mean)) +
geom_smooth(method = "loess",
se = FALSE, col = "black",
linetype= "solid") +
geom_smooth(data = mic.df, aes(ID*1000, Q025),
method = "loess",
se = FALSE, col = "grey",
linetype= "dashed") +
geom_smooth(data = mic.df, aes(ID*1000, Q975),
method = "loess",
se = FALSE, col = "grey",
linetype= "dashed") +
geom_hline(yintercept = 0,
linetype = "dotted",
col = "red",
size = 1) +
geom_vline(xintercept = 0,
linetype = "dotted",
col = "red",
size = 1) +
xlim(c(0,2000)) +
xlab("Distance to Nearest Neighbor (km)") +
ylab("Logit") +
ggtitle("B") +
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(hjust = 0.5))
grid.arrange(Full.p, Zoom.p, ncol=1)
Threshold Distance
###Max Distance for dispersal #1585
FindZero(mic.df)[1]*1000## [1] 1585Fixed Effects
ggplotmargin Function
Function for plotting the density of the marginal terms.
ggplotmargin = function(x, type, effect, xlab, ylab = "Posterior Density",
int.value = c(value = 0, 5, 95),
color = c("red", "gray", "gray")){
xx = as.data.frame(inla.smarginal(x[[paste("marginals", type, sep=".")]][[effect]]))
out = ggplot(xx, aes(x, y)) + geom_line(size = 1) + ylab(ylab) + xlab(xlab) + 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(hjust = 0.5)) +
if(length(int.value) == 0) int.value = 0
int.value = lapply(int.value, function(x) if(is.character(x))
type.convert(x, as.is = TRUE) else x)
int.value = lapply(int.value, function(x) if(is.character(x))
lapply(strsplit(x, "=")[[1]], type.convert, as.is = TRUE) else x)
nx = names(int.value)
if(!is.null(nx))
for(i in which(nx != "")) int.value[[i]] = list(nx[i], int.value[[i]])
int.value = sapply(int.value, function(x) {
if(is.numeric(x)) xx$x[which.max(cumsum(xx$y)/sum(xx$y) >= as.numeric(x/100))]
else switch(x[[1]],
mean = sum(xx$y*xx$x)/sum(xx$y),
median = xx$x[which.max(cumsum(xx$y)/sum(xx$y) >=.5)],
mode = xx$x[which.max(xx$y)],
value = x[[2]],
zero = 0)})
if(length(color) <= length(int.value)) color = rep(color, length = length(int.value))
for(i in 1:length(int.value)) out = out + geom_vline(xintercept = int.value[i], color = color[i])
out
}Intercept and Hyperparameters
Posterior marginal distribution of intercept, posterior marginal distribution to nominal variance of the random field, and posterior marginal distribution of the practical range.
plotI = ggplotmargin(Model5, type = "fixed",
effect = "intercept0", xlab = "Intercept")
mRF.df = as.data.frame(Model5$marginals.hy[[1]])
plotH1 = ggplot(mRF.df, aes(x, y)) +
geom_smooth(method = "loess",
se = FALSE, col = "black",
span = 0.15,
linetype= "solid") +
xlab(expression(phi)) +
ylab("Posterior Density") +
theme_classic() +
theme(axis.title.y = element_text(size = 18),
axis.title.x = element_text(size = 18))
result.field = inla.spde.result(Model5, "field0", spde)
mRF.df = as.data.frame(result.field$marginals.variance.nominal[[1]])
plotVN = ggplot(mRF.df, aes(x, y)) +
geom_smooth(method = "loess",
se = FALSE, col = "black",
span = 0.05,
linetype= "solid") +
xlab(expression(sigma[x]^2)) +
ylab("Posterior Density") +
theme_classic() +
theme(axis.title.y = element_text(size = 18),
axis.title.x = element_text(size = 18))
mRF.df = as.data.frame(result.field$marginals.range.nominal[[1]])
plotRng = ggplot(mRF.df, aes(x, y)) +
geom_smooth(method = "loess",
se = FALSE, col = "black",
span = 0.05,
linetype= "solid") +
xlab("Practical Range") +
ylab("Posterior Density") +
theme_classic() +
theme(axis.title.y = element_text(size = 18),
axis.title.x = element_text(size = 18))
grid.arrange(plotI, plotH1, plotVN, plotRng, ncol = 2)
Spatial Random Effect
Point estimates and 95% credible intervals for the course random spatial effect.
mRF.df = as.data.frame(Model0$summary.random$field0)[,1:6]
names(mRF.df) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975")
mRF.df = arrange(mRF.df, Mean)
mRF.df$Index = 1:nrow(mRF.df)
SpCor1 = ggplot(mRF.df, aes(Index, Mean)) +
geom_smooth(method = "loess",
col = "black",
span = 0.25,
se = FALSE) +
geom_line(data = mRF.df,
aes(Index, Q025),
col = "grey") +
geom_line(data = mRF.df,
aes(Index, Q975),
col = "grey") +
geom_hline(yintercept = 0,
linetype = "dotted",
col = "red",
size = 1) +
scale_y_continuous(breaks=c(-15,-10,-5, 0, 5, 10, 15),
labels=c("-15","-10","-5", "0", "5", "10", "15")) +
#ylim(-16, 16) +
xlab(NULL) +
ylab("Spatial Field") +
ggtitle("A") +
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(hjust = 0.5))
mRF.df = as.data.frame(Model5$summary.random$field0)[,1:6]
names(mRF.df) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975")
mRF.df = arrange(mRF.df, Mean)
mRF.df$Index = 1:nrow(mRF.df)
SpCor2 = ggplot(mRF.df, aes(Index, Mean)) +
geom_smooth(method = "loess",
col = "black",
span = 0.25,
se = FALSE) +
geom_line(data = mRF.df,
aes(Index, Q025),
col = "grey") +
geom_line(data = mRF.df,
aes(Index, Q975),
col = "grey") +
geom_hline(yintercept = 0,
linetype = "dotted",
col = "red",
size = 1) +
xlab("Index") +
ylab("Spatial Field") +
ggtitle("B") +
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(hjust = 0.5))
grid.arrange(SpCor1, SpCor2, ncol=1)
Density of the Marginal Terms
(Fixed Effects)
results = Model5
results$marginals.fixed$Pop[, 1] = results$marginals.fixed$Pop[, 1]
results$marginals.fixed$CTI[, 1] = results$marginals.fixed$CTI[, 1]
results$marginals.fixed$Pwq[, 1] = results$marginals.fixed$Pwq[, 1]
results$marginals.fixed$CoH[, 1] = results$marginals.fixed$CoH[, 1]
plotI = ggplotmargin(Model5, type = "fixed",
effect = "intercept0", xlab = "Intercept")
plot1 = ggplotmargin(results, type = "fixed",
effect = "Pop", xlab = "Pd")
plot2 = ggplotmargin(results, type = "fixed",
effect = "CTI", xlab = "CTI")
plot3 = ggplotmargin(results, type = "fixed",
effect = "Pwq", xlab = "Pq")
plot4 = ggplotmargin(results, type = "fixed",
effect = "CoH", xlab = "Ic")
Top.grd = grid.arrange(plot1, plot2, plot3, plot4, ncol=2)
Prediction
Global Prediction
Predicting the probability of occurrence globally using “Model1”.
Model.pred = inla(Frm5, #Formula for selcted model
data = inla.stack.data(CStack.global, spde=spde), #
family = "binomial",
verbose = TRUE,
control.fixed = list(prec = 0.1, prec.intercept = 0.1),
control.predictor = list(
A = inla.stack.A(CStack.global),
compute = TRUE,
link = 1),
control.mode = list(restart = TRUE, theta = theta5),
control.results = list(return.marginals.random = FALSE, #Turning off un-needed calculations
return.marginals.predictor = FALSE),
control.compute=list(dic = FALSE, cpo = FALSE, waic = FALSE)) summary(Model.pred) ##
## Call:
## c("inla(formula = Frm5, family = \"binomial\", data = inla.stack.data(CStack.global, ", " spde = spde), verbose = TRUE, control.compute = list(dic = FALSE, ", " cpo = FALSE, waic = FALSE), control.predictor = list(A = inla.stack.A(CStack.global), ", " compute = TRUE, link = 1), control.results = list(return.marginals.random = FALSE, ", " return.marginals.predictor = FALSE), control.fixed = list(prec = 0.1, ", " prec.intercept = 0.1), control.mode = list(restart = TRUE, ", " theta = theta5))")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 1.4447 35909.4810 1.6790 35912.6047
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept0 -9.2230 1.5337 -12.2905 -9.2053 -6.2571 -9.1709 0
## Pop 0.3328 0.1556 0.0441 0.3266 0.6555 0.3138 0
## CTI -0.8209 0.3281 -1.4689 -0.8196 -0.1806 -0.8170 0
## Pwq 0.9616 0.2332 0.5004 0.9626 1.4168 0.9645 0
## CoH 0.7396 0.0951 0.5569 0.7381 0.9307 0.7353 0
##
## Random effects:
## Name Model
## field0 SPDE2 model
## mTcm RW1 model
## NN RW1 model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant 0.975quant mode
## Range for field0 0.4347 0.1536 0.2160 0.4072 0.8102 0.3578
## Stdev for field0 0.9570 0.1937 0.6376 0.9361 1.3949 0.8949
## Precision for mTcm 0.8272 0.4519 0.2684 0.7250 1.9844 0.5591
## Precision for NN 0.5369 0.1186 0.3453 0.5226 0.8089 0.4946
##
## Expected number of effective parameters(std dev): 62.94(8.886)
## Number of equivalent replicates : 413.99
##
## Marginal log-Likelihood: -9393.83
## Posterior marginals for linear predictor and fitted values computedPredicted Distribution
Under Current Conditions(1950 - 2000)
I0 = inla.stack.index(CStack.global, "pred0") #Index to pull results
E0 = Model.pred$summary.linear.predictor[I0$data,] #pull results
Pred.pnts = LygPA
Pred.pnts$Cp = plogis(E0[,"mean"]) #exponentiate
M0p = rasterize(Pred.pnts, Domain.r, "Cp", #Rasterize predicted values
background = NA)
M0p = sum(M0p, Domain.r)
rng = seq(0, 1, 0.001)
cr0 = rev(colorRampPalette(c("maroon", "yellow", "dodgerblue"))(n = 299))
cr = colorRampPalette(c("tan", cr0),
bias = 0.6, space = "rgb")
Cp = levelplot(M0p,
margin = FALSE,
#sub = " Probability of Occurrence",
xlab = " ", ylab = NULL,
col.regions = cr, at = rng,
colorkey = list(at=seq(0, 1, 0.001),
labels=list(cex=1.5),
labels=list(at=c(0, 0.25, 0.5, 0.75, 1),
labels=c("0%", "25%", "50%", "75%", "100%")),
space = "bottom"),
par.settings = list(axis.line = list(col = "black"),
strip.background = list(col = 'transparent'),
strip.border = list(col = 'transparent')),
scales = list(cex = 2))
Cp +
latticeExtra::layer(sp.polygons(Countries, col = "dimgray", cex=2, lwd = 1))
Sum of Linear Components
As an alternative to “analytic” prediction, we explore results from a simple sum of linear components, to include the spatial random field and non-linear random walk results.
Pred.pnts = LygPA #Copy Points
pLoc = cbind(Pred.pnts$Long, Pred.pnts$Lat)
pLoc = inla.mesh.map(pLoc,
projection = "longlat",
inverse = TRUE)
Ap = inla.spde.make.A(mesh1, loc=pLoc) #project new points to mesh to obtain random field
Pred.pnts$RF = drop(Ap %*% Model5$summary.random$field0$mean) #Record the random field at each location
#Get the intercept
Pred.pnts$Int = Model5$summary.fix[1,1]
mydf = as.data.frame(Model5$summary.fixed) #Write results to data frame
#Find closest match for the mTcm RW1 results
Rw.lu = as.numeric(Model5$summary.random$mTcm[,"ID"])
Pred.pnts$RW = sapply(Pred.pnts$mTcm, function(x)which.min(abs(x - Rw.lu)))
#Write rw1 value to look up table and note rank/position
Rw.lu = as.data.frame(Model5$summary.random$mTcm[,"mean"])
names(Rw.lu) = "Mean"
Rw.lu$POS = as.integer(seq(1, 766, 1))
#Use rank/position to pull nearest closest effect size
Pred.pnts$mTcm_lu = as.numeric(with(Rw.lu,
Mean[match(Pred.pnts$RW, POS)]))
#Same as above, but for the NN RW1 result
Rw.lu = as.numeric(Model5$summary.random$NN[,"ID"])
Pred.pnts$RW = sapply(Pred.pnts$NN, function(x)which.min(abs(x - Rw.lu)))
#Write rw1 value to look up table and note rank/position
Rw.lu = as.data.frame(Model5$summary.random$NN[,"mean"])
names(Rw.lu) = "Mean"
Rw.lu$POS = as.integer(seq(1, 2189, 1))
#Use rank/position to pull nearest neighbor effect size
Pred.pnts$NN_lu = as.numeric(with(Rw.lu,
Mean[match(Pred.pnts$RW, POS)]))
#Sum intercept, random effect, random walk variables, and fixed effects.
Pred.pnts@data = Pred.pnts@data %>%
mutate(LP = Int + RF + mTcm_lu + NN_lu +
Pop*Model5$summary.fix[2,1] +
CTI*Model5$summary.fix[3,1] +
Pwq*Model5$summary.fix[4,1] +
CoH*Model5$summary.fix[5,1],
LP = ifelse(is.na(LP), 0, LP),
LPe = plogis(LP)) Compare Predictions
Analytic and Summed Linear Component predictions are compared. First visually, then by correlation test.
Interpretation: Results are comparable. Use of the Summed Linear Components is a reasonable alternative to analytic computation, which can be computationally expensive and prone to numerical issues (because of combined spatial and dispersal random effects).
Global_LP = rasterize(Pred.pnts, Domain.r, "LPe",
background = NA)
Global_LP = sum(Global_LP, Domain.r)
CompPred.stk = stack(M0p, Global_LP)
names(CompPred.stk) = c("A", "B")
rng = seq(0, 1, 0.001)
cr0 = rev(colorRampPalette(c("maroon", "yellow", "dodgerblue"))(n = 299))
cr = colorRampPalette(c("tan", cr0),
bias = 0.7, space = "rgb")
Cp = levelplot(M0p,
margin = FALSE,
sub = " Probability of Occurrence",
xlab = NULL, ylab = NULL,
col.regions = cr, at = rng,
colorkey = list(at=seq(0, 1, 0.001),
labels=list(at=c(0, 0.25, 0.5, 0.75, 1),
labels=c("0%", "25%", "50%", "75%", "100%")),
panel = panel.levelplot, space = "right",
col = cr, par.settings = list(fontsize = list(text = 18))))
Cp +
latticeExtra::layer(sp.polygons(Countries, col = "black", lwd = 1))
Test for Correlation
Testing for correlation between prediction methods.
cor.test(values(M0p), values(Global_LP))##
## Pearson's product-moment correlation
##
## data: values(M0p) and values(Global_LP)
## t = 8220.1, df = 24019, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.9998178 0.9998268
## sample estimates:
## cor
## 0.9998223Prediction for Southeastern United States
Pred.pnts = SECP.pnts #Copy Points
pLoc = cbind(Pred.pnts$Long, Pred.pnts$Lat)
pLoc = inla.mesh.map(pLoc,
projection = "longlat",
inverse = TRUE)
Ap = inla.spde.make.A(mesh1, loc=pLoc)
Pred.pnts$RF = drop(Ap %*% Model5$summary.random$field0$mean)
#Get CoH covariate
Pred.pnts$CoH = drop(Ap %*% ModelR$summary.random$fieldR$mean) #RElative density of Invasive Cohort Richness
Pred.pnts$Int = Model5$summary.fix[1,1]
mydf = as.data.frame(Model5$summary.fixed)
Rw.lu = as.numeric(Model5$summary.random$mTcm[,"ID"])
Pred.pnts$RW = sapply(Pred.pnts$mTcm, function(x)which.min(abs(x - Rw.lu)))
Rw.lu = as.data.frame(Model5$summary.random$mTcm[,"mean"])
names(Rw.lu) = "Mean"
Rw.lu$POS = as.integer(seq(1, 766, 1))
Pred.pnts$mTcm_lu = as.numeric(with(Rw.lu,
Mean[match(Pred.pnts$RW, POS)]))
Rw.lu = as.numeric(Model5$summary.random$NN[,"ID"])
Pred.pnts$RW = sapply(Pred.pnts$NN, function(x)which.min(abs(x - Rw.lu)))
Rw.lu = as.data.frame(Model5$summary.random$NN[,"mean"])
names(Rw.lu) = "Mean"
Rw.lu$POS = as.integer(seq(1, 2189, 1))
Pred.pnts$NN_lu = as.numeric(with(Rw.lu,
Mean[match(Pred.pnts$RW, POS)]))
Pred.pnts@data = Pred.pnts@data %>%
mutate(LP = Int + RF + mTcm_lu + NN_lu +
Pop*Model5$summary.fix[2,1] +
CTI*Model5$summary.fix[3,1] +
Pwq*Model5$summary.fix[4,1] +
CoH*Model5$summary.fix[5,1],
LP = ifelse(is.na(LP), 0, LP),
LPe = plogis(LP))
#Unconstrained Dispersal
UnDisp = Model5$summary.random$NN[1, "mean"]
#Cuurent condition, unconstrained dispersal
Pred.pnts@data = Pred.pnts@data %>%
mutate(LP = Int + RF + mTcm_lu + UnDisp +
Pop*Model5$summary.fix[2,1] +
CTI*Model5$summary.fix[3,1] +
Pwq*Model5$summary.fix[4,1] +
CoH*Model5$summary.fix[5,1],
LP = ifelse(is.na(LP), 0, LP),
LPeC = plogis(LP))
#Changing to Project Climate (Changing mTcm and Pwq)
Rw.lu = as.numeric(Model5$summary.random$mTcm[,"ID"])
Pred.pnts$RWF = sapply(Pred.pnts$mTcm70, function(x)which.min(abs(x - Rw.lu)))
Rw.lu = as.data.frame(Model5$summary.random$mTcm[,"mean"])
names(Rw.lu) = "Mean"
Rw.lu$POS = as.integer(seq(1, 766, 1))
Pred.pnts$FmTcm_lu = as.numeric(with(Rw.lu,
Mean[match(Pred.pnts$RWF, POS)]))
#Constrained Dispersal
Pred.pnts@data = Pred.pnts@data %>%
mutate(LP = Int + RF + FmTcm_lu + NN_lu +
Pop*Model5$summary.fix[2,1] +
CTI*Model5$summary.fix[3,1] +
Pwq70*Model5$summary.fix[4,1] +
CoH*Model5$summary.fix[5,1],
FLP = ifelse(is.na(LP), 0, LP),
FLPe = plogis(LP))
Pred.pnts@data = Pred.pnts@data %>%
mutate(LP = Int + RF + FmTcm_lu + UnDisp +
Pop*Model5$summary.fix[2,1] +
CTI*Model5$summary.fix[3,1] +
Pwq70*Model5$summary.fix[4,1] +
CoH*Model5$summary.fix[5,1],
FLP = ifelse(is.na(LP), 0, LP),
FLPeD = plogis(LP))Predicted Occurence
Current Conditions 1950-2000.
SEUS_LP = rasterize(Pred.pnts, SEStates.r, "LPe",
background = NA)
SEUS_LP = sum(SEUS_LP, SEStates.r)
SEUS_LPb = rasterize(Pred.pnts, SEStates.r, "LPeC",
background = NA)
SEUS_LPb = sum(SEUS_LPb, SEStates.r)
SEUS8570 = rasterize(Pred.pnts, SEStates.r, "FLPe",
background = NA)
SEUS8570 = sum(SEUS8570, SEStates.r)
SEUS8570b = rasterize(Pred.pnts, SEStates.r, "FLPeD",
background = NA)
SEUS8570b = sum(SEUS8570b, SEStates.r)
SEPred.stk = stack(SEUS_LP, SEUS8570, SEUS_LPb, SEUS8570b)
names(SEPred.stk) = c("A", "B", "C", "D")
rng = seq(0, 1, 0.1)
cr0 = rev(colorRampPalette(c("maroon", "yellow", "dodgerblue"))(n = 299))
cr = colorRampPalette(c("tan", cr0),
space = "rgb")
Cp = levelplot(SEPred.stk,
margin = FALSE,
layout=c(2, 2),
#sub = " Probability of Occurrence",
xlab = " ", ylab = NULL,
col.regions = cr, at = rng,
colorkey = list(at=seq(0, 1, 0.001),
labels=list(cex=1.5),
labels=list(at=c(0, 0.25, 0.5, 0.75, 1),
labels=c("0%", "25%", "50%", "75%", "100%")),
space = "bottom"),
par.settings = list(axis.line = list(col = "black"),
strip.background = list(col = 'transparent'),
strip.border = list(col = 'transparent')),
scales = list(cex = 1.5))
Cp +
latticeExtra::layer(sp.polygons(StatesP, col = "dimgray", lwd = 1))
Compare Prediction to Presence Locations
SE.obs = subset(LygPA, OBS == 1)
Cp.us = levelplot(SEUS_LP,
margin = FALSE,
#sub = " Probability of Occurrence",
xlab = " ", ylab = NULL,
col.regions = cr, at = rng,
colorkey = list(at=seq(0, 1, 0.001),
labels=list(cex=1.5),
labels=list(at=c(0, 0.25, 0.5, 0.75, 1),
labels=c("0%", "25%", "50%", "75%", "100%")),
space = "bottom"),
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(StatesP, col = "dimgray", lwd = 1)) +
latticeExtra::layer(sp.polygons(SE.obs, col = "black", pch = 16, cex = 0.5))
Cp.us
Prediction for South East Asia and Oceania
Pred.pnts = SEAsia.pnts #Copy Points
pLoc = cbind(Pred.pnts$Long, Pred.pnts$Lat)
pLoc = inla.mesh.map(pLoc,
projection = "longlat",
inverse = TRUE)
Ap = inla.spde.make.A(mesh1, loc=pLoc)
Pred.pnts$RF = drop(Ap %*% Model5$summary.random$field0$mean)
#Get CoH covariate
Pred.pnts$CoH = drop(Ap %*% ModelR$summary.random$fieldR$mean)
Pred.pnts$Int = Model5$summary.fix[1,1]
mydf = as.data.frame(Model5$summary.fixed)
Rw.lu = as.numeric(Model5$summary.random$mTcm[,"ID"])
Pred.pnts$RW = sapply(Pred.pnts$mTcm, function(x)which.min(abs(x - Rw.lu)))
Rw.lu = as.data.frame(Model5$summary.random$mTcm[,"mean"])
names(Rw.lu) = "Mean"
Rw.lu$POS = as.integer(seq(1, 766, 1))
Pred.pnts$mTcm_lu = as.numeric(with(Rw.lu,
Mean[match(Pred.pnts$RW, POS)]))
Rw.lu = as.numeric(Model5$summary.random$NN[,"ID"])
Pred.pnts$RW = sapply(Pred.pnts$NN, function(x)which.min(abs(x - Rw.lu)))
Rw.lu = as.data.frame(Model5$summary.random$NN[,"mean"])
names(Rw.lu) = "Mean"
Rw.lu$POS = as.integer(seq(1, 2189, 1))
Pred.pnts$NN_lu = as.numeric(with(Rw.lu,
Mean[match(Pred.pnts$RW, POS)]))
Pred.pnts@data = Pred.pnts@data %>%
mutate(LP = Int + RF + mTcm_lu + NN_lu +
Pop*Model5$summary.fix[2,1] +
CTI*Model5$summary.fix[3,1] +
Pwq*Model5$summary.fix[4,1] +
CoH*Model5$summary.fix[5,1],
LP = ifelse(is.na(LP), 0, LP),
LPe = plogis(LP))
Pred.pnts@data = Pred.pnts@data %>%
mutate(LP = Int + RF + mTcm_lu + UnDisp +
Pop*Model5$summary.fix[2,1] +
CTI*Model5$summary.fix[3,1] +
Pwq*Model5$summary.fix[4,1] +
CoH*Model5$summary.fix[5,1],
LP = ifelse(is.na(LP), 0, LP),
LPeC = plogis(LP))
#Changing to Project Climate (Changing mTcm and Pwq)
Rw.lu = as.numeric(Model5$summary.random$mTcm[,"ID"])
Pred.pnts$RWF = sapply(Pred.pnts$mTcm70, function(x)which.min(abs(x - Rw.lu)))
Rw.lu = as.data.frame(Model5$summary.random$mTcm[,"mean"])
names(Rw.lu) = "Mean"
Rw.lu$POS = as.integer(seq(1, 766, 1))
Pred.pnts$FmTcm_lu = as.numeric(with(Rw.lu,
Mean[match(Pred.pnts$RWF, POS)]))
Pred.pnts@data = Pred.pnts@data %>%
mutate(LP = Int + RF + FmTcm_lu + NN_lu +
Pop*Model5$summary.fix[2,1] +
CTI*Model5$summary.fix[3,1] +
Pwq70*Model5$summary.fix[4,1] +
CoH*Model5$summary.fix[5,1],
FLP = ifelse(is.na(LP), 0, LP),
FLPe = plogis(LP))
#Unconstrained dispersal
Pred.pnts@data = Pred.pnts@data %>%
mutate(LP = Int + RF + FmTcm_lu + UnDisp +
Pop*Model5$summary.fix[2,1] +
CTI*Model5$summary.fix[3,1] +
Pwq70*Model5$summary.fix[4,1] +
CoH*Model5$summary.fix[5,1],
FLP = ifelse(is.na(LP), 0, LP),
FLPeD = plogis(LP))Predicted Occurence
SEA_LP = rasterize(Pred.pnts, SEAsia.r, "LPe",
background = NA)
SEA_LP = sum(SEA_LP, SEAsia.r)
SEA_LPb = rasterize(Pred.pnts, SEAsia.r, "LPeC",
background = NA)
SEA_LPb = sum(SEA_LPb, SEAsia.r)
SEA8570 = rasterize(Pred.pnts, SEAsia.r, "FLPe",
background = NA)
SEA8570 = sum(SEA8570, SEAsia.r)
SEA8570b = rasterize(Pred.pnts, SEAsia.r, "FLPeD",
background = NA)
SEA8570b = sum(SEA8570b, SEAsia.r)
SEAPred.stk = stack(SEA_LP, SEA8570, SEA_LPb, SEA8570b)
names(SEAPred.stk) = c("A", "B", "C", "D")
rng = seq(0, 1, 0.1)
cr0 = rev(colorRampPalette(c("maroon", "yellow", "dodgerblue"))(n = 299))
cr = colorRampPalette(c("tan", cr0),
space = "rgb")
Cp = levelplot(SEAPred.stk,
margin = FALSE,
layout=c(2, 2),
#sub = " Probability of Occurrence",
xlab = " ", ylab = NULL,
col.regions = cr, at = rng,
colorkey = list(at=seq(0, 1, 0.001),
labels=list(cex=1.5),
labels=list(at=c(0, 0.25, 0.5, 0.75, 1),
labels=c("0%", "25%", "50%", "75%", "100%")),
space = "bottom"),
par.settings = list(axis.line = list(col = "black"),
strip.background = list(col = 'transparent'),
strip.border = list(col = 'transparent')),
scales = list(cex = 1.25))
Cp +
latticeExtra::layer(sp.polygons(SEAsia.sp, col = "dimgray", lwd = 1))
Compare Prediction to Presence Locations
SE.obs = subset(LygPA, OBS == 1)
Cp.A = levelplot(SEA_LP,
margin = FALSE,
#sub = " Probability of Occurrence",
xlab = " ", ylab = NULL,
col.regions = cr, at = rng,
colorkey = list(at=seq(0, 1, 0.001),
labels=list(cex=1.5),
labels=list(at=c(0, 0.25, 0.5, 0.75, 1),
labels=c("0%", "25%", "50%", "75%", "100%")),
space = "bottom"),
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(SEAsia.sp, col = "dimgray", lwd = 1)) +
latticeExtra::layer(sp.polygons(SE.obs, col = "black", pch = 16, cex = 0.5))
Cp.A
Prediction for Mexico, Central and South America
Pred.pnts = MC.pnts #Copy Points
pLoc = cbind(Pred.pnts$Long, Pred.pnts$Lat)
pLoc = inla.mesh.map(pLoc,
projection = "longlat",
inverse = TRUE)
Ap = inla.spde.make.A(mesh1, loc=pLoc)
Pred.pnts$RF = drop(Ap %*% Model5$summary.random$field0$mean)
#Get CoH covariate
Pred.pnts$CoH = drop(Ap %*% ModelR$summary.random$fieldR$mean)
Pred.pnts$Int = Model5$summary.fix[1,1]
mydf = as.data.frame(Model5$summary.fixed)
Rw.lu = as.numeric(Model5$summary.random$mTcm[,"ID"])
Pred.pnts$RW = sapply(Pred.pnts$mTcm, function(x)which.min(abs(x - Rw.lu)))
Rw.lu = as.data.frame(Model5$summary.random$mTcm[,"mean"])
names(Rw.lu) = "Mean"
Rw.lu$POS = as.integer(seq(1, 766, 1))
Pred.pnts$mTcm_lu = as.numeric(with(Rw.lu,
Mean[match(Pred.pnts$RW, POS)]))
Rw.lu = as.numeric(Model5$summary.random$NN[,"ID"])
Pred.pnts$RW = sapply(Pred.pnts$NN, function(x)which.min(abs(x - Rw.lu)))
Rw.lu = as.data.frame(Model5$summary.random$NN[,"mean"])
names(Rw.lu) = "Mean"
Rw.lu$POS = as.integer(seq(1, 2189, 1))
Pred.pnts$NN_lu = as.numeric(with(Rw.lu,
Mean[match(Pred.pnts$RW, POS)]))
Pred.pnts@data = Pred.pnts@data %>%
mutate(LP = Int + RF + mTcm_lu + NN_lu +
Pop*Model5$summary.fix[2,1] +
CTI*Model5$summary.fix[3,1] +
Pwq*Model5$summary.fix[4,1] +
CoH*Model5$summary.fix[5,1],
LP = ifelse(is.na(LP), 0, LP),
LPe = plogis(LP))
Pred.pnts@data = Pred.pnts@data %>%
mutate(LP = Int + RF + mTcm_lu + UnDisp +
Pop*Model5$summary.fix[2,1] +
CTI*Model5$summary.fix[3,1] +
Pwq*Model5$summary.fix[4,1] +
CoH*Model5$summary.fix[5,1],
LP = ifelse(is.na(LP), 0, LP),
LPeC = plogis(LP))
#Changing to Project Climate (Changing mTcm and Pwq)
Rw.lu = as.numeric(Model5$summary.random$mTcm[,"ID"])
Pred.pnts$RWF = sapply(Pred.pnts$mTcm70, function(x)which.min(abs(x - Rw.lu)))
Rw.lu = as.data.frame(Model5$summary.random$mTcm[,"mean"])
names(Rw.lu) = "Mean"
Rw.lu$POS = as.integer(seq(1, 766, 1))
Pred.pnts$FmTcm_lu = as.numeric(with(Rw.lu,
Mean[match(Pred.pnts$RWF, POS)]))
Pred.pnts@data = Pred.pnts@data %>%
mutate(LP = Int + RF + FmTcm_lu + NN_lu +
Pop*Model5$summary.fix[2,1] +
CTI*Model5$summary.fix[3,1] +
Pwq70*Model5$summary.fix[4,1] +
CoH*Model5$summary.fix[5,1],
FLP = ifelse(is.na(LP), 0, LP),
FLPe = plogis(LP))
#Unconstrained dispersal
Pred.pnts@data = Pred.pnts@data %>%
mutate(LP = Int + RF + FmTcm_lu + UnDisp +
Pop*Model5$summary.fix[2,1] +
CTI*Model5$summary.fix[3,1] +
Pwq70*Model5$summary.fix[4,1] +
CoH*Model5$summary.fix[5,1],
FLP = ifelse(is.na(LP), 0, LP),
FLPeD = plogis(LP))Predicted Occurence
MC_LP = rasterize(Pred.pnts, MC.r, "LPe",
background = NA)
MC_LP = sum(MC_LP, MC.r)
MC_LPb = rasterize(Pred.pnts, MC.r, "LPeC",
background = NA)
MC_LPb = sum(MC_LPb, MC.r)
MC8570 = rasterize(Pred.pnts, MC.r, "FLPe",
background = NA)
MC8570 = sum(MC8570, MC.r)
MC8570b = rasterize(Pred.pnts, MC.r, "FLPeD",
background = NA)
MC8570b = sum(MC8570b, MC.r)
MCPred.stk = stack(MC_LP, MC8570, MC_LPb, MC8570b)
names(MCPred.stk) = c("A", "B", "C", "D")
rng = seq(0, 1, 0.1)
cr0 = rev(colorRampPalette(c("maroon", "yellow", "dodgerblue"))(n = 299))
cr = colorRampPalette(c("tan", cr0),
space = "rgb")
Cp = levelplot(MCPred.stk,
margin = FALSE,
layout=c(2, 2),
#sub = " Probability of Occurrence",
xlab = " ", ylab = NULL,
col.regions = cr, at = rng,
colorkey = list(at=seq(0, 1, 0.001),
labels=list(cex=1.5),
labels=list(at=c(0, 0.25, 0.5, 0.75, 1),
labels=c("0%", "25%", "50%", "75%", "100%")),
space = "bottom"),
par.settings = list(axis.line = list(col = "black"),
strip.background = list(col = 'transparent'),
strip.border = list(col = 'transparent')),
scales = list(cex = 1.25))
Cp +
latticeExtra::layer(sp.polygons(MC.sp, col = "dimgray", lwd = 1))
Session Info:
sessionInfo()## R version 3.3.2 (2016-10-31)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 8.1 x64 (build 9600)
##
## 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] xtable_1.8-2 threejs_0.2.2 dplyr_0.5.0
## [4] dismo_1.1-1 gridExtra_2.2.1 GISTools_0.7-4
## [7] MASS_7.3-45 ggplot2_2.2.0 spdep_0.6-8
## [10] mapproj_1.2-4 maptools_0.8-40 Thermimage_2.2.3
## [13] rgeos_0.3-21 rasterVis_0.41 latticeExtra_0.6-28
## [16] RColorBrewer_1.1-2 lattice_0.20-34 fields_8.4-1
## [19] maps_3.1.1 spam_1.4-0 raster_2.5-8
## [22] reshape2_1.4.2 rgdal_1.2-4 INLA_0.0-1482340637
## [25] Matrix_1.2-7.1 sp_1.2-3 rgl_0.96.0
## [28] knitr_1.15.1
##
## loaded via a namespace (and not attached):
## [1] Rcpp_0.12.8 deldir_0.1-12 zoo_1.7-13
## [4] gtools_3.5.0 assertthat_0.1 rprojroot_1.1
## [7] digest_0.6.10 mime_0.5 R6_2.2.0
## [10] plyr_1.8.4 backports_1.0.4 evaluate_0.10
## [13] coda_0.19-1 lazyeval_0.2.0 gdata_2.17.0
## [16] hexbin_1.27.1 gmodels_2.16.2 rmarkdown_1.2
## [19] labeling_0.3 splines_3.3.2 stringr_1.1.0
## [22] foreign_0.8-67 htmlwidgets_0.8 munsell_0.4.3
## [25] shiny_0.14.2 numDeriv_2016.8-1 httpuv_1.3.3
## [28] base64enc_0.1-3 htmltools_0.3.5 tibble_1.2
## [31] viridisLite_0.1.3 nlme_3.1-128 jsonlite_1.1
## [34] gtable_0.2.0 DBI_0.5-1 magrittr_1.5
## [37] scales_0.4.1 stringi_1.1.2 LearnBayes_2.15
## [40] boot_1.3-18 tools_3.3.2 markdown_0.7.7
## [43] yaml_2.1.14 parallel_3.3.2 colorspace_1.3-2