Regional Hare Occupancy
Prior Versions:
1. Original model (no validation data) here: https://rpubs.com/JMHumphreys/Hare1
2. Original with validation data here: https://rpubs.com/JMHumphreys/Hare2
3. Revised version of original with validation here: https://rpubs.com/JMHumphreys/Hare3
4. Large domain version here: http://rpubs.com/JMHumphreys/HareL
5. Reduced domain version here: https://rpubs.com/JMHumphreys/HareR
6. 15km Buffer version here: https://rpubs.com/JMHumphreys/Hare15k
1 Data Prep
1.0.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.1 Get 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)
#Point grid version
Grd.pnt = rasterToPoints(Domain.r, spatial = TRUE)
Grd.pnt@data = Grd.pnt@data %>%
mutate(Long = Grd.pnt@coords[,1],
Lat = Grd.pnt@coords[,2],
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 and Define Physical Barriers
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 = 8,
min.angle = 25,
offset = c(max.edge, bound.outer))mesh$n #number of nodes## [1] 5551plot(mesh, lwd=0.5)
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
Red areas are 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)
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))Identify 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, "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.1 Seperate Datasets
Hare.mod = subset(All.pnts, Spp == "Hare") #Observations only
Mesh.mod = subset(All.pnts, Spp == "Mesh" & is.na(Domain) == FALSE) #Mesh locations excluding buffer extension
HareMesh.mod = subset(All.pnts, Spp != "Grid" & Spp != "Valid") #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 4739
## 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)
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
## 1 1.000 0.059 0.048 0.214 0.064
## 2 1.547 0.832 0.003 0.786 0.001
## 3 3.394 0.109 0.949 0.000 0.9346 Model Set up
Define Regions (Three Regions)
Converting each region to an integer value to be used to internally “replicate” estimates by region.
HR.df = HareMesh.mod@data
levels(factor(HR.df$StateUP)) #Level names## [1] "All" "Mich.UP" "Michigan" "Wisconsin"HR.df$StateUPWI = as.integer(as.factor(HR.df$StateUP)) #convert to integer
levels(factor(HR.df$StateUPWI)) #levels as integer## [1] "1" "2" "3" "4"#Keep U.P. Hare observations as "1", set other regions to "0"
MupSet = HR.df
MupSet$OBS2 = ifelse(MupSet$StateUPWI == "2", MupSet$OBS, 0)
#Keep Lower MI Hare observations as "1", set other regions to "0"
MSet = HR.df
MSet$OBS2 = ifelse(MSet$StateUPWI == "3", MSet$OBS, 0)
#Keep Wisconsin Hare observations as "1", set other regions to "0"
WSet = HR.df
WSet$OBS2 = ifelse(WSet$StateUPWI == "4", WSet$OBS, 0)
MupSet$Rep = 1 #Renumbering from 1-3 instead of 2-4
MSet$Rep = 2
WSet$Rep = 3
HR.df = rbind(MupSet, MSet, WSet) #Join data
HR.df %>% #Count of hare observations by region
group_by(Rep) %>%
summarise(Cnt = sum(OBS2))## # A tibble: 3 x 2
## Rep Cnt
## <dbl> <dbl>
## 1 1 38
## 2 2 38
## 3 3 536.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
locs = cbind(HR.df[,"Long"], HR.df[,"Lat"])
A.env = inla.spde.make.A(mesh,
alpha = 2,
loc=locs)
#Prior for barrier model
barrier.model = inla.barrier.pcmatern(mesh, #mesh
barrier.triangles = barrier.triangles, #Lake boundaries
prior.range = c(1000, 0.5), #0.5 probabilty of effect within 1000 km
prior.sigma = c(1, 0.01))
#Same as above, but without barriers
spde = inla.spde2.pcmatern(mesh,
prior.range = c(1000, 0.5),
prior.sigma = c(1, 0.01))
#Create index to track locations of mesh nodes
field.det = inla.spde.make.index("field.det", spde$n.spde) #index for detection level
field.det.c = inla.spde.make.index("field.det.c", spde$n.spde) #copy of above to pass to covariate level
field.env = inla.spde.make.index("field.env", spde$n.spde) #index for covariate level6.2 Seperate Set for Base Model
Base.set = subset(HareMesh.mod, Source != "Grid")
locs = cbind(Base.set@coords[,1], Base.set@coords[,2])
A.base = inla.spde.make.A(mesh,
alpha = 2,
loc=locs)
field.base = inla.spde.make.index("field.base", spde$n.spde)
base.lst = list(c(field.base,
list(intercept3 = 1)),
list(XX = Base.set@data[,"Long"])) #Placeholder only
base.stk = inla.stack(data = list(Y = Base.set$OBS), #Standard model for comparison
A = list(A.base, 1),
effects = base.lst,
tag = "base.0")7 Model Runs
7.1 Three Region Model
Three Region Data Prep
Need to use list objects rather than data frames.
#Detection level
DT.df = Hare.mod@data
DT.lst = list(c(field.det, #Spatial index
list(intercept1 = 1)), #Intercept
list(XX = DT.df[,"Lat"], #List of variables/covariates (placeholder for detection)
Site = DT.df[,"Site2"])) #Site identifier (to allow sites to idependently vary)
detect.stk = inla.stack(data = list(Y = DT.df$Counts, #Number of successes
Field.trials = DT.df$Trials), #trials
A = list(A.det, 1), #Projection matrix
effects = DT.lst, #Spatial index and covariates
tag = "det.0") #just a name/tag
jdetect.stk = inla.stack(data = list(Y = cbind(DT.df$Counts, NA), #"NA" = space for next model level
Field.trials = DT.df$Trials), #as above
A = list(A.det, 1),
effects = DT.lst,
tag = "jdet.0")
#Environment/covariate level
HR.lst = list(c(field.env, #index for covariate level
field.det.c, #copy of field from detection level
list(intercept2 = 1)), #intercept
list(mSnw5yr = round(HR.df[,"mSnow5yrE"], 3), #Snow weeks
mMxTemp = round(HR.df[,"mxTempE"], 3), #Temperature
Forest = HR.df[,"Forest1kmE"], #Forest
Region = HR.df[,"Rep"])) #Region identifier
env.stk = inla.stack(data = list(Y = HR.df$OBS, #Standard model for comparison
Field.trials = rep(1, dim(HR.df)[1])), #1 trial
A = list(A.env, 1),
effects = HR.lst,
tag = "env.0")
jenv.stk = inla.stack(data = list(Y = cbind(NA, HR.df$OBS2), #NA = space for detection level
Field.trials = rep(1, dim(HR.df)[1])),
A = list(A.env, 1),
effects = HR.lst,
tag = "jenv1.0")
#Combine detection and covariate levels
Joint.stk = inla.stack(jdetect.stk, jenv.stk)
#Save data to run in HPCC
#save(list=c("Joint.stk", "barrier.model", "spde"), file="./HPC/Feb10/Comb_0210.RData") #File for HPCCThree Regions Model The joint models are computationaly expensive, so there were ran on the HPCC.
Site.prior = list(theta=list(prior = "normal", param=c(0, 3))) #Prior for site effect
Reg.prior = list(theta=list(prior = "normal", param=c(0, 3))) #prior for region
JFrm0 = Y ~ -1 + intercept1 + #intercept (detection level)
intercept2 + #Intercept (covariate level)
f(field.det, #spatial index (detection)
model=barrier.model) + #change "barrier.model" to "spde" for no-barrier version.
f(field.det.c, #Shared spatial field to account for correlation between model levels
copy = "field.det", #copied from detection level of model
fixed = FALSE) +
f(field.env, #spatial index for covariate level
model=barrier.model) +
f(Site, #Site identifier
model="iid",
hyper=Site.prior) + #site prior
f(mMxTemp, #temperature
model="rw1", #random walk of order 1
replicate = Region, #replicate for each region
hyper = Reg.prior) + #prior
f(mSnw5yr, #weekly snowfall
model="rw1",
replicate = Region,
hyper = Reg.prior) +
f(Forest, #Forest
model="iid",
replicate = Region,
hyper = Reg.prior)
#thetaJ = JModel.full$internal.summary.hyperpar$mean #from initial run (to speed up model)
thetaJ = c(1.1935975, 3.7979780, -2.1432853, 6.6863289, -0.6710575, -0.1090578, -0.9891725, 0.6587151, 1)
JModel.full = inla(JFrm0, #formula
data = inla.stack.data(Joint.stk), #data list object
family = c("binomial","binomial"), #families for detection and covariate levels
verbose = FALSE, #Show running process
Ntrials = inla.stack.data(Joint.stk)$Field.trials, #number trials (varies for detection level)
control.fixed = list(prec = 0.001, #priors for intercept
prec.intercept = 0.0001),
control.predictor = list(
A = inla.stack.A(Joint.stk), #data again
compute = TRUE, #estimate fitted values
link = 1), #transform fitted values from logit
control.mode = list(restart = TRUE, theta = thetaJ), #to speed up
control.inla = list(strategy="gaussian", #to speed up
int.strategy = "eb"),
control.results = list(return.marginals.random = TRUE, #results to report
return.marginals.predictor = TRUE),
control.compute=list(dic = TRUE, cpo = TRUE, waic = TRUE)) #calculate indicies for comparisonLoad Results from HPCC
load("~/Michigan_State/Sean/HPC/Feb10/Full_0210_RES.RData")
GetMets(JModel.full)[,c(1,3)]## Metric Tier.2
## 1 DIC 1207.9115
## 2 WAIC 1094.2853
## 3 lCPO 0.03137.2 Two Region Models
Treating Wisconsin and the UP as a single/combined region.
HR2.df = HareMesh.mod@data
levels(factor(HR2.df$StateUP)) #Level names## [1] "All" "Mich.UP" "Michigan" "Wisconsin"HR2.df$StateUPWI = as.integer(as.factor(HR2.df$StateUP)) #convert to integer
levels(factor(HR2.df$StateUPWI)) #levels as integer## [1] "1" "2" "3" "4"#Keep U.P. and WI Hare observations as "1", set Michigan (Lower) locations to "0"
UP_WI_Set = HR2.df
UP_WI_Set$OBS3 = ifelse(UP_WI_Set$StateUPWI == "2" |
UP_WI_Set$StateUPWI == "4", UP_WI_Set$OBS, 0)
#Keep Lower MI Hare observations as "1", set other regions to "0"
MI_Set = HR2.df
MI_Set$OBS3 = ifelse(MI_Set$StateUPWI == "3", MI_Set$OBS, 0)
UP_WI_Set$Rep = 1 #Renumbering from 1-3 instead of 2-4
MI_Set$Rep = 2
HR2.df = rbind(UP_WI_Set, MI_Set) #Join data
HR2.df %>% #Count of hare observations by region
group_by(Rep) %>%
summarise(Cnt = sum(OBS3))## # A tibble: 2 x 2
## Rep Cnt
## <dbl> <dbl>
## 1 1 91
## 2 2 38Two Regions Model Set-up As above.
#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
locs = cbind(HR2.df[,"Long"], HR2.df[,"Lat"])
A.env = inla.spde.make.A(mesh,
alpha = 2,
loc=locs)
#Prior for barrier model
barrier.model = inla.barrier.pcmatern(mesh, #mesh
barrier.triangles = barrier.triangles, #Lake boundaries
prior.range = c(1000, 0.5), #0.5 probabilty of effect within 1000 km
prior.sigma = c(1, 0.01))
#Same as above, but without barriers
spde = inla.spde2.pcmatern(mesh,
prior.range = c(1000, 0.5),
prior.sigma = c(1, 0.01))
#Create index to track locations of mesh nodes
field.det = inla.spde.make.index("field.det", spde$n.spde) #index for detection level
field.det.c = inla.spde.make.index("field.det.c", spde$n.spde) #copy of above to pass to covariate level
field.env = inla.spde.make.index("field.env", spde$n.spde) #index for covariate levelOrganize Two Regions Data for Fitting
#Detection level
DT.df = Hare.mod@data
DT.lst = list(c(field.det, #Spatial index
list(intercept1 = 1)), #Intercept
list(XX = DT.df[,"Lat"], #List of variables/covariates (placeholder for detection)
Site = DT.df[,"Site2"])) #Site identifier (to allow sites to idependently vary)
jdetect.stk = inla.stack(data = list(Y = cbind(DT.df$Counts, NA), #"NA" = space for next model level
Field.trials = DT.df$Trials), #as above
A = list(A.det, 1),
effects = DT.lst,
tag = "jdet.0")
#Environment/covariate level
HR.lst = list(c(field.env, #index for covariate level
field.det.c, #copy of field from detection level
list(intercept2 = 1)), #intercept
list(mSnw5yr = round(HR2.df[,"mSnow5yrE"], 3), #Snow weeks
mMxTemp = round(HR2.df[,"mxTempE"], 3), #Temperature
Forest = HR2.df[,"Forest1kmE"], #Forest
Region = HR2.df[,"Rep"])) #Region identifier
Env.stk = inla.stack(data = list(Y = HR2.df$OBS3),
A = list(A.env, 1),
effects = HR.lst,
tag = "env.0")
jenv.stk = inla.stack(data = list(Y = cbind(NA, HR2.df$OBS3), #NA = space for detection level
Field.trials = rep(1, dim(HR2.df)[1])),
A = list(A.env, 1),
effects = HR.lst,
tag = "jenv.0")
#Combine detection and covariate levels
Joint2.stk = inla.stack(jdetect.stk, jenv.stk)
#Keep copy for in sample validation
Keep.valid.df = HR2.df
#Save data to run in HPCC
#save(list=c("Joint2.stk", "barrier.model", "spde"), file="./HPC/Feb25/TwoReg.RData") #File for HPCC
#save(list=c("Joint2.stk", "barrier.model", "spde"), file="./HPC/Mar10/TwoReg2.RData") #File for HPCC
#save(list=c("Env.stk", "barrier.model"), file="./HPC/June19/Env.RData") #File for HPCCTwo Regions (Non Spatial w/ Covariates)
Site.prior = list(theta=list(prior = "normal", param=c(0, 3))) #Prior for site effect
Reg.prior = list(theta=list(prior = "normal", param=c(0, 3))) #prior for region
JFrm0 = Y ~ -1 + intercept2 + #Intercept (covariate level)
f(mMxTemp, #temperature
model="rw1", #random walk of order 1
replicate = Region, #replicate for each region
hyper = Reg.prior) + #prior
f(mSnw5yr, #weekly snowfall
model="rw1",
replicate = Region,
hyper = Reg.prior) +
f(Forest, #Forest
model="iid",
replicate = Region,
hyper = Reg.prior)
JModel.2Reg.ns = inla(JFrm0, #formula
data = inla.stack.data(Env.stk), #data list object
family = "binomial", #families for detection and covariate levels
verbose = FALSE, #Show running process
#Ntrials = inla.stack.data(Joint2.stk)$Field.trials, #number trials
control.fixed = list(prec = 0.001, #priors for intercept
prec.intercept = 0.0001),
control.predictor = list(
A = inla.stack.A(Env.stk), #data again
compute = TRUE, #estimate fitted values
link = 1), #transform fitted values from logit
#control.mode = list(restart = TRUE, theta = thetaJ), #to speed up
control.inla = list(strategy="gaussian", #to speed up
int.strategy = "eb"),
control.results = list(return.marginals.random = TRUE, #results to report
return.marginals.predictor = TRUE),
control.compute=list(dic = TRUE, cpo = TRUE, waic = TRUE)) #calculate indicies for comparisonload("~/Michigan_State/Sean/HPC/June19/JModel_2Reg_ns.RData")
GetMets(JModel.2Reg.ns)## Metric Tier.1
## 1 DIC 1893.4959
## 2 WAIC 1837.0123
## 3 lCPO 0.0049Two Regions (Spatial, Covariates, no detection)
Site.prior = list(theta=list(prior = "normal", param=c(0, 3))) #Prior for site effect
Reg.prior = list(theta=list(prior = "normal", param=c(0, 3))) #prior for region
JFrm0 = Y ~ -1 + intercept2 + #Intercept (covariate level)
f(field.env,
model=barrier.model) +
f(mMxTemp, #temperature
model="rw1", #random walk of order 1
replicate = Region, #replicate for each region
hyper = Reg.prior) + #prior
f(mSnw5yr, #weekly snowfall
model="rw1",
replicate = Region,
hyper = Reg.prior) +
f(Forest, #Forest
model="iid",
replicate = Region,
hyper = Reg.prior)
JModel.2Reg.nd = inla(JFrm0, #formula
data = inla.stack.data(Env.stk), #data list object
family = "binomial", #families for detection and covariate levels
verbose = FALSE, #Show running process
#Ntrials = inla.stack.data(Joint2.stk)$Field.trials, #number trials
control.fixed = list(prec = 0.001, #priors for intercept
prec.intercept = 0.0001),
control.predictor = list(
A = inla.stack.A(Env.stk), #data again
compute = TRUE, #estimate fitted values
link = 1), #transform fitted values from logit
#control.mode = list(restart = TRUE, theta = thetaJ), #to speed up
control.inla = list(strategy="gaussian", #to speed up
int.strategy = "eb"),
control.results = list(return.marginals.random = TRUE, #results to report
return.marginals.predictor = TRUE),
control.compute=list(dic = TRUE, cpo = TRUE, waic = TRUE)) #calculate indicies for comparisonload("~/Michigan_State/Sean/HPC/June19/Env_061919_RES.RData")
GetMets(JModel.2Reg.nd)## Metric Tier.1
## 1 DIC 1368.2536
## 2 WAIC 1318.4894
## 3 lCPO 0.0036Two Regions (No Covariates)
Site.prior = list(theta=list(prior = "normal", param=c(0, 3)))
Reg.prior = list(theta=list(prior = "normal", param=c(0, 3)))
JFrm0 = Y ~ -1 + intercept1 +
intercept2 +
f(field.det,
model=barrier.model) +
f(field.det.c,
copy = "field.det",
fixed = FALSE) +
f(field.env,
model=barrier.model) +
f(Site,
model="iid",
hyper=Site.prior)
JModel.2Regst = inla(JFrm0,
#num.threads = 8,
data = inla.stack.data(Joint2.stk),
family = c("binomial","binomial"),
verbose = FALSE,
Ntrials = inla.stack.data(Joint2.stk)$Field.trials, #trials
control.fixed = list(prec = 0.001,
prec.intercept = 0.0001),
control.predictor = list(
A = inla.stack.A(Joint2.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)) load("~/Michigan_State/Sean/HPC/June19/TwoReg_0619st_RES.RData")
GetMets(JModel.2Regst)[,c(1,3)]## Metric Tier.2
## 1 DIC 1358.7232
## 2 WAIC 1231.6249
## 3 lCPO 0.0563Two Regions Model (full)
Site.prior = list(theta=list(prior = "normal", param=c(0, 3))) #Prior for site effect
Reg.prior = list(theta=list(prior = "normal", param=c(0, 3))) #prior for region
JFrm0 = Y ~ -1 + intercept1 + #intercept (detection level)
intercept2 + #Intercept (covariate level)
f(field.det, #spatial index (detection)
model=barrier.model) + #change "barrier.model" to "spde" for no-barrier version.
f(field.det.c, #Shared spatial field to account for correlation between model levels
copy = "field.det", #copied from detection level of model
fixed = FALSE) +
f(field.env, #spatial index for covariate level
model=barrier.model) +
f(Site, #Site identifier
model="iid",
hyper=Site.prior) + #site prior
f(mMxTemp, #temperature
model="rw1", #random walk of order 1
replicate = Region, #replicate for each region
hyper = Reg.prior) + #prior
f(mSnw5yr, #weekly snowfall
model="rw1",
replicate = Region,
hyper = Reg.prior) +
f(Forest, #Forest
model="iid",
replicate = Region,
hyper = Reg.prior)
#thetaJ = JModel.2Reg$internal.summary.hyperpar$mean #from initial run (to speed up model)
thetaJ = c(1.1410600, 3.9267409, -1.9521533, 6.4758302, -0.9262989, -0.4777043, 0.2826159, 0.4211844, 0.486779)
JModel.2Reg = inla(JFrm0, #formula
data = inla.stack.data(Joint2.stk), #data list object
family = c("binomial","binomial"), #families for detection and covariate levels
verbose = FALSE, #Show running process
Ntrials = inla.stack.data(Joint2.stk)$Field.trials, #number trials
control.fixed = list(prec = 0.001, #priors for intercept
prec.intercept = 0.0001),
control.predictor = list(
A = inla.stack.A(Joint2.stk), #data again
compute = TRUE, #estimate fitted values
link = 1), #transform fitted values from logit
control.mode = list(restart = TRUE, theta = thetaJ), #to speed up
control.inla = list(strategy="gaussian", #to speed up
int.strategy = "eb"),
control.results = list(return.marginals.random = TRUE, #results to report
return.marginals.predictor = TRUE),
control.compute=list(dic = TRUE, cpo = TRUE, waic = TRUE)) #calculate indicies for comparisonload("~/Michigan_State/Sean/HPC/Mar10/TwoReg_0310_RES.RData")
GetMets(JModel.2Reg)[,c(1,3)]## Metric Tier.2
## 1 DIC 1127.2103
## 2 WAIC 1024.0180
## 3 lCPO 0.04387.3 Combined Model
Site.prior = list(theta=list(prior = "normal", param=c(0, 3))) #Prior for site effect
Reg.prior = list(theta=list(prior = "normal", param=c(0, 3))) #prior for region
JFrm0 = Y ~ -1 + intercept1 +
intercept2 +
f(field.det,
model=barrier.model) +
f(field.det.c,
copy = "field.det",
fixed = FALSE) +
f(field.env,
model=barrier.model) +
f(Site,
model="iid",
hyper=Site.prior) +
f(mMxTemp,
model="rw1",
#replicate = Region,
hyper = Reg.prior) +
f(mSnw5yr,
model="rw1",
#replicate = Region,
hyper = Reg.prior) +
f(Forest,
model="iid",
#replicate = Region,
hyper = Reg.prior)
thetaJ = c(1.1410600, 3.9267409, -1.9521533, 6.4758302, -0.9262989, -0.4777043, 0.2826159, 0.4211844, 0.486779)
JModel.comb = inla(JFrm0, #formula
data = inla.stack.data(Joint2.stk),
family = c("binomial","binomial"),
verbose = FALSE,
Ntrials = inla.stack.data(Joint2.stk)$Field.trials, #number trials
control.fixed = list(prec = 0.001,
prec.intercept = 0.0001),
control.predictor = list(
A = inla.stack.A(Joint2.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)) load("~/Michigan_State/Sean/HPC/June19/Comb_062019B_RES.RData")
GetMets(JModel.comb)## Metric Tier.1 Tier.2
## 1 DIC 844.7262 1251.3275
## 2 WAIC 1246.0023 1138.4055
## 3 lCPO 16.9696 0.04867.4 Base Models
Model Spatial Fields for Comparison.
Spatial Field w/ Barriers
JFrmB = Y ~ -1 + intercept3 +
f(field.base,
model=barrier.model)
#thetaJ = JModel.base$internal.summary.hyperpar$mean #from initial run (to speed up model)
thetaJ = c(1.444903, 4.073292)
JModel.base = inla(JFrmB,
data = inla.stack.data(base.stk),
family = "binomial",
verbose = FALSE,
control.fixed = list(prec = 0.001,
prec.intercept = 0.0001),
control.predictor = list(
A = inla.stack.A(base.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))Spatial Field w/out Barriers
JFrmNB = Y ~ -1 + intercept3 +
f(field.base,
model=spde)
thetaJ = c(1.444903, 4.073292)
JModel.base.nb = inla(JFrmNB,
data = inla.stack.data(base.stk),
family = "binomial",
verbose = FALSE,
control.fixed = list(prec = 0.001,
prec.intercept = 0.0001),
control.predictor = list(
A = inla.stack.A(base.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)) 8 Model Results
Map Spatial Fields
SF.Ext = extent(c(-93.173415, -81.932841, 41.217132, 47.740649))
SF.Ext = as(SF.Ext, 'SpatialPolygons')
proj4string(SF.Ext) = proj4string(Domain)
SF.r = raster(res = 0.02, ext = extent(SF.Ext),
crs = proj4string(SF.Ext))
SF.domain.r = rasterize(SF.Ext, SF.r,
field = 0,
background = NA)
SF.pnt = rasterToPoints(SF.domain.r, spatial = TRUE)
SF.pnt@data = SF.pnt@data %>%
mutate(Long = SF.pnt@coords[,1],
Lat = SF.pnt@coords[,2],
Spp = "SFGrid") %>%
select(-layer)
Map.pnts = SF.pnt
ModResult1 = JModel.base #base, barrier
ModResult2 = JModel.base.nb # base, no barrier
ModResult3 = JModel.2Reg
Map.pntsp = spTransform(Map.pnts, proj4string(HareMesh.mod))
pLoc = cbind(Map.pntsp@coords[,1], Map.pntsp@coords[,2])
Ap = inla.spde.make.A(mesh, loc=pLoc)
Map.pnts$wBarrier = drop(Ap %*% ModResult1$summary.random$field.base$mean)
Map.pnts$noBarrier = drop(Ap %*% ModResult2$summary.random$field.base$mean)
Map.pnts$Reg2SF = drop(Ap %*% ModResult3$summary.random$field.env$mean)
#Create rasters
SF.barrier.r = rasterize(spTransform(
Map.pnts,
CRS(proj4string(SF.domain.r))),
SF.domain.r,
"wBarrier",
background = NA)
SF.nobarrier.r = rasterize(spTransform(
Map.pnts,
CRS(proj4string(SF.domain.r))),
SF.domain.r,
"noBarrier",
background = NA)
SF.Reg2.r = rasterize(spTransform(
Map.pnts,
CRS(proj4string(SF.domain.r))),
SF.domain.r,
"Reg2SF",
background = NA)
Barrier.comp = stack(SF.nobarrier.r, SF.barrier.r, SF.Reg2.r)
names(Barrier.comp) = c("NoBarrier", "Barrier", "Region2")
Domain2 = subset(Domain, ID != 21 & ID != 48)
rng = seq(-5.5, 5.6, 0.01)
mCols = rev(inferno(25))
cr = colorRampPalette(c(mCols),
bias = 0.3)
levelplot(Barrier.comp, #[[1:2]],
layout = c(3,1),
margin = FALSE,
xlab = NULL,
ylab = NULL,
#main = "Occurrence Probabilty",
names.attr= c("A", "B", "C"),
maxpixels = 1e5,
col.regions = cr, at = rng,
colorkey = list(labels=list(
fontface='bold', cex=1.5),
labels=list(cex=16),
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 = "black", fill="transparent", lwd = 1)) +
latticeExtra::layer(sp.polygons(LakesLL , fill = "transparent", alpha = 1, lwd=1, col = "black")) +
latticeExtra::layer({SpatialPolygonsRescale(layout.north.arrow(),
offset = c(-83.5, 46.3),
scale = c(1.9,1.3))})
Model Comparison Table
#Environment Level
Compare.tab = as.data.frame(matrix(ncol = 3, nrow= 6))
names(Compare.tab) = c("Model", "DIC", "WAIC")
Compare.tab[1, 1] = "Two Region (Non-spatial)"
Compare.tab[1, 2] = round(GetMets(JModel.2Reg.ns)[1,2],2)
Compare.tab[1, 3] = round(GetMets(JModel.2Reg.ns)[2,2],2)
Compare.tab[2, 1] = "Two Region (Spatial only)"
Compare.tab[2, 2] = round(GetMets(JModel.2Regst)[1,3],2)
Compare.tab[2, 3] = round(GetMets(JModel.2Regst)[2,3],2)
Compare.tab[3, 1] = "Two Region (No detection)"
Compare.tab[3, 2] = round(GetMets(JModel.2Reg.nd)[1,2],2)
Compare.tab[3, 3] = round(GetMets(JModel.2Reg.nd)[2,2],2)
Compare.tab[4, 1] = "Two Region"
Compare.tab[4, 2] = round(GetMets(JModel.2Reg)[1,3],2)
Compare.tab[4, 3] = round(GetMets(JModel.2Reg)[2,3],2)
Compare.tab[5, 1] = "Three Region"
Compare.tab[5, 2] = round(GetMets(JModel.full)[1,3],2)
Compare.tab[5, 3] = round(GetMets(JModel.full)[2,3],2)
Compare.tab[6, 1] = "Combined"
Compare.tab[6, 2] = round(GetMets(JModel.comb)[1,3],2)
Compare.tab[6, 3] = round(GetMets(JModel.comb)[2,3],2)
kable(Compare.tab, caption = "Model Comparison") %>%
kable_styling("striped", full_width = F) %>%
row_spec(0, font_size = 20) %>%
column_spec(1, bold = T)| Model | DIC | WAIC |
|---|---|---|
| Two Region (Non-spatial) | 1893.50 | 1837.01 |
| Two Region (Spatial only) | 1358.72 | 1231.62 |
| Two Region (No detection) | 1368.25 | 1318.49 |
| Two Region | 1127.21 | 1024.02 |
| Three Region | 1207.91 | 1094.29 |
| Combined | 1251.33 | 1138.41 |
#SelMod.fx = xtable(Compare.tab)
#print(SelMod.fx, include.rownames = F)Detection Site Effect
PhySig.df2 = as.data.frame(JModel.2Reg$summary.random$Site)
names(PhySig.df2) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975", "mode", "kld")
PhySig.df2$PSigL = ifelse(PhySig.df2$Q025>0 & PhySig.df2$Q975>0, 1, 0)
PhySig.df2$PSigH = ifelse(PhySig.df2$Q025<0 & PhySig.df2$Q975<0, 1, 0)
PhySig.df3 = PhySig.df2 %>%
filter(PSigL == 1 | PSigH == 1)
RM.m2 = ggplot(PhySig.df2, aes(x=ID, y=Mean)) +
geom_point(size=2, pch=1, col = "gray75") +
geom_linerange(aes(ymin=Q025, ymax=Q975), colour="gray75") +
geom_point(data=PhySig.df3, aes(x=ID, y=Mean),
size=2, pch=19, col = "red") +
geom_linerange(data=PhySig.df3, aes(ymin=Q025, ymax=Q975), colour="black") +
geom_hline(yintercept = 0,
linetype = "solid",
col = "red",
size = 0.5) +
ylim(-5,5) +
xlab("Survey Site Location") +
ylab("Detection Probabilty (logit)") +
theme_classic() +
theme(plot.margin = unit(c(1,1,1,1),"cm"),
panel.grid.minor = element_blank(),
panel.grid.major = element_blank(),
panel.background = element_blank(),
plot.background = element_blank(),
panel.border = element_blank(),
legend.direction = "vertical",
#legend.position = "bottom",
strip.text = element_text(size=16, face="bold"),
strip.background = element_blank(),
legend.position=c(0.74, 0.25),
legend.key.size = unit(3,"line"),
axis.ticks.x=element_blank(),
legend.text = element_text(size=14, face="bold"),
legend.title = element_text(size=18, face="bold"),
axis.text.x = element_blank(),
axis.text.y =element_text(size=18, face="bold"),
axis.title.x = element_text(size=22, face="bold"),
axis.title.y =element_text(size=22, face="bold"),
plot.title = element_blank())
RM.m2
Average Maximum Temperature
Comb.plt.df = as.data.frame(JModel.2Reg$summary.random$mMxTemp)[,1:6]
names(Comb.plt.df) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975") #Clean labels
Comb.plt.df$Region = rep(c("Wisconsin + Upper", "Lower"), each = 3502) #Add Region ID
range(Hare.mod$mxTempE)## [1] 14.23175 18.58800Plt.df = Comb.plt.df
Plt.df$Region = as.factor(Plt.df$Region)
RAnge = range(Plt.df$ID)
MySpan = 0.2
#Set 1
g1 = Plt.df %>% filter(Region == "Lower")
g1 = ggplot(g1) +
stat_smooth(aes(x = ID, y = Q025), span = MySpan, method = "loess", fullrange=TRUE, se = FALSE) +
stat_smooth(aes(x = ID, y = Q975), span = MySpan, method = "loess", fullrange=TRUE, se = FALSE)
gg1 = ggplot_build(g1)
df1 = data.frame(x = gg1$data[[1]]$x,
ymin = gg1$data[[1]]$y,
ymax = gg1$data[[2]]$y)
df1$Region = "Lower Peninsula"
#Set 2
g2 = Plt.df %>% filter(Region == "Wisconsin + Upper")
g2 = ggplot(g2) +
stat_smooth(aes(x = ID, y = Q025), span = MySpan, method = "loess", fullrange=TRUE, se = FALSE) +
stat_smooth(aes(x = ID, y = Q975), span = MySpan, method = "loess", fullrange=TRUE, se = FALSE)
gg2 = ggplot_build(g2)
df2 = data.frame(x = gg2$data[[1]]$x,
ymin = gg2$data[[1]]$y,
ymax = gg2$data[[2]]$y)
df2$Region = "Wisconsin + Upper Peninsula"
Plt.df$Region = ifelse(Plt.df$Region == "Lower", "Lower Peninsula", "Wisconsin + Upper Peninsula")
#Combine
Plt.dft = rbind(df1, df2)
ggplot() +
geom_smooth(data = Plt.df,
aes(ID, Mean, group = Region,
col = Region,
linetype= Region),
fullrange=TRUE,
method = "loess",
span = MySpan,
se = FALSE,
lwd = 1) +
geom_ribbon(data = Plt.dft,
aes(x = x, ymin=ymin, ymax=ymax, group = Region),
alpha=0.3, #transparency
linetype=1, #solid, dashed or other line types
colour="transparent", #border line color
size=1, #border line size
fill="darkgray") + #fill color
scale_linetype_manual(values=c("solid", "dotted")) +
scale_colour_manual(values=c("black", "black")) +
geom_hline(yintercept = 0,
linetype = "solid",
col = "red",
size = 0.5) +
scale_x_continuous(expand=c(0,0), limits=c(13.5,19.5)) +
scale_y_continuous(expand=c(0,0), limits=c(-5.9,5)) +
xlab("Maximum Temperature (°C)") +
ylab("Occurrence Probability (logit)") +
theme_classic() +
theme(plot.margin = unit(c(1,1,1,1),"cm"),
panel.grid.minor = element_blank(),
panel.grid.major = element_blank(),
panel.background = element_blank(),
plot.background = element_blank(),
panel.border = element_blank(),
legend.direction = "vertical",
#legend.position = "bottom",
strip.text = element_text(size=16, face="bold"),
strip.background = element_blank(),
legend.position=c(0.74, 0.25),
legend.key.size = unit(3,"line"),
legend.text = element_text(size=14, face="bold"),
legend.title = element_text(size=18, face="bold"),
axis.text.x = element_text(size=18, face="bold"),
axis.text.y =element_text(size=18, face="bold"),
axis.title.x = element_text(size=22, face="bold"),
axis.title.y =element_text(size=22, face="bold"),
plot.title = element_blank()) +
coord_cartesian(xlim=c(14,19), ylim=c(-5.5,4))
Snow Weeks
Comb.plt.df = as.data.frame(JModel.2Reg$summary.random$mSnw5yr)[,1:6]
names(Comb.plt.df) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975")
Comb.plt.df$Region = rep(c("Wisconsin + Upper Peninsula", "Lower Peninsula"), each = 90)
SPlt.df = Comb.plt.df
range(Hare.mod$mSnow5yrE)## [1] 5.0 21.6SPlt.df$Region = as.factor(SPlt.df$Region)
MySpan = 0.2
#Set 1
g1 = SPlt.df %>% filter(Region == "Lower Peninsula")
g1 = ggplot(g1) +
stat_smooth(aes(x = ID, y = Q025), span = MySpan, method = "loess", fullrange=TRUE, se = FALSE) +
stat_smooth(aes(x = ID, y = Q975), span = MySpan, method = "loess", fullrange=TRUE, se = FALSE)
gg1 = ggplot_build(g1)
df1 = data.frame(x = gg1$data[[1]]$x,
ymin = gg1$data[[1]]$y,
ymax = gg1$data[[2]]$y)
df1$Region = "Lower Peninsula"
#Set 2
g2 = SPlt.df %>% filter(Region == "Wisconsin + Upper Peninsula")
g2 = ggplot(g2) +
stat_smooth(aes(x = ID, y = Q025), span = MySpan, method = "loess", fullrange=TRUE, se = FALSE) +
stat_smooth(aes(x = ID, y = Q975), span = MySpan, method = "loess", fullrange=TRUE, se = FALSE)
gg2 = ggplot_build(g2)
df2 = data.frame(x = gg2$data[[1]]$x,
ymin = gg2$data[[1]]$y,
ymax = gg2$data[[2]]$y)
df2$Region = "Wisconsin + Upper Peninsula"
#Combine
SPlt.dft = rbind(df1, df2)
ggplot() +
geom_smooth(data = SPlt.df,
aes(ID, Mean, group = Region,
col = Region,
linetype= Region),
fullrange=TRUE,
method = "loess",
span = MySpan,
se = FALSE,
lwd = 1) +
geom_ribbon(data = SPlt.dft,
aes(x = x, ymin=ymin, ymax=ymax, group = Region),
alpha=0.3, #transparency
linetype=1, #solid, dashed or other line types
colour="transparent", #border line color
size=1, #border line size
fill="darkgray") + #fill color
scale_linetype_manual(values=c("solid", "dotted")) +
scale_colour_manual(values=c("black", "black")) +
geom_hline(yintercept = 0,
linetype = "solid",
col = "red",
size = 0.5) +
scale_x_continuous(expand=c(0,0), limits=c(5,22)) +
scale_y_continuous(expand=c(0,0), limits=c(-5.5,5.5)) +
xlab("Snow Weeks") +
ylab("Occurrence Probability (logit)") +
theme_classic() +
theme(plot.margin = unit(c(1,1,1,1),"cm"),
panel.grid.minor = element_blank(),
panel.grid.major = element_blank(),
panel.background = element_blank(),
plot.background = element_blank(),
panel.border = element_blank(),
legend.direction = "vertical",
#legend.position = "bottom",
strip.text = element_text(size=16, face="bold"),
strip.background = element_blank(),
legend.position=c(0.73, 0.21),
legend.key.size = unit(3,"line"),
legend.text = element_text(size=14, face="bold"),
legend.title = element_text(size=18, face="bold"),
axis.text.x = element_text(size=18, face="bold"),
axis.text.y =element_text(size=18, face="bold"),
axis.title.x = element_text(size=22, face="bold"),
axis.title.y =element_text(size=22, face="bold"),
plot.title = element_blank()) +
coord_cartesian(xlim=c(5,22), ylim=c(-5,5))## Warning: Removed 12 rows containing non-finite values (stat_smooth).
Forest Cover
Comb.plt.df = as.data.frame(JModel.2Reg$summary.random$Forest)[,1:6]
names(Comb.plt.df) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975")
Comb.plt.df$Region = rep(c("Wisconsin + Upper Peninsula", "Lower Peninsula"), each = 110)
FPlt.df = Comb.plt.df
range(Hare.mod$Forest1kmE)## [1] 2.3 10.9FPlt.df$Region = as.factor(FPlt.df$Region)
MySpan = 0.2
#Set 1
g1 = FPlt.df %>% filter(Region == "Lower Peninsula")
g1 = ggplot(g1) +
stat_smooth(aes(x = ID, y = Q025), span = MySpan, method = "loess", fullrange=TRUE, se = FALSE) +
stat_smooth(aes(x = ID, y = Q975), span = MySpan, method = "loess", fullrange=TRUE, se = FALSE)
gg1 = ggplot_build(g1)
df1 = data.frame(x = gg1$data[[1]]$x,
ymin = gg1$data[[1]]$y,
ymax = gg1$data[[2]]$y)
df1$Region = "Lower Peninsula"
#Set 2
g2 = FPlt.df %>% filter(Region == "Wisconsin + Upper Peninsula")
g2 = ggplot(g2) +
stat_smooth(aes(x = ID, y = Q025), span = MySpan, method = "loess", fullrange=TRUE, se = FALSE) +
stat_smooth(aes(x = ID, y = Q975), span = MySpan, method = "loess", fullrange=TRUE, se = FALSE)
gg2 = ggplot_build(g2)
df2 = data.frame(x = gg2$data[[1]]$x,
ymin = gg2$data[[1]]$y,
ymax = gg2$data[[2]]$y)
df2$Region = "Wisconsin + Upper Peninsula"
#Combine
FPlt.dft = rbind(df1, df2)
FPlt.dft$x = Scale(FPlt.dft$x)*100
FPlt.df$ID = Scale(FPlt.df$ID)*100
ggplot() +
geom_smooth(data = FPlt.df,
aes(ID, Mean, group = Region,
col = Region,
linetype= Region),
fullrange=TRUE,
method = "loess",
span = MySpan,
se = FALSE,
lwd = 1) +
geom_ribbon(data = FPlt.dft,
aes(x = x, ymin=ymin, ymax=ymax, group = Region),
alpha=0.3, #transparency
linetype=1, #solid, dashed or other line types
colour="transparent", #border line color
size=1, #border line size
fill="darkgray") + #fill color
scale_linetype_manual(values=c("solid", "dotted")) +
scale_colour_manual(values=c("black", "black")) +
geom_hline(yintercept = 0,
linetype = "solid",
col = "red",
size = 0.5) +
scale_x_continuous(expand=c(0,0), limits=c(0,100)) +
scale_y_continuous(expand=c(0,0), limits=c(-4,4)) +
xlab("Forest (% Cover)") +
ylab("Occurrence Probability (logit)") +
theme_classic() +
theme(plot.margin = unit(c(1,1,1,1),"cm"),
panel.grid.minor = element_blank(),
panel.grid.major = element_blank(),
panel.background = element_blank(),
plot.background = element_blank(),
panel.border = element_blank(),
legend.direction = "vertical",
#legend.position = "bottom",
strip.text = element_text(size=16, face="bold"),
strip.background = element_blank(),
legend.position=c(0.8, 0.22),
legend.key.size = unit(3,"line"),
legend.text = element_text(size=14, face="bold"),
legend.title = element_text(size=18, face="bold"),
axis.text.x = element_text(size=18, face="bold"),
axis.text.y =element_text(size=18, face="bold"),
axis.title.x = element_text(size=22, face="bold"),
axis.title.y =element_text(size=22, face="bold"),
plot.title = element_blank()) +
coord_cartesian(xlim=c(0,100), ylim=c(-4,4))
Mapped Spatial Predictions Two Region Model Based on covariates only.
ModResult = JModel.2Reg
#Update Region ID for prediction locations (MI and WI only)
MI_WI_Domain = subset(Domain, Name == "Michigan" | Name == "Wisconsin")
Pred.pnts = Grd.pnts
MI_WI_DomainP = spTransform(MI_WI_Domain, proj4string(Pred.pnts))
Pred.pnts$nDom = over(Pred.pnts, MI_WI_DomainP)[,1]
Pred.pntsN = subset(Pred.pnts, is.na(nDom) == FALSE)
Pred.pntsN$Region = ifelse(Pred.pntsN@data$nDom == 21 & Pred.pntsN@data$UP == "FALSE", "Upper", #FALSE indicates UP
ifelse(Pred.pntsN@data$nDom == 21 & Pred.pntsN@data$UP == "TRUE", "Lower",
ifelse(Pred.pntsN@data$nDom == 48, "Wisconsin", NA)))
Pred.pntsN$Region = ifelse(Pred.pntsN$Region == "Upper" | Pred.pntsN$Region == "Wisconsin", "Wisconsin + UP", Pred.pntsN$Region)
#Get RF
pLoc = cbind(Pred.pntsN@coords[,1], Pred.pntsN@coords[,2])
Ap = inla.spde.make.A(mesh, loc=pLoc)
Pred.pntsN$REg2det.rf = drop(Ap %*% ModResult$summary.random$field.det$mean)
Pred.pntsN$REg2detC.rf = drop(Ap %*% ModResult$summary.random$field.det.c$mean)
Pred.pntsN$REg2.rf = drop(Ap %*% ModResult$summary.random$field.env$mean)
#Get Temperature
Comb.plt.df = as.data.frame(JModel.2Reg$summary.random$mMxTemp)[,1:6] #Get effect estimates (Full model)
names(Comb.plt.df) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975") #Clean labels
Comb.plt.df$Region = rep(c("Wisconsin + Upper", "Lower"), each = 3502)
#Region-sepcific estimates
LW.lu = Comb.plt.df %>% filter(Region == "Lower")
WI.lu = Comb.plt.df %>% filter(Region == "Wisconsin + Upper")
LW.lu$LW.POS = 1:dim(LW.lu)[1]
WI.lu$WI.POS = 1:dim(WI.lu)[1]
#Look up Lower Temp Estimates
Pred.pntsN$LW.temp = sapply(Pred.pntsN$mxTempE, function(x)which.min(abs(x - LW.lu$ID)))
Pred.pntsN$LW.TmpEst = as.numeric(with(LW.lu,
Mean[match(Pred.pntsN$LW.temp,
LW.POS)]))
#Look up WI + UP Temp Estimates
Pred.pntsN$WI.temp = sapply(Pred.pntsN$mxTempE, function(x)which.min(abs(x - WI.lu$ID)))
Pred.pntsN$WI.TmpEst = as.numeric(with(WI.lu,
Mean[match(Pred.pntsN$WI.temp,
WI.POS)]))
#Match to region
Pred.pntsN$Temp.Full = ifelse(Pred.pntsN$Region == "Wisconsin + UP", Pred.pntsN$WI.TmpEst,
ifelse(Pred.pntsN$Region == "Lower", Pred.pntsN$LW.TmpEst, NA))
#Get Snow Week (same process as above, overwriting a few lables)
Comb.plt.df = as.data.frame(JModel.2Reg$summary.random$mSnw5yr)[,1:6]
names(Comb.plt.df) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975")
Comb.plt.df$Region = rep(c("Wisconsin + Upper", "Lower"), each = 90)
LW.lu = Comb.plt.df %>% filter(Region == "Lower")
WI.lu = Comb.plt.df %>% filter(Region == "Wisconsin + Upper")
LW.lu$LW.POS = 1:dim(LW.lu)[1]
WI.lu$WI.POS = 1:dim(WI.lu)[1]
Pred.pntsN$LW.temp = sapply(Pred.pntsN$mSnow5yrE, function(x)which.min(abs(x - LW.lu$ID)))
Pred.pntsN$LW.TmpEst = as.numeric(with(LW.lu,
Mean[match(Pred.pntsN$LW.temp,
LW.POS)]))
Pred.pntsN$WI.temp = sapply(Pred.pntsN$mSnow5yrE, function(x)which.min(abs(x - WI.lu$ID)))
Pred.pntsN$WI.TmpEst = as.numeric(with(WI.lu,
Mean[match(Pred.pntsN$WI.temp,
WI.POS)]))
#Match to region
Pred.pntsN$Snow.Full = ifelse(Pred.pntsN$Region == "Lower", Pred.pntsN$LW.TmpEst,
ifelse(Pred.pntsN$Region == "Wisconsin + UP", Pred.pntsN$WI.TmpEst, NA))
#Get Forest (same process as above, overwriting a few lables)
Comb.plt.df = as.data.frame(JModel.2Reg$summary.random$Forest)[,1:6]
names(Comb.plt.df) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975")
Comb.plt.df$Region = rep(c("Wisconsin + Upper", "Lower"), each = 110)
LW.lu = Comb.plt.df %>% filter(Region == "Lower")
WI.lu = Comb.plt.df %>% filter(Region == "Wisconsin + Upper")
LW.lu$LW.POS = 1:dim(LW.lu)[1]
WI.lu$WI.POS = 1:dim(WI.lu)[1]
Pred.pntsN$LW.temp = sapply(Pred.pntsN$Forest1kmE, function(x)which.min(abs(x - LW.lu$ID)))
Pred.pntsN$LW.TmpEst = as.numeric(with(LW.lu,
Mean[match(Pred.pntsN$LW.temp,
LW.POS)]))
Pred.pntsN$WI.temp = sapply(Pred.pntsN$Forest1kmE, function(x)which.min(abs(x - WI.lu$ID)))
Pred.pntsN$WI.TmpEst = as.numeric(with(WI.lu,
Mean[match(Pred.pntsN$WI.temp,
WI.POS)]))
#Match to region
Pred.pntsN$Forest.Full = ifelse(Pred.pntsN$Region == "Lower", Pred.pntsN$LW.TmpEst,
ifelse(Pred.pntsN$Region == "Wisconsin + UP", Pred.pntsN$WI.TmpEst, NA))
Pred.pntsN$Temp.Full.low = Pred.pntsN$Temp.Full - 0.2766
Pred.pntsN$Temp.Full.high = Pred.pntsN$Temp.Full + 0.2766
Pred.pntsN$Snow.Full.low = Pred.pntsN$Snow.Full - 0.6193
Pred.pntsN$Snow.Full.high = Pred.pntsN$Snow.Full + 0.6193
Pred.pntsN$Forest.Full.low = Pred.pntsN$Forest.Full - 0.6945
Pred.pntsN$Forest.Full.high = Pred.pntsN$Forest.Full + 0.6945
Pred.pntsN@data = Pred.pntsN@data %>%
mutate(Pred.Full = plogis(Temp.Full + Snow.Full + Forest.Full))
Pred.pntsN@data = Pred.pntsN@data %>%
mutate(LPred.Full = plogis(Temp.Full.low + Snow.Full.low + Forest.Full.low))
Pred.pntsN@data = Pred.pntsN@data %>%
mutate(HPred.Full = plogis(Temp.Full.high + Snow.Full.high + Forest.Full.high))
Reg2.pred.r = rasterize(spTransform(
Pred.pntsN,
CRS(proj4string(Domain.r))),
Domain.r,
"Pred.Full",
background = NA)
Reg2.Lpred.r = rasterize(spTransform(
Pred.pntsN,
CRS(proj4string(Domain.r))),
Domain.r,
"LPred.Full",
background = NA)
Reg2.Hpred.r = rasterize(spTransform(
Pred.pntsN,
CRS(proj4string(Domain.r))),
Domain.r,
"HPred.Full",
background = NA)Pred.stk = stack(Reg2.Lpred.r, Reg2.pred.r, Reg2.Hpred.r)
names(Pred.stk) = c("Low", "Mean", "High")
rng = seq(0, 1, 0.01)
mCols = brewer.pal(9, "YlOrBr")
cr = colorRampPalette(mCols)(n = 500)
cr = colorRampPalette(cr,
bias = 0.4, space = "rgb")
Pred.stk2 = Pred.stk
Pred.stk2[is.na(Pred.stk2)] = 0
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.2)
levelplot(Pred.stk2, #Trial.comp,
layout = c(3,1),
margin = FALSE,
xlab = NULL,
ylab = NULL,
#main = "Occurrence Probabilty",
names.attr= c("A", "B", "C"),
maxpixels = 1e5,
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))})
9 Validation
Out of Sample prediction for validation
#Two Region
ModResult = JModel.2Reg
Obs.set = Validation.set2
#Update Region ID for prediction locations (MI and WI only)
MI_WI_Domain = subset(Domain, Name == "Michigan" | Name == "Wisconsin")
MI_WI_DomainP = spTransform(MI_WI_Domain, proj4string(Obs.set))
Obs.set$nDom = over(Obs.set, MI_WI_DomainP)[,1]
Obs.setN = subset(Obs.set, is.na(nDom) == FALSE)
Obs.setN$Region = ifelse(Obs.setN@data$nDom == 21 & Obs.setN@data$UP == "FALSE", "Upper", #FALSE indicates UP
ifelse(Obs.setN@data$nDom == 21 & Obs.setN@data$UP == "TRUE", "Lower",
ifelse(Obs.setN@data$nDom == 48, "Wisconsin", NA)))
Obs.setN$Region = ifelse(Obs.setN$Region == "Upper" | Obs.setN$Region == "Wisconsin", "Wisconsin + UP", Obs.setN$Region)
#Get RF
pLoc = cbind(Obs.setN@coords[,1], Obs.setN@coords[,2])
Ap = inla.spde.make.A(mesh, loc=pLoc)
Obs.setN$REg2.rf = drop(Ap %*% ModResult$summary.random$field.env$mean)
#Get Temperature
Comb.plt.df = as.data.frame(JModel.2Reg$summary.random$mMxTemp)[,1:6] #Get effect estimates (Full model)
names(Comb.plt.df) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975") #Clean labels
Comb.plt.df$Region = rep(c("Wisconsin + Upper", "Lower"), each = 3502)
#Region-sepcific estimates
LW.lu = Comb.plt.df %>% filter(Region == "Lower")
WI.lu = Comb.plt.df %>% filter(Region == "Wisconsin + Upper")
LW.lu$LW.POS = 1:dim(LW.lu)[1]
WI.lu$WI.POS = 1:dim(WI.lu)[1]
#Look up Lower Temp Estimates
Obs.setN$LW.temp = sapply(Obs.setN$mxTempE, function(x)which.min(abs(x - LW.lu$ID)))
Obs.setN$LW.TmpEst = as.numeric(with(LW.lu,
Mean[match(Obs.setN$LW.temp,
LW.POS)]))
#Look up WI + UP Temp Estimates
Obs.setN$WI.temp = sapply(Obs.setN$mxTempE, function(x)which.min(abs(x - WI.lu$ID)))
Obs.setN$WI.TmpEst = as.numeric(with(WI.lu,
Mean[match(Obs.setN$WI.temp,
WI.POS)]))
#Match to region
Obs.setN$Temp.Full = ifelse(Obs.setN$Region == "Wisconsin + UP", Obs.setN$WI.TmpEst,
ifelse(Obs.setN$Region == "Lower", Obs.setN$LW.TmpEst, NA))
#Get Snow Week (same process as above, overwriting a few lables)
Comb.plt.df = as.data.frame(JModel.2Reg$summary.random$mSnw5yr)[,1:6]
names(Comb.plt.df) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975")
Comb.plt.df$Region = rep(c("Wisconsin + Upper", "Lower"), each = 90)
LW.lu = Comb.plt.df %>% filter(Region == "Lower")
WI.lu = Comb.plt.df %>% filter(Region == "Wisconsin + Upper")
LW.lu$LW.POS = 1:dim(LW.lu)[1]
WI.lu$WI.POS = 1:dim(WI.lu)[1]
Obs.setN$LW.temp = sapply(Obs.setN$mSnow5yrE, function(x)which.min(abs(x - LW.lu$ID)))
Obs.setN$LW.TmpEst = as.numeric(with(LW.lu,
Mean[match(Obs.setN$LW.temp,
LW.POS)]))
Obs.setN$WI.temp = sapply(Obs.setN$mSnow5yrE, function(x)which.min(abs(x - WI.lu$ID)))
Obs.setN$WI.TmpEst = as.numeric(with(WI.lu,
Mean[match(Obs.setN$WI.temp,
WI.POS)]))
#Match to region
Obs.setN$Snow.Full = ifelse(Obs.setN$Region == "Lower", Obs.setN$LW.TmpEst,
ifelse(Obs.setN$Region == "Wisconsin + UP", Obs.setN$WI.TmpEst, NA))
#Get Forest (same process as above, overwriting a few lables)
Comb.plt.df = as.data.frame(JModel.2Reg$summary.random$Forest)[,1:6]
names(Comb.plt.df) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975")
Comb.plt.df$Region = rep(c("Wisconsin + Upper", "Lower"), each = 110)
LW.lu = Comb.plt.df %>% filter(Region == "Lower")
WI.lu = Comb.plt.df %>% filter(Region == "Wisconsin + Upper")
LW.lu$LW.POS = 1:dim(LW.lu)[1]
WI.lu$WI.POS = 1:dim(WI.lu)[1]
Obs.setN$LW.temp = sapply(Obs.setN$Forest1kmE, function(x)which.min(abs(x - LW.lu$ID)))
Obs.setN$LW.TmpEst = as.numeric(with(LW.lu,
Mean[match(Obs.setN$LW.temp,
LW.POS)]))
Obs.setN$WI.temp = sapply(Obs.setN$Forest1kmE, function(x)which.min(abs(x - WI.lu$ID)))
Obs.setN$WI.TmpEst = as.numeric(with(WI.lu,
Mean[match(Obs.setN$WI.temp,
WI.POS)]))
#Match to region
Obs.setN$Forest.Full = ifelse(Obs.setN$Region == "Lower", Obs.setN$LW.TmpEst,
ifelse(Obs.setN$Region == "Wisconsin + UP", Obs.setN$WI.TmpEst, NA))
Obs.setN@data = Obs.setN@data %>%
mutate(Pred.Full = plogis(Temp.Full + Snow.Full + Forest.Full))
Obs.set.R2.keep = Obs.setN
Reg2.PA = as.data.frame(cbind(1:dim(Obs.setN)[1], Obs.setN$OBS, Obs.setN$Pred.Full))
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.95
## 2 PredPrev=Obs 0.97Reg2.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
OutofSample = Reg2.PA.outIn Sample prediction for validation
#Two Region
ModResult = JModel.2Reg
Obs.set = Keep.valid.df
idat = inla.stack.index(Joint2.stk, "jenv.0")$data
Obs.set$Est = JModel.2Reg$summary.fitted.values$mean[idat]
Obs.set = Obs.set %>% filter(Data == "Fit")
Obs.set = SpatialPointsDataFrame(Obs.set[,c("Long","Lat")], Obs.set)
proj4string(Obs.set) = proj4string(Validation.set2)
#Update Region ID for prediction locations (MI and WI only)
MI_WI_Domain = subset(Domain, Name == "Michigan" | Name == "Wisconsin")
MI_WI_DomainP = spTransform(MI_WI_Domain, proj4string(Obs.set))
Obs.set$nDom = over(Obs.set, MI_WI_DomainP)[,1]
Obs.setN = subset(Obs.set, is.na(nDom) == FALSE)
Obs.setN$Region = ifelse(Obs.setN@data$nDom == 21 & Obs.setN@data$UP == "FALSE", "Upper", #FALSE indicates UP
ifelse(Obs.setN@data$nDom == 21 & Obs.setN@data$UP == "TRUE", "Lower",
ifelse(Obs.setN@data$nDom == 48, "Wisconsin", NA)))
Obs.setN$Region = ifelse(Obs.setN$Region == "Upper" | Obs.setN$Region == "Wisconsin", "Wisconsin + UP", Obs.setN$Region)
#Get RF
pLoc = cbind(Obs.setN@coords[,1], Obs.setN@coords[,2])
Ap = inla.spde.make.A(mesh, loc=pLoc)
Obs.setN$REg2.rf = drop(Ap %*% ModResult$summary.random$field.env$mean)
#Get Temperature
Comb.plt.df = as.data.frame(JModel.2Reg$summary.random$mMxTemp)[,1:6] #Get effect estimates (Full model)
names(Comb.plt.df) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975") #Clean labels
Comb.plt.df$Region = rep(c("Wisconsin + Upper", "Lower"), each = 3502)
#Region-sepcific estimates
LW.lu = Comb.plt.df %>% filter(Region == "Lower")
WI.lu = Comb.plt.df %>% filter(Region == "Wisconsin + Upper")
LW.lu$LW.POS = 1:dim(LW.lu)[1]
WI.lu$WI.POS = 1:dim(WI.lu)[1]
#Look up Lower Temp Estimates
Obs.setN$LW.temp = sapply(Obs.setN$mxTempE, function(x)which.min(abs(x - LW.lu$ID)))
Obs.setN$LW.TmpEst = as.numeric(with(LW.lu,
Mean[match(Obs.setN$LW.temp,
LW.POS)]))
#Look up WI + UP Temp Estimates
Obs.setN$WI.temp = sapply(Obs.setN$mxTempE, function(x)which.min(abs(x - WI.lu$ID)))
Obs.setN$WI.TmpEst = as.numeric(with(WI.lu,
Mean[match(Obs.setN$WI.temp,
WI.POS)]))
#Match to region
Obs.setN$Temp.Full = ifelse(Obs.setN$Region == "Wisconsin + UP", Obs.setN$WI.TmpEst,
ifelse(Obs.setN$Region == "Lower", Obs.setN$LW.TmpEst, NA))
#Get Snow Week (same process as above, overwriting a few lables)
Comb.plt.df = as.data.frame(JModel.2Reg$summary.random$mSnw5yr)[,1:6]
names(Comb.plt.df) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975")
Comb.plt.df$Region = rep(c("Wisconsin + Upper", "Lower"), each = 90)
LW.lu = Comb.plt.df %>% filter(Region == "Lower")
WI.lu = Comb.plt.df %>% filter(Region == "Wisconsin + Upper")
LW.lu$LW.POS = 1:dim(LW.lu)[1]
WI.lu$WI.POS = 1:dim(WI.lu)[1]
Obs.setN$LW.temp = sapply(Obs.setN$mSnow5yrE, function(x)which.min(abs(x - LW.lu$ID)))
Obs.setN$LW.TmpEst = as.numeric(with(LW.lu,
Mean[match(Obs.setN$LW.temp,
LW.POS)]))
Obs.setN$WI.temp = sapply(Obs.setN$mSnow5yrE, function(x)which.min(abs(x - WI.lu$ID)))
Obs.setN$WI.TmpEst = as.numeric(with(WI.lu,
Mean[match(Obs.setN$WI.temp,
WI.POS)]))
#Match to region
Obs.setN$Snow.Full = ifelse(Obs.setN$Region == "Lower", Obs.setN$LW.TmpEst,
ifelse(Obs.setN$Region == "Wisconsin + UP", Obs.setN$WI.TmpEst, NA))
#Get Forest (same process as above, overwriting a few lables)
Comb.plt.df = as.data.frame(JModel.2Reg$summary.random$Forest)[,1:6]
names(Comb.plt.df) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975")
Comb.plt.df$Region = rep(c("Wisconsin + Upper", "Lower"), each = 110)
LW.lu = Comb.plt.df %>% filter(Region == "Lower")
WI.lu = Comb.plt.df %>% filter(Region == "Wisconsin + Upper")
LW.lu$LW.POS = 1:dim(LW.lu)[1]
WI.lu$WI.POS = 1:dim(WI.lu)[1]
Obs.setN$LW.temp = sapply(Obs.setN$Forest1kmE, function(x)which.min(abs(x - LW.lu$ID)))
Obs.setN$LW.TmpEst = as.numeric(with(LW.lu,
Mean[match(Obs.setN$LW.temp,
LW.POS)]))
Obs.setN$WI.temp = sapply(Obs.setN$Forest1kmE, function(x)which.min(abs(x - WI.lu$ID)))
Obs.setN$WI.TmpEst = as.numeric(with(WI.lu,
Mean[match(Obs.setN$WI.temp,
WI.POS)]))
#Match to region
Obs.setN$Forest.Full = ifelse(Obs.setN$Region == "Lower", Obs.setN$LW.TmpEst,
ifelse(Obs.setN$Region == "Wisconsin + UP", Obs.setN$WI.TmpEst, NA))
Obs.setN@data = Obs.setN@data %>%
mutate(Pred.Full = plogis(Temp.Full + Snow.Full + Forest.Full))
Obs.set.R2.keep = Obs.setN
Reg2.PA = as.data.frame(cbind(1:dim(Obs.setN)[1], Obs.setN$OBS, Obs.setN$Pred.Full))
names(Reg2.PA) = c("ID", "OBS", "Pred")
Thresh = optimal.thresholds(Reg2.PA, opt.methods = c("MaxSens+Spec", "PredPrev=Obs"))
Reg2.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
InSample = Reg2.PA.outValidation Table
InSample$Set = "In Sample"
OutofSample$Set = "Out of Sample"
Valid.table = rbind(InSample, OutofSample)
Valid.table = Valid.table %>%
select(Set, threshold, PCC, sensitivity, specificity, AUC, TSS)
names(Valid.table) = c("Type", "Threshold", "PCC", "Sensitivity", "Specificity", "AUC", "TSS")
kable(Valid.table, caption = "Model Validation") %>%
kable_styling("striped", full_width = F) %>%
row_spec(0, font_size = 20) %>%
column_spec(1, bold = T)| Type | Threshold | PCC | Sensitivity | Specificity | AUC | TSS |
|---|---|---|---|---|---|---|
| In Sample | 0.94 | 0.6527861 | 0.9147287 | 0.6422562 | 0.8676953 | 0.5569848 |
| In Sample | 0.50 | 0.3511084 | 0.9922481 | 0.3253350 | 0.8676953 | 0.3175831 |
| Out of Sample | 0.95 | 0.7674476 | 0.9717349 | 0.6982502 | 0.8218957 | 0.6699851 |
| Out of Sample | 0.50 | 0.5094945 | 1.0000000 | 0.3433476 | 0.8218957 | 0.3433476 |
#SelMod.fx = xtable(Valid.table)
#print(SelMod.fx, include.rownames = F)