Pig Spatiotemporal Movement Model
1 Data Preperation
1.1 Political Boundaries
States = map("state",
fill = TRUE,
plot = FALSE)
IDs = sapply(strsplit(States$names, ":"), function(x) x[1])
LL84 = "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
StatesP = map2SpatialPolygons(States, IDs = IDs,
proj4string = CRS(projection(LL84)))
#Add a dataframe
pid = sapply(slot(StatesP, "polygons"),
function(x) slot(x, "ID"))
p.df = data.frame( ID=1:length(StatesP),
row.names = pid)
States = SpatialPolygonsDataFrame(StatesP, p.df)
States = spTransform(States, LL84)
Domain = subset(States, ID == 21)
Domain$Name = str_cap_words(rownames(Domain@data))
DomLLU = gUnaryUnion(Domain)
#Rasterized version
Ras = raster(res = 0.02, ext = extent(DomLLU),
crs = proj4string(DomLLU))
Domain.r = rasterize(DomLLU, Ras,
field = 0,
background = NA)
#Point grid version
Grd.pnt = rasterToPoints(Domain.r, spatial = TRUE)
Grd.pnt@data = Grd.pnt@data %>%
mutate(Long = Grd.pnt@coords[,1],
Lat = Grd.pnt@coords[,2]) %>%
select(-layer)1.2 Load and Prep Telemetry Data
wpmove = read.csv("./ST_060119/MovementModeling_Final4.csv", header=T)
wpmove = wpmove %>%
filter(id != "23192" & id != "20131") %>% #Remove UP pigs
select(tspan, Dvlpd_A, Rprn_Ar, Open_Ar, Hrdms_A, PCI_2, PCI_5, Edge1_3, HOURLYDRYBULBTEMPC, HOURLYStationPressure,
cohort, hour, day, jdate, month, year, id, x, y, state, dt)
wpmove %>%
group_by(state) %>%
summarise(AvgT = median(tspan, na.rm=T))## # A tibble: 3 x 2
## state AvgT
## <int> <dbl>
## 1 1 203.
## 2 2 234.
## 3 3 253.Scaled.set = as.data.frame(scale(wpmove[,1:10], scale=T, center=T))
apply(Scaled.set, 2, range)## tspan Dvlpd_A Rprn_Ar Open_Ar Hrdms_A PCI_2
## [1,] -1.235169 -0.1977294 -1.075746 -0.4295997 -0.09612832 -0.475362
## [2,] 9.655976 17.3070993 1.203094 3.7447104 21.18692598 2.839339
## PCI_5 Edge1_3 HOURLYDRYBULBTEMPC HOURLYStationPressure
## [1,] -0.3766494 -0.4633794 -3.702357 -4.393409
## [2,] 3.8110604 3.8251009 2.272631 2.991898names(Scaled.set) = paste("s_", names(wpmove[,1:10]), sep="")
wpmove = cbind(wpmove, Scaled.set)
### Create training and testing dataset - 20% from each pig ############
# Partition data for training and testing
set.seed(123)
wpmovetest = wpmove %>%
group_by(id) %>%
sample_frac(0.20, replace=F)
wpmovetest = as.data.frame(wpmovetest)
dim(wpmovetest) ## [1] 4693 31wpmovetrain = anti_join(wpmove,wpmovetest) ## Joining, by = c("tspan", "Dvlpd_A", "Rprn_Ar", "Open_Ar", "Hrdms_A", "PCI_2", "PCI_5", "Edge1_3", "HOURLYDRYBULBTEMPC", "HOURLYStationPressure", "cohort", "hour", "day", "jdate", "month", "year", "id", "x", "y", "state", "dt", "s_tspan", "s_Dvlpd_A", "s_Rprn_Ar", "s_Open_Ar", "s_Hrdms_A", "s_PCI_2", "s_PCI_5", "s_Edge1_3", "s_HOURLYDRYBULBTEMPC", "s_HOURLYStationPressure")dim(wpmovetrain) ## [1] 18780 311.3 Set Time Steps
wpmovetrain %>% #Observations per pig
group_by(id) %>%
summarise(Count = length(id))## # A tibble: 9 x 2
## id Count
## <fct> <int>
## 1 36635 2210
## 2 38303 864
## 3 38306 5711
## 4 38307 874
## 5 38308 806
## 6 38309 95
## 7 38312 2154
## 8 38313 1566
## 9 midland 4500wpmovetrain$Date = as.POSIXct(as.character(wpmovetrain$dt), #Convert to date
tz = "America/New_York",
format = "%m/%d/%Y %H:%M")
Pig.lvls = unique(wpmovetrain$id)
for(i in 1:length(Pig.lvls)){ #For all pigs
tmp.pig = wpmovetrain %>% filter(id == Pig.lvls[i]) #take one at a time
tmp.pig = arrange(tmp.pig, Date) #ensure that date stamps are ordered
tmp.pig$TStep = 1:nrow(tmp.pig) #Number the time steps as intergers
tmp.pig$Pig = i #Keep track of loop
if(i == 1){Out.pig = tmp.pig} #reassemble dataset
else{Out.pig = rbind(Out.pig, tmp.pig)}
}
wpmovetrain = Out.pig1.4 Locations Map
nProj = "+proj=utm +zone=16 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"
nProj2 = "+proj=utm +zone=16 +datum=NAD83 +units=km +no_defs +ellps=GRS80 +towgs84=0,0,0"
Pnts = SpatialPointsDataFrame(wpmovetrain[,c("x","y")], wpmovetrain)
proj4string(Pnts) = nProj
Pnts = spTransform(Pnts, nProj2)
DomainP = spTransform(Domain, proj4string(Pnts))
plot(DomainP)
plot(Pnts, add=T)
1.5 Nearest Neighbor
Find each pigs nearest neighbor by movement state.
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.60-1 (nickname: 'Swinging Sixties')
## For an introduction to spatstat, type 'beginner'##
## Attaching package: 'spatstat'## The following object is masked _by_ '.GlobalEnv':
##
## im## The following object is masked from 'package:data.table':
##
## shift## The following object is masked from 'package:dismo':
##
## domain## The following object is masked from 'package:PresenceAbsence':
##
## auc## The following object is masked from 'package:MASS':
##
## area## The following object is masked from 'package:lattice':
##
## panel.histogram## The following objects are masked from 'package:raster':
##
## area, rotate, shiftNN.pnt = Pnts
By = as.factor(NN.pnt$state)
P.pp = as.ppp(NN.pnt)
NN = as.data.frame(nndist(P.pp, by=By, k = 1))
Pnts$nS1 = round(NN[,"1"], 3)
Pnts$nS2 = round(NN[,"2"], 3)
Pnts$nS3 = round(NN[,"3"], 3)1.6 Spatial Mesh
Define spatial domain.
max.edge = 1 #Make the outer edge length (km)
bound.outer = 10 #Outer extension (km)
Hull = gBuffer(
gConvexHull(Pnts),
width = 5)
plot(Hull)
bdry = inla.sp2segment(Hull) #Formatting boundary
mesh = inla.mesh.2d(boundary = bdry,
loc = Pnts, #Fit to pig locations
max.edge = c(1, 2)*max.edge,
cutoff = 0.3,
min.angle = 30,
offset = c(6, bound.outer))Result
mesh$n #number of nodes## [1] 5643plot(mesh, rgl=F)
With Observations overlaid.
mesh$n #number of nodes## [1] 5643plot(mesh, rgl=F)
plot(Pnts, col="red", add=T)
1.7 Projection Matricies
Relating observations to their corresponding loction on the mesh.
locs = cbind(Pnts@coords[,1], Pnts@coords[,2]) #point locations
A.mx1 = inla.spde.make.A(mesh, #the mesh
alpha = 2, #default setting
loc=locs) #our locations
A.mx3 = A.mx2 = A.mx1
spde = inla.spde2.pcmatern(mesh,
prior.range = c(10, 0.5),
prior.sigma = c(1, 0.5))
#Create index to track locations of mesh nodes in the model.
field1 = inla.spde.make.index("field1", spde$n.spde)
field2 = inla.spde.make.index("field2", spde$n.spde)
field3 = inla.spde.make.index("field3", spde$n.spde)
field1x = inla.spde.make.index("field1x", spde$n.spde)
field1xx = inla.spde.make.index("field1xx", spde$n.spde)
field2x = inla.spde.make.index("field2x", spde$n.spde) 1.8 Organize Data
Creating a unique list object for each movement state.
pnts.df = Pnts@data
pnts.df$Point.IID = 1:nrow(pnts.df)
#Movement State 1
HR1.lst = list(c(field1, #Spatial index for State 1
list(intercept1 = 1)), #intercept
list(Dvlpd_A1 = round(pnts.df[,"s_Dvlpd_A"], 3), #Covarites and effects below
Rprn_Ar1 = round(pnts.df[,"s_Rprn_Ar"], 3),
Open_Ar1 = round(pnts.df[,"s_Open_Ar"], 3),
Hrdms_A1 = round(pnts.df[,"s_Hrdms_A"], 3),
PCI_21 = round(pnts.df[,"s_PCI_2"], 3),
PCI_51 = round(pnts.df[,"s_PCI_5"], 3),
Edge1_31 = round(pnts.df[,"s_Edge1_3"], 3),
HOURLYDRYBULBTEMPC1 = round(pnts.df[,"s_HOURLYDRYBULBTEMPC"], 3),
HOURLYStationPressure1 = round(pnts.df[,"s_HOURLYStationPressure"], 3),
NN = pnts.df[,"nS1"],
hour1 = pnts.df[,"hour"],
day1 = pnts.df[,"jdate"],
month1 = pnts.df[,"month"],
year = pnts.df[,"year"],
cohort1 = pnts.df[,"cohort"],
State = as.integer(pnts.df[,"state"]),
TStep1 = as.integer(pnts.df[,"TStep"]),
pig1 = pnts.df[,"Pig"],
pig = pnts.df[,"Pig"],
Point.IID = pnts.df[,"Point.IID"]))
pnts.df$State1 = ifelse(pnts.df$state == 1, 1, 0) #if state 1, then 1, else zero (probabilty of State 1)
sum(pnts.df$State1)## [1] 4531#Test version
TS1.stk = inla.stack(data = list(Y = pnts.df$State1), #Single response to model State 1 by itself
A = list(A.mx1, 1), #Spatial matrix
effects = HR1.lst,
tag = "T1.0")
#Joint model version
S1.stk = inla.stack(data = list(Y = cbind(pnts.df$State1, NA, NA)), #NA are placeholders for State 2 and State 3 responses.
A = list(A.mx1, 1),
effects = HR1.lst,
tag = "s1.0")
###State 2
HR2.lst = list(c(field2, #Spatial index for State 2
field1x, #Copy of State 1 spatial effect estimated in Model Tier 1
list(intercept2 = 1)), #intercept
list(Dvlpd_A2 = round(pnts.df[,"s_Dvlpd_A"], 3),
Rprn_Ar2 = round(pnts.df[,"s_Rprn_Ar"], 3),
Open_Ar2 = round(pnts.df[,"s_Open_Ar"], 3),
Hrdms_A2 = round(pnts.df[,"s_Hrdms_A"], 3),
PCI_22 = round(pnts.df[,"s_PCI_2"], 3),
PCI_52 = round(pnts.df[,"s_PCI_5"], 3),
Edge1_32 = round(pnts.df[,"s_Edge1_3"], 3),
HOURLYDRYBULBTEMPC2 = round(pnts.df[,"s_HOURLYDRYBULBTEMPC"], 3),
HOURLYStationPressure2 = round(pnts.df[,"s_HOURLYStationPressure"], 3),
NN = pnts.df[,"nS2"],
hour2 = pnts.df[,"hour"],
day2 = pnts.df[,"jdate"],
month2 = pnts.df[,"month"],
year = pnts.df[,"year"],
cohort2 = pnts.df[,"cohort"],
State = as.integer(pnts.df[,"state"]),
TStep2 = as.integer(pnts.df[,"TStep"]),
pig2 = pnts.df[,"Pig"],
pig = pnts.df[,"Pig"],
Point.IID = pnts.df[,"Point.IID"]))
pnts.df$State2 = ifelse(pnts.df$state == 2, 1, 0) #State 2 set to 1, States 1 and 3 as zero.
sum(pnts.df$State2)## [1] 9279#Test version
TS2.stk = inla.stack(data = list(Y = pnts.df$State2),
A = list(A.mx2, 1),
effects = HR2.lst,
tag = "T2.0")
S2.stk = inla.stack(data = list(Y = cbind(NA, pnts.df$State2, NA)), #NA's for State 1 and 3 placeholders
A = list(A.mx2, 1),
effects = HR2.lst,
tag = "s2.0")
###State 3
HR3.lst = list(c(field3, #Spatial index for State 3
field1xx, #Copy of field estimated for State 1
field2x, #Copy of spatial field estimated for State 2
list(intercept3 = 1)),
list(Dvlpd_A3 = round(pnts.df[,"s_Dvlpd_A"], 3),
Rprn_Ar3 = round(pnts.df[,"s_Rprn_Ar"], 3),
Open_Ar3 = round(pnts.df[,"s_Open_Ar"], 3),
Hrdms_A3 = round(pnts.df[,"s_Hrdms_A"], 3),
PCI_23 = round(pnts.df[,"s_PCI_2"], 3),
PCI_53 = round(pnts.df[,"s_PCI_5"], 3),
Edge1_33 = round(pnts.df[,"s_Edge1_3"], 3),
HOURLYDRYBULBTEMPC3 = round(pnts.df[,"s_HOURLYDRYBULBTEMPC"], 3),
HOURLYStationPressure3 = round(pnts.df[,"s_HOURLYStationPressure"], 3),
NN = pnts.df[,"nS3"],
hour3 = pnts.df[,"hour"],
day3 = pnts.df[,"jdate"],
month3 = pnts.df[,"month"],
year = pnts.df[,"year"],
cohort3 = pnts.df[,"cohort"],
State = as.integer(pnts.df[,"state"]),
TStep3 = as.integer(pnts.df[,"TStep"]),
pig3 = pnts.df[,"Pig"],
pig = pnts.df[,"Pig"],
Point.IID = pnts.df[,"Point.IID"]))
pnts.df$State3 = ifelse(pnts.df$state == 3, 1, 0) #State 3 set to 1, States 1 and 2 as zero
sum(pnts.df$State3) ## [1] 4970#Test version
TS3.stk = inla.stack(data = list(Y = pnts.df$State3),
A = list(A.mx3, 1),
effects = HR3.lst,
tag = "T3.0")
S3.stk = inla.stack(data = list(Y = cbind(NA, NA, pnts.df$State3)),
A = list(A.mx3, 1),
effects = HR3.lst,
tag = "s3.0")1.9 Combine Stacks
Data set for joint model, combines the unique stack created for each movement state above. Note that the covariates have different names, so model covariates are model level (movement state) specific.
Joint.stk = inla.stack(S1.stk, S2.stk, S3.stk)
#Copies of data o run on HPCC
#save(list=c("Joint.stk", "spde"), file="./ST_060119/HPC/PigM_060519.RData") #bundled for HPCC
#save(list=c("Joint.stk", "spde"), file="./ST_060119/HPC_061719/PigM_061719.RData")
#save(list=c("TS1.stk", "TS2.stk","TS3.stk", "spde"), file="./ST_060119/HPC_061719/BaseMod_061919.RData")2 Base Models
Modeling each movement state seperately for later comparison.
2.1 State 1 (Base Model)
pcprior1 = list(prec = list(prior="pc.prec", param = c(1, 0.001)))
pcprior2 = list(prec = list(prior="pc.prec", param = c(3, 0.1)))
LCPrior = list(theta=list(prior = "normal", param=c(0, 5)))
TPrior = list(theta=list(prior = "normal", param=c(0, 5)))
LCPrior2 = list(theta=list(prior = "normal", param=c(0, 3)))
h.hyper = list(theta = list(prior="pc.cor1", param=c(0, 3)))
Frm1 = Y ~ -1 + intercept1 +
f(field1,
model = spde) +
f(NN,
model="rw1",
constr=TRUE,
scale.model = TRUE,
hyper=TPrior) +
f(pig,
constr=TRUE,
model="iid",
hyper=LCPrior) +
f(year,
constr=TRUE,
model="iid",
hyper=LCPrior) +
f(cohort1,
constr=TRUE,
model="iid",
hyper=LCPrior) +
f(day1,
model="iid",
constr=TRUE,
hyper=LCPrior) +
f(hour1,
model="rw1",
constr=TRUE,
scale.model = TRUE,
hyper=TPrior) +
f(TStep1,
model="ar1",
constr=TRUE,
replicate = pig1,
hyper = h.hyper) +
Dvlpd_A1 +
Rprn_Ar1 +
Open_Ar1 +
Hrdms_A1 +
PCI_21 +
PCI_51 +
Edge1_31 +
HOURLYDRYBULBTEMPC1 +
HOURLYStationPressure1
Model.S1 = inla(Frm1,
data = inla.stack.data(TS1.stk, spde=spde),
family = "binomial",
verbose = TRUE,
control.fixed = list(prec = 0.001,
prec.intercept = 0.0001),
control.predictor = list(
A = inla.stack.A(TS1.stk),
compute = TRUE,
link = 1),
control.inla = list(strategy="gaussian",
int.strategy = "eb"),
#control.mode = list(restart = TRUE, theta = theta1),
control.results = list(return.marginals.random = TRUE,
return.marginals.predictor = TRUE),
control.compute=list(dic = TRUE, cpo = TRUE, waic = TRUE))2.2 State 2 (Base Model)
Frm2 = Y ~ -1 + intercept2 +
f(field2,
model = spde) +
f(NN,
model="rw1",
constr=TRUE,
scale.model = TRUE,
hyper=TPrior) +
f(pig,
constr=TRUE,
model="iid",
hyper=LCPrior) +
f(year,
constr=TRUE,
model="iid",
hyper=LCPrior) +
f(cohort2,
constr=TRUE,
model="iid",
hyper=LCPrior) +
f(day2,
model="iid",
constr=TRUE,
hyper=LCPrior) +
f(hour2,
model="rw1",
constr=TRUE,
scale.model = TRUE,
hyper=TPrior) +
f(TStep2,
model="ar1",
constr=TRUE,
replicate = pig2,
hyper = h.hyper) +
Dvlpd_A2 +
Rprn_Ar2 +
Open_Ar2 +
Hrdms_A2 +
PCI_22 +
PCI_52 +
Edge1_32 +
HOURLYDRYBULBTEMPC2 +
HOURLYStationPressure2
Model.S2 = inla(Frm2,
data = inla.stack.data(TS2.stk, spde=spde),
family = "binomial",
verbose = TRUE,
control.fixed = list(prec = 0.001,
prec.intercept = 0.0001),
control.predictor = list(
A = inla.stack.A(TS2.stk),
compute = TRUE,
link = 1),
control.inla = list(strategy="gaussian",
int.strategy = "eb"),
#control.mode = list(restart = TRUE, theta = theta1),
control.results = list(return.marginals.random = TRUE,
return.marginals.predictor = TRUE),
control.compute=list(dic = TRUE, cpo = TRUE, waic = TRUE))2.3 State 3 (Base Model)
Frm3 = Y ~ -1 + intercept3 +
f(field3,
model = spde) +
f(NN,
model="rw1",
constr=TRUE,
scale.model = TRUE,
hyper=TPrior) +
f(pig,
constr=TRUE,
model="iid",
hyper=LCPrior) +
f(year,
constr=TRUE,
model="iid",
hyper=LCPrior) +
f(cohort3,
constr=TRUE,
model="iid",
hyper=LCPrior) +
f(day3,
model="iid",
constr=TRUE,
hyper=LCPrior) +
f(hour3,
model="rw1",
constr=TRUE,
scale.model = TRUE,
hyper=TPrior) +
f(TStep3,
model="ar1",
constr=TRUE,
replicate = pig3,
hyper = h.hyper) +
Dvlpd_A3 +
Rprn_Ar3 +
Open_Ar3 +
Hrdms_A3 +
PCI_23 +
PCI_53 +
Edge1_33 +
HOURLYDRYBULBTEMPC3 +
HOURLYStationPressure3
Model.S3 = inla(Frm3,
data = inla.stack.data(TS3.stk, spde=spde),
family = "binomial",
verbose = TRUE,
control.fixed = list(prec = 0.001,
prec.intercept = 0.0001),
control.predictor = list(
A = inla.stack.A(TS3.stk),
compute = TRUE,
link = 1),
control.inla = list(strategy="gaussian",
int.strategy = "eb"),
#control.mode = list(restart = TRUE, theta = theta1),
control.results = list(return.marginals.random = TRUE,
return.marginals.predictor = TRUE),
control.compute=list(dic = TRUE, cpo = TRUE, waic = TRUE))2.4 Check Base models
load("~/Michigan_State/Steve/ST_060119/HPC_061719/Base_RES_061919B.RData")
GetMets(Model.S1)## Metric Tier.1
## 1 DIC 6272.1211
## 2 WAIC 6263.3147
## 3 lCPO 0.9073GetMets(Model.S2)## Metric Tier.1
## 1 DIC 11600.4019
## 2 WAIC 11036.4863
## 3 lCPO 0.7645GetMets(Model.S3)## Metric Tier.1
## 1 DIC 13369.3451
## 2 WAIC 12915.5064
## 3 lCPO 0.36983 Joint Models
3.1 Random Effect Model
Joint model for all movemnet states, but, with random effects only.
pcprior1 = list(prec = list(prior="pc.prec", param = c(1, 0.001))) #Different prior explored.
pcprior2 = list(prec = list(prior="pc.prec", param = c(3, 0.1)))
LCPrior = list(theta=list(prior = "normal", param=c(0, 5)))
TPrior = list(theta=list(prior = "normal", param=c(0, 5)))
LCPrior2 = list(theta=list(prior = "normal", param=c(0, 3)))
h.hyper = list(theta = list(prior="pc.cor1", param=c(0, 3)))
#List all level specifc covarites
Frm1 = Y ~ -1 + intercept1 + #Intercepts
intercept2 +
intercept3 +
f(field1, #State 1 spatial
model = spde) +
f(field2, #State 2 spatial
model = spde) +
f(field3, #State 3 spatial
model = spde) +
f(field1x, #Copy of State 1 spatial for State 2
copy = "field1",
fixed = FALSE) +
f(field1xx, #Copy of State 1 spatial for State 3
copy = "field1",
fixed = FALSE) +
f(field2x, #Copy of State 2 spatial for State 3
copy = "field2",
fixed = FALSE) +
f(NN, #Nearest neighbor
model="rw1",
constr=TRUE,
scale.model = TRUE,
hyper=TPrior) +
f(pig, #Variation due to individual differences (all states)
constr=TRUE,
model="iid",
hyper=LCPrior) +
f(year, #Variation due to years (all states)
constr=TRUE,
model="iid",
hyper=LCPrior) +
f(TStep1, #Time steps for State 1
model="ar1", #Autoregressive Order 1
constr=TRUE, #Force the effect to have a zero mean average (helps avoid confounding/colinearity)
replicate = pig1, #Repeat analysis seperatly for each unique pig
hyper = h.hyper) + #Flat prior
f(TStep2, #As baove but for State 2
model="ar1",
constr=TRUE,
replicate = pig2,
hyper = h.hyper) +
f(TStep3, #As above but for State 3
model="ar1",
constr=TRUE,
replicate = pig3,
hyper = h.hyper)
Model1.st = inla(Frm1, #above formula
num.threads = 8, #Number of cores ran on HPC
data = inla.stack.data(Joint.stk, spde=spde), #Data to use
family = c("binomial","binomial","binomial"), #Lieklihoods for each level
verbose = TRUE,
control.fixed = list(prec = 0.001, #intercept/fixed effect prior (flat/default)
prec.intercept = 0.0001),
control.predictor = list( #calculate fitted values
A = inla.stack.A(Joint.stk),
compute = TRUE,
link = 1),
control.inla = list(strategy="gaussian", #Use mode instead of mean (speeds things up)
int.strategy = "eb"),
control.results = list(return.marginals.random = TRUE,
return.marginals.predictor = TRUE),
control.compute=list(dic = TRUE, cpo = TRUE, waic = TRUE))Load Results from HPCC
load("~/Michigan_State/Steve/ST_060119/HPC_061719/PigM_stRES_061719.RData")3.2 Full Joint Model1
With Autoregressive term replicated for each individual.
pcprior1 = list(prec = list(prior="pc.prec", param = c(1, 0.001)))
pcprior2 = list(prec = list(prior="pc.prec", param = c(3, 0.1)))
LCPrior = list(theta=list(prior = "normal", param=c(0, 5)))
TPrior = list(theta=list(prior = "normal", param=c(0, 5)))
LCPrior2 = list(theta=list(prior = "normal", param=c(0, 3)))
h.hyper = list(theta = list(prior="pc.cor1", param=c(0, 3)))
Frm1 = Y ~ -1 + intercept1 +
intercept2 +
intercept3 +
f(field1,
model = spde) +
f(field2,
model = spde) +
f(field3,
model = spde) +
f(field1x,
copy = "field1",
fixed = FALSE) +
f(field1xx,
copy = "field1",
fixed = FALSE) +
f(field2x,
copy = "field2",
fixed = FALSE) +
f(NN,
model="rw1",
constr=TRUE,
scale.model = TRUE,
hyper=TPrior) +
f(pig,
constr=TRUE,
model="iid",
hyper=LCPrior) +
f(year,
constr=TRUE,
model="iid",
hyper=LCPrior) +
f(cohort1,
constr=TRUE,
model="iid",
hyper=LCPrior) +
f(cohort2,
constr=TRUE,
model="iid",
hyper=LCPrior) +
f(cohort3,
constr=TRUE,
model="iid",
hyper=LCPrior) +
f(day1,
model="iid",
constr=TRUE,
hyper=LCPrior) +
f(hour1,
model="rw1",
constr=TRUE,
scale.model = TRUE,
hyper=TPrior) +
f(day2,
model="iid",
constr=TRUE,
hyper=LCPrior) +
f(hour2,
model="rw1",
constr=TRUE,
scale.model = TRUE,
hyper=TPrior) +
f(day3,
model="iid",
constr=TRUE,
hyper=LCPrior) +
f(hour3,
model="rw1",
constr=TRUE,
scale.model = TRUE,
hyper=TPrior) +
f(TStep1,
model="ar1",
constr=TRUE,
replicate = pig1,
hyper = h.hyper) +
f(TStep2,
model="ar1",
constr=TRUE,
replicate = pig2,
hyper = h.hyper) +
f(TStep3,
model="ar1",
constr=TRUE,
replicate = pig3,
hyper = h.hyper) +
Dvlpd_A1 + Dvlpd_A2 + Dvlpd_A3 +
Rprn_Ar1 + Rprn_Ar2 + Rprn_Ar3 +
Open_Ar1 + Open_Ar2 + Open_Ar3 +
Hrdms_A1 + Hrdms_A2 + Hrdms_A3 +
PCI_21 + PCI_22 + PCI_23 +
PCI_51 + PCI_52 + PCI_53 +
Edge1_31 + Edge1_32 + Edge1_33 +
HOURLYDRYBULBTEMPC1 + HOURLYDRYBULBTEMPC2 + HOURLYDRYBULBTEMPC3 +
HOURLYStationPressure1 + HOURLYStationPressure2 + HOURLYStationPressure3
#theta1 = Model1$internal.summary.hyperpar$mean
theta1 = c(1.12260417, 2.58477356, 1.35086993, 2.38248830, 2.70405057, 1.04526546, 0.33116481, 0.13586773,
-0.85066598, -0.86655625, 0.12612075, 0.11837675, 0.04338175, 0.48530462, 0.48937685, 1.42749601,
0.44789598, -3.44024949, 4.16518162, -2.10516022, 2.80414702, -1.61869425, 2.09393767, 0.69071527,
-0.84613170, -0.39377607)
Model1 = inla(Frm1,
data = inla.stack.data(Joint.stk, spde=spde),
family = c("binomial","binomial","binomial"),
verbose = TRUE,
control.fixed = list(prec = 0.001,
prec.intercept = 0.0001),
control.predictor = list(
A = inla.stack.A(Joint.stk),
compute = TRUE,
link = 1),
control.inla = list(strategy="gaussian",
int.strategy = "eb"),
#control.mode = list(restart = TRUE, theta = theta1),
control.results = list(return.marginals.random = TRUE,
return.marginals.predictor = TRUE),
control.compute=list(dic = TRUE, cpo = TRUE, waic = TRUE))Load Results from HPCC
#load("~/Michigan_State/Steve/ST_060119/HPC/PigM_RES_060519.RData")
load("~/Michigan_State/Steve/ST_060119/HPC_061719/PigM_RES_061719.RData") #Updated3.3 Full Joint Model (OU)
Implementing ab Ornstein-Uhlenbeck process term for time correlation (see, https://en.wikipedia.org/wiki/Ornstein%E2%80%93Uhlenbeck_process).
pcprior1 = list(prec = list(prior="pc.prec", param = c(1, 0.001)))
pcprior2 = list(prec = list(prior="pc.prec", param = c(3, 0.1)))
LCPrior = list(theta=list(prior = "normal", param=c(0, 5)))
TPrior = list(theta=list(prior = "normal", param=c(0, 5)))
LCPrior2 = list(theta=list(prior = "normal", param=c(0, 3)))
h.hyper = list(theta = list(prior="pc.cor1", param=c(0, 3)))
Frm1 = Y ~ -1 + intercept1 +
intercept2 +
intercept3 +
f(field1,
model = spde) +
f(field2,
model = spde) +
f(field3,
model = spde) +
f(field1x,
copy = "field1",
fixed = FALSE) +
f(field1xx,
copy = "field1",
fixed = FALSE) +
f(field2x,
copy = "field2",
fixed = FALSE) +
f(NN,
model="rw1",
constr=TRUE,
scale.model = TRUE,
hyper=TPrior) +
f(pig,
constr=TRUE,
model="iid",
hyper=LCPrior) +
f(year,
constr=TRUE,
model="iid",
hyper=LCPrior) +
f(cohort1,
constr=TRUE,
model="iid",
hyper=LCPrior) +
f(cohort2,
constr=TRUE,
model="iid",
hyper=LCPrior) +
f(cohort3,
constr=TRUE,
model="iid",
hyper=LCPrior) +
f(day1,
model="iid",
constr=TRUE,
hyper=LCPrior) +
f(hour1,
model="rw1",
constr=TRUE,
scale.model = TRUE,
hyper=TPrior) +
f(day2,
model="iid",
constr=TRUE,
hyper=LCPrior) +
f(hour2,
model="rw1",
constr=TRUE,
scale.model = TRUE,
hyper=TPrior) +
f(day3,
model="iid",
constr=TRUE,
hyper=LCPrior) +
f(hour3,
model="ou",
constr=TRUE,
scale.model = TRUE,
hyper=TPrior) +
f(TStep1,
model="ou",
constr=TRUE,
replicate = pig1,
hyper = h.hyper) +
f(TStep2,
model="ou",
constr=TRUE,
replicate = pig2,
hyper = h.hyper) +
f(TStep3,
model="ar1",
constr=TRUE,
replicate = pig3,
hyper = h.hyper) +
Dvlpd_A1 + Dvlpd_A2 + Dvlpd_A3 +
Rprn_Ar1 + Rprn_Ar2 + Rprn_Ar3 +
Open_Ar1 + Open_Ar2 + Open_Ar3 +
Hrdms_A1 + Hrdms_A2 + Hrdms_A3 +
PCI_21 + PCI_22 + PCI_23 +
PCI_51 + PCI_52 + PCI_53 +
Edge1_31 + Edge1_32 + Edge1_33 +
HOURLYDRYBULBTEMPC1 + HOURLYDRYBULBTEMPC2 + HOURLYDRYBULBTEMPC3 +
HOURLYStationPressure1 + HOURLYStationPressure2 + HOURLYStationPressure3
#theta1 = Model1$internal.summary.hyperpar$mean
theta1 = c(1.12260417, 2.58477356, 1.35086993, 2.38248830, 2.70405057, 1.04526546, 0.33116481, 0.13586773,
-0.85066598, -0.86655625, 0.12612075, 0.11837675, 0.04338175, 0.48530462, 0.48937685, 1.42749601,
0.44789598, -3.44024949, 4.16518162, -2.10516022, 2.80414702, -1.61869425, 2.09393767, 0.69071527,
-0.84613170, -0.39377607)
Model1.ou = inla(Frm1,
data = inla.stack.data(Joint.stk, spde=spde),
family = c("binomial","binomial","binomial"),
verbose = TRUE,
control.fixed = list(prec = 0.001,
prec.intercept = 0.0001),
control.predictor = list(
A = inla.stack.A(Joint.stk),
compute = TRUE,
link = 1),
control.inla = list(strategy="gaussian",
int.strategy = "eb"),
#control.mode = list(restart = TRUE, theta = theta1),
control.results = list(return.marginals.random = TRUE,
return.marginals.predictor = TRUE),
control.compute=list(dic = TRUE, cpo = TRUE, waic = TRUE))Load Results from HPCC
load("~/Michigan_State/Steve/ST_060119/HPC_061719/PigM_RESOU_061719.RData")3.4 Full Joint Model2
Non spatiotemporal version.
pcprior1 = list(prec = list(prior="pc.prec", param = c(1, 0.001)))
pcprior2 = list(prec = list(prior="pc.prec", param = c(3, 0.1)))
LCPrior = list(theta=list(prior = "normal", param=c(0, 5)))
TPrior = list(theta=list(prior = "normal", param=c(0, 5)))
LCPrior2 = list(theta=list(prior = "normal", param=c(0, 3)))
h.hyper = list(theta = list(prior="pc.cor1", param=c(0, 3)))
Frm1 = Y ~ -1 + intercept1 +
intercept2 +
intercept3 +
# f(field1,
# model = spde) +
# f(field2,
# model = spde) +
# f(field3,
# model = spde) +
# f(field1x,
# copy = "field1",
# fixed = FALSE) +
# f(field1xx,
# copy = "field1",
# fixed = FALSE) +
# f(field2x,
# copy = "field2",
# fixed = FALSE) +
# f(NN,
# model="rw1",
# constr=TRUE,
# scale.model = TRUE,
# hyper=TPrior) +
f(pig,
constr=TRUE,
model="iid",
hyper=LCPrior) +
f(year,
constr=TRUE,
model="iid",
hyper=LCPrior) +
f(cohort1,
constr=TRUE,
model="iid",
hyper=LCPrior) +
f(cohort2,
constr=TRUE,
model="iid",
hyper=LCPrior) +
f(cohort3,
constr=TRUE,
model="iid",
hyper=LCPrior) +
f(day1,
model="iid",
constr=TRUE,
hyper=LCPrior) +
f(hour1,
model="rw1",
constr=TRUE,
scale.model = TRUE,
hyper=TPrior) +
f(day2,
model="iid",
constr=TRUE,
hyper=LCPrior) +
f(hour2,
model="rw1",
constr=TRUE,
scale.model = TRUE,
hyper=TPrior) +
f(day3,
model="iid",
constr=TRUE,
hyper=LCPrior) +
f(hour3,
model="rw1",
constr=TRUE,
scale.model = TRUE,
hyper=TPrior) +
# f(TStep1,
# model="ar1",
## constr=TRUE,
# replicate = pig1,
# hyper = h.hyper) +
# f(TStep2,
# model="ar1",
# constr=TRUE,
# replicate = pig2,
# hyper = h.hyper) +
# f(TStep3,
# model="ar1",
# constr=TRUE,
# replicate = pig3,
# hyper = h.hyper) +
Dvlpd_A1 + Dvlpd_A2 + Dvlpd_A3 +
Rprn_Ar1 + Rprn_Ar2 + Rprn_Ar3 +
Open_Ar1 + Open_Ar2 + Open_Ar3 +
Hrdms_A1 + Hrdms_A2 + Hrdms_A3 +
PCI_21 + PCI_22 + PCI_23 +
PCI_51 + PCI_52 + PCI_53 +
Edge1_31 + Edge1_32 + Edge1_33 +
HOURLYDRYBULBTEMPC1 + HOURLYDRYBULBTEMPC2 + HOURLYDRYBULBTEMPC3 +
HOURLYStationPressure1 + HOURLYStationPressure2 + HOURLYStationPressure3
Model2 = inla(Frm1,
data = inla.stack.data(Joint.stk, spde=spde),
family = c("binomial","binomial","binomial"),
verbose = TRUE,
control.fixed = list(prec = 0.001,
prec.intercept = 0.0001),
control.predictor = list(
A = inla.stack.A(Joint.stk),
compute = TRUE,
link = 1),
control.inla = list(strategy="gaussian",
int.strategy = "eb"),
#control.mode = list(restart = TRUE, theta = theta1),
control.results = list(return.marginals.random = TRUE,
return.marginals.predictor = TRUE),
control.compute=list(dic = TRUE, cpo = TRUE, waic = TRUE))
#save(list=c("Model2"), file="./ST_060119/HPC/Model2_ns_082019.RData") Load Results from prior run.
load("~/Michigan_State/Steve/ST_060119/HPC/Model2_ns_082019.RData")4 Results
4.1 Parsimony Comparison
Compare all models.
Models = c("Model1", "Model2", "Model3",
"Model4", "Model5", "Model6",
"Model7")
Type = c("Base Model (RM-LF)",
"Base Model (LM-MF)",
"Base Model (MM-HF)",
"Joint Random Effects",
"Joint (Full)",
"Joint (Ornstein-Uhlenbeck)",
"Joint (Non Spatiotemporal)")
#Deviance information criteria
S1.DICs = c(round(GetMets(Model.S1)[1,"Tier.1"], 2),
"-",
"-",
round(GetMets(Model1.st)[1,"Tier.1"], 2),
round(GetMets(Model1)[1,"Tier.1"], 2),
round(GetMets(Model1.ou)[1,"Tier.1"], 2),
round(GetMets(Model2)[1,"Tier.1"], 2))
S2.DICs = c("-",
round(GetMets(Model.S2)[1,"Tier.1"], 2),
"-",
round(GetMets(Model1.st)[1,"Tier.2"], 2),
round(GetMets(Model1)[1,"Tier.2"], 2),
round(GetMets(Model1.ou)[1,"Tier.2"], 2),
round(GetMets(Model2)[1,"Tier.2"], 2))
S3.DICs = c("-",
"-",
round(GetMets(Model.S3)[1,"Tier.1"], 2),
round(GetMets(Model1.st)[1,"Tier.3"], 2),
round(GetMets(Model1)[1,"Tier.3"], 2),
round(GetMets(Model1.ou)[1,"Tier.3"], 2),
round(GetMets(Model2)[1,"Tier.3"], 2))
#Watanabe-Akaike information criteria
S1.WAICs = c(round(GetMets(Model.S1)[2,"Tier.1"], 2),
"-",
"-",
round(GetMets(Model1.st)[2,"Tier.1"], 2),
round(GetMets(Model1)[2,"Tier.1"], 2),
round(GetMets(Model1.ou)[2,"Tier.1"], 2),
round(GetMets(Model2)[2,"Tier.1"], 2))
S2.WAICs = c("-",
round(GetMets(Model.S2)[2,"Tier.1"], 2),
"-",
round(GetMets(Model1.st)[2,"Tier.2"], 2),
round(GetMets(Model1)[2,"Tier.2"], 2),
round(GetMets(Model1.ou)[2,"Tier.2"], 2),
round(GetMets(Model2)[2,"Tier.2"], 2))
S3.WAICs = c("-",
"-",
round(GetMets(Model.S3)[2,"Tier.1"], 2),
round(GetMets(Model1.st)[2,"Tier.3"], 2),
round(GetMets(Model1)[2,"Tier.3"], 2),
round(GetMets(Model1.ou)[2,"Tier.3"], 2),
round(GetMets(Model2)[2,"Tier.3"], 2))
Model_mets2 = as.data.frame(cbind(Model = Models,
DIC1 = S1.DICs,
WAIC1 = S1.WAICs,
DIC2 = S2.DICs,
WAIC2 = S2.WAICs,
DIC3 = S3.DICs,
WAIC3 = S3.WAICs,
Type = Type))kable(Model_mets2) %>%
kable_styling("striped", full_width = F) %>%
row_spec(0, font_size = 20) | Model | DIC1 | WAIC1 | DIC2 | WAIC2 | DIC3 | WAIC3 | Type |
|---|---|---|---|---|---|---|---|
| Model1 | 6272.12 | 6263.31 |
|
|
|
|
Base Model (RM-LF) |
| Model2 |
|
|
11600.4 | 11036.49 |
|
|
Base Model (LM-MF) |
| Model3 |
|
|
|
|
13369.35 | 12915.51 | Base Model (MM-HF) |
| Model4 | 5106.56 | 4857.62 | 10077.37 | 9723.6 | 21235.58 | 27294.45 | Joint Random Effects |
| Model5 | 5077.18 | 4766.05 | 10019.93 | 9687.76 | 18299.66 | 19876.19 | Joint (Full) |
| Model6 | 4881.74 | 4827.58 | 10318.47 | 10253.8 | 21952.87 | 27375.4 | Joint (Ornstein-Uhlenbeck) |
| Model7 | 9073.67 | 8904.25 | 17470.98 | 17447.3 | 17625.67 | 17630.07 | Joint (Non Spatiotemporal) |
4.2 Weekly Activity Proportion
idat = inla.stack.index(Joint.stk, "s1.0")$data
pnts.df$S1_pr = Model1$summary.fitted.values$mean[idat]
idat = inla.stack.index(Joint.stk, "s2.0")$data
pnts.df$S2_pr = Model1$summary.fitted.values$mean[idat]
idat = inla.stack.index(Joint.stk, "s3.0")$data
pnts.df$S3_pr = Model1$summary.fitted.values$mean[idat]
#Quick correlation check
cor(pnts.df[,c("S1_pr","S2_pr","S3_pr")])## S1_pr S2_pr S3_pr
## S1_pr 1.0000000 -0.6641648 -0.3403381
## S2_pr -0.6641648 1.0000000 -0.4208792
## S3_pr -0.3403381 -0.4208792 1.0000000Weekly Movement State Probability
pnts.df$NState = ifelse(pnts.df$S1_pr > pnts.df$S2_pr & pnts.df$S1_pr > pnts.df$S3_pr, "State1",
ifelse(pnts.df$S2_pr > pnts.df$S1_pr & pnts.df$S2_pr > pnts.df$S3_pr, "State2",
ifelse(pnts.df$S3_pr > pnts.df$S1_pr & pnts.df$S3_pr > pnts.df$S2_pr, "State3", NA)))
which(is.na(pnts.df$NState))## integer(0)unique(pnts.df$NState)## [1] "State2" "State3" "State1"pnts.df$Week = week(pnts.df$Date)
Plot.df = pnts.df %>%
select(Week, NState) %>%
group_by(Week, NState) %>%
summarise(Count = length(NState))
Plot.df1 = Plot.df %>%
group_by(Week) %>%
summarise(Sum = sum(Count))
Plot.df2 = Plot.df %>%
group_by(NState) %>%
summarise(StateT = sum(Count))
Plot.df$TWeek = with(Plot.df1,
Sum[match(
Plot.df$Week,
Week)])
Plot.df$TState = with(Plot.df2,
StateT[match(
Plot.df$NState,
NState)])
Plot.df$StateProp = Plot.df$Count/Plot.df$TState
Plot.df3 = Plot.df %>%
group_by(Week) %>%
summarise(PropT = sum(StateProp))
Plot.df$WeekProp = with(Plot.df3,
PropT[match(
Plot.df$Week,
Week)])
Plot.df$IndState = Plot.df$StateProp/Plot.df$WeekProp
colScale = scale_fill_manual(name = "Activity",
values = c("lightgray", "gray65", "gray50"),
labels = c("RM-LF", "LM-MF", "MM-HF"))
ggplot() +
geom_bar(data = Plot.df,
aes(x = Plot.df$Week, y = Plot.df$IndState,
fill = Plot.df$NState), stat="identity",
width = 1,
position="stack",
col ="transparent") +
colScale +
theme_classic() +
xlab("Week of Year") +
ylab("Activity Proportion") +
scale_x_continuous(breaks = seq(0, 54, 13),
labels = c("1", seq(13, 52, 13)),
limits =c(0,52)) +
theme(plot.title = element_text(hjust = 0.5),
legend.position = "bottom",
legend.title=element_text(size=22),
legend.text=element_text(size=18),
strip.background = element_blank(),
title = element_text(face="bold", size=22, hjust=0.5),
strip.text = element_text(face="bold", size = 22),
axis.title.y = element_text(face="bold", size = 24),
axis.text.x = element_text(face="bold", size=16, vjust=0.5,
hjust=0.5, angle=0),
axis.text.y = element_text(face="bold", size=20),
axis.title.x = element_text(face="bold", size = 24)) 
4.3 Cohort Effect
Comb.plt.df1 = as.data.frame(Model1$summary.random$cohort1)[,1:6] #Get effect estimates (Full model)
names(Comb.plt.df1) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975") #Clean labels
Comb.plt.df1$State = "RM-LF"
Comb.plt.df2 = as.data.frame(Model1$summary.random$cohort2)[,1:6] #Get effect estimates (Full model)
names(Comb.plt.df2) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975") #Clean labels
Comb.plt.df2$State = "LM-MF"
Comb.plt.df3 = as.data.frame(Model1$summary.random$cohort3)[,1:6] #Get effect estimates (Full model)
names(Comb.plt.df3) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975") #Clean labels
Comb.plt.df3$State = "MM-HF"
Comb.plt.df = rbind(Comb.plt.df1, Comb.plt.df2, Comb.plt.df3)
Plt.all = ggplot(Comb.plt.df , aes(ID, Mean)) +
geom_hline(yintercept = 0,
linetype = "solid",
col = "lightgray",
size = 0.5) +
geom_point(size=2, pch=19, col = "black") +
geom_linerange(aes(ymin=Q025, ymax=Q975), colour="black") +
geom_hline(yintercept = 0,
linetype = "dotted",
colour = "darkgray",
size = 0.5) +
ylim(-4,4.2) +
facet_wrap(~State, ncol=1) +
xlab("Cohort Group") +
ylab("State Probability (logit)") +
theme_classic() +
theme(plot.margin = unit(c(1,1,1,1),"cm"),
axis.text=element_text(size=18),
legend.title = element_text(size=18, face="bold"),
legend.key = element_rect(fill = "gray80", linetype=0),
strip.background = element_blank(),
legend.background = element_rect(fill = "gray80", linetype=0),
strip.text = element_text(face="bold", size = 20),
axis.title.y = element_text(face="bold", size = 22),
axis.title.x = element_text(face="bold", size = 22, vjust=-2))
Plt.all
4.4 State Transition Probabilities
Using the markovchain package to estimate transition probabilities based on model fitted values. Looping through by month to look for temporal patterns.
library(markovchain) ## Package: markovchain
## Version: 0.7.0
## Date: 2019-08-13
## BugReport: http://github.com/spedygiorgio/markovchain/issues#Overall
Transitions.overall = markovchainFit(data=pnts.df$NState)
Transitions.overall$estimate## MLE Fit
## A 3 - dimensional discrete Markov Chain defined by the following states:
## State1, State2, State3
## The transition matrix (by rows) is defined as follows:
## State1 State2 State3
## State1 0.925909784 0.00544781 0.06864241
## State2 0.002422328 0.91363876 0.08393892
## State3 0.067731629 0.16911608 0.76315229for(i in 1:12){
tmp.wk = pnts.df %>% filter(month == i)
tmp.trans = markovchainFit(data=tmp.wk$NState)
S1xS1 = tmp.trans$estimate["State1"][1]
S1xS2 = tmp.trans$estimate["State1"][2]
S1xS3 = tmp.trans$estimate["State1"][3]
S2xS1 = tmp.trans$estimate["State2"][1]
S2xS2 = tmp.trans$estimate["State2"][2]
S2xS3 = tmp.trans$estimate["State2"][3]
S3xS1 = tmp.trans$estimate["State3"][1]
S3xS2 = tmp.trans$estimate["State3"][2]
S3xS3 = tmp.trans$estimate["State3"][3]
tmp.trans.df = as.data.frame(
cbind(S1xS2, S1xS3,
S2xS3, S2xS1,
S3xS1, S3xS2))
tmp.trans.df$Month = i
if(i == 1){State1.loop = tmp.trans.df
} else{State1.loop = rbind(State1.loop, tmp.trans.df)}
}
State1.loop.mlt = melt(State1.loop, "Month")
State1.loop.mlt$variable = gsub("x", ":", State1.loop.mlt$variable)
State1.loop.mlt$Label = month.abb[State1.loop.mlt$Month]
State1.loop.mlt$Label = ordered(factor(State1.loop.mlt$Label), levels = month.abb)
State1.loop.mlt$variable2 = State1.loop.mlt$variable
State1.loop.mlt$variable2[State1.loop.mlt$variable2 == "S1:S2"] = "RM-LF:LM-MF"
State1.loop.mlt$variable2[State1.loop.mlt$variable2 == "S1:S3"] = "RM-LF:MM-HF"
State1.loop.mlt$variable2[State1.loop.mlt$variable2 == "S2:S3"] = "LM-MF:MM-HF"
State1.loop.mlt$variable2[State1.loop.mlt$variable2 == "S2:S1"] = "LM-MF:RM-LF"
State1.loop.mlt$variable2[State1.loop.mlt$variable2 == "S3:S1"] = "MM-HF:RM-LF"
State1.loop.mlt$variable2[State1.loop.mlt$variable2 == "S3:S2"] = "MM-HF:LM-MF"
ggplot(State1.loop.mlt, aes(factor(Label), value)) +
geom_bar(stat = "identity",
col = "darkgray",
fill = "lightgray",
linetype= "solid",
lwd = 0.2) +
ylab("State Transition Probability") +
xlab("Month of Year") +
facet_wrap(~variable2, ncol=2, scales ="free") +
theme_classic() +
#scale_y_continuous(breaks = seq(0, 0.5, 0.1),
# limits = c(0, 0.5)) +
#scale_x_discrete(breaks = seq(1, 12, 1),
# labels = month.abb,
# limits =c(1, 12)) +
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),
strip.background = element_blank(),
axis.title.y = element_text(face="bold", size = 23),
axis.text.x = element_text(face="bold", size=16, vjust=0.5,
hjust=0.95, angle=90),
axis.title.x = element_text(face="bold", size = 23))
4.5 Random Effects Table
FE.table = Model1$summary.hyperpar[, c(1:3, 5)]
names(FE.table) = c("Mean", "sd", "Q025", "Q975")
right = function (string, char){
substr(string,nchar(string)-(char-1),nchar(string))
}
FE.table$Covariate = substr(rownames(FE.table), 1, nchar(rownames(FE.table))-1)
FE.table$State = paste("State", right(rownames(FE.table), 1), sep=" ")
FE.table = FE.table %>%
select(Covariate, Mean, sd, Q025, Q975)
#SelMod.fx = xtable(FE.table )
#print(SelMod.fx, include.rownames = F)
kable(FE.table) %>%
kable_styling("striped", full_width = F) %>%
row_spec(0, font_size = 20) | Covariate | Mean | sd | Q025 | Q975 | |
|---|---|---|---|---|---|
| Range for field1 | Range for field | 3.2587474 | 0.4223770 | 2.5189431 | 4.1769130 |
| Stdev for field1 | Stdev for field | 6.0211619 | 0.5481103 | 5.0200143 | 7.1725607 |
| Range for field2 | Range for field | 3.4436076 | 0.5700988 | 2.5131559 | 4.7442700 |
| Stdev for field2 | Stdev for field | 3.5055135 | 0.3022570 | 2.9250788 | 4.1085095 |
| Range for field3 | Range for field | 19.6327805 | 7.2547087 | 8.8792265 | 37.0712287 |
| Stdev for field3 | Stdev for field | 2.2904984 | 0.5717713 | 1.4357924 | 3.6650889 |
| Precision for NN | Precision for N | 0.0189401 | 0.0032128 | 0.0132516 | 0.0258353 |
| Precision for pig | Precision for pi | 1.0723591 | 0.2880699 | 0.5888864 | 1.7241024 |
| Precision for year | Precision for yea | 1.1206685 | 0.3456492 | 0.5650801 | 1.9217497 |
| Precision for cohort1 | Precision for cohort | 0.6180443 | 0.1991991 | 0.3067747 | 1.0824203 |
| Precision for cohort2 | Precision for cohort | 0.7504572 | 0.2444314 | 0.3649799 | 1.3185499 |
| Precision for cohort3 | Precision for cohort | 1.1401212 | 0.3539134 | 0.5763347 | 1.9666626 |
| Precision for day1 | Precision for day | 1.0813514 | 0.3365788 | 0.5696912 | 1.8774277 |
| Precision for hour1 | Precision for hour | 1.7791871 | 0.5208115 | 0.9766455 | 3.0095204 |
| Precision for day2 | Precision for day | 3.2157452 | 0.6858512 | 2.0273846 | 4.7200324 |
| Precision for hour2 | Precision for hour | 2.1088858 | 0.5720289 | 1.1721126 | 3.4079348 |
| Precision for day3 | Precision for day | 2.3529206 | 0.4496542 | 1.5342066 | 3.2981755 |
| Precision for hour3 | Precision for hour | 0.8481031 | 0.1674493 | 0.5489370 | 1.2049255 |
| Precision for TStep1 | Precision for TStep | 0.0572201 | 0.0070622 | 0.0438889 | 0.0716756 |
| Rho for TStep1 | Rho for TStep | 0.9872225 | 0.0019324 | 0.9832598 | 0.9908314 |
| Precision for TStep2 | Precision for TStep | 0.4397445 | 0.0365464 | 0.3736638 | 0.5177768 |
| Rho for TStep2 | Rho for TStep | 0.9366355 | 0.0080606 | 0.9191139 | 0.9507792 |
| Precision for TStep3 | Precision for TStep | 0.0825132 | 0.0142331 | 0.0562984 | 0.1126557 |
| Rho for TStep3 | Rho for TStep | 0.6641886 | 0.0141731 | 0.6336971 | 0.6891912 |
| Beta for field1x | Beta for field1 | 0.1769045 | 0.0666443 | 0.0495452 | 0.3117547 |
| Beta for field1xx | Beta for field1x | -1.2394266 | 0.0984295 | -1.4465569 | -1.0597126 |
| Beta for field2x | Beta for field2 | -1.3284747 | 0.1014085 | -1.5412968 | -1.1402709 |
4.6 Hourly Movement States
mic.d1 = as.data.frame(Model1$summary.random$hour1)
names(mic.d1) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975", "mode", "kld")
mic.d1$State = "RM-LF"
mic.d2 = as.data.frame(Model1$summary.random$hour2)
names(mic.d2) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975", "mode", "kld")
mic.d2$State = "LM-MF"
mic.d3 = as.data.frame(Model1$summary.random$hour3)
names(mic.d3) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975", "mode", "kld")
mic.d3$State = "MM-HF"
mic.d = rbind(mic.d1, mic.d2, mic.d3)
mic.d$State = ordered(factor(mic.d$State), levels = c("RM-LF", "LM-MF", "MM-HF"))
Raw.plt2 = ggplot(mic.d, aes(ID, Mean)) +
geom_hline(yintercept = 0,
linetype = "solid",
col = "darkgray",
size = 0.5) +
geom_smooth(col = "black",
linetype= "solid",
method = "loess",
span = 0.3,
se = FALSE,
lwd = 1) +
geom_smooth(data = mic.d, aes(ID, Q025),
col = "grey",
method = "loess",
span = 0.3,
size = 0.5,
se = FALSE,
linetype= "dashed") +
geom_smooth(data = mic.d, aes(ID, Q975),
col = "grey",
method = "loess",
span = 0.3,
size = 0.5,
se = FALSE,
linetype= "dashed") +
ylab("Movement State Probability (logit)") +
xlab("Hour of Day") +
#ggtitle("Aerial Counts (COMPRECNTDEN)") +
theme_classic() +
scale_y_continuous(breaks = seq(-1.5, 2, 0.5)) +
#limits = c(-1, 1)) +
scale_x_continuous(breaks = seq(0, 24, 6),
labels = c(seq(0, 18, 6), "23"),
limits =c(0,23)) +
facet_grid(rows = vars(State)) +
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 = 22),
strip.background=element_rect(size=0.5),
axis.title.y = element_text(face="bold", size = 22),
axis.text.x = element_text(face="bold", size=18, vjust=0.5,
hjust=0.5, angle=0),
axis.title.x = element_text(face="bold", size = 22))
Raw.plt2
4.7 Fixed Effects Table
FE.table = Model1$summary.fixed[, c(1:3, 5)]
names(FE.table) = c("Mean", "sd", "Q025", "Q975")
right = function (string, char){
substr(string,nchar(string)-(char-1),nchar(string))
}
FE.table$Covariate = substr(rownames(FE.table), 1, nchar(rownames(FE.table))-1)
FE.table$State = paste("State", right(rownames(FE.table), 1), sep=" ")
FE.table = FE.table %>%
select(Covariate, State, Mean, sd, Q025, Q975)
for(i in 3:6){
FE.table[,i] = round(FE.table[,i], 2)
}
#SelMod.fx = xtable(FE.table )
#print(SelMod.fx, include.rownames = F)
kable(FE.table) %>%
kable_styling("striped", full_width = F) %>%
row_spec(0, font_size = 20) | Covariate | State | Mean | sd | Q025 | Q975 | |
|---|---|---|---|---|---|---|
| intercept1 | intercept | State 1 | -4.47 | 2.42 | -9.22 | 0.28 |
| intercept2 | intercept | State 2 | -6.85 | 2.99 | -12.72 | -0.99 |
| intercept3 | intercept | State 3 | 0.64 | 2.37 | -4.02 | 5.30 |
| Dvlpd_A1 | Dvlpd_A | State 1 | 0.01 | 0.07 | -0.12 | 0.15 |
| Dvlpd_A2 | Dvlpd_A | State 2 | 0.00 | 0.05 | -0.09 | 0.09 |
| Dvlpd_A3 | Dvlpd_A | State 3 | 0.02 | 0.04 | -0.05 | 0.09 |
| Rprn_Ar1 | Rprn_Ar | State 1 | 0.35 | 0.14 | 0.07 | 0.63 |
| Rprn_Ar2 | Rprn_Ar | State 2 | 0.24 | 0.10 | 0.04 | 0.44 |
| Rprn_Ar3 | Rprn_Ar | State 3 | -0.45 | 0.08 | -0.61 | -0.30 |
| Open_Ar1 | Open_Ar | State 1 | 0.46 | 0.10 | 0.27 | 0.65 |
| Open_Ar2 | Open_Ar | State 2 | -0.04 | 0.09 | -0.22 | 0.15 |
| Open_Ar3 | Open_Ar | State 3 | -0.35 | 0.06 | -0.47 | -0.23 |
| Hrdms_A1 | Hrdms_A | State 1 | -0.17 | 0.08 | -0.32 | -0.02 |
| Hrdms_A2 | Hrdms_A | State 2 | -0.06 | 0.05 | -0.15 | 0.03 |
| Hrdms_A3 | Hrdms_A | State 3 | 0.08 | 0.04 | 0.01 | 0.16 |
| PCI_21 | PCI_2 | State 1 | 0.34 | 0.19 | -0.04 | 0.72 |
| PCI_22 | PCI_2 | State 2 | 0.41 | 0.10 | 0.22 | 0.60 |
| PCI_23 | PCI_2 | State 3 | -0.41 | 0.08 | -0.58 | -0.25 |
| PCI_51 | PCI_5 | State 1 | 0.20 | 0.10 | 0.01 | 0.40 |
| PCI_52 | PCI_5 | State 2 | -0.02 | 0.08 | -0.18 | 0.13 |
| PCI_53 | PCI_5 | State 3 | -0.11 | 0.06 | -0.22 | 0.00 |
| Edge1_31 | Edge1_3 | State 1 | -0.13 | 0.08 | -0.28 | 0.02 |
| Edge1_32 | Edge1_3 | State 2 | 0.04 | 0.05 | -0.07 | 0.14 |
| Edge1_33 | Edge1_3 | State 3 | 0.06 | 0.04 | -0.02 | 0.14 |
| HOURLYDRYBULBTEMPC1 | HOURLYDRYBULBTEMPC | State 1 | 0.85 | 0.31 | 0.23 | 1.46 |
| HOURLYDRYBULBTEMPC2 | HOURLYDRYBULBTEMPC | State 2 | 0.08 | 0.14 | -0.19 | 0.35 |
| HOURLYDRYBULBTEMPC3 | HOURLYDRYBULBTEMPC | State 3 | -0.15 | 0.10 | -0.34 | 0.04 |
| HOURLYStationPressure1 | HOURLYStationPressure | State 1 | 0.16 | 0.25 | -0.33 | 0.66 |
| HOURLYStationPressure2 | HOURLYStationPressure | State 2 | 0.10 | 0.11 | -0.12 | 0.32 |
| HOURLYStationPressure3 | HOURLYStationPressure | State 3 | -0.10 | 0.07 | -0.23 | 0.04 |
4.8 Landscape Proportions
Proportion Emergent
range(Pnts$Open_Ar)## [1] 0 100myFixed = Model1$summary.fixed[,c(1:3,5)]
names(myFixed) = c("Mean", "SD", "Q0.025", "Q0.975")
PlotFrame = Pnts@data %>% select(state, Open_Ar, s_Open_Ar)
PlotFrame$ExtRipM = plogis(PlotFrame$s_Open_Ar*myFixed["Open_Ar1","Mean"])
PlotFrame$ExtRipL = plogis(PlotFrame$s_Open_Ar*myFixed["Open_Ar1","Q0.025"])
PlotFrame$ExtRipH = plogis(PlotFrame$s_Open_Ar*myFixed["Open_Ar1","Q0.975"])
PlotFrame$State = "RM-LF"
PlotFrame2 = Pnts@data %>% select(state, Open_Ar, s_Open_Ar)
PlotFrame2$ExtRipM = plogis(PlotFrame2$s_Open_Ar*myFixed["Open_Ar2","Mean"])
PlotFrame2$ExtRipL = plogis(PlotFrame2$s_Open_Ar*myFixed["Open_Ar2","Q0.025"])
PlotFrame2$ExtRipH = plogis(PlotFrame2$s_Open_Ar*myFixed["Open_Ar2","Q0.975"])
PlotFrame2$State = "LM-MF"
PlotFrame3 = Pnts@data %>% select(state, Open_Ar, s_Open_Ar)
PlotFrame3$ExtRipM = plogis(PlotFrame3$s_Open_Ar*myFixed["Open_Ar3","Mean"])
PlotFrame3$ExtRipL = plogis(PlotFrame3$s_Open_Ar*myFixed["Open_Ar3","Q0.025"])
PlotFrame3$ExtRipH = plogis(PlotFrame3$s_Open_Ar*myFixed["Open_Ar3","Q0.975"])
PlotFrame3$State = "MM-HF"
CFrn = rbind(PlotFrame, PlotFrame2, PlotFrame3)
CFrn$State = factor(CFrn$State)
ggplot(CFrn, aes(Open_Ar, ExtRipM, group=State)) +
geom_line(aes(col = State,
linetype= State),
lwd = 1) +
scale_linetype_manual(values=c("solid", "longdash", "dotted")) +
scale_colour_manual(values=c("black", "black", "black")) +
geom_ribbon(aes(ymin=ExtRipL,
ymax=ExtRipH),
alpha=0.1, #transparency
linetype=1, #solid, dashed or other line types
colour="transparent", #border line color
size=1, #border line size
fill="gray30") +
geom_hline(yintercept = 0,
linetype = "solid",
col = "darkgray",
size = 0.5) +
ylab("Movement Probability") +
xlab("Proportion Emergent/Herbaceous") +
#ggtitle("Aerial Counts (COMPRECNTDEN)") +
theme_classic() +
scale_y_continuous(breaks = seq(0, 1, 0.2),
limits = c(0, 1)) +
scale_x_continuous(breaks = seq(0, 100, 20),
limits =c(0,100)) +
#facet_grid(rows = vars(State)) +
theme(plot.title = element_text(hjust = 0.5),
title = element_text(face="bold", size=18, hjust=0.5),
legend.text = element_text(size=14, face="bold"),
legend.direction = "vertical",
legend.position=c(0.2, 0.2),
legend.key.size = unit(3,"line"),
legend.title = element_text(size=18, face="bold"),
axis.text=element_text(size=12),
strip.text = element_text(face="bold", size = 20),
axis.title.y = element_text(face="bold", size = 22),
axis.text.x = element_text(face="bold", size=18, vjust=0.5,
hjust=0.5, angle=0),
axis.text.y = element_text(face="bold", size=18, vjust=0.5,
hjust=0.5, angle=0),
axis.title.x = element_text(face="bold", size = 22))
CFrn$Covariate = "c"
CFrn1 = CFrnProportion Riparian
myFixed = Model1$summary.fixed[,c(1:3,5)]
names(myFixed) = c("Mean", "SD", "Q0.025", "Q0.975")
PlotFrame = Pnts@data %>% select(state, Rprn_Ar, s_Rprn_Ar)
PlotFrame$ExtRipM = plogis(PlotFrame$s_Rprn_Ar*myFixed["Rprn_Ar1","Mean"])
PlotFrame$ExtRipL = plogis(PlotFrame$s_Rprn_Ar*myFixed["Rprn_Ar1","Q0.025"])
PlotFrame$ExtRipH = plogis(PlotFrame$s_Rprn_Ar*myFixed["Rprn_Ar1","Q0.975"])
PlotFrame$State = "RM-LF"
PlotFrame2 = Pnts@data %>% select(state, Rprn_Ar, s_Rprn_Ar)
PlotFrame2$ExtRipM = plogis(PlotFrame2$s_Rprn_Ar*myFixed["Rprn_Ar2","Mean"])
PlotFrame2$ExtRipL = plogis(PlotFrame2$s_Rprn_Ar*myFixed["Rprn_Ar2","Q0.025"])
PlotFrame2$ExtRipH = plogis(PlotFrame2$s_Rprn_Ar*myFixed["Rprn_Ar2","Q0.975"])
PlotFrame2$State = "LM-MF"
PlotFrame3 = Pnts@data %>% select(state, Rprn_Ar, s_Rprn_Ar)
PlotFrame3$ExtRipM = plogis(PlotFrame3$s_Rprn_Ar*myFixed["Rprn_Ar3","Mean"])
PlotFrame3$ExtRipL = plogis(PlotFrame3$s_Rprn_Ar*myFixed["Rprn_Ar3","Q0.025"])
PlotFrame3$ExtRipH = plogis(PlotFrame3$s_Rprn_Ar*myFixed["Rprn_Ar3","Q0.975"])
PlotFrame3$State = "MM-HF"
CFrn = rbind(PlotFrame, PlotFrame2, PlotFrame3)
CFrn$State = factor(CFrn$State)
ggplot(CFrn, aes(Rprn_Ar, ExtRipM, group=State)) +
geom_line(aes(col = State,
linetype= State),
lwd = 1) +
scale_linetype_manual(values=c("solid", "longdash", "dotted")) +
scale_colour_manual(values=c("black", "black", "black")) +
geom_ribbon(aes(ymin=ExtRipL,
ymax=ExtRipH),
alpha=0.1, #transparency
linetype=1, #solid, dashed or other line types
colour="transparent", #border line color
size=1, #border line size
fill="gray30") +
geom_hline(yintercept = 0,
linetype = "solid",
col = "darkgray",
size = 0.5) +
ylab("Movement Probability") +
xlab("Proportion Riparian") +
#ggtitle("Aerial Counts (COMPRECNTDEN)") +
theme_classic() +
scale_y_continuous(breaks = seq(0, 1, 0.2),
limits = c(0, 1)) +
scale_x_continuous(breaks = seq(0, 100, 20),
limits =c(0,100)) +
#facet_grid(rows = vars(State)) +
theme(plot.title = element_text(hjust = 0.5),
title = element_text(face="bold", size=18, hjust=0.5),
legend.text = element_text(size=14, face="bold"),
legend.direction = "vertical",
legend.position=c(0.2, 0.2),
legend.key.size = unit(3,"line"),
legend.title = element_text(size=18, face="bold"),
axis.text=element_text(size=12),
strip.text = element_text(face="bold", size = 20),
axis.title.y = element_text(face="bold", size = 22),
axis.text.x = element_text(face="bold", size=18, vjust=0.5,
hjust=0.5, angle=0),
axis.text.y = element_text(face="bold", size=18, vjust=0.5,
hjust=0.5, angle=0),
axis.title.x = element_text(face="bold", size = 22))
CFrn$Covariate = "b"
CFrn2 = CFrnProportion Developed
myFixed = Model1$summary.fixed[,c(1:3,5)]
names(myFixed) = c("Mean", "SD", "Q0.025", "Q0.975")
PlotFrame = Pnts@data %>% select(state, Dvlpd_A, s_Dvlpd_A)
PlotFrame$ExtRipM = plogis(PlotFrame$s_Dvlpd_A*myFixed["Dvlpd_A1","Mean"])
PlotFrame$ExtRipL = plogis(PlotFrame$s_Dvlpd_A*myFixed["Dvlpd_A1","Q0.025"])
PlotFrame$ExtRipH = plogis(PlotFrame$s_Dvlpd_A*myFixed["Dvlpd_A1","Q0.975"])
PlotFrame$State = "RM-LF"
PlotFrame2 = Pnts@data %>% select(state, Dvlpd_A, s_Dvlpd_A)
PlotFrame2$ExtRipM = plogis(PlotFrame2$s_Dvlpd_A*myFixed["Dvlpd_A2","Mean"])
PlotFrame2$ExtRipL = plogis(PlotFrame2$s_Dvlpd_A*myFixed["Dvlpd_A2","Q0.025"])
PlotFrame2$ExtRipH = plogis(PlotFrame2$s_Dvlpd_A*myFixed["Dvlpd_A2","Q0.975"])
PlotFrame2$State = "LM-MF"
PlotFrame3 = Pnts@data %>% select(state, Dvlpd_A, s_Dvlpd_A)
PlotFrame3$ExtRipM = plogis(PlotFrame3$s_Dvlpd_A*myFixed["Dvlpd_A3","Mean"])
PlotFrame3$ExtRipL = plogis(PlotFrame3$s_Dvlpd_A*myFixed["Dvlpd_A3","Q0.025"])
PlotFrame3$ExtRipH = plogis(PlotFrame3$s_Dvlpd_A*myFixed["Dvlpd_A3","Q0.975"])
PlotFrame3$State = "MM-HF"
CFrn = rbind(PlotFrame, PlotFrame2, PlotFrame3)
CFrn$State = factor(CFrn$State)
ggplot(CFrn, aes(Dvlpd_A, ExtRipM, group=State)) +
geom_line(aes(col = State,
linetype= State),
lwd = 1) +
scale_linetype_manual(values=c("solid", "longdash", "dotted")) +
scale_colour_manual(values=c("black", "black", "black")) +
geom_ribbon(aes(ymin=ExtRipL,
ymax=ExtRipH),
alpha=0.1, #transparency
linetype=1, #solid, dashed or other line types
colour="transparent", #border line color
size=1, #border line size
fill="gray30") +
geom_hline(yintercept = 0,
linetype = "solid",
col = "darkgray",
size = 0.5) +
ylab("Movement Probability") +
xlab("Proportion Developed") +
#ggtitle("Aerial Counts (COMPRECNTDEN)") +
theme_classic() +
scale_y_continuous(breaks = seq(0, 1, 0.2),
limits = c(0, 1)) +
scale_x_continuous(breaks = seq(0, 100, 20),
limits =c(0,100)) +
#facet_grid(rows = vars(State)) +
theme(plot.title = element_text(hjust = 0.5),
title = element_text(face="bold", size=18, hjust=0.5),
legend.text = element_text(size=14, face="bold"),
legend.direction = "vertical",
legend.position=c(0.2, 0.2),
legend.key.size = unit(3,"line"),
legend.title = element_text(size=18, face="bold"),
axis.text=element_text(size=12),
strip.text = element_text(face="bold", size = 20),
axis.title.y = element_text(face="bold", size = 22),
axis.text.x = element_text(face="bold", size=18, vjust=0.5,
hjust=0.5, angle=0),
axis.text.y = element_text(face="bold", size=18, vjust=0.5,
hjust=0.5, angle=0),
axis.title.x = element_text(face="bold", size = 22))
CFrn$Covariate = "a"
CFrn3 = CFrnProportion Hard Mast
myFixed = Model1$summary.fixed[,c(1:3,5)]
names(myFixed) = c("Mean", "SD", "Q0.025", "Q0.975")
PlotFrame = Pnts@data %>% select(state, Hrdms_A, s_Hrdms_A)
PlotFrame$ExtRipM = plogis(PlotFrame$s_Hrdms_A*myFixed["Hrdms_A1","Mean"])
PlotFrame$ExtRipL = plogis(PlotFrame$s_Hrdms_A*myFixed["Hrdms_A1","Q0.025"])
PlotFrame$ExtRipH = plogis(PlotFrame$s_Hrdms_A*myFixed["Hrdms_A1","Q0.975"])
PlotFrame$State = "RM-LF"
PlotFrame2 = Pnts@data %>% select(state, Hrdms_A, s_Hrdms_A)
PlotFrame2$ExtRipM = plogis(PlotFrame2$s_Hrdms_A*myFixed["Hrdms_A2","Mean"])
PlotFrame2$ExtRipL = plogis(PlotFrame2$s_Hrdms_A*myFixed["Hrdms_A2","Q0.025"])
PlotFrame2$ExtRipH = plogis(PlotFrame2$s_Hrdms_A*myFixed["Hrdms_A2","Q0.975"])
PlotFrame2$State = "LM-MF"
PlotFrame3 = Pnts@data %>% select(state, Hrdms_A, s_Hrdms_A)
PlotFrame3$ExtRipM = plogis(PlotFrame3$s_Hrdms_A*myFixed["Hrdms_A3","Mean"])
PlotFrame3$ExtRipL = plogis(PlotFrame3$s_Hrdms_A*myFixed["Hrdms_A3","Q0.025"])
PlotFrame3$ExtRipH = plogis(PlotFrame3$s_Hrdms_A*myFixed["Hrdms_A3","Q0.975"])
PlotFrame3$State = "MM-HF"
CFrn = rbind(PlotFrame, PlotFrame2, PlotFrame3)
CFrn$State = factor(CFrn$State)
ggplot(CFrn, aes(Hrdms_A, ExtRipM, group=State)) +
geom_line(aes(col = State,
linetype= State),
lwd = 1) +
scale_linetype_manual(values=c("solid", "longdash", "dotted")) +
scale_colour_manual(values=c("black", "black", "black")) +
geom_ribbon(aes(ymin=ExtRipL,
ymax=ExtRipH),
alpha=0.1, #transparency
linetype=1, #solid, dashed or other line types
colour="transparent", #border line color
size=1, #border line size
fill="gray30") +
geom_hline(yintercept = 0,
linetype = "solid",
col = "darkgray",
size = 0.5) +
ylab("Movement Probability") +
xlab("Proportion Hard Mast") +
#ggtitle("Aerial Counts (COMPRECNTDEN)") +
theme_classic() +
scale_y_continuous(breaks = seq(0, 1, 0.2),
limits = c(0, 1)) +
scale_x_continuous(breaks = seq(0, 100, 20),
limits =c(0,100)) +
#facet_grid(rows = vars(State)) +
theme(plot.title = element_text(hjust = 0.5),
title = element_text(face="bold", size=18, hjust=0.5),
legend.text = element_text(size=14, face="bold"),
legend.direction = "vertical",
legend.position=c(0.15, 0.8),
legend.key.size = unit(3,"line"),
legend.title = element_text(size=18, face="bold"),
axis.text=element_text(size=12),
strip.text = element_text(face="bold", size = 20),
axis.title.y = element_text(face="bold", size = 22),
axis.text.x = element_text(face="bold", size=18, vjust=0.5,
hjust=0.5, angle=0),
axis.text.y = element_text(face="bold", size=18, vjust=0.5,
hjust=0.5, angle=0),
axis.title.x = element_text(face="bold", size = 22))
CFrn$Covariate = "d"
CFrn4 = CFrnProportions Combined
names(CFrn1)[2:3] = c("Value", "s_Value")
names(CFrn2)[2:3] = c("Value", "s_Value")
names(CFrn3)[2:3] = c("Value", "s_Value")
names(CFrn4)[2:3] = c("Value", "s_Value")
Comb.set = rbind(CFrn1, CFrn2, CFrn3, CFrn4)
ggplot(Comb.set , aes(Value, ExtRipM, group=State)) +
geom_line(aes(col = State,
linetype= State),
lwd = 1) +
scale_linetype_manual(values=c("solid", "longdash", "dotted")) +
scale_colour_manual(values=c("black", "black", "black")) +
geom_ribbon(aes(ymin=ExtRipL,
ymax=ExtRipH),
alpha=0.1, #transparency
linetype=1, #solid, dashed or other line types
colour="transparent", #border line color
size=1, #border line size
fill="gray30") +
geom_hline(yintercept = 0,
linetype = "solid",
col = "darkgray",
size = 0.5) +
ylab("Movement Probability") +
xlab("Proportion") +
#ggtitle("Aerial Counts (COMPRECNTDEN)") +
theme_classic() +
scale_y_continuous(breaks = seq(0, 1, 0.2),
limits = c(0, 1)) +
scale_x_continuous(breaks = seq(0, 100, 20),
limits =c(0,100)) +
facet_wrap(~Covariate, ncol=2) +
theme(plot.title = element_text(hjust = 0.5),
title = element_text(face="bold", size=18, hjust=0.5),
legend.text = element_text(size=16, face="bold"),
legend.direction = "horizontal",
#legend.position=c(0.15, 0.8),
legend.position="bottom",
legend.key.size = unit(3,"line"),
legend.title = element_text(size=22, face="bold"),
axis.text=element_text(size=12),
strip.text = element_text(face="bold", size = 20),
strip.background = element_blank(),
axis.title.y = element_text(face="bold", size = 22),
axis.text.x = element_text(face="bold", size=18, vjust=0.5,
hjust=0.5, angle=0),
axis.text.y = element_text(face="bold", size=18, vjust=0.5,
hjust=0.5, angle=0),
axis.title.x = element_text(face="bold", size = 22))
4.9 PCs and Climate
Agriculture PCI
range(Pnts$PCI_2)## [1] 0.000000 4.444444myFixed = Model1$summary.fixed[,c(1:3,5)]
names(myFixed) = c("Mean", "SD", "Q0.025", "Q0.975")
PlotFrame = Pnts@data %>% select(state, PCI_2, s_PCI_2)
PlotFrame$ExtRipM = plogis(PlotFrame$s_PCI_2*myFixed["PCI_21","Mean"])
PlotFrame$ExtRipL = plogis(PlotFrame$s_PCI_2*myFixed["PCI_21","Q0.025"])
PlotFrame$ExtRipH = plogis(PlotFrame$s_PCI_2*myFixed["PCI_21","Q0.975"])
PlotFrame$State = "RM-LF"
PlotFrame2 = Pnts@data %>% select(state, PCI_2, s_PCI_2)
PlotFrame2$ExtRipM = plogis(PlotFrame2$s_PCI_2*myFixed["PCI_22","Mean"])
PlotFrame2$ExtRipL = plogis(PlotFrame2$s_PCI_2*myFixed["PCI_22","Q0.025"])
PlotFrame2$ExtRipH = plogis(PlotFrame2$s_PCI_2*myFixed["PCI_22","Q0.975"])
PlotFrame2$State = "LM-MF"
PlotFrame3 = Pnts@data %>% select(state, PCI_2, s_PCI_2)
PlotFrame3$ExtRipM = plogis(PlotFrame3$s_PCI_2*myFixed["PCI_23","Mean"])
PlotFrame3$ExtRipL = plogis(PlotFrame3$s_PCI_2*myFixed["PCI_23","Q0.025"])
PlotFrame3$ExtRipH = plogis(PlotFrame3$s_PCI_2*myFixed["PCI_23","Q0.975"])
PlotFrame3$State = "MM-HF"
CFrn = rbind(PlotFrame, PlotFrame2, PlotFrame3)
CFrn$State = factor(CFrn$State)
ggplot(CFrn, aes(PCI_2, ExtRipM, group=State)) +
geom_line(aes(col = State,
linetype= State),
lwd = 1) +
scale_linetype_manual(values=c("solid", "longdash", "dotted")) +
scale_colour_manual(values=c("black", "black", "black")) +
geom_ribbon(aes(ymin=ExtRipL,
ymax=ExtRipH),
alpha=0.1, #transparency
linetype=1, #solid, dashed or other line types
colour="transparent", #border line color
size=1, #border line size
fill="gray30") +
geom_hline(yintercept = 0,
linetype = "solid",
col = "darkgray",
size = 0.5) +
ylab("Movement Probability") +
xlab("Agriculture (PCI)") +
#ggtitle("Aerial Counts (COMPRECNTDEN)") +
theme_classic() +
scale_y_continuous(breaks = seq(0, 1, 0.2),
limits = c(0, 1)) +
scale_x_continuous(breaks = seq(0, 5, 1),
limits =c(0,5)) +
#facet_grid(rows = vars(State)) +
theme(plot.title = element_text(hjust = 0.5),
title = element_text(face="bold", size=18, hjust=0.5),
legend.text = element_text(size=14, face="bold"),
legend.direction = "vertical",
legend.position=c(0.15, 0.25),
legend.key.size = unit(3,"line"),
legend.title = element_text(size=18, face="bold"),
axis.text=element_text(size=12),
strip.text = element_text(face="bold", size = 20),
axis.title.y = element_text(face="bold", size = 22),
axis.text.x = element_text(face="bold", size=18, vjust=0.5,
hjust=0.5, angle=0),
axis.text.y = element_text(face="bold", size=18, vjust=0.5,
hjust=0.5, angle=0),
axis.title.x = element_text(face="bold", size = 22))
CFrn$Covariate = "a"
CFrn1 = CFrnShrub PCI
range(Pnts$PCI_5)## [1] 0.000000 4.309644myFixed = Model1$summary.fixed[,c(1:3,5)]
names(myFixed) = c("Mean", "SD", "Q0.025", "Q0.975")
PlotFrame = Pnts@data %>% select(state, PCI_5, s_PCI_5)
PlotFrame$ExtRipM = plogis(PlotFrame$s_PCI_5*myFixed["PCI_51","Mean"])
PlotFrame$ExtRipL = plogis(PlotFrame$s_PCI_5*myFixed["PCI_51","Q0.025"])
PlotFrame$ExtRipH = plogis(PlotFrame$s_PCI_5*myFixed["PCI_51","Q0.975"])
PlotFrame$State = "RM-LF"
PlotFrame2 = Pnts@data %>% select(state, PCI_5, s_PCI_5)
PlotFrame2$ExtRipM = plogis(PlotFrame2$s_PCI_5*myFixed["PCI_52","Mean"])
PlotFrame2$ExtRipL = plogis(PlotFrame2$s_PCI_5*myFixed["PCI_52","Q0.025"])
PlotFrame2$ExtRipH = plogis(PlotFrame2$s_PCI_5*myFixed["PCI_52","Q0.975"])
PlotFrame2$State = "LM-MF"
PlotFrame3 = Pnts@data %>% select(state, PCI_5, s_PCI_5)
PlotFrame3$ExtRipM = plogis(PlotFrame3$s_PCI_5*myFixed["PCI_53","Mean"])
PlotFrame3$ExtRipL = plogis(PlotFrame3$s_PCI_5*myFixed["PCI_53","Q0.025"])
PlotFrame3$ExtRipH = plogis(PlotFrame3$s_PCI_5*myFixed["PCI_53","Q0.975"])
PlotFrame3$State = "MM-HF"
CFrn = rbind(PlotFrame, PlotFrame2, PlotFrame3)
CFrn$State = factor(CFrn$State)
ggplot(CFrn, aes(PCI_5, ExtRipM, group=State)) +
geom_line(aes(col = State,
linetype= State),
lwd = 1) +
scale_linetype_manual(values=c("solid", "longdash", "dotted")) +
scale_colour_manual(values=c("black", "black", "black")) +
geom_ribbon(aes(ymin=ExtRipL,
ymax=ExtRipH),
alpha=0.1, #transparency
linetype=1, #solid, dashed or other line types
colour="transparent", #border line color
size=1, #border line size
fill="gray30") +
geom_hline(yintercept = 0,
linetype = "solid",
col = "darkgray",
size = 0.5) +
ylab("Movement Probability") +
xlab("Shrub (PCI)") +
#ggtitle("Aerial Counts (COMPRECNTDEN)") +
theme_classic() +
scale_y_continuous(breaks = seq(0, 1, 0.2),
limits = c(0, 1)) +
scale_x_continuous(breaks = seq(0, 5, 1),
limits =c(0,5)) +
#facet_grid(rows = vars(State)) +
theme(plot.title = element_text(hjust = 0.5),
title = element_text(face="bold", size=18, hjust=0.5),
legend.text = element_text(size=14, face="bold"),
legend.direction = "vertical",
legend.position=c(0.15, 0.25),
legend.key.size = unit(3,"line"),
legend.title = element_text(size=18, face="bold"),
axis.text=element_text(size=12),
strip.text = element_text(face="bold", size = 20),
axis.title.y = element_text(face="bold", size = 22),
axis.text.x = element_text(face="bold", size=18, vjust=0.5,
hjust=0.5, angle=0),
axis.text.y = element_text(face="bold", size=18, vjust=0.5,
hjust=0.5, angle=0),
axis.title.x = element_text(face="bold", size = 22))
CFrn$Covariate = "b"
CFrn2 = CFrnOpen-Forest Edge
range(Pnts$Edge1_3)## [1] 0 26myFixed = Model1$summary.fixed[,c(1:3,5)]
names(myFixed) = c("Mean", "SD", "Q0.025", "Q0.975")
PlotFrame = Pnts@data %>% select(state, Edge1_3, s_Edge1_3)
PlotFrame$ExtRipM = plogis(PlotFrame$s_Edge1_3*myFixed["Edge1_31","Mean"])
PlotFrame$ExtRipL = plogis(PlotFrame$s_Edge1_3*myFixed["Edge1_31","Q0.025"])
PlotFrame$ExtRipH = plogis(PlotFrame$s_Edge1_3*myFixed["Edge1_31","Q0.975"])
PlotFrame$State = "RM-LF"
PlotFrame2 = Pnts@data %>% select(state, Edge1_3, s_Edge1_3)
PlotFrame2$ExtRipM = plogis(PlotFrame2$s_Edge1_3*myFixed["Edge1_32","Mean"])
PlotFrame2$ExtRipL = plogis(PlotFrame2$s_Edge1_3*myFixed["Edge1_32","Q0.025"])
PlotFrame2$ExtRipH = plogis(PlotFrame2$s_Edge1_3*myFixed["Edge1_32","Q0.975"])
PlotFrame2$State = "LM-MF"
PlotFrame3 = Pnts@data %>% select(state, Edge1_3, s_Edge1_3)
PlotFrame3$ExtRipM = plogis(PlotFrame3$s_Edge1_3*myFixed["Edge1_33","Mean"])
PlotFrame3$ExtRipL = plogis(PlotFrame3$s_Edge1_3*myFixed["Edge1_33","Q0.025"])
PlotFrame3$ExtRipH = plogis(PlotFrame3$s_Edge1_3*myFixed["Edge1_33","Q0.975"])
PlotFrame3$State = "MM-HF"
CFrn = rbind(PlotFrame, PlotFrame2, PlotFrame3)
CFrn$State = factor(CFrn$State)
ggplot(CFrn, aes(Edge1_3, ExtRipM, group=State)) +
geom_line(aes(col = State,
linetype= State),
lwd = 1) +
scale_linetype_manual(values=c("solid", "longdash", "dotted")) +
scale_colour_manual(values=c("black", "black", "black")) +
geom_ribbon(aes(ymin=ExtRipL,
ymax=ExtRipH),
alpha=0.1, #transparency
linetype=1, #solid, dashed or other line types
colour="transparent", #border line color
size=1, #border line size
fill="gray30") +
geom_hline(yintercept = 0,
linetype = "solid",
col = "darkgray",
size = 0.5) +
ylab("Movement Probability") +
xlab("Open-Forest Edge") +
#ggtitle("Aerial Counts (COMPRECNTDEN)") +
theme_classic() +
scale_y_continuous(breaks = seq(0, 1, 0.2),
limits = c(0, 1)) +
scale_x_continuous(breaks = seq(0, 26, 2),
limits =c(0,26)) +
#facet_grid(rows = vars(State)) +
theme(plot.title = element_text(hjust = 0.5),
title = element_text(face="bold", size=18, hjust=0.5),
legend.text = element_text(size=14, face="bold"),
legend.direction = "vertical",
legend.position=c(0.15, 0.25),
legend.key.size = unit(3,"line"),
legend.title = element_text(size=18, face="bold"),
axis.text=element_text(size=12),
strip.text = element_text(face="bold", size = 20),
axis.title.y = element_text(face="bold", size = 22),
axis.text.x = element_text(face="bold", size=18, vjust=0.5,
hjust=0.5, angle=0),
axis.text.y = element_text(face="bold", size=18, vjust=0.5,
hjust=0.5, angle=0),
axis.title.x = element_text(face="bold", size = 22))
CFrn$Covariate = "c"
CFrn3 = CFrnProportions Combined
names(CFrn1)[2:3] = c("Value", "s_Value")
names(CFrn2)[2:3] = c("Value", "s_Value")
names(CFrn3)[2:3] = c("Value", "s_Value")
Comb.set = rbind(CFrn1, CFrn2, CFrn3)
Cmb.plt = ggplot(Comb.set , aes(Value, ExtRipM, group=State)) +
geom_line(aes(col = State,
linetype= State),
lwd = 1) +
scale_linetype_manual(values=c("solid", "longdash", "dotted")) +
scale_colour_manual(values=c("black", "black", "black")) +
geom_ribbon(aes(ymin=ExtRipL,
ymax=ExtRipH),
alpha=0.1, #transparency
linetype=1, #solid, dashed or other line types
colour="transparent", #border line color
size=1, #border line size
fill="gray30") +
geom_hline(yintercept = 0,
linetype = "solid",
col = "darkgray",
size = 0.5) +
ylab("Movement Probability") +
xlab(" ") +
#ggtitle("Aerial Counts (COMPRECNTDEN)") +
theme_classic() +
scale_y_continuous(breaks = seq(0, 1, 0.2),
limits = c(0, 1)) +
scale_x_continuous(breaks = scales::pretty_breaks(5), limits = c(0, NA)) +
#scale_x_continuous(breaks = seq(0, 100, 20),
# limits =c(0,100)) +
facet_wrap(~Covariate, ncol=1, scales="free") +
theme(plot.title = element_text(hjust = 0.5),
title = element_text(face="bold", size=18, hjust=0.5),
legend.text = element_text(size=14, face="bold"),
legend.direction = "horizontal",
#legend.position=c(0.15, 0.8),
legend.position="bottom",
legend.key.size = unit(3,"line"),
legend.title = element_text(size=20, face="bold"),
axis.text=element_text(size=12),
strip.text = element_text(face="bold", size = 22),
strip.background = element_blank(),
axis.title.y = element_text(face="bold", size = 22),
axis.text.x = element_text(face="bold", size=18, vjust=0.5,
hjust=0.5, angle=0),
axis.text.y = element_text(face="bold", size=18, vjust=0.5,
hjust=0.5, angle=0),
axis.title.x = element_text(face="bold", size = 22))
Cmb.plt
4.10 Atmospheric
Temperature
range(Pnts$HOURLYDRYBULBTEMPC)## [1] -33.3 32.8myFixed = Model1$summary.fixed[,c(1:3,5)]
names(myFixed) = c("Mean", "SD", "Q0.025", "Q0.975")
PlotFrame = Pnts@data %>% select(state, HOURLYDRYBULBTEMPC, s_HOURLYDRYBULBTEMPC)
PlotFrame$ExtRipM = plogis(PlotFrame$s_HOURLYDRYBULBTEMPC*myFixed["HOURLYDRYBULBTEMPC1","Mean"])
PlotFrame$ExtRipL = plogis(PlotFrame$s_HOURLYDRYBULBTEMPC*myFixed["HOURLYDRYBULBTEMPC1","Q0.025"])
PlotFrame$ExtRipH = plogis(PlotFrame$s_HOURLYDRYBULBTEMPC*myFixed["HOURLYDRYBULBTEMPC1","Q0.975"])
PlotFrame$State = "RM-LF"
Find.temp = PlotFrame %>%
filter(ExtRipM >= 0.49 &
ExtRipM <= 0.51 )
mean(Find.temp$HOURLYDRYBULBTEMPC)## [1] 7.515716PlotFrame2 = Pnts@data %>% select(state, HOURLYDRYBULBTEMPC, s_HOURLYDRYBULBTEMPC)
PlotFrame2$ExtRipM = plogis(PlotFrame2$s_HOURLYDRYBULBTEMPC*myFixed["HOURLYDRYBULBTEMPC2","Mean"])
PlotFrame2$ExtRipL = plogis(PlotFrame2$s_HOURLYDRYBULBTEMPC*myFixed["HOURLYDRYBULBTEMPC2","Q0.025"])
PlotFrame2$ExtRipH = plogis(PlotFrame2$s_HOURLYDRYBULBTEMPC*myFixed["HOURLYDRYBULBTEMPC2","Q0.975"])
PlotFrame2$State = "LM-MF"
PlotFrame3 = Pnts@data %>% select(state, HOURLYDRYBULBTEMPC, s_HOURLYDRYBULBTEMPC)
PlotFrame3$ExtRipM = plogis(PlotFrame3$s_HOURLYDRYBULBTEMPC*myFixed["HOURLYDRYBULBTEMPC3","Mean"])
PlotFrame3$ExtRipL = plogis(PlotFrame3$s_HOURLYDRYBULBTEMPC*myFixed["HOURLYDRYBULBTEMPC3","Q0.025"])
PlotFrame3$ExtRipH = plogis(PlotFrame3$s_HOURLYDRYBULBTEMPC*myFixed["HOURLYDRYBULBTEMPC3","Q0.975"])
PlotFrame3$State = "MM-HF"
CFrn = rbind(PlotFrame, PlotFrame2, PlotFrame3)
CFrn$State = factor(CFrn$State)
ggplot(CFrn, aes(HOURLYDRYBULBTEMPC, ExtRipM, group=State)) +
geom_line(aes(col = State,
linetype= State),
lwd = 1) +
scale_linetype_manual(values=c("solid", "longdash", "dotted")) +
scale_colour_manual(values=c("black", "black", "black")) +
geom_ribbon(aes(ymin=ExtRipL,
ymax=ExtRipH),
alpha=0.1, #transparency
linetype=1, #solid, dashed or other line types
colour="transparent", #border line color
size=1, #border line size
fill="gray30") +
geom_hline(yintercept = 0,
linetype = "solid",
col = "darkgray",
size = 0.5) +
ylab("Movement Probability") +
xlab("Temperature (戼㸰C)") +
#ggtitle("Aerial Counts (COMPRECNTDEN)") +
theme_classic() +
scale_y_continuous(breaks = seq(0, 1, 0.2),
limits = c(0, 1)) +
scale_x_continuous(breaks = seq(-35, 35, 5),
limits =c(-35,35)) +
#facet_grid(rows = vars(State)) +
theme(plot.title = element_text(hjust = 0.5),
title = element_text(face="bold", size=18, hjust=0.5),
legend.text = element_text(size=14, face="bold"),
legend.direction = "vertical",
legend.position=c(0.7, 0.25),
legend.key.size = unit(3,"line"),
legend.title = element_text(size=18, face="bold"),
axis.text=element_text(size=12),
strip.text = element_text(face="bold", size = 20),
axis.title.y = element_text(face="bold", size = 22),
axis.text.x = element_text(face="bold", size=18, vjust=0.5,
hjust=0.5, angle=0),
axis.text.y = element_text(face="bold", size=18, vjust=0.5,
hjust=0.5, angle=0),
axis.title.x = element_text(face="bold", size = 22))
CFrn$Covariate = "a"
CFrn2 = CFrnPressure
range(Pnts$HOURLYStationPressure)## [1] 27.76 29.48myFixed = Model1$summary.fixed[,c(1:3,5)]
names(myFixed) = c("Mean", "SD", "Q0.025", "Q0.975")
PlotFrame = Pnts@data %>% select(state, HOURLYStationPressure, s_HOURLYStationPressure)
PlotFrame$ExtRipM = plogis(PlotFrame$s_HOURLYStationPressure*myFixed["HOURLYStationPressure1","Mean"])
PlotFrame$ExtRipL = plogis(PlotFrame$s_HOURLYStationPressure*myFixed["HOURLYStationPressure1","Q0.025"])
PlotFrame$ExtRipH = plogis(PlotFrame$s_HOURLYStationPressure*myFixed["HOURLYStationPressure1","Q0.975"])
PlotFrame$State = "RM-LF"
Find.temp = PlotFrame %>%
filter(ExtRipM >= 0.49 &
ExtRipM <= 0.51 )
mean(Find.temp$HOURLYStationPressure)## [1] 28.78259PlotFrame2 = Pnts@data %>% select(state, HOURLYStationPressure, s_HOURLYStationPressure)
PlotFrame2$ExtRipM = plogis(PlotFrame2$s_HOURLYStationPressure*myFixed["HOURLYStationPressure2","Mean"])
PlotFrame2$ExtRipL = plogis(PlotFrame2$s_HOURLYStationPressure*myFixed["HOURLYStationPressure2","Q0.025"])
PlotFrame2$ExtRipH = plogis(PlotFrame2$s_HOURLYStationPressure*myFixed["HOURLYStationPressure2","Q0.975"])
PlotFrame2$State = "LM-MF"
PlotFrame3 = Pnts@data %>% select(state, HOURLYStationPressure, s_HOURLYStationPressure)
PlotFrame3$ExtRipM = plogis(PlotFrame3$s_HOURLYStationPressure*myFixed["HOURLYStationPressure3","Mean"])
PlotFrame3$ExtRipL = plogis(PlotFrame3$s_HOURLYStationPressure*myFixed["HOURLYStationPressure3","Q0.025"])
PlotFrame3$ExtRipH = plogis(PlotFrame3$s_HOURLYStationPressure*myFixed["HOURLYStationPressure3","Q0.975"])
PlotFrame3$State = "MM-HF"
CFrn = rbind(PlotFrame, PlotFrame2, PlotFrame3)
CFrn$State = factor(CFrn$State)
ggplot(CFrn, aes(HOURLYStationPressure, ExtRipM, group=State)) +
geom_line(aes(col = State,
linetype= State),
lwd = 1) +
scale_linetype_manual(values=c("solid", "longdash", "dotted")) +
scale_colour_manual(values=c("black", "black", "black")) +
geom_ribbon(aes(ymin=ExtRipL,
ymax=ExtRipH),
alpha=0.1, #transparency
linetype=1, #solid, dashed or other line types
colour="transparent", #border line color
size=1, #border line size
fill="gray30") +
geom_hline(yintercept = 0,
linetype = "solid",
col = "darkgray",
size = 0.5) +
ylab("Movement Probability") +
xlab("Pressure (in)") +
#ggtitle("Aerial Counts (COMPRECNTDEN)") +
theme_classic() +
scale_y_continuous(breaks = seq(0, 1, 0.2),
limits = c(0, 1)) +
scale_x_continuous(breaks = seq(27.75, 29.5, 0.5),
limits =c(27.75,29.5)) +
#facet_grid(rows = vars(State)) +
theme(plot.title = element_text(hjust = 0.5),
title = element_text(face="bold", size=18, hjust=0.5),
legend.text = element_text(size=14, face="bold"),
legend.direction = "vertical",
legend.position=c(0.6, 0.25),
legend.key.size = unit(3,"line"),
legend.title = element_text(size=18, face="bold"),
axis.text=element_text(size=12),
strip.text = element_text(face="bold", size = 20),
axis.title.y = element_text(face="bold", size = 22),
axis.text.x = element_text(face="bold", size=18, vjust=0.5,
hjust=0.5, angle=0),
axis.text.y = element_text(face="bold", size=18, vjust=0.5,
hjust=0.5, angle=0),
axis.title.x = element_text(face="bold", size = 22))
CFrn$Covariate = "b"
CFrn2 = CFrnProportions Combined
names(CFrn1)[2:3] = c("Value", "s_Value")
names(CFrn2)[2:3] = c("Value", "s_Value")
Comb.set = rbind(CFrn1, CFrn2)
Cmb.plt = ggplot(Comb.set , aes(Value, ExtRipM, group=State)) +
geom_line(aes(col = State,
linetype= State),
lwd = 1) +
scale_linetype_manual(values=c("solid", "longdash", "dotted")) +
scale_colour_manual(values=c("black", "black", "black")) +
geom_ribbon(aes(ymin=ExtRipL,
ymax=ExtRipH),
alpha=0.1, #transparency
linetype=1, #solid, dashed or other line types
colour="transparent", #border line color
size=1, #border line size
fill="gray30") +
geom_hline(yintercept = 0,
linetype = "solid",
col = "darkgray",
size = 0.5) +
ylab("Movement Probability") +
xlab(" ") +
theme_classic() +
scale_y_continuous(breaks = seq(0, 1, 0.2),
limits = c(0, 1)) +
scale_x_continuous(breaks = scales::pretty_breaks(5), limits = c(NA, NA)) +
facet_wrap(~Covariate, ncol=1, scales="free") +
theme(plot.title = element_text(hjust = 0.5),
title = element_text(face="bold", size=18, hjust=0.5),
legend.text = element_text(size=14, face="bold"),
legend.direction = "horizontal",
#legend.position=c(0.15, 0.8),
legend.position="bottom",
legend.key.size = unit(3,"line"),
legend.title = element_text(size=20, face="bold"),
axis.text=element_text(size=12),
strip.text = element_text(face="bold", size = 22),
strip.background = element_blank(),
axis.title.y = element_text(face="bold", size = 22),
axis.text.x = element_text(face="bold", size=18, vjust=0.5,
hjust=0.5, angle=0),
axis.text.y = element_text(face="bold", size=18, vjust=0.5,
hjust=0.5, angle=0),
axis.title.x = element_text(face="bold", size = 22))
Cmb.plt
4.11 Validation
Predict the probabilty for each state based on testing set variables. This is true prediction as the spatial and temporal error is not added in to the predictions. However, random effects year, day, hour, and cohort are considered for prediction as these may reflect behavior.
ModResult = Model1
Pred.pnts = SpatialPointsDataFrame(wpmovetest[, c("x","y")], wpmovetest)
proj4string(Pred.pnts) = nProj
Pred.pnts = spTransform(Pred.pnts, proj4string(DomainP))
#pLoc = cbind(Pred.pnts@coords[,1], Pred.pnts@coords[,2])
#At = inla.spde.make.A(mesh, loc=pLoc)
mydf = ModResult$summary.fix
###############################################################Predict Probability of State 1
Pred.pnts$Int1 = ModResult$summary.fix[1,1]
#Find Year
Year.tab = ModResult$summary.random$year
Pred.pnts$Year.cov = with(Year.tab,
mean[match(
Pred.pnts$year,
ID)])
#Find day
LU.tab = ModResult$summary.random$day1
Pred.pnts$Day.cov = with(LU.tab,
mean[match(
Pred.pnts$jdate,
ID)])
#Find hour
LU.tab = ModResult$summary.random$hour1
Pred.pnts$Hour.cov = with(LU.tab,
mean[match(
Pred.pnts$hour,
ID)])
#Find cohort
LU.tab = ModResult$summary.random$cohort1
Pred.pnts$Cohort.cov = with(LU.tab,
mean[match(
Pred.pnts$cohort,
ID)])
Pred.pnts@data = Pred.pnts@data %>%
mutate(LP = Year.cov + Day.cov + Hour.cov + Cohort.cov +
s_Dvlpd_A*mydf["Dvlpd_A1","mean"] +
s_Rprn_Ar*mydf["Rprn_Ar1","mean"] +
s_Open_Ar*mydf["Open_Ar1","mean"] +
s_Hrdms_A*mydf["Hrdms_A1","mean"] +
s_PCI_2*mydf["PCI_21","mean"] +
s_PCI_5*mydf["PCI_51","mean"] +
s_Edge1_3*mydf["Edge1_31","mean"] +
s_HOURLYDRYBULBTEMPC*mydf["HOURLYDRYBULBTEMPC1","mean"] +
s_HOURLYStationPressure*mydf["HOURLYStationPressure1","mean"],
Pred.S1 = plogis(LP))
##########################################################Predict Probability of State 2
Pred.pnts$Int2 = ModResult$summary.fix[2,1]
#Find Year
Year.tab = ModResult$summary.random$year
Pred.pnts$Year.cov = with(Year.tab,
mean[match(
Pred.pnts$year,
ID)])
#Find day
LU.tab = ModResult$summary.random$day2
Pred.pnts$Day.cov = with(LU.tab,
mean[match(
Pred.pnts$jdate,
ID)])
#Find hour
LU.tab = ModResult$summary.random$hour2
Pred.pnts$Hour.cov = with(LU.tab,
mean[match(
Pred.pnts$hour,
ID)])
#Find cohort
LU.tab = ModResult$summary.random$cohort2
Pred.pnts$Cohort.cov = with(LU.tab,
mean[match(
Pred.pnts$cohort,
ID)])
Pred.pnts@data = Pred.pnts@data %>%
mutate(LP = Year.cov + Day.cov + Hour.cov + Cohort.cov +
s_Dvlpd_A*mydf["Dvlpd_A2","mean"] +
s_Rprn_Ar*mydf["Rprn_Ar2","mean"] +
s_Open_Ar*mydf["Open_Ar2","mean"] +
s_Hrdms_A*mydf["Hrdms_A2","mean"] +
s_PCI_2*mydf["PCI_22","mean"] +
s_PCI_5*mydf["PCI_52","mean"] +
s_Edge1_3*mydf["Edge1_32","mean"] +
s_HOURLYDRYBULBTEMPC*mydf["HOURLYDRYBULBTEMPC2","mean"] +
s_HOURLYStationPressure*mydf["HOURLYStationPressure2","mean"],
Pred.S2 = plogis(LP))
##########################################################Predict Probability of State 3
Pred.pnts$Int3 = ModResult$summary.fix[3,1]
#Find Year
Year.tab = ModResult$summary.random$year
Pred.pnts$Year.cov = with(Year.tab,
mean[match(
Pred.pnts$year,
ID)])
#Find day
LU.tab = ModResult$summary.random$day3
Pred.pnts$Day.cov = with(LU.tab,
mean[match(
Pred.pnts$jdate,
ID)])
#Find hour
LU.tab = ModResult$summary.random$hour3
Pred.pnts$Hour.cov = with(LU.tab,
mean[match(
Pred.pnts$hour,
ID)])
#Find cohort
LU.tab = ModResult$summary.random$cohort3
Pred.pnts$Cohort.cov = with(LU.tab,
mean[match(
Pred.pnts$cohort,
ID)])
Pred.pnts@data = Pred.pnts@data %>%
mutate(LP = Year.cov + Day.cov + Hour.cov + Cohort.cov +
s_Dvlpd_A*mydf["Dvlpd_A3","mean"] +
s_Rprn_Ar*mydf["Rprn_Ar3","mean"] +
s_Open_Ar*mydf["Open_Ar3","mean"] +
s_Hrdms_A*mydf["Hrdms_A3","mean"] +
s_PCI_2*mydf["PCI_23","mean"] +
s_PCI_5*mydf["PCI_53","mean"] +
s_Edge1_3*mydf["Edge1_33","mean"] +
s_HOURLYDRYBULBTEMPC*mydf["HOURLYDRYBULBTEMPC3","mean"] +
s_HOURLYStationPressure*mydf["HOURLYStationPressure3","mean"],
Pred.S3 = plogis(LP))
###########
#AUC State 1
Pred.pnts$S1.OBS = ifelse(Pred.pnts$state == 1, 1, 0)
Test.state1 = as.data.frame(cbind(1:dim(Pred.pnts@data)[1], Pred.pnts$S1.OBS, Pred.pnts$Pred.S1))
names(Test.state1) = c("ID", "OBS", "Pred")
Thresh = optimal.thresholds(Test.state1, opt.methods = c("MaxSens+Spec"))
Thresh## Method Pred
## 1 MaxSens+Spec 0.47Test.state1.out = presence.absence.accuracy(Test.state1, threshold = Thresh[1,2])
Test.state1.out$Count = sum(Pred.pnts$S1.OBS)
Test.state1.out$State = "RM-LF"
#AUC State 2
Pred.pnts$S2.OBS = ifelse(Pred.pnts$state == 2, 1, 0)
Test.state2 = as.data.frame(cbind(1:dim(Pred.pnts@data)[1], Pred.pnts$S2.OBS, Pred.pnts$Pred.S2))
names(Test.state2) = c("ID", "OBS", "Pred")
Thresh = optimal.thresholds(Test.state2, opt.methods = c("MaxSens+Spec"))
Thresh## Method Pred
## 1 MaxSens+Spec 0.13Test.state2.out = presence.absence.accuracy(Test.state2, threshold = Thresh[1,2])
Test.state2.out$Count = sum(Pred.pnts$S2.OBS)
Test.state2.out$State = "LM-MF"
#AUC State 3
Pred.pnts$S3.OBS = ifelse(Pred.pnts$state == 3, 1, 0)
Test.state3 = as.data.frame(cbind(1:dim(Pred.pnts@data)[1], Pred.pnts$S3.OBS, Pred.pnts$Pred.S3))
names(Test.state3) = c("ID", "OBS", "Pred")
Thresh = optimal.thresholds(Test.state3, opt.methods = c("MaxSens+Spec"))
Thresh## Method Pred
## 1 MaxSens+Spec 0.56Test.state3.out = presence.absence.accuracy(Test.state3, threshold = Thresh[1,2])
Test.state3.out$Count = sum(Pred.pnts$S3.OBS)
Test.state3.out$State = "MM-HF"Validation Table
Validation.tab = rbind(Test.state1.out, Test.state2.out, Test.state3.out) %>%
select(State, Count, threshold, PCC, sensitivity, specificity, AUC)
names(Validation.tab) = c("State", "Count", "Threshold", "PCC", "Sensitivity", "specificity", "AUC")
for(i in 4:7){
Validation.tab[,i] = round(Validation.tab[,i], 2)
}
#SelMod.fx = xtable(Validation.tab)
#print(SelMod.fx, include.rownames = F)
kable(Validation.tab) %>%
kable_styling("striped", full_width = F) %>%
row_spec(0, font_size = 20) | State | Count | Threshold | PCC | Sensitivity | specificity | AUC |
|---|---|---|---|---|---|---|
| RM-LF | 1106 | 0.47 | 0.60 | 0.78 | 0.54 | 0.68 |
| LM-MF | 2300 | 0.13 | 0.66 | 0.86 | 0.47 | 0.71 |
| MM-HF | 1287 | 0.56 | 0.71 | 0.57 | 0.76 | 0.70 |