Phyllotis osilae
Bayesian Hierarchical Modeling
Joint-modeling of environmental, morphometric, and phylogenetic data.
John: jmh09r@my.fsu.edu
date() ## [1] "Sat Oct 07 16:29:44 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(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") # See "PhyllotisPrePro072117.Rmd" for source (data pre-processing)View Mesh
(click, drag, and zoom to view)
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 Model
PC Prior and spatial index.
FE.df = Mod.pnts@data
spde1 = inla.spde2.pcmatern(mesh.dom, alpha = 2,
prior.range=c(0.5, 0.9),
prior.sigma=c(1, 0.5),
constr = TRUE)
field1 = inla.spde.make.index("field1",
spde1$n.spde) #spatial indexError Precision
error.prec = list(hyper=list(theta=list(prior = "pc.prec", param=c(1, 0.01))))Construct Data Stack
Selecting P. osilae.
FE.df$OBS2 = ifelse(FE.df$nRvTax == "osilae" | FE.df$nRvTax == "phaeus" | FE.df$nRvTax == "tucumanus", 1, 0)
osi.pnt = subset(Mod.pnts, nRvTax == "osilae" | nRvTax == "phaeus" | nRvTax == "tucumanus")
sum(FE.df$OBS2)## [1] 74FE.lst = list(c(field1,
list(intercept1 = 1)),
list(bio1 = FE.df[,"bio1"],
bio2 = FE.df[,"bio2"],
bio3 = FE.df[,"bio3"],
bio4 = FE.df[,"bio4"],
bio5 = FE.df[,"bio5"],
bio6 = FE.df[,"bio6"],
bio7 = FE.df[,"bio7"],
bio8 = FE.df[,"bio8"],
bio9 = FE.df[,"bio9"],
bio10 = FE.df[,"bio10"],
bio11 = FE.df[,"bio11"],
bio12 = FE.df[,"bio12"],
bio13 = FE.df[,"bio13"],
bio14 = FE.df[,"bio14"],
bio15 = FE.df[,"bio15"],
bio16 = FE.df[,"bio16"],
bio17 = FE.df[,"bio17"],
bio18 = FE.df[,"bio18"],
bio19 = FE.df[,"bio19"],
DEM = FE.df[,"DEM"],
Slope = FE.df[,"Slope"],
Aspect = FE.df[,"Aspect"],
TWet = FE.df[,"TWet"],
TRI = FE.df[,"TRI"],
aPET = FE.df[,"aPET"],
AIT = FE.df[,"AIT"],
CMI = FE.df[,"CMI"],
Cont = FE.df[,"Cont"],
EmbQ = FE.df[,"EmbQ"],
GD0 = FE.df[,"GD0"],
GD5 = FE.df[,"GD5"],
mTc = FE.df[,"mTc"],
mTw = FE.df[,"mTw"],
mC10 = FE.df[,"mC10"],
PETcq = FE.df[,"PETcq"],
PETdq = FE.df[,"PETdq"],
sPET = FE.df[,"sPET"],
PETwq = FE.df[,"PETwq"],
PETwetq = FE.df[,"PETwetq"],
ThrmI = FE.df[,"ThrmI"],
Pop = FE.df[,"Pop"],
Nrd = FE.df[,"Nrd"]))
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)
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.0001),
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(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 = 1e-04))" )
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 1.4130 134.4151 0.7872 136.6152
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept1 -12.6341 2.013 -16.859 -12.4125 -9.4092 -12.3415 0
##
## Random effects:
## Name Model
## field1 SPDE2 model
##
## Model hyperparameters:
## mean sd
## zero-probability parameter for zero-inflated binomial_1 0.3722 0.0812
## Range for field1 0.0966 0.0234
## Stdev for field1 4.0029 0.6307
## 0.025quant
## zero-probability parameter for zero-inflated binomial_1 0.2171
## Range for field1 0.0601
## Stdev for field1 2.9124
## 0.5quant
## zero-probability parameter for zero-inflated binomial_1 0.3720
## Range for field1 0.0934
## Stdev for field1 3.9515
## 0.975quant mode
## zero-probability parameter for zero-inflated binomial_1 0.5323 0.3761
## Range for field1 0.1513 0.0872
## Stdev for field1 5.3882 3.8495
##
## Expected number of effective parameters(std dev): 61.28(7.479)
## Number of equivalent replicates : 63.74
##
## Deviance Information Criterion (DIC) ...: 297.87
## Effective number of parameters .........: 31.90
##
## Watanabe-Akaike information criterion (WAIC) ...: 284.72
## Effective number of parameters .................: 16.82
##
## Marginal log-Likelihood: -198.65
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computedPlot Gaussian Random Field
Pred.pnts = Pred.Grid #Copy Points
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) #project new points to mesh to obtain random field
#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(-0.57, 14.99, 0.01)
mCols = c("#A50026", "#D73027", "#F46D43", "#FDAE61", "#E0F3F8", "#ABD9E9", "#74ADD1", "#4575B4","#313695")
#brewer.pal(11, "RdBu")
cr = colorRampPalette(c(rev(mCols)),
bias = 4.5, 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(FE.df)[1]
n.covariates = 5
X = cbind(rep(1,n.data),
FE.df$DEM,
FE.df$TWet,
FE.df$bio3,
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.int = list(A = as.matrix(t(Q)%*%A.osi),
e = rep(0, n.covariates))) ## Note: method with signature 'matrix#sparseMatrix' chosen for function '%*%',
## target signature 'matrix#dgTMatrix'.
## "ANY#TsparseMatrix" would also be validBase Model (Model1)
Including all fixed and random effects.
Frm1 = OBS ~ -1 + intercept1 +
f(field1, model=spde) +
DEM + TWet + bio3 + Pop
#theta1 = Model1$internal.summary.hyperpar$mean
theta1 = c(-0.2907427, -2.4372514, 1.4966095) #previous run
#osilae
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), ", " control.mode = list(restart = TRUE, theta = theta1))")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 1.3721 220.7073 0.9128 222.9922
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept1 -11.5990 1.6484 -15.1349 -11.4474 -8.8384 -11.0867 0
## DEM 2.3069 0.4348 1.5307 2.2802 3.2336 2.2252 0
## TWet -1.0692 0.4616 -2.0236 -1.0530 -0.2075 -1.0211 0
## bio3 -0.8505 0.8146 -2.4892 -0.8402 0.7323 -0.8190 0
## Pop 0.0697 0.0914 -0.1403 0.0829 0.2139 0.1153 0
##
## Random effects:
## Name Model
## field1 SPDE2 model
##
## Model hyperparameters:
## mean sd
## zero-probability parameter for zero-inflated binomial_1 0.3587 0.0670
## Range for field1 0.0782 0.0232
## Stdev for field1 3.7568 0.7036
## 0.025quant
## zero-probability parameter for zero-inflated binomial_1 0.2292
## Range for field1 0.0448
## Stdev for field1 2.5756
## 0.5quant
## zero-probability parameter for zero-inflated binomial_1 0.359
## Range for field1 0.074
## Stdev for field1 3.687
## 0.975quant mode
## zero-probability parameter for zero-inflated binomial_1 0.4901 0.3639
## Range for field1 0.1348 0.0662
## Stdev for field1 5.3365 3.5483
##
## Expected number of effective parameters(std dev): 61.23(6.393)
## Number of equivalent replicates : 63.79
##
## Deviance Information Criterion (DIC) ...: 286.22
## Effective number of parameters .........: 34.26
##
## Watanabe-Akaike information criterion (WAIC) ...: 273.67
## Effective number of parameters .................: 19.15
##
## Marginal log-Likelihood: -194.06
## 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]])
DEM.df = as.data.frame(Model1$marginals.fix[[2]])
TWet.df = as.data.frame(Model1$marginals.fix[[3]])
bio3.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))
DEM.plt = ggplot(DEM.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))
TWet.plt = ggplot(TWet.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))
bio3.plt = ggplot(bio3.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, DEM.plt, TWet.plt, bio3.plt, Pop.plt, ncol = 2)
Sum of Linear Components
Probability of occurrence for Phyllotis osilae.
Pred.pnts = Pred.Grid #Copy Points
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 +
DEM*ModResult$summary.fix[2,1] +
TWet*ModResult$summary.fix[3,1] +
bio3*ModResult$summary.fix[4,1] +
Pop*ModResult$summary.fix[5,1],
LP = ifelse(is.na(LP), 0, LP),
LPe = plogis(LP)) Plot Probabilty of Occurence
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=c(-15,-10,-5, 0, 5, 10, 15, 20, 25),
labels=c("-15","-10","-5", "0", "5", "10", "15", "20", "25")) +
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(-4.55, 14.99, 0.01)
mCols = c("#A50026", "#D73027", "#F46D43", "#FDAE61", "#E0F3F8", "#ABD9E9", "#74ADD1", "#4575B4","#313695")
#brewer.pal(11, "RdBu")
cr = colorRampPalette(c(rev(mCols)),
bias = 2.2, 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 +
DEM + TWet + bio3 + Pop
#osilae
Model2 = inla(Frm2,
data = inla.stack.data(Stack1, spde=spde1),
family = "binomial",
verbose = TRUE,
control.fixed = list(prec = 0.1, prec.intercept = 0.0001),
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 = 1e-04))")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 1.9104 5.9342 0.7275 8.5721
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept1 -5.5417 0.2741 -6.1128 -5.5302 -5.0357 -5.5072 0
## DEM 1.3732 0.1547 1.0791 1.3698 1.6868 1.3631 0
## TWet -0.5786 0.2075 -0.9927 -0.5764 -0.1776 -0.5718 0
## bio3 -0.9061 0.2563 -1.4220 -0.9016 -0.4154 -0.8927 0
## Pop 0.0349 0.0295 -0.0325 0.0389 0.0820 0.0481 0
##
## The model has no random effects
##
## The model has no hyperparameters
##
## Expected number of effective parameters(std dev): 4.98(0.00)
## Number of equivalent replicates : 784.36
##
## Deviance Information Criterion (DIC) ...: 532.63
## Effective number of parameters .........: 4.986
##
## Watanabe-Akaike information criterion (WAIC) ...: 531.35
## Effective number of parameters .................: 3.554
##
## Marginal log-Likelihood: -278.99
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computedPlot Non-Spatial Model
Pred.pnts = Pred.Grid #Copy Points
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 +
DEM*ModResult$summary.fix[2,1] +
TWet*ModResult$summary.fix[3,1] +
bio3*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 = 2, 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
Combining Model1 above with two additional models to perform joint-modeling.
Additional Spatial Indices
Spatial indices are constructed for two additional models and two shared fields.
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)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")
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")
names(myCol) = levels(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"))
ggplot(wmap_df, aes(long,lat, group=group)) +
geom_polygon(fill = "tan", col="black") +
geom_point(data=Pnts.df,
aes(Long, Lat,
group=NULL, fill=NULL,
col = Pnts.df$Species), size=1.5) +
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_blank(),
panel.border = element_blank(),
legend.title = element_text(size=12, face="bold"),
axis.title.x = element_text(size=16, face="bold"),
axis.title.y = element_text(size=16, face="bold"),
plot.title = element_blank())
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 #copy2007 Variance Covariance Matrix
Loading and labeling matrix for model inclusion. Derived from Steppan et al, 2007.
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
Modeling PCA2 from DAPC Analysis. Including index for variance co-variance matrix (“VcVx”). The model response (“Y”) includes three elements; PCA2 and PCA3 from morphological analysis (all species) and presence/absence vector for the target species.
FE2.lst = list(c(field2,
list(intercept2 = 1)),
list(GD5 = FE2.df[,"GD5"],
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$Morph5),
A = list(Aphy, 1),
effects = FE2.lst,
tag = "phy0")
Phy.stk = inla.stack(data = list(Y = cbind(FE2.df$Morph5, NA, NA)), #three elements
A = list(Aphy, 1),
effects = FE2.lst,
tag = "phy0")Morphological Stack 2
Modeling PCA3 from DAPC Analysis. Including index for variance co-variance matrix (“VcVx”).
FE3.lst = list(c(field3,
list(intercept3 = 1)),
list(bio19 = FE2.df[,"bio19"],
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$Morph3),
A = list(Aphy2, 1),
effects = FE3.lst,
tag = "phy2")
Phy.stk2 = inla.stack(data = list(Y = cbind(NA, FE2.df$Morph3, NA)), #three elements
A = list(Aphy2, 1),
effects = FE3.lst,
tag = "phy2")Modify Data Stack from Model1
Model1 included a single vector response (presence/absence, 1/0). The stack is modified to create slots for the two morphological vectors and the shared spatial fields (field1x and field2x).
FE.df = Mod.pnts@data
FE.df$OBS2 = ifelse(FE.df$nRvTax == "osilae" | FE.df$nRvTax == "phaeus" | FE.df$nRvTax == "tucumanus", 1, 0)
FE.lst = list(c(field1, field1x, field2x,
list(intercept1 = 1)),
list(bio3 = FE.df[,"bio3"],
DEM = FE.df[,"DEM"],
TWet = FE.df[,"TWet"],
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")Update PC Prior
n.data = dim(FE.df)[1]
n.covariates = 5
X = cbind(rep(1,n.data),
FE.df$DEM,
FE.df$TWet,
FE.df$bio3,
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.int = list(A = as.matrix(t(Q)%*%A.osi),
e = rep(0, n.covariates))) Combine Three Data Stacks
J.Stack = inla.stack(Phy.stk, Phy.stk2, Stack1)Model for Morphological 5th PCA
(Spatial Effect Only)
Frm.M2A = M2 ~ -1 + intercept2 +
f(field2, model=spde1)
#theta2A = Model.M2A$internal.summary.hyperpar$mean
theta2A = c(-0.09491031, -3.83812689, 0.04044648)
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.3665 19.3932 0.5941 21.3538
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept2 0.1902 0.1254 -0.0143 0.1747 0.4764 0.1413 0
##
## Random effects:
## Name Model
## field2 SPDE2 model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant
## Precision for the Gaussian observations 3.0347 0.0896 2.8617 3.0336
## Range for field2 0.0852 0.0233 0.0502 0.0814
## Stdev for field2 0.8607 0.1406 0.6281 0.8451
## 0.975quant mode
## Precision for the Gaussian observations 3.2143 3.0319
## Range for field2 0.1411 0.0742
## Stdev for field2 1.1790 0.8114
##
## Expected number of effective parameters(std dev): 115.07(7.334)
## Number of equivalent replicates : 21.96
##
## Deviance Information Criterion (DIC) ...: 4487.81
## Effective number of parameters .........: 117.75
##
## Watanabe-Akaike information criterion (WAIC) ...: 4481.95
## Effective number of parameters .................: 103.21
##
## Marginal log-Likelihood: -2342.62
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computedAdding covariates
hc1 = list(theta=list(prior = "normal", param=c(0, 10))) #prior for "iid"
Frm.M2B = M2 ~ -1 + intercept2 +
f(field2, model=spde1) +
f(VcVx, model="generic0",
Cmatrix = VcV,
hyper = hc1,
rankdef=1, constr=TRUE,
diagonal=1e-05) +
f(Sex, model="iid",
hyper = hc1) +
f(Age, model="iid",
hyper = hc1) +
f(Clust, model="iid",
hyper = hc1) +
GD5
#theta2B = Model.M2B$internal.summary.hyperpar$mean
theta2B = c(1.42171300, -1.04822645, -0.46586519, 0.51644606, 0.09567485, 0.13234126, 2.14607270)
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), ", " control.mode = list(restart = TRUE, theta = theta2B))")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 1.4770 39.9903 1.0003 42.4676
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept2 0.4036 0.7597 -1.0939 0.4031 1.9014 0.4020 0
## GD5 0.0910 0.0410 0.0106 0.0909 0.1717 0.0907 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 4.1452 0.1283 3.8975 4.1436
## Range for field2 0.4573 0.3448 0.1055 0.3631
## Stdev for field2 0.7306 0.3515 0.2770 0.6552
## Precision for VcVx 1.7665 0.5657 0.9122 1.6817
## Precision for Sex 1.1606 0.3767 0.5939 1.1034
## Precision for Age 1.2038 0.3887 0.6183 1.1451
## Precision for Clust 8.6161 1.0137 6.8174 8.5486
## 0.975quant mode
## Precision for the Gaussian observations 4.403 4.1410
## Range for field2 1.368 0.2381
## Stdev for field2 1.621 0.5318
## Precision for VcVx 3.112 1.5249
## Precision for Sex 2.058 0.9980
## Precision for Age 2.130 1.0367
## Precision for Clust 10.788 8.4074
##
## Expected number of effective parameters(std dev): 315.41(12.41)
## Number of equivalent replicates : 8.012
##
## Deviance Information Criterion (DIC) ...: 3899.04
## Effective number of parameters .........: 317.81
##
## Watanabe-Akaike information criterion (WAIC) ...: 3902.27
## Effective number of parameters .................: 279.96
##
## Marginal log-Likelihood: -2121.62
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computedCompare Morph5 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(-1.6, 0.99, 0.001)
mCols = c("#A50026", "#D73027", "#F46D43", "#FDAE61",
"#E0F3F8", "#ABD9E9", "#74ADD1", "#4575B4","#313695")
cr = colorRampPalette(c(rev(mCols)),
bias=0.7,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")
PhySig.df2$Spp = CNames
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 PCA-2") +
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=1, angle=90))
Repeated Measures Effect (Morph5)
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)
ggplot(PhySig.df2, aes(x=ID, y=Mean)) +
geom_point(size=2, pch=1, col = "darkgray") +
geom_linerange(aes(ymin=Q025, ymax=Q975), colour="gray") +
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 PCA-2") +
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())
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")
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 PCA-5") +
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=1, angle=90))
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")
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 PCA-5") +
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=1, angle=90))
Map of PCA5
Pred.Grid$M2C = drop(Ap %*% Model.M2B$summary.random$field2$mean) +
Model.M2B$summary.fixed[1,1] +
Pred.Grid$GD5*Model.M2B$summary.fixed[2,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=2.5, space = "rgb")
rng = seq(-0.5, 0.98, 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"),
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 3rd PCA
(Spatial Effect Only)
Frm.M3A = M3 ~ -1 + intercept3 +
f(field3, model=spde1)
#theta3A = Model.M3A$internal.summary.hyperpar$mean
theta3A = c(0.016777155, -2.808926794, 0.005255357)
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
## 2.0569 14.9529 0.5730 17.5827
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept3 0.7805 0.2658 0.3386 0.751 1.3813 0.6868 0
##
## Random effects:
## Name Model
## field3 SPDE2 model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant
## Precision for the Gaussian observations 0.5441 0.0159 0.5132 0.5439
## Range for field3 0.0873 0.0242 0.0500 0.0838
## Stdev for field3 1.8936 0.3136 1.3603 1.8643
## 0.975quant mode
## Precision for the Gaussian observations 0.5760 0.5436
## Range for field3 0.1443 0.0773
## Stdev for field3 2.5904 1.8045
##
## Expected number of effective parameters(std dev): 109.52(7.29)
## Number of equivalent replicates : 23.07
##
## Deviance Information Criterion (DIC) ...: 8828.16
## Effective number of parameters .........: 111.91
##
## Watanabe-Akaike information criterion (WAIC) ...: 8822.90
## Effective number of parameters .................: 98.60
##
## Marginal log-Likelihood: -4507.90
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computedModel for Morphological 3rd PCA
Covariates added.
Frm.M3B = M3 ~ -1 + intercept3 +
f(field3, model=spde1) +
f(VcVx, model="generic0",
Cmatrix = VcV,
hyper = hc1,
rankdef=1, constr=TRUE,
diagonal=1e-05) +
f(Sex, model="iid",
hyper = hc1) +
f(Age, model="iid",
hyper = hc1) +
f(Clust, model="iid",
hyper = hc1) +
bio19
#theta3B = Model.M3B$internal.summary.hyperpar$mean
theta3B = c(-0.22616501, -2.70967238, -0.51123069, 0.09317573, 0.10334072, 0.09823951, 0.79876206)
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.mode = list(restart = TRUE, theta = theta3B),
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.4356 110.8929 0.8941 113.2227
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept3 0.006 0.7978 -1.5708 0.0071 1.5744 0.0093 0
## bio19 -0.723 0.2585 -1.2370 -0.7219 -0.2152 -0.7199 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 0.7980 0.0249 0.7502 0.7976
## Range for field3 0.0799 0.0580 0.0246 0.0630
## Stdev for field3 0.6141 0.1375 0.4038 0.5931
## Precision for VcVx 1.1391 0.3177 0.6485 1.0939
## Precision for Sex 1.1658 0.3735 0.5827 1.1176
## Precision for Age 1.1578 0.3660 0.5916 1.1079
## Precision for Clust 2.2425 0.3001 1.7181 2.2194
## 0.975quant mode
## Precision for the Gaussian observations 0.8481 0.7968
## Range for field3 0.2334 0.0434
## Stdev for field3 0.9387 0.5502
## Precision for VcVx 1.8860 1.0093
## Precision for Sex 2.0334 1.0248
## Precision for Age 2.0159 1.0139
## Precision for Clust 2.8946 2.1711
##
## Expected number of effective parameters(std dev): 290.78(13.55)
## Number of equivalent replicates : 8.69
##
## Deviance Information Criterion (DIC) ...: 8042.41
## Effective number of parameters .........: 293.49
##
## Watanabe-Akaike information criterion (WAIC) ...: 8055.76
## Effective number of parameters .................: 269.55
##
## Marginal log-Likelihood: -4159.23
## 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$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(-3.56, 2.49, 0.001)
mCols = c("#A50026", "#D73027", "#F46D43", "#FDAE61",
"#E0F3F8", "#ABD9E9", "#74ADD1", "#4575B4","#313695")
cr = colorRampPalette(c(rev(mCols)),
bias=0.75,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 Morph3
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
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 PCA-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=1, angle=90))
Repeated Measures Effect (Morph3)
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)
ggplot(PhySig.df2, aes(x=ID, y=Mean)) +
geom_point(size=2, pch=1, col = "darkgray") +
geom_linerange(aes(ymin=Q025, ymax=Q975), colour="gray") +
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 PCA-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())
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")
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 PCA-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=1, angle=90))
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")
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 PCA-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=1, angle=90))
Map of PCA3
Pred.Grid$M3C = drop(Ap %*% Model.M3B$summary.random$field3$mean) +
Model.M3B$summary.fixed[1,1] +
Pred.Grid$bio19*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 = rev(colorRampPalette(mCols)(n = 500))
cr = colorRampPalette(cr0, bias=0.22, space = "rgb")
rng = seq(-6.38, 1.52, 0.001)
M2.p = levelplot(M3C.r,
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.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
Formula for Joint-Model
hc3 = hc2 = list(theta=list(prior = "normal", param=c(0, 3))) #prior for shared spatial fields
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, hyper=hc2) +
f(field2x, copy = "field3",
fixed = FALSE, hyper=hc3) +
f(VcVx, model="generic0",
Cmatrix = VcV,
hyper =hc1,
rankdef=1, constr=TRUE,
diagonal=1e-05) +
f(Clust, model="iid",
hyper = hc1) +
DEM + TWet + bio3 + Pop + GD5 + bio19Execute Modelkeep4x
#thetaJ = ModelJ$internal.summary.hyperpar$mean
thetaJ = c(1.33673058, -0.59656711, -0.59715131, -3.90055312, -1.12287332, -2.78505805, 0.34285537,
-2.58614577, 1.27763381, 4.17915393, 2.98391653, 0.06797478, 0.61024328)
ModelJ2 = 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(ModelJ2)##
## 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
## 3.4576 1173.1349 1.8314 1178.4240
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept1 -11.2881 0.7030 -12.7813 -11.2479 -10.0174 -11.1662 0
## intercept2 0.2227 0.0522 0.1202 0.2227 0.3252 0.2227 0
## intercept3 -0.2020 0.2494 -0.6915 -0.2020 0.2873 -0.2020 0
## DEM 2.3749 0.3717 1.6781 2.3630 3.1394 2.3393 0
## TWet -1.0689 0.4449 -1.9603 -1.0626 -0.2127 -1.0503 0
## bio3 -1.0301 0.8077 -2.6587 -1.0151 0.5140 -0.9853 0
## Pop 0.0698 0.0829 -0.1233 0.0834 0.1943 0.1187 0
## GD5 0.0774 0.0401 -0.0014 0.0774 0.1562 0.0774 0
## bio19 -1.2418 0.3146 -1.8594 -1.2419 -0.6247 -1.2419 0
##
## Random effects:
## Name Model
## field2 SPDE2 model
## field3 SPDE2 model
## field1 SPDE2 model
## VcVx Generic0 model
## Clust IID model
## field1x Copy
## field2x Copy
##
## Model hyperparameters:
## mean sd
## Precision for the Gaussian observations 3.8733 0.1226
## Precision for the Gaussian observations[2] 0.5510 0.0168
## zero-probability parameter for zero-inflated binomial_1[3] 0.3594 0.0672
## Range for field2 0.0230 0.0085
## Stdev for field2 0.3056 0.0392
## Range for field3 0.0625 0.0244
## Stdev for field3 1.3995 0.2261
## Range for field1 0.0783 0.0248
## Stdev for field1 3.5766 0.6792
## Precision for VcVx 1.8062 0.5761
## Precision for Clust 11.2440 1.4355
## Beta for field1x 0.0471 0.5351
## Beta for field2x 0.6149 0.3476
## 0.025quant
## Precision for the Gaussian observations 3.6366
## Precision for the Gaussian observations[2] 0.5187
## zero-probability parameter for zero-inflated binomial_1[3] 0.2306
## Range for field2 0.0107
## Stdev for field2 0.2362
## Range for field3 0.0291
## Stdev for field3 1.0183
## Range for field1 0.0430
## Stdev for field1 2.4314
## Precision for VcVx 0.9143
## Precision for Clust 8.6688
## Beta for field1x -0.9987
## Beta for field2x -0.0635
## 0.5quant
## Precision for the Gaussian observations 3.8718
## Precision for the Gaussian observations[2] 0.5508
## zero-probability parameter for zero-inflated binomial_1[3] 0.3591
## Range for field2 0.0216
## Stdev for field2 0.3029
## Range for field3 0.0578
## Stdev for field3 1.3768
## Range for field1 0.0737
## Stdev for field1 3.5118
## Precision for VcVx 1.7287
## Precision for Clust 11.1597
## Beta for field1x 0.0456
## Beta for field2x 0.6133
## 0.975quant
## Precision for the Gaussian observations 4.1192
## Precision for the Gaussian observations[2] 0.5847
## zero-probability parameter for zero-inflated binomial_1[3] 0.4922
## Range for field2 0.0436
## Stdev for field2 0.3901
## Range for field3 0.1232
## Stdev for field3 1.9065
## Range for field1 0.1390
## Stdev for field1 5.0957
## Precision for VcVx 3.1527
## Precision for Clust 14.3128
## Beta for field1x 1.1034
## Beta for field2x 1.3022
## mode
## Precision for the Gaussian observations 3.8695
## Precision for the Gaussian observations[2] 0.5505
## zero-probability parameter for zero-inflated binomial_1[3] 0.3620
## Range for field2 0.0190
## Stdev for field2 0.2976
## Range for field3 0.0497
## Stdev for field3 1.3289
## Range for field1 0.0654
## Stdev for field1 3.3852
## Precision for VcVx 1.5817
## Precision for Clust 11.0002
## Beta for field1x 0.0404
## Beta for field2x 0.6077
##
## Expected number of effective parameters(std dev): 477.04(0.00)
## Number of equivalent replicates : 18.78
##
## Deviance Information Criterion (DIC) ...: 13133.35
## Effective number of parameters .........: 449.42
##
## Watanabe-Akaike information criterion (WAIC) ...: 13127.39
## Effective number of parameters .................: 396.03
##
## Marginal log-Likelihood: -6870.00
## 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.
A) Intercept
B) Elevation (DEM)
C) Growing Days (5-degree threshold) (GD5)
D) Topographic wetness index (TWet)
E) Isothermality (BIO2/BIO7)*100 (Bio3)
F) Human population density (Pop)
Intj.df = as.data.frame(ModelJ$marginals.fix[[1]])
DEMj.df = as.data.frame(ModelJ$marginals.fix[[4]])
TWetj.df = as.data.frame(ModelJ$marginals.fix[[5]])
bio3j.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=Intj.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))
DEM.plt = ggplot(DEMj.df , aes(x, y)) +
geom_line(size = 1, col = "black") +
geom_line(data=DEM.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(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))
TWet.plt = ggplot(TWetj.df , aes(x, y)) +
geom_line(size = 1, col = "black") +
geom_line(data=TWet.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("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))
bio3.plt = ggplot(bio3j.df , aes(x, y)) +
geom_line(size = 1, col = "black") +
geom_line(data=bio3.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(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(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("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, DEM.plt, TWet.plt, bio3.plt, Pop.plt, ncol = 2)
Sum of Linear Components
(Joint Model)
Pred.pnts = Pred.Grid #Copy Points
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)
#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]
mydf = as.data.frame(ModResult$summary.fixed)
Pred.pnts@data = Pred.pnts@data %>%
mutate(LP = Int + RF + Int2 + RF2 + Int3 + RF3 + RF4 + RF5 +
DEM*ModResult$summary.fix[4,1] +
TWet*ModResult$summary.fix[5,1] +
bio3*ModResult$summary.fix[6,1] +
Pop*ModResult$summary.fix[7,1] +
GD5*ModResult$summary.fix[8,1] +
bio19*ModResult$summary.fix[9,1],
LP = ifelse(is.na(LP), 0, LP),
LPe = plogis(LP))
RF_osi1J = rasterize(Pred.pnts, Domain.r, "RF",
background = NA)Compare Predicted Distributions
- 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(-20, 20, 5), limits = c(-20, 20)) +
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"))
#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(-20, 20, 5), limits = c(-20, 20)) +
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 Effect Only (Model0)
- Spatial Effect and Covariates (Model1)
- Joint Model (ModelJ)
RF1.stk = stack(GRF_osi0, RF_osi1, RF_osi1J)
names(RF1.stk) = c("A", "B", "C")
rng = seq(-4.7, 14.99, 0.01)
mCols = c("#A50026", "#D73027", "#F46D43", "#FDAE61", "#E0F3F8", "#ABD9E9", "#74ADD1", "#4575B4","#313695")
#brewer.pal(11, "RdBu")
cr = colorRampPalette(c(rev(mCols)),
bias = 2, 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
(Joint model has slightly greater parsimony)
#Model1
Model1$dic$dic## [1] 286.2182#ModelJ
tapply(ModelJ$dic$local.dic, ModelJ$dic$family, sum)[3]## 3
## 282.0807Model Comparison
- MaxEnt (MaxEnt, Dismo)
- Random Forest (RandomForest)
- Gaussian random field species distribution models (GRaF package)
- Hierarchical Bayesian species distribution models (hSDM)
- Model1
- Joint model (ModelJ)
MaxEnt
Pred.pnts = Mod.pnts
xFE.df = Pred.pnts@data #Copy point data to dataframe
#Presence locations
xFE.df$OBS2 = ifelse(xFE.df$nRvTax == "osilae" | xFE.df$nRvTax == "phaeus" | xFE.df$nRvTax == "tucumanus", 1, 0)
sum(xFE.df$OBS2)
#Select environmental variables
select.var = c("DEM" , "TWet", "bio3", "Pop")
SDMvar.df = xFE.df[, select.var]
#call MaxEnt to fit data
MaXEnt_dar = maxent(x = SDMvar.df, p = xFE.df$OBS2)
#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 ~ DEM + TWet + bio3 + Pop, data=xFE.df)
rf.mod
#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
Gaussian random field species distribution models (GRaF package)
suppressMessages(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 = ~DEM + TWet + bio3 + 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 = ~DEM + TWet + bio3 + 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 +
DEM*ModResult$summary.fix[2,1] +
TWet*ModResult$summary.fix[3,1] +
bio3*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 = 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)
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)
#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]
mydf = as.data.frame(ModResult$summary.fixed)
Pred.pnts@data = Pred.pnts@data %>%
mutate(LP = Int + RF + Int2 + RF2 + Int3 + RF3 + RF4 + RF5 +
DEM*ModResult$summary.fix[4,1] +
TWet*ModResult$summary.fix[5,1] +
bio3*ModResult$summary.fix[6,1] +
Pop*ModResult$summary.fix[7,1] +
GD5*ModResult$summary.fix[8,1] +
bio19*ModResult$summary.fix[9,1],
LP = ifelse(is.na(LP), 0, LP),
LPe = plogis(LP))
#Evalute 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 +
DEM*ModResult$summary.fix[2,1] +
TWet*ModResult$summary.fix[3,1] +
bio3*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
TSS = True Skill Statistic
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)
colnames(ModResults.p) = c("Model", "Sensitivity", "Specificity", "AUC", "TSS")
ModResults.p## Model Sensitivity Specificity AUC TSS
## 1 Mod1 0.7567568 0.9919102 0.9934213 0.748666986
## 2 JMod 0.8783784 0.9861691 0.9934143 0.864547481
## 3 NSMod 0.0000000 0.9992171 0.9296148 -0.000782881
## 4 GRaF 0.0000000 1.0000000 0.9539017 0.000000000
## 5 MaxEnt 0.7972973 0.9618998 0.9767393 0.759197089
## 6 RF 0.7162162 0.9994781 0.9993758 0.715694296
## 7 hSDM 0.0000000 0.9989562 0.9288178 -0.001043841Compare All Model Predictions
- Model1
- Joint Model (ModelJ)
- Non-Spatial Model (Model2)
- GRaF Model
- MaxEnt
- Random Forrest
- hSDM
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 = 2, 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%", "25%", "50%", "75%", "100%"), 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
#save.image("D:/Phyllotis/PhyProj/PhyloProj/Osi_091617.RData")
date()## [1] "Sat Oct 07 16:31:27 2017"sessionInfo()## R version 3.4.1 (2017-06-30)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 8.1 x64 (build 9600)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=English_United States.1252
## [2] LC_CTYPE=English_United States.1252
## [3] LC_MONETARY=English_United States.1252
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United States.1252
##
## attached base packages:
## [1] grid stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] bindrcpp_0.2 dplyr_0.7.2 dismo_1.1-4
## [4] PresenceAbsence_1.1.9 gridExtra_2.2.1 GISTools_0.7-4
## [7] MASS_7.3-47 ggthemes_3.4.0 spdep_0.6-13
## [10] mapproj_1.2-5 maptools_0.9-2 rgeos_0.3-23
## [13] rasterVis_0.41 latticeExtra_0.6-28 RColorBrewer_1.1-2
## [16] lattice_0.20-35 corrplot_0.77 factoextra_1.0.5
## [19] ggplot2_2.2.1 adegenet_2.0.1 ade4_1.7-8
## [22] FactoMineR_1.36 fields_9.0 maps_3.2.0
## [25] spam_2.1-1 dotCall64_0.9-04 raster_2.5-8
## [28] reshape2_1.4.2 rgdal_1.2-8 INLA_17.08.19
## [31] Matrix_1.2-10 sp_1.2-5 rgl_0.98.1
## [34] 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.1 backports_1.1.0 R6_2.2.2
## [7] vegan_2.4-3 lazyeval_0.2.0 mgcv_1.8-17
## [10] colorspace_1.3-2 permute_0.9-4 compiler_3.4.1
## [13] flashClust_1.01-2 expm_0.999-2 labeling_0.3
## [16] GRaF_0.1-12 scales_0.4.1 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.4
## [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.12 munsell_0.4.3 ape_4.1
## [37] scatterplot3d_0.3-40 stringi_1.1.5 plyr_1.8.4
## [40] parallel_3.4.1 gdata_2.18.0 ggrepel_0.6.5
## [43] deldir_0.1-14 splines_3.4.1 igraph_1.1.2
## [46] boot_1.3-19 markdown_0.8 seqinr_3.4-5
## [49] LearnBayes_2.15 glue_1.1.1 evaluate_0.10.1
## [52] httpuv_1.3.5 gtable_0.2.0 assertthat_0.2.0
## [55] mime_0.5 xtable_1.8-2 coda_0.19-1
## [58] viridisLite_0.2.0 tibble_1.3.4 cluster_2.0.6