Marten Scale (02-15-20)
1 Prepare Data
Set WD
#setwd("C:/Users/roloff/OneDrive - Michigan State University/Sault Tribe Projects/Sault Tribe Marten/GrayAnalysis2019/NewFeb2020/RRunFeb6")
setwd("E:/Marten_2")Load Libraries
Needed.packages = c("rgdal", "raster", "corrplot", "spdep",
"rgeos","maptools","mapproj","GISTools",
"sp","ggplot2","ggthemes","plyr","dplyr",
"kableExtra", "INLABMA", "sf", "reshape2")
for(p in Needed.packages){
if(!require(p, character.only = TRUE)) install.packages(p, repos = "http://cran.us.r-project.org")
suppressMessages(library(p, character.only = TRUE))
}## Loading required package: rgdal## Loading required package: sp## rgdal: version: 1.4-8, (SVN revision 845)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 2.2.3, released 2017/11/20
## Path to GDAL shared files: C:/Program Files/R/R-3.5.3/library/rgdal/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
## Path to PROJ.4 shared files: C:/Program Files/R/R-3.5.3/library/rgdal/proj
## Linking to sp version: 1.3-2## Loading required package: raster## Loading required package: corrplot## corrplot 0.84 loaded## Loading required package: spdep## Loading required package: spData## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`## Loading required package: sf## Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3## Loading required package: rgeos## rgeos version: 0.5-2, (SVN revision 621)
## GEOS runtime version: 3.6.1-CAPI-1.10.1
## Linking to sp version: 1.3-1
## Polygon checking: TRUE## Loading required package: maptools## Checking rgeos availability: TRUE## Loading required package: mapproj## Loading required package: maps## Loading required package: GISTools## Loading required package: RColorBrewer## Loading required package: MASS##
## Attaching package: 'MASS'## The following objects are masked from 'package:raster':
##
## area, select##
## Attaching package: 'GISTools'## The following object is masked from 'package:maps':
##
## map.scale## Loading required package: ggplot2## Loading required package: ggthemes## Loading required package: plyr##
## Attaching package: 'plyr'## The following object is masked from 'package:maps':
##
## ozone## Loading required package: dplyr##
## Attaching package: 'dplyr'## The following objects are masked from 'package:plyr':
##
## arrange, count, desc, failwith, id, mutate, rename, summarise,
## summarize## The following object is masked from 'package:MASS':
##
## select## The following objects are masked from 'package:rgeos':
##
## intersect, setdiff, union## The following objects are masked from 'package:raster':
##
## intersect, select, union## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union## Loading required package: kableExtra##
## Attaching package: 'kableExtra'## The following object is masked from 'package:dplyr':
##
## group_rows## Loading required package: INLABMA## Loading required package: parallel## Loading required package: reshape2#INLA not on Cran
if(!require("INLA", character.only = TRUE)) install.packages("INLA", repos=c(getOption("repos"), INLA="https://inla.r-inla-download.org/R/stable"), dep=TRUE)## Loading required package: INLA## Loading required package: Matrix## This is INLA_18.07.12 built 2018-07-12 11:05:18 UTC.
## See www.r-inla.org/contact-us for how to get help.suppressMessages(library("INLA", character.only = TRUE))Load Data
#martenscl <- read.csv("Marten_Scale.csv", header=T)
martenscl = read.csv("MartenDataNew.csv", header=T)
head(martenscl)## FID Field1 Duration DOP Satellites id Longitude Latitude RASTERVALU
## 1 0 11007 26 2.8 4 7004m 676032.0 5141529 0.6205056
## 2 1 21003 2 2.6 4 7004m 675910.1 5141611 0.6644616
## 3 2 3569 2 2.4 4 7004m 675943.2 5141910 0.5873879
## 4 3 4108 9 3.4 4 7004m 675836.6 5141944 0.5672649
## 5 4 5108 15 6.0 4 7004m 675631.0 5142643 0.3531531
## 6 5 6108 4 9.0 4 7004m 675811.5 5143019 0.4251542
## hardfor90 confor90 hardtall90 contall90 hardsh90 consh90 herbopen90
## 1 0.03125 0.59375 0.03125 0.28125 0 0.31250 0
## 2 0.03125 0.75000 0.03125 0.75000 0 0.00000 0
## 3 0.03125 0.81250 0.03125 0.71875 0 0.09375 0
## 4 0.25000 0.40625 0.25000 0.40625 0 0.00000 0
## 5 0.09375 0.81250 0.09375 0.81250 0 0.00000 0
## 6 0.00000 0.78125 0.00000 0.78125 0 0.00000 0
## riparian90 hardfor120 confor120 hardtal120 contall120 hardsh120 consh120
## 1 0 0.01923077 0.7692308 0.01923077 0.4615385 0 0.30769231
## 2 0 0.01923077 0.6346154 0.01923077 0.6346154 0 0.00000000
## 3 0 0.07692308 0.8269231 0.07692308 0.7692308 0 0.05769231
## 4 0 0.21153846 0.3461538 0.21153846 0.3461538 0 0.00000000
## 5 0 0.13461539 0.8269231 0.13461539 0.8269231 0 0.00000000
## 6 0 0.00000000 0.8461538 0.00000000 0.8461538 0 0.00000000
## herbop120 ripar120 hardfor150 confor150 hardtal150 contall150 hardsh150
## 1 0 0 0.0125 0.8250 0.0125 0.5375 0
## 2 0 0 0.0625 0.4750 0.0625 0.4625 0
## 3 0 0 0.1375 0.8000 0.1375 0.7625 0
## 4 0 0 0.1000 0.3750 0.1000 0.3750 0
## 5 0 0 0.1000 0.9000 0.1000 0.9000 0
## 6 0 0 0.0000 0.8375 0.0000 0.8250 0
## consh150 herbop150 rip150 Paved_m Unpaved_m hydro_m
## 1 0.2875 0.0000 0 1488 943 3377
## 2 0.0125 0.0000 0 1480 1061 3479
## 3 0.0375 0.0000 0 1184 1015 3767
## 4 0.0000 0.0000 0 1145 1120 3820
## 5 0.0000 0.0000 0 436 1294 3970
## 6 0.0125 0.0125 0 70 1098 3864length(unique(martenscl$id))## [1] 13unique(martenscl$id)## [1] 7004m 6977m 5516m 1114m 5509m 5529m 5517m 8217m houdini
## [10] 5538m 5064f 5496m 2254m
## 13 Levels: 1114m 2254m 5064f 5496m 5509m 5516m 5517m 5529m 5538m ... houdiniExtract month and year from GMT_Time field
# Julian date variable
mart_time <- read.csv("mart_TimeFormat.csv", header=T)
colnames(mart_time)[1] <- "Field1"
marten <- merge(martenscl, mart_time, by=c("Field1"))
marten$moda <- as.character(substr(marten$GMT.Time,1,5))
marten$moda <- as.Date(marten$moda, format="%m/%d")
marten$jdate <- as.numeric(strftime(marten$moda, format="%j"))
# Year for random effect
marten$year <- as.character(substr(marten$GMT.Time,7,10))
unique(marten$year)## [1] "2016" "2017"Scale and Center environmental variables (columns 10:36)
Cov.scale = marten[,c(10:36)]
head(Cov.scale)## hardfor90 confor90 hardtall90 contall90 hardsh90 consh90 herbopen90
## 1 0.28125 0.09375 0.28125 0.09375 0 0 0
## 2 0.06250 0.06250 0.06250 0.06250 0 0 0
## 3 0.00000 0.00000 0.00000 0.00000 0 0 0
## 4 0.00000 0.00000 0.00000 0.00000 0 0 0
## 5 0.03125 0.00000 0.03125 0.00000 0 0 0
## 6 0.65625 0.25000 0.65625 0.25000 0 0 0
## riparian90 hardfor120 confor120 hardtal120 contall120 hardsh120 consh120
## 1 0.00000 0.26923077 0.07692308 0.26923077 0.07692308 0 0
## 2 0.12500 0.03846154 0.03846154 0.03846154 0.03846154 0 0
## 3 0.00000 0.00000000 0.00000000 0.00000000 0.00000000 0 0
## 4 0.00000 0.00000000 0.00000000 0.00000000 0.00000000 0 0
## 5 0.00000 0.01923077 0.00000000 0.01923077 0.00000000 0 0
## 6 0.03125 0.61538461 0.32692308 0.61538461 0.32692308 0 0
## herbop120 ripar120 hardfor150 confor150 hardtal150 contall150 hardsh150
## 1 0 0.03846154 0.2250 0.0750 0.2250 0.0750 0
## 2 0 0.03846154 0.0125 0.0125 0.0125 0.0125 0
## 3 0 0.00000000 0.0000 0.0000 0.0000 0.0000 0
## 4 0 0.00000000 0.0000 0.0000 0.0000 0.0000 0
## 5 0 0.00000000 0.0000 0.0000 0.0000 0.0000 0
## 6 0 0.05769231 0.5000 0.3500 0.5000 0.3500 0
## consh150 herbop150 rip150 Paved_m Unpaved_m hydro_m
## 1 0 0 0.025 10 185 553
## 2 0 0 0.025 70 130 485
## 3 0 0 0.000 136 59 428
## 4 0 0 0.000 136 60 427
## 5 0 0 0.000 134 59 437
## 6 0 0 0.125 145 339 870for(i in 1:ncol(Cov.scale)){
Cov.scale[,i] = as.numeric(scale(Cov.scale[,i], scale = T, center = T))
}
names(Cov.scale) = paste("s_", names(Cov.scale), sep="")
#Cov.scale$s_jdate = as.numeric(scale(marten$jdate, scale = T, center = T))
#Add scaled version to data frame with an "s" in front of column name
Mod.cov.df = cbind(marten, Cov.scale)
head(Mod.cov.df)## Field1 FID Duration.x DOP.x Satellites.x id.x Longitude.x Latitude.x
## 1 18 7987 59 2.4 4 2254m 675505.0 5136966
## 2 19 7988 6 1.8 4 2254m 675562.8 5137027
## 3 20 7989 3 1.8 4 2254m 675631.1 5136975
## 4 21 7990 2 2.2 4 2254m 675630.5 5136985
## 5 22 7991 2 2.8 4 2254m 675630.1 5136948
## 6 23 7992 9 3.4 4 2254m 675362.8 5136545
## RASTERVALU hardfor90 confor90 hardtall90 contall90 hardsh90 consh90
## 1 0.12328015 0.28125 0.09375 0.28125 0.09375 0 0
## 2 0.19138694 0.06250 0.06250 0.06250 0.06250 0 0
## 3 0.21412432 0.00000 0.00000 0.00000 0.00000 0 0
## 4 0.21412432 0.00000 0.00000 0.00000 0.00000 0 0
## 5 0.21412432 0.03125 0.00000 0.03125 0.00000 0 0
## 6 0.02193207 0.65625 0.25000 0.65625 0.25000 0 0
## herbopen90 riparian90 hardfor120 confor120 hardtal120 contall120 hardsh120
## 1 0 0.00000 0.26923077 0.07692308 0.26923077 0.07692308 0
## 2 0 0.12500 0.03846154 0.03846154 0.03846154 0.03846154 0
## 3 0 0.00000 0.00000000 0.00000000 0.00000000 0.00000000 0
## 4 0 0.00000 0.00000000 0.00000000 0.00000000 0.00000000 0
## 5 0 0.00000 0.01923077 0.00000000 0.01923077 0.00000000 0
## 6 0 0.03125 0.61538461 0.32692308 0.61538461 0.32692308 0
## consh120 herbop120 ripar120 hardfor150 confor150 hardtal150 contall150
## 1 0 0 0.03846154 0.2250 0.0750 0.2250 0.0750
## 2 0 0 0.03846154 0.0125 0.0125 0.0125 0.0125
## 3 0 0 0.00000000 0.0000 0.0000 0.0000 0.0000
## 4 0 0 0.00000000 0.0000 0.0000 0.0000 0.0000
## 5 0 0 0.00000000 0.0000 0.0000 0.0000 0.0000
## 6 0 0 0.05769231 0.5000 0.3500 0.5000 0.3500
## hardsh150 consh150 herbop150 rip150 Paved_m Unpaved_m hydro_m
## 1 0 0 0 0.025 10 185 553
## 2 0 0 0 0.025 70 130 485
## 3 0 0 0 0.000 136 59 428
## 4 0 0 0 0.000 136 60 427
## 5 0 0 0 0.000 134 59 437
## 6 0 0 0 0.125 145 339 870
## GMT.Time Duration.y DOP.y Satellites.y id.y Longitude.y Latitude.y
## 1 03/05/2016 19:30 59 2.4 4 2254m 675505.0 5136966
## 2 03/05/2016 19:45 6 1.8 4 2254m 675562.8 5137027
## 3 03/05/2016 20:00 3 1.8 4 2254m 675631.1 5136975
## 4 03/05/2016 20:15 2 2.2 4 2254m 675630.5 5136985
## 5 03/05/2016 20:30 2 2.8 4 2254m 675630.1 5136948
## 6 03/05/2016 20:45 9 3.4 4 2254m 675362.8 5136545
## moda jdate year s_hardfor90 s_confor90 s_hardtall90 s_contall90
## 1 2020-03-05 65 2016 1.45620350 0.06917282 1.52546893 0.14368603
## 2 2020-03-05 65 2016 0.03757526 -0.09513455 0.06185808 -0.03794139
## 3 2020-03-05 65 2016 -0.36774710 -0.42374930 -0.35631645 -0.40119623
## 4 2020-03-05 65 2016 -0.36774710 -0.42374930 -0.35631645 -0.40119623
## 5 2020-03-05 65 2016 -0.16508592 -0.42374930 -0.14722918 -0.40119623
## 6 2020-03-05 65 2016 3.88813763 0.89070969 4.03451610 1.05182313
## s_hardsh90 s_consh90 s_herbopen90 s_riparian90 s_hardfor120 s_confor120
## 1 -0.1089675 -0.1766424 -0.1242558 -0.26375740 1.61567196 0.05854051
## 2 -0.1089675 -0.1766424 -0.1242558 0.87941736 -0.05509775 -0.16030494
## 3 -0.1089675 -0.1766424 -0.1242558 -0.26375740 -0.33355936 -0.37915037
## 4 -0.1089675 -0.1766424 -0.1242558 -0.26375740 -0.33355936 -0.37915037
## 5 -0.1089675 -0.1766424 -0.1242558 -0.26375740 -0.19432855 -0.37915037
## 6 -0.1089675 -0.1766424 -0.1242558 0.02203629 4.12182651 1.48103586
## s_hardtal120 s_contall120 s_hardsh120 s_consh120 s_herbop120 s_ripar120
## 1 1.68966047 0.1294811 -0.1063301 -0.1682362 -0.1259219 0.1410059
## 2 -0.03512312 -0.1158332 -0.1063301 -0.1682362 -0.1259219 0.1410059
## 3 -0.32258704 -0.3611475 -0.1063301 -0.1682362 -0.1259219 -0.2423773
## 4 -0.32258704 -0.3611475 -0.1063301 -0.1682362 -0.1259219 -0.2423773
## 5 -0.17885508 -0.3611475 -0.1063301 -0.1682362 -0.1259219 -0.2423773
## 6 4.27683585 1.7240241 -0.1063301 -0.1682362 -0.1259219 0.3326975
## s_hardfor150 s_confor150 s_hardtal150 s_contall150 s_hardsh150 s_consh150
## 1 1.6021914 0.1460311 1.6761332 0.2260010 -0.1007484 -0.1587382
## 2 -0.1980222 -0.2623019 -0.1842321 -0.2369178 -0.1007484 -0.1587382
## 3 -0.3039172 -0.3439685 -0.2936654 -0.3295015 -0.1007484 -0.1587382
## 4 -0.3039172 -0.3439685 -0.2936654 -0.3295015 -0.1007484 -0.1587382
## 5 -0.3039172 -0.3439685 -0.2936654 -0.3295015 -0.1007484 -0.1587382
## 6 3.9318797 1.9426962 4.0836649 2.2628437 -0.1007484 -0.1587382
## s_herbop150 s_rip150 s_Paved_m s_Unpaved_m s_hydro_m
## 1 -0.1244583 0.05927052 -1.1602269 -1.190281 -0.31487316
## 2 -0.1244583 0.05927052 -1.0151395 -1.202578 -0.39481041
## 3 -0.1244583 -0.22632565 -0.8555433 -1.218453 -0.46181662
## 4 -0.1244583 -0.22632565 -0.8555433 -1.218229 -0.46299217
## 5 -0.1244583 -0.22632565 -0.8603795 -1.218453 -0.45123669
## 6 -0.1244583 1.20165522 -0.8337802 -1.155849 0.05777545View Correlation Table
Cov.cor = cor(Cov.scale)
corrplot(Cov.cor,
tl.cex=0.7,
number.cex = 0.6,
method = "number")
Convert to Spatial Points
nProj = "+proj=utm +zone=16 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"
Marten.pnts = SpatialPointsDataFrame(Mod.cov.df[,c("Longitude.x", "Latitude.x")], Mod.cov.df)
proj4string(Marten.pnts) = nProjLoad County Boundaries
Counties = readOGR(dsn = "./Counties_v17a",
layer = "Counties_v17a",
stringsAsFactors = FALSE)## OGR data source with driver: ESRI Shapefile
## Source: "E:\Marten_2\Counties_v17a", layer: "Counties_v17a"
## with 83 features
## It has 15 fields
## Integer64 fields read as strings: OBJECTID FIPSNUMCounties = spTransform(Counties, proj4string(Marten.pnts))Plot Points
wmap_df = fortify(Counties)## Regions defined for each Polygonstmp.df = Marten.pnts@data
ggplot(wmap_df, aes(long,lat, group=group)) +
geom_polygon(fill = "lightgray", col="darkgray") +
geom_point(data=tmp.df,
aes(Longitude.x, Latitude.x,
group=NULL,
fill=NULL,
col = "red")) +
xlab("Longitude") +
ylab("Latitude") +
ggtitle("Marten Locations") +
coord_equal() +
theme(panel.grid.minor = element_blank(),
panel.grid.major = element_blank(),
panel.background = element_blank(),
panel.border = element_blank(),
legend.title = element_blank(),
legend.position = "none",
axis.title.x = element_text(size=22, face="bold"),
axis.title.y = element_text(size=22, face="bold"),
plot.title = element_text(size=24, face="bold", hjust = 0.5))
Calculate Nearest Neighbor Distance by Marten ID
library(spatstat)## Loading required package: spatstat.data## Loading required package: nlme##
## Attaching package: 'nlme'## The following object is masked from 'package:dplyr':
##
## collapse## The following object is masked from 'package:raster':
##
## getData## Loading required package: rpart##
## spatstat 1.63-0 (nickname: 'Space camouflage')
## For an introduction to spatstat, type 'beginner'##
## Attaching package: 'spatstat'## The following object is masked from 'package:MASS':
##
## area## The following objects are masked from 'package:raster':
##
## area, rotate, shiftMart.ids = levels(as.factor(Marten.pnts$id.x))
for(i in 1:length(Mart.ids)){ #for all marten IDs
Mart.tmp = subset(Marten.pnts, id.x == Mart.ids[i]) #choose one marten at a time
P.pp = as.ppp(Mart.tmp) #convert to ppp object
Mart.tmp@data$NN = as.numeric(nndist(P.pp, k = 1)) #find nearest neighboring point
if(i == 1){Marten.df = Mart.tmp@data #re-assemble data frame
} else{Marten.df = rbind(Marten.df, Mart.tmp@data)}
}
#large range
range(Marten.df$NN)## [1] 0.000 1984.528#Scale neighbor distance
Marten.df$NNs = Marten.df$NN/100
#Convert back to points
Marten.pnts = SpatialPointsDataFrame(Marten.df[,c("Longitude.x", "Latitude.x")], Marten.df)
proj4string(Marten.pnts) = nProj
detach("package:spatstat", unload=TRUE)2 Individual Covariate Models
Looping through each covariate one at a time with marten ID as a grouping random effect and nearest neigbor distance (within individual marten) to account for spatial bias.
3 Response Transformation
Logit Transformaion (https://besjournals.onlinelibrary.wiley.com/doi/full/10.1111/2041-210X.13234)
Marten.pnts$t.RASTERVALU = log(Marten.pnts$RASTERVALU/(1-Marten.pnts$RASTERVALU))
plot(density(Marten.pnts$t.RASTERVALU))
Covariates.list = levels(factor(names(Marten.pnts@data[,47:73]))) #Get names of all covariates
Covariates.list
Marten.data = Marten.pnts@data #copy of the data
for(i in 1:length(Covariates.list)){ #for each covariate
Cov.lst = list(list(intercept1 = rep(1, dim(Marten.data)[1])), #intercept = 1
list(Covariate.X = Marten.data[,Covariates.list[i]], #Choose 1 covariate at a time
NN = Marten.data[,"NNs"], #Scaled neighbor distance
ID = Marten.data[,"id.x"], #Marten ID
JDATE = as.factor(as.character(Marten.data[,"jdate"])), #Day
YEAR = as.factor(as.character(Marten.data[,"year"])))) #Year
Data.stk = inla.stack(data = list( #Organizing data as a list
RASTERVALU = Marten.data$t.RASTERVALU), #Response/dependent variable
A = list(1,1), #Not used, empty space
effects = Cov.lst, #Individual covariate, neighbor, ID, Year
tag = "M1.0") #Just a label
nPrior = list(theta=list(prior = "normal", param=c(0, 5))) #flat prior
Formula = RASTERVALU ~ -1 + intercept1 + #response and intercept
f(ID, #grouping effect for individual marten
model="iid",
constr = TRUE, #enforce a zero mean (average)
hyper = nPrior) +
f(inla.group(NN, #distance as a smooth effect
n = 20, #binning distances because some distances are very close
method = "cut"),
model="rw1",
constr = TRUE,
scale.model = TRUE,
hyper = nPrior) +
f(YEAR, #Year as grouping variable
model="iid",
constr = TRUE,
hyper = nPrior) +
f(JDATE, #DOY as grouping variable
model="iid",
constr = TRUE,
hyper = nPrior) +
Covariate.X #Individual covariate of interest
Model.X = inla(Formula, #Run model iteration
data = inla.stack.data(Data.stk),
family = "gaussian",
verbose = TRUE,
control.fixed = list(prec = 0.001,
prec.intercept = 0.0001),
control.predictor = list(
A = inla.stack.A(Data.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.compute=list(dic = TRUE, cpo = TRUE, waic = TRUE))
#idat = inla.stack.index(Data.stk, "M1.0")$data #Checking estimates
#Marten.data$Fit.v = Model.X$summary.fitted.values$mean[idat] #Fitted values
#cor(Marten.data$RASTERVALU, plogis(Marten.data$Fit.v))
#Create a data frame to store results
Fixed.effects = Model.X$summary.fixed[,c(1:3,5)]
names(Fixed.effects) = c("Mean", "sd", "Q.025", "Q.975")
#Add comparison metrics DIC and WAIC
Fixed.effects$DIC = Model.X$dic$dic
Fixed.effects$WAIC = Model.X$waic$waic
#Add label to identify actual covariate name
Fixed.effects$Covariate = Covariates.list[i]
#Get estimates for ID effect
ID.effect = Model.X$summary.random$ID[,c(1:4,6)]
names(ID.effect) = c("ID", "Mean", "sd", "Q.025", "Q.975")
ID.effect$Covariate = Covariates.list[i]
#Get estimates for NN distance effect
NN.effect = Model.X$summary.random$`inla.group(NN, n = 20, method = "cut")`[,c(1:4,6)]
names(NN.effect) = c("Distance", "Mean", "sd", "Q.025", "Q.975")
NN.effect$Covariate = Covariates.list[i]
#Get estimates for JDATE effect
JDate.effect = Model.X$summary.random$JDATE[,c(1:4,6)]
names(JDate.effect) = c("JDATE", "Mean", "sd", "Q.025", "Q.975")
JDate.effect$Covariate = Covariates.list[i]
#Get estimates for YEAR effect
Year.effect = Model.X$summary.random$YEAR[,c(1:4,6)]
names(Year.effect) = c("YEAR", "Mean", "sd", "Q.025", "Q.975")
Year.effect$Covariate = Covariates.list[i]
#Combine data for all iterations
if(i == 1){
Fixed.effects.all = Fixed.effects
ID.effect.all = ID.effect
NN.effect.all = NN.effect
JDate.effect.all = JDate.effect
Year.effect.all = Year.effect
} else{
Fixed.effects.all = rbind(Fixed.effects.all, Fixed.effects)
ID.effect.all = rbind(ID.effect.all, ID.effect)
NN.effect.all = rbind(NN.effect.all, NN.effect)
JDate.effect.all = rbind(JDate.effect.all,JDate.effect)
Year.effect.all = rbind(Year.effect.all, Year.effect)}
}
#Determine which are signficant as judged by credible interval
Fixed.effects.all$Credible = ifelse(Fixed.effects.all$Q.025 >= 0 & Fixed.effects.all$Q.975 >= 0, "Inside",
ifelse(Fixed.effects.all$Q.025 <= 0 & Fixed.effects.all$Q.975 <= 0, "Inside",
"Outside"))
#Are row values the covariate or intercept?
Fixed.effects.all$Type = paste(substr(rownames(Fixed.effects.all), 1, 3))Save copies to WD
write.csv(Fixed.effects.all, "./Ind_fixed_effects.csv") #fixed effects
write.csv(ID.effect.all, "./Ind_ID_effect.csv")
write.csv(NN.effect.all, "./Ind_Distance_effect.csv")
write.csv(Year.effect.all, "./Year_effect.csv")Covariate Selection
#Remove those "outside" credible interval and intercepts
Fixed.effects.all2 = Fixed.effects.all %>% filter(Credible == "Inside" & Type == "Cov")
#How many dropped
dim(Fixed.effects.all)[1] - dim(Fixed.effects.all2)[1]## [1] 30Fixed.effects.all2$Var = gsub("\\d", "", Fixed.effects.all2$Covariate)
Fixed.effects.all2$Var = gsub("s_Area_", "", Fixed.effects.all2$Var)
Fixed.effects.all2 = arrange(Fixed.effects.all2, (DIC))
Fixed.effects.all2$Dup = duplicated(Fixed.effects.all2[,"Var"])
Fixed.effects.all3 = Fixed.effects.all2 %>%
filter(Dup != TRUE)
#Top of each type
Fixed.effects.all3## Mean sd Q.025 Q.975 DIC WAIC
## 1 0.20667635 0.01915107 0.1690763527 0.24424497 29797.50 29803.07
## 2 -0.48824294 0.04538085 -0.5773408448 -0.39921938 29797.52 29802.69
## 3 0.19662727 0.01853118 0.1602443189 0.23297986 29802.08 29807.35
## 4 0.17228950 0.01728078 0.1383614997 0.20618918 29819.91 29826.28
## 5 0.16542823 0.01724706 0.1315664326 0.19926177 29827.12 29833.49
## 6 0.15981253 0.01791238 0.1246444954 0.19495122 29837.70 29842.42
## 7 0.15914260 0.01803003 0.1237435839 0.19451208 29839.64 29844.04
## 8 0.15274821 0.01724275 0.1188948841 0.18657328 29839.89 29845.90
## 9 0.12863010 0.01835400 0.0925950005 0.16463512 29867.87 29872.31
## 10 0.04924124 0.01649194 0.0168619947 0.08159346 29905.31 29910.82
## 11 0.08260723 0.02551086 0.0325208111 0.13265185 29905.52 29910.98
## 12 0.04756556 0.01657257 0.0150280123 0.08007595 29906.18 29911.87
## 13 0.03594293 0.01657555 0.0033995390 0.06845916 29909.46 29914.46
## 14 0.20352965 0.10318081 0.0009509321 0.40593931 29910.81 29915.79
## Covariate Credible Type Var Dup
## 1 s_confor90 Inside Cov s_confor FALSE
## 2 s_hydro_m Inside Cov s_hydro_m FALSE
## 3 s_contall90 Inside Cov s_contall FALSE
## 4 s_rip150 Inside Cov s_rip FALSE
## 5 s_ripar120 Inside Cov s_ripar FALSE
## 6 s_hardtal150 Inside Cov s_hardtal FALSE
## 7 s_hardfor150 Inside Cov s_hardfor FALSE
## 8 s_riparian90 Inside Cov s_riparian FALSE
## 9 s_hardtall90 Inside Cov s_hardtall FALSE
## 10 s_herbopen90 Inside Cov s_herbopen FALSE
## 11 s_Paved_m Inside Cov s_Paved_m FALSE
## 12 s_herbop120 Inside Cov s_herbop FALSE
## 13 s_consh120 Inside Cov s_consh FALSE
## 14 s_Unpaved_m Inside Cov s_Unpaved_m FALSEwrite.csv(Fixed.effects.all3, "./Top_sig_effects.csv")4 Global Models
Check correlation between selected covariates
Best.Cov = Marten.pnts@data %>%
dplyr::select(s_confor90, s_hydro_m, s_contall90, s_rip150, s_ripar120, s_hardtal150, s_hardfor150, s_riparian90,
s_hardtall90, s_herbopen90, s_Paved_m, s_herbop120, s_consh120, s_Unpaved_m)
Cov.cor2 = cor(Best.Cov)
corrplot(Cov.cor2,
tl.cex=0.7,
number.cex = 0.6,
method = "number")
#The AC, AH, MC, and MH variables look highly correlated. Let's test with a different method:
library(perturb)##
## Attaching package: 'perturb'## The following object is masked from 'package:raster':
##
## reclassifyCI = colldiag(Cov.cor2)
CI #Last value in the Index column needs to be less than 30; it's still too high## Condition
## Index Variance Decomposition Proportions
## intercept s_confor90 s_hydro_m s_contall90 s_rip150 s_ripar120
## 1 1.000 0.000 0.000 0.000 0.000 0.000 0.000
## 2 1.218 0.000 0.000 0.000 0.000 0.000 0.000
## 3 1.391 0.000 0.000 0.000 0.000 0.000 0.000
## 4 1.712 0.000 0.000 0.002 0.001 0.000 0.000
## 5 2.103 0.000 0.000 0.002 0.001 0.000 0.000
## 6 2.424 0.000 0.000 0.002 0.000 0.000 0.000
## 7 2.918 0.000 0.000 0.008 0.000 0.000 0.000
## 8 3.000 0.000 0.001 0.001 0.006 0.000 0.000
## 9 24.842 0.007 0.000 0.009 0.000 0.000 0.000
## 10 36.204 0.009 0.000 0.010 0.000 0.018 0.000
## 11 58.185 0.520 0.076 0.523 0.001 0.000 0.000
## 12 65.199 0.460 0.089 0.436 0.000 0.000 0.000
## 13 203.334 0.000 0.019 0.001 0.013 0.000 0.000
## 14 397.780 0.003 0.813 0.007 0.977 0.982 0.999
## s_hardtal150 s_hardfor150 s_riparian90 s_hardtall90 s_herbopen90 s_Paved_m
## 1 0.000 0.000 0.000 0.000 0.000 0.000
## 2 0.000 0.000 0.000 0.000 0.000 0.000
## 3 0.000 0.000 0.000 0.001 0.000 0.000
## 4 0.000 0.000 0.000 0.000 0.000 0.000
## 5 0.000 0.000 0.000 0.000 0.000 0.006
## 6 0.000 0.000 0.000 0.000 0.000 0.024
## 7 0.000 0.000 0.000 0.000 0.000 0.001
## 8 0.000 0.000 0.000 0.000 0.000 0.000
## 9 0.004 0.007 0.000 0.465 0.000 0.011
## 10 0.000 0.000 0.032 0.001 0.000 0.017
## 11 0.000 0.000 0.005 0.162 0.210 0.489
## 12 0.000 0.000 0.003 0.154 0.788 0.445
## 13 0.996 0.993 0.000 0.210 0.001 0.004
## 14 0.000 0.000 0.960 0.007 0.001 0.002
## s_herbop120 s_consh120 s_Unpaved_m
## 1 0.000 0.000 0.000
## 2 0.000 0.002 0.000
## 3 0.000 0.000 0.000
## 4 0.000 0.000 0.001
## 5 0.000 0.002 0.005
## 6 0.000 0.001 0.003
## 7 0.000 0.011 0.020
## 8 0.000 0.052 0.001
## 9 0.000 0.002 0.012
## 10 0.001 0.006 0.023
## 11 0.527 0.302 0.500
## 12 0.470 0.327 0.432
## 13 0.001 0.031 0.003
## 14 0.000 0.264 0.000Remove highly correlated variables that are redundant and reassess.
# dropping s_confor90, s_ripar120, s_ripar150 and re-assessing
Best.Cov2 = Best.Cov %>% dplyr::select(-c("s_confor90","s_ripar120","s_rip150","s_hardfor150","s_hardtal150","s_herbop120")) #Have the high individual values
Cov.cor3 = cor(Best.Cov2)
CI = colldiag(Cov.cor3)
CI #Looks good (all values <3.0). Nested subsets of cover types highly correlated (as expected)## Condition
## Index Variance Decomposition Proportions
## intercept s_hydro_m s_contall90 s_riparian90 s_hardtall90 s_herbopen90
## 1 1.000 0.548 0.000 0.034 0.030 0.031 0.032
## 2 1.139 0.019 0.146 0.047 0.054 0.012 0.033
## 3 1.458 0.054 0.159 0.014 0.000 0.198 0.159
## 4 1.646 0.005 0.032 0.127 0.411 0.174 0.089
## 5 1.833 0.011 0.007 0.006 0.004 0.443 0.010
## 6 1.974 0.086 0.101 0.390 0.311 0.005 0.319
## 7 2.020 0.011 0.010 0.212 0.065 0.050 0.348
## 8 2.166 0.265 0.545 0.169 0.126 0.087 0.010
## s_Paved_m s_consh120 s_Unpaved_m
## 1 0.028 0.031 0.034
## 2 0.051 0.052 0.068
## 3 0.074 0.093 0.065
## 4 0.110 0.032 0.006
## 5 0.673 0.044 0.068
## 6 0.007 0.061 0.133
## 7 0.055 0.653 0.003
## 8 0.000 0.035 0.625Organize data
Marten.data = Marten.pnts@data
Global.lst = list(list(intercept1 = rep(1, dim(Marten.data)[1])),
list(s_confor90 = Marten.data[,"s_confor90"],
s_hydro_m = Marten.data[,"s_hydro_m"],
s_contall90 = Marten.data[,"s_contall90"],
s_rip150 = Marten.data[,"s_rip150"],
s_ripar120 = Marten.data[,"s_ripar120"],
s_hardtal150 = Marten.data[,"s_hardtal150"],
s_hardfor150 = Marten.data[,"s_hardfor150"],
s_riparian90 = Marten.data[,"s_riparian90"],
s_hardtall90 = Marten.data[,"s_hardtall90"],
s_herbopen90 = Marten.data[,"s_herbopen90"],
s_Paved_m = Marten.data[,"s_Paved_m"],
s_herbop120 = Marten.data[,"s_herbop120"],
s_consh120 = Marten.data[,"s_consh120"],
s_Unpaved_m = Marten.data[,"s_Unpaved_m"],
JDATE = Marten.data[,"jdate"],
NN = round(Marten.data[,"NNs"],3),
ID = Marten.data[,"id.x"],
YEAR = Marten.data[,"year"]))
Global.stk = inla.stack(data = list(RASTERVALU = Marten.data$t.RASTERVALU),
A = list(1,1),
effects = Global.lst,
tag = "G1.0") Run Full Model 1 All Covariates Included (some highly correlated)
nPrior = list(theta=list(prior = "normal", param=c(0, 5))) #flat prior
Formula = RASTERVALU ~ -1 + intercept1 + #response and intercept
f(ID, #grouping effect for individual marten
model="iid",
constr = TRUE,
hyper = nPrior) +
f(inla.group(NN, #distance as a smooth effect
n = 20, #
method = "cut"),
model="rw1",
constr = TRUE,
scale.model = TRUE,
hyper = nPrior ) +
f(JDATE,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(YEAR,
model="iid",
constr = TRUE,
hyper = nPrior) +
s_confor90 + s_hydro_m + s_contall90 + s_rip150 + s_ripar120 + s_hardtal150 + s_hardfor150 + s_riparian90 + s_hardtall90 + s_herbopen90 + s_Paved_m + s_herbop120 + s_consh120 + s_Unpaved_m
GModel.1 = inla(Formula,
data = inla.stack.data(Global.stk),
family = "gaussian",
verbose = TRUE,
control.fixed = list(prec = 0.001,
prec.intercept = 0.0001),
control.predictor = list(
A = inla.stack.A(Global.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.compute=list(dic = TRUE, cpo = TRUE, waic = TRUE))summary(GModel.1)##
## Call:
## c("inla(formula = Formula, family = \"gaussian\", data = inla.stack.data(Global.stk), ", " verbose = TRUE, control.compute = list(dic = TRUE, cpo = TRUE, ", " waic = TRUE), control.predictor = list(A = inla.stack.A(Global.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
## 2.1597 47.5138 0.9678 50.6413
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept1 0.7691 0.2798 0.2197 0.7691 1.3180 0.7691 0
## s_confor90 0.2064 0.1517 -0.0915 0.2064 0.5039 0.2064 0
## s_hydro_m -0.4861 0.0449 -0.5743 -0.4861 -0.3980 -0.4861 0
## s_contall90 0.0191 0.1366 -0.2490 0.0191 0.2870 0.0191 0
## s_rip150 0.1467 0.0811 -0.0125 0.1467 0.3058 0.1467 0
## s_ripar120 -0.0573 0.1265 -0.3057 -0.0573 0.1909 -0.0573 0
## s_hardtal150 0.1168 0.0862 -0.0524 0.1168 0.2858 0.1168 0
## s_hardfor150 0.0611 0.0806 -0.0972 0.0611 0.2193 0.0611 0
## s_riparian90 0.0810 0.0671 -0.0508 0.0810 0.2126 0.0810 0
## s_hardtall90 0.0110 0.0367 -0.0610 0.0110 0.0830 0.0110 0
## s_herbopen90 0.1032 0.0487 0.0075 0.1032 0.1988 0.1032 0
## s_Paved_m 0.0746 0.0250 0.0255 0.0746 0.1237 0.0746 0
## s_herbop120 -0.0272 0.0491 -0.1235 -0.0272 0.0691 -0.0272 0
## s_consh120 -0.0301 0.0522 -0.1326 -0.0301 0.0724 -0.0301 0
## s_Unpaved_m 0.2016 0.1030 -0.0006 0.2016 0.4036 0.2016 0
##
## Random effects:
## Name Model
## ID IID model
## inla.group(NN, n = 20, method = "cut") RW1 model
## JDATE IID model
## YEAR IID model
##
## Model hyperparameters:
## mean sd 0.025quant
## Precision for the Gaussian observations 0.5264 0.0082 0.5104
## Precision for ID 0.6628 0.1941 0.3512
## Precision for inla.group(NN, n = 20, method = "cut") 0.8738 0.3307 0.3959
## Precision for JDATE 0.4648 0.0592 0.3609
## Precision for YEAR 1.1849 0.5434 0.4422
## 0.5quant 0.975quant mode
## Precision for the Gaussian observations 0.5264 0.5426 0.5264
## Precision for ID 0.6405 1.1074 0.5976
## Precision for inla.group(NN, n = 20, method = "cut") 0.8179 1.6781 0.7166
## Precision for JDATE 0.4604 0.5930 0.4512
## Precision for YEAR 1.0822 2.5298 0.8957
##
## Expected number of effective parameters(std dev): 147.62(0.00)
## Number of equivalent replicates : 57.11
##
## Deviance Information Criterion (DIC) ...............: 29482.87
## Deviance Information Criterion (DIC, saturated) ....: 8581.99
## Effective number of parameters .....................: 147.62
##
## Watanabe-Akaike information criterion (WAIC) ...: 29493.52
## Effective number of parameters .................: 153.60
##
## Marginal log-Likelihood: -15077.70
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computedRun Global Model 2 Removing strongly correlated covariates. (s_confor90, s_rip150, s_ripar120, s_hardtal150, s_hardfor150, s_herbop120, s_consh120, s_riparian90)
Recheck for colinearity
Cov.cor3 = as.data.frame(Cov.cor2) %>% select(-c(s_confor90, s_rip150, s_ripar120, s_hardtal150, s_hardfor150, s_herbop120, s_consh120, s_riparian90))
CI = colldiag(Cov.cor3)
CI # Condition index <5.0 for remaining variables.## Condition
## Index Variance Decomposition Proportions
## intercept s_hydro_m s_contall90 s_hardtall90 s_herbopen90 s_Paved_m
## 1 1.000 0.029 0.000 0.031 0.039 0.023 0.028
## 2 1.175 0.002 0.123 0.064 0.003 0.021 0.044
## 3 1.475 0.000 0.047 0.002 0.119 0.208 0.089
## 4 1.779 0.000 0.000 0.265 0.058 0.290 0.022
## 5 1.814 0.000 0.004 0.011 0.303 0.007 0.488
## 6 2.150 0.015 0.393 0.317 0.036 0.004 0.000
## 7 4.857 0.954 0.433 0.310 0.443 0.447 0.328
## s_Unpaved_m
## 1 0.024
## 2 0.051
## 3 0.053
## 4 0.118
## 5 0.018
## 6 0.228
## 7 0.5075 Model at 90m resolution
Marten.data = Marten.pnts@data
Res90.lst = list(list(intercept1 = rep(1, dim(Marten.data)[1])),
list(s_hydro_m = Marten.data[,"s_hydro_m"],
s_contall = Marten.data[,"s_contall90"],
s_hardtall = Marten.data[,"s_hardtall90"],
s_herbopen = Marten.data[,"s_herbopen90"],
s_Paved_m = Marten.data[,"s_Paved_m"],
s_Unpaved_m = Marten.data[,"s_Unpaved_m"],
JDATE = Marten.data[,"jdate"],
NN = round(Marten.data[,"NNs"],3),
ID = Marten.data[,"id.x"],
YEAR = Marten.data[,"year"]))
Res90.stk = inla.stack(data = list(RASTERVALU = Marten.data$t.RASTERVALU),
A = list(1,1),
effects = Res90.lst,
tag = "res.90") nPrior = list(theta=list(prior = "normal", param=c(0, 5))) #flat prior
Formula = RASTERVALU ~ -1 + intercept1 + #response and intercept
f(ID, #grouping effect for individual marten
model="iid",
constr = TRUE,
hyper = nPrior) +
f(inla.group(NN, #distance as a smooth effect
n = 20, #
method = "cut"),
model="rw1",
constr = TRUE,
scale.model = TRUE,
hyper = nPrior ) +
f(JDATE,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(YEAR,
model="iid",
constr = TRUE,
hyper = nPrior) +
s_hydro_m + s_contall +
s_hardtall + s_herbopen + s_Paved_m + s_Unpaved_m
Model.90 = inla(Formula,
data = inla.stack.data(Res90.stk),
family = "gaussian",
verbose = TRUE,
control.fixed = list(prec = 0.001,
prec.intercept = 0.0001),
control.predictor = list(
A = inla.stack.A(Res90.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.compute=list(dic = TRUE, cpo = TRUE, waic = TRUE))summary(Model.90)##
## Call:
## c("inla(formula = Formula, family = \"gaussian\", data = inla.stack.data(Res90.stk), ", " verbose = TRUE, control.compute = list(dic = TRUE, cpo = TRUE, ", " waic = TRUE), control.predictor = list(A = inla.stack.A(Res90.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
## 2.2281 39.0685 1.1769 42.4735
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept1 1.2172 0.2780 0.6714 1.2172 1.7625 1.2172 0
## s_hydro_m -0.5169 0.0451 -0.6054 -0.5169 -0.4285 -0.5169 0
## s_contall 0.2107 0.0184 0.1746 0.2107 0.2468 0.2107 0
## s_hardtall 0.1584 0.0182 0.1227 0.1584 0.1940 0.1584 0
## s_herbopen 0.0759 0.0163 0.0440 0.0759 0.1079 0.0759 0
## s_Paved_m 0.0715 0.0252 0.0220 0.0715 0.1209 0.0715 0
## s_Unpaved_m 0.1136 0.1029 -0.0884 0.1136 0.3154 0.1136 0
##
## Random effects:
## Name Model
## ID IID model
## inla.group(NN, n = 20, method = "cut") RW1 model
## JDATE IID model
## YEAR IID model
##
## Model hyperparameters:
## mean sd 0.025quant
## Precision for the Gaussian observations 0.5189 0.0081 0.5035
## Precision for ID 0.6729 0.2005 0.3389
## Precision for inla.group(NN, n = 20, method = "cut") 0.7440 0.2802 0.3436
## Precision for JDATE 0.4539 0.0584 0.3526
## Precision for YEAR 1.1877 0.5427 0.4394
## 0.5quant 0.975quant mode
## Precision for the Gaussian observations 0.5187 0.5356 0.5179
## Precision for ID 0.6555 1.1161 0.6195
## Precision for inla.group(NN, n = 20, method = "cut") 0.6952 1.4267 0.6082
## Precision for JDATE 0.4491 0.5820 0.4384
## Precision for YEAR 1.0880 2.5203 0.9022
##
## Expected number of effective parameters(std dev): 139.87(0.00)
## Number of equivalent replicates : 60.27
##
## Deviance Information Criterion (DIC) ...............: 29591.44
## Deviance Information Criterion (DIC, saturated) ....: 8552.90
## Effective number of parameters .....................: 139.87
##
## Watanabe-Akaike information criterion (WAIC) ...: 29599.50
## Effective number of parameters .................: 143.84
##
## Marginal log-Likelihood: -15089.69
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computed6 Model at 120m resolution
Marten.data = Marten.pnts@data
Res120.lst = list(list(intercept1 = rep(1, dim(Marten.data)[1])),
list(s_hydro_m = Marten.data[,"s_hydro_m"],
s_contall = Marten.data[,"s_contall120"],
s_hardtall = Marten.data[,"s_hardtal120"],
s_herbopen = Marten.data[,"s_herbop120"],
s_Paved_m = Marten.data[,"s_Paved_m"],
s_Unpaved_m = Marten.data[,"s_Unpaved_m"],
JDATE = Marten.data[,"jdate"],
NN = round(Marten.data[,"NNs"],3),
ID = Marten.data[,"id.x"],
YEAR = Marten.data[,"year"]))
Res120.stk = inla.stack(data = list(RASTERVALU = Marten.data$t.RASTERVALU),
A = list(1,1),
effects = Res120.lst,
tag = "res.120") nPrior = list(theta=list(prior = "normal", param=c(0, 5))) #flat prior
Formula = RASTERVALU ~ -1 + intercept1 + #response and intercept
f(ID, #grouping effect for individual marten
model="iid",
constr = TRUE,
hyper = nPrior) +
f(inla.group(NN, #distance as a smooth effect
n = 20, #
method = "cut"),
model="rw1",
constr = TRUE,
scale.model = TRUE,
hyper = nPrior ) +
f(JDATE,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(YEAR,
model="iid",
constr = TRUE,
hyper = nPrior) +
s_hydro_m + s_contall +
s_hardtall + s_herbopen + s_Paved_m + s_Unpaved_m
Model.120 = inla(Formula,
data = inla.stack.data(Res120.stk),
family = "gaussian",
verbose = TRUE,
control.fixed = list(prec = 0.001,
prec.intercept = 0.0001),
control.predictor = list(
A = inla.stack.A(Res120.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.compute=list(dic = TRUE, cpo = TRUE, waic = TRUE))summary(Model.120)##
## Call:
## c("inla(formula = Formula, family = \"gaussian\", data = inla.stack.data(Res120.stk), ", " verbose = TRUE, control.compute = list(dic = TRUE, cpo = TRUE, ", " waic = TRUE), control.predictor = list(A = inla.stack.A(Res120.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.9339 38.7563 1.1125 41.8026
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept1 1.1292 0.2793 0.5808 1.1292 1.6771 1.1292 0
## s_hydro_m -0.5189 0.0451 -0.6075 -0.5189 -0.4305 -0.5189 0
## s_contall 0.1995 0.0183 0.1636 0.1995 0.2354 0.1995 0
## s_hardtall 0.1706 0.0179 0.1355 0.1706 0.2057 0.1706 0
## s_herbopen 0.0659 0.0163 0.0339 0.0659 0.0979 0.0659 0
## s_Paved_m 0.0695 0.0252 0.0200 0.0695 0.1189 0.0695 0
## s_Unpaved_m 0.1038 0.1032 -0.0988 0.1038 0.3063 0.1038 0
##
## Random effects:
## Name Model
## ID IID model
## inla.group(NN, n = 20, method = "cut") RW1 model
## JDATE IID model
## YEAR IID model
##
## Model hyperparameters:
## mean sd 0.025quant
## Precision for the Gaussian observations 0.5193 0.0082 0.5040
## Precision for ID 0.6587 0.1924 0.3505
## Precision for inla.group(NN, n = 20, method = "cut") 0.7740 0.2891 0.3307
## Precision for JDATE 0.4579 0.0592 0.3441
## Precision for YEAR 1.1658 0.5582 0.4485
## 0.5quant 0.975quant mode
## Precision for the Gaussian observations 0.5190 0.5363 0.5179
## Precision for ID 0.6364 1.0994 0.5935
## Precision for inla.group(NN, n = 20, method = "cut") 0.7355 1.4490 0.6553
## Precision for JDATE 0.4578 0.5755 0.4614
## Precision for YEAR 1.0447 2.5866 0.8483
##
## Expected number of effective parameters(std dev): 139.54(0.00)
## Number of equivalent replicates : 60.41
##
## Deviance Information Criterion (DIC) ...............: 29584.13
## Deviance Information Criterion (DIC, saturated) ....: 8544.55
## Effective number of parameters .....................: 139.54
##
## Watanabe-Akaike information criterion (WAIC) ...: 29592.39
## Effective number of parameters .................: 143.75
##
## Marginal log-Likelihood: -15085.97
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computed7 Model at 120m resolution
Marten.data = Marten.pnts@data
Res150.lst = list(list(intercept1 = rep(1, dim(Marten.data)[1])),
list(s_hydro_m = Marten.data[,"s_hydro_m"],
s_contall = Marten.data[,"s_contall150"],
s_hardtall = Marten.data[,"s_hardtal150"],
s_herbopen = Marten.data[,"s_herbop150"],
s_Paved_m = Marten.data[,"s_Paved_m"],
s_Unpaved_m = Marten.data[,"s_Unpaved_m"],
JDATE = Marten.data[,"jdate"],
NN = round(Marten.data[,"NNs"],3),
ID = Marten.data[,"id.x"],
YEAR = Marten.data[,"year"]))
Res150.stk = inla.stack(data = list(RASTERVALU = Marten.data$t.RASTERVALU),
A = list(1,1),
effects = Res150.lst,
tag = "res.150") nPrior = list(theta=list(prior = "normal", param=c(0, 5))) #flat prior
Formula = RASTERVALU ~ -1 + intercept1 + #response and intercept
f(ID, #grouping effect for individual marten
model="iid",
constr = TRUE,
hyper = nPrior) +
f(inla.group(NN, #distance as a smooth effect
n = 20, #
method = "cut"),
model="rw1",
constr = TRUE,
scale.model = TRUE,
hyper = nPrior ) +
f(JDATE,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(YEAR,
model="iid",
constr = TRUE,
hyper = nPrior) +
s_hydro_m + s_contall +
s_hardtall + s_herbopen + s_Paved_m + s_Unpaved_m
Model.150 = inla(Formula,
data = inla.stack.data(Res150.stk),
family = "gaussian",
verbose = TRUE,
control.fixed = list(prec = 0.001,
prec.intercept = 0.0001),
control.predictor = list(
A = inla.stack.A(Res150.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.compute=list(dic = TRUE, cpo = TRUE, waic = TRUE))summary(Model.120)##
## Call:
## c("inla(formula = Formula, family = \"gaussian\", data = inla.stack.data(Res120.stk), ", " verbose = TRUE, control.compute = list(dic = TRUE, cpo = TRUE, ", " waic = TRUE), control.predictor = list(A = inla.stack.A(Res120.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.9339 38.7563 1.1125 41.8026
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept1 1.1292 0.2793 0.5808 1.1292 1.6771 1.1292 0
## s_hydro_m -0.5189 0.0451 -0.6075 -0.5189 -0.4305 -0.5189 0
## s_contall 0.1995 0.0183 0.1636 0.1995 0.2354 0.1995 0
## s_hardtall 0.1706 0.0179 0.1355 0.1706 0.2057 0.1706 0
## s_herbopen 0.0659 0.0163 0.0339 0.0659 0.0979 0.0659 0
## s_Paved_m 0.0695 0.0252 0.0200 0.0695 0.1189 0.0695 0
## s_Unpaved_m 0.1038 0.1032 -0.0988 0.1038 0.3063 0.1038 0
##
## Random effects:
## Name Model
## ID IID model
## inla.group(NN, n = 20, method = "cut") RW1 model
## JDATE IID model
## YEAR IID model
##
## Model hyperparameters:
## mean sd 0.025quant
## Precision for the Gaussian observations 0.5193 0.0082 0.5040
## Precision for ID 0.6587 0.1924 0.3505
## Precision for inla.group(NN, n = 20, method = "cut") 0.7740 0.2891 0.3307
## Precision for JDATE 0.4579 0.0592 0.3441
## Precision for YEAR 1.1658 0.5582 0.4485
## 0.5quant 0.975quant mode
## Precision for the Gaussian observations 0.5190 0.5363 0.5179
## Precision for ID 0.6364 1.0994 0.5935
## Precision for inla.group(NN, n = 20, method = "cut") 0.7355 1.4490 0.6553
## Precision for JDATE 0.4578 0.5755 0.4614
## Precision for YEAR 1.0447 2.5866 0.8483
##
## Expected number of effective parameters(std dev): 139.54(0.00)
## Number of equivalent replicates : 60.41
##
## Deviance Information Criterion (DIC) ...............: 29584.13
## Deviance Information Criterion (DIC, saturated) ....: 8544.55
## Effective number of parameters .....................: 139.54
##
## Watanabe-Akaike information criterion (WAIC) ...: 29592.39
## Effective number of parameters .................: 143.75
##
## Marginal log-Likelihood: -15085.97
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computed8 Bayesian Model Averaging
myMods = list(Model.90, Model.120, Model.150)
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)
library(INLABMA)
#Average Fixed Effects
BMA.fx = INLABMA:::fitmatrixBMA(myMods, rep(1/length(myMods), length(myMods)), "summary.fixed")
BMA.fx## mean sd 0.025quant 0.5quant 0.975quant
## intercept1 1.14581602 0.28007428 0.59593578 1.14580814 1.69523737
## s_hydro_m -0.51789792 0.04507271 -0.60639085 -0.51789919 -0.42947884
## s_contall 0.19837919 0.01829473 0.16246047 0.19837868 0.23426794
## s_hardtall 0.16761234 0.01792427 0.13242094 0.16761183 0.20277436
## s_herbopen 0.06577665 0.01631447 0.03374585 0.06577619 0.09778072
## s_Paved_m 0.06983082 0.02518721 0.02037983 0.06983012 0.11924055
## s_Unpaved_m 0.10691475 0.10320815 -0.09571764 0.10691185 0.30937804
## mode kld
## intercept1 1.14581602 4.743971e-16
## s_hydro_m -0.51789792 6.586167e-30
## s_contall 0.19837919 3.409785e-15
## s_hardtall 0.16761234 3.557147e-15
## s_herbopen 0.06577665 8.178390e-16
## s_Paved_m 0.06983082 1.138960e-16
## s_Unpaved_m 0.10691475 2.731798e-17#Average Random effects
BMA.hyp = INLABMA:::fitmatrixBMA(myMods, probs, "summary.hyperpar")
BMA.hyp## mean sd
## Precision for the Gaussian observations 0.01390839 0.0002196675
## Precision for ID 0.01755927 0.0051752788
## Precision for inla.group(NN, n = 20, method = "cut") 0.02069374 0.0078026735
## Precision for JDATE 0.01222633 0.0015968253
## Precision for YEAR 0.03123040 0.0149480406
## 0.025quant 0.5quant
## Precision for the Gaussian observations 0.013496825 0.01390052
## Precision for ID 0.009381869 0.01691404
## Precision for inla.group(NN, n = 20, method = "cut") 0.008978089 0.01955820
## Precision for JDATE 0.009262773 0.01218068
## Precision for YEAR 0.012012068 0.02799172
## 0.975quant mode
## Precision for the Gaussian observations 0.01435842 0.01387526
## Precision for ID 0.02952591 0.01568847
## Precision for inla.group(NN, n = 20, method = "cut") 0.03914410 0.01730898
## Precision for JDATE 0.01550207 0.01215005
## Precision for YEAR 0.06927108 0.022732639 Comparison Models
Here we re-run without the spatial (NN) effect.
Run Non Spatial Model Dropping NN.
nPrior = list(theta=list(prior = "normal", param=c(0, 5)))
Formula = RASTERVALU ~ -1 + intercept1 + #response and intercept
f(ID, #grouping effect for individual marten
model="iid",
constr = TRUE,
hyper = nPrior) +
f(JDATE,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(YEAR,
model="iid",
constr = TRUE,
hyper = nPrior) +
s_hydro_m + s_contall90 +
s_hardtall90 + s_herbopen90 + s_Paved_m + s_Unpaved_m
Non.Spatial.mod = inla(Formula,
data = inla.stack.data(Global.stk),
family = "gaussian",
verbose = TRUE,
control.fixed = list(prec = 0.001,
prec.intercept = 0.0001),
control.predictor = list(
A = inla.stack.A(Global.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.compute=list(dic = TRUE, cpo = TRUE, waic = TRUE)) summary(Non.Spatial.mod)##
## Call:
## c("inla(formula = Formula, family = \"gaussian\", data = inla.stack.data(Global.stk), ", " verbose = TRUE, control.compute = list(dic = TRUE, cpo = TRUE, ", " waic = TRUE), control.predictor = list(A = inla.stack.A(Global.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.5380 39.3110 1.1289 41.9779
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept1 -1.4496 0.1401 -1.7247 -1.4496 -1.1747 -1.4496 0
## s_hydro_m -0.5151 0.0455 -0.6044 -0.5151 -0.4258 -0.5151 0
## s_contall90 0.2688 0.0171 0.2351 0.2688 0.3024 0.2688 0
## s_hardtall90 0.2123 0.0173 0.1783 0.2123 0.2463 0.2123 0
## s_herbopen90 0.1120 0.0160 0.0805 0.1120 0.1434 0.1120 0
## s_Paved_m 0.0707 0.0255 0.0207 0.0707 0.1206 0.0707 0
## s_Unpaved_m 0.1144 0.1044 -0.0905 0.1144 0.3191 0.1144 0
##
## Random effects:
## Name Model
## ID IID model
## JDATE IID model
## YEAR IID model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant
## Precision for the Gaussian observations 0.5060 0.0079 0.4907 0.5060
## Precision for ID 0.6657 0.1990 0.3582 0.6381
## Precision for JDATE 0.4424 0.0554 0.3420 0.4395
## Precision for YEAR 1.1779 0.5438 0.4445 1.0714
## 0.975quant mode
## Precision for the Gaussian observations 0.5216 0.5059
## Precision for ID 1.1340 0.5863
## Precision for JDATE 0.5598 0.4344
## Precision for YEAR 2.5373 0.8845
##
## Expected number of effective parameters(std dev): 133.98(0.00)
## Number of equivalent replicates : 62.92
##
## Deviance Information Criterion (DIC) ...............: 29801.66
## Deviance Information Criterion (DIC, saturated) ....: 8566.83
## Effective number of parameters .....................: 133.98
##
## Watanabe-Akaike information criterion (WAIC) ...: 29806.36
## Effective number of parameters .................: 135.72
##
## Marginal log-Likelihood: -15183.97
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computedRun spatial model without grouping by ID Day, or Year effects
nPrior = list(theta=list(prior = "normal", param=c(0, 5))) #flat prior
Formula = RASTERVALU ~ -1 + intercept1 +
f(inla.group(NN, #distance as a smooth effect
n = 20, #
method = "cut"),
model="rw1",
constr = TRUE,
scale.model = TRUE,
hyper = nPrior ) +
s_hydro_m + s_contall90 +
s_hardtall90 + s_herbopen90 + s_Paved_m + s_Unpaved_m
No.ID.mod = inla(Formula,
data = inla.stack.data(Global.stk),
family = "gaussian",
verbose = TRUE,
control.fixed = list(prec = 0.001,
prec.intercept = 0.0001),
control.predictor = list(
A = inla.stack.A(Global.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.compute=list(dic = TRUE, cpo = TRUE, waic = TRUE)) summary(No.ID.mod)##
## Call:
## c("inla(formula = Formula, family = \"gaussian\", data = inla.stack.data(Global.stk), ", " verbose = TRUE, control.compute = list(dic = TRUE, cpo = TRUE, ", " waic = TRUE), control.predictor = list(A = inla.stack.A(Global.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.5189 71.7807 0.7413 74.0410
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept1 0.8575 0.3851 0.1014 0.8575 1.6130 0.8575 0
## s_hydro_m 0.5103 0.0272 0.4569 0.5103 0.5637 0.5103 0
## s_contall90 0.6417 0.0290 0.5847 0.6417 0.6986 0.6417 0
## s_hardtall90 0.3686 0.0281 0.3134 0.3686 0.4238 0.3686 0
## s_herbopen90 0.1639 0.0269 0.1111 0.1639 0.2167 0.1639 0
## s_Paved_m 0.8124 0.0265 0.7604 0.8124 0.8644 0.8124 0
## s_Unpaved_m 0.9465 0.0270 0.8935 0.9465 0.9994 0.9465 0
##
## Random effects:
## Name Model
## inla.group(NN, n = 20, method = "cut") RW1 model
##
## Model hyperparameters:
## mean sd 0.025quant
## Precision for the Gaussian observations 0.1723 0.0015 0.1687
## Precision for inla.group(NN, n = 20, method = "cut") 0.7184 0.2671 0.3355
## 0.5quant 0.975quant mode
## Precision for the Gaussian observations 0.1726 0.1743 0.1741
## Precision for inla.group(NN, n = 20, method = "cut") 0.6722 1.3674 0.5896
##
## Expected number of effective parameters(std dev): 11.54(0.00)
## Number of equivalent replicates : 730.36
##
## Deviance Information Criterion (DIC) ...............: 39495.46
## Deviance Information Criterion (DIC, saturated) ....: 9275.18
## Effective number of parameters .....................: 11.54
##
## Watanabe-Akaike information criterion (WAIC) ...: 39494.99
## Effective number of parameters .................: 10.76
##
## Marginal log-Likelihood: -19814.94
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computedRun non-spatial model without grouping by ID or Year effects (just a linear model)
nPrior = list(theta=list(prior = "normal", param=c(0, 5)))
Formula = RASTERVALU ~ -1 + intercept1 +
s_hydro_m + s_contall90 +
s_hardtall90 + s_herbopen90 + s_Paved_m + s_Unpaved_m
No.Space.ID.mod = inla(Formula,
data = inla.stack.data(Global.stk),
family = "gaussian",
verbose = TRUE,
control.fixed = list(prec = 0.001,
prec.intercept = 0.0001),
control.predictor = list(
A = inla.stack.A(Global.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.compute=list(dic = TRUE, cpo = TRUE, waic = TRUE))summary(No.Space.ID.mod)##
## Call:
## c("inla(formula = Formula, family = \"gaussian\", data = inla.stack.data(Global.stk), ", " verbose = TRUE, control.compute = list(dic = TRUE, cpo = TRUE, ", " waic = TRUE), control.predictor = list(A = inla.stack.A(Global.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.1499 15.0471 1.1081 17.3051
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept1 -2.1054 0.0276 -2.1597 -2.1054 -2.0512 -2.1054 0
## s_hydro_m 0.5017 0.0287 0.4454 0.5017 0.5580 0.5017 0
## s_contall90 0.7207 0.0280 0.6658 0.7207 0.7756 0.7207 0
## s_hardtall90 0.4444 0.0278 0.3898 0.4444 0.4990 0.4444 0
## s_herbopen90 0.2082 0.0277 0.1538 0.2082 0.2625 0.2082 0
## s_Paved_m 0.8102 0.0279 0.7554 0.8102 0.8650 0.8102 0
## s_Unpaved_m 0.9383 0.0284 0.8825 0.9383 0.9941 0.9383 0
##
## The model has no random effects
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant
## Precision for the Gaussian observations 0.1566 0.0024 0.1519 0.1566
## 0.975quant mode
## Precision for the Gaussian observations 0.1614 0.1565
##
## Expected number of effective parameters(std dev): 7.008(0.00)
## Number of equivalent replicates : 1202.83
##
## Deviance Information Criterion (DIC) ...............: 39562.96
## Deviance Information Criterion (DIC, saturated) ....: 8437.82
## Effective number of parameters .....................: 7.008
##
## Watanabe-Akaike information criterion (WAIC) ...: 39561.82
## Effective number of parameters .................: 5.848
##
## Marginal log-Likelihood: -19838.74
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computed10 Results
Compare Model Metrics. The two global models look comparable. Global 1 has lowest DIC/WAIC, but contains mutliple correlated variables. Global 2, with uncorrelated variables and all the random effects seems to be the way to go. The Linear model without any spatial or grouping effects performs the worst.
Models = c("Global 1",
"Model.90",
"Model.120",
"Model.150",
"Non-Spatial",
"No ID", "Linear Model")
#Deviance information criteria
S1.DICs = c(GModel.1$dic$dic,
Model.90$dic$dic,
Model.120$dic$dic,
Model.150$dic$dic,
Non.Spatial.mod$dic$dic,
No.ID.mod$dic$dic, No.Space.ID.mod$dic$dic)
S1.WAICs = c(GModel.1$waic$waic,
Model.90$waic$waic,
Model.120$waic$waic,
Model.150$waic$waic,
Non.Spatial.mod$waic$waic,
No.ID.mod$waic$waic, No.Space.ID.mod$waic$waic)
Model_out = as.data.frame(cbind(Model = Models,
DIC = round(S1.DICs, 2),
WAIC = round(S1.WAICs, 2)))
kable(Model_out, caption = "Model Comparison") %>%
kable_styling("striped", full_width = F) %>%
row_spec(0, font_size = 20) | Model | DIC | WAIC |
|---|---|---|
| Global 1 | 29482.87 | 29493.52 |
| Model.90 | 29591.44 | 29599.5 |
| Model.120 | 29584.13 | 29592.39 |
| Model.150 | 29586.97 | 29595.59 |
| Non-Spatial | 29801.66 | 29806.36 |
| No ID | 39495.46 | 39494.99 |
| Linear Model | 39562.96 | 39561.82 |
Global Model 2 (Model.90) Fixed Effect Estimates
FE.table = Model.90$summary.fixed[,c(1:3,5)]
names(FE.table) = c("Mean", "sd", "Q0.025", "Q0.975")
kable(FE.table, caption = "Fixed Effects") %>%
kable_styling("striped", full_width = F) %>%
row_spec(0, font_size = 20) | Mean | sd | Q0.025 | Q0.975 | |
|---|---|---|---|---|
| intercept1 | 1.2171650 | 0.2779907 | 0.6713754 | 1.7624990 |
| s_hydro_m | -0.5168758 | 0.0450746 | -0.6053725 | -0.4284529 |
| s_contall | 0.2107209 | 0.0184165 | 0.1745632 | 0.2468484 |
| s_hardtall | 0.1583540 | 0.0181610 | 0.1226979 | 0.1939804 |
| s_herbopen | 0.0759395 | 0.0162831 | 0.0439703 | 0.1078819 |
| s_Paved_m | 0.0714674 | 0.0251916 | 0.0220078 | 0.1208856 |
| s_Unpaved_m | 0.1135738 | 0.1028763 | -0.0884070 | 0.3153860 |
Plot Fixed Effects (Model.90)
myNames = names(Model.90$marginals.fix)
myNames## [1] "intercept1" "s_hydro_m" "s_contall" "s_hardtall" "s_herbopen"
## [6] "s_Paved_m" "s_Unpaved_m"myLevels = length(names(Model.90$marginals.fix))
myLevels## [1] 7for(i in 1:myLevels){
tmp.df = as.data.frame(Model.90$marginals.fix[[i]])
tmp.df$Name = myNames[i]
if(i == 1){fixed.df = tmp.df
} else{fixed.df = rbind(fixed.df, tmp.df)}
}
CI.df = as.data.frame(Model.90$summary.fixed)[,c(3,5)]
names(CI.df) = c("Low","High")
CI.df$Name = row.names(CI.df)
fixed.df = fixed.df %>% filter(Name != "intercept1")
CI.df = CI.df %>% filter(Name != "intercept1")
fixed.df$Name = as.factor(fixed.df$Name)
fixed.df$Name = mapvalues(fixed.df$Name,
from = c("s_contall", "s_hardtall", "s_herbopen", "s_hydro_m", "s_Paved_m", "s_Unpaved_m"),
to = c("Conifer", "Hardwood", "Open Herbaceous", "Hydrology", "Paved", "Unpaved"))
CI.df$Name = as.factor(CI.df$Name)
CI.df$Name = mapvalues(CI.df$Name,
from = c("s_contall", "s_hardtall", "s_herbopen", "s_hydro_m", "s_Paved_m", "s_Unpaved_m"),
to = c("Conifer", "Hardwood", "Open Herbaceous", "Hydrology", "Paved", "Unpaved"))
ggplot(fixed.df, aes(x, y)) +
facet_wrap(~Name, ncol = 2, scales = "free") +
geom_line(size = 1, col = "black") +
geom_vline(data = CI.df,
aes(xintercept = Low), col = "gray") +
geom_vline(data = CI.df,
aes(xintercept = High), col = "gray") +
geom_vline(xintercept = 0, col = "red") +
xlab(NULL) +
ylab("Density") +
ggtitle("Fixed Effects") +
theme_classic() +
theme(legend.title = element_text(size=16, face="bold"),
legend.key = element_rect(fill = "white", linetype=0),
legend.background = element_rect(fill = "white", linetype=0),
legend.position = "bottom",
legend.text = element_text(size=12),
plot.title = element_text(hjust = 0.5),
axis.title.y = element_text(face="bold", size=22),
axis.title.x = element_text(face="bold", size=18),
title = element_text(face="bold", size=24, hjust=0.5),
strip.text.x = element_text(face="bold", size = 14, colour = "black"),
strip.background = element_blank(),
panel.border = element_blank(),
axis.text.y = element_text(face="bold", size=18),
axis.text.x = element_text(face="bold", size=18, vjust=0.5,
hjust=0.5, angle=0))
Plot Fixed Effects (All Resolutions)
myNames = names(Model.90$marginals.fix)
myNames## [1] "intercept1" "s_hydro_m" "s_contall" "s_hardtall" "s_herbopen"
## [6] "s_Paved_m" "s_Unpaved_m"myLevels = length(names(Model.90$marginals.fix))
myLevels## [1] 7for(i in 1:myLevels){
tmp.df = as.data.frame(Model.90$marginals.fix[[i]])
tmp.df$Name = myNames[i]
tmp.df$Model = "90 m"
tmp.df2 = as.data.frame(Model.120$marginals.fix[[i]])
tmp.df2$Name = myNames[i]
tmp.df2$Model = "120 m"
tmp.df3 = as.data.frame(Model.150$marginals.fix[[i]])
tmp.df3$Name = myNames[i]
tmp.df3$Model = "150m "
tmp.df = rbind(tmp.df, tmp.df2, tmp.df3)
if(i == 1){fixed2.df = tmp.df
} else{fixed2.df = rbind(fixed2.df, tmp.df)}
}
CI.df = as.data.frame(Model.90$summary.fixed)[,c(3,5)]
names(CI.df) = c("Low","High")
CI.df$Name = row.names(CI.df)
fixed2.df = fixed2.df %>% filter(Name != "intercept1")
CI.df = CI.df %>% filter(Name != "intercept1")
fixed2.df$Name = as.factor(fixed2.df$Name)
fixed2.df$Name = mapvalues(fixed2.df$Name,
from = c("s_contall", "s_hardtall", "s_herbopen", "s_hydro_m", "s_Paved_m", "s_Unpaved_m"),
to = c("Conifer", "Hardwood", "Open Herbaceous", "Hydrology", "Paved", "Unpaved"))
CI.df$Name = as.factor(CI.df$Name)
CI.df$Name = mapvalues(CI.df$Name,
from = c("s_contall", "s_hardtall", "s_herbopen", "s_hydro_m", "s_Paved_m", "s_Unpaved_m"),
to = c("Conifer", "Hardwood", "Open Herbaceous", "Hydrology", "Paved", "Unpaved"))
ggplot(fixed2.df, aes(x, y, group=Model,
linetype = Model,
col = Model)) +
facet_wrap(~Name, ncol = 2, scales = "free") +
geom_line(size = 0.75, col = "black") +
geom_vline(data = CI.df,
aes(xintercept = Low), col = "gray") +
geom_vline(data = CI.df,
aes(xintercept = High), col = "gray") +
geom_vline(xintercept = 0, col = "red") +
scale_linetype_manual(values=c("solid", "dashed", "dotted")) +
scale_colour_manual(values=c("red", "gray", "black")) +
xlab(NULL) +
ylab("Density") +
ggtitle(" ") +
theme_classic() +
theme(legend.title = element_text(size=16, face="bold"),
legend.key = element_rect(fill = "white", linetype=0),
legend.background = element_rect(fill = "white", linetype=0),
legend.position = "bottom",
legend.text = element_text(size=12),
plot.title = element_text(hjust = 0.5),
axis.title.y = element_text(face="bold", size=22),
axis.title.x = element_text(face="bold", size=18),
title = element_text(face="bold", size=24, hjust=0.5),
strip.text.x = element_text(face="bold", size = 14, colour = "black"),
strip.background = element_blank(),
panel.border = element_blank(),
axis.text.y = element_text(face="bold", size=18),
axis.text.x = element_text(face="bold", size=18, vjust=0.5,
hjust=0.5, angle=0))
ID Effect by Marten ID
mic.d = as.data.frame(Model.90$summary.random$ID)
names(mic.d) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975", "mode", "kld")
ggplot(mic.d, aes(factor(ID), y=Mean)) +
geom_point(size=2, pch=19, col = "red") +
geom_linerange(aes(ymin=Q025, ymax=Q975), colour="black") +
geom_hline(yintercept = 0,
linetype = "solid",
col = "darkgray",
size = 0.5) +
theme_classic() +
xlab("Marten ID") +
ylab("RASTERVALU (log)") +
ggtitle("Marten ID Effect") +
theme(plot.title = element_text(hjust = 0.5),
axis.title.y = element_text(face="bold", size=18),
axis.title.x = element_text(face="bold", size=18),
title = element_text(face="bold", size=18, hjust=0.5),
strip.text.x = element_text(face="bold", size = 14, colour = "black"),
axis.text.y = element_text(face="bold", size=14),
axis.text.x = element_text(face="bold", size=14, vjust=0.5,
hjust=0.5, angle=90))
Effect by Year
mic.d = as.data.frame(Model.90$summary.random$YEAR)
names(mic.d) = c("YEAR", "Mean", "sd", "Q025", "Q50", "Q975", "mode", "kld")
ggplot(mic.d, aes(factor(YEAR), y=Mean)) +
geom_point(size=2, pch=19, col = "red") +
geom_linerange(aes(ymin=Q025, ymax=Q975), colour="black") +
geom_hline(yintercept = 0,
linetype = "solid",
col = "darkgray",
size = 0.5) +
theme_classic() +
xlab("Year") +
ylab("RASTERVALU (log)") +
ggtitle("Year Effect") +
theme(plot.title = element_text(hjust = 0.5),
axis.title.y = element_text(face="bold", size=18),
axis.title.x = element_text(face="bold", size=18),
title = element_text(face="bold", size=18, hjust=0.5),
strip.text.x = element_text(face="bold", size = 14, colour = "black"),
axis.text.y = element_text(face="bold", size=14),
axis.text.x = element_text(face="bold", size=14, vjust=0.5,
hjust=0.5, angle=90))
JDate Effect
mic.d = as.data.frame(Model.90$summary.random$JDATE)
names(mic.d) = c("Day", "Mean", "sd", "Q025", "Q50", "Q975", "mode", "kld")
ggplot(mic.d, aes(Day, y=Mean)) +
geom_point(size=2, pch=19, col = "red") +
geom_linerange(aes(ymin=Q025, ymax=Q975), colour="black") +
geom_hline(yintercept = 0,
linetype = "solid",
col = "darkgray",
size = 0.5) +
scale_x_continuous(breaks = seq(64, 325, 5), labels=seq(64, 325, 5)) +
theme_classic() +
xlab("JDATE") +
ylab("RASTERVALU (log)") +
ggtitle("JDATE Effect") +
theme(plot.title = element_text(hjust = 0.5),
axis.title.y = element_text(face="bold", size=18),
axis.title.x = element_text(face="bold", size=18),
title = element_text(face="bold", size=18, hjust=0.5),
strip.text.x = element_text(face="bold", size = 14, colour = "black"),
axis.text.y = element_text(face="bold", size=14),
axis.text.x = element_text(face="bold", size=14, vjust=0.5,
hjust=0.5, angle=90))
Spatial Effect (Neighbor distance)
mic.d = as.data.frame(Model.90$summary.random$`inla.group(NN, n = 20, method = "cut")`)
names(mic.d) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975", "mode", "kld")
mic.d$ID = mic.d$ID*100 #Back to original scale, meters
ggplot(mic.d, aes(ID, Mean), lab) +
geom_smooth(col = "black",
linetype= "solid",
method = "loess",
span = 0.9,
se = FALSE,
lwd = 1) +
geom_smooth(data = mic.d, aes(ID, Q025),
col = "grey",
method = "loess",
span = 0.9,
se = FALSE,
linetype= "dashed") +
geom_smooth(data = mic.d, aes(ID, Q975),
col = "grey",
method = "loess",
span = 0.9,
se = FALSE,
linetype= "dashed") +
geom_hline(yintercept = 0,
linetype = "solid",
col = "darkgray",
size = 0.5) +
ylab("RASTERVALU (log)") +
xlab("Distance (meters)") +
ggtitle("Spatial Effect") +
theme_classic() +
theme(plot.title = element_text(hjust = 0.5),
title = element_text(face="bold", size=18, hjust=0.5),
axis.text=element_text(size=16),
strip.text = element_text(face="bold", size = 20),
axis.title.y = element_text(face="bold", size = 20),
axis.text.x = element_text(face="bold", size=14, vjust=0.5, angle=90),
axis.title.x = element_text(face="bold", size = 20))
11 Validation
Select records from each marten.
#Testing
set.seed(1976)
test.data = as.data.frame(
Marten.pnts@data %>%
group_by(id.x) %>%
sample_frac(0.20, replace=FALSE))
dim(test.data) ## [1] 1686 76#Training set 80%
train.data = anti_join(Marten.pnts@data, test.data) ## Joining, by = c("Field1", "FID", "Duration.x", "DOP.x", "Satellites.x", "id.x", "Longitude.x", "Latitude.x", "RASTERVALU", "hardfor90", "confor90", "hardtall90", "contall90", "hardsh90", "consh90", "herbopen90", "riparian90", "hardfor120", "confor120", "hardtal120", "contall120", "hardsh120", "consh120", "herbop120", "ripar120", "hardfor150", "confor150", "hardtal150", "contall150", "hardsh150", "consh150", "herbop150", "rip150", "Paved_m", "Unpaved_m", "hydro_m", "GMT.Time", "Duration.y", "DOP.y", "Satellites.y", "id.y", "Longitude.y", "Latitude.y", "moda", "jdate", "year", "s_hardfor90", "s_confor90", "s_hardtall90", "s_contall90", "s_hardsh90", "s_consh90", "s_herbopen90", "s_riparian90", "s_hardfor120", "s_confor120", "s_hardtal120", "s_contall120", "s_hardsh120", "s_consh120", "s_herbop120", "s_ripar120", "s_hardfor150", "s_confor150", "s_hardtal150", "s_contall150", "s_hardsh150", "s_consh150", "s_herbop150", "s_rip150", "s_Paved_m", "s_Unpaved_m", "s_hydro_m", "NN", "NNs", "t.RASTERVALU")dim(train.data)## [1] 6744 76Organize training data to fit the predictive model. Using 80 percent to create a data set called Train.stk.
Train.lst = list(list(intercept1 = rep(1, dim(train.data)[1])),
list(s_hydro_m = train.data[,"s_hydro_m"],
s_contall90 = train.data[,"s_contall90"],
s_hardtall90 = train.data[,"s_hardtall90"],
s_herbopen90 = train.data[,"s_herbopen90"],
s_Paved_m = train.data[,"s_Paved_m"],
s_Unpaved_m = train.data[,"s_Unpaved_m"],
JDATE = train.data[,"jdate"],
NN = round(train.data[,"NNs"],3),
ID = train.data[,"id.x"],
YEAR = train.data[,"year"]))
Train.stk = inla.stack(data = list(RASTERVALU = train.data$t.RASTERVALU),
A = list(1,1),
effects = Train.lst,
tag = "train.0") Organize testing data to be predicted by trained model model. Using 20 percent to create a data set called Test.stk.
Test.lst = list(list(intercept1 = rep(1, dim(test.data)[1])),
list(s_confor90 = test.data[,"s_confor90"],
s_hydro_m = test.data[,"s_hydro_m"],
s_contall90 = test.data[,"s_contall90"],
s_hardtall90 = test.data[,"s_hardtall90"],
s_herbopen90 = test.data[,"s_herbopen90"],
s_Paved_m = test.data[,"s_Paved_m"],
s_Unpaved_m = test.data[,"s_Unpaved_m"],
JDATE = test.data[,"jdate"],
NN = round(test.data[,"NNs"],3),
ID = test.data[,"id.x"],
YEAR = test.data[,"year"]))
test.data$NA_RASTERVALU = NA #Remove the RASTERVALU in the testing set so it can't be used; these NA values will be predicted.
Test.stk = inla.stack(data = list(RASTERVALU = test.data$NA_RASTERVALU), #NA values to be predicted
A = list(1,1),
effects = Test.lst,
tag = "test.0") #This label will be used to pull results laterJoin training and testing data to feed into model.
Validation.stk = inla.stack(Train.stk, Test.stk)Run formula from best Global model from above.
nPrior = list(theta=list(prior = "normal", param=c(0, 5))) #flat prior
Formula = RASTERVALU ~ -1 + intercept1 + #response and intercept
f(ID, #grouping effect for individual marten
model="iid",
constr = TRUE,
hyper = nPrior) +
f(inla.group(NN, #distance as a smooth effect
n = 20, #
method = "cut"),
model="rw1",
constr = TRUE,
scale.model = TRUE,
hyper = nPrior ) +
f(JDATE,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(YEAR,
model="iid",
constr = TRUE,
hyper = nPrior) +
s_hydro_m + s_contall90 +
s_hardtall90 + s_herbopen90 + s_Paved_m + s_Unpaved_m
Validation.mod = inla(Formula,
data = inla.stack.data(Validation.stk), #Validation data: 80% to train, 20% to be predicted.
family = "gaussian",
verbose = TRUE,
control.fixed = list(prec = 0.001,
prec.intercept = 0.0001),
control.predictor = list(
A = inla.stack.A(Validation.stk), #Need to provide the data twice
compute = TRUE,
link = 1),
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)) Put predicted values in new column to compare to RASTERVALU
idat = inla.stack.index(Validation.stk, "test.0")$data #Index to identify testing data locations from full data set
test.data$Prediction = plogis(Validation.mod$summary.fitted.values$mean[idat]) #Backtransform fitted values for tetsing data only
test.data$Low = plogis(Validation.mod$summary.fitted.values$`0.025quant`[idat])
test.data$High = plogis(Validation.mod$summary.fitted.values$`0.975quant`[idat])Compare predicted to actual RASTERVALU
#Correlation test
cor.test(test.data$Prediction, test.data$RASTERVALU)##
## Pearson's product-moment correlation
##
## data: test.data$Prediction and test.data$RASTERVALU
## t = 40.474, df = 1684, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.6771698 0.7256226
## sample estimates:
## cor
## 0.7022083#Calculate error
test.data$Error = test.data$RASTERVALU - test.data$Prediction
#Root Mean Squared Error
rmse = function(error)
{
sqrt(mean(error^2))
}
rmse(test.data$Error)## [1] 0.2089019#Mean Absolute Error
mae = function(error)
{
mean(abs(error))
}
mae(test.data$Error)## [1] 0.1381893#Plot
range(test.data$RASTERVALU)## [1] 0.000324589 0.988708973range(test.data$Prediction) ## [1] 0.0007083507 0.9858150639#length(which(test.data$Prediction > 1)) # No longer issue
ggplot(test.data, aes(RASTERVALU, Prediction)) +
geom_point(shape = 19,
col="lightgray",
size = 1) +
geom_smooth(col = "black",
linetype= "solid",
method = "lm",
se = FALSE,
lwd = 1) +
xlim(0,1) +
ylim(0,1) +
ylab("Actual") +
xlab("Predicted") +
ggtitle(" ") +
theme_classic() +
theme(plot.title = element_text(hjust = 0.5),
title = element_text(face="bold", size=18, hjust=0.5),
axis.text=element_text(size=16),
strip.text = element_text(face="bold", size = 20),
axis.title.y = element_text(face="bold", size = 20),
axis.text.x = element_text(face="bold", size=14, vjust=0.5, angle=90),
axis.title.x = element_text(face="bold", size = 20))