Pig Spatiotemporal Movement Model
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)2 Load and Prep 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.2, replace=F)
wpmovetest = as.data.frame(wpmovetest)
dim(wpmovetest) #7981## [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) #31933## [1] 18780 313 Set Time steps
wpmovetrain %>%
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),
tz = "America/New_York",
format = "%m/%d/%Y %H:%M")
Pig.lvls = unique(wpmovetrain$id)
for(i in 1:length(Pig.lvls)){
tmp.pig = wpmovetrain %>% filter(id == Pig.lvls[i])
tmp.pig = arrange(tmp.pig, Date)
tmp.pig$TStep = 1:nrow(tmp.pig)
tmp.pig$Pig = i
if(i == 1){Out.pig = tmp.pig}
else{Out.pig = rbind(Out.pig, tmp.pig)}
}
wpmovetrain = Out.pig4 Quick plot
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)
5 Domain
max.edge = 1 #Make the outer edge length km
bound.outer = 15 #Outer extension km
Hull = gBuffer(
gConvexHull(Pnts),
width = 5)
plot(Hull)
Hull = gBuffer(Hull, width = 10)
plot(Hull)
bdry = inla.sp2segment(Hull) #Formatting boundary
mesh = inla.mesh.2d(boundary = bdry,
loc = Pnts, #Fit to point locations
max.edge = c(1, 4)*max.edge,
cutoff = 0.75,
min.angle = 23,
offset = c(max.edge, bound.outer))
mesh$n #number of nodes## [1] 5016plot(mesh)
plot(Pnts, col="red", add=T)
6 Projection Matricies
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(100, 0.5),
prior.sigma = c(1, 0.5))
#Create index to track locations of mesh nodes
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) 7 Organize Data (Full Model)
pnts.df = Pnts@data
pnts.df$Point.IID = 1:nrow(pnts.df)
#Movement State 1
HR1.lst = list(c(field1,
list(intercept1 = 1)), #intercept
list(Dvlpd_A1 = round(pnts.df[,"s_Dvlpd_A"], 3),
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),
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)
sum(pnts.df$State1)## [1] 4531#Test version
TS1.stk = inla.stack(data = list(Y = pnts.df$State1),
A = list(A.mx1, 1),
effects = HR1.lst,
tag = "T1.0")
S1.stk = inla.stack(data = list(Y = cbind(pnts.df$State1, NA, NA)),
A = list(A.mx1, 1),
effects = HR1.lst,
tag = "s1.0")
###State 2
HR2.lst = list(c(field2,
field1x,
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),
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)
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)),
A = list(A.mx2, 1),
effects = HR2.lst,
tag = "s2.0")
###State 3
HR3.lst = list(c(field3,
field1xx,
field2x,
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),
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)
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")
Joint.stk = inla.stack(S1.stk, S2.stk, S3.stk)
#save(list=c("Joint.stk", "spde"), file="./ST_060119/HPC/PigM_060519.RData") #bundled for HPCC8 Full 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 +
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(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
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 = theta.m1),
control.results = list(return.marginals.random = TRUE,
return.marginals.predictor = TRUE),
control.compute=list(dic = TRUE, cpo = TRUE, waic = TRUE))9 Load Results from HPCC
load("~/Michigan_State/Steve/ST_060119/HPC/PigM_RES_060519.RData")Get Fit values bi tier
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]
cor(pnts.df[,c("S1_pr","S2_pr","S3_pr")])## S1_pr S2_pr S3_pr
## S1_pr 1.0000000 -0.6667208 -0.3624334
## S2_pr -0.6667208 1.0000000 -0.4483404
## S3_pr -0.3624334 -0.4483404 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("State 1", "State 2", "State 3"))
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 = 22)) 
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 = "State1"
mic.d2 = as.data.frame(Model1$summary.random$hour2)
names(mic.d2) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975", "mode", "kld")
mic.d2$State = "State2"
mic.d3 = as.data.frame(Model1$summary.random$hour3)
names(mic.d3) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975", "mode", "kld")
mic.d3$State = "State3"
mic.d = rbind(mic.d1, mic.d2, mic.d3)
mic.d$State = ordered(factor(mic.d$State), levels = c("State1", "State2", "State3"))
Raw.plt2 = ggplot(mic.d, aes(ID, Mean)) +
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") +
geom_hline(yintercept = 0,
linetype = "solid",
col = "darkgray",
size = 0.5) +
ylab("Movement State Probability (logit)") +
xlab("Hour of Day") +
#ggtitle("Aerial Counts (COMPRECNTDEN)") +
theme_classic() +
scale_y_continuous(breaks = seq(-1, 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 = 20),
axis.title.y = element_text(face="bold", size = 20),
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 = 20))
Raw.plt2
10 Cohort
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 = "State 1"
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 = "State 2"
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 = "State 3"
Comb.plt.df = rbind(Comb.plt.df1, Comb.plt.df2, Comb.plt.df3)
Plt.all = ggplot(Comb.plt.df , aes(ID, Mean)) +
geom_point(size=2, pch=19, col = "red") +
geom_linerange(aes(ymin=Q025, ymax=Q975), colour="black") +
geom_hline(yintercept = 0,
linetype = "dotted",
colour = "darkgray",
size = 0.5) +
geom_hline(yintercept = 0,
linetype = "solid",
col = "lightgray",
size = 0.5) +
facet_wrap(~State, ncol=1) +
xlab("Cohort Group") +
ylab("State Probability (logit)") +
theme_classic() +
theme(axis.text=element_text(size=16),
legend.title = element_text(size=16, face="bold"),
legend.key = element_rect(fill = "gray80", linetype=0),
legend.background = element_rect(fill = "gray80", linetype=0),
strip.text = element_text(face="bold", size = 20),
axis.title.y = element_text(face="bold", size = 20),
axis.title.x = element_text(face="bold", size = 20, vjust=-2))
Plt.all
Covariate Tables
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)
FE.table## Covariate State Mean
## intercept1 intercept State 1 -4.46830553
## intercept2 intercept State 2 -6.85475532
## intercept3 intercept State 3 0.64227457
## Dvlpd_A1 Dvlpd_A State 1 0.01490958
## Dvlpd_A2 Dvlpd_A State 2 0.00238311
## Dvlpd_A3 Dvlpd_A State 3 0.01974940
## Rprn_Ar1 Rprn_Ar State 1 0.35175384
## Rprn_Ar2 Rprn_Ar State 2 0.24100972
## Rprn_Ar3 Rprn_Ar State 3 -0.45422603
## Open_Ar1 Open_Ar State 1 0.46187742
## Open_Ar2 Open_Ar State 2 -0.03504699
## Open_Ar3 Open_Ar State 3 -0.35144676
## Hrdms_A1 Hrdms_A State 1 -0.16876482
## Hrdms_A2 Hrdms_A State 2 -0.06252680
## Hrdms_A3 Hrdms_A State 3 0.08445231
## PCI_21 PCI_2 State 1 0.34098407
## PCI_22 PCI_2 State 2 0.41244997
## PCI_23 PCI_2 State 3 -0.41454738
## PCI_51 PCI_5 State 1 0.20332383
## PCI_52 PCI_5 State 2 -0.02283126
## PCI_53 PCI_5 State 3 -0.11052205
## Edge1_31 Edge1_3 State 1 -0.12853918
## Edge1_32 Edge1_3 State 2 0.03537606
## Edge1_33 Edge1_3 State 3 0.05945077
## HOURLYDRYBULBTEMPC1 HOURLYDRYBULBTEMPC State 1 0.84565593
## HOURLYDRYBULBTEMPC2 HOURLYDRYBULBTEMPC State 2 0.07948688
## HOURLYDRYBULBTEMPC3 HOURLYDRYBULBTEMPC State 3 -0.15114193
## HOURLYStationPressure1 HOURLYStationPressure State 1 0.16217375
## HOURLYStationPressure2 HOURLYStationPressure State 2 0.09858855
## HOURLYStationPressure3 HOURLYStationPressure State 3 -0.09672927
## sd Q025 Q975
## intercept1 2.42235370 -9.224202639 0.2836225524
## intercept2 2.98861327 -12.722410899 -0.9919965939
## intercept3 2.37286022 -4.016450153 5.2971113439
## Dvlpd_A1 0.06836815 -0.119320160 0.1490272995
## Dvlpd_A2 0.04558388 -0.087113407 0.0918049380
## Dvlpd_A3 0.03645204 -0.051818243 0.0912573251
## Rprn_Ar1 0.14166125 0.073625037 0.6296505384
## Rprn_Ar2 0.10062947 0.043440138 0.4384144262
## Rprn_Ar3 0.07760500 -0.606590814 -0.3019884094
## Open_Ar1 0.09615774 0.273087339 0.6505099367
## Open_Ar2 0.09482482 -0.221220078 0.1509707300
## Open_Ar3 0.06130920 -0.471817399 -0.2311765792
## Hrdms_A1 0.07767469 -0.321266427 -0.0163904735
## Hrdms_A2 0.04703719 -0.154876655 0.0297459930
## Hrdms_A3 0.03706168 0.011687726 0.1571561595
## PCI_21 0.19172646 -0.035439617 0.7170936205
## PCI_22 0.09725517 0.221505275 0.6032353211
## PCI_23 0.08269579 -0.576907093 -0.2523231681
## PCI_51 0.09948732 0.007996663 0.3984879772
## PCI_52 0.07784755 -0.175672247 0.1298821821
## PCI_53 0.05682338 -0.222085500 0.0009483042
## Edge1_31 0.07636023 -0.278460052 0.0212565721
## Edge1_32 0.05311142 -0.068899565 0.1395646679
## Edge1_33 0.04210585 -0.023217211 0.1420497533
## HOURLYDRYBULBTEMPC1 0.31197519 0.233143436 1.4576572551
## HOURLYDRYBULBTEMPC2 0.13949412 -0.194387106 0.3531323000
## HOURLYDRYBULBTEMPC3 0.09555317 -0.338745026 0.0363046008
## HOURLYStationPressure1 0.25183304 -0.332259426 0.6561942926
## HOURLYStationPressure2 0.11077561 -0.118901330 0.3158969208
## HOURLYStationPressure3 0.06972872 -0.233630254 0.0400574692#Random
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)
FE.table## Covariate Mean sd
## Range for field1 Range for field 3.11301286 0.509522753
## Stdev for field1 Stdev for field 13.36639556 1.699828251
## Range for field2 Range for field 3.93447353 0.765107274
## Stdev for field2 Stdev for field 10.95012036 1.610804643
## Range for field3 Range for field 16.99140216 9.821446081
## Stdev for field3 Stdev for field 3.01926077 1.101838319
## Precision for pig Precision for pi 1.51463939 0.649941824
## Precision for year Precision for yea 1.24347016 0.517243614
## Precision for cohort1 Precision for cohort 0.45836106 0.181021415
## Precision for cohort2 Precision for cohort 0.44878166 0.168280772
## Precision for cohort3 Precision for cohort 1.23450030 0.530543232
## Precision for day1 Precision for day 1.20296032 0.446170328
## Precision for hour1 Precision for hour 1.09486136 0.348474040
## Precision for day2 Precision for day 1.69246010 0.493644150
## Precision for hour2 Precision for hour 1.69772818 0.478784265
## Precision for day3 Precision for day 4.31021971 1.123828882
## Precision for hour3 Precision for hour 1.62850697 0.461608915
## Precision for TStep1 Precision for TStep 0.03233787 0.004314471
## Rho for TStep1 Rho for TStep 0.96916888 0.004010085
## Precision for TStep2 Precision for TStep 0.12244218 0.012503031
## Rho for TStep2 Rho for TStep 0.88548624 0.008759658
## Precision for TStep3 Precision for TStep 0.19885495 0.016614916
## Rho for TStep3 Rho for TStep 0.78040416 0.010555067
## Beta for field1x Beta for field1 0.69071527 0.121176378
## Beta for field1xx Beta for field1x -0.84613170 0.067890392
## Beta for field2x Beta for field2 -0.39377607 0.054458272
## Q025 Q975
## Range for field1 2.27393762 4.26359971
## Stdev for field1 10.39517755 17.06767062
## Range for field2 2.59283479 5.59054999
## Stdev for field2 8.00959469 14.33086063
## Range for field3 6.29331460 42.84199086
## Stdev for field3 1.50814670 5.76461095
## Precision for pig 0.62456880 3.13281393
## Precision for year 0.50251270 2.49935702
## Precision for cohort1 0.20960749 0.90860346
## Precision for cohort2 0.20823729 0.85735943
## Precision for cohort3 0.50568978 2.55000079
## Precision for day1 0.53444398 2.26311748
## Precision for hour1 0.58193279 1.93612754
## Precision for day2 0.92475091 2.84853229
## Precision for hour2 0.90281641 2.76303094
## Precision for day3 2.46566132 6.84702555
## Precision for hour3 0.87684679 2.67323657
## Precision for TStep1 0.02492093 0.04186309
## Rho for TStep1 0.96042780 0.97616393
## Precision for TStep2 0.10248290 0.15111207
## Rho for TStep2 0.86714485 0.90160384
## Precision for TStep3 0.16676308 0.23197039
## Rho for TStep3 0.75816279 0.79964453
## Beta for field1x 0.44672568 0.92290415
## Beta for field1xx -0.98148559 -0.71429655
## Beta for field2x -0.50648055 -0.29292481Proportion Emergent
range(Pnts$Open_Ar)## [1] 0 100PlotFrame = Pnts@data %>% select(state, Open_Ar, s_Open_Ar)
PlotFrame$ExtRipM = plogis(PlotFrame$s_Open_Ar*0.46187742)
PlotFrame$ExtRipL = plogis(PlotFrame$s_Open_Ar*(0.46187742 - 0.09615774))
PlotFrame$ExtRipH = plogis(PlotFrame$s_Open_Ar*(0.46187742 + 0.09615774))
PlotFrame$State = "State 1"
PlotFrame2 = Pnts@data %>% select(state, Open_Ar, s_Open_Ar)
PlotFrame2$ExtRipM = plogis(PlotFrame2$s_Open_Ar* 0.03504699)
PlotFrame2$ExtRipL = plogis(PlotFrame2$s_Open_Ar*(0.03504699 - 0.09482482))
PlotFrame2$ExtRipH = plogis(PlotFrame2$s_Open_Ar*(0.03504699 + 0.09482482))
PlotFrame2$State = "State 2"
PlotFrame3 = Pnts@data %>% select(state, Open_Ar, s_Open_Ar)
PlotFrame3$ExtRipM = 1 - plogis(PlotFrame3$s_Open_Ar*0.35144676)
PlotFrame3$ExtRipL = 1 - plogis(PlotFrame3$s_Open_Ar*(0.35144676 - 0.06130920))
PlotFrame3$ExtRipH = 1 - plogis(PlotFrame3$s_Open_Ar*(0.35144676 + 0.06130920))
PlotFrame3$State = "State 3"
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),
#method = "lm",
#span = 0.3,
#se = FALSE,
lwd = 1) +
scale_linetype_manual(values=c("solid", "longdash", "dotted")) +
scale_colour_manual(values=c("red", "blue", "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="darkgray") +
geom_hline(yintercept = 0,
linetype = "solid",
col = "darkgray",
size = 0.5) +
ylab("Movement Probability") +
xlab("Proportion Emergent/Herbaceous (Open_Ar)") +
#ggtitle("Aerial Counts (COMPRECNTDEN)") +
theme_classic() +
scale_y_continuous(breaks = seq(0, 1, 0.1),
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),
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=18, vjust=0.5,
hjust=0.5, angle=0),
axis.title.x = element_text(face="bold", size = 20))
Proportion Riparian
range(Pnts$Rprn_Ar)## [1] 0 100PlotFrame = Pnts@data %>% select(state, Rprn_Ar, s_Rprn_Ar)
PlotFrame$ExtRipM = plogis(PlotFrame$s_Rprn_Ar*0.35175384)
PlotFrame$ExtRipL = plogis(PlotFrame$s_Rprn_Ar*(0.35175384 - 0.14166125))
PlotFrame$ExtRipH = plogis(PlotFrame$s_Rprn_Ar*(0.35175384 + 0.14166125))
PlotFrame$State = "State 1"
PlotFrame2 = Pnts@data %>% select(state, Rprn_Ar, s_Rprn_Ar)
PlotFrame2$ExtRipM = plogis(PlotFrame2$s_Rprn_Ar* 0.24100972)
PlotFrame2$ExtRipL = plogis(PlotFrame2$s_Rprn_Ar*(0.24100972 - 0.10062947))
PlotFrame2$ExtRipH = plogis(PlotFrame2$s_Rprn_Ar*(0.24100972 + 0.10062947))
PlotFrame2$State = "State 2"
PlotFrame3 = Pnts@data %>% select(state, Rprn_Ar, s_Rprn_Ar)
PlotFrame3$ExtRipM = 1 - plogis(PlotFrame3$s_Rprn_Ar*0.45422603)
PlotFrame3$ExtRipL = 1 - plogis(PlotFrame3$s_Rprn_Ar*(0.45422603 - 0.07760500))
PlotFrame3$ExtRipH = 1 - plogis(PlotFrame3$s_Rprn_Ar*(0.45422603 + 0.07760500))
PlotFrame3$State = "State 3"
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),
#method = "lm",
#span = 0.3,
#se = FALSE,
lwd = 1) +
scale_linetype_manual(values=c("solid", "longdash", "dotted")) +
scale_colour_manual(values=c("red", "blue", "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="darkgray") +
geom_hline(yintercept = 0,
linetype = "solid",
col = "darkgray",
size = 0.5) +
ylab("Movement Probability") +
xlab("Proportion Riparian (Rprn_Ar)") +
#ggtitle("Aerial Counts (COMPRECNTDEN)") +
theme_classic() +
scale_y_continuous(breaks = seq(0, 1, 0.1),
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),
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=18, vjust=0.5,
hjust=0.5, angle=0),
axis.title.x = element_text(face="bold", size = 20))
Proportion Developed
range(Pnts$Dvlpd_A)## [1] 0.0000 99.3719PlotFrame = Pnts@data %>% select(state, Dvlpd_A, s_Dvlpd_A)
PlotFrame$ExtRipM = plogis(PlotFrame$s_Dvlpd_A*0.01490958)
PlotFrame$ExtRipL = plogis(PlotFrame$s_Dvlpd_A*(0.01490958 - 0.06836815))
PlotFrame$ExtRipH = plogis(PlotFrame$s_Dvlpd_A*(0.01490958 + 0.06836815))
PlotFrame$State = "State 1"
PlotFrame2 = Pnts@data %>% select(state, Dvlpd_A, s_Dvlpd_A)
PlotFrame2$ExtRipM = plogis(PlotFrame2$s_Dvlpd_A* 0.00238311)
PlotFrame2$ExtRipL = plogis(PlotFrame2$s_Dvlpd_A*(0.00238311 - 0.04558388))
PlotFrame2$ExtRipH = plogis(PlotFrame2$s_Dvlpd_A*(0.00238311 + 0.04558388))
PlotFrame2$State = "State 2"
PlotFrame3 = Pnts@data %>% select(state, Dvlpd_A, s_Dvlpd_A)
PlotFrame3$ExtRipM = 1 - plogis(PlotFrame3$s_Dvlpd_A*0.01974940)
PlotFrame3$ExtRipL = 1 - plogis(PlotFrame3$s_Dvlpd_A*(0.01974940 - 0.03645204))
PlotFrame3$ExtRipH = 1 - plogis(PlotFrame3$s_Dvlpd_A*(0.01974940 + 0.03645204))
PlotFrame3$State = "State 3"
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),
#method = "lm",
#span = 0.3,
#se = FALSE,
lwd = 1) +
scale_linetype_manual(values=c("solid", "longdash", "dotted")) +
scale_colour_manual(values=c("red", "blue", "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="darkgray") +
geom_hline(yintercept = 0,
linetype = "solid",
col = "darkgray",
size = 0.5) +
ylab("Movement Probability") +
xlab("Proportion Development (Dvlpd_A)") +
#ggtitle("Aerial Counts (COMPRECNTDEN)") +
theme_classic() +
scale_y_continuous(breaks = seq(0, 1, 0.1),
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),
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=18, vjust=0.5,
hjust=0.5, angle=0),
axis.title.x = element_text(face="bold", size = 20))
Proportion Hard Mast
range(Pnts$Hrdms_A)## [1] 0 100PlotFrame = Pnts@data %>% select(state, Hrdms_A, s_Hrdms_A)
PlotFrame$ExtRipM = plogis(PlotFrame$s_Hrdms_A*0.16876482)
PlotFrame$ExtRipL = plogis(PlotFrame$s_Hrdms_A*(0.16876482 - 0.07767469))
PlotFrame$ExtRipH = plogis(PlotFrame$s_Hrdms_A*(0.16876482 + 0.07767469))
PlotFrame$State = "State 1"
PlotFrame2 = Pnts@data %>% select(state, Hrdms_A, s_Hrdms_A)
PlotFrame2$ExtRipM = plogis(PlotFrame2$s_Hrdms_A* 0.06252680)
PlotFrame2$ExtRipL = plogis(PlotFrame2$s_Hrdms_A*(0.06252680 - 0.04703719))
PlotFrame2$ExtRipH = plogis(PlotFrame2$s_Hrdms_A*(0.06252680 + 0.04703719))
PlotFrame2$State = "State 2"
PlotFrame3 = Pnts@data %>% select(state, Hrdms_A, s_Hrdms_A)
PlotFrame3$ExtRipM = 1 - plogis(PlotFrame3$s_Hrdms_A*0.08445231)
PlotFrame3$ExtRipL = 1 - plogis(PlotFrame3$s_Hrdms_A*(0.08445231 - 0.03706168))
PlotFrame3$ExtRipH = 1 - plogis(PlotFrame3$s_Hrdms_A*(0.08445231 + 0.03706168))
PlotFrame3$State = "State 3"
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),
#method = "lm",
#span = 0.3,
#se = FALSE,
lwd = 1) +
scale_linetype_manual(values=c("solid", "longdash", "dotted")) +
scale_colour_manual(values=c("red", "blue", "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="darkgray") +
geom_hline(yintercept = 0,
linetype = "solid",
col = "darkgray",
size = 0.5) +
ylab("Movement Probability") +
xlab("Proportion Hard Mast (Hrdms_A)") +
#ggtitle("Aerial Counts (COMPRECNTDEN)") +
theme_classic() +
scale_y_continuous(breaks = seq(0, 1, 0.1),
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),
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=18, vjust=0.5,
hjust=0.5, angle=0),
axis.title.x = element_text(face="bold", size = 20))