UP Marten (Martes americana)
Spatial-temporal
date() ## [1] "Tue Dec 11 23:19:44 2018"Options
library(knitr)
library(rgl)
opts_knit$set(verbose = TRUE)
knit_hooks$set(webgl = hook_webgl)Set Directory
setwd("~/Michigan_State/Talesha")Load Packages:
Get County Boundaries
Counties = readOGR(dsn = "C:/Users/humph173/Documents/Michigan_State/SLP_Beam_Diam/Counties_v17a",
layer = "Counties_v18a",
stringsAsFactors = FALSE)## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\humph173\Documents\Michigan_State\SLP_Beam_Diam\Counties_v17a", layer: "Counties_v18a"
## with 83 features
## It has 15 fields
## Integer64 fields read as strings: OBJECTID FIPSNUMUP = subset(Counties, PENINSULA == "upper")
UP.union = gUnaryUnion(UP)
#Rasterized version
Ras = raster(res = 0.01, ext = extent(UP.union),
crs = proj4string(UP.union))
Domain.r = rasterize(UP.union, Ras,
field = 0,
background = NA)
#Point grid version
Grd.pnt = rasterToPoints(Domain.r, spatial = TRUE)
Grd.pnt@data = Grd.pnt@data %>%
mutate(Long = Grd.pnt@coords[,1],
Lat = Grd.pnt@coords[,2]) %>%
select(-layer)Marten Harvest Data
Mart.df = read.csv("./Marten_cmb.csv",
header = TRUE, sep=",") %>%
mutate(Long = Y, #east & northing
Lat = X,
Year = Season_Yea,
Spp = "Marten",
OBS = 1,
Source = "MDNR") %>%
select(Spp, OBS, Year, Long, Lat, Source)
dim(Mart.df)## [1] 2382 6Mart.df %>%
group_by(Year) %>%
summarise(Count = length(Year))## # A tibble: 11 x 2
## Year Count
## <int> <int>
## 1 2006 181
## 2 2007 256
## 3 2008 244
## 4 2009 222
## 5 2010 233
## 6 2011 159
## 7 2012 251
## 8 2013 244
## 9 2014 247
## 10 2015 232
## 11 2016 113Convert to Points
Marten.pnt = SpatialPointsDataFrame(Mart.df[,c("Long", "Lat")], Mart.df)
proj4string(Marten.pnt) = proj4string(UP)
#Keep only those in the UP
Marten.pnt$UP = is.na(over(Marten.pnt, UP.union))
Marten.pnt = subset(Marten.pnt, UP == FALSE)Quick Plot
UP_map = fortify(UP)## Regions defined for each Polygonstmp.df = Marten.pnt@data %>%
select(Long, Lat, Year)
ggplot(UP_map, aes(long,lat, group=group)) +
geom_polygon(fill = "tan", col="black") +
geom_point(data=tmp.df,
aes(Long, Lat,
group=NULL, fill=NULL),
col = "red", size=1) +
xlab("Longitude") +
ylab("Latitude") +
ggtitle("Marten Harvest") +
guides(color=guide_legend(title = "Year")) +
facet_wrap(~Year, ncol = 3) +
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_blank(),
axis.title.x = element_text(size=16, face="bold"),
axis.title.y = element_text(size=16, face="bold"),
plot.title = element_text(size=20, face="bold", hjust = 0.5))
Change Projection
nProj = "+proj=utm +zone=16 +datum=NAD83 +units=km +no_defs +ellps=GRS80 +towgs84=0,0,0"
UP.prj = spTransform(UP, nProj)
UP.union.prj = spTransform(UP.union, nProj)
Marten.pnt.prj = spTransform(Marten.pnt, nProj)
Marten.pnt.prj@data = Marten.pnt.prj@data %>%
mutate(Longp = Marten.pnt.prj@coords[,1],
Latp = Marten.pnt.prj@coords[,2])Construct Mesh
max.edge = 2
bound.outer = 30
bdry = inla.sp2segment(UP.union.prj)
mesh = inla.mesh.2d(boundary = bdry,
loc=cbind(Marten.pnt.prj),
max.edge = c(1, 12)*max.edge,
cutoff = 5,
min.angle = 23,
offset = c(max.edge, bound.outer))
plot(mesh, lwd=0.5, main=" ",
rgl=F, draw.vertices = TRUE,
vertex.color = "blue")
mesh$n## [1] 1855Get Mesh Nodes
dd = as.data.frame(cbind(mesh$loc[,1],
mesh$loc[,2]))
names(dd) = c("Longp", "Latp")
dd = dd %>%
mutate(OBS = 0,
Source = "Mesh")
##Make a copy for each year
Year.lvls = levels(factor(Marten.pnt.prj$Year))
for(i in 1:length(Year.lvls)){
tmp.dd = dd %>%
mutate(Year = Year.lvls[i])
if(i==1){dd.df = tmp.dd}
else{dd.df = rbind(dd.df, tmp.dd)}}
dim(dd.df)## [1] 20405 5levels(factor(dd.df$Year))## [1] "2006" "2007" "2008" "2009" "2010" "2011" "2012" "2013" "2014" "2015"
## [11] "2016"head(dd.df)## Longp Latp OBS Source Year
## 1 522.7310 5050.779 0 Mesh 2006
## 2 700.1876 5066.286 0 Mesh 2006
## 3 489.1956 5057.737 0 Mesh 2006
## 4 539.7490 5065.947 0 Mesh 2006
## 5 512.3788 5064.130 0 Mesh 2006
## 6 526.6298 5064.886 0 Mesh 2006Combine Marten and Nodes
Mod.pnts = Marten.pnt.prj@data %>%
select(Longp, Latp, OBS, Source, Year)
Mod.pnts = rbind(Mod.pnts, dd.df)Construct Model
Select Time Steps
This model using an autoregressive term for temporal correlation.
k = length(levels(factor(Mod.pnts$Year)))
#Change years to integers
Mod.pnts$Step = as.integer(as.factor(Mod.pnts$Year))
range(Mod.pnts$Step)## [1] 1 11Spatial Indexing
Relating points to the mesh.
locs = cbind(Mod.pnts$Longp, Mod.pnts$Latp)
A.mart = inla.spde.make.A(mesh,
alpha = 2,
loc=locs,
group = Mod.pnts$Step)Set Spatial Prior
FE.df = Mod.pnts
spde0 = inla.spde2.pcmatern(mesh,
alpha = 2,
prior.range=c(500, 0.01),
prior.sigma=c(1, 0.01),
constr = TRUE)
field1 = inla.spde.make.index("field1",
spde0$n.spde,
n.group=k) Construct Data Stack
Organizing data for model fitting.
FE.lst = list(c(field1,
list(intercept1 = 1)),
list(Year = FE.df[,"Year"])) #Place holder for future covariates
Stack1 = inla.stack(data = list(PA = FE.df$OBS),
A = list(A.mart, 1),
effects = FE.lst,
tag = "mart.0")
#save(list=c("Stack1", "spde0"), file="./HPC/Mart_121118.RData") Setting priors & Constructing model formula
h.spec = list(theta=list(prior='pccor1', param=c(0, 0.9)))
YearPrior = list(prec = list(prior="pc.prec", param = c(1, 0.0001)))
Frm0 = PA ~ -1 + intercept1 +
f(field1,
model=spde,
group = field1.group,
control.group=list(model="ar1"),
hyper=h.spec) +
f(Year,
model = "rw1",
hyper = YearPrior) Run Model
#theta0 = Model0$internal.summary.hyperpar$mean
theta0 = c(5.019970, 1.272142, 7.275585, 2.745718)
Model0 = inla(Frm0,
data = inla.stack.data(Stack1, spde=spde0),
family = "binomial",
verbose = TRUE,
control.fixed = list(prec = 0.01,
prec.intercept = 0.001),
control.predictor = list(
A = inla.stack.A(Stack1),
compute = TRUE,
link = 1),
control.mode = list(restart = TRUE, theta = theta0),
control.inla = list(int.strategy = "eb"),
control.results = list(return.marginals.random = TRUE,
return.marginals.predictor = TRUE),
control.compute=list(dic = TRUE, cpo = TRUE, waic = TRUE)) summary(Model0)##
## Call:
## c("inla(formula = Frm0, family = \"binomial\", data = inla.stack.data(Stack1, ", " spde = spde0), verbose = TRUE, control.compute = list(dic = TRUE, ", " cpo = TRUE, waic = TRUE), control.predictor = list(A = inla.stack.A(Stack1), ", " compute = TRUE, link = 1), control.inla = list(int.strategy = \"eb\"), ", " control.results = list(return.marginals.random = TRUE, return.marginals.predictor = TRUE), ", " control.fixed = list(prec = 0.01, prec.intercept = 0.001), ", " control.mode = list(restart = TRUE, theta = theta0), num.threads = 8)" )
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 2.1450 6755.1342 4.0183 6761.2974
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept1 -6.3995 0.302 -7.0397 -6.3826 -5.8517 -6.3476 0
##
## Random effects:
## Name Model
## field1 SPDE2 model
## Year RW1 model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant 0.975quant
## Range for field1 153.0273 22.4951 114.0324 151.2096 202.3125
## Stdev for field1 3.6033 0.5052 2.7220 3.5641 4.7057
## GroupRho for field1 0.9984 0.0008 0.9965 0.9986 0.9995
## Precision for Year 17.7930 9.8353 5.6559 15.5661 42.8895
## mode
## Range for field1 147.5045
## Stdev for field1 3.4840
## GroupRho for field1 0.9989
## Precision for Year 11.9356
##
## Expected number of effective parameters(std dev): 340.25(0.00)
## Number of equivalent replicates : 66.89
##
## Deviance Information Criterion (DIC) ...............: 11359.31
## Deviance Information Criterion (DIC, saturated) ....: NULL
## Effective number of parameters .....................: 326.40
##
## Watanabe-Akaike information criterion (WAIC) ...: 11302.69
## Effective number of parameters .................: 261.20
##
## Marginal log-Likelihood: -5883.67
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computedYear Effect
mic.df = as.data.frame(Model0$summary.random$Year)
names(mic.df) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975")
ggplot(mic.df, aes(x=ID, y=Mean)) +
geom_point(size=2, pch=19, col = "red") +
geom_linerange(aes(ymin=Q025, ymax=Q975), colour="black") +
geom_hline(yintercept = 0,
linetype = "dotted",
colour = "red",
size = 1) +
theme_classic() +
xlab("Year Effect") +
ylab("Effect (logit)") +
theme(axis.title.y = element_text(face="bold", size=16),
axis.title.x = element_text(face="bold", size=14),
axis.text.y = element_text(face="bold", size=16),
axis.text.x = element_text(face="bold", size=10,
vjust=1, hjust=1, angle=45),
strip.text = element_text(face="bold", size = 20))
Spatial Fields
F.s1 = cbind(Model0$summary.random$field1$mean,
field1$field1.group)
s1xmean = list()
s1xmean = split(F.s1[,1], F.s1[,2])
Pred.pnts = spTransform(Grd.pnt, nProj)
Grd_A = inla.spde.make.A(mesh, loc=Pred.pnts@coords)
for(i in 1:k){
Pred.pnts$m1RE = drop(Grd_A %*% s1xmean[[i]])
tmp.r = rasterize(spTransform(Pred.pnts, proj4string(UP)),
Domain.r,
"m1RE",
background = NA)
if(i == 1){M0.stk = tmp.r}
else{M0.stk = stack(M0.stk, tmp.r)}
}
names(M0.stk) = Year.lvlsPlot Spatial Fields
rng = seq(-4.5, 7.1, 0.01)
mCols = c("#A50026", "#D73027", "#F46D43", "#FDAE61", "#E0F3F8", "#ABD9E9", "#74ADD1", "#4575B4","#313695")
#brewer.pal(11, "RdBu")
cr = colorRampPalette(c(rev(mCols)),
bias = 1.4, space = "rgb")
RF0osi = levelplot(M0.stk,
margin = FALSE,
layout=c(4, 4),
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(UP, col = "black", lwd = 1))
RF0osi
`