SVC Hare Occupancy
1 Data Prep
1.1 Load and Prep Observation Data
Michigan Loading raw Michigan data, dropping previously acquired environmental covariates and summing trial successes.
MI.df = read.csv("./Data_102418/mi_data10-24-18.csv",
header = TRUE, stringsAsFactors = FALSE, sep=",")
MI.df = MI.df %>%
mutate(Long = Lat, #these appear reversed in original data
Lat = Lon,
State = "Michigan", #transect state
Spp = "Hare",
OBS = Hare, #overall presence (at any trial)
Counts = T.1 + T.2 + T.3 + T.4 + T.5 + T.6 + T.7 + T.8 + T.9, #Sum across trials
Trials = 9, #No NAs in this dataset all have 9 trials
Data = "Fit",
Site2 = paste("M", 1:nrow(MI.df), sep=".")) %>%
select(Site2, Long, Lat, State, OBS, Counts, Trials, Data, Spp)
dim(MI.df) #Check data## [1] 125 9head(MI.df)## Site2 Long Lat State OBS Counts Trials Data Spp
## 1 M.1 -88.24388 47.41681 Michigan 1 6 9 Fit Hare
## 2 M.2 -84.58309 45.53829 Michigan 1 6 9 Fit Hare
## 3 M.3 -84.14892 45.46854 Michigan 1 6 9 Fit Hare
## 4 M.4 -90.07851 46.35078 Michigan 1 1 9 Fit Hare
## 5 M.5 -86.56820 46.37911 Michigan 1 1 9 Fit Hare
## 6 M.6 -86.75846 46.25573 Michigan 1 9 9 Fit HareWisconsin Loading data as above, but trials vary by transect, therefore, they are individually summed.
WI.df = read.csv("./Data_102418/wi_data10-24-18.csv",
header = TRUE, stringsAsFactors = FALSE, sep=",")
WI.Trial.df = WI.df[2:21] #pull-out transects
WI.Trial.df[WI.Trial.df >= 1] = 1 #set track counts to 1
WI.Trial.df$Counts = rowSums(WI.Trial.df, na.rm=T) #sum total successes by transect
WI.Trial.df$Trials = rowSums(is.na(WI.df[2:21])==FALSE) #count number of trials
WI.df = WI.df %>% #Basically the same as MI above
mutate(Long = Lat,
Lat = Lon,
State = "Wisconsin",
OBS = Hare,
Spp = "Hare",
Data = "Fit",
Counts = WI.Trial.df$Counts, #Sum across trials
Trials = WI.Trial.df$Trials, #trials vary by transect
Site2 = paste("W", 1:nrow(WI.df), sep=".")) %>% #Site identifier
select(Site2, Long, Lat, State, OBS, Counts, Trials, Data, Spp)
dim(WI.df)## [1] 195 9head(WI.df)## Site2 Long Lat State OBS Counts Trials Data Spp
## 1 W.1 -90.72816 44.36486 Wisconsin 1 0 8 Fit Hare
## 2 W.2 -90.84052 44.27223 Wisconsin 0 0 6 Fit Hare
## 3 W.3 -90.82058 44.45853 Wisconsin 1 5 8 Fit Hare
## 4 W.4 -92.11831 45.44736 Wisconsin 0 0 6 Fit Hare
## 5 W.5 -91.22783 45.02406 Wisconsin 0 0 6 Fit Hare
## 6 W.6 -90.01462 44.21396 Wisconsin 0 0 8 Fit HareCombine & Convert to Points
Fit.df = rbind(MI.df, WI.df)
Fit.pnt = SpatialPointsDataFrame(Fit.df[, c("Long","Lat")], Fit.df)
proj4string(Fit.pnt) = "+proj=longlat +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0"1.2 Domain Boundaries
Defining domain extent, downloading state boundaries, and then creating raster versions for later plotting. Note that the “Lakes” file is used to better identify lake boundaries. A copy of this shapefile is in the folder with this script. The adress to the shapefile will need to be updated in this chunk.
World = map("world",
fill = TRUE,
plot = FALSE)
IDs = sapply(strsplit(World$names, ":"), function(x) x[1])
LL84 = "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
WorldP = map2SpatialPolygons(World, IDs = IDs,
proj4string = CRS(projection(LL84)))
#Add a dataframe
pid = sapply(slot(WorldP, "polygons"),
function(x) slot(x, "ID"))
p.df = data.frame( ID=1:length(WorldP),
row.names = pid)
World = SpatialPolygonsDataFrame(WorldP, p.df)
World = spTransform(World, LL84)
World = gBuffer(World, width = 0, byid = F)
#States
States = map("state",
fill = TRUE,
plot = FALSE)
IDs = sapply(strsplit(States$names, ":"), function(x) x[1])
LL84 = "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
StatesP = map2SpatialPolygons(States, IDs = IDs,
proj4string = CRS(projection(LL84)))
#Add a dataframe
pid = sapply(slot(StatesP, "polygons"),
function(x) slot(x, "ID"))
p.df = data.frame( ID=1:length(StatesP),
row.names = pid)
States = SpatialPolygonsDataFrame(StatesP, p.df)
States = spTransform(States, LL84)
Lakes = readOGR(dsn = "C:/Users/humph173/Documents/Michigan_State/Marten/ArcWork/Lakes",
layer = "Lake_2ks",
stringsAsFactors = FALSE)## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\humph173\Documents\Michigan_State\Marten\ArcWork\Lakes", layer: "Lake_2ks"
## with 1 features
## It has 15 fields
## Integer64 fields read as strings: OBJECTID Id area InPoly_FID SimPgnFlagLakesLL = spTransform(Lakes, proj4string(States))
Ext = c(-93.173415, -81.932841, 41.217132, 47.740649)
Domain = crop(States, Ext)
Domain$Name = str_cap_words(rownames(Domain@data))
Water0 = as(extent(Domain), "SpatialPolygons")
p.df = data.frame(ID=1:length(Water0))
Water0 = SpatialPolygonsDataFrame(Water0, p.df, match.ID = F)
proj4string(Water0) = proj4string(Domain)
Water1 = gDifference(Water0, spTransform(Lakes, proj4string(Water0)))
DomLLU = gUnaryUnion(Domain)
#Rasterized version
Ras = raster(res = 0.02, ext = extent(DomLLU),
crs = proj4string(DomLLU))
Domain.r = rasterize(DomLLU, Ras,
field = 0,
background = NA)
Domain.full = rasterize(Water0, Ras,
field = 0,
background = NA)
#Point grid version
Grd.pnt = rasterToPoints(Domain.full, spatial = TRUE)
Grd.pnt@data = Grd.pnt@data %>%
mutate(Long = Grd.pnt@coords[,1],
Lat = Grd.pnt@coords[,2],
Spp = "Grid") %>%
select(-layer)
#Dense for hi-res prediction
Zext = as(raster::extent(-92, -91, 45, 46), "SpatialPolygons")
proj4string(Zext) = proj4string(States)
Z_samp = disaggregate(crop(Domain.r, Zext), fact=6)
zPred.pnts = rasterToPoints(Z_samp, sp=T)
zPred.pnts@data = zPred.pnts@data %>%
mutate(Long = zPred.pnts@coords[,1],
Lat = zPred.pnts@coords[,2],
Spp = "Zoom") %>%
select(-layer)
Grd.pnt.df = rbind(Grd.pnt@data, zPred.pnts@data)
Grd.pnt = SpatialPointsDataFrame(Grd.pnt.df[,c("Long","Lat")], Grd.pnt.df)
proj4string(Grd.pnt) = proj4string(States)Quick Plot to Check Data Locations used to fit Model.
rng = seq(0, 255, 1)
mCols = brewer.pal(11, "RdYlBu")[-6]
cr0 = rev(colorRampPalette((mCols))(n = 256))
cr = colorRampPalette(c("tan", cr0),
bias = 1, space = "rgb")
MyMatrix = matrix(nrow=7, ncol=2)
rownames(MyMatrix) = rownames(coordinates(Domain))
MyMatrix[,1] = c(-89.07686, -86.11260, -91.91021, -83.82014, -91.31258, -83.46125, -89.51171)
MyMatrix[,2] = c(41.80000, 41.36514, 42.28710, 43.0, 47.5041, 41.61285, 43.63285)
Hare.LLpnt = spTransform(Fit.pnt, proj4string(States))
levelplot(Domain.r,
margin = FALSE,
xlab = NULL,
ylab = NULL,
maxpixels = 1e5,
col.regions = cr, at = rng,
colorkey = FALSE, par.strip.text = list(fontface='bold', cex=1.5),
par.settings = list(axis.line = list(col = "black"),
strip.background = list(col = 'transparent'),
strip.border = list(col = 'transparent')),
scales = list(cex = 1.25)) +
latticeExtra::layer(sp.polygons(Domain, col = "black", lwd = 0.5)) +
latticeExtra::layer(sp.polygons(LakesLL, fill = "lightblue", col = "transparent", lwd = 0.5)) +
latticeExtra::layer(sp.polygons(Hare.LLpnt , col = "red", pch=factor(Hare.LLpnt$OBS), cex = 1)) +
latticeExtra::layer(sp.text(MyMatrix, txt = Domain$Name,
pos =c(1,3,3,2,2,1,1),
col="black",font=list(face="bold"), cex=1)) +
latticeExtra::layer({SpatialPolygonsRescale(layout.north.arrow(),
offset = c(-83, 45.5),
scale = 2)})
Project Everything to KM
nProj = "+proj=utm +zone=16 +datum=NAD83 +units=km +no_defs +ellps=GRS80 +towgs84=0,0,0"
LakesKM = spTransform(Lakes, nProj)
FocalKM = spTransform(Water1, nProj)
DomP = spTransform(Domain, nProj)
DomPU = gUnaryUnion(DomP)
Fit.pntP = spTransform(Fit.pnt, nProj)
Grd.pntP = spTransform(Grd.pnt, nProj)2 Construct Mesh
These values are scaled to correspond to the geographic projection (kilometers)
max.edge = 8 #Make the outer edge length 8km
bound.outer = 75 #Outer extension can be 75km
Fit.pntPX = subset(Fit.pntP, Spp == "Hare" & Data == "Fit") #Only training data
bdry = inla.sp2segment(DomPU) #Formatting boundary for r-INLA
mesh = inla.mesh.2d(boundary = bdry, #Boundary
loc = Fit.pntPX, #Fit to point locations
max.edge = c(1, 5)*max.edge, #mesh size specifications
cutoff = 10,#8
min.angle = 25,
offset = c(max.edge, bound.outer))mesh$n #number of nodes## [1] 3792plot(mesh, lwd=0.5)
3 Validation Data
Wisconsin
WIV.df = read.csv("./WI_HareValid2000-18.csv",
header = TRUE, stringsAsFactors = FALSE, sep=",")
WIV.df = WIV.df %>%
select(X, Y, HARESA_12_, HARESA_13_, HARESA_14_)
WIV.df[is.na(WIV.df)] = 0
WIV.df = WIV.df %>%
mutate(Site2 = paste("WIV_", 1:dim(WIV.df)[1]),
Long = X,
Lat = Y,
Source = "WI",
Spp = "Valid",
State = "Wisconsin",
Data = "Valid",
Counts = (HARESA_12_ + HARESA_13_ + HARESA_14_),
OBS = ifelse(Counts >= 1, 1, 0),
Trials = 1) %>%
select(Site2, Long, Lat, State, OBS, Counts, Trials, Data, Spp, Source)
dim(WIV.df)## [1] 51 10head(WIV.df)## Site2 Long Lat State OBS Counts Trials Data Spp Source
## 1 WIV_ 1 459786.6 640333.1 Wisconsin 0 0 1 Valid Valid WI
## 2 WIV_ 2 459119.3 623022.5 Wisconsin 0 0 1 Valid Valid WI
## 3 WIV_ 3 448003.0 645595.7 Wisconsin 1 20 1 Valid Valid WI
## 4 WIV_ 4 544112.4 409977.9 Wisconsin 1 40 1 Valid Valid WI
## 5 WIV_ 5 429676.8 692700.4 Wisconsin 1 16 1 Valid Valid WI
## 6 WIV_ 6 323853.0 603638.8 Wisconsin 1 58 1 Valid Valid WI#NAD_1983_HARN_Wisconsin_TM
WIval.nproj = "+proj=tmerc +lat_0=0 +lon_0=-90 +k=0.9996 +x_0=520000 +y_0=-4480000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
WIv.pnt = SpatialPointsDataFrame(WIV.df[, c("Long","Lat")], WIV.df)
proj4string(WIv.pnt) = WIval.nproj
WIv.pnt = spTransform(WIv.pnt, nProj)Hare (Sault Tribe)
library("readxl")
files = list.files(path="C:/Users/humph173/Documents/Michigan_State/Eric/Pellets",
pattern="*.xlsx", full.names=T, recursive=FALSE)
for(i in 1:length(files)) {
tmp.ex = as.data.frame(read_excel(files[i]))[,1:2]
names(tmp.ex) = c("X", "Y")
if(i == 1){Pellet.df = tmp.ex}
else{Pellet.df = rbind(Pellet.df, tmp.ex)}
}## New names:
## * `` -> ...3
## New names:
## * `` -> ...3STV.df = Pellet.df %>%
mutate(Site2 = paste("STV_", 1:dim(Pellet.df)[1]),
Long = X,
Lat = Y,
Source = "ST",
Spp = "Valid",
State = "Michigan",
Data = "Valid",
Counts = 1,
OBS = 1,
Trials = 1) %>%
select(Site2, Long, Lat, State, OBS, Counts, Trials, Data, Spp, Source)
head(STV.df)## Site2 Long Lat State OBS Counts Trials Data Spp Source
## 1 STV_ 1 649079.4 5120467 Michigan 1 1 1 Valid Valid ST
## 2 STV_ 2 661239.2 5135238 Michigan 1 1 1 Valid Valid ST
## 3 STV_ 3 661243.7 5135231 Michigan 1 1 1 Valid Valid ST
## 4 STV_ 4 661246.0 5135232 Michigan 1 1 1 Valid Valid ST
## 5 STV_ 5 661253.6 5135214 Michigan 1 1 1 Valid Valid ST
## 6 STV_ 6 661249.1 5135212 Michigan 1 1 1 Valid Valid STdim(STV.df)## [1] 847 10MnProj = "+proj=utm +zone=16 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"
STV.pnt = SpatialPointsDataFrame(STV.df[, c("Long","Lat")], STV.df)
proj4string(STV.pnt) = MnProj
STV.pnt = spTransform(STV.pnt, nProj)Predator Prey Data
PPV.df = read.csv("C:/Users/humph173/Documents/Michigan_State/Marten/Data/PredPrey/Pred_prey_Comb_121218.csv",
header = TRUE, sep=",")
PPV.df = PPV.df %>%
mutate(Site2 = paste("PPV_", 1:dim(PPV.df)[1]),
Long = Easting,
Lat = Northing,
Source = "PP",
Spp = "Valid",
State = "Michigan",
Data = "Valid",
Counts = 1,
OBS = 1,
Trials = 1) %>%
select(Site2, Long, Lat, State, OBS, Counts, Trials, Data, Spp, Source)
head(PPV.df)## Site2 Long Lat State OBS Counts Trials Data Spp Source
## 1 PPV_ 1 399299 5117030 Michigan 1 1 1 Valid Valid PP
## 2 PPV_ 2 399299 5117030 Michigan 1 1 1 Valid Valid PP
## 3 PPV_ 3 411292 5117450 Michigan 1 1 1 Valid Valid PP
## 4 PPV_ 4 400838 5119454 Michigan 1 1 1 Valid Valid PP
## 5 PPV_ 5 400838 5119454 Michigan 1 1 1 Valid Valid PP
## 6 PPV_ 6 396858 5127841 Michigan 1 1 1 Valid Valid PPdim(PPV.df)## [1] 25939 10PPV.pnt = SpatialPointsDataFrame(PPV.df[, c("Long","Lat")], PPV.df)
proj4string(PPV.pnt) = MnProj
PPV.pnt = spTransform(PPV.pnt, nProj)Lower Survey (Pellet )
LPS.df = read.csv("./2019_pellet_sites.csv",
header = TRUE, sep=",")
LPS.df = LPS.df %>%
mutate(Site2 = paste("PPV_", 1:dim(LPS.df)[1]),
Long = Latitude, #Reversed in data
Lat = Longitude,
Source = "LPP",
Spp = "Valid",
State = "Michigan",
Data = "Valid",
Counts = 1,
OBS = CurrentOcc,
Trials = 1) %>%
filter(plot_ID != "main") %>% #Remove Original Surveys
select(Site2, Long, Lat, State, OBS, Counts, Trials, Data, Spp, Source)
head(LPS.df)## Site2 Long Lat State OBS Counts Trials Data Spp Source
## 1 PPV_ 1 619035.8 519897.4 Michigan 1 1 1 Valid Valid LPP
## 2 PPV_ 2 617098.3 528769.7 Michigan 1 1 1 Valid Valid LPP
## 3 PPV_ 4 620157.4 513569.1 Michigan 1 1 1 Valid Valid LPP
## 4 PPV_ 5 609769.3 522842.3 Michigan 0 1 1 Valid Valid LPP
## 5 PPV_ 6 608475.4 518285.8 Michigan 1 1 1 Valid Valid LPP
## 6 PPV_ 7 613170.1 511702.0 Michigan 0 1 1 Valid Valid LPPdim(LPS.df)## [1] 31 10#Long & Lat Reversed & not UTM
#NAD_1983_2011_Michigan_GeoRef_Meters
StateGeoref.Proj = "+proj=omerc +lat_0=45.30916666666666 +lonc=-86 +alpha=337.25556 +k=0.9996 +x_0=2546731.496 +y_0=-4354009.816 +no_uoff +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
LPS.pnt = SpatialPointsDataFrame(LPS.df[, c("Long","Lat")], LPS.df)
proj4string(LPS.pnt) =StateGeoref.Proj
LPS.pnt = spTransform(LPS.pnt, nProj)Update Coordinates and Combine
Fit.pnt.tab = Fit.pntP@data %>%
mutate(Source = "Fit",
Long = Fit.pntP@coords[,1],
Lat = Fit.pntP@coords[,2])
WIv.pnt.tab = WIv.pnt@data %>%
mutate(Long = WIv.pnt@coords[,1],
Lat = WIv.pnt@coords[,2])
PPV.pnt.tab = PPV.pnt@data %>%
mutate(Long = PPV.pnt@coords[,1],
Lat = PPV.pnt@coords[,2])
STV.pnt.tab = STV.pnt@data %>%
mutate(Long = STV.pnt@coords[,1],
Lat = STV.pnt@coords[,2])
LPS.pnt.tab = LPS.pnt@data %>%
mutate(Long = LPS.pnt@coords[,1],
Lat = LPS.pnt@coords[,2])Combine All Data Join data to common dataframe.
hare.df = rbind(Fit.pnt.tab, WIv.pnt.tab, STV.pnt.tab, PPV.pnt.tab, LPS.pnt.tab)
dim(hare.df)## [1] 27188 10hare.df %>%
group_by(State, Data) %>%
summarise(Count = length(Data))## # A tibble: 4 x 3
## # Groups: State [2]
## State Data Count
## <chr> <chr> <int>
## 1 Michigan Fit 125
## 2 Michigan Valid 26817
## 3 Wisconsin Fit 195
## 4 Wisconsin Valid 51Convert to Points
Hare.pnt = SpatialPointsDataFrame(hare.df[, c("Long","Lat")], hare.df)
proj4string(Hare.pnt) = nProjQuick Plot to Check Validation Data Locations used for model validation.
rng = seq(0, 255, 1)
mCols = brewer.pal(11, "RdYlBu")[-6]
cr0 = rev(colorRampPalette((mCols))(n = 256))
cr = colorRampPalette(c("tan", cr0),
bias = 1, space = "rgb")
MyMatrix = matrix(nrow=7, ncol=2)
rownames(MyMatrix) = rownames(coordinates(Domain))
MyMatrix[,1] = c(-89.07686, -86.11260, -91.91021, -83.82014, -91.31258, -83.46125, -89.51171)
MyMatrix[,2] = c(41.80000, 41.36514, 42.28710, 43.0, 47.5041, 41.61285, 43.63285)
Hare.LLpnt = spTransform(Hare.pnt, proj4string(States))
Hare.LLpnt = subset(Hare.LLpnt, Data == "Valid")
levelplot(Domain.r,
margin = FALSE,
xlab = NULL,
ylab = NULL,
maxpixels = 1e5,
col.regions = cr, at = rng,
colorkey = FALSE, par.strip.text = list(fontface='bold', cex=1.5),
par.settings = list(axis.line = list(col = "black"),
strip.background = list(col = 'transparent'),
strip.border = list(col = 'transparent')),
scales = list(cex = 1.25)) +
latticeExtra::layer(sp.polygons(Domain, col = "black", lwd = 0.5)) +
latticeExtra::layer(sp.polygons(LakesLL, fill = "lightblue", col = "transparent", lwd = 0.5)) +
latticeExtra::layer(sp.polygons(Hare.LLpnt , col = "red", pch=factor(Hare.pnt$Data), cex = 1)) +
latticeExtra::layer(sp.text(MyMatrix, txt = Domain$Name,
pos =c(1,3,3,2,2,1,1),
col="black",font=list(face="bold"), cex=1)) +
latticeExtra::layer({SpatialPolygonsRescale(layout.north.arrow(),
offset = c(-83, 45.5),
scale = 2)})
Combine Observations, Mesh Nodes, and Grid Points
#Node coordinates
dd = as.data.frame(cbind(mesh$loc[,1],
mesh$loc[,2]))
names(dd) = c("Long", "Lat") #name coordinates
dd$OBS = 0 #no hare at these locations
dd$Site2 = paste("N", 1:nrow(dd), sep = ".") #to match with observation data
dd$State = "All"
dd$Spp = "Mesh"
dd$Counts = 0
dd$Trials = 0
dd$Source = "Mesh"
dd$Data = "Fit"
#Hare Obs
hare.set = Hare.pnt@data %>%
mutate(Long = Hare.pnt@coords[,1],
Lat = Hare.pnt@coords[,2])
#Grid
grid.set = Grd.pntP@data %>%
mutate(Long = Grd.pntP@coords[,1],
Lat = Grd.pntP@coords[,2],
Spp = Spp,
OBS = 0,
State = "All",
Counts = 0,
Trials = 0,
Data = "grid",
Source = "Grid",
Site2 = paste("G", 1:nrow(Grd.pntP@data), sep = "."))
All.pnts = rbind(hare.set, dd, grid.set)
All.pnts = SpatialPointsDataFrame(All.pnts[, c("Long","Lat")], All.pnts)
proj4string(All.pnts) = nProj4 Extract Covariates
Forest1km = raster("C:/Users/humph173/Documents/Michigan_State/Sean/Loop_020819/Forest1km.grd")
mSnow5Yr = raster("./Hare1/Mean5yrSnow.tif")
mxTemp = raster("./Hare1/meanTMAX.tif")
All.pnts$mSnow5yrE = extract(mSnow5Yr,
spTransform(All.pnts,
CRS(proj4string(mSnow5Yr))),
method="simple")
All.pnts$mSnow5yrE[is.na(All.pnts$mSnow5yrE)] = mean(All.pnts$mSnow5yrE, na.rm=T)
All.pnts$mxTempE = extract(mxTemp,
spTransform(All.pnts,
CRS(proj4string(mxTemp))),
method="simple")
All.pnts$mxTempE[is.na(All.pnts$mxTempE)] = mean(All.pnts$mxTempE, na.rm=T)
All.pnts$Forest1kmE = extract(Forest1km,
spTransform(All.pnts,
CRS(proj4string(Forest1km))),
method="simple")
All.pnts$Forest1kmE[is.na(All.pnts$Forest1kmE)] = mean(All.pnts$Forest1kmE, na.rm=T)
All.pnts$Forest1kmE = round(All.pnts$Forest1kmE/100, 1)
# 2019 Lower ###################
mSnow5Yr2019 = raster("./Mean5yrlp2019.tif")
mxTemp2019 = raster("./meanTMAXLP2019.tif")
All.pnts$mSnow5yrE2 = extract(mSnow5Yr2019,
spTransform(All.pnts,
CRS(proj4string(mSnow5Yr2019))),
method="simple")
All.pnts$mSnow5yrE2[is.na(All.pnts$mSnow5yrE2)] = mean(All.pnts$mSnow5yrE2, na.rm=T)
All.pnts$mxTempE2 = extract(mxTemp2019,
spTransform(All.pnts,
CRS(proj4string(mxTemp2019))),
method="simple")
All.pnts$mxTempE2[is.na(All.pnts$mxTempE2)] = mean(All.pnts$mxTempE2, na.rm=T)
#Select Variable Year
All.pnts$mSnow5yrE = ifelse(All.pnts$Source == "LPP", All.pnts$mSnow5yrE2, All.pnts$mSnow5yrE)
All.pnts$mxTempE2 = ifelse(All.pnts$Source == "LPP", All.pnts$mxTempE2, All.pnts$mxTempE)
#Clean up frame
All.pnts@data = All.pnts@data %>% select(-c(mSnow5yrE2, mxTempE2))4.1 Earth Environment Variables
Land Cover, Full Conncensus verision: http://www.earthenv.org/landcover
EE.files = list.files(path="E:/EarthEnv",
pattern="*.tif", full.names=T, recursive=T)
EE.stack = stack(EE.files)
EE.df = as.data.frame(
extract(EE.stack,
spTransform(All.pnts,
proj4string(EE.stack)),
method = "simple"))
names(EE.df) = paste("Type", 1:12, sep=".")
EE.dfs = cbind(All.pnts@data %>% select(mSnow5yrE, mxTempE), EE.df)
for(i in 1:dim(EE.dfs)[2]){
EE.dfs[,i] = scale(EE.dfs[,i], scale=T, center=T)
}
names(EE.dfs) = paste("s", names(EE.dfs), sep=".")
All.pnts@data = cbind(All.pnts@data, EE.dfs)
#Lables
EE.names = read.csv("E:/EarthEnv/EE_lables.csv",
header = TRUE, sep=",")
EE.names## Type NAME
## 1 1 Evergreen/Deciduous Needleleaf Trees
## 2 2 Evergreen Broadleaf Trees
## 3 3 Deciduous Broadleaf Trees
## 4 4 Mixed/Other Trees
## 5 5 Shrubs
## 6 6 Herbaceous Vegetation
## 7 7 Cultivated and Managed Vegetation
## 8 8 Regularly Flooded Vegetation
## 9 9 Urban/Built-up
## 10 10 Snow/Ice
## 11 11 Barren
## 12 12 Open WaterIdentify the UP Identify locations from the U.P. based on a well-defined county boundaries
UP = readOGR(dsn = "C:/Users/humph173/Documents/Michigan_State/SLP_Beam_Diam/Counties_v17a",
layer = "MI_UP",
stringsAsFactors = FALSE)## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\humph173\Documents\Michigan_State\SLP_Beam_Diam\Counties_v17a", layer: "MI_UP"
## with 15 features
## It has 15 fields
## Integer64 fields read as strings: OBJECTID FIPSNUMXX = readOGR(dsn = "C:/Users/humph173/Documents/Michigan_State/Sean/Hare1/CoordTest",
layer = "Test7b",
stringsAsFactors = FALSE)## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\humph173\Documents\Michigan_State\Sean\Hare1\CoordTest", layer: "Test7b"
## with 35 features
## It has 13 fields
## Integer64 fields read as strings: FID_1 Site CurrentOcc Site_ID BUFF_DISTUPp = gUnaryUnion(spTransform(UP, proj4string(All.pnts)))
All.pnts$UP = is.na(over(All.pnts, UPp))Update State Designation to Include UP Cleaning up labels
All.pnts$StateUP = ifelse(All.pnts$UP == FALSE & All.pnts$Spp == "Hare", "Mich.UP", All.pnts$State)
All.pnts$StateUPW = ifelse(All.pnts$StateUP == "Mich.UP", "Wisconsin", All.pnts$StateUP)
All.pnts$Domain = over(All.pnts, DomPU)
levels(factor(All.pnts$StateUP))## [1] "All" "Mich.UP" "Michigan" "Wisconsin"levels(factor(All.pnts$StateUPW))## [1] "All" "Michigan" "Wisconsin"4.2 Seperate Datasets
Hare.mod = subset(All.pnts, Spp == "Hare") #Observations only
OneTier.mod = subset(All.pnts, Spp == "Hare" | Spp == "Mesh" & is.na(Domain) == FALSE)
Mesh.mod = subset(All.pnts, Spp == "Mesh" & is.na(Domain) == FALSE) #Mesh locations excluding buffer extension
HareMesh.mod = subset(All.pnts, Data == "Fit") #Observations excluding Grid points for prediction/plotting
HareMesh.mod$ID = 1:nrow(HareMesh.mod@data)
Validation.set = subset(All.pnts, Data == "Valid")
Grd.pnts = subset(All.pnts, Spp == "Grid") #Grid loactions for prediction/plotting
zGrd.pnts = subset(All.pnts, Spp == "Zoom") #zoomed locations for prediction/plottingAdd Background Points to Validation Set
Validation.set2 = spRbind(Validation.set, Mesh.mod)
Validation.set2@data %>%
group_by(OBS) %>%
summarise(Cnt =length(OBS))## # A tibble: 2 x 2
## OBS Cnt
## <dbl> <int>
## 1 0 3061
## 2 1 26857Validation.set2$dups = duplicated(Validation.set2@data[,c("Long","Lat")])
Validation.set2 = subset(Validation.set2, dups == FALSE)5 Colinearty Check
Index value needs to be below 30. This suggest colineartity will not be an issue.
library(perturb)##
## Attaching package: 'perturb'## The following object is masked from 'package:raster':
##
## reclassifyColin.df = HareMesh.mod@data %>% #select covariates
select(mSnow5yrE, mxTempE, Forest1kmE, paste("s.Type", 1:12, sep="."))
CorCov = cor(Colin.df) #calculate correlation
corrplot(CorCov) #view correlation table
CI = colldiag(CorCov) #Apply metric
CI## Condition
## Index Variance Decomposition Proportions
## intercept mSnow5yrE mxTempE Forest1kmE s.Type.1 s.Type.2 s.Type.3
## 1 1.000 0.000 0.000 0.001 0.000 0.004 0.000 0.002
## 2 1.564 0.000 0.001 0.000 0.000 0.003 0.001 0.001
## 3 1.745 0.070 0.000 0.000 0.000 0.011 0.060 0.000
## 4 2.243 0.019 0.000 0.000 0.000 0.000 0.202 0.000
## 5 2.362 0.007 0.000 0.000 0.000 0.019 0.024 0.001
## 6 2.507 0.008 0.000 0.000 0.000 0.008 0.001 0.002
## 7 2.632 0.001 0.000 0.000 0.000 0.005 0.001 0.000
## 8 2.932 0.008 0.000 0.000 0.000 0.058 0.429 0.002
## 9 3.243 0.008 0.000 0.000 0.000 0.263 0.096 0.023
## 10 4.008 0.516 0.001 0.001 0.000 0.006 0.103 0.001
## 11 5.906 0.025 0.001 0.017 0.000 0.370 0.008 0.099
## 12 7.926 0.044 0.000 0.012 0.000 0.205 0.003 0.146
## 13 8.983 0.069 0.021 0.116 0.003 0.032 0.038 0.044
## 14 23.054 0.046 0.058 0.001 0.963 0.005 0.006 0.421
## 15 32.157 0.179 0.918 0.852 0.033 0.011 0.028 0.259
## s.Type.4 s.Type.5 s.Type.6 s.Type.7 s.Type.8 s.Type.9 s.Type.10 s.Type.11
## 1 0.002 0.000 0.002 0.001 0.000 0.000 0.000 0.000
## 2 0.000 0.003 0.004 0.004 0.026 0.005 0.000 0.024
## 3 0.000 0.043 0.001 0.000 0.014 0.042 0.012 0.010
## 4 0.000 0.021 0.001 0.001 0.058 0.115 0.101 0.076
## 5 0.000 0.205 0.001 0.000 0.021 0.006 0.431 0.004
## 6 0.000 0.158 0.000 0.000 0.009 0.001 0.200 0.411
## 7 0.000 0.095 0.000 0.000 0.589 0.000 0.004 0.190
## 8 0.001 0.056 0.000 0.000 0.000 0.304 0.005 0.000
## 9 0.005 0.194 0.003 0.001 0.001 0.138 0.001 0.001
## 10 0.001 0.092 0.000 0.003 0.172 0.029 0.183 0.248
## 11 0.112 0.002 0.119 0.016 0.002 0.012 0.014 0.002
## 12 0.079 0.005 0.641 0.070 0.005 0.002 0.013 0.000
## 13 0.245 0.067 0.055 0.126 0.002 0.033 0.010 0.003
## 14 0.350 0.013 0.172 0.190 0.012 0.269 0.009 0.017
## 15 0.204 0.046 0.001 0.587 0.088 0.045 0.016 0.014
## s.Type.12
## 1 0.000
## 2 0.004
## 3 0.001
## 4 0.000
## 5 0.000
## 6 0.000
## 7 0.000
## 8 0.001
## 9 0.000
## 10 0.004
## 11 0.001
## 12 0.004
## 13 0.005
## 14 0.017
## 15 0.963#Remove Forest1kmE
Colin.df = HareMesh.mod@data %>% #select covariates
select(mSnow5yrE, mxTempE, paste("s.Type", 1:12, sep="."))
CorCov = cor(Colin.df) #calculate correlation
corrplot(CorCov) #view correlation table
CI = colldiag(CorCov) #Apply metric
CI## Condition
## Index Variance Decomposition Proportions
## intercept mSnow5yrE mxTempE s.Type.1 s.Type.2 s.Type.3 s.Type.4
## 1 1.000 0.000 0.000 0.001 0.008 0.000 0.005 0.006
## 2 1.391 0.000 0.001 0.000 0.000 0.001 0.004 0.001
## 3 1.531 0.074 0.000 0.000 0.013 0.054 0.001 0.000
## 4 1.996 0.017 0.000 0.000 0.001 0.209 0.000 0.000
## 5 2.103 0.009 0.000 0.000 0.021 0.017 0.004 0.000
## 6 2.243 0.008 0.000 0.000 0.008 0.000 0.008 0.000
## 7 2.349 0.001 0.000 0.000 0.005 0.001 0.002 0.000
## 8 2.597 0.009 0.000 0.001 0.076 0.390 0.011 0.003
## 9 2.813 0.005 0.000 0.000 0.199 0.144 0.056 0.010
## 10 3.611 0.565 0.001 0.000 0.007 0.106 0.006 0.001
## 11 4.937 0.035 0.001 0.016 0.430 0.011 0.207 0.180
## 12 6.705 0.056 0.000 0.014 0.192 0.004 0.244 0.110
## 13 7.570 0.055 0.019 0.113 0.031 0.036 0.133 0.476
## 14 26.851 0.166 0.977 0.855 0.009 0.025 0.319 0.212
## s.Type.5 s.Type.6 s.Type.7 s.Type.8 s.Type.9 s.Type.10 s.Type.11 s.Type.12
## 1 0.001 0.005 0.003 0.000 0.000 0.000 0.000 0.000
## 2 0.001 0.003 0.003 0.024 0.005 0.000 0.024 0.004
## 3 0.043 0.001 0.000 0.016 0.051 0.011 0.011 0.001
## 4 0.030 0.001 0.001 0.057 0.157 0.087 0.067 0.000
## 5 0.183 0.001 0.000 0.019 0.004 0.459 0.005 0.000
## 6 0.146 0.000 0.000 0.006 0.001 0.177 0.452 0.001
## 7 0.099 0.000 0.000 0.613 0.000 0.005 0.170 0.000
## 8 0.083 0.000 0.000 0.001 0.329 0.007 0.000 0.001
## 9 0.201 0.004 0.001 0.001 0.213 0.000 0.001 0.000
## 10 0.105 0.000 0.004 0.175 0.047 0.197 0.254 0.004
## 11 0.001 0.112 0.013 0.003 0.017 0.018 0.003 0.001
## 12 0.004 0.819 0.070 0.006 0.005 0.015 0.001 0.005
## 13 0.065 0.052 0.122 0.001 0.052 0.008 0.002 0.007
## 14 0.039 0.002 0.782 0.079 0.119 0.015 0.010 0.9776 Model Setup
6.1 Spatial Priors and Indicies
Relating observations to mesh locations as a matrix and defining flat spatial priors.
#Relate mesh for detection level
locs = cbind(Hare.mod@coords[,1], Hare.mod@coords[,2]) #point locations
A.det = inla.spde.make.A(mesh, #the mesh
alpha = 2, #default setting
loc=locs) #our locations
#Relate mesh for covariate level
HR.df = OneTier.mod@data
HR.df$Trials = ifelse(HR.df$Spp == "Mesh", 0, HR.df$Trials) #Zero trials for background points
locs = cbind(HR.df[,"Long"], HR.df[,"Lat"])
A.obs = inla.spde.make.A(mesh, #F
alpha = 2,
loc=locs)
#These two include covariates
A.tem = inla.spde.make.A(mesh,
alpha = 2,
loc=locs,
weights = HR.df[,"mxTempE"]) #covariate to model as SVC
A.sno = inla.spde.make.A(mesh,
alpha = 2,
loc=locs,
weights = HR.df[,"mSnow5yrE"]) #covariate to model as SVC
spde = inla.spde2.pcmatern(mesh,
prior.range = c(1500, 0.5),
prior.sigma = c(5, 0.01))
#Create index to track locations of mesh nodes
field.env = inla.spde.make.index("field.env", spde$n.spde) #index for observations
field.tem = inla.spde.make.index("field.tem", spde$n.spde) #index for Temperature SVC
field.sno = inla.spde.make.index("field.sno", spde$n.spde) #index for snow SVC6.2 Organize data
#Data for full model
HR.lst = list(field.env,
field.tem,
field.sno,
list(intercept1 = rep(1,dim(HR.df)[1]),
mxTempE = HR.df[,"s.mxTempE"],
mSnow5yrE = HR.df[,"s.mSnow5yrE"],
Type.1 = HR.df[,"s.Type.1"],
Type.2 = HR.df[,"s.Type.2"],
Type.3 = HR.df[,"s.Type.3"],
Type.4 = HR.df[,"s.Type.4"],
Type.5 = HR.df[,"s.Type.5"],
Type.6 = HR.df[,"s.Type.6"],
Type.7 = HR.df[,"s.Type.7"],
Type.8 = HR.df[,"s.Type.8"],
Type.9 = HR.df[,"s.Type.9"],
Type.10 = HR.df[,"s.Type.10"],
Type.11 = HR.df[,"s.Type.11"],
Type.12 = HR.df[,"s.Type.12"]))
Env.stk = inla.stack(data = list(Y = HR.df$Counts,
Field.trials = HR.df$Trials),
A = list(A.obs, A.tem, A.sno, 1), #spatial matricies, three with covariate weights
effects = HR.lst,
tag = "env.0")
#Non Background version
#Data for just snow svc model (no temp)
sno.lst = list(field.env,
#field.tem,
field.sno,
list(intercept1 = rep(1,dim(HR.df)[1]),
mxTempE = HR.df[,"s.mxTempE"],
mSnow5yrE = HR.df[,"s.mSnow5yrE"],
Type.1 = HR.df[,"s.Type.1"],
Type.2 = HR.df[,"s.Type.2"],
Type.3 = HR.df[,"s.Type.3"],
Type.4 = HR.df[,"s.Type.4"],
Type.5 = HR.df[,"s.Type.5"],
Type.6 = HR.df[,"s.Type.6"],
Type.7 = HR.df[,"s.Type.7"],
Type.8 = HR.df[,"s.Type.8"],
Type.9 = HR.df[,"s.Type.9"],
Type.10 = HR.df[,"s.Type.10"],
Type.11 = HR.df[,"s.Type.11"],
Type.12 = HR.df[,"s.Type.12"]))
Sno.stk = inla.stack(data = list(Y = HR.df$Counts,
Field.trials = HR.df$Trials),
A = list(A.obs, A.sno, 1),
effects = sno.lst,
tag = "sno.0")
#Data for just temp svc model (no snow)
temp.lst = list(field.env,
field.tem,
#field.sno,
list(intercept1 = rep(1,dim(HR.df)[1]),
mxTempE = HR.df[,"s.mxTempE"],
mSnow5yrE = HR.df[,"s.mSnow5yrE"],
Type.1 = HR.df[,"s.Type.1"],
Type.2 = HR.df[,"s.Type.2"],
Type.3 = HR.df[,"s.Type.3"],
Type.4 = HR.df[,"s.Type.4"],
Type.5 = HR.df[,"s.Type.5"],
Type.6 = HR.df[,"s.Type.6"],
Type.7 = HR.df[,"s.Type.7"],
Type.8 = HR.df[,"s.Type.8"],
Type.9 = HR.df[,"s.Type.9"],
Type.10 = HR.df[,"s.Type.10"],
Type.11 = HR.df[,"s.Type.11"],
Type.12 = HR.df[,"s.Type.12"]))
Tem.stk = inla.stack(data = list(Y = HR.df$Counts,
Field.trials = HR.df$Trials),
A = list(A.obs, A.tem, 1), #
effects = temp.lst,
tag = "tem.0")
#Set for non-SVC model
HR2.lst = list(field.env,
# field.tem,
#field.sno,
list(intercept1 = rep(1,dim(HR.df)[1]),
mxTempE = HR.df[,"s.mxTempE"],
mSnow5yrE = HR.df[,"s.mSnow5yrE"],
Type.1 = HR.df[,"s.Type.1"],
Type.2 = HR.df[,"s.Type.2"],
Type.3 = HR.df[,"s.Type.3"],
Type.4 = HR.df[,"s.Type.4"],
Type.5 = HR.df[,"s.Type.5"],
Type.6 = HR.df[,"s.Type.6"],
Type.7 = HR.df[,"s.Type.7"],
Type.8 = HR.df[,"s.Type.8"],
Type.9 = HR.df[,"s.Type.9"],
Type.10 = HR.df[,"s.Type.10"],
Type.11 = HR.df[,"s.Type.11"],
Type.12 = HR.df[,"s.Type.12"]))
FE.stk = inla.stack(data = list(Y = HR.df$Counts,
Field.trials = HR.df$Trials),
A = list(A.obs, 1),
effects = HR2.lst,
tag = "fe.0")6.3 Site Summary
Sites.df = as.data.frame(
HR.df %>%
filter(Spp == "Hare") %>%
group_by(StateUP) %>%
summarise(Count = length(StateUP),
Detection = sum(Counts, na.rm=T),
Transects = sum(Trials, na.rm=T),
Snow = mean(mSnow5yrE),
minSnow = min(mSnow5yrE),
mxSnow = max(mSnow5yrE),
Temp = mean(mxTempE),
mnTemp = min(mxTempE),
mxTemp = max(mxTempE)))
Sites.df## StateUP Count Detection Transects Snow minSnow mxSnow Temp mnTemp
## 1 Mich.UP 51 218 459 18.94902 14.4 21.6 15.67417 14.23175
## 2 Michigan 74 135 666 15.95135 12.8 18.4 17.26791 15.66925
## 3 Wisconsin 195 240 1885 15.36205 5.0 18.6 17.06331 15.76150
## mxTemp
## 1 16.79450
## 2 18.58800
## 3 18.23675#Transects = Trials? Some have NAs...
sum(Sites.df$Transects)## [1] 3010NonDetect = length(which(Fit.df$Counts == 0))
Detect = length(which(Fit.df$Counts > 0))
#Mean detection rate
Detect/NonDetect## [1] 0.6494845Det.sites = Fit.df %>% filter(Counts > 0)
Total.det.trans = sum(Det.sites$Counts)
1- dim(Det.sites)[1]/Total.det.trans## [1] 0.7875211#Landcover summary
LC.sum.pnts = subset(OneTier.mod, Spp == "Hare")
EE2.df = as.data.frame(
extract(EE.stack,
spTransform(LC.sum.pnts,
proj4string(EE.stack)),
method = "simple"))
names(EE2.df) = paste("Type", 1:12, sep=".")
EE2.df = cbind(LC.sum.pnts@data, EE2.df)
Sites2.df = as.data.frame(
EE2.df %>%
group_by(StateUP) %>%
summarise(T1= mean(Type.1, na.rm=T),
T1.min = min(Type.1, na.rm=T),
T1.max = max(Type.1, na.rm=T),
T2= mean(Type.2, na.rm=T),
T2.min = min(Type.2, na.rm=T),
T2.max = max(Type.2, na.rm=T),
T3= mean(Type.3, na.rm=T),
T3.min = min(Type.3, na.rm=T),
T3.max = max(Type.3, na.rm=T),
T4= mean(Type.4, na.rm=T),
T4.min = min(Type.4, na.rm=T),
T4.max = max(Type.4, na.rm=T),
T6= mean(Type.6, na.rm=T),
T6.min = min(Type.6, na.rm=T),
T6.max = max(Type.6, na.rm=T)))
Sites2.df[,2:16] = round(Sites2.df[,2:16], 2)
Sites2.df$Type1 = paste(Sites2.df$T1, " (", Sites2.df$T1.min, ", ", Sites2.df$T1.max,")", sep="")
Sites2.df$Type2 = paste(Sites2.df$T2, " (", Sites2.df$T2.min, ", ", Sites2.df$T2.max,")", sep="")
Sites2.df$Type3 = paste(Sites2.df$T3, " (", Sites2.df$T3.min, ", ", Sites2.df$T3.max,")", sep="")
Sites2.df$Type4 = paste(Sites2.df$T4, " (", Sites2.df$T4.min, ", ", Sites2.df$T4.max,")", sep="")
Sites2.df$Type6 = paste(Sites2.df$T6, " (", Sites2.df$T6.min, ", ", Sites2.df$T6.max,")", sep="")
Sites2.df = Sites2.df %>%
select(StateUP, Type1, Type2, Type3, Type4, Type6)
names(Sites2.df) = c("Type", "Needleleaf", "Evergreen Broadleaf", "Deciduous Broadleaf", "Mixed", "Herbaceous")
Sites2.df = as.data.frame(t(Sites2.df[,2:6]))
names(Sites2.df) = c("Mich.UP", "Michigan", "Wisconsin")
kable(Sites2.df, caption = "Forest", digits=2) %>%
kable_styling("striped", full_width = F) %>%
row_spec(0, font_size = 20) %>%
column_spec(1, bold = T)| Mich.UP | Michigan | Wisconsin | |
|---|---|---|---|
| Needleleaf | 6.69 (0, 53) | 6.31 (0, 44) | 4.52 (0, 44) |
| Evergreen Broadleaf | 0 (0, 0) | 0 (0, 0) | 0.05 (0, 6) |
| Deciduous Broadleaf | 19.71 (0, 67) | 36.3 (0, 100) | 31.47 (0, 100) |
| Mixed | 64.9 (0, 100) | 47.69 (0, 100) | 21.05 (0, 100) |
| Herbaceous | 0.47 (0, 6) | 2.23 (0, 14) | 2.47 (0, 17) |
6.4 Define Spatial Barriers
Article (Bakka 2019): https://www.sciencedirect.com/science/article/pii/S221167531830099X See tutorial: https://haakonbakka.bitbucket.io/btopic107.html
tl = length(mesh$graph$tv[,1])
posTri = matrix(0, tl, 2)
for (t in 1:tl){
temp = mesh$loc[mesh$graph$tv[t, ], ]
posTri[t,] = colMeans(temp)[c(1,2)]
}
posTri = SpatialPoints(posTri)
proj4string(posTri) = proj4string(FocalKM)
normal = over(DomPU, posTri, returnList=T)
normal = unlist(normal)
barrier.triangles = setdiff(1:tl, normal)
poly.barrier = inla.barrier.polygon(mesh, barrier.triangles)Plot Mesh with Barriers
plot(poly.barrier, col='lightgray')
plot(mesh, main="Mesh", lwd=1, edge.color = "darkgray", add=T)
plot(Fit.pntPX, col="red", pch=19, add=T, cex=0.75)
6.5 Spatial Index and priors
barrier.model1 = inla.barrier.pcmatern(mesh, #Observations
barrier.triangles = barrier.triangles, #boundaries
prior.range = c(400, 0.5),
prior.sigma = c(5, 0.01))
barrier.model2 = inla.barrier.pcmatern(mesh, #for temperature svc
barrier.triangles = barrier.triangles,
prior.range = c(1000, 0.5),
prior.sigma = c(1, 0.01))
barrier.model3 = inla.barrier.pcmatern(mesh, #for snow svc
barrier.triangles = barrier.triangles,
prior.range = c(400, 0.5),
prior.sigma = c(1, 0.01))
spde = inla.spde2.pcmatern(mesh, #for comparison to barrier models
prior.range = c(400, 0.5),
prior.sigma = c(5, 0.01))
#save(list=c("Env.stk", "Tem.stk", "Sno.stk", "spde", "barrier.model1", "barrier.model2", "barrier.model3"), file="C:/Users/humph173/Documents/Michigan_State/Sean/HPC/DEC1419/Barrier.RData")7 Run Models
7.1 Full Model
Climate variables as SVC with forest fixed effects and barriers.
#Previously run
load("~/Michigan_State/Sean/HPC/Dec1419/Barrier_RES121419.RData")
load("~/Michigan_State/Sean/HPC/Dec1419/Barrier2_RES121419.RData")
load("~/Michigan_State/Sean/HPC/Dec1419/Barrier_intRES121419A.RData")JFrm0 = Y ~ -1 + intercept1 +
f(field.env,
model=barrier.model1) +
f(field.tem,
model=barrier.model2) +
f(field.sno,
model=barrier.model3) +
mxTempE + mSnow5yrE +
Type.1 + Type.2 + Type.3 + Type.4 + Type.6
#thetJ = SVC.mod$internal.summary.hyperpar$mean
thetaJ = c(7.4849618, -0.4471538, 8.4521755, -1.3595618, 8.5256066, -1.5158039)
SVC.Full = inla(JFrm0, #formula
data = inla.stack.data(Env.stk),
family = "binomial",
verbose = TRUE,
Ntrials = inla.stack.data(Env.stk)$Field.trials, #number trials
control.fixed = list(prec = 0.001,
prec.intercept = 0.0001),
control.predictor = list(
A = inla.stack.A(Env.stk),
compute = TRUE,
link = 1),
#control.mode = list(restart = TRUE, theta = thetaJ),
control.inla = list(strategy="gaussian",
int.strategy = "eb"),
control.results = list(return.marginals.random = TRUE,
return.marginals.predictor = TRUE),
control.compute=list(dic = TRUE, cpo = TRUE, waic = TRUE)) summary(SVC.Full)##
## Call:
## c("inla(formula = JFrm0, family = \"binomial\", data = inla.stack.data(Env.stk), ", " Ntrials = inla.stack.data(Env.stk)$Field.trials, verbose = TRUE, ", " control.compute = list(dic = TRUE, cpo = TRUE, waic = TRUE), ", " control.predictor = list(A = inla.stack.A(Env.stk), compute = TRUE, ", " link = 1), control.inla = list(strategy = \"gaussian\", ", " int.strategy = \"eb\"), control.results = list(return.marginals.random = TRUE, ", " return.marginals.predictor = TRUE), control.fixed = list(prec = 0.1, ", " prec.intercept = 0.01), num.threads = 8)")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 2.6560 337.1147 1.8531 341.6238
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept1 -3.7326 0.4505 -4.6170 -3.7326 -2.8489 -3.7326 0
## mxTempE -0.5828 0.2913 -1.1546 -0.5828 -0.0114 -0.5828 0
## mSnow5yrE 0.4333 0.2996 -0.1549 0.4332 1.0209 0.4333 0
## Type.1 0.1552 0.0665 0.0248 0.1552 0.2856 0.1552 0
## Type.2 0.1211 0.2202 -0.3113 0.1210 0.5531 0.1211 0
## Type.3 0.3174 0.1487 0.0254 0.3174 0.6092 0.3174 0
## Type.4 0.2686 0.1363 0.0009 0.2686 0.5360 0.2686 0
## Type.6 0.1437 0.2257 -0.2994 0.1437 0.5865 0.1437 0
##
## Random effects:
## Name Model
## field.env RGeneric2
## field.tem RGeneric2
## field.sno RGeneric2
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant 0.975quant mode
## Theta1 for field.env -1.368 0.6457 -2.678 -1.352 -0.1316 -1.296
## Theta2 for field.env 5.934 0.9120 4.330 5.862 7.8835 5.599
## Theta1 for field.tem -3.923 0.7912 -5.469 -3.929 -2.3472 -3.952
## Theta2 for field.tem 6.982 1.0795 5.029 6.913 9.2603 6.664
## Theta1 for field.sno -4.067 0.6749 -5.388 -4.072 -2.7251 -4.088
## Theta2 for field.sno 6.013 1.0352 4.105 5.964 8.1659 5.788
##
## Expected number of effective parameters(std dev): 15.53(0.00)
## Number of equivalent replicates : 216.97
##
## Deviance Information Criterion (DIC) ...............: 456.27
## Deviance Information Criterion (DIC, saturated) ....: NULL
## Effective number of parameters .....................: 16.19
##
## Watanabe-Akaike information criterion (WAIC) ...: 452.51
## Effective number of parameters .................: 10.17
##
## Marginal log-Likelihood: -250.71
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computed7.2 Full Model + Interaction
JFrm0 = Y ~ -1 + intercept1 +
f(field.env,
model=barrier.model1) +
f(field.tem,
model=barrier.model2) +
f(field.sno,
model=barrier.model3) +
mxTempE + mSnow5yrE + mxTempE*mSnow5yrE +
Type.1 + Type.2 + Type.3 + Type.4 + Type.6
#thetJ = SVC.mod$internal.summary.hyperpar$mean
thetaJ = c(7.4849618, -0.4471538, 8.4521755, -1.3595618, 8.5256066, -1.5158039)
SVCI.Full = inla(JFrm0, #formula
data = inla.stack.data(Env.stk),
family = "binomial",
verbose = TRUE,
Ntrials = inla.stack.data(Env.stk)$Field.trials, #number trials
control.fixed = list(prec = 0.001,
prec.intercept = 0.0001),
control.predictor = list(
A = inla.stack.A(Env.stk),
compute = TRUE,
link = 1),
control.mode = list(restart = TRUE, theta = thetaJ),
control.inla = list(strategy="gaussian",
int.strategy = "eb"),
control.results = list(return.marginals.random = TRUE,
return.marginals.predictor = TRUE),
control.compute=list(dic = TRUE, cpo = TRUE, waic = TRUE)) summary(SVCI.Full)##
## Call:
## c("inla(formula = JFrm0, family = \"binomial\", data = inla.stack.data(Env.stk), ", " Ntrials = inla.stack.data(Env.stk)$Field.trials, verbose = TRUE, ", " control.compute = list(dic = TRUE, cpo = TRUE, waic = TRUE), ", " control.predictor = list(A = inla.stack.A(Env.stk), compute = TRUE, ", " link = 1), control.inla = list(strategy = \"gaussian\", ", " int.strategy = \"eb\"), control.results = list(return.marginals.random = TRUE, ", " return.marginals.predictor = TRUE), control.fixed = list(prec = 0.1, ", " prec.intercept = 0.01), num.threads = 8)")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 3.1060 500.2534 2.1881 505.5475
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept1 -3.7300 0.4574 -4.6280 -3.7300 -2.8328 -3.7300 0
## mxTempE -0.8497 0.3182 -1.4744 -0.8497 -0.2255 -0.8497 0
## mSnow5yrE 0.8216 0.3634 0.1081 0.8216 1.5346 0.8216 0
## Type.1 0.1510 0.0666 0.0202 0.1510 0.2817 0.1510 0
## Type.2 0.1216 0.2217 -0.3137 0.1216 0.5566 0.1216 0
## Type.3 0.3288 0.1484 0.0375 0.3288 0.6200 0.3288 0
## Type.4 0.3032 0.1373 0.0336 0.3032 0.5726 0.3032 0
## Type.6 0.1495 0.2250 -0.2923 0.1495 0.5909 0.1495 0
## mxTempE:mSnow5yrE 0.7333 0.3727 0.0016 0.7333 1.4644 0.7333 0
##
## Random effects:
## Name Model
## field.env RGeneric2
## field.tem RGeneric2
## field.sno RGeneric2
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant 0.975quant mode
## Theta1 for field.env -1.229 0.6121 -2.435 -1.230 -0.0245 -1.231
## Theta2 for field.env 5.808 0.8353 4.277 5.761 7.5530 5.595
## Theta1 for field.tem -3.812 0.8063 -5.347 -3.831 -2.1822 -3.902
## Theta2 for field.tem 7.013 1.0790 5.113 6.925 9.3271 6.605
## Theta1 for field.sno -4.081 0.6354 -5.334 -4.081 -2.8339 -4.082
## Theta2 for field.sno 5.946 0.9430 4.211 5.897 7.9056 5.721
##
## Expected number of effective parameters(std dev): 16.70(0.00)
## Number of equivalent replicates : 201.83
##
## Deviance Information Criterion (DIC) ...............: 453.68
## Deviance Information Criterion (DIC, saturated) ....: NULL
## Effective number of parameters .....................: 17.43
##
## Watanabe-Akaike information criterion (WAIC) ...: 449.16
## Effective number of parameters .................: 10.60
##
## Marginal log-Likelihood: -251.00
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computed7.3 Temperature SVC
Including landcover fixed effects
JFrm0 = Y ~ -1 + intercept1 +
f(field.env,
model=barrier.model1) +
f(field.tem,
model=barrier.model2) +
mxTempE +
Type.1 + Type.2 + Type.3 + Type.4 + Type.6
#thetJ = SVC.mod$internal.summary.hyperpar$mean
thetaJ = c(-0.9392654, 7.1265513, -3.1710490, 7.6613096, -2.3975088, 8.6556755)
SVC.temp = inla(JFrm0, #formula
num.threads=8,
data = inla.stack.data(Tem.stk),
family = "binomial",
verbose = TRUE,
Ntrials = inla.stack.data(Tem.stk)$Field.trials, #number trials
control.fixed = list(prec = 0.1,
prec.intercept = 0.01),
control.predictor = list(
A = inla.stack.A(Tem.stk),
compute = TRUE,
link = 1),
#control.mode = list(restart = TRUE, theta = thetaJ),
control.inla = list(strategy="gaussian",
int.strategy = "eb"),
control.results = list(return.marginals.random = TRUE,
return.marginals.predictor = TRUE),
control.compute=list(dic = TRUE, cpo = TRUE, waic = TRUE)) summary(SVC.temp)##
## Call:
## c("inla(formula = JFrm0, family = \"binomial\", data = inla.stack.data(Tem.stk), ", " Ntrials = inla.stack.data(Tem.stk)$Field.trials, verbose = TRUE, ", " control.compute = list(dic = TRUE, cpo = TRUE, waic = TRUE), ", " control.predictor = list(A = inla.stack.A(Tem.stk), compute = TRUE, ", " link = 1), control.inla = list(strategy = \"gaussian\", ", " int.strategy = \"eb\"), control.results = list(return.marginals.random = TRUE, ", " return.marginals.predictor = TRUE), control.fixed = list(prec = 0.1, ", " prec.intercept = 0.01), num.threads = 8)")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 1.5278 112.9460 1.8173 116.2911
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept1 -3.6001 0.4003 -4.3861 -3.6001 -2.8147 -3.6001 0
## mxTempE -0.7840 0.2290 -1.2335 -0.7840 -0.3348 -0.7840 0
## Type.1 0.1495 0.0660 0.0200 0.1495 0.2789 0.1495 0
## Type.2 0.1159 0.2185 -0.3132 0.1159 0.5446 0.1159 0
## Type.3 0.3213 0.1461 0.0345 0.3213 0.6080 0.3213 0
## Type.4 0.2981 0.1319 0.0392 0.2981 0.5568 0.2981 0
## Type.6 0.1500 0.2241 -0.2901 0.1500 0.5896 0.1500 0
##
## Random effects:
## Name Model
## field.env RGeneric2
## field.tem RGeneric2
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant 0.975quant mode
## Theta1 for field.env -1.718 0.6874 -3.239 -1.640 -0.5803 -1.346
## Theta2 for field.env 5.729 0.8511 4.200 5.672 7.5286 5.466
## Theta1 for field.tem -3.892 0.7462 -5.288 -3.922 -2.3538 -4.031
## Theta2 for field.tem 7.548 1.1895 5.677 7.388 10.2384 6.773
##
## Expected number of effective parameters(std dev): 12.72(0.00)
## Number of equivalent replicates : 265.02
##
## Deviance Information Criterion (DIC) ...............: 455.80
## Deviance Information Criterion (DIC, saturated) ....: NULL
## Effective number of parameters .....................: 13.25
##
## Watanabe-Akaike information criterion (WAIC) ...: 453.68
## Effective number of parameters .................: 9.001
##
## Marginal log-Likelihood: -247.22
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computed7.4 Snow SVC
Including landcover fixed effects
JFrm0 = Y ~ -1 + intercept1 +
f(field.env,
model=barrier.model1) +
#f(field.tem,
# model=barrier.model2) +
f(field.sno,
model=barrier.model3) +
mSnow5yrE +
Type.1 + Type.2 + Type.3 + Type.4 + Type.6
#thetJ = SVC.mod$internal.summary.hyperpar$mean
thetaJ = c(-0.9392654, 7.1265513, -3.1710490, 7.6613096, -2.3975088, 8.6556755)
SVC.sno = inla(JFrm0, #formula
num.threads=8,
data = inla.stack.data(Sno.stk),
family = "binomial",
verbose = TRUE,
Ntrials = inla.stack.data(Sno.stk)$Field.trials, #number trials
control.fixed = list(prec = 0.1,
prec.intercept = 0.01),
control.predictor = list(
A = inla.stack.A(Sno.stk),
compute = TRUE,
link = 1),
#control.mode = list(restart = TRUE, theta = thetaJ),
control.inla = list(strategy="gaussian",
int.strategy = "eb"),
control.results = list(return.marginals.random = TRUE,
return.marginals.predictor = TRUE),
control.compute=list(dic = TRUE, cpo = TRUE, waic = TRUE)) summary(SVC.sno)##
## Call:
## c("inla(formula = JFrm0, family = \"binomial\", data = inla.stack.data(Sno.stk), ", " Ntrials = inla.stack.data(Sno.stk)$Field.trials, verbose = TRUE, ", " control.compute = list(dic = TRUE, cpo = TRUE, waic = TRUE), ", " control.predictor = list(A = inla.stack.A(Sno.stk), compute = TRUE, ", " link = 1), control.inla = list(strategy = \"gaussian\", ", " int.strategy = \"eb\"), control.results = list(return.marginals.random = TRUE, ", " return.marginals.predictor = TRUE), control.fixed = list(prec = 0.1, ", " prec.intercept = 0.01), num.threads = 8)")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 0.8680 108.3998 1.1357 110.4035
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept1 -3.7507 0.2614 -4.2640 -3.7507 -3.2378 -3.7507 0
## mSnow5yrE 0.7771 0.2542 0.2781 0.7771 1.2756 0.7771 0
## Type.1 0.1410 0.0660 0.0114 0.1410 0.2706 0.1410 0
## Type.2 0.1227 0.2221 -0.3134 0.1227 0.5585 0.1227 0
## Type.3 0.3100 0.1470 0.0215 0.3100 0.5983 0.3100 0
## Type.4 0.2760 0.1326 0.0158 0.2760 0.5361 0.2760 0
## Type.6 0.1070 0.2246 -0.3339 0.1070 0.5475 0.1070 0
##
## Random effects:
## Name Model
## field.env RGeneric2
## field.sno RGeneric2
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant 0.975quant mode
## Theta1 for field.env -1.380 0.6238 -2.638 -1.367 -0.1834 -1.323
## Theta2 for field.env 5.602 0.7762 4.181 5.559 7.2215 5.402
## Theta1 for field.sno -4.253 0.6863 -5.617 -4.246 -2.9164 -4.223
## Theta2 for field.sno 5.888 0.8454 4.425 5.811 7.7149 5.529
##
## Expected number of effective parameters(std dev): 13.99(0.00)
## Number of equivalent replicates : 240.82
##
## Deviance Information Criterion (DIC) ...............: 457.37
## Deviance Information Criterion (DIC, saturated) ....: NULL
## Effective number of parameters .....................: 14.59
##
## Watanabe-Akaike information criterion (WAIC) ...: 454.32
## Effective number of parameters .................: 9.429
##
## Marginal log-Likelihood: -248.46
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computed7.5 No climate
Barrier model with landcover fixed effects only
JFrm0 = Y ~ -1 + intercept1 +
f(field.env,
model=barrier.model1) +
#f(field.tem,
# model=barrier.model2) +
#f(field.sno,
# model=barrier.model3) +
#mxTempE + mSnow5yrE +
Type.1 + Type.2 + Type.3 + Type.4 + Type.6
#thetJ = SVC.mod$internal.summary.hyperpar$mean
thetaJ = c(-0.9392654, 7.1265513, -3.1710490, 7.6613096, -2.3975088, 8.6556755)
No.Clim.mod = inla(JFrm0, #formula
num.threads=8,
data = inla.stack.data(Env.stk),
family = "binomial",
verbose = TRUE,
Ntrials = inla.stack.data(Env.stk)$Field.trials, #number trials
control.fixed = list(prec = 0.1,
prec.intercept = 0.01),
control.predictor = list(
A = inla.stack.A(Env.stk),
compute = TRUE,
link = 1),
#control.mode = list(restart = TRUE, theta = thetaJ),
control.inla = list(strategy="gaussian",
int.strategy = "eb"),
control.results = list(return.marginals.random = TRUE,
return.marginals.predictor = TRUE),
control.compute=list(dic = TRUE, cpo = TRUE, waic = TRUE)) summary(No.Clim.mod)##
## Call:
## c("inla(formula = JFrm0, family = \"binomial\", data = inla.stack.data(Env.stk), ", " Ntrials = inla.stack.data(Env.stk)$Field.trials, verbose = TRUE, ", " control.compute = list(dic = TRUE, cpo = TRUE, waic = TRUE), ", " control.predictor = list(A = inla.stack.A(Env.stk), compute = TRUE, ", " link = 1), control.inla = list(strategy = \"gaussian\", ", " int.strategy = \"eb\"), control.results = list(return.marginals.random = TRUE, ", " return.marginals.predictor = TRUE), control.fixed = list(prec = 0.1, ", " prec.intercept = 0.01), num.threads = 8)")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 1.2546 41.0931 1.5536 43.9013
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept1 -3.5615 0.3496 -4.2478 -3.5615 -2.8757 -3.5615 0
## Type.1 0.1278 0.0661 -0.0020 0.1278 0.2574 0.1278 0
## Type.2 0.1248 0.2202 -0.3074 0.1248 0.5567 0.1248 0
## Type.3 0.3282 0.1503 0.0331 0.3282 0.6230 0.3282 0
## Type.4 0.3378 0.1363 0.0702 0.3378 0.6053 0.3378 0
## Type.6 0.0884 0.2242 -0.3517 0.0884 0.5281 0.0884 0
##
## Random effects:
## Name Model
## field.env RGeneric2
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant 0.975quant mode
## Theta1 for field.env -0.6338 0.4216 -1.472 -0.6307 0.189 -0.6193
## Theta2 for field.env 5.8184 0.6853 4.580 5.7744 7.262 5.6153
##
## Expected number of effective parameters(std dev): 16.36(0.00)
## Number of equivalent replicates : 206.02
##
## Deviance Information Criterion (DIC) ...............: 461.50
## Deviance Information Criterion (DIC, saturated) ....: NULL
## Effective number of parameters .....................: 17.05
##
## Watanabe-Akaike information criterion (WAIC) ...: 457.47
## Effective number of parameters .................: 10.72
##
## Marginal log-Likelihood: -248.86
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computed7.6 Full No SVC
Barrier Model with climate and landcover as fixed effects.
JFrm0 = Y ~ -1 + intercept1 +
f(field.env,
model=barrier.model1) +
#f(field.tem,
# model=barrier.model2) +
#f(field.sno,
# model=barrier.model3) +
mxTempE + mSnow5yrE +
Type.1 + Type.2 + Type.3 + Type.4 + Type.6
#thetJ = SVC.mod$internal.summary.hyperpar$mean
thetaJ = c(-0.9392654, 7.1265513, -3.1710490, 7.6613096, -2.3975088, 8.6556755)
Full.noSVC = inla(JFrm0, #formula
num.threads=8,
data = inla.stack.data(Env.stk),
family = "binomial",
verbose = TRUE,
Ntrials = inla.stack.data(Env.stk)$Field.trials, #number trials
control.fixed = list(prec = 0.1,
prec.intercept = 0.01),
control.predictor = list(
A = inla.stack.A(Env.stk),
compute = TRUE,
link = 1),
#control.mode = list(restart = TRUE, theta = thetaJ),
control.inla = list(strategy="gaussian",
int.strategy = "eb"),
control.results = list(return.marginals.random = TRUE,
return.marginals.predictor = TRUE),
control.compute=list(dic = TRUE, cpo = TRUE, waic = TRUE)) summary(Full.noSVC)##
## Call:
## c("inla(formula = JFrm0, family = \"binomial\", data = inla.stack.data(Env.stk), ", " Ntrials = inla.stack.data(Env.stk)$Field.trials, verbose = TRUE, ", " control.compute = list(dic = TRUE, cpo = TRUE, waic = TRUE), ", " control.predictor = list(A = inla.stack.A(Env.stk), compute = TRUE, ", " link = 1), control.inla = list(strategy = \"gaussian\", ", " int.strategy = \"eb\"), control.results = list(return.marginals.random = TRUE, ", " return.marginals.predictor = TRUE), control.fixed = list(prec = 0.1, ", " prec.intercept = 0.01), num.threads = 8)")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 1.6676 40.2546 2.1394 44.0617
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept1 -3.7037 0.2416 -4.1782 -3.7037 -3.2297 -3.7037 0
## mxTempE -0.5344 0.2689 -1.0624 -0.5344 -0.0069 -0.5344 0
## mSnow5yrE 0.4207 0.2858 -0.1404 0.4207 0.9812 0.4207 0
## Type.1 0.1541 0.0659 0.0247 0.1541 0.2833 0.1541 0
## Type.2 0.1137 0.2207 -0.3196 0.1137 0.5467 0.1137 0
## Type.3 0.3188 0.1455 0.0331 0.3188 0.6042 0.3188 0
## Type.4 0.2819 0.1308 0.0251 0.2819 0.5386 0.2819 0
## Type.6 0.1571 0.2246 -0.2839 0.1571 0.5978 0.1571 0
##
## Random effects:
## Name Model
## field.env RGeneric2
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant 0.975quant mode
## Theta1 for field.env -1.442 0.5731 -2.671 -1.396 -0.4332 -1.228
## Theta2 for field.env 6.127 0.9993 4.531 5.998 8.3670 5.496
##
## Expected number of effective parameters(std dev): 13.37(0.00)
## Number of equivalent replicates : 252.12
##
## Deviance Information Criterion (DIC) ...............: 454.92
## Deviance Information Criterion (DIC, saturated) ....: NULL
## Effective number of parameters .....................: 13.94
##
## Watanabe-Akaike information criterion (WAIC) ...: 452.22
## Effective number of parameters .................: 9.128
##
## Marginal log-Likelihood: -246.74
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computed7.7 Temperature SVC (no LC)
Excluding landcover and snow
JFrm0 = Y ~ -1 + intercept1 +
f(field.env,
model=barrier.model1) +
f(field.tem,
model=barrier.model2) +
#f(field.sno,
# model=barrier.model3) +
mxTempE
#Type.1 + Type.2 + Type.3 + Type.4 + Type.6
#thetJ = SVC.mod$internal.summary.hyperpar$mean
thetaJ = c(-0.9392654, 7.1265513, -3.1710490, 7.6613096, -2.3975088, 8.6556755)
No.LC.temp = inla(JFrm0, #formula
num.threads=8,
data = inla.stack.data(Tem.stk),
family = "binomial",
verbose = TRUE,
Ntrials = inla.stack.data(Tem.stk)$Field.trials, #number trials
control.fixed = list(prec = 0.1,
prec.intercept = 0.01),
control.predictor = list(
A = inla.stack.A(Tem.stk),
compute = TRUE,
link = 1),
#control.mode = list(restart = TRUE, theta = thetaJ),
control.inla = list(strategy="gaussian",
int.strategy = "eb"),
control.results = list(return.marginals.random = TRUE,
return.marginals.predictor = TRUE),
control.compute=list(dic = TRUE, cpo = TRUE, waic = TRUE)) summary(No.LC.temp)##
## Call:
## c("inla(formula = JFrm0, family = \"binomial\", data = inla.stack.data(Tem.stk), ", " Ntrials = inla.stack.data(Tem.stk)$Field.trials, verbose = TRUE, ", " control.compute = list(dic = TRUE, cpo = TRUE, waic = TRUE), ", " control.predictor = list(A = inla.stack.A(Tem.stk), compute = TRUE, ", " link = 1), control.inla = list(strategy = \"gaussian\", ", " int.strategy = \"eb\"), control.results = list(return.marginals.random = TRUE, ", " return.marginals.predictor = TRUE), control.fixed = list(prec = 0.1, ", " prec.intercept = 0.01), num.threads = 8)")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 0.8329 127.8330 0.7279 129.3937
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept1 -3.1851 0.5919 -4.3472 -3.1851 -2.0239 -3.1851 0
## mxTempE -0.7884 0.2291 -1.2382 -0.7884 -0.3391 -0.7884 0
##
## Random effects:
## Name Model
## field.env RGeneric2
## field.tem RGeneric2
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant 0.975quant mode
## Theta1 for field.env -1.665 0.6481 -3.036 -1.623 -0.4925 -1.470
## Theta2 for field.env 6.287 1.0625 4.492 6.176 8.6160 5.761
## Theta1 for field.tem -3.602 0.5475 -4.624 -3.627 -2.4604 -3.718
## Theta2 for field.tem 7.383 0.8839 5.828 7.316 9.2926 7.074
##
## Expected number of effective parameters(std dev): 6.931(0.00)
## Number of equivalent replicates : 486.24
##
## Deviance Information Criterion (DIC) ...............: 454.84
## Deviance Information Criterion (DIC, saturated) ....: NULL
## Effective number of parameters .....................: 7.002
##
## Watanabe-Akaike information criterion (WAIC) ...: 451.47
## Effective number of parameters .................: 3.523
##
## Marginal log-Likelihood: -235.59
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computed7.8 Snow SVC (no LC)
Excluding landcover and temperature
JFrm0 = Y ~ -1 + intercept1 +
f(field.env,
model=barrier.model1) +
#f(field.tem,
# model=barrier.model2) +
f(field.sno,
model=barrier.model3) +
mSnow5yrE
#Type.1 + Type.2 + Type.3 + Type.4 + Type.6
#thetJ = SVC.mod$internal.summary.hyperpar$mean
thetaJ = c(-0.9392654, 7.1265513, -3.1710490, 7.6613096, -2.3975088, 8.6556755)
No.LC.sno = inla(JFrm0, #formula
num.threads=8,
data = inla.stack.data(Sno.stk),
family = "binomial",
verbose = TRUE,
Ntrials = inla.stack.data(Sno.stk)$Field.trials, #number trials
control.fixed = list(prec = 0.1,
prec.intercept = 0.01),
control.predictor = list(
A = inla.stack.A(Sno.stk),
compute = TRUE,
link = 1),
#control.mode = list(restart = TRUE, theta = thetaJ),
control.inla = list(strategy="gaussian",
int.strategy = "eb"),
control.results = list(return.marginals.random = TRUE,
return.marginals.predictor = TRUE),
control.compute=list(dic = TRUE, cpo = TRUE, waic = TRUE))summary(No.LC.sno)##
## Call:
## c("inla(formula = JFrm0, family = \"binomial\", data = inla.stack.data(Sno.stk), ", " Ntrials = inla.stack.data(Sno.stk)$Field.trials, verbose = TRUE, ", " control.compute = list(dic = TRUE, cpo = TRUE, waic = TRUE), ", " control.predictor = list(A = inla.stack.A(Sno.stk), compute = TRUE, ", " link = 1), control.inla = list(strategy = \"gaussian\", ", " int.strategy = \"eb\"), control.results = list(return.marginals.random = TRUE, ", " return.marginals.predictor = TRUE), control.fixed = list(prec = 0.1, ", " prec.intercept = 0.01), num.threads = 8)")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 0.8629 92.3857 1.2083 94.4569
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept1 -3.4516 0.2722 -3.9860 -3.4516 -2.9177 -3.4516 0
## mSnow5yrE 0.8464 0.2507 0.3542 0.8464 1.3382 0.8464 0
##
## Random effects:
## Name Model
## field.env RGeneric2
## field.sno RGeneric2
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant 0.975quant mode
## Theta1 for field.env -1.557 0.6656 -2.972 -1.510 -0.3671 -1.340
## Theta2 for field.env 5.876 0.9060 4.151 5.853 7.7061 5.770
## Theta1 for field.sno -4.056 0.6288 -5.248 -4.077 -2.7715 -4.150
## Theta2 for field.sno 6.571 1.1263 4.664 6.459 9.0401 6.043
##
## Expected number of effective parameters(std dev): 7.864(0.00)
## Number of equivalent replicates : 428.52
##
## Deviance Information Criterion (DIC) ...............: 455.21
## Deviance Information Criterion (DIC, saturated) ....: NULL
## Effective number of parameters .....................: 7.961
##
## Watanabe-Akaike information criterion (WAIC) ...: 451.17
## Effective number of parameters .................: 3.812
##
## Marginal log-Likelihood: -236.07
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computed7.9 SPDE Random field
Random field, no bariers
JFrm0 = Y ~ -1 + intercept1 +
f(field.env,
model=spde)
#thetJ = SVC.mod$internal.summary.hyperpar$mean
thetaJ = c(-0.9392654, 7.1265513, -3.1710490, 7.6613096, -2.3975088, 8.6556755)
spde.RF.mod = inla(JFrm0, #formula
num.threads=8,
data = inla.stack.data(Env.stk),
family = "binomial",
verbose = TRUE,
Ntrials = inla.stack.data(Env.stk)$Field.trials, #number trials
control.fixed = list(prec = 0.1,
prec.intercept = 0.01),
control.predictor = list(
A = inla.stack.A(Env.stk),
compute = TRUE,
link = 1),
#control.mode = list(restart = TRUE, theta = thetaJ),
control.inla = list(strategy="gaussian",
int.strategy = "eb"),
control.results = list(return.marginals.random = TRUE,
return.marginals.predictor = TRUE),
control.compute=list(dic = TRUE, cpo = TRUE, waic = TRUE)) summary(spde.RF.mod)##
## Call:
## c("inla(formula = JFrm0, family = \"binomial\", data = inla.stack.data(Env.stk), ", " Ntrials = inla.stack.data(Env.stk)$Field.trials, verbose = TRUE, ", " control.compute = list(dic = TRUE, cpo = TRUE, waic = TRUE), ", " control.predictor = list(A = inla.stack.A(Env.stk), compute = TRUE, ", " link = 1), control.inla = list(strategy = \"gaussian\", ", " int.strategy = \"eb\"), control.results = list(return.marginals.random = TRUE, ", " return.marginals.predictor = TRUE), control.fixed = list(prec = 0.1, ", " prec.intercept = 0.01), num.threads = 8)")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 0.4826 14.9804 1.2919 16.7548
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept1 -3.3244 0.776 -4.848 -3.3244 -1.802 -3.3244 0
##
## Random effects:
## Name Model
## field.env SPDE2 model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant 0.975quant mode
## Range for field.env 989.525 851.4802 231.9230 738.8574 3242.110 468.9798
## Stdev for field.env 1.023 0.4891 0.4036 0.9141 2.274 0.7416
##
## Expected number of effective parameters(std dev): 10.07(0.00)
## Number of equivalent replicates : 334.56
##
## Deviance Information Criterion (DIC) ...............: 455.41
## Deviance Information Criterion (DIC, saturated) ....: NULL
## Effective number of parameters .....................: 10.23
##
## Watanabe-Akaike information criterion (WAIC) ...: 450.26
## Effective number of parameters .................: 4.904
##
## Marginal log-Likelihood: -235.83
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computed7.10 Barrier Random Field
JFrm0 = Y ~ -1 + intercept1 +
f(field.env,
model=barrier.model1)
#thetJ = SVC.mod$internal.summary.hyperpar$mean
thetaJ = c(-0.9392654, 7.1265513, -3.1710490, 7.6613096, -2.3975088, 8.6556755)
barrier.RF.mod = inla(JFrm0, #formula
num.threads=8,
data = inla.stack.data(Env.stk),
family = "binomial",
verbose = TRUE,
Ntrials = inla.stack.data(Env.stk)$Field.trials,
control.fixed = list(prec = 0.1,
prec.intercept = 0.01),
control.predictor = list(
A = inla.stack.A(Env.stk),
compute = TRUE,
link = 1),
#control.mode = list(restart = TRUE, theta = thetaJ),
control.inla = list(strategy="gaussian",
int.strategy = "eb"),
control.results = list(return.marginals.random = TRUE,
return.marginals.predictor = TRUE),
control.compute=list(dic = TRUE, cpo = TRUE, waic = TRUE)) summary(barrier.RF.mod)##
## Call:
## c("inla(formula = JFrm0, family = \"binomial\", data = inla.stack.data(Env.stk), ", " Ntrials = inla.stack.data(Env.stk)$Field.trials, verbose = TRUE, ", " control.compute = list(dic = TRUE, cpo = TRUE, waic = TRUE), ", " control.predictor = list(A = inla.stack.A(Env.stk), compute = TRUE, ", " link = 1), control.inla = list(strategy = \"gaussian\", ", " int.strategy = \"eb\"), control.results = list(return.marginals.random = TRUE, ", " return.marginals.predictor = TRUE), control.fixed = list(prec = 0.1, ", " prec.intercept = 0.01), num.threads = 8)")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 0.4760 42.5353 1.0889 44.1001
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept1 -3.2327 0.5054 -4.2249 -3.2328 -2.2414 -3.2327 0
##
## Random effects:
## Name Model
## field.env RGeneric2
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant 0.975quant mode
## Theta1 for field.env -0.3273 0.4027 -1.089 -0.3401 0.4945 -0.3863
## Theta2 for field.env 6.3626 0.6773 5.176 6.3079 7.8171 6.1061
##
## Expected number of effective parameters(std dev): 10.63(0.00)
## Number of equivalent replicates : 316.89
##
## Deviance Information Criterion (DIC) ...............: 458.80
## Deviance Information Criterion (DIC, saturated) ....: NULL
## Effective number of parameters .....................: 10.82
##
## Watanabe-Akaike information criterion (WAIC) ...: 453.53
## Effective number of parameters .................: 5.309
##
## Marginal log-Likelihood: -237.63
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computed7.11 Non Spatial
No spatial effect, No SVC, simple linear regression
JFrm0 = Y ~ -1 + intercept1 +
mxTempE + mSnow5yrE +
Type.1 + Type.2 + Type.3 + Type.4 + Type.6
#thetJ = NonSVC.mod$internal.summary.hyperpar$mean
thetaJ = c(7.4867414, -0.4521236, 8.4357741, -1.4313818, 8.4668390, -1.4568577)
Linear.mod = inla(JFrm0, #formula
data = inla.stack.data(Env.stk),
family = "binomial",
verbose = TRUE,
Ntrials = inla.stack.data(Env.stk)$Field.trials, #number trials
control.fixed = list(prec = 0.001,
prec.intercept = 0.0001),
control.predictor = list(
A = inla.stack.A(Env.stk),
compute = TRUE,
link = 1),
#control.mode = list(restart = TRUE, theta = thetaJ),
control.inla = list(strategy="gaussian",
int.strategy = "eb"),
control.results = list(return.marginals.random = TRUE,
return.marginals.predictor = TRUE),
control.compute=list(dic = TRUE, cpo = TRUE, waic = TRUE)) summary(Linear.mod)##
## Call:
## c("inla(formula = JFrm0, family = \"binomial\", data = inla.stack.data(Env.stk), ", " Ntrials = inla.stack.data(Env.stk)$Field.trials, verbose = TRUE, ", " control.compute = list(dic = TRUE, cpo = TRUE, waic = TRUE), ", " control.predictor = list(A = inla.stack.A(Env.stk), compute = TRUE, ", " link = 1), control.inla = list(strategy = \"gaussian\", ", " int.strategy = \"eb\"), control.results = list(return.marginals.random = TRUE, ", " return.marginals.predictor = TRUE), control.fixed = list(prec = 0.001, ", " prec.intercept = 1e-04))")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 1.2397 4.3741 1.2374 6.8512
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept1 -3.6670 0.1700 -4.0008 -3.6670 -3.3334 -3.6670 0
## mxTempE -0.4302 0.2438 -0.9089 -0.4302 0.0481 -0.4302 0
## mSnow5yrE 0.4070 0.2719 -0.1269 0.4070 0.9405 0.4070 0
## Type.1 0.1551 0.0646 0.0283 0.1551 0.2818 0.1551 0
## Type.2 0.0976 0.2231 -0.3403 0.0976 0.5352 0.0976 0
## Type.3 0.3373 0.1390 0.0645 0.3373 0.6099 0.3373 0
## Type.4 0.3334 0.1180 0.1018 0.3334 0.5648 0.3334 0
## Type.6 0.2109 0.2231 -0.2270 0.2109 0.6485 0.2109 0
##
## The model has no random effects
##
## The model has no hyperparameters
##
## Expected number of effective parameters(std dev): 8.306(0.00)
## Number of equivalent replicates : 405.72
##
## Deviance Information Criterion (DIC) ...............: 456.58
## Deviance Information Criterion (DIC, saturated) ....: 213.54
## Effective number of parameters .....................: 8.761
##
## Watanabe-Akaike information criterion (WAIC) ...: 456.54
## Effective number of parameters .................: 6.815
##
## Marginal log-Likelihood: -263.34
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computed8 Results
8.1 Compare Parsimony
Models = c("barrier.RF.mod", "Linear.mod", "No.LC.sno",
"SVC.sno", "No.LC.temp", "SVC.temp",
"No.Clim.mod", "Full.noSVC", "SVC.Full", "SVCI.Full")
Type = c("Random Field",
"Non-spatial Regression",
"Snow SVC",
"Snow SVC + Landcover",
"Temperature SVC",
"Temperature SVC + Landcover",
"Landcover Only (spatial)",
"All Covariates (no SVC)",
"Full SVC Model",
"Full SVC + Interaction")
#Deviance information criteria
DICs = c(barrier.RF.mod$dic$dic,
Linear.mod$dic$dic,
No.LC.sno$dic$dic,
SVC.sno$dic$dic,
No.LC.temp$dic$dic,
SVC.temp$dic$dic,
No.Clim.mod$dic$dic,
Full.noSVC$dic$dic,
SVC.Full$dic$dic,
SVCI.Full$dic$dic)
#Watanabe-Akaike information criteria
WAICs = c(barrier.RF.mod$waic$waic,
Linear.mod$waic$waic,
No.LC.sno$waic$waic,
SVC.sno$waic$waic,
No.LC.temp$waic$waic,
SVC.temp$waic$waic,
No.Clim.mod$waic$waic,
Full.noSVC$waic$waic,
SVC.Full$waic$waic,
SVCI.Full$waic$waic)
lCPOs = c(-mean(log(barrier.RF.mod$cpo$cpo)),
-mean(log(Linear.mod$cpo$cpo)),
-mean(log(No.LC.sno$cpo$cpo)),
-mean(log(SVC.sno$cpo$cpo)),
-mean(log(No.LC.temp$cpo$cpo)),
-mean(log(SVC.temp$cpo$cpo)),
-mean(log(No.Clim.mod$cpo$cpo)),
-mean(log(Full.noSVC$cpo$cpo)),
-mean(log(SVC.Full$cpo$cpo)),
-mean(log(SVCI.Full$cpo$cpo)))
Model_mets = as.data.frame(cbind(Model = Models,
DIC = round(DICs, 2),
WAIC = round(WAICs, 2),
lCPO = round(lCPOs, 2),
Type = Type))
kable(Model_mets) %>%
kable_styling("striped", full_width = F) %>%
row_spec(0, font_size = 20) | Model | DIC | WAIC | lCPO | Type |
|---|---|---|---|---|
| barrier.RF.mod | 458.8 | 453.53 | 0.07 | Random Field |
| Linear.mod | 456.58 | 456.54 | 0.07 | Non-spatial Regression |
| No.LC.sno | 455.21 | 451.17 | 0.07 | Snow SVC |
| SVC.sno | 457.37 | 454.32 | 0.07 | Snow SVC + Landcover |
| No.LC.temp | 454.84 | 451.47 | 0.07 | Temperature SVC |
| SVC.temp | 455.8 | 453.68 | 0.07 | Temperature SVC + Landcover |
| No.Clim.mod | 461.5 | 457.47 | 0.07 | Landcover Only (spatial) |
| Full.noSVC | 454.92 | 452.22 | 0.07 | All Covariates (no SVC) |
| SVC.Full | 456.27 | 452.51 | 0.07 | Full SVC Model |
| SVCI.Full | 453.68 | 449.16 | 0.07 | Full SVC + Interaction |
8.2 Covariates
Cov.fe = SVCI.Full$summary.fixed[,c(1:3,5)]
names(Cov.fe) = c("Mean", "sd", "Q0.025", "Q0.975")
kable(Cov.fe, caption = "Fixed Effects", digits=2) %>%
kable_styling("striped", full_width = F) %>%
row_spec(0, font_size = 20) %>%
column_spec(1, bold = T)| Mean | sd | Q0.025 | Q0.975 | |
|---|---|---|---|---|
| intercept1 | -3.73 | 0.46 | -4.63 | -2.83 |
| mxTempE | -0.85 | 0.32 | -1.47 | -0.23 |
| mSnow5yrE | 0.82 | 0.36 | 0.11 | 1.53 |
| Type.1 | 0.15 | 0.07 | 0.02 | 0.28 |
| Type.2 | 0.12 | 0.22 | -0.31 | 0.56 |
| Type.3 | 0.33 | 0.15 | 0.04 | 0.62 |
| Type.4 | 0.30 | 0.14 | 0.03 | 0.57 |
| Type.6 | 0.15 | 0.23 | -0.29 | 0.59 |
| mxTempE:mSnow5yrE | 0.73 | 0.37 | 0.00 | 1.46 |
#Landcover covariates
kable(EE.names, caption = "Landcover Variables", digits=2) %>%
kable_styling("striped", full_width = F) %>%
row_spec(0, font_size = 20) %>%
column_spec(1, bold = T)| Type | NAME |
|---|---|
| 1 | Evergreen/Deciduous Needleleaf Trees |
| 2 | Evergreen Broadleaf Trees |
| 3 | Deciduous Broadleaf Trees |
| 4 | Mixed/Other Trees |
| 5 | Shrubs |
| 6 | Herbaceous Vegetation |
| 7 | Cultivated and Managed Vegetation |
| 8 | Regularly Flooded Vegetation |
| 9 | Urban/Built-up |
| 10 | Snow/Ice |
| 11 | Barren |
| 12 | Open Water |
8.3 Compare Spatial Fields
Pred.pntsN = Grd.pnts
#Get RF
pLoc = cbind(Pred.pntsN@coords[,1], Pred.pntsN@coords[,2])
Ap = inla.spde.make.A(mesh, loc=pLoc)
Pred.pntsN$Barrier.RF = drop(Ap %*% barrier.RF.mod$summary.random$field.env$mean)
Pred.pntsN$SPDE.RF = drop(Ap %*% spde.RF.mod$summary.random$field.env$mean)
Barrier.RF = rasterize(spTransform(
Pred.pntsN,
CRS(proj4string(Domain.r))),
Domain.r,
"Barrier.RF",
background = NA)
SPDE.RF = rasterize(spTransform(
Pred.pntsN,
CRS(proj4string(Domain.r))),
Domain.r,
"SPDE.RF",
background = NA)
RF.compare = stack(SPDE.RF, Barrier.RF)
names(RF.compare) = c("No.Barrer","Barrier")
Domain2 = subset(Domain, ID != 21 & ID != 48)
rng = seq(-1.1, 1.1, 0.001)
mCols = mCols = c("#A50026", "#D73027", "#F46D43", "#FDAE61", "#E0F3F8", "#ABD9E9", "#74ADD1", "#4575B4","#313695")
cr = colorRampPalette(rev(mCols))(n = 500)
cr = colorRampPalette(cr,
bias = 1, space = "rgb")
levelplot(RF.compare,
layout = c(2,1),
margin = FALSE,
xlab = NULL,
ylab = NULL,
names.attr= c("No Barrier", "Barrier"),
maxpixels = 1e5,
col.regions = cr, at = rng,
colorkey = list(labels=list(at=c(seq(-1,1, 0.5)),
labels=c(paste(seq(-1,1, 0.5))),
fontface='bold', cex=1.5),
labels=list(cex=18),
space = "bottom"),
par.strip.text = list(fontface='bold', cex=1.9),
par.settings = list(axis.line = list(col = "black"),
strip.background = list(col = 'transparent'),
strip.border = list(col = 'transparent')),
scales = list(cex = 1.25)) +
latticeExtra::layer(sp.polygons(Domain2, col = "grey29", fill="transparent", lwd = 0.5)) +
latticeExtra::layer(sp.polygons(States , fill = "transparent", lwd=0.5, col = "grey29")) +
latticeExtra::layer({SpatialPolygonsRescale(layout.north.arrow(),
offset = c(-83.5, 46.3),
scale = c(1.9,1.3))})
8.4 Bayesian Model Averaging
myMods = list(SVC.Full, SVCI.Full)
prior.beta <- function(x, mu = mean(d.mis$bmi, na.rm = TRUE),
sigma = 2*sd(d.mis$bmi, na.rm = TRUE), log = TRUE) {
res <- dnorm(x, mean = mu, sd= sigma, log = log)
if(log) {
return(sum(res))
} else {
return(prod(res))
}
}
mliks = sapply(myMods, function(X){ X$mlik})
probs = exp(mliks - min(mliks))
probs = probs/sum(probs)
#BMA.fx = INLABMA:::fitmatrixBMA(myMods, rep(1/length(myMods), length(myMods)), "summary.fixed")
#BMA.fx
BMA.hyp = INLABMA:::fitmatrixBMA(myMods, probs, "summary.hyperpar")
BMA.hyp## mean sd 0.025quant 0.5quant 0.975quant
## Theta1 for field.env -0.7013524 0.3492621 -1.389836 -0.7017027 -0.01415414
## Theta2 for field.env 3.3135059 0.4766803 2.439882 3.2870524 4.30960649
## Theta1 for field.tem -2.1747651 0.4599868 -3.050548 -2.1860452 -1.24522495
## Theta2 for field.tem 4.0006275 0.6155832 2.916754 3.9506550 5.32108714
## Theta1 for field.sno -2.3279717 0.3625721 -3.043152 -2.3282568 -1.61655982
## Theta2 for field.sno 3.3925361 0.5381687 2.402214 3.3646359 4.51066105
## mode
## Theta1 for field.env -0.7026234
## Theta2 for field.env 3.1919232
## Theta1 for field.tem -2.2263206
## Theta2 for field.tem 3.7683038
## Theta1 for field.sno -2.3289191
## Theta2 for field.sno 3.2641809#marg.hyperpar = INLABMA:::fitmargBMA2(myMods, rep(1/length(myMods), length(myMods)), "marginals.hyperpar")
#marg.hyperpar8.5 Mapped Predictions
ModResult = SVCI.Full
Pred.pntsN = Grd.pnts
#Get RF
pLoc = cbind(Pred.pntsN@coords[,1], Pred.pntsN@coords[,2])
Ap = inla.spde.make.A(mesh, loc=pLoc)
Pred.pntsN$RF = drop(Ap %*% ModResult$summary.random$field.env$mean)
#Intercept
Fixed.df = ModResult$summary.fixed
Pred.pntsN$Sno = Pred.pntsN$s.mSnow5yrE*drop(Ap %*% ModResult$summary.random$field.sno$mean)
Pred.pntsN$Tem = Pred.pntsN$s.mxTempE*drop(Ap %*% ModResult$summary.random$field.tem$mean)
Pred.pntsN$TemEf = drop(Ap %*% ModResult$summary.random$field.tem$mean) + Fixed.df[2,1]
Pred.pntsN$SnoEf = drop(Ap %*% ModResult$summary.random$field.sno$mean) + Fixed.df[3,1]
Pred.pntsN@data = Pred.pntsN@data %>%
mutate(E.LP = RF + Sno + Tem +
Fixed.df[1,1] +
s.mxTempE*Fixed.df[2,1] +
s.mSnow5yrE*Fixed.df[3,1] +
s.Type.1*Fixed.df[4,1] +
s.Type.2*Fixed.df[5,1] +
s.Type.3*Fixed.df[6,1] +
s.Type.4*Fixed.df[7,1] +
s.Type.6*Fixed.df[8,1] +
(s.mxTempE*s.mSnow5yrE)*Fixed.df[9,1],
P.LP = Sno + Tem +
s.mxTempE*Fixed.df[2,1] +
s.mSnow5yrE*Fixed.df[3,1] +
s.Type.1*Fixed.df[4,1] +
s.Type.2*Fixed.df[5,1] +
s.Type.3*Fixed.df[6,1] +
s.Type.4*Fixed.df[7,1] +
s.Type.6*Fixed.df[8,1] +
(s.mxTempE*s.mSnow5yrE)*Fixed.df[9,1],
Est = plogis(E.LP),
Pred = plogis(P.LP))
Test.pred = rasterize(spTransform(
Pred.pntsN,
CRS(proj4string(Domain.r))),
Domain.r,
"Pred",
background = NA)
Test.pred = Test.pred + Domain.r #remove areas outside domain
Test.pred = disaggregate(Test.pred, fact=2)
Test.pred = focal(Test.pred, w=matrix(1, 3, 3), mean) #smooth image
Test.est = rasterize(spTransform(
Pred.pntsN,
CRS(proj4string(Domain.r))),
Domain.r,
"Est",
background = NA)
Test.sno= rasterize(spTransform(
Pred.pntsN,
CRS(proj4string(Domain.r))),
Domain.r,
"SnoEf",
background = NA)
Test.tem= rasterize(spTransform(
Pred.pntsN,
CRS(proj4string(Domain.r))),
Domain.r,
"TemEf",
background = NA)
Test.rf = rasterize(spTransform(
Pred.pntsN,
CRS(proj4string(Domain.r))),
Domain.r,
"RF",
background = NA)rng = seq(0, 1, 0.01)
mCols = brewer.pal(9, "YlOrBr")
cr = colorRampPalette(mCols)(n = 500)
cr = colorRampPalette(cr,
bias = 0.75, space = "rgb")
Domain2 = subset(Domain, ID != 21 & ID != 48)
rng = seq(0, 1, 0.01)
mCols = rev(inferno(10))
#mCols[1] = "lightgray"
cr = colorRampPalette(c(mCols),
bias = 0.3)
levelplot(Test.pred, #Trial.comp,
#layout = c(2,1),
margin = FALSE,
xlab = NULL,
ylab = NULL,
#main = "Occurrence Probabilty",
#names.attr= c("A", "B", "C"),
maxpixels = 1e6,
col.regions = cr, at = rng,
colorkey = list(labels=list(at=c(0.02, 0.25, 0.50, 0.75, 0.98),
labels=c("0.00", "0.25", "0.50", "0.75", "1.00"),
fontface='bold', cex=1.5),
labels=list(cex=18),
space = "bottom"),
par.strip.text = list(fontface='bold', cex=1.9),
par.settings = list(axis.line = list(col = "black"),
strip.background = list(col = 'transparent'),
strip.border = list(col = 'transparent')),
scales = list(cex = 1.25)) +
latticeExtra::layer(sp.polygons(Domain2, col = "grey26", fill="lightgray", lwd = 1)) +
latticeExtra::layer(sp.polygons(LakesLL , fill = "lightblue", alpha = 1, lwd=2, col = "grey29")) +
latticeExtra::layer({SpatialPolygonsRescale(layout.north.arrow(),
offset = c(-83.5, 46.3),
scale = c(1.9,1.3))})
#Save grids
#writeRaster(Test.pred, "C:/Users/humph173/Documents/Michigan_State/Sean/Dec1419_figs/pred_map.tif")8.6 SVC Covariates
Cov.stk = stack(Test.sno, Test.tem)
names(Cov.stk) = c("Snow", "Temp")
Cov.stk = Cov.stk+Domain.r
Domain2 = subset(Domain, ID != 21 & ID != 48)
rng = seq(0.80, 0.83, 0.0001)
mCols = rev(viridis(100))
cr = colorRampPalette(c(mCols),
bias = 0.5)
levelplot(Cov.stk[[1]], #Trial.comp,
#layout = c(2,1),
margin = FALSE,
xlab = NULL,
ylab = NULL,
main = "Snow Weeks Effect Size \n(mean = 0.82)",
maxpixels = 1e5,
col.regions = cr, at = rng,
colorkey = list(labels=list(at=c(seq(0.80, 0.83, 0.01)),
labels=c(paste(round(seq(0.80, 0.83, 0.01),2))),
fontface='bold', cex=1.5),
labels=list(cex=18),
space = "bottom"),
par.strip.text = list(fontface='bold', cex=1.9),
par.settings = list(axis.line = list(col = "black"),
strip.background = list(col = 'transparent'),
strip.border = list(col = 'transparent')),
scales = list(cex = 1.25)) +
latticeExtra::layer(sp.polygons(Domain2, col = "grey26", fill="lightgray", lwd = 1)) +
latticeExtra::layer(sp.polygons(LakesLL , fill = "lightblue", alpha = 1, lwd=2, col = "grey29")) +
latticeExtra::layer({SpatialPolygonsRescale(layout.north.arrow(),
offset = c(-83.5, 46.3),
scale = c(1.9,1.3))})
Cov.stk = stack(Test.sno, Test.tem)
names(Cov.stk) = c("Snow", "Temp")
Cov.stk = Cov.stk+Domain.r
Domain2 = subset(Domain, ID != 21 & ID != 48)
rng = seq(-0.84, -0.87, -0.0001)
mCols = rev(viridis(100))
cr = colorRampPalette(c(mCols),
bias = 0.5)
levelplot(Cov.stk[[2]], #Trial.comp,
#layout = c(2,1),
margin = FALSE,
xlab = NULL,
ylab = NULL,
main = "Temperature Effect Size \n(mean = -0.85)",
#names.attr= c("Snow Weeks (mean = 0.55)"),
maxpixels = 1e5,
col.regions = cr, at = rng,
colorkey = list(labels=list(at=c(seq(-0.84, -0.87, -0.01)),
labels=c(paste(seq(-0.84, -0.87, -0.01))),
fontface='bold', cex=1.5),
labels=list(cex=18),
space = "bottom"),
par.strip.text = list(fontface='bold', cex=1.9),
par.settings = list(axis.line = list(col = "black"),
strip.background = list(col = 'transparent'),
strip.border = list(col = 'transparent')),
scales = list(cex = 1.25)) +
latticeExtra::layer(sp.polygons(Domain2, col = "grey26", fill="lightgray", lwd = 1)) +
latticeExtra::layer(sp.polygons(LakesLL , fill = "lightblue", alpha = 1, lwd=2, col = "grey29")) +
latticeExtra::layer({SpatialPolygonsRescale(layout.north.arrow(),
offset = c(-83.5, 46.3),
scale = c(1.9,1.3))})
8.7 SVC Correlation
cor.test(values(Cov.stk[[1]]), values(Cov.stk[[2]]))##
## Pearson's product-moment correlation
##
## data: values(Cov.stk[[1]]) and values(Cov.stk[[2]])
## t = 497.02, df = 124350, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.8137119 0.8174340
## sample estimates:
## cor
## 0.8155814library(spatialEco)##
## Attaching package: 'spatialEco'## The following object is masked from 'package:data.table':
##
## shift## The following object is masked from 'package:spData':
##
## elev## The following object is masked from 'package:raster':
##
## shiftr.cor = rasterCorrelation(Cov.stk[[1]], Cov.stk[[2]], s = 3, type = "pearson")
rng = seq(-1, 1, 0.001)
mCols = mCols = c("#A50026", "#D73027", "#F46D43", "#FDAE61", "#E0F3F8", "#ABD9E9", "#74ADD1", "#4575B4","#313695")
cr = colorRampPalette(rev(mCols))(n = 500)
cr = colorRampPalette(cr,
bias = 1, space = "rgb")
levelplot(r.cor, #Trial.comp,
margin = FALSE,
xlab = NULL,
ylab = NULL,
main = "Snow and Temp SVC Correlation",
#names.attr= c("Snow Weeks (mean = 0.55)"),
maxpixels = 1e5,
col.regions = cr, at = rng,
colorkey = list(labels=list(at=c(seq(-1, 1, 0.25)),
labels=c(paste(paste(seq(-1, 1, 0.25)))),
fontface='bold', cex=1.5),
labels=list(cex=18),
space = "bottom"),
par.strip.text = list(fontface='bold', cex=1.9),
par.settings = list(axis.line = list(col = "black"),
strip.background = list(col = 'transparent'),
strip.border = list(col = 'transparent')),
scales = list(cex = 1.25)) +
latticeExtra::layer(sp.polygons(Domain2, col = "grey26", fill="lightgray", lwd = 1)) +
latticeExtra::layer(sp.polygons(LakesLL , fill = "lightblue", alpha = 1, lwd=2, col = "grey29")) +
latticeExtra::layer({SpatialPolygonsRescale(layout.north.arrow(),
offset = c(-83.5, 46.3),
scale = c(1.9,1.3))})
8.8 Compare Region Coefficients
ModResult = SVCI.Full
Pred.pntsN = subset(OneTier.mod, Spp == "Hare")
#Get LOcations
pLoc = cbind(Pred.pntsN@coords[,1], Pred.pntsN@coords[,2])
Ap = inla.spde.make.A(mesh, loc=pLoc)
Pred.pntsN$Snow = drop(Ap %*% ModResult$summary.random$field.sno$mean) + ModResult$summary.fixed[3,1]
Pred.pntsN$Temp = drop(Ap %*% ModResult$summary.random$field.tem$mean) + ModResult$summary.fixed[2,1]
plot.df = Pred.pntsN@data
plot.df$SiteX = 1:nrow(plot.df)
plot.df$StateUP[plot.df$StateUP == "Mich.UP"] = "Michigan (Upper)"
plot.df$StateUP[plot.df$StateUP == "Michigan"] = "Michigan (Lower)"
plot.df$StateUP[plot.df$StateUP == "Wisconsin"] = "Wisconsin"8.9 Temperature
sd(HR.df$mxTempE)## [1] 1.639839mean(HR.df$mxTempE)## [1] 17.61804plot.df$Temp.obs = plogis(c(-1*plot.df$Temp))
ggplot(plot.df, aes(StateUP, -1*Temp.obs, fill = StateUP)) +
geom_dotplot(binaxis='y', stackdir='center',
stackratio=0.2, dotsize=1,
position = position_jitter(0, 0.0005), col="gray90") +
geom_violin(trim = FALSE, fill = "gray95", col = "gray30", alpha=0.6) +
scale_fill_manual(name="Survey Sites", values=c("#999999", "#E69F00", "#56B4E9")) +
scale_y_continuous(breaks = seq(-0.704, -0.698, 0.001),
labels = scales::percent_format(accuracy = 0.1),
limits = c(-0.704, -0.698)) +
theme_classic() +
xlab(" ") +
ylab("Chance of Occurrence") +
#ylim(-0.865, -0.840) +
theme(plot.title = element_text(hjust = 0.5),
legend.position = "bottom",
plot.margin = unit(c(1,1,4,1),"cm"),
legend.title=element_text(size=22),
legend.text=element_text(size=18),
#panel.grid.major = element_blank(),
#panel.grid.minor = element_blank(),
strip.background = element_blank(),
#panel.border = element_rect(colour = "black"),
title = element_text(face="bold", size=22, hjust=0.5),
strip.text = element_text(face="bold", size = 22),
axis.title.y = element_text(face="bold", size = 24),
axis.ticks.x = element_blank(),
axis.text.x = element_blank(),
axis.text.y = element_text(face="bold", size=20),
axis.title.x = element_text(face="bold", size = 22))## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.
8.10 Snow Cover
sd(HR.df$mSnow5yrE)## [1] 4.822815mean(HR.df$mSnow5yrE)## [1] 13.14616plot.df$Snow.obs = plogis(plot.df$Snow)
ggplot(plot.df, aes(StateUP, Snow.obs, fill = StateUP)) +
geom_dotplot(binaxis='y', stackdir='center',
stackratio=0.3, dotsize=1,
position = position_jitter(0, 0.0005), col="gray90") +
geom_violin(trim = FALSE, fill = "gray95", col = "gray30", alpha=0.6) +
scale_fill_manual(name="Survey Sites", values=c("#999999", "#E69F00", "#56B4E9")) +
scale_y_continuous(breaks = seq(0.691, 0.697, 0.001),
labels = scales::percent_format(accuracy = 0.1),
limits = c(0.691, 0.697)) +
theme_classic() +
xlab(" ") +
ylab("Chance of Occurrence") +
#ylim(0.691, 0.697) +
theme(plot.title = element_text(hjust = 0.5),
legend.position = "bottom",
legend.title=element_text(size=22),
legend.text=element_text(size=18),
plot.margin = unit(c(1,1,3,1),"cm"),
#panel.grid.major = element_blank(),
#panel.grid.minor = element_blank(),
strip.background = element_blank(),
#panel.border = element_rect(colour = "black"),
title = element_text(face="bold", size=22, hjust=0.5),
strip.text = element_text(face="bold", size = 22),
axis.title.y = element_text(face="bold", size = 24),
axis.ticks.x = element_blank(),
axis.text.x = element_blank(),
axis.text.y = element_text(face="bold", size=20),
axis.title.x = element_text(face="bold", size = 22))## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.
9 Validation
Full SVC Model (Out of Sample)
ModResult = SVCI.Full
Pred.pntsN = Validation.set2
#Get RF
pLoc = cbind(Pred.pntsN@coords[,1], Pred.pntsN@coords[,2])
Ap = inla.spde.make.A(mesh, loc=pLoc)
Pred.pntsN$RF = drop(Ap %*% ModResult$summary.random$field.env$mean)
#Intercept
Fixed.df = ModResult$summary.fixed
Pred.pntsN$Sno = drop(Ap %*% ModResult$summary.random$field.sno$mean) + ModResult$summary.fixed[3,1]
Pred.pntsN$Tem = drop(Ap %*% ModResult$summary.random$field.tem$mean) + ModResult$summary.fixed[2,1]
Pred.pntsN@data = Pred.pntsN@data %>%
mutate(E.LP = RF + Sno + Tem +
Fixed.df[1,1] +
s.mxTempE*Fixed.df[2,1] +
s.mSnow5yrE*Fixed.df[3,1] +
s.Type.1*Fixed.df[4,1] +
s.Type.2*Fixed.df[5,1] +
s.Type.3*Fixed.df[6,1] +
s.Type.4*Fixed.df[7,1] +
s.Type.6*Fixed.df[8,1] +
(s.mxTempE*s.mSnow5yrE)*Fixed.df[9,1],
P.LP = Sno + Tem +
s.mxTempE*Fixed.df[2,1] +
s.mSnow5yrE*Fixed.df[3,1] +
s.Type.1*Fixed.df[4,1] +
s.Type.2*Fixed.df[5,1] +
s.Type.3*Fixed.df[6,1] +
s.Type.4*Fixed.df[7,1] +
s.Type.6*Fixed.df[8,1] +
(s.mxTempE*s.mSnow5yrE)*Fixed.df[9,1],
Est = plogis(E.LP),
Pred = plogis(P.LP))Reg2.PA = as.data.frame(cbind(1:dim(Pred.pntsN@data)[1], Pred.pntsN$OBS, Pred.pntsN$Pred))
names(Reg2.PA) = c("ID", "OBS", "Pred")
Thresh = optimal.thresholds(Reg2.PA, opt.methods = c("MaxSens+Spec", "PredPrev=Obs"))
Thresh## Method Pred
## 1 MaxSens+Spec 0.77
## 2 PredPrev=Obs 0.81Reg2.PA.out1 = presence.absence.accuracy(Reg2.PA, threshold = Thresh[1,2])
Reg2.PA.out2 = presence.absence.accuracy(Reg2.PA, threshold = 0.5)
Reg2.PA.out = rbind(Reg2.PA.out1, Reg2.PA.out2)
#Calculate TSS
Reg2.PA.out$TSS = Reg2.PA.out$sensitivity + Reg2.PA.out$specificity - 1
Reg2.PA.out## model threshold PCC sensitivity specificity Kappa AUC
## 1 Pred 0.77 0.8580866 0.9688109 0.8209735 0.6766628 0.8909207
## 2 Pred 0.50 0.7080989 0.9980507 0.6109115 0.4396149 0.8909207
## PCC.sd sensitivity.sd specificity.sd Kappa.sd AUC.sd TSS
## 1 0.005459188 0.005429488 0.006930467 0.01189092 0.005121451 0.7897845
## 2 0.007112392 0.001377703 0.008813586 0.01132580 0.005121451 0.6089621Validation.tab = Reg2.PA.out %>% select(threshold, PCC, sensitivity, specificity, AUC, TSS)
kable(Validation.tab, caption = "Out of Sample (SVC)", digits=2) %>%
kable_styling("striped", full_width = F) %>%
row_spec(0, font_size = 20) %>%
column_spec(1, bold = T)| threshold | PCC | sensitivity | specificity | AUC | TSS |
|---|---|---|---|---|---|
| 0.77 | 0.86 | 0.97 | 0.82 | 0.89 | 0.79 |
| 0.50 | 0.71 | 1.00 | 0.61 | 0.89 | 0.61 |
9.1 One Sample per Cell
Thinned Out of Sample (SVC)
Thin.grid = mxTemp
Thin.grid = projectRaster(Thin.grid, crs = proj4string(Pred.pntsN))
#Select one set of coordinates per cell
set.seed(1976)
Grid.Samp = as.data.frame(gridSample(Pred.pntsN@data[,c("Long", "Lat")], Thin.grid, n = 1))
names(Grid.Samp) = c("Long", "Lat")
#Sample size
dim(Grid.Samp)[1] ## [1] 3178#Match to data
Grid.Samp$Sample = "Keep"
temp.sample = right_join(Pred.pntsN@data, Grid.Samp, by = c("Long", "Lat"))
#Scores
Reg2.PA = as.data.frame(cbind(1:dim(temp.sample)[1], temp.sample$OBS, temp.sample$Pred))
names(Reg2.PA) = c("ID", "OBS", "Pred")
Thresh = optimal.thresholds(Reg2.PA, opt.methods = c("MaxSens+Spec", "PredPrev=Obs"))
Thresh## Method Pred
## 1 MaxSens+Spec 0.64
## 2 PredPrev=Obs 0.83Reg2.PA.out1 = presence.absence.accuracy(Reg2.PA, threshold = Thresh[1,2])
Reg2.PA.out2 = presence.absence.accuracy(Reg2.PA, threshold = 0.5)
Reg2.PA.out = rbind(Reg2.PA.out1, Reg2.PA.out2)
#Calculate TSS
Reg2.PA.out$TSS = Reg2.PA.out$sensitivity + Reg2.PA.out$specificity - 1
Reg2.PA.out## model threshold PCC sensitivity specificity Kappa AUC
## 1 Pred 0.64 0.7193203 0.9692308 0.7086614 0.1593171 0.883028
## 2 Pred 0.50 0.6286973 0.9846154 0.6135171 0.1122165 0.883028
## PCC.sd sensitivity.sd specificity.sd Kappa.sd AUC.sd TSS
## 1 0.007971830 0.01520467 0.008231560 0.013076840 0.008300148 0.6778922
## 2 0.008571883 0.01083632 0.008821496 0.009519187 0.008300148 0.5981324Validation.tab = Reg2.PA.out %>% select(threshold, PCC, sensitivity, specificity, AUC, TSS)
kable(Validation.tab, caption = "Thinned Out of Sample (SVC)", digits=2) %>%
kable_styling("striped", full_width = F) %>%
row_spec(0, font_size = 20) %>%
column_spec(1, bold = T)| threshold | PCC | sensitivity | specificity | AUC | TSS |
|---|---|---|---|---|---|
| 0.64 | 0.72 | 0.97 | 0.71 | 0.88 | 0.68 |
| 0.50 | 0.63 | 0.98 | 0.61 | 0.88 | 0.60 |
No SVC (Out of Sample)
ModResult = Full.noSVC
Pred.pntsN = Validation.set2
#Get RF
pLoc = cbind(Pred.pntsN@coords[,1], Pred.pntsN@coords[,2])
Ap = inla.spde.make.A(mesh, loc=pLoc)
Pred.pntsN$RF = drop(Ap %*% ModResult$summary.random$field.env$mean)
#Intercept
Fixed.df = ModResult$summary.fixed
#Pred.pntsN$Sno = drop(Ap %*% ModResult$summary.random$field.sno$mean) + ModResult$summary.fixed[3,1]
#Pred.pntsN$Tem = drop(Ap %*% ModResult$summary.random$field.tem$mean) + ModResult$summary.fixed[2,1]
Pred.pntsN@data = Pred.pntsN@data %>%
mutate(E.LP = RF + #Sno + Tem +
Fixed.df[1,1] +
s.mxTempE*Fixed.df[2,1] +
s.mSnow5yrE*Fixed.df[3,1] +
s.Type.1*Fixed.df[4,1] +
s.Type.2*Fixed.df[5,1] +
s.Type.3*Fixed.df[6,1] +
s.Type.4*Fixed.df[7,1] +
s.Type.6*Fixed.df[8,1],
P.LP = #Sno + Tem +
s.mxTempE*Fixed.df[2,1] +
s.mSnow5yrE*Fixed.df[3,1] +
s.Type.1*Fixed.df[4,1] +
s.Type.2*Fixed.df[5,1] +
s.Type.3*Fixed.df[6,1] +
s.Type.4*Fixed.df[7,1] +
s.Type.6*Fixed.df[8,1],
Est = plogis(E.LP),
Pred = plogis(P.LP))Reg2.PA = as.data.frame(cbind(1:dim(Pred.pntsN@data)[1], Pred.pntsN$OBS, Pred.pntsN$Pred))
names(Reg2.PA) = c("ID", "OBS", "Pred")
Thresh = optimal.thresholds(Reg2.PA, opt.methods = c("MaxSens+Spec", "PredPrev=Obs"))
Thresh## Method Pred
## 1 MaxSens+Spec 0.77
## 2 PredPrev=Obs 0.79Reg2.PA.out1 = presence.absence.accuracy(Reg2.PA, threshold = Thresh[1,2])
Reg2.PA.out2 = presence.absence.accuracy(Reg2.PA, threshold = 0.5)
Reg2.PA.out = rbind(Reg2.PA.out1, Reg2.PA.out2)
#Calculate TSS
Reg2.PA.out$TSS = Reg2.PA.out$sensitivity + Reg2.PA.out$specificity - 1
Reg2.PA.out## model threshold PCC sensitivity specificity Kappa AUC
## 1 Pred 0.77 0.8705652 0.9541910 0.8425351 0.6983110 0.8797982
## 2 Pred 0.50 0.7320773 0.9990253 0.6426005 0.4738342 0.8797982
## PCC.sd sensitivity.sd specificity.sd Kappa.sd AUC.sd TSS
## 1 0.005251423 0.0065302676 0.006584536 0.01182564 0.005404061 0.7967262
## 2 0.006928417 0.0009746589 0.008663369 0.01157761 0.005404061 0.6416258Validation.tab = Reg2.PA.out %>% select(threshold, PCC, sensitivity, specificity, AUC, TSS)
kable(Validation.tab, caption = "Out of Sample (No SVC)", digits=2) %>%
kable_styling("striped", full_width = F) %>%
row_spec(0, font_size = 20) %>%
column_spec(1, bold = T)| threshold | PCC | sensitivity | specificity | AUC | TSS |
|---|---|---|---|---|---|
| 0.77 | 0.87 | 0.95 | 0.84 | 0.88 | 0.80 |
| 0.50 | 0.73 | 1.00 | 0.64 | 0.88 | 0.64 |
9.2 One Sample per Cell
Thinned Out of Sample (No SVC)
#Select one set of coordinates per cell
set.seed(1976)
Grid.Samp = as.data.frame(gridSample(Pred.pntsN@data[,c("Long", "Lat")], Thin.grid, n = 1))
names(Grid.Samp) = c("Long", "Lat")
#Sample size
dim(Grid.Samp)[1] ## [1] 3178#Match to data
Grid.Samp$Sample = "Keep"
temp.sample = right_join(Pred.pntsN@data, Grid.Samp, by = c("Long", "Lat"))
#Scores
Reg2.PA = as.data.frame(cbind(1:dim(temp.sample)[1], temp.sample$OBS, temp.sample$Pred))
names(Reg2.PA) = c("ID", "OBS", "Pred")
Thresh = optimal.thresholds(Reg2.PA, opt.methods = c("MaxSens+Spec", "PredPrev=Obs"))
Thresh## Method Pred
## 1 MaxSens+Spec 0.59
## 2 PredPrev=Obs 0.86Reg2.PA.out1 = presence.absence.accuracy(Reg2.PA, threshold = Thresh[1,2])
Reg2.PA.out2 = presence.absence.accuracy(Reg2.PA, threshold = 0.5)
Reg2.PA.out = rbind(Reg2.PA.out1, Reg2.PA.out2)
#Calculate TSS
Reg2.PA.out$TSS = Reg2.PA.out$sensitivity + Reg2.PA.out$specificity - 1
Reg2.PA.out## model threshold PCC sensitivity specificity Kappa AUC
## 1 Pred 0.59 0.7240403 0.9692308 0.7135827 0.1625841 0.8757167
## 2 Pred 0.50 0.6595343 0.9923077 0.6453412 0.1281273 0.8757167
## PCC.sd sensitivity.sd specificity.sd Kappa.sd AUC.sd TSS
## 1 0.007930409 0.015204672 0.008190031 0.01329871 0.008406231 0.6828134
## 2 0.008407111 0.007692308 0.008666901 0.01058299 0.008406231 0.6376489Validation.tab = Reg2.PA.out %>% select(threshold, PCC, sensitivity, specificity, AUC, TSS)
kable(Validation.tab, caption = "Thinned Out of Sample (No SVC)", digits=2) %>%
kable_styling("striped", full_width = F) %>%
row_spec(0, font_size = 20) %>%
column_spec(1, bold = T)| threshold | PCC | sensitivity | specificity | AUC | TSS |
|---|---|---|---|---|---|
| 0.59 | 0.72 | 0.97 | 0.71 | 0.88 | 0.68 |
| 0.50 | 0.66 | 0.99 | 0.65 | 0.88 | 0.64 |
Full SVC Model (In Sample)
ModResult = SVCI.Full
Pred.pntsN = subset(OneTier.mod, Spp == "Hare")
#Get RF
pLoc = cbind(Pred.pntsN@coords[,1], Pred.pntsN@coords[,2])
Ap = inla.spde.make.A(mesh, loc=pLoc)
Pred.pntsN$RF = drop(Ap %*% ModResult$summary.random$field.env$mean)
#Intercept
Fixed.df = ModResult$summary.fixed
Pred.pntsN$Sno = drop(Ap %*% ModResult$summary.random$field.sno$mean) + ModResult$summary.fixed[3,1]
Pred.pntsN$Tem = drop(Ap %*% ModResult$summary.random$field.tem$mean) + ModResult$summary.fixed[2,1]
Pred.pntsN@data = Pred.pntsN@data %>%
mutate(E.LP = RF + Sno + Tem +
Fixed.df[1,1] +
s.mxTempE*Fixed.df[2,1] +
s.mSnow5yrE*Fixed.df[3,1] +
s.Type.1*Fixed.df[4,1] +
s.Type.2*Fixed.df[5,1] +
s.Type.3*Fixed.df[6,1] +
s.Type.4*Fixed.df[7,1] +
s.Type.6*Fixed.df[8,1] +
(s.mxTempE*s.mSnow5yrE)*Fixed.df[9,1],
P.LP = Sno + Tem +
s.mxTempE*Fixed.df[2,1] +
s.mSnow5yrE*Fixed.df[3,1] +
s.Type.1*Fixed.df[4,1] +
s.Type.2*Fixed.df[5,1] +
s.Type.3*Fixed.df[6,1] +
s.Type.4*Fixed.df[7,1] +
s.Type.6*Fixed.df[8,1] +
(s.mxTempE*s.mSnow5yrE)*Fixed.df[9,1],
Est = plogis(E.LP),
Pred = plogis(P.LP))Reg2.PA = as.data.frame(cbind(1:dim(Pred.pntsN@data)[1], Pred.pntsN$OBS, Pred.pntsN$Pred))
names(Reg2.PA) = c("ID", "OBS", "Pred")
Thresh = optimal.thresholds(Reg2.PA, opt.methods = c("MaxSens+Spec", "PredPrev=Obs"))
Thresh## Method Pred
## 1 MaxSens+Spec 0.71
## 2 PredPrev=Obs 0.73Reg2.PA.out1 = presence.absence.accuracy(Reg2.PA, threshold = Thresh[1,2])
Reg2.PA.out2 = presence.absence.accuracy(Reg2.PA, threshold = 0.5)
Reg2.PA.out = rbind(Reg2.PA.out1, Reg2.PA.out2)
#Calculate TSS
Reg2.PA.out$TSS = Reg2.PA.out$sensitivity + Reg2.PA.out$specificity - 1
Reg2.PA.out## model threshold PCC sensitivity specificity Kappa AUC
## 1 Pred 0.71 0.71875 0.7209302 0.7172775 0.4285034 0.7805917
## 2 Pred 0.50 0.51250 0.9767442 0.1989529 0.1478029 0.7805917
## PCC.sd sensitivity.sd specificity.sd Kappa.sd AUC.sd TSS
## 1 0.02517328 0.03964588 0.03266984 0.05056029 0.02613074 0.4382077
## 2 0.02798588 0.01332144 0.02896192 0.02864680 0.02613074 0.1756971Validation.tab = Reg2.PA.out %>% select(threshold, PCC, sensitivity, specificity, AUC, TSS)
kable(Validation.tab, caption = "In Sample (SVC)", digits=2) %>%
kable_styling("striped", full_width = F) %>%
row_spec(0, font_size = 20) %>%
column_spec(1, bold = T)| threshold | PCC | sensitivity | specificity | AUC | TSS |
|---|---|---|---|---|---|
| 0.71 | 0.72 | 0.72 | 0.72 | 0.78 | 0.44 |
| 0.50 | 0.51 | 0.98 | 0.20 | 0.78 | 0.18 |
9.3 One Sample per Cell
Thinned In Sample (SVC)
#Select one set of coordinates per cell
set.seed(1976)
Grid.Samp = as.data.frame(gridSample(Pred.pntsN@data[,c("Long", "Lat")], Thin.grid, n = 1))
names(Grid.Samp) = c("Long", "Lat")
#Sample size
dim(Grid.Samp)[1] ## [1] 314#Match to data
Grid.Samp$Sample = "Keep"
temp.sample = right_join(Pred.pntsN@data, Grid.Samp, by = c("Long", "Lat"))
#Scores
Reg2.PA = as.data.frame(cbind(1:dim(temp.sample)[1], temp.sample$OBS, temp.sample$Pred))
names(Reg2.PA) = c("ID", "OBS", "Pred")
Thresh = optimal.thresholds(Reg2.PA, opt.methods = c("MaxSens+Spec", "PredPrev=Obs"))
Thresh## Method Pred
## 1 MaxSens+Spec 0.710
## 2 PredPrev=Obs 0.725Reg2.PA.out1 = presence.absence.accuracy(Reg2.PA, threshold = Thresh[1,2])
Reg2.PA.out2 = presence.absence.accuracy(Reg2.PA, threshold = 0.5)
Reg2.PA.out = rbind(Reg2.PA.out1, Reg2.PA.out2)
#Calculate TSS
Reg2.PA.out$TSS = Reg2.PA.out$sensitivity + Reg2.PA.out$specificity - 1
Reg2.PA.out## model threshold PCC sensitivity specificity Kappa AUC
## 1 Pred 0.71 0.7261146 0.7244094 0.7272727 0.4427339 0.7849173
## 2 Pred 0.50 0.5127389 0.9763780 0.1978610 0.1469512 0.7849173
## PCC.sd sensitivity.sd specificity.sd Kappa.sd AUC.sd TSS
## 1 0.02520663 0.03980512 0.03265551 0.05074806 0.02624454 0.4516822
## 2 0.02825250 0.01352952 0.02921113 0.02895548 0.02624454 0.1742389Validation.tab = Reg2.PA.out %>% select(threshold, PCC, sensitivity, specificity, AUC, TSS)
kable(Validation.tab, caption = "Thinned In Sample (SVC)", digits=2) %>%
kable_styling("striped", full_width = F) %>%
row_spec(0, font_size = 20) %>%
column_spec(1, bold = T)| threshold | PCC | sensitivity | specificity | AUC | TSS |
|---|---|---|---|---|---|
| 0.71 | 0.73 | 0.72 | 0.73 | 0.78 | 0.45 |
| 0.50 | 0.51 | 0.98 | 0.20 | 0.78 | 0.17 |
No SVC (In Sample)
ModResult = Full.noSVC
Pred.pntsN = subset(OneTier.mod, Spp == "Hare")
#Get RF
pLoc = cbind(Pred.pntsN@coords[,1], Pred.pntsN@coords[,2])
Ap = inla.spde.make.A(mesh, loc=pLoc)
Pred.pntsN$RF = drop(Ap %*% ModResult$summary.random$field.env$mean)
#Intercept
Fixed.df = ModResult$summary.fixed
#Pred.pntsN$Sno = drop(Ap %*% ModResult$summary.random$field.sno$mean) + ModResult$summary.fixed[3,1]
#Pred.pntsN$Tem = drop(Ap %*% ModResult$summary.random$field.tem$mean) + ModResult$summary.fixed[2,1]
Pred.pntsN@data = Pred.pntsN@data %>%
mutate(E.LP = RF + #Sno + Tem +
Fixed.df[1,1] +
s.mxTempE*Fixed.df[2,1] +
s.mSnow5yrE*Fixed.df[3,1] +
s.Type.1*Fixed.df[4,1] +
s.Type.2*Fixed.df[5,1] +
s.Type.3*Fixed.df[6,1] +
s.Type.4*Fixed.df[7,1] +
s.Type.6*Fixed.df[8,1],
P.LP = #Sno + Tem +
s.mxTempE*Fixed.df[2,1] +
s.mSnow5yrE*Fixed.df[3,1] +
s.Type.1*Fixed.df[4,1] +
s.Type.2*Fixed.df[5,1] +
s.Type.3*Fixed.df[6,1] +
s.Type.4*Fixed.df[7,1] +
s.Type.6*Fixed.df[8,1],
Est = plogis(E.LP),
Pred = plogis(P.LP))Reg2.PA = as.data.frame(cbind(1:dim(Pred.pntsN@data)[1], Pred.pntsN$OBS, Pred.pntsN$Pred))
names(Reg2.PA) = c("ID", "OBS", "Pred")
Thresh = optimal.thresholds(Reg2.PA, opt.methods = c("MaxSens+Spec", "PredPrev=Obs"))
Thresh## Method Pred
## 1 MaxSens+Spec 0.68
## 2 PredPrev=Obs 0.68Reg2.PA.out1 = presence.absence.accuracy(Reg2.PA, threshold = Thresh[1,2])
Reg2.PA.out2 = presence.absence.accuracy(Reg2.PA, threshold = 0.5)
Reg2.PA.out = rbind(Reg2.PA.out1, Reg2.PA.out2)
#Calculate TSS
Reg2.PA.out$TSS = Reg2.PA.out$sensitivity + Reg2.PA.out$specificity - 1
Reg2.PA.out## model threshold PCC sensitivity specificity Kappa AUC
## 1 Pred 0.68 0.734375 0.6744186 0.7748691 0.4487231 0.7794959
## 2 Pred 0.50 0.537500 0.9457364 0.2617801 0.1775779 0.7794959
## PCC.sd sensitivity.sd specificity.sd Kappa.sd AUC.sd TSS
## 1 0.02472852 0.04141804 0.03030086 0.05085353 0.02628688 0.4492877
## 2 0.02791578 0.02002323 0.03189219 0.03411622 0.02628688 0.2075165Validation.tab = Reg2.PA.out %>% select(threshold, PCC, sensitivity, specificity, AUC, TSS)
kable(Validation.tab, caption = "In Sample (No SVC)", digits=2) %>%
kable_styling("striped", full_width = F) %>%
row_spec(0, font_size = 20) %>%
column_spec(1, bold = T)| threshold | PCC | sensitivity | specificity | AUC | TSS |
|---|---|---|---|---|---|
| 0.68 | 0.73 | 0.67 | 0.77 | 0.78 | 0.45 |
| 0.50 | 0.54 | 0.95 | 0.26 | 0.78 | 0.21 |
9.4 One Sample per Cell
Thinned In Sample (No SVC)
#Select one set of coordinates per cell
set.seed(1976)
Grid.Samp = as.data.frame(gridSample(Pred.pntsN@data[,c("Long", "Lat")], Thin.grid, n = 1))
names(Grid.Samp) = c("Long", "Lat")
#Sample size
dim(Grid.Samp)[1] ## [1] 314#Match to data
Grid.Samp$Sample = "Keep"
temp.sample = right_join(Pred.pntsN@data, Grid.Samp, by = c("Long", "Lat"))
#Scores
Reg2.PA = as.data.frame(cbind(1:dim(temp.sample)[1], temp.sample$OBS, temp.sample$Pred))
names(Reg2.PA) = c("ID", "OBS", "Pred")
Thresh = optimal.thresholds(Reg2.PA, opt.methods = c("MaxSens+Spec", "PredPrev=Obs"))
Thresh## Method Pred
## 1 MaxSens+Spec 0.68
## 2 PredPrev=Obs 0.68Reg2.PA.out1 = presence.absence.accuracy(Reg2.PA, threshold = Thresh[1,2])
Reg2.PA.out2 = presence.absence.accuracy(Reg2.PA, threshold = 0.5)
Reg2.PA.out = rbind(Reg2.PA.out1, Reg2.PA.out2)
#Calculate TSS
Reg2.PA.out$TSS = Reg2.PA.out$sensitivity + Reg2.PA.out$specificity - 1
Reg2.PA.out## model threshold PCC sensitivity specificity Kappa AUC
## 1 Pred 0.68 0.7420382 0.6771654 0.7860963 0.4638475 0.7844962
## 2 Pred 0.50 0.5382166 0.9448819 0.2620321 0.1775353 0.7844962
## PCC.sd sensitivity.sd specificity.sd Kappa.sd AUC.sd TSS
## 1 0.02472969 0.04165356 0.03006703 0.05093924 0.02638573 0.4632616
## 2 0.02817900 0.02033062 0.03224330 0.03458085 0.02638573 0.2069140Validation.tab = Reg2.PA.out %>% select(threshold, PCC, sensitivity, specificity, AUC, TSS)
kable(Validation.tab, caption = "Thinned In Sample (No SVC)", digits=2) %>%
kable_styling("striped", full_width = F) %>%
row_spec(0, font_size = 20) %>%
column_spec(1, bold = T)| threshold | PCC | sensitivity | specificity | AUC | TSS |
|---|---|---|---|---|---|
| 0.68 | 0.74 | 0.68 | 0.79 | 0.78 | 0.46 |
| 0.50 | 0.54 | 0.94 | 0.26 | 0.78 | 0.21 |
10 Variable Plots
10.1 Snow Weeks
rng = seq(0, 25, 1)
mCols = brewer.pal(9, "GnBu")
cr0 = (colorRampPalette((mCols))(n = 256))
cr = colorRampPalette(c(cr0),
bias = 0.7, space = "rgb")
MyMatrix = matrix(nrow=7, ncol=2)
rownames(MyMatrix) = rownames(coordinates(Domain))
MyMatrix[,1] = c(-89.07686, -86.11260, -91.91021, -83.82014, -91.31258, -83.46125, -89.51171)
MyMatrix[,2] = c(41.80000, 41.36514, 42.28710, 43.0, 47.5041, 41.61285, 44.63285)
Hare.LLpnt = spTransform(Fit.pnt, proj4string(States))
mSnow5Yr.crp = crop(mSnow5Yr, Domain)
levelplot(mSnow5Yr.crp,
margin = FALSE,
xlab = NULL,
ylab = NULL,
maxpixels = 1e5,
sub= "Snow Weeks",
col.regions = cr, at = rng,
colorkey = list(space = "bottom"),
par.strip.text = list(fontface='bold', cex=4),
par.settings = list(axis.line = list(col = "black"),
strip.background = list(col = 'transparent'),
strip.border = list(col = 'transparent')),
scales = list(cex = 1.25)) +
latticeExtra::layer(sp.polygons(Domain, col = "black", lwd = 0.5)) +
latticeExtra::layer(sp.polygons(LakesLL, fill = "lightblue", col = "transparent", lwd = 0.5)) +
latticeExtra::layer(sp.text(MyMatrix, txt = Domain$Name,
pos =c(1,3,3,2,2,1,1),
col="black",font=list(face="bold"), cex=1)) +
latticeExtra::layer({SpatialPolygonsRescale(layout.north.arrow(),
offset = c(-83, 45.5),
scale = 2)})
10.2 Maximum Temperature
rng = seq(10, 25, 1)
mCols = brewer.pal(9, "YlOrRd")
cr0 = (colorRampPalette((mCols))(n = 256))
cr = colorRampPalette(c(cr0),
bias = 0.7, space = "rgb")
MyMatrix = matrix(nrow=7, ncol=2)
rownames(MyMatrix) = rownames(coordinates(Domain))
MyMatrix[,1] = c(-89.07686, -86.11260, -91.91021, -83.82014, -91.31258, -83.46125, -89.51171)
MyMatrix[,2] = c(41.80000, 41.36514, 42.28710, 43.0, 47.5041, 41.61285, 44.63285)
Hare.LLpnt = spTransform(Fit.pnt, proj4string(States))
mxTemp.crp = crop(mxTemp, Domain)
levelplot(mxTemp.crp,
margin = FALSE,
xlab = NULL,
ylab = NULL,
maxpixels = 1e6,
sub= "Maximum Temperature",
col.regions = cr, at = rng,
colorkey = list(space = "bottom"),
par.strip.text = list(fontface='bold', cex=4),
par.settings = list(axis.line = list(col = "black"),
strip.background = list(col = 'transparent'),
strip.border = list(col = 'transparent')),
scales = list(cex = 1.25)) +
latticeExtra::layer(sp.polygons(Domain, col = "black", lwd = 0.5)) +
latticeExtra::layer(sp.polygons(LakesLL, fill = "lightblue", col = "transparent", lwd = 0.5)) +
latticeExtra::layer(sp.text(MyMatrix, txt = Domain$Name,
pos =c(1,3,3,2,2,1,1),
col="black",font=list(face="bold"), cex=1)) +
latticeExtra::layer({SpatialPolygonsRescale(layout.north.arrow(),
offset = c(-83, 45.5),
scale = 2)})