Marten Scale
1 Prepare Data
Set WD
#setwd("C:/Users/roloff/Desktop/Projects/SaultTribeMarten")
setwd("E:/Marten_2")
#setwd("C:/Users/roloff/Desktop/Temp/RRunFeb6")Load Libraries
Needed.packages = c("rgdal", "raster", "corrplot", "spdep",
"rgeos","maptools","mapproj","GISTools",
"sp","ggplot2","ggthemes","plyr","dplyr",
"kableExtra")
for(p in Needed.packages){
if(!require(p, character.only = TRUE)) install.packages(p)
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#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
## 1 0 0.01923077 0.7692308 0.01923077 0.4615385 0
## 2 0 0.01923077 0.6346154 0.01923077 0.6346154 0
## 3 0 0.07692308 0.8269231 0.07692308 0.7692308 0
## 4 0 0.21153846 0.3461538 0.21153846 0.3461538 0
## 5 0 0.13461539 0.8269231 0.13461539 0.8269231 0
## 6 0 0.00000000 0.8461538 0.00000000 0.8461538 0
## consh120 herbop120 ripar120 hardfor150 confor150 hardtal150 contall150
## 1 0.30769231 0 0 0.0125 0.8250 0.0125 0.5375
## 2 0.00000000 0 0 0.0625 0.4750 0.0625 0.4625
## 3 0.05769231 0 0 0.1375 0.8000 0.1375 0.7625
## 4 0.00000000 0 0 0.1000 0.3750 0.1000 0.3750
## 5 0.00000000 0 0 0.1000 0.9000 0.1000 0.9000
## 6 0.00000000 0 0 0.0000 0.8375 0.0000 0.8250
## hardsh150 consh150 herbop150 rip150 Paved_m Unpaved_m hydro_m
## 1 0 0.2875 0.0000 0 1488 943 3377
## 2 0 0.0125 0.0000 0 1480 1061 3479
## 3 0 0.0375 0.0000 0 1184 1015 3767
## 4 0 0.0000 0.0000 0 1145 1120 3820
## 5 0 0.0000 0.0000 0 436 1294 3970
## 6 0 0.0125 0.0125 0 70 1098 3864length(unique(martenscl$id))## [1] 13unique(martenscl$id)## [1] 7004m 6977m 5516m 1114m 5509m 5529m 5517m 8217m
## [9] houdini 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
## 1 0.00000 0.26923077 0.07692308 0.26923077 0.07692308 0
## 2 0.12500 0.03846154 0.03846154 0.03846154 0.03846154 0
## 3 0.00000 0.00000000 0.00000000 0.00000000 0.00000000 0
## 4 0.00000 0.00000000 0.00000000 0.00000000 0.00000000 0
## 5 0.00000 0.01923077 0.00000000 0.01923077 0.00000000 0
## 6 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 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
## 1 0 0.00000 0.26923077 0.07692308 0.26923077 0.07692308
## 2 0 0.12500 0.03846154 0.03846154 0.03846154 0.03846154
## 3 0 0.00000 0.00000000 0.00000000 0.00000000 0.00000000
## 4 0 0.00000 0.00000000 0.00000000 0.00000000 0.00000000
## 5 0 0.00000 0.01923077 0.00000000 0.01923077 0.00000000
## 6 0 0.03125 0.61538461 0.32692308 0.61538461 0.32692308
## hardsh120 consh120 herbop120 ripar120 hardfor150 confor150 hardtal150
## 1 0 0 0 0.03846154 0.2250 0.0750 0.2250
## 2 0 0 0 0.03846154 0.0125 0.0125 0.0125
## 3 0 0 0 0.00000000 0.0000 0.0000 0.0000
## 4 0 0 0 0.00000000 0.0000 0.0000 0.0000
## 5 0 0 0 0.00000000 0.0000 0.0000 0.0000
## 6 0 0 0 0.05769231 0.5000 0.3500 0.5000
## contall150 hardsh150 consh150 herbop150 rip150 Paved_m Unpaved_m hydro_m
## 1 0.0750 0 0 0 0.025 10 185 553
## 2 0.0125 0 0 0 0.025 70 130 485
## 3 0.0000 0 0 0 0.000 136 59 428
## 4 0.0000 0 0 0 0.000 136 60 427
## 5 0.0000 0 0 0 0.000 134 59 437
## 6 0.3500 0 0 0 0.125 145 339 870
## GMT.Time Duration.y DOP.y Satellites.y id.y Longitude.y
## 1 03/05/2016 19:30 59 2.4 4 2254m 675505.0
## 2 03/05/2016 19:45 6 1.8 4 2254m 675562.8
## 3 03/05/2016 20:00 3 1.8 4 2254m 675631.1
## 4 03/05/2016 20:15 2 2.2 4 2254m 675630.5
## 5 03/05/2016 20:30 2 2.8 4 2254m 675630.1
## 6 03/05/2016 20:45 9 3.4 4 2254m 675362.8
## Latitude.y moda jdate year s_hardfor90 s_confor90 s_hardtall90
## 1 5136966 2020-03-05 65 2016 1.45620350 0.06917282 1.52546893
## 2 5137027 2020-03-05 65 2016 0.03757526 -0.09513455 0.06185808
## 3 5136975 2020-03-05 65 2016 -0.36774710 -0.42374930 -0.35631645
## 4 5136985 2020-03-05 65 2016 -0.36774710 -0.42374930 -0.35631645
## 5 5136948 2020-03-05 65 2016 -0.16508592 -0.42374930 -0.14722918
## 6 5136545 2020-03-05 65 2016 3.88813763 0.89070969 4.03451610
## s_contall90 s_hardsh90 s_consh90 s_herbopen90 s_riparian90 s_hardfor120
## 1 0.14368603 -0.1089675 -0.1766424 -0.1242558 -0.26375740 1.61567196
## 2 -0.03794139 -0.1089675 -0.1766424 -0.1242558 0.87941736 -0.05509775
## 3 -0.40119623 -0.1089675 -0.1766424 -0.1242558 -0.26375740 -0.33355936
## 4 -0.40119623 -0.1089675 -0.1766424 -0.1242558 -0.26375740 -0.33355936
## 5 -0.40119623 -0.1089675 -0.1766424 -0.1242558 -0.26375740 -0.19432855
## 6 1.05182313 -0.1089675 -0.1766424 -0.1242558 0.02203629 4.12182651
## s_confor120 s_hardtal120 s_contall120 s_hardsh120 s_consh120 s_herbop120
## 1 0.05854051 1.68966047 0.1294811 -0.1063301 -0.1682362 -0.1259219
## 2 -0.16030494 -0.03512312 -0.1158332 -0.1063301 -0.1682362 -0.1259219
## 3 -0.37915037 -0.32258704 -0.3611475 -0.1063301 -0.1682362 -0.1259219
## 4 -0.37915037 -0.32258704 -0.3611475 -0.1063301 -0.1682362 -0.1259219
## 5 -0.37915037 -0.17885508 -0.3611475 -0.1063301 -0.1682362 -0.1259219
## 6 1.48103586 4.27683585 1.7240241 -0.1063301 -0.1682362 -0.1259219
## s_ripar120 s_hardfor150 s_confor150 s_hardtal150 s_contall150
## 1 0.1410059 1.6021914 0.1460311 1.6761332 0.2260010
## 2 0.1410059 -0.1980222 -0.2623019 -0.1842321 -0.2369178
## 3 -0.2423773 -0.3039172 -0.3439685 -0.2936654 -0.3295015
## 4 -0.2423773 -0.3039172 -0.3439685 -0.2936654 -0.3295015
## 5 -0.2423773 -0.3039172 -0.3439685 -0.2936654 -0.3295015
## 6 0.3326975 3.9318797 1.9426962 4.0836649 2.2628437
## s_hardsh150 s_consh150 s_herbop150 s_rip150 s_Paved_m s_Unpaved_m
## 1 -0.1007484 -0.1587382 -0.1244583 0.05927052 -1.1602269 -1.190281
## 2 -0.1007484 -0.1587382 -0.1244583 0.05927052 -1.0151395 -1.202578
## 3 -0.1007484 -0.1587382 -0.1244583 -0.22632565 -0.8555433 -1.218453
## 4 -0.1007484 -0.1587382 -0.1244583 -0.22632565 -0.8555433 -1.218229
## 5 -0.1007484 -0.1587382 -0.1244583 -0.22632565 -0.8603795 -1.218453
## 6 -0.1007484 -0.1587382 -0.1244583 1.20165522 -0.8337802 -1.155849
## s_hydro_m
## 1 -0.31487316
## 2 -0.39481041
## 3 -0.46181662
## 4 -0.46299217
## 5 -0.45123669
## 6 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.61-0 (nickname: 'Puppy zoomies')
## For an introduction to spatstat, type 'beginner'##
## Note: R version 3.5.3 (2019-03-11) is more than 9 months old; we strongly recommend upgrading to the latest version##
## 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
Looks like s_Area_AC150, s_Area_AH120, s_Area_AM150, s_Area_HO150, s_Area_MC150, s_Area_MH120, s_Dist_Paved, s_Dist_River, s_Dist_Unpav are the winners with good agreement between the DIC and WAIC.
#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
## 1 0.000 0.000 0.000 0.000 0.000
## 2 0.000 0.000 0.000 0.000 0.000
## 3 0.000 0.000 0.000 0.001 0.000
## 4 0.000 0.000 0.000 0.000 0.000
## 5 0.000 0.000 0.000 0.000 0.000
## 6 0.000 0.000 0.000 0.000 0.000
## 7 0.000 0.000 0.000 0.000 0.000
## 8 0.000 0.000 0.000 0.000 0.000
## 9 0.004 0.007 0.000 0.465 0.000
## 10 0.000 0.000 0.032 0.001 0.000
## 11 0.000 0.000 0.005 0.162 0.210
## 12 0.000 0.000 0.003 0.154 0.788
## 13 0.996 0.993 0.000 0.210 0.001
## 14 0.000 0.000 0.960 0.007 0.001
## s_Paved_m s_herbop120 s_consh120 s_Unpaved_m
## 1 0.000 0.000 0.000 0.000
## 2 0.000 0.000 0.002 0.000
## 3 0.000 0.000 0.000 0.000
## 4 0.000 0.000 0.000 0.001
## 5 0.006 0.000 0.002 0.005
## 6 0.024 0.000 0.001 0.003
## 7 0.001 0.000 0.011 0.020
## 8 0.000 0.000 0.052 0.001
## 9 0.011 0.000 0.002 0.012
## 10 0.017 0.001 0.006 0.023
## 11 0.489 0.527 0.302 0.500
## 12 0.445 0.470 0.327 0.432
## 13 0.004 0.001 0.031 0.003
## 14 0.002 0.000 0.264 0.000#dropping the MH150 and re-assessing
#Best.Cov2 = Best.Cov %>% dplyr::select(-c("s_rip150", "s_Area_MC120")) #Have the high individual values
#Cov.cor3 = cor(Best.Cov2)
#CI = colldiag(Cov.cor3)
#CI #Now less than 30Organize 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 are strongly 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.3068 50.1461 1.1360 53.5889
##
## 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
## Precision for the Gaussian observations 0.5264 0.0082
## Precision for ID 0.6628 0.1941
## Precision for inla.group(NN, n = 20, method = "cut") 0.8738 0.3307
## Precision for JDATE 0.4648 0.0592
## Precision for YEAR 1.1849 0.5434
## 0.025quant 0.5quant
## Precision for the Gaussian observations 0.5104 0.5264
## Precision for ID 0.3512 0.6405
## Precision for inla.group(NN, n = 20, method = "cut") 0.3959 0.8179
## Precision for JDATE 0.3609 0.4604
## Precision for YEAR 0.4422 1.0822
## 0.975quant mode
## Precision for the Gaussian observations 0.5426 0.5264
## Precision for ID 1.1074 0.5976
## Precision for inla.group(NN, n = 20, method = "cut") 1.6781 0.7166
## Precision for JDATE 0.5930 0.4512
## Precision for YEAR 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_contall90, s_rip150, s_ripar120, s_riparian90, s_hardfor150)
Recheck for colinearity
Cov.cor3 = as.data.frame(Cov.cor2) %>% select(-c(s_confor90, s_contall90, s_rip150, s_ripar120, s_riparian90, s_hardfor150, s_herbop120))
CI = colldiag(Cov.cor3)
CI ## Condition
## Index Variance Decomposition Proportions
## intercept s_hydro_m s_hardtal150 s_hardtall90 s_herbopen90
## 1 1.000 0.014 0.001 0.001 0.001 0.008
## 2 1.285 0.007 0.049 0.000 0.000 0.085
## 3 1.423 0.001 0.090 0.001 0.001 0.033
## 4 1.854 0.006 0.054 0.001 0.001 0.034
## 5 2.128 0.009 0.098 0.000 0.000 0.351
## 6 2.239 0.011 0.124 0.000 0.000 0.163
## 7 5.145 0.776 0.503 0.001 0.005 0.310
## 8 23.269 0.176 0.081 0.997 0.992 0.016
## s_Paved_m s_consh120 s_Unpaved_m
## 1 0.017 0.006 0.011
## 2 0.022 0.113 0.009
## 3 0.011 0.014 0.096
## 4 0.469 0.001 0.013
## 5 0.071 0.035 0.221
## 6 0.003 0.580 0.064
## 7 0.301 0.242 0.502
## 8 0.106 0.010 0.083nPrior = 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_hardtal150 +
s_hardtall90 + s_herbopen90 + s_Paved_m + s_consh120 + s_Unpaved_m
GModel.2 = 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.2)##
## 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.0246 51.4428 0.6935 54.1609
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept1 1.3549 0.2804 0.8044 1.3549 1.9051 1.3549 0
## s_hydro_m -0.5149 0.0453 -0.6039 -0.5149 -0.4261 -0.5149 0
## s_hardtal150 0.2137 0.0361 0.1429 0.2137 0.2845 0.2137 0
## s_hardtall90 -0.0444 0.0369 -0.1168 -0.0444 0.0280 -0.0444 0
## s_herbopen90 0.0596 0.0165 0.0272 0.0596 0.0920 0.0596 0
## s_Paved_m 0.0747 0.0253 0.0250 0.0747 0.1243 0.0747 0
## s_consh120 0.0300 0.0165 -0.0025 0.0300 0.0624 0.0300 0
## s_Unpaved_m 0.1629 0.1036 -0.0404 0.1629 0.3661 0.1629 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
## Precision for the Gaussian observations 0.5134 0.0080
## Precision for ID 0.6461 0.1883
## Precision for inla.group(NN, n = 20, method = "cut") 0.7553 0.2770
## Precision for JDATE 0.4559 0.0571
## Precision for YEAR 1.1680 0.5560
## 0.025quant 0.5quant
## Precision for the Gaussian observations 0.4978 0.5134
## Precision for ID 0.3444 0.6242
## Precision for inla.group(NN, n = 20, method = "cut") 0.3415 0.7140
## Precision for JDATE 0.3517 0.4533
## Precision for YEAR 0.4488 1.0490
## 0.975quant mode
## Precision for the Gaussian observations 0.5293 0.5134
## Precision for ID 1.0774 0.5823
## Precision for inla.group(NN, n = 20, method = "cut") 1.4131 0.6346
## Precision for JDATE 0.5756 0.4492
## Precision for YEAR 2.5790 0.8541
##
## Expected number of effective parameters(std dev): 140.69(0.00)
## Number of equivalent replicates : 59.92
##
## Deviance Information Criterion (DIC) ...............: 29685.93
## Deviance Information Criterion (DIC, saturated) ....: 8574.11
## Effective number of parameters .....................: 140.69
##
## Watanabe-Akaike information criterion (WAIC) ...: 29694.16
## Effective number of parameters .................: 144.79
##
## Marginal log-Likelihood: -15142.55
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computed5 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_hardtal150 +
s_hardtall90 + s_herbopen90 + s_Paved_m + s_consh120 + 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.7463 36.4247 0.8557 39.0268
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept1 -1.4338 0.1424 -1.7134 -1.4338 -1.1545 -1.4338 0
## s_hydro_m -0.5035 0.0459 -0.5937 -0.5035 -0.4134 -0.5035 0
## s_hardtal150 0.3122 0.0359 0.2418 0.3122 0.3826 0.3122 0
## s_hardtall90 -0.0537 0.0370 -0.1263 -0.0537 0.0188 -0.0537 0
## s_herbopen90 0.1034 0.0164 0.0711 0.1034 0.1356 0.1034 0
## s_Paved_m 0.0791 0.0257 0.0287 0.0791 0.1295 0.0791 0
## s_consh120 0.0559 0.0166 0.0233 0.0559 0.0885 0.0559 0
## s_Unpaved_m 0.1859 0.1053 -0.0208 0.1859 0.3925 0.1859 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.4967 0.0077 0.4817 0.4967
## Precision for ID 0.6390 0.1893 0.3457 0.6129
## Precision for JDATE 0.4373 0.0547 0.3378 0.4346
## Precision for YEAR 1.1740 0.5479 0.4487 1.0625
## 0.975quant mode
## Precision for the Gaussian observations 0.5121 0.4966
## Precision for ID 1.0839 0.5640
## Precision for JDATE 0.5528 0.4300
## Precision for YEAR 2.5483 0.8730
##
## Expected number of effective parameters(std dev): 134.96(0.00)
## Number of equivalent replicates : 62.46
##
## Deviance Information Criterion (DIC) ...............: 29958.72
## Deviance Information Criterion (DIC, saturated) ....: 8567.02
## Effective number of parameters .....................: 134.96
##
## Watanabe-Akaike information criterion (WAIC) ...: 29963.40
## Effective number of parameters .................: 136.64
##
## Marginal log-Likelihood: -15268.41
## 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_hardtal150 +
s_hardtall90 + s_herbopen90 + s_Paved_m + s_consh120 + 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
## 2.1026 54.8417 0.6275 57.5719
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept1 1.4343 0.4431 0.5643 1.4343 2.3037 1.4343 0
## s_hydro_m 0.5901 0.0322 0.5268 0.5901 0.6533 0.5901 0
## s_hardtal150 -0.0043 0.0719 -0.1455 -0.0043 0.1367 -0.0043 0
## s_hardtall90 0.3232 0.0716 0.1828 0.3232 0.4636 0.3232 0
## s_herbopen90 0.1364 0.0325 0.0725 0.1364 0.2002 0.1364 0
## s_Paved_m 0.8129 0.0317 0.7507 0.8129 0.8751 0.8129 0
## s_consh120 0.1863 0.0323 0.1228 0.1863 0.2497 0.1863 0
## s_Unpaved_m 0.9655 0.0323 0.9022 0.9655 1.0288 0.9655 0
##
## Random effects:
## Name Model
## inla.group(NN, n = 20, method = "cut") RW1 model
##
## Model hyperparameters:
## mean sd
## Precision for the Gaussian observations 0.1226 0.0006
## Precision for inla.group(NN, n = 20, method = "cut") 0.6585 0.2288
## 0.025quant 0.5quant
## Precision for the Gaussian observations 0.1217 0.1225
## Precision for inla.group(NN, n = 20, method = "cut") 0.3014 0.6300
## 0.975quant mode
## Precision for the Gaussian observations 0.1242 0.1219
## Precision for inla.group(NN, n = 20, method = "cut") 1.1892 0.5718
##
## Expected number of effective parameters(std dev): 12.09(0.00)
## Number of equivalent replicates : 697.32
##
## Deviance Information Criterion (DIC) ...............: 40039.72
## Deviance Information Criterion (DIC, saturated) ....: 6795.12
## Effective number of parameters .....................: 12.09
##
## Watanabe-Akaike information criterion (WAIC) ...: 40036.52
## Effective number of parameters .................: 8.614
##
## Marginal log-Likelihood: -20091.96
## 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_hardtal150 +
s_hardtall90 + s_herbopen90 + s_Paved_m + s_consh120 + 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.1880 15.8937 0.6762 17.7580
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept1 -2.1054 0.0285 -2.1614 -2.1054 -2.0495 -2.1054 0
## s_hydro_m 0.5989 0.0293 0.5413 0.5989 0.6564 0.5989 0
## s_hardtal150 0.1502 0.0642 0.0242 0.1502 0.2761 0.1502 0
## s_hardtall90 0.3433 0.0643 0.2170 0.3433 0.4695 0.3433 0
## s_herbopen90 0.2234 0.0290 0.1665 0.2234 0.2803 0.2234 0
## s_Paved_m 0.8156 0.0288 0.7590 0.8156 0.8722 0.8156 0
## s_consh120 0.2554 0.0290 0.1984 0.2553 0.3122 0.2554 0
## s_Unpaved_m 0.9605 0.0294 0.9029 0.9605 1.0181 0.9605 0
##
## The model has no random effects
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant
## Precision for the Gaussian observations 0.1467 0.0023 0.1421 0.1467
## 0.975quant mode
## Precision for the Gaussian observations 0.151 0.1469
##
## Expected number of effective parameters(std dev): 8.008(0.00)
## Number of equivalent replicates : 1052.70
##
## Deviance Information Criterion (DIC) ...............: 40120.36
## Deviance Information Criterion (DIC, saturated) ....: 8459.18
## Effective number of parameters .....................: 8.008
##
## Watanabe-Akaike information criterion (WAIC) ...: 40119.51
## Effective number of parameters .................: 7.122
##
## Marginal log-Likelihood: -20122.50
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computed6 Results
Compare Model Metrics. The two global models look comparable, but, dropping the one non-significant covariate increases both the DIC and WAIC slightly. The Linear model without any spatial or grouping effects performs the worst.
Models = c("Global 1", "Global 2",
"Non-Spatial",
"No ID", "Linear Model")
#Deviance information criteria
S1.DICs = c(GModel.1$dic$dic, GModel.2$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, GModel.2$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 |
| Global 2 | 29685.93 | 29694.16 |
| Non-Spatial | 29958.72 | 29963.4 |
| No ID | 40039.72 | 40036.52 |
| Linear Model | 40120.36 | 40119.51 |
Global Model 1 (GModel.2) Fixed Effect Estimates
FE.table = GModel.1$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 | 0.7690793 | 0.2798142 | 0.2197096 | 1.3179905 |
| s_confor90 | 0.2063560 | 0.1516977 | -0.0914778 | 0.5039412 |
| s_hydro_m | -0.4861302 | 0.0449093 | -0.5743023 | -0.3980317 |
| s_contall90 | 0.0191078 | 0.1365799 | -0.2490446 | 0.2870364 |
| s_rip150 | 0.1467073 | 0.0810817 | -0.0124833 | 0.3057651 |
| s_ripar120 | -0.0572668 | 0.1265230 | -0.3056742 | 0.1909331 |
| s_hardtal150 | 0.1167549 | 0.0861590 | -0.0524044 | 0.2857729 |
| s_hardfor150 | 0.0611222 | 0.0806178 | -0.0971577 | 0.2192700 |
| s_riparian90 | 0.0809545 | 0.0671145 | -0.0508139 | 0.2126129 |
| s_hardtall90 | 0.0110246 | 0.0366971 | -0.0610241 | 0.0830132 |
| s_herbopen90 | 0.1031978 | 0.0487429 | 0.0074990 | 0.1988166 |
| s_Paved_m | 0.0746001 | 0.0250144 | 0.0254883 | 0.1236708 |
| s_herbop120 | -0.0271608 | 0.0490716 | -0.1235049 | 0.0691028 |
| s_consh120 | -0.0300504 | 0.0522328 | -0.1326009 | 0.0724146 |
| s_Unpaved_m | 0.2015993 | 0.1029677 | -0.0005610 | 0.4035910 |
ID Effect by Marten ID
mic.d = as.data.frame(GModel.1$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(GModel.1$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(GModel.1$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(GModel.1$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))
7 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_confor90 = train.data[,"s_confor90"],
s_hydro_m = train.data[,"s_hydro_m"],
s_contall90 = train.data[,"s_contall90"],
s_rip150 = train.data[,"s_rip150"],
s_ripar120 = train.data[,"s_ripar120"],
s_hardtal150 = train.data[,"s_hardtal150"],
s_hardfor150 = train.data[,"s_hardfor150"],
s_riparian90 = train.data[,"s_riparian90"],
s_hardtall90 = train.data[,"s_hardtall90"],
s_herbopen90 = train.data[,"s_herbopen90"],
s_Paved_m = train.data[,"s_Paved_m"],
s_herbop120 = train.data[,"s_herbop120"],
s_consh120 = train.data[,"s_consh120"],
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_rip150 = test.data[,"s_rip150"],
s_ripar120 = test.data[,"s_ripar120"],
s_hardtal150 = test.data[,"s_hardtal150"],
s_hardfor150 = test.data[,"s_hardfor150"],
s_riparian90 = test.data[,"s_riparian90"],
s_hardtall90 = test.data[,"s_hardtall90"],
s_herbopen90 = test.data[,"s_herbopen90"],
s_Paved_m = test.data[,"s_Paved_m"],
s_herbop120 = test.data[,"s_herbop120"],
s_consh120 = test.data[,"s_consh120"],
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_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
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.668, df = 1684, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.6789913 0.7272152
## sample estimates:
## cor
## 0.7039135#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.2084358#Mean Absolute Error
mae = function(error)
{
mean(abs(error))
}
mae(test.data$Error)## [1] 0.1383627#Plot
range(test.data$RASTERVALU)## [1] 0.000324589 0.988708973range(test.data$Prediction) ## [1] 0.000724005 0.980441886#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))