Phyllotis amicus
Bayesian Hierarchical Modeling
Joint-modeling of environmental, morphometric, and phylogenetic data.
John: jmh09r@my.fsu.edu
date() ## [1] "Tue Oct 24 12:56:53 2017"Options
library(knitr)
library(rgl)
opts_knit$set(verbose = TRUE)
knit_hooks$set(webgl = hook_webgl)Load Packages:
suppressMessages(library(INLA))
suppressMessages(library(rgdal))
suppressMessages(library(reshape2))
suppressMessages(library(raster))
suppressMessages(library(fields))
suppressMessages(library(FactoMineR))
suppressMessages(library(adegenet))
suppressMessages(library(factoextra))
suppressMessages(library(corrplot))
suppressMessages(library(rasterVis))
suppressMessages(library(rgeos))
suppressMessages(library(maptools))
suppressMessages(library(mapproj))
suppressMessages(library(spdep))
suppressMessages(library(ggplot2))
suppressMessages(library(ggthemes))
suppressMessages(library(GISTools))
suppressMessages(library(lattice))
suppressMessages(library(gridExtra))
suppressMessages(library(xtable))
suppressMessages(library(PresenceAbsence))
suppressMessages(library(dismo))
suppressMessages(library(dplyr))
source("RScripts.R")Set Directory
setwd("D:/Phyllotis/PhyProj/PhyloProj")Load Pre-processed Data
load("PrePro082917_8Ld.RData")View Mesh
(Interactive plot: click, drag, and zoom to view)
Constructed during pre-processing.
O.pnts$Ocean = is.na(over(O.pnts, SA.dom2)) #Land Vs Ocean
cr = colorRampPalette(c("tan", "steelblue"),
space = "rgb")
plot(mesh.dom, rgl = TRUE,
col = O.pnts$Ocean,
color.palette = cr,
draw.edges = TRUE,
draw.segments = TRUE,
draw.vertices =FALSE)You must enable Javascript to view this page properly.
Convert to Cartesian Coordinates
Separate point sets for Poisson and Binomial Likelihoods.
locs = cbind(Mod.pnts$Long, Mod.pnts$Lat)
locs = inla.mesh.map(locs,
projection = "longlat",
inverse = TRUE)
A.osi = inla.spde.make.A(mesh.dom, loc=locs)Prior for SPDE Model0
PC Prior.
FE.df = Mod.pnts@data
n.data = dim(Mod.pnts@data)[1]
n.covariates = 5
X = cbind(rep(1,n.data),
FE.df$bio1,
FE.df$bio15,
FE.df$aPET,
FE.df$Pop)
Q = qr.Q(qr(X))
spde1 = inla.spde2.pcmatern(mesh.dom, alpha = 2,
prior.range=c(0.5, 0.9),
prior.sigma=c(1, 0.5),
constr = TRUE,
extraconstr = list(A = as.matrix(t(Q)%*%A.osi),
e = rep(0, n.covariates)))
field1 = inla.spde.make.index("field1",
spde0$n.spde) #spatial indexError Precision
PC Prior.
error.prec = list(hyper=list(theta=list(prior = "pc.prec", param=c(1, 0.01))))Construct Data Stack
FE.df$OBS2 = ifelse(FE.df$nRvTax == "amicus", 1, 0)
osi.pnt = subset(Mod.pnts, FE.df$OBS2 == 1)
sum(FE.df$OBS2)## [1] 24FE.lst = list(c(field1,
list(intercept1 = 1)),
list(bio1 = FE.df[,"bio1"],
bio15 = FE.df[,"bio15"],
aPET = FE.df[,"aPET"],
Pop = FE.df[,"Pop"]))
Stack1 = inla.stack(data = list(OBS = FE.df$OBS2),
A = list(A.osi, 1),
effects = FE.lst,
tag = "osi.1")Spatial Effects Only (Model0)
Including only spatial effect.
Frm0 = OBS ~ -1 + intercept1 +
f(field1, model=spde)
#theta0 = Model0$internal.summary.hyperpar$mean
theta0 = c(0.793335, -1.710971, 1.481912) #previous run
#osilae
Model0 = inla(Frm0,
data = inla.stack.data(Stack1, spde=spde1),
family = "zeroinflatedbinomial1",
verbose = TRUE,
control.family = list(error.prec),
control.fixed = list(prec = 0.1, prec.intercept = 0.001),
control.predictor = list(
A = inla.stack.A(Stack1),
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 = \"zeroinflatedbinomial1\", data = inla.stack.data(Stack1, ", " spde = spde1), verbose = TRUE, control.compute = list(dic = TRUE, ", " cpo = TRUE, waic = TRUE), control.predictor = list(A = inla.stack.A(Stack1), ", " compute = TRUE, link = 1), control.family = list(error.prec), ", " control.results = list(return.marginals.random = TRUE, return.marginals.predictor = TRUE), ", " control.fixed = list(prec = 0.1, prec.intercept = 0.001), ", " control.mode = list(restart = TRUE, theta = theta0))")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 1.4036 186.0160 0.8463 188.2659
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept1 -8.5559 1.6779 -11.7747 -8.597 -5.9472 -9.0413 0
##
## Random effects:
## Name Model
## field1 SPDE2 model
##
## Model hyperparameters:
## mean sd
## zero-probability parameter for zero-inflated binomial_1 0.7213 0.1204
## Range for field1 0.0805 0.0301
## Stdev for field1 3.0931 0.6133
## 0.025quant
## zero-probability parameter for zero-inflated binomial_1 0.4402
## Range for field1 0.0381
## Stdev for field1 2.0553
## 0.5quant
## zero-probability parameter for zero-inflated binomial_1 0.740
## Range for field1 0.075
## Stdev for field1 3.038
## 0.975quant mode
## zero-probability parameter for zero-inflated binomial_1 0.9033 0.7848
## Range for field1 0.1547 0.0654
## Stdev for field1 4.4534 2.9317
##
## Expected number of effective parameters(std dev): 54.49(9.707)
## Number of equivalent replicates : 71.69
##
## Deviance Information Criterion (DIC) ...: 182.45
## Effective number of parameters .........: 24.53
##
## Watanabe-Akaike information criterion (WAIC) ...: 171.63
## Effective number of parameters .................: 12.24
##
## Marginal log-Likelihood: -124.10
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computedPlot Gaussian Random Field
Pred.pnts = Pred.Grid
ModResult = Model0
ModMesh = mesh.dom
pLoc = cbind(Pred.pnts$Long, Pred.pnts$Lat)
pLoc = inla.mesh.map(pLoc,
projection = "longlat",
inverse = TRUE)
Ap = inla.spde.make.A(ModMesh, loc=pLoc)
#Get the random field
Pred.pnts$RF0 = drop(Ap %*% ModResult$summary.random$field1$mean)
GRF_osi0 = rasterize(Pred.pnts, Domain.r, "RF0",
background = NA)
rng = seq(-12.1, 8.94, 0.01)
mCols = c("#A50026", "#D73027", "#F46D43", "#FDAE61", "#E0F3F8", "#ABD9E9", "#74ADD1", "#4575B4","#313695")
#brewer.pal(11, "RdBu")
cr = colorRampPalette(c(rev(mCols)),
bias = 0.8, space = "rgb")
RF0osi = levelplot(GRF_osi0,
margin = FALSE,
xlab = " ", ylab = NULL,
col.regions = cr, at = rng,
colorkey = list(labels=list(cex=1.5),
space = "right"),
par.settings = list(axis.line = list(col = "black"),
strip.background = list(col = 'transparent'),
strip.border = list(col = 'transparent')),
scales = list(cex = 1.5)) +
latticeExtra::layer(sp.polygons(SpineAndes, col = "black", lwd = 1))
RF0osi
Update PC Prior
n.data = dim(Mod.pnts@data)[1]
n.covariates = 5
X = cbind(rep(1,n.data),
FE.df$bio1,
FE.df$bio15,
FE.df$aPET,
FE.df$Pop)
Q = qr.Q(qr(X))
spde1 = inla.spde2.pcmatern(mesh.dom, alpha = 2,
prior.range=c(0.5, 0.9),
prior.sigma=c(1, 0.5),
constr = TRUE,
extraconstr = list(A = as.matrix(t(Q)%*%A.osi),
e = rep(0, n.covariates)))Evaluate Possible Multi-Colinearity
library(perturb)##
## Attaching package: 'perturb'## The following object is masked from 'package:raster':
##
## reclassifyColin.df = FE.df[,94:135] %>%
select(bio1, bio15, aPET, Pop)
CorCov = cor(Colin.df)
CI = colldiag(CorCov)
CI## Condition
## Index Variance Decomposition Proportions
## intercept bio1 bio15 aPET Pop
## 1 1.000 0.002 0.038 0.002 0.003 0.000
## 2 1.303 0.005 0.007 0.025 0.000 0.021
## 3 1.730 0.001 0.001 0.050 0.001 0.060
## 4 16.588 0.992 0.955 0.923 0.996 0.919Base Model (Model1)
Including all fixed and random effects.
Frm1 = OBS ~ -1 + intercept1 +
f(field1, model=spde) +
bio1 + bio15 + aPET + Pop
#theta1 = Model1$internal.summary.hyperpar$mean
theta1 = c(-0.2018038, -0.9674923, 1.3514013)
Model1 = inla(Frm1,
data = inla.stack.data(Stack1, spde=spde1),
family = "zeroinflatedbinomial1",
verbose = TRUE,
control.family = list(error.prec),
control.fixed = list(prec = 0.1, prec.intercept = 0.0001),
control.predictor = list(
A = inla.stack.A(Stack1),
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 = \"zeroinflatedbinomial1\", data = inla.stack.data(Stack1, ", " spde = spde1), verbose = TRUE, control.compute = list(dic = TRUE, ", " cpo = TRUE, waic = TRUE), control.predictor = list(A = inla.stack.A(Stack1), ", " compute = TRUE, link = 1), control.family = list(error.prec), ", " control.results = list(return.marginals.random = TRUE, return.marginals.predictor = TRUE), ", " control.fixed = list(prec = 0.1, prec.intercept = 1e-04))" )
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 14.6864 730.4262 1.0411 746.1537
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept1 -8.1780 0.5435 -9.3199 -8.1516 -7.1852 -8.0984 0
## bio1 2.3544 0.6956 1.0131 2.3463 3.7427 2.3302 0
## bio15 1.4600 0.2597 0.9720 1.4521 1.9929 1.4364 0
## aPET -2.1075 0.7208 -3.4990 -2.1154 -0.6705 -2.1312 0
## Pop 0.0476 0.0252 0.0056 0.0449 0.1041 0.0389 0
##
## Random effects:
## Name Model
## field1 SPDE2 model
##
## Model hyperparameters:
## mean sd
## zero-probability parameter for zero-inflated binomial_1 0.4499 0.0310
## Range for field1 0.3808 0.0236
## Stdev for field1 3.8673 0.1874
## 0.025quant
## zero-probability parameter for zero-inflated binomial_1 0.3817
## Range for field1 0.3434
## Stdev for field1 3.5394
## 0.5quant
## zero-probability parameter for zero-inflated binomial_1 0.4535
## Range for field1 0.3775
## Stdev for field1 3.8521
## 0.975quant mode
## zero-probability parameter for zero-inflated binomial_1 0.5009 0.468
## Range for field1 0.4347 0.367
## Stdev for field1 4.2716 3.806
##
## Expected number of effective parameters(std dev): 32.09(1.27)
## Number of equivalent replicates : 121.72
##
## Deviance Information Criterion (DIC) ...: 198.74
## Effective number of parameters .........: 25.39
##
## Watanabe-Akaike information criterion (WAIC) ...: 184.91
## Effective number of parameters .................: 10.53
##
## Marginal log-Likelihood: -134.85
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computedDensity of the Marginal Terms
(Fixed Effects)
Int.df = as.data.frame(Model1$marginals.fix[[1]])
bio1.df = as.data.frame(Model1$marginals.fix[[2]])
bio15.df = as.data.frame(Model1$marginals.fix[[3]])
aPET.df = as.data.frame(Model1$marginals.fix[[4]])
Pop.df = as.data.frame(Model1$marginals.fix[[5]])
CI.df = as.data.frame(Model1$summary.fixed)[,c(3,5)]
In.plt = ggplot(Int.df , aes(x, y)) +
geom_line(size = 1, col = "black") +
geom_vline(xintercept = CI.df[1,1], col = "gray") +
geom_vline(xintercept = CI.df[1,2], col = "gray") +
geom_vline(xintercept = 0, col = "red") +
xlab(NULL) +
ylab("Density") +
ggtitle("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(size = 20, hjust = 0.5))
bio1.plt = ggplot(bio1.df , aes(x, y)) +
geom_line(size = 1, col = "black") +
geom_vline(xintercept = CI.df[2,1], col = "gray") +
geom_vline(xintercept = CI.df[2,2], col = "gray") +
geom_vline(xintercept = 0, col = "red") +
xlab(NULL) +
ylab(NULL) +
ggtitle("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(size = 20, hjust = 0.5))
bio15.plt = ggplot(bio15.df , aes(x, y)) +
geom_line(size = 1, col = "black") +
geom_vline(xintercept = CI.df[3,1], col = "gray") +
geom_vline(xintercept = CI.df[3,2], col = "gray") +
geom_vline(xintercept = 0, col = "red") +
xlab(NULL) +
ylab("Density") +
ggtitle("C") +
theme_classic() +
theme(axis.text=element_text(size=16),
axis.title.y = element_text(size = 20),
axis.title.x = element_text(size = 20, vjust=-2),
plot.title = element_text(size = 20, hjust = 0.5))
aPET.plt = ggplot(aPET.df , aes(x, y)) +
geom_line(size = 1, col = "black") +
geom_vline(xintercept = CI.df[4,1], col = "gray") +
geom_vline(xintercept = CI.df[4,2], col = "gray") +
geom_vline(xintercept = 0, col = "red") +
xlab(NULL) +
ylab(NULL) +
ggtitle("D") +
theme_classic() +
theme(axis.text=element_text(size=16),
axis.title.y = element_text(size = 20),
axis.title.x = element_text(size = 20, vjust=-2),
plot.title = element_text(size = 20, hjust = 0.5))
Pop.plt = ggplot(Pop.df , aes(x, y)) +
geom_line(size = 1, col = "black") +
geom_vline(xintercept = CI.df[5,1], col = "gray") +
geom_vline(xintercept = CI.df[5,2], col = "gray") +
geom_vline(xintercept = 0, col = "red") +
xlab(NULL) +
ylab("Density") +
ggtitle("E") +
theme_classic() +
theme(axis.text=element_text(size=16),
axis.title.y = element_text(size = 20),
axis.title.x = element_text(size = 20, vjust=-2),
plot.title = element_text(size = 20, hjust = 0.5))
grid.arrange(In.plt, bio1.plt, bio15.plt, aPET.plt, Pop.plt, ncol = 2)
Sum of Linear Components
Pred.pnts = Pred.Grid
ModResult = Model1
ModMesh = mesh.dom
pLoc = cbind(Pred.pnts$Long, Pred.pnts$Lat)
pLoc = inla.mesh.map(pLoc,
projection = "longlat",
inverse = TRUE)
Ap = inla.spde.make.A(ModMesh, loc=pLoc) #project new points to mesh to obtain random field
#Get the random field
Pred.pnts$RF = drop(Ap %*% ModResult$summary.random$field1$mean)
#Get the intercept
Pred.pnts$Int = ModResult$summary.fix[1,1]
#Get fixed effects
mydf = as.data.frame(ModResult$summary.fixed) #Write results to data frame
#Sum intercept, random effect, and fixed effects.
Pred.pnts@data = Pred.pnts@data %>%
mutate(LP = Int + RF +
bio1*ModResult$summary.fix[2,1] +
bio15*ModResult$summary.fix[3,1] +
aPET*ModResult$summary.fix[4,1] +
Pop*ModResult$summary.fix[5,1],
LP = ifelse(is.na(LP), 0, LP),
LPe = plogis(LP)) Plot Predicted Distribution
Pred_osi1 = rasterize(Pred.pnts, Domain.r, "LPe",
background = NA)
rng = seq(0, 1, 0.01)
mCols = brewer.pal(9, "YlOrRd")[2:9]
cr0 = rev(colorRampPalette(rev(mCols))(n = 299))
cr = colorRampPalette(c("tan", cr0),
bias = 2, space = "rgb")
M1oL = levelplot(Pred_osi1,
margin = FALSE,
xlab = " ", ylab = NULL,
col.regions = cr, at = rng,
colorkey = list(labels=list(at=c(0, 0.25, 0.5, 0.75, 1),
labels=c("0%", "25%", "50%", "75%", "100%"), cex=1.5),
labels=list(cex=12),
space = "right"),
par.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(SpineAndes, col = "black", lwd = 1)) +
latticeExtra::layer(sp.polygons(osi.pnt, col = "black", lwd = 0.01, pch=19))
M1op = levelplot(Pred_osi1,
margin = FALSE,
xlab = " ", ylab = NULL,
col.regions = cr, at = rng,
colorkey = list(labels=list(at=c(0, 0.25, 0.5, 0.75, 1),
labels=c("0%", "25%", "50%", "75%", "100%"), cex=1.5),
labels=list(cex=12),
space = "right"),
par.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(SpineAndes, col = "black", lwd = 1))
M1op
M1oL
Spatial Random Effect
Point estimates and 95% credible intervals for the course random spatial effect.
mRF.df = as.data.frame(Model1$summary.random$field1)[,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)
SpCor = 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("Spatial Index") +
ylab("Spatial Field") +
scale_x_continuous(breaks=seq(0, 3400, 680), limits = c(0, 3400)) +
scale_y_continuous(breaks=seq(-15, 15, 5), limits = c(-15, 15)) +
theme_classic() +
theme(axis.text=element_text(size=14),
axis.title.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
plot.title = element_text(hjust = 0.5))
SpCor
Plot Random Field
RF_osi1 = rasterize(Pred.pnts, Domain.r, "RF",
background = NA)
XStack = stack(GRF_osi0, RF_osi1)
names(XStack) = c("A", "B")
rng = seq(-12.1, 8.94, 0.01)
mCols = c("#A50026", "#D73027", "#F46D43", "#FDAE61", "#E0F3F8", "#ABD9E9", "#74ADD1", "#4575B4","#313695")
#brewer.pal(11, "RdBu")
cr = colorRampPalette(c(rev(mCols)),
bias = 0.8, space = "rgb")
RFosi = levelplot(XStack,
margin = FALSE,
xlab = " ", ylab = NULL,
col.regions = cr, at = rng,
colorkey = list(labels=list(cex=1.5),
space = "right"),
par.settings = list(axis.line = list(col = "black"),
strip.background = list(col = 'transparent'),
strip.border = list(col = 'transparent')),
scales = list(cex = 1.5)) +
latticeExtra::layer(sp.polygons(SpineAndes, col = "black", lwd = 1))
RFosi
Non-spatial Model (Model2)
For comparison only.
Frm2 = OBS ~ -1 + intercept1 +
bio1 + bio15 + aPET + Pop
Model2 = inla(Frm2,
data = inla.stack.data(Stack1, spde=spde1),
family = "binomial",
verbose = TRUE,
control.fixed = list(prec = 0.1, prec.intercept = 0.001),
control.predictor = list(
A = inla.stack.A(Stack1),
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(Model2)##
## Call:
## c("inla(formula = Frm2, family = \"binomial\", data = inla.stack.data(Stack1, ", " spde = spde1), verbose = TRUE, control.compute = list(dic = TRUE, ", " cpo = TRUE, waic = TRUE), control.predictor = list(A = inla.stack.A(Stack1), ", " compute = TRUE, link = 1), control.results = list(return.marginals.random = TRUE, ", " return.marginals.predictor = TRUE), control.fixed = list(prec = 0.1, ", " prec.intercept = 0.001))")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 0.3162 4.4422 0.9136 5.6719
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept1 -5.9208 0.3008 -6.5507 -5.9071 -5.3682 -5.8795 0
## bio1 1.5831 0.4725 0.6284 1.5926 2.4841 1.6116 0
## bio15 1.2356 0.1658 0.9092 1.2359 1.5599 1.2365 0
## aPET -1.2577 0.5707 -2.3255 -1.2767 -0.0819 -1.3148 0
## Pop 0.0569 0.0171 0.0213 0.0576 0.0886 0.0590 0
##
## The model has no random effects
##
## The model has no hyperparameters
##
## Expected number of effective parameters(std dev): 4.933(0.00)
## Number of equivalent replicates : 791.74
##
## Deviance Information Criterion (DIC) ...: 250.50
## Effective number of parameters .........: 4.706
##
## Watanabe-Akaike information criterion (WAIC) ...: 250.24
## Effective number of parameters .................: 4.063
##
## Marginal log-Likelihood: -137.78
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computedPlot Non-Spatial Model
Pred.pnts = Pred.Grid
ModResult = Model2
#Get the intercept
Pred.pnts$Int = ModResult$summary.fix[1,1]
#Get fixed effects
mydf = as.data.frame(ModResult$summary.fixed)
Pred.pnts@data = Pred.pnts@data %>%
mutate(LP = Int +
bio1*ModResult$summary.fix[2,1] +
bio15*ModResult$summary.fix[3,1] +
aPET*ModResult$summary.fix[4,1] +
Pop*ModResult$summary.fix[5,1],
LP = ifelse(is.na(LP), 0, LP),
LPe = plogis(LP))
NS_M2 = rasterize(Pred.pnts, Domain.r, "LPe",
background = NA)
rng = seq(0, 1, 0.01)
mCols = brewer.pal(9, "YlOrRd")[2:9]
cr0 = rev(colorRampPalette(rev(mCols))(n = 299))
cr = colorRampPalette(c("tan", cr0),
bias = 1, space = "rgb")
M2ns = levelplot(NS_M2,
margin = FALSE,
xlab = " ", ylab = NULL,
col.regions = cr, at = rng,
colorkey = list(labels=list(at=c(0, 0.25, 0.5, 0.75, 1),
labels=c("0%", "25%", "50%", "75%", "100%"), cex=1.5),
labels=list(cex=12),
space = "right"),
par.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(SpineAndes, col = "black", lwd = 1))
M2ns 
Joint Model
Additional Spatial Indices
Spatial indices are constructed.
field1x = inla.spde.make.index("field1x",
spde1$n.spde)
field2 = inla.spde.make.index("field2",
spde1$n.spde)
field2x = inla.spde.make.index("field2x",
spde1$n.spde)
field3 = inla.spde.make.index("field3",
spde1$n.spde)
field3x = inla.spde.make.index("field3x",
spde1$n.spde)Get Locations with Morphological Data
OBS.pnts2 = subset(OBS.pnts, Source != "phylo")
OBS.pnts2@data = OBS.pnts2@data %>%
filter(nRvTax != "alisosiensis",
nRvTax != "definitus",
nRvTax != "ricardulus")
MapCent = gCentroid(SpineAndes, byid=T)
cnt.df = as.data.frame(MapCent@coords)
cnt.df$id = rownames(cnt.df)
wmap_df = fortify(SpineAndes)## Regions defined for each PolygonsPnts.df = OBS.pnts2@data %>%
select(Long, Lat, nSpp)
names(Pnts.df) = c("Long","Lat","Species")
#myCol = c("darkblue", "darkgreen", "red", "black", "purple", "orange", "maroon", "lightgreen")
myCol = c("black", "red", "green3", "blue", "cyan", "magenta", "orange", "brown")
names(myCol) = levels(factor(Pnts.df$Species))
myShapes = c(16:19, 16:19)
names(myShapes) = levels(factor(Pnts.df$Species))
colScale = scale_colour_manual(name = "Species", values = myCol,
labels = c("P. amicus", "P. andium", "P.darwini", "P. gerbillus",
"P. limatus", "P. magister", "P. osilae","P. xanthopygus"))
ShpScale = scale_shape_manual(name = "Species", values = myShapes,
labels = c("P. amicus", "P. andium", "P.darwini", "P. gerbillus",
"P. limatus", "P. magister", "P. osilae","P. xanthopygus"))
ggplot(wmap_df, aes(long,lat, group=group)) +
geom_polygon(fill = "tan", col="gray21") +
geom_point(data=Pnts.df,
aes(Long, Lat,
group=Pnts.df$Species, fill=NULL,
col = Pnts.df$Species, shape=Pnts.df$Species), size=1) +
geom_text(data=cnt.df, aes(label = id, x = x, y = y, group=id), cex=5, col = "gray21") +
ShpScale +
colScale +
xlab("Longitude") +
ylab("Latitude") +
scale_x_longitude(xmin=-80, xmax=-55, step=10) +
scale_y_latitude(ymin=-50, ymax=10, step=10) +
coord_equal() +
theme(panel.grid.minor = element_blank(),
panel.grid.major = element_blank(),
panel.background = element_blank(),
plot.background = element_rect(fill = "gray80", linetype="solid"),
panel.border = element_blank(),
legend.title = element_text(size=16, face="bold"),
legend.key = element_rect(fill = "gray80", linetype=0),
legend.background = element_rect(fill = "gray80", linetype=0),
axis.title.x = element_text(size=18, face="bold"),
axis.title.y = element_text(size=18, face="bold"),
#axis.text.x = element_text(size=10, face="bold"),
plot.title = element_blank())
#grid.arrange(Map.plt, ncol=2)Project Morphological Specimens to Mesh
locs = cbind(OBS.pnts2$Long, OBS.pnts2$Lat)
locs = inla.mesh.map(locs,
projection = "longlat",
inverse = TRUE)
Aphy = inla.spde.make.A(mesh.dom, loc=locs)
Aphy2 = Aphy #copyLoad and Plot Phylogenetic Tree
Steppan et al., 2007
suppressMessages(library(ape))
suppressMessages(library(ggtree))
MolReap.tree = read.nexus("D:/Phyllotis/RawData/reaprisal_2007/Festschrift_cytb_amicus.unix ML.tre")[[1]]
p1 = ggtree(MolReap.tree, layout="circular", branch.length="none", size=1,
ladderize=TRUE) +
geom_tiplab(aes(subset=(abs(angle) < 90), angle=angle), cex = 3) +
geom_tiplab(aes(subset=(abs(angle) >= 90), angle=angle+180), cex = 3, hjust=1) +
geom_hilight(127, fill="orange", alpha = 0.75) +
geom_hilight(144, fill="black") +
geom_hilight(146, fill="blue", alpha = 0.75) +
geom_hilight(136, fill="red", alpha = 0.75) +
geom_hilight(149, fill="magenta", alpha = 0.75) +
geom_hilight(154, fill="green3", alpha = 0.75) +
geom_hilight(160, fill="brown", alpha = 0.5) +
geom_hilight(176, fill="cyan", alpha = 0.75)
p1 = rotate_tree(p1, 150)
p1
VCV Matrix
Creating a variance covariance matrix from phylogent=etic tree.
VCV2007.mat = vcv.phylo(MolReap.tree, model = "Brownian", corr = FALSE)
VCV2007.df = as.data.frame(VCV2007.mat)
New.testSet = VCV2007.df[!grepl("pseudogene",row.names(VCV2007.df)),
!grepl("pseudogene", colnames(VCV2007.df))]
TSppName = c("amicus", "andium", "chilensis", "darwini", "gerbillus", "limatus", "magister",
"osilae", "posticalis", "rupestris", "vaccarum", "xanthopygus")
VcVr.mat = as.data.frame(matrix(ncol=length(TSppName), nrow=length(TSppName)))
for(i in 1:length(TSppName)){
Tmp1 = colMeans(VCV2007.df[grep(TSppName[i], colnames(VCV2007.df)),])
if( i == 1){Tmp.rows = Tmp1}
else {Tmp.rows = rbind(Tmp.rows, Tmp1)}
}
Tmp.rows = as.data.frame(Tmp.rows)
rownames(Tmp.rows) = TSppName
for(i in 1:length(TSppName)){
Tmp2 = as.data.frame(Tmp.rows[,grep(TSppName[i], colnames(Tmp.rows))])
if(dim(Tmp2)[2] > 1){Tmp2 = rowMeans(Tmp.rows[,grep(TSppName[i], colnames(Tmp.rows))])}
else{Tmp2 = Tmp2}
if( i == 1){Tmp.cols = Tmp2}
else {Tmp.cols = cbind(Tmp.cols, Tmp2)}
}
colnames(Tmp.cols) = TSppName
VcV2007.mat = Tmp.cols
#write.csv(VcV2007.mat, "D:/Phyllotis/PhyProj/PhyloProj/VcV07_082117.csv") #To be used during model fitting Pruned Tree for Figures
KeepSpp = c("P._x._xanthopygus_AY275128", "P._limatus_107615_5C", "P._magister_1804_1813_1894", "P._darwini_LCM2511" ,
"P._gerbillus_ORB58E","P._amicus_10787", "P._andium_5B", "P._osilae_ORB98E")
pruned.tree = drop.tip(MolReap.tree, MolReap.tree$tip.label[-match(KeepSpp, MolReap.tree$tip.label)])
#write.nexus(pruned.tree, file="D:/Phyllotis/PhyProj/PhyloProj/Exec_Mod/RedTRee/Reduced2007.nex")
dd = as.data.frame(pruned.tree$tip.label)
names(dd) = "Label"
dd$Species = 0
dd$Species[grep("P._osilae_ORB98E", dd$Label)] = "P. osilae"
dd$Species[grep("P._andium_5B", dd$Label)] = "P. andium"
dd$Species[grep("P._amicus_10787", dd$Label)] = "P. amicus"
dd$Species[grep("P._gerbillus_ORB58E", dd$Label)] = "P. gerbillus"
dd$Species[grep("P._magister_1804_1813_1894", dd$Label)] = "P. magister"
dd$Species[grep("P._darwini_LCM2511", dd$Label)] = "P. darwini"
dd$Species[grep("P._x._xanthopygus_AY275128", dd$Label)] = "P. xanthopygus"
dd$Species[grep("P._limatus_107615_5C", dd$Label)] = "P. limatus"Plot Pruned Tree
pruned.tree$tip.label = dd$Species
myCol = c("orange", "red", "black", "blue", "magenta", "green3", "brown", "cyan")
RedTree = ggtree(pruned.tree , ladderize=TRUE) +
geom_tiplab(cex = 4, hjust = -.2) + ggplot2::xlim(0, 0.26)
RedTree
2007 Variance Covariance Matrix
Loading and labeling matrix for model inclusion.
FE2.df = OBS.pnts2@data
VCV.mat2 = read.csv("VcV07_082117.csv", header = TRUE, sep=",", row.names = 1)
CNames = colnames(VCV.mat2)
colnames(VCV.mat2) = 1:12
rownames(VCV.mat2) = 1:12
VCV.mat2 = as.matrix(VCV.mat2)
FE2.df$nSpp = ifelse(FE2.df$nRvTax == "fulvescens", "darwini",
ifelse(FE2.df$nRvTax == "tucumanus", "osilae",
FE2.df$nRvTax))
LuTable = as.data.frame(cbind(1:12, CNames))
FE2.df$VcVx = with(LuTable,
V1[match(FE2.df$nSpp,
CNames)])
FE2.df$VcVx = as.numeric(as.character(FE2.df$VcVx))Sex and Age
Cleaning up the age and sex designations.
FE2.df$Sex2 = factor(
as.character(
ifelse(FE2.df$sex == "f" | FE2.df$sex == "F", "Female",
ifelse(FE2.df$sex == "m" | FE2.df$sex == "M", "Male", "Unk"))))
FE2.df$Sex2[is.na(FE2.df$Sex2)] = "Unk"
FE2.df$Age2 = as.factor(FE2.df$age)Repeated Measures Effect
In many instances, multiple field specimens are collected from the same geographic location. Sometimes they are of the same species, other times they are different species. Here, we identify all unique location and species combinations in the study domain.
FE2.df = transform(FE2.df, Cluster_ID = as.numeric(interaction(Long, Lat, nRvTax, drop=TRUE)))
FE2.df$Cluster_ID = as.factor(FE2.df$Cluster_ID)Morphological Stack 1
FE2.lst = list(c(field2,
list(intercept2 = 1)),
list(bio3 = FE2.df[,"bio3"],
bio19 = FE2.df[,"bio19"],
VcVx = FE2.df[,"VcVx"],
Sex = FE2.df[,"Sex2"],
Age = FE2.df[,"Age2"],
Clust = FE2.df[,"Cluster_ID"]))
M2.stk = inla.stack(data = list(M2 = FE2.df$Morph3),
A = list(Aphy, 1),
effects = FE2.lst,
tag = "phy0")
Phy.stk = inla.stack(data = list(Y = cbind(FE2.df$Morph3, NA, NA)), #three elements
A = list(Aphy, 1),
effects = FE2.lst,
tag = "phy0")Morphological Stack 2
FE3.lst = list(c(field3, field3x,
list(intercept3 = 1)),
list(bio2 = FE2.df[,"bio2"],
VcVx = FE2.df[,"VcVx"],
Sex = FE2.df[,"Sex2"],
Age = FE2.df[,"Age2"],
Clust = FE2.df[,"Cluster_ID"]))
M3.stk = inla.stack(data = list(M3 = FE2.df$Morph4),
A = list(Aphy2, 1),
effects = FE3.lst,
tag = "phy2")
Phy.stk2 = inla.stack(data = list(Y = cbind(NA, FE2.df$Morph4, NA)), #three elements
A = list(Aphy2, 1),
effects = FE3.lst,
tag = "phy2")Modify Data Stack from Model1
FE.df = Mod.pnts@data
FE.df$OBS2 = ifelse(FE.df$nRvTax == "amicus", 1, 0)
FE.lst = list(c(field1, field1x, field2x,
list(intercept1 = 1)),
list(bio1 = FE.df[,"bio1"],
bio15 = FE.df[,"bio15"],
aPET = FE.df[,"aPET"],
Pop = FE.df[,"Pop"]))
Stack1 = inla.stack(data = list(Y = cbind(NA, NA, FE.df$OBS2)), #three elements
A = list(A.osi, 1),
effects = FE.lst,
tag = "osi.1j")Combine Three Data Stacks
J.Stack = inla.stack(Phy.stk, Phy.stk2, Stack1)Update PC Prior
n.data = dim(FE.df)[1]
n.covariates = 5
X = cbind(rep(1,n.data),
FE.df$bio1,
FE.df$bio15,
FE.df$aPET,
FE.df$Pop)
Q = qr.Q(qr(X))
spde1 = inla.spde2.pcmatern(mesh.dom,
prior.range=c(0.5, 0.9),
prior.sigma=c(1, 0.5),
constr = TRUE,
extraconstr = list(A = as.matrix(t(Q)%*%A.osi),
e = rep(0, n.covariates))) Model for Morphological 3rd PCA
(Spatial Effect Only)
Frm.M2A = M2 ~ -1 + intercept2 +
f(field2, model=spde1)
#theta2A = Model.M2A$internal.summary.hyperpar$mean
theta2A = c(-0.6060693, -2.8078510, 0.5193457)
Model.M2A = inla(Frm.M2A,
data = inla.stack.data(M2.stk,
spde1=spde1,
VcV = VCV.mat2),
family = "gaussian",
verbose = TRUE,
control.family = list(error.prec),
control.fixed = list(prec = 0.001, prec.intercept = 0.0001),
control.predictor = list(
A = inla.stack.A(M2.stk),
compute = TRUE,
link = 1),
control.mode = list(restart = TRUE, theta = theta2A),
control.results = list(return.marginals.random = TRUE,
return.marginals.predictor = TRUE),
control.compute=list(dic = TRUE, cpo = TRUE, waic = TRUE)) summary(Model.M2A)##
## Call:
## c("inla(formula = Frm.M2A, family = \"gaussian\", data = inla.stack.data(M2.stk, ", " spde1 = spde1, VcV = VCV.mat2), verbose = TRUE, control.compute = list(dic = TRUE, ", " cpo = TRUE, waic = TRUE), control.predictor = list(A = inla.stack.A(M2.stk), ", " compute = TRUE, link = 1), control.family = list(error.prec), ", " 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 = theta2A))")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 1.5331 8.6141 0.6907 10.8379
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept2 0.2028 0.0986 0.0147 0.2002 0.4067 0.1961 0
##
## Random effects:
## Name Model
## field2 SPDE2 model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant
## Precision for the Gaussian observations 0.5451 0.0160 0.5146 0.5447
## Range for field2 0.0664 0.0193 0.0383 0.0631
## Stdev for field2 1.7585 0.2719 1.3113 1.7270
## 0.975quant mode
## Precision for the Gaussian observations 0.5774 0.5438
## Range for field2 0.1131 0.0568
## Stdev for field2 2.3743 1.6560
##
## Expected number of effective parameters(std dev): 119.14(6.856)
## Number of equivalent replicates : 21.21
##
## Deviance Information Criterion (DIC) ...: 8830.66
## Effective number of parameters .........: 121.14
##
## Watanabe-Akaike information criterion (WAIC) ...: 8825.15
## Effective number of parameters .................: 105.97
##
## Marginal log-Likelihood: -4516.07
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computedAdding covariates
MorphPC = list(theta=list(prior = "normal", param=c(0, 10)))
MorphPC2 = list(theta=list(prior = "normal", param=c(1, 70)))
MorphPC3 = list(theta=list(prior = "normal", param=c(0, 20)))
Frm.M2B = M2 ~ -1 + intercept2 +
f(field2, model=spde1) +
f(VcVx, model="generic0",
Cmatrix = VcV,
hyper = MorphPC3,
rankdef=1, constr=TRUE,
diagonal=1e-05) +
f(Sex, model="iid",
hyper = MorphPC) +
f(Age, model="iid",
hyper = MorphPC) +
f(Clust, model="iid",
hyper = MorphPC2) +
bio3 + bio19
#theta2B = Model.M2B$internal.summary.hyperpar$mean
theta2B = c(-0.22683311, -2.06702227, -0.20763410, 0.13541314, 0.08782475, 0.08322901, 0.80024507)
Model.M2B = inla(Frm.M2B,
data = inla.stack.data(M2.stk,
spde1=spde1,
VcV = VCV.mat2),
family = "gaussian",
verbose = TRUE,
control.family = list(error.prec),
control.fixed = list(prec = 0.001, prec.intercept = 0.0001),
control.predictor = list(
A = inla.stack.A(M2.stk),
compute = TRUE,
link = 1),
#control.mode = list(restart = TRUE, theta = theta2B),
control.results = list(return.marginals.random = TRUE,
return.marginals.predictor = TRUE),
control.compute=list(dic = TRUE, cpo = TRUE, waic = TRUE)) summary(Model.M2B)##
## Call:
## c("inla(formula = Frm.M2B, family = \"gaussian\", data = inla.stack.data(M2.stk, ", " spde1 = spde1, VcV = VCV.mat2), verbose = TRUE, control.compute = list(dic = TRUE, ", " cpo = TRUE, waic = TRUE), control.predictor = list(A = inla.stack.A(M2.stk), ", " compute = TRUE, link = 1), control.family = list(error.prec), ", " 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.7644 103.0648 1.0937 105.9229
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept2 -0.0391 0.7985 -1.6173 -0.0379 1.5304 -0.0354 0
## bio3 0.3268 0.1422 0.0452 0.3275 0.6039 0.3292 0
## bio19 -0.6531 0.2741 -1.1982 -0.6519 -0.1157 -0.6500 0
##
## Random effects:
## Name Model
## field2 SPDE2 model
## VcVx Generic0 model
## Sex IID model
## Age IID model
## Clust IID model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant
## Precision for the Gaussian observations 0.7901 0.0244 0.7430 0.7898
## Range for field2 0.1331 0.0866 0.0420 0.1095
## Stdev for field2 0.8210 0.2528 0.4676 0.7726
## Precision for VcVx 1.0935 0.2287 0.7136 1.0701
## Precision for Sex 1.1585 0.3728 0.5887 1.1050
## Precision for Age 1.1535 0.3658 0.5926 1.1015
## Precision for Clust 2.7067 0.2569 2.2421 2.6923
## 0.975quant mode
## Precision for the Gaussian observations 0.8392 0.7893
## Range for field2 0.3644 0.0790
## Stdev for field2 1.4508 0.6831
## Precision for VcVx 1.6089 1.0247
## Precision for Sex 2.0411 1.0052
## Precision for Age 2.0180 1.0045
## Precision for Clust 3.2502 2.6613
##
## Expected number of effective parameters(std dev): 269.53(10.56)
## Number of equivalent replicates : 9.376
##
## Deviance Information Criterion (DIC) ...: 8044.89
## Effective number of parameters .........: 272.09
##
## Watanabe-Akaike information criterion (WAIC) ...: 8060.01
## Effective number of parameters .................: 254.50
##
## Marginal log-Likelihood: -4155.71
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computedCompare Morph3 GRFs
pLocs = cbind(Pred.Grid$Long, Pred.Grid$Lat)
pLocs = inla.mesh.map(pLocs,
projection = "longlat",
inverse = TRUE)
Ap = inla.spde.make.A(mesh.dom, loc=pLocs)
Pred.Grid$M2A = drop(Ap %*% Model.M2A$summary.random$field2$mean)
Pred.Grid$M2B = drop(Ap %*% Model.M2B$summary.random$field2$mean)
M2A.r = rasterize(Pred.Grid, Domain.r, "M2A",
background = NA)
M2B.r = rasterize(Pred.Grid, Domain.r, "M2B",
background = NA)
M2.rstk = stack(M2A.r, M2B.r)
names(M2.rstk) = c("A","B")
rng = seq(-2.76, 3.05, 0.001)
mCols = c("#A50026", "#D73027", "#F46D43", "#FDAE61",
"#E0F3F8", "#ABD9E9", "#74ADD1", "#4575B4","#313695")
cr = colorRampPalette(c(rev(mCols)),
bias=1.1,space = "rgb")
M2.p = levelplot(M2.rstk,
margin = FALSE,
xlab = " ", ylab = NULL,
col.regions = cr, at = rng,
colorkey = list(labels=list(cex=1.5),
space = "right"),
par.strip.text = list(fontface='bold', cex=1),
par.settings = list(axis.line = list(col = "black"),
strip.background = list(col = 'transparent'),
strip.border = list(col = 'transparent')),
scales = list(cex = 1.5)) +
latticeExtra::layer(sp.polygons(SpineAndes, col = "black", lwd = 1))
M2.p
Phylogenetic Effect
Labs = CNames
PhySig.df2 = as.data.frame(Model.M2B$summary.random$VcVx)
names(PhySig.df2) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975", "mode", "kld")
SppLabs = c("P. amicus", "P. andium", "P. x. chilensis", "P.darwini", "P. gerbillus",
"P. limatus", "P. magister", "P. osilae", "P. x. posticalis",
"P. x. rupestris", "P. x. vaccarum", "P. x. xanthopygus")
PhySig.df2$Spp = CNames
PE.m2 = ggplot(PhySig.df2, aes(x=Spp, y=Mean)) +
geom_point(size=2, pch=19, col = "red") +
geom_linerange(aes(ymin=Q025, ymax=Q975), colour="black") +
geom_hline(yintercept = 0,
linetype = "dotted",
colour = "red",
size = 1) +
scale_x_discrete(labels=SppLabs) +
theme_classic() +
xlab(NULL) +
ylab("Morphology PC-3") +
theme(axis.title.y = element_text(face="bold", size=14),
axis.text.y = element_text(face="bold", size=12),
axis.text.x = element_text(face="bold", size=14, vjust=1, hjust=1, angle=45))
PE.m2
Repeated Measures Effect
PhySig.df2 = as.data.frame(Model.M2B$summary.random$Clust)
names(PhySig.df2) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975", "mode", "kld")
PhySig.df2$PSigL = ifelse(PhySig.df2$Q025>0 & PhySig.df2$Q975>0, 1, 0)
PhySig.df2$PSigH = ifelse(PhySig.df2$Q025<0 & PhySig.df2$Q975<0, 1, 0)
PhySig.df3 = PhySig.df2 %>%
filter(PSigL == 1 | PSigH == 1)
RM.m2 = ggplot(PhySig.df2, aes(x=ID, y=Mean)) +
geom_point(size=2, pch=1, col = "gray75") +
geom_linerange(aes(ymin=Q025, ymax=Q975), colour="gray75") +
geom_point(data=PhySig.df3, aes(x=ID, y=Mean),
size=2, pch=19, col = "red") +
geom_linerange(data=PhySig.df3, aes(ymin=Q025, ymax=Q975), colour="black") +
geom_hline(yintercept = 0,
linetype = "dotted",
colour = "red",
size = 0.75) +
theme_classic() +
xlab("Location") +
ylab("Morphology PC-3") +
theme(axis.title.y = element_text(face="bold", size=14),
axis.title.x = element_text(face="bold", size=14),
axis.text.y = element_text(face="bold", size=12),
axis.ticks.x=element_blank(),
axis.text.x = element_blank())
RM.m2
Sex Effect
PhySig.df2 = as.data.frame(Model.M2B$summary.random$Sex)
names(PhySig.df2) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975", "mode", "kld")
PhySig.df2$Spp = c("Female", "Male", "Unknown")
SE.m2 = ggplot(PhySig.df2, aes(x=Spp, y=Mean)) +
geom_point(size=2, pch=19, col = "red") +
geom_linerange(aes(ymin=Q025, ymax=Q975), colour="black") +
geom_hline(yintercept = 0,
linetype = "dotted",
colour = "red",
size = 1) +
theme_classic() +
xlab(NULL) +
ylab("Morphology PC-3") +
theme(axis.title.y = element_text(face="bold", size=14),
axis.text.y = element_text(face="bold", size=12),
axis.text.x = element_text(face="bold", size=14, vjust=0.5, hjust=0.5))
SE.m2
Age Effect
PhySig.df2 = as.data.frame(Model.M2B$summary.random$Age)
names(PhySig.df2) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975", "mode", "kld")
PhySig.df2$Spp = c("2 yrs", "2.5 yrs", "3 yrs", "4 yrs")
AE.m2 = ggplot(PhySig.df2, aes(x=Spp, y=Mean)) +
geom_point(size=2, pch=19, col = "red") +
geom_linerange(aes(ymin=Q025, ymax=Q975), colour="black") +
geom_hline(yintercept = 0,
linetype = "dotted",
colour = "red",
size = 1) +
theme_classic() +
xlab(NULL) +
ylab("Morphology PC-3") +
theme(axis.title.y = element_text(face="bold", size=14),
axis.text.y = element_text(face="bold", size=12),
axis.text.x = element_text(face="bold", size=14))
AE.m2
Map of PCA3
Pred.Grid$M2C = drop(Ap %*% Model.M2B$summary.random$field2$mean) +
Model.M2B$summary.fixed[1,1] +
Pred.Grid$bio3*Model.M2B$summary.fixed[2,1] +
Pred.Grid$bio19*Model.M2B$summary.fixed[3,1] +
mean(Model.M2B$summary.random$VcVx$mean) +
mean(Model.M2B$summary.random$Sex$mean) +
mean(Model.M2B$summary.random$Age$mean) +
mean(Model.M2B$summary.random$Clust$mean)
M2C.r = rasterize(Pred.Grid, Domain.r, "M2C",
background = NA)
mCols = brewer.pal(9, "BrBG")[-c(5,6,7)]
cr0 = rev(colorRampPalette(rev(mCols))(n = 500))
cr = colorRampPalette(cr0, bias=0.7, space = "rgb")
rng = seq(-4.51, 1.77, 0.001)
M2.p = levelplot(M2C.r,
margin = FALSE,
xlab = " ", ylab = NULL,
col.regions = cr, at = rng,
colorkey = list(labels=list(cex=1.5),
space = "right"),
main=list("A", side=1, line=0.5),
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.5)) +
latticeExtra::layer(sp.polygons(SpineAndes, col = "black", lwd = 1))
M2.p
Model for Morphological 4th PCA
(Spatial Effect Only)
Frm.M3A = M3 ~ -1 + intercept3 +
f(field3, model=spde1)
#theta3A = Model.M3A$internal.summary.hyperpar$mean
theta3A = c(-0.02949533, -3.59886349, -0.31498436)
Model.M3A = inla(Frm.M3A,
data = inla.stack.data(M3.stk,
spde1=spde1,
VcV = VCV.mat2),
family = "gaussian",
verbose = TRUE,
control.family = list(error.prec),
control.fixed = list(prec = 0.001, prec.intercept = 0.0001),
control.predictor = list(
A = inla.stack.A(M3.stk),
compute = TRUE,
link = 1),
control.mode = list(restart = TRUE, theta = theta3A),
control.results = list(return.marginals.random = TRUE,
return.marginals.predictor = TRUE),
control.compute=list(dic = TRUE, cpo = TRUE, waic = TRUE)) summary(Model.M3A)##
## Call:
## c("inla(formula = Frm.M3A, family = \"gaussian\", data = inla.stack.data(M3.stk, ", " spde1 = spde1, VcV = VCV.mat2), verbose = TRUE, control.compute = list(dic = TRUE, ", " cpo = TRUE, waic = TRUE), control.predictor = list(A = inla.stack.A(M3.stk), ", " compute = TRUE, link = 1), control.family = list(error.prec), ", " 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 = theta3A))")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 1.2882 19.3338 0.7319 21.3539
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept3 -0.0405 0.0424 -0.1244 -0.0403 0.0424 -0.04 0
##
## Random effects:
## Name Model
## field3 SPDE2 model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant
## Precision for the Gaussian observations 0.9714 0.0285 0.9162 0.9711
## Range for field3 0.0284 0.0066 0.0178 0.0276
## Stdev for field3 0.7335 0.0675 0.6101 0.7303
## 0.975quant mode
## Precision for the Gaussian observations 1.0285 0.9708
## Range for field3 0.0434 0.0261
## Stdev for field3 0.8753 0.7239
##
## Expected number of effective parameters(std dev): 110.68(7.23)
## Number of equivalent replicates : 22.83
##
## Deviance Information Criterion (DIC) ...: 7364.40
## Effective number of parameters .........: 113.02
##
## Watanabe-Akaike information criterion (WAIC) ...: 7371.16
## Effective number of parameters .................: 110.31
##
## Marginal log-Likelihood: -3756.91
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computedCovariates added.
MorphPC = list(theta=list(prior = "normal", param=c(0, 10)))
MorphPC2 = list(theta=list(prior = "normal", param=c(1, 70)))
MorphPC3 = list(theta=list(prior = "normal", param=c(0, 20)))
error.prec = list(hyper=list(theta=list(prior='pc.prec', param=c(1, 0.01))))
Frm.M3B = M3 ~ -1 + intercept3 +
f(field3, model=spde1) +
f(VcVx, model="generic0",
Cmatrix = VcV,
hyper = MorphPC3,
rankdef=1, constr=TRUE,
diagonal=1e-05) +
f(Sex, model="iid",
hyper = MorphPC) +
f(Age, model="iid",
hyper = MorphPC) +
f(Clust, model="iid",
hyper = MorphPC2) + bio2
Model.M3B = inla(Frm.M3B,
data = inla.stack.data(M3.stk,
spde1=spde1,
VcV = VCV.mat2),
family = "gaussian",
verbose = TRUE,
control.family = list(error.prec),
control.fixed = list(prec = 0.001, prec.intercept = 0.0001),
control.predictor = list(
A = inla.stack.A(M3.stk),
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(Model.M3B)##
## Call:
## c("inla(formula = Frm.M3B, family = \"gaussian\", data = inla.stack.data(M3.stk, ", " spde1 = spde1, VcV = VCV.mat2), verbose = TRUE, control.compute = list(dic = TRUE, ", " cpo = TRUE, waic = TRUE), control.predictor = list(A = inla.stack.A(M3.stk), ", " compute = TRUE, link = 1), control.family = list(error.prec), ", " 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.6750 83.1344 1.2184 86.0279
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept3 0.2613 0.7759 -1.2717 0.2618 1.7899 0.2630 0
## bio2 -0.1082 0.0550 -0.2175 -0.1078 -0.0013 -0.1067 0
##
## Random effects:
## Name Model
## field3 SPDE2 model
## VcVx Generic0 model
## Sex IID model
## Age IID model
## Clust IID model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant
## Precision for the Gaussian observations 1.330 0.0411 1.2496 1.3295
## Range for field3 0.545 0.5836 0.0371 0.3694
## Stdev for field3 1.227 0.9910 0.1587 0.9693
## Precision for VcVx 1.321 0.2946 0.8275 1.2928
## Precision for Sex 1.147 0.3837 0.5934 1.0793
## Precision for Age 1.153 0.3729 0.6019 1.0921
## Precision for Clust 2.861 0.2426 2.4279 2.8440
## 0.975quant mode
## Precision for the Gaussian observations 1.412 1.3296
## Range for field3 2.103 0.1007
## Stdev for field3 3.819 0.4576
## Precision for VcVx 1.980 1.2397
## Precision for Sex 2.081 0.9589
## Precision for Age 2.052 0.9814
## Precision for Clust 3.381 2.8038
##
## Expected number of effective parameters(std dev): 309.21(10.07)
## Number of equivalent replicates : 8.172
##
## Deviance Information Criterion (DIC) ...: 6769.25
## Effective number of parameters .........: 311.84
##
## Watanabe-Akaike information criterion (WAIC) ...: 6781.89
## Effective number of parameters .................: 283.03
##
## Marginal log-Likelihood: -3535.59
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computedCompare Morph4 GRFs
pLocs = cbind(Pred.Grid$Long, Pred.Grid$Lat)
pLocs = inla.mesh.map(pLocs,
projection = "longlat",
inverse = TRUE)
Ap = inla.spde.make.A(mesh.dom, loc=pLocs)
Pred.Grid$M3A = drop(Ap %*% Model.M3A$summary.random$field3$mean)
Pred.Grid$M3B = drop(Ap %*% Model.M3B$summary.random$field3$mean)
M3A.r = rasterize(Pred.Grid, Domain.r, "M3A",
background = NA)
M3B.r = rasterize(Pred.Grid, Domain.r, "M3B",
background = NA)
M3.rstk = stack(M3A.r, M3B.r)
names(M3.rstk) = c("A","B")
rng = seq(-1.73, 1.77, 0.001)
mCols = c("#A50026", "#D73027", "#F46D43", "#FDAE61",
"#E0F3F8", "#ABD9E9", "#74ADD1", "#4575B4","#313695")
cr = colorRampPalette(c(rev(mCols)),
bias=1.1,space = "rgb")
M3.p = levelplot(M3.rstk,
margin = FALSE,
xlab = " ", ylab = NULL,
col.regions = cr, at = rng,
colorkey = list(labels=list(cex=1.5),
space = "right"),
par.strip.text = list(fontface='bold', cex=1),
par.settings = list(axis.line = list(col = "black"),
strip.background = list(col = 'transparent'),
strip.border = list(col = 'transparent')),
scales = list(cex = 1.5)) +
latticeExtra::layer(sp.polygons(SpineAndes, col = "black", lwd = 1))
M3.p
Phylogenetic Effect for Morph4
Labs = CNames
PhySig.df2 = as.data.frame(Model.M3B$summary.random$VcVx)
names(PhySig.df2) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975", "mode", "kld")
PhySig.df2$Spp = CNames
PE.m3 = ggplot(PhySig.df2, aes(x=Spp, y=Mean)) +
geom_point(size=2, pch=19, col = "red") +
geom_linerange(aes(ymin=Q025, ymax=Q975), colour="black") +
geom_hline(yintercept = 0,
linetype = "dotted",
colour = "red",
size = 1) +
scale_x_discrete(labels=SppLabs) +
theme_classic() +
xlab(NULL) +
ylab("Morphology PC-4") +
theme(axis.title.y = element_text(face="bold", size=14),
axis.text.y = element_text(face="bold", size=12),
axis.text.x = element_text(face="bold", size=14, vjust=1, hjust=1, angle=45))
PE.m3
#grid.arrange(PE.m2, PE.m3, ncol=1)Repeated Measures Effect (Morph4)
PhySig.df2 = as.data.frame(Model.M3B$summary.random$Clust)
names(PhySig.df2) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975", "mode", "kld")
PhySig.df2$PSigL = ifelse(PhySig.df2$Q025>0 & PhySig.df2$Q975>0, 1, 0)
PhySig.df2$PSigH = ifelse(PhySig.df2$Q025<0 & PhySig.df2$Q975<0, 1, 0)
PhySig.df3 = PhySig.df2 %>%
filter(PSigL == 1 | PSigH == 1)
RM.m3 = ggplot(PhySig.df2, aes(x=ID, y=Mean)) +
geom_point(size=2, pch=1, col = "gray75") +
geom_linerange(aes(ymin=Q025, ymax=Q975), colour="gray75", size = 0.5) +
geom_point(data=PhySig.df3, aes(x=ID, y=Mean),
size=2, pch=19, col = "red") +
geom_linerange(data=PhySig.df3, aes(ymin=Q025, ymax=Q975), colour="black") +
geom_hline(yintercept = 0,
linetype = "dotted",
colour = "red",
size = 0.75) +
theme_classic() +
xlab("Location") +
ylab("Morphology PC-4") +
theme(axis.title.y = element_text(face="bold", size=14),
axis.title.x = element_text(face="bold", size=14),
axis.text.y = element_text(face="bold", size=12),
axis.ticks.x=element_blank(),
axis.text.x = element_blank())
RM.m3
#grid.arrange(RM.m2, RM.m3, ncol=1)Sex Effect
PhySig.df2 = as.data.frame(Model.M3B$summary.random$Sex)
names(PhySig.df2) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975", "mode", "kld")
PhySig.df2$Spp = c("Female", "Male", "Unknown")
SE.m3 = ggplot(PhySig.df2, aes(x=Spp, y=Mean)) +
geom_point(size=2, pch=19, col = "red") +
geom_linerange(aes(ymin=Q025, ymax=Q975), colour="black") +
geom_hline(yintercept = 0,
linetype = "dotted",
colour = "red",
size = 1) +
theme_classic() +
xlab(NULL) +
ylab("Morphology PC-4") +
theme(axis.title.y = element_text(face="bold", size=14),
axis.text.y = element_text(face="bold", size=12),
axis.text.x = element_text(face="bold", size=14))
SE.m3
#grid.arrange(SE.m2, SE.m3, ncol=1)Age Effect
PhySig.df2 = as.data.frame(Model.M3B$summary.random$Age)
names(PhySig.df2) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975", "mode", "kld")
PhySig.df2$Spp = c("2 yrs", "2.5 yrs", "3 yrs", "4 yrs")
AE.m3 = ggplot(PhySig.df2, aes(x=Spp, y=Mean)) +
geom_point(size=2, pch=19, col = "red") +
geom_linerange(aes(ymin=Q025, ymax=Q975), colour="black") +
geom_hline(yintercept = 0,
linetype = "dotted",
colour = "red",
size = 1) +
theme_classic() +
xlab(NULL) +
ylab("Morphology PC-4") +
theme(axis.title.y = element_text(face="bold", size=14),
axis.text.y = element_text(face="bold", size=12),
axis.text.x = element_text(face="bold", size=14))
AE.m3
#grid.arrange(AE.m2, AE.m3, ncol=1)Map of PCA4
Pred.Grid$M3C = drop(Ap %*% Model.M3B$summary.random$field3$mean) +
Model.M3B$summary.fixed[1,1] +
Model.M3B$summary.fixed[2,1] +
mean(Model.M3B$summary.random$VcVx$mean) +
mean(Model.M3B$summary.random$Sex$mean) +
mean(Model.M3B$summary.random$Age$mean) +
mean(Model.M3B$summary.random$Clust$mean)
M3C.r = rasterize(Pred.Grid, Domain.r, "M3C",
background = NA)
mCols = brewer.pal(9, "BrBG")[-c(5,6,7)]
cr0 = (colorRampPalette(mCols)(n = 500))
cr = colorRampPalette(cr0, bias=2, space = "rgb")
rng = seq(-0.49, 0.722, 0.001)
names(M3C.r) = "PC 4"
M3.p = levelplot(M3C.r,
margin = FALSE,
xlab = " ", ylab = NULL,
col.regions = cr, at = rng,
colorkey = list(labels=list(cex=1.5),
space = "right"),
main=list("B", side=1, line=0.5),
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.5)) +
latticeExtra::layer(sp.polygons(SpineAndes, col = "black", lwd = 1))
M3.p
#grid.arrange(M2.p, M3.p, ncol=2)Formula for Joint-Model
MorphPC = list(theta=list(prior = "normal", param=c(0, 10)))
MorphPC2 = list(theta=list(prior = "normal", param=c(1, 70)))
MorphPC3 = list(theta=list(prior = "normal", param=c(0, 20)))
Frm.J = Y ~ -1 + intercept1 +
intercept2 +
intercept3 +
f(field2, model = spde1) +
f(field3, model = spde1) +
f(field1, model = spde1) +
f(field1x, copy = "field2",
fixed = FALSE) +
f(field2x, copy = "field3",
fixed = FALSE) +
f(field3x, copy = "field2",
fixed = FALSE) +
f(VcVx, model="generic0",
Cmatrix = VcV,
hyper = MorphPC3,
rankdef=1, constr=TRUE,
diagonal=1e-05) +
f(Sex, model="iid",
hyper = MorphPC) +
f(Age, model="iid",
hyper = MorphPC) +
f(Clust, model="iid",
hyper = MorphPC2) +
bio1 + bio15 + aPET + Pop + bio2 + bio3 + bio19Execute Model
#thetaJ = ModelJ$internal.summary.hyperpar$mean
thetaJ = c(-0.52024760, 0.11733760, 1.20819307, -4.15500234, 0.16643314, -2.95315759, 1, -0.59823897,
-2.50261798, 1.13007045, 0.21531003, 0.09058275, 0.11792903, 1.41772367, 0.69690391, 0.06779955)
ModelJ = inla(Frm.J,
data = inla.stack.data(J.Stack,
spde1=spde1,
VcV = VCV.mat2),
family = c("gaussian", "gaussian", "zeroinflatedbinomial1"),
verbose = TRUE,
control.family = list(error.prec, error.prec, error.prec),
control.fixed = list(prec = 0.1, prec.intercept = 0.001),
control.predictor = list(
A = inla.stack.A(J.Stack),
compute = TRUE,
link = 1),
control.mode = list(restart = TRUE, theta = thetaJ),
control.inla = list(int.strategy = "eb"),
control.results = list(return.marginals.random = TRUE,
return.marginals.predictor = TRUE),
control.compute=list(dic = TRUE, cpo = TRUE, waic = TRUE)) summary(ModelJ)##
## Call:
## c("inla(formula = Frm.J, family = c(\"gaussian\", \"gaussian\", \"zeroinflatedbinomial1\"), ", " data = inla.stack.data(J.Stack, spde1 = spde1, VcV = VCV.mat2), ", " verbose = TRUE, control.compute = list(dic = TRUE, cpo = TRUE, ", " waic = TRUE), control.predictor = list(A = inla.stack.A(J.Stack), ", " compute = TRUE, link = 1), control.family = list(error.prec, ", " error.prec, error.prec), control.inla = list(int.strategy = \"eb\"), ", " control.results = list(return.marginals.random = TRUE, return.marginals.predictor = TRUE), ", " control.fixed = list(prec = 0.1, prec.intercept = 0.001), ", " control.mode = list(restart = TRUE, theta = thetaJ))")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 4.0848 14132.2067 2.2861 14138.5777
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept1 -9.2893 0.4998 -10.3505 -9.2609 -8.3853 -9.2032 0
## intercept2 -0.3541 0.7161 -1.7602 -0.3541 1.0507 -0.3541 0
## intercept3 0.1544 0.7071 -1.2338 0.1544 1.5415 0.1544 0
## bio1 2.6863 0.8535 1.0219 2.6824 4.3708 2.6746 0
## bio15 1.8641 0.3451 1.2179 1.8528 2.5745 1.8301 0
## aPET -2.4293 0.9366 -4.2117 -2.4489 -0.5331 -2.4880 0
## Pop 0.0826 0.0365 0.0246 0.0776 0.1661 0.0657 0
## bio2 -0.0193 0.0519 -0.1212 -0.0193 0.0826 -0.0193 0
## bio3 0.2780 0.1369 0.0091 0.2780 0.5466 0.2780 0
## bio19 -1.1515 0.2356 -1.6139 -1.1515 -0.6893 -1.1516 0
##
## Random effects:
## Name Model
## field2 SPDE2 model
## field3 SPDE2 model
## field1 SPDE2 model
## VcVx Generic0 model
## Sex IID model
## Age IID model
## Clust IID model
## field1x Copy
## field2x Copy
## field3x Copy
##
## Model hyperparameters:
## mean sd
## Precision for the Gaussian observations 0.6015 0.0196
## Precision for the Gaussian observations[2] 1.1282 0.0379
## zero-probability parameter for zero-inflated binomial_1[3] 0.7634 0.0424
## Range for field2 0.0131 0.0033
## Stdev for field2 1.1231 0.1229
## Range for field3 0.6703 0.6147
## Stdev for field3 1.2940 0.7221
## Range for field1 0.0889 0.0290
## Stdev for field1 3.0757 0.4132
## Precision for VcVx 1.2676 0.2798
## Precision for Sex 1.1625 0.3744
## Precision for Age 1.1905 0.3825
## Precision for Clust 4.1408 0.3570
## Beta for field1x 1.3672 0.1936
## Beta for field2x 0.9723 0.2675
## Beta for field3x -0.4394 0.0746
## 0.025quant
## Precision for the Gaussian observations 0.5637
## Precision for the Gaussian observations[2] 1.0553
## zero-probability parameter for zero-inflated binomial_1[3] 0.7040
## Range for field2 0.0076
## Stdev for field2 0.9179
## Range for field3 0.1023
## Stdev for field3 0.3960
## Range for field1 0.0564
## Stdev for field1 2.5729
## Precision for VcVx 0.8100
## Precision for Sex 0.5938
## Precision for Age 0.6121
## Precision for Clust 3.4907
## Beta for field1x 1.0892
## Beta for field2x 0.4757
## Beta for field3x -0.5853
## 0.5quant
## Precision for the Gaussian observations 0.6012
## Precision for the Gaussian observations[2] 1.1277
## zero-probability parameter for zero-inflated binomial_1[3] 0.7572
## Range for field2 0.0128
## Stdev for field2 1.1093
## Range for field3 0.4944
## Stdev for field3 1.1327
## Range for field1 0.0811
## Stdev for field1 2.9790
## Precision for VcVx 1.2362
## Precision for Sex 1.1079
## Precision for Age 1.1334
## Precision for Clust 4.1223
## Beta for field1x 1.3348
## Beta for field2x 0.9600
## Beta for field3x -0.4396
## 0.975quant
## Precision for the Gaussian observations 0.6410
## Precision for the Gaussian observations[2] 1.2046
## zero-probability parameter for zero-inflated binomial_1[3] 0.8577
## Range for field2 0.0206
## Stdev for field2 1.4002
## Range for field3 2.3114
## Stdev for field3 3.1524
## Range for field1 0.1654
## Stdev for field1 4.1219
## Precision for VcVx 1.9051
## Precision for Sex 2.0475
## Precision for Age 2.1009
## Precision for Clust 4.8919
## Beta for field1x 1.8173
## Beta for field2x 1.5275
## Beta for field3x -0.2923
## mode
## Precision for the Gaussian observations 0.6008
## Precision for the Gaussian observations[2] 1.1266
## zero-probability parameter for zero-inflated binomial_1[3] 0.7159
## Range for field2 0.0122
## Stdev for field2 1.0744
## Range for field3 0.2650
## Stdev for field3 0.8634
## Range for field1 0.0647
## Stdev for field1 2.6840
## Precision for VcVx 1.1752
## Precision for Sex 1.0064
## Precision for Age 1.0276
## Precision for Clust 4.0820
## Beta for field1x 1.2036
## Beta for field2x 0.9155
## Beta for field3x -0.4401
##
## Expected number of effective parameters(std dev): 476.48(0.00)
## Number of equivalent replicates : 18.80
##
## Deviance Information Criterion (DIC) ...: 15917.15
## Effective number of parameters .........: 449.06
##
## Watanabe-Akaike information criterion (WAIC) ...: 15923.88
## Effective number of parameters .................: 409.64
##
## Marginal log-Likelihood: -8232.44
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computedCompare Density of the Marginal Terms
Red dotted line is from Model1.
Intj.df = as.data.frame(ModelJ$marginals.fix[[1]])
bio1j.df = as.data.frame(ModelJ$marginals.fix[[4]])
bio15j.df = as.data.frame(ModelJ$marginals.fix[[5]])
aPETj.df = as.data.frame(ModelJ$marginals.fix[[6]])
Popj.df = as.data.frame(ModelJ$marginals.fix[[7]])
CI2.df = as.data.frame(ModelJ$summary.fixed)[,c(3,5)]
In.plt = ggplot(Intj.df , aes(x, y)) +
geom_line(size = 1, col = "black") +
geom_line(data=Int.df,aes(x, y), size = 1,
col = "red", linetype = "dotted") +
geom_vline(xintercept = CI2.df[1,1], col = "gray") +
geom_vline(xintercept = CI2.df[1,2], col = "gray") +
geom_vline(xintercept = 0, col = "red") +
xlab(NULL) +
ylab("Density") +
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(size = 20, hjust = 0.5))
bio1.plt = ggplot(bio1j.df , aes(x, y)) +
geom_line(size = 1, col = "black") +
geom_line(data=bio1.df,aes(x, y), size = 1,
col = "red", linetype = "dotted") +
geom_vline(xintercept = CI2.df[4,1], col = "gray") +
geom_vline(xintercept = CI2.df[4,2], col = "gray") +
geom_vline(xintercept = 0, col = "red") +
xlab(NULL) +
ylab("Density") +
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(size = 20, hjust = 0.5))
bio15.plt = ggplot(bio15j.df , aes(x, y)) +
geom_line(size = 1, col = "black") +
geom_line(data=bio15.df,aes(x, y), size = 1,
col = "red", linetype = "dotted") +
geom_vline(xintercept = CI2.df[5,1], col = "gray") +
geom_vline(xintercept = CI2.df[5,2], col = "gray") +
geom_vline(xintercept = 0, col = "red") +
xlab(NULL) +
ylab(NULL) +
ggtitle("C") +
theme_classic() +
theme(axis.text=element_text(size=16),
axis.title.y = element_text(size = 20),
axis.title.x = element_text(size = 20, vjust=-2),
plot.title = element_text(size = 20, hjust = 0.5))
aPET.plt = ggplot(aPETj.df , aes(x, y)) +
geom_line(size = 1, col = "black") +
geom_line(data=aPET.df,aes(x, y), size = 1,
col = "red", linetype = "dotted") +
geom_vline(xintercept = CI2.df[6,1], col = "gray") +
geom_vline(xintercept = CI2.df[6,2], col = "gray") +
geom_vline(xintercept = 0, col = "red") +
xlab(NULL) +
ylab("Density") +
ggtitle("D") +
theme_classic() +
theme(axis.text=element_text(size=16),
axis.title.y = element_text(size = 20),
axis.title.x = element_text(size = 20, vjust=-2),
plot.title = element_text(size = 20, hjust = 0.5))
Pop.plt = ggplot(Popj.df , aes(x, y)) +
geom_line(size = 1, col = "black") +
geom_line(data=Pop.df,aes(x, y), size = 1,
col = "red", linetype = "dotted") +
geom_vline(xintercept = CI2.df[7,1], col = "gray") +
geom_vline(xintercept = CI2.df[7,2], col = "gray") +
geom_vline(xintercept = 0, col = "red") +
xlab(NULL) +
ylab(NULL) +
ggtitle("E") +
theme_classic() +
theme(axis.text=element_text(size=16),
axis.title.y = element_text(size = 20),
axis.title.x = element_text(size = 20, vjust=-2),
plot.title = element_text(size = 20, hjust = 0.5))
grid.arrange(In.plt, bio1.plt, bio15.plt, aPET.plt, Pop.plt, ncol = 2)
#Top.rw = grid.arrange(bio1.plt, bio15.plt, ncol = 2)
#Bot.row = grid.arrange(aPET.plt, Pop.plt, ncol = 2)
#grid.arrange(In.plt, Top.rw, Bot.row, ncol = 1)Sum of Linear Components
Pred.pnts = Pred.Grid
ModResult = ModelJ
ModMesh = mesh.dom
pLoc = cbind(Pred.pnts$Long, Pred.pnts$Lat)
pLoc = inla.mesh.map(pLoc,
projection = "longlat",
inverse = TRUE)
Ap = inla.spde.make.A(ModMesh, loc=pLoc) #project new points to mesh to obtain random field
Pred.pnts$RF = drop(Ap %*% ModResult$summary.random$field1$mean)
Pred.pnts$RF2 = drop(Ap %*% ModResult$summary.random$field2$mean)
Pred.pnts$RF3 = drop(Ap %*% ModResult$summary.random$field3$mean)
Pred.pnts$RF4 = drop(Ap %*% ModResult$summary.random$field1x$mean)
Pred.pnts$RF5 = drop(Ap %*% ModResult$summary.random$field2x$mean)
Pred.pnts$RF6 = drop(Ap %*% ModResult$summary.random$field3x$mean)
#Get the intercept
Pred.pnts$Int = ModResult$summary.fix[1,1]
Pred.pnts$Int2 = ModResult$summary.fix[2,1]
Pred.pnts$Int3 = ModResult$summary.fix[3,1]
#Additional effects
#Using the mean effect
Phylo.e = mean(ModResult$summary.random$VcVx$mean)
Sex.e = mean(ModResult$summary.random$Sex$mean)
Age.e = mean(ModResult$summary.random$Age$mean)
MM.e = mean(ModResult$summary.random$Clust$mean)
mydf = as.data.frame(ModResult$summary.fixed)
#Sum linear components
Pred.pnts@data = Pred.pnts@data %>%
mutate(LP = Int + RF + Int2 + RF2 + Int3 + RF3 + RF4 + RF5 + RF6 +
Phylo.e + Sex.e + Age.e + MM.e +
bio1*ModResult$summary.fix[4,1] +
bio15*ModResult$summary.fix[5,1] +
aPET*ModResult$summary.fix[6,1] +
Pop*ModResult$summary.fix[7,1] +
bio2*ModResult$summary.fix[8,1] +
bio3*ModResult$summary.fix[9,1] +
bio19*ModResult$summary.fix[10,1],
LP = ifelse(is.na(LP), 0, LP),
LPe = plogis(LP)) Compare Predicted Distribution
- Model1 B)Joint-model
Pred_osi1J = rasterize(Pred.pnts, Domain.r, "LPe",
background = NA)
Pred.stk = stack(Pred_osi1, Pred_osi1J)
names(Pred.stk) = c("A", "B")
rng = seq(0, 1, 0.01)
mCols = brewer.pal(9, "YlOrRd")[2:9]
cr0 = rev(colorRampPalette(rev(mCols))(n = 299))
cr = colorRampPalette(c("tan", cr0),
bias = 2, space = "rgb")
M1Jop = levelplot(Pred.stk,
margin = FALSE,
xlab = " ", ylab = NULL,
col.regions = cr, at = rng,
colorkey = list(labels=list(at=c(0, 0.25, 0.5, 0.75, 1),
labels=c("0%", "25%", "50%", "75%", "100%"), cex=1.5),
labels=list(cex=12),
space = "right"),
par.strip.text = list(fontface='bold', cex=1),
par.settings = list(axis.line = list(col = "black"),
strip.background = list(col = 'transparent'),
strip.border = list(col = 'transparent')),
scales = list(cex = 1.5)) +
latticeExtra::layer(sp.polygons(SpineAndes, col = "black", lwd = 1))
M1Jop
Compare Spatial Random Effect
- Model1 B)Joint-model
#Base Model
mRF.df = as.data.frame(Model1$summary.random$field1)[,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)
SpCor = 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(NULL) +
ylab("Spatial Field") +
scale_x_continuous(breaks=seq(0, 3400, 680), limits = c(0, 3400)) +
scale_y_continuous(breaks=seq(-15, 15, 5), limits = c(-15, 15)) +
theme_classic() +
ggtitle("A") +
theme(axis.text=element_text(size=14),
axis.title.y = element_text(size = 16),
axis.title.x = element_text(size = 16),
plot.title = element_text(hjust = 0.5, size = 16, face ="bold"))
SpCor
#Joint-model
mRF.df = as.data.frame(ModelJ$summary.random$field1)[,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("Spatial Index") +
ylab("Spatial Field") +
scale_x_continuous(breaks=seq(0, 3400, 680), limits = c(0, 3400)) +
scale_y_continuous(breaks=seq(-15, 15, 5), limits = c(-15, 15)) +
theme_classic() +
ggtitle("B") +
theme(axis.text=element_text(size=14),
axis.title.y = element_text(size = 16),
axis.title.x = element_text(size = 16),
plot.title = element_text(hjust = 0.5, size = 16, face ="bold"))
grid.arrange(SpCor, SpCor2, ncol=1)
Compare Random Field
Spatial residuals
RF_osi1J = rasterize(Pred.pnts, Domain.r, "RF",
background = NA)
RF1.stk = stack(GRF_osi0, RF_osi1, RF_osi1J)
names(RF1.stk) = c("A", "B", "C")
rng = seq(-12.1,12, 0.01)
mCols = c("#A50026", "#D73027", "#F46D43", "#FDAE61", "#E0F3F8", "#ABD9E9", "#74ADD1", "#4575B4","#313695")
#brewer.pal(11, "RdBu")
cr = colorRampPalette(c(rev(mCols)),
bias = 1, space = "rgb")
RFc = levelplot(RF1.stk,
margin = FALSE,
xlab = " ", ylab = NULL,
col.regions = cr, at = rng,
colorkey = list(labels=list(cex=1.5),
space = "right"),
par.strip.text = list(fontface='bold', cex=1),
par.settings = list(axis.line = list(col = "black"),
strip.background = list(col = 'transparent'),
strip.border = list(col = 'transparent')),
scales = list(cex = 1.5)) +
latticeExtra::layer(sp.polygons(SpineAndes, col = "black", lwd = 1))
RFc
Compare Deviance Information Criteria
#Model2 (non-spatial)
Model2$dic$dic## [1] 250.5033#Model0 (spatial effect only)
Model0$dic$dic## [1] 182.45#Model1 (spatial + environment)
Model1$dic$dic## [1] 198.7375#ModelJ (joint model)
tapply(ModelJ$dic$local.dic, ModelJ$dic$family, sum)[3]## 3
## 154.091Model Comparison
- MaxEnt (MaxEnt, Dismo)
- Random Forest (RandomForest)
- Gaussian random field species distribution models (GRaF package)
- Hierarchical Bayesian species distribution models (hSDM)
- Model1 (INLA/SPDE)
- Joint model (ModelJ, INLA/SPDE)
- Non-spatial model (Model2, INLA/SPDE)
MaxEnt
library(PresenceAbsence)
Pred.pnts = Mod.pnts
xFE.df = Pred.pnts@data #Copy point data to dataframe
#Presence locations
xFE.df$OBS2 = ifelse(FE.df$nRvTax == "amicus", 1, 0)
sum(xFE.df$OBS2)## [1] 24#Select environmental variables
select.var = c("bio1" , "bio15", "aPET", "Pop")
SDMvar.df = xFE.df[, select.var]
#call MaxEnt to fit data
MaXEnt_dar = maxent(x = SDMvar.df, p = xFE.df$OBS2)## Loading required namespace: rJava#Predict to all domain locations
PredGrid_ME_dar = Pred.Grid@data[, select.var]
Pred.Grid$Dar_ME = predict(MaXEnt_dar, PredGrid_ME_dar)
#Rasterize prediction
MaxEnt.ras = rasterize(Pred.Grid,
Domain.r, "Dar_ME",
background = NA)
#Predict at observations locations for model comparison
SDMvar.df$Dar_MEp = predict(MaXEnt_dar, SDMvar.df)
#Evalute performance
MaxEnt.df = as.data.frame(cbind(xFE.df$Idx, xFE.df$OBS2, SDMvar.df$Dar_MEp))
MaxEnt.out = presence.absence.accuracy(MaxEnt.df, threshold = 0.50)
#Calculate TSS
MaxEnt.out$TSS = MaxEnt.out$sensitivity + MaxEnt.out$specificity - 1
MaxEnt.out$MSE = mean((SDMvar.df$Dar_MEp - xFE.df$OBS2)^2)
MaxEnt.out$ModLab = "MaxEnt"Random Forest (RF)
suppressMessages(library(randomForest))
#fit model
rf.mod = randomForest(OBS2 ~ bio1 + bio15 + aPET + Pop, data=xFE.df, importance=T)
#print(rf.mod)
#round(importance(rf.mod), 2)
#Predict at observations locations for model comparison
xFE.df$RF.pred = as.numeric(predict(rf.mod, xFE.df))
#Evalute performance
RF.df = as.data.frame(cbind(xFE.df$Idx, xFE.df$OBS2, xFE.df$RF.pred))
RF.out = presence.absence.accuracy(RF.df, threshold = 0.50)
#Calculate TSS
RF.out$TSS = RF.out$sensitivity + RF.out$specificity - 1
RF.out$MSE = mean((xFE.df$RF.pred - xFE.df$OBS2)^2)
RF.out$ModLab = "RF"
#Predict to all domain locations
PredGrid_rf_dar = Pred.Grid@data[, select.var]
Pred.Grid$RF_pred = predict(rf.mod, PredGrid_rf_dar)
#Rasterize prediction
RF.ras = rasterize(Pred.Grid,
Domain.r, "RF_pred",
background = NA)GRaF
library(GRaF)
GRaF.cov = xFE.df[,select.var]
# fit a model
GRaF.mod = graf(xFE.df$OBS2, GRaF.cov)
#Predict at observations locations for model comparison
GRaF.fout = predict(GRaF.mod, GRaF.cov, CI = NULL)
GRaF.cov$graf.pred = as.numeric(GRaF.fout[,1])
#Evalute performance
graf.rdf = as.data.frame(cbind(xFE.df$Idx, xFE.df$OBS2, GRaF.cov$graf.pred))
GRaF.out = presence.absence.accuracy(graf.rdf, threshold = 0.50)
#Calculate TSS
GRaF.out$TSS = GRaF.out$sensitivity + GRaF.out$specificity - 1
GRaF.out$MSE = mean((GRaF.cov$graf.pred - xFE.df$OBS2)^2)
GRaF.out$ModLab = "GRaF"
#Predict to all domain locations
graf.map = Pred.Grid@data[, select.var]
Pred.grd2 = Pred.Grid
GRaF_pred = predict(GRaF.mod, graf.map, CI = NULL)
Pred.grd2$GRaF_pred = as.numeric(GRaF_pred[,1])
#Rasterize prediction
GRaF.ras = rasterize(Pred.grd2,
Domain.r, "GRaF_pred",
background = NA)hSDM
Hierarchical Bayesian species distribution models (hSDM Package)
suppressMessages(library(hSDM))
#Fit model
hSDM.mod = hSDM.binomial(presences = xFE.df$OBS2,
trials = rep(1, length(xFE.df$OBS2)),
suitability = ~bio1 + bio15 + aPET + Pop,
data = xFE.df,
burnin = 1000, mcmc = 1000, thin = 1,
beta.start = 0,
mubeta = 0, Vbeta = 1.0E6,
seed = 1976, verbose = 1, save.p = 0)
#summary(hSDM.mod$mcmc)
#Prediction for observation loactions
xFE.df$hSDM.pred = hSDM.mod$theta.pred
#Evalute performance
hSDM.df = as.data.frame(cbind(xFE.df$Idx, xFE.df$OBS2, xFE.df$hSDM.pred))
hSDM.out = presence.absence.accuracy(hSDM.df, threshold = 0.50)
#Calculate TSS
hSDM.out$TSS = hSDM.out$sensitivity + hSDM.out$specificity - 1
hSDM.out$MSE = mean((xFE.df$hSDM.pred - xFE.df$OBS2)^2)
hSDM.out$ModLab = "hSDM"
#Predict to all domain locations
data.pred = Pred.Grid@data
hSDM.dom = hSDM.binomial(presences = xFE.df$OBS2,
trials = rep(1, length(xFE.df$OBS2)),
suitability = ~bio1 + bio15 + aPET + Pop,
data = xFE.df,
suitability.pred = data.pred,
burnin = 1000, mcmc = 1000, thin = 1,
beta.start = 0,
mubeta = 0, Vbeta = 1.0E6,
seed = 1976, verbose = 1, save.p = 0)
Pred.Grid$hSDM.p = hSDM.dom$theta.pred
#Rasterize prediction
hSDM.ras = rasterize(Pred.Grid,
Domain.r, "hSDM.p",
background = NA)Undertake Observation location Predictions using Base Model (Model1)
ModResult = Model1
ModMesh = mesh.dom
pLoc = cbind(Pred.pnts$Long, Pred.pnts$Lat)
pLoc = inla.mesh.map(pLoc,
projection = "longlat",
inverse = TRUE)
Ap = inla.spde.make.A(ModMesh, loc=pLoc)
#Get the random field
Pred.pnts$RF = drop(Ap %*% ModResult$summary.random$field1$mean)
#Get the intercept
Pred.pnts$Int = ModResult$summary.fix[1,1]
#Get fixed effects
mydf = as.data.frame(ModResult$summary.fixed) #Write results to data frame
#Sum intercept, random effect, and fixed effects.
Pred.pnts@data = Pred.pnts@data %>%
mutate(LP = Int + RF +
bio1*ModResult$summary.fix[2,1] +
bio15*ModResult$summary.fix[3,1] +
aPET*ModResult$summary.fix[4,1] +
Pop*ModResult$summary.fix[5,1],
LP = ifelse(is.na(LP), 0, LP),
LPe = plogis(LP))
#Evalute performance
Mod1.df = as.data.frame(cbind(xFE.df$Idx, xFE.df$OBS2, Pred.pnts$LPe))
Mod1.out = presence.absence.accuracy(Mod1.df, threshold = 0.50)
#Calculate TSS
Mod1.out$TSS = Mod1.out$sensitivity + Mod1.out$specificity - 1
Mod1.out$MSE = mean((Pred.pnts$LPe - xFE.df$OBS2)^2)
Mod1.out$ModLab = "Mod1"Undertake Observation location Predictions using Joint-Model (ModelJ)
ModResult = ModelJB2
ModMesh = mesh.dom
pLoc = cbind(Pred.pnts$Long, Pred.pnts$Lat)
pLoc = inla.mesh.map(pLoc,
projection = "longlat",
inverse = TRUE)
Ap = inla.spde.make.A(ModMesh, loc=pLoc) #project new points to mesh to obtain random field
Pred.pnts$RF = drop(Ap %*% ModResult$summary.random$field1$mean)
Pred.pnts$RF2 = drop(Ap %*% ModResult$summary.random$field2$mean)
Pred.pnts$RF3 = drop(Ap %*% ModResult$summary.random$field3$mean)
Pred.pnts$RF4 = drop(Ap %*% ModResult$summary.random$field1x$mean)
Pred.pnts$RF5 = drop(Ap %*% ModResult$summary.random$field2x$mean)
Pred.pnts$RF6 = drop(Ap %*% ModResult$summary.random$field3x$mean)
#Get the intercept
Pred.pnts$Int = ModResult$summary.fix[1,1]
Pred.pnts$Int2 = ModResult$summary.fix[2,1]
Pred.pnts$Int3 = ModResult$summary.fix[3,1]
#Additional effects
#Using the mean effect
Phylo.e = mean(ModResult$summary.random$VcVx$mean)
Sex.e = mean(ModResult$summary.random$Sex$mean)
Age.e = mean(ModResult$summary.random$Age$mean)
MM.e = mean(ModResult$summary.random$Clust$mean)
mydf = as.data.frame(ModResult$summary.fixed)
#Sum linear components
Pred.pnts@data = Pred.pnts@data %>%
mutate(LP = Int + RF + Int2 + RF2 + Int3 +
RF3 + RF4 + RF5 + RF6 +
Phylo.e + Sex.e + Age.e + MM.e +
bio1*ModResult$summary.fix[4,1] +
bio15*ModResult$summary.fix[5,1] +
aPET*ModResult$summary.fix[6,1] +
Pop*ModResult$summary.fix[7,1] +
bio2*ModResult$summary.fix[8,1] +
bio3*ModResult$summary.fix[9,1] +
bio19*ModResult$summary.fix[10,1],
LP = ifelse(is.na(LP), 0, LP),
LPe = plogis(LP))
#Evaluate performance
JMod.df = as.data.frame(cbind(xFE.df$Idx, xFE.df$OBS2, Pred.pnts$LPe))
JMod.out = presence.absence.accuracy(JMod.df, threshold = 0.50)
#Calculate TSS
JMod.out$TSS = JMod.out$sensitivity + JMod.out$specificity - 1
JMod.out$MSE = mean((Pred.pnts$LPe - xFE.df$OBS2)^2)
JMod.out$ModLab = "JMod"Undertake Observation location Predictions using Non-spatial Model (Model2)
ModResult = Model2
#Get the intercept
Pred.pnts$Int = ModResult$summary.fix[1,1]
#Get fixed effects
mydf = as.data.frame(ModResult$summary.fixed)
Pred.pnts@data = Pred.pnts@data %>%
mutate(LP = Int +
bio1*ModResult$summary.fix[2,1] +
bio15*ModResult$summary.fix[3,1] +
aPET*ModResult$summary.fix[4,1] +
Pop*ModResult$summary.fix[5,1],
LP = ifelse(is.na(LP), 0, LP),
LPe = plogis(LP))
#Evalute performance
NSMod.df = as.data.frame(cbind(xFE.df$Idx, xFE.df$OBS2, Pred.pnts$LPe))
NSMod.out = presence.absence.accuracy(NSMod.df, threshold = 0.50)
#Calculate TSS
NSMod.out$TSS = NSMod.out$sensitivity + NSMod.out$specificity - 1
NSMod.out$MSE = mean((Pred.pnts$LPe - xFE.df$OBS2)^2)
NSMod.out$ModLab = "NSMod"Compare Model Metrics
ModResults = rbind(Mod1.out, JMod.out, NSMod.out, GRaF.out, MaxEnt.out, RF.out, hSDM.out)
ModResults.p = ModResults %>%
select(ModLab, sensitivity, specificity, AUC, TSS)
ModResults.p## ModLab sensitivity specificity AUC TSS
## 1 Mod1 0.2083333 0.9979392 0.9892560 0.2062725
## 2 JMod 0.9166667 0.9894384 0.9937962 0.9061051
## 3 NSMod 0.0000000 1.0000000 0.7892463 0.0000000
## 4 GRaF 0.0000000 1.0000000 0.8804686 0.0000000
## 5 MaxEnt 0.8750000 0.9289026 0.9263320 0.8039026
## 6 RF 0.3750000 0.9997424 0.9990608 0.3747424
## 7 hSDM 0.0000000 1.0000000 0.7850281 0.0000000#SelMod.tab2 = xtable(ModResults.p)
#print(SelMod.tab2, include.rownames = F)Compare Al Model Predictions
Mod.stk = stack(Pred_osi1, Pred_osi1J, NS_M2, GRaF.ras, MaxEnt.ras, RF.ras, hSDM.ras)
names(Mod.stk) = c("A", "B", "C", "D", "E", "F", "G")
#Random Forest produced negative predictions; setting these to zero for visualization
Mod.stk[Mod.stk<0] = 0
rng = seq(0, 1, 0.01)
mCols = brewer.pal(9, "YlOrRd")[2:9]
cr0 = rev(colorRampPalette(rev(mCols))(n = 299))
cr = colorRampPalette(c("tan", cr0),
bias = 1, space = "rgb")
AMods = levelplot(Mod.stk,
layout=c(7, 1),
margin = FALSE,
xlab = " ", ylab = NULL,
col.regions = cr, at = rng,
colorkey = list(labels=list(at=c(0, 0.25, 0.5, 0.75, 1),
labels=c("0.0", "0.25", "0.50", "0.75", "1.0"), cex=1.5),
labels=list(cex=12),
space = "bottom"),
par.strip.text = list(fontface='bold', cex=1),
par.settings = list(axis.line = list(col = "black"),
strip.background = list(col = 'transparent'),
strip.border = list(col = 'transparent')),
scales = list(cex = 1.5)) +
latticeExtra::layer(sp.polygons(SpineAndes, col = "black", lwd = 1))
AMods
Session Info
date()## [1] "Tue Oct 24 12:58:55 2017"sessionInfo()## R version 3.4.2 (2017-09-28)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 8.1 x64 (build 9600)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=English_United States.1252
## [2] LC_CTYPE=English_United States.1252
## [3] LC_MONETARY=English_United States.1252
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United States.1252
##
## attached base packages:
## [1] grid stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] ggtree_1.8.2 treeio_1.0.2 ape_4.1
## [4] bindrcpp_0.2 perturb_2.05 dplyr_0.7.4
## [7] dismo_1.1-4 PresenceAbsence_1.1.9 xtable_1.8-2
## [10] gridExtra_2.3 GISTools_0.7-4 MASS_7.3-47
## [13] ggthemes_3.4.0 spdep_0.6-15 mapproj_1.2-5
## [16] maptools_0.9-2 rgeos_0.3-25 rasterVis_0.41
## [19] latticeExtra_0.6-28 RColorBrewer_1.1-2 lattice_0.20-35
## [22] corrplot_0.77 factoextra_1.0.5 ggplot2_2.2.1
## [25] adegenet_2.1.0 ade4_1.7-8 FactoMineR_1.38
## [28] fields_9.0 maps_3.2.0 spam_2.1-1
## [31] dotCall64_0.9-04 raster_2.5-8 reshape2_1.4.2
## [34] rgdal_1.2-13 INLA_17.08.19 Matrix_1.2-11
## [37] sp_1.2-5 rgl_0.98.1 knitr_1.17
##
## loaded via a namespace (and not attached):
## [1] nlme_3.1-131 gmodels_2.16.2 rprojroot_1.2
## [4] tools_3.4.2 backports_1.1.1 R6_2.2.2
## [7] vegan_2.4-4 lazyeval_0.2.0 mgcv_1.8-20
## [10] colorspace_1.3-2 permute_0.9-4 compiler_3.4.2
## [13] flashClust_1.01-2 expm_0.999-2 labeling_0.3
## [16] GRaF_0.1-12 scales_0.5.0 hexbin_1.27.1
## [19] stringr_1.2.0 digest_0.6.12 foreign_0.8-69
## [22] rmarkdown_1.6 pkgconfig_2.0.1 htmltools_0.3.6
## [25] htmlwidgets_0.9 rlang_0.1.2 shiny_1.0.5
## [28] bindr_0.1 zoo_1.8-0 jsonlite_1.5
## [31] gtools_3.5.0 magrittr_1.5 leaps_3.0
## [34] Rcpp_0.12.13 munsell_0.4.3 scatterplot3d_0.3-40
## [37] stringi_1.1.5 plyr_1.8.4 parallel_3.4.2
## [40] gdata_2.18.0 ggrepel_0.7.0 deldir_0.1-14
## [43] splines_3.4.2 igraph_1.1.2 boot_1.3-20
## [46] markdown_0.8 seqinr_3.4-5 LearnBayes_2.15
## [49] glue_1.1.1 evaluate_0.10.1 httpuv_1.3.5
## [52] purrr_0.2.3 tidyr_0.7.1 gtable_0.2.0
## [55] assertthat_0.2.0 mime_0.5 coda_0.19-1
## [58] viridisLite_0.2.0 tibble_1.3.4 rJava_0.9-9
## [61] rvcheck_0.0.9 cluster_2.0.6