Chicago Deer Models
Demographics of intensively managed urban deer
Contact: John (humph173@msu.edu)
Load Preprocessed Data
load("~/Michigan_State/Chicago/Run_Feb072019.RData")Model for Live Weight
#Model formula
LiveWeight.frm = Weight ~ -1 + intercept1 +
f(COMPRECNTDEN, #density
model="rw1",
constr = TRUE,
scale.model = TRUE,
hyper = nPrior) +
f(SexAgeGroup.f, #sex-age class
model="iid",
constr = TRUE,
hyper = nPrior) +
f(Month, #Month
model="iid",
constr = TRUE,
hyper = nPrior) +
f(Int.DATE_YR, #year
model="rw1",
constr = TRUE,
scale.model = TRUE,
hyper = nPrior ) +
f(GrowSeas,#grow season
model="rw1",
constr = TRUE,
scale.model = TRUE,
hyper = nPrior ) +
f(SumDur,#summer duration
model="rw1",
constr = TRUE,
scale.model = TRUE,
hyper = nPrior ) +
f(PREG, #pregnancy status
model="iid",
constr = TRUE,
hyper = nPrior) +
DSubZero
pcprior = list(prec = list(prior="pc.prec", param = c(1, 0.001)))
pcprior2 = list(prec = list(prior="pc.prec", param = c(3, 0.001)))
nPrior = list(theta=list(prior = "normal", param=c(0, 10)))
theta1 = c(-3.99299354, -0.55776308, -2.91058370, -1.14850504, -0.15903006, -0.04047979, 0.08543686, -0.02441974)
LiveWeight.mod = inla(LiveWeight.frm,
data = inla.stack.data(Sing.stk),
family = "gaussian",
verbose = TRUE,
control.fixed = list(prec = 0.001,
prec.intercept = 0.0001),
control.predictor = list(
A = inla.stack.A(Sing.stk),
compute = TRUE,
link = 1),
control.mode = list(restart = TRUE, theta = theta1),
control.inla = list(strategy="gaussian",
int.strategy = "eb"),
control.results = list(return.marginals.random = TRUE,
return.marginals.predictor = TRUE),
control.compute=list(dic = TRUE, cpo = TRUE, waic = TRUE))
summary(LiveWeight.mod)##
## Call:
## c("inla(formula = LiveWeight.frm, family = \"gaussian\", data = inla.stack.data(Sing.stk), ", " verbose = TRUE, control.compute = list(dic = TRUE, cpo = TRUE, ", " waic = TRUE), control.predictor = list(A = inla.stack.A(Sing.stk), ", " compute = TRUE, link = 1), control.inla = list(strategy = \"gaussian\", ", " int.strategy = \"eb\"), control.results = list(return.marginals.random = TRUE, ", " return.marginals.predictor = TRUE), control.fixed = list(prec = 0.001, ", " prec.intercept = 1e-04), control.mode = list(restart = TRUE, ", " theta = theta1))")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 2.2642 4.9627 0.4099 7.6368
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept1 55.7483 0.2236 55.3093 55.7483 56.1869 55.7483 0
## DSubZero -1.0227 0.3623 -1.7340 -1.0227 -0.3120 -1.0227 0
##
## Random effects:
## Name Model
## COMPRECNTDEN RW1 model
## SexAgeGroup.f IID model
## Month IID model
## Int.DATE_YR RW1 model
## GrowSeas RW1 model
## SumDur RW1 model
## PREG IID model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant
## Precision for the Gaussian observations 0.0184 0.0004 0.0176 0.0184
## Precision for COMPRECNTDEN 0.5831 0.1713 0.3116 0.5619
## Precision for SexAgeGroup.f 0.0553 0.0081 0.0411 0.0546
## Precision for Month 0.3432 0.0680 0.2226 0.3394
## Precision for Int.DATE_YR 0.8863 0.2441 0.5034 0.8542
## Precision for GrowSeas 1.0757 0.2993 0.5950 1.0411
## Precision for SumDur 1.1203 0.3000 0.6415 1.0835
## Precision for PREG 1.0049 0.2840 0.5567 0.9687
## 0.975quant mode
## Precision for the Gaussian observations 0.0193 0.0184
## Precision for COMPRECNTDEN 0.9798 0.5216
## Precision for SexAgeGroup.f 0.0731 0.0533
## Precision for Month 0.4881 0.3332
## Precision for Int.DATE_YR 1.4551 0.7935
## Precision for GrowSeas 1.7635 0.9751
## Precision for SumDur 1.8129 1.0137
## Precision for PREG 1.6639 0.9002
##
## Expected number of effective parameters(std dev): 27.15(0.00)
## Number of equivalent replicates : 140.02
##
## Deviance Information Criterion (DIC) ...............: 26007.37
## Deviance Information Criterion (DIC, saturated) ....: 3838.40
## Effective number of parameters .....................: 27.15
##
## Watanabe-Akaike information criterion (WAIC) ...: 26008.09
## Effective number of parameters .................: 27.65
##
## Marginal log-Likelihood: -13211.06
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computedLiveweight Curve
mic.d = as.data.frame(LiveWeight.mod$summary.random$COMPRECNTDEN[,c(1:6)])
names(mic.d) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975")
Raw.plt = ggplot(mic.d, aes(ID, Mean), lab) +
geom_point() +
geom_line(col = "black",
linetype= "solid",
lwd = 1) +
geom_ribbon(aes(ymin=Q025, ymax=Q975),
alpha=0.1,
linetype=1,
colour="transparent",
size=1,
fill="gray") +
geom_hline(yintercept = 0,
linetype = "solid",
col = "darkgray",
size = 0.5) +
#geom_point(data=FE.df1,
# aes(COMPREHARKM, LiveWtKg)) +
ylab("Liveweight (log)") +
xlab(expression(Deer~km^2)) +
ggtitle("Liveweight Vs Density (COMPRECNTDEN)") +
theme_classic() +
scale_y_continuous(breaks = c(-5, -4, -3, -1, -2, 0, 1, 2), limits = c(-5, 2)) +
scale_x_continuous(breaks = seq(4, 32, 8), limits =c(4,32)) +
theme(plot.title = element_text(hjust = 0.5),
title = element_text(face="bold", size=18, hjust=0.5),
axis.text=element_text(size=16),
strip.text = element_text(face="bold", size = 20),
axis.title.y = element_text(face="bold", size = 20),
axis.text.x = element_text(face="bold", size=14, vjust=0.5,
hjust=0.5, angle=0),
axis.title.x = element_text(face="bold", size = 20))
#Raw.plt
FindZero(mic.d$Mean, mic.d$ID)## [1] 9.25MaxX = arrange(mic.d, desc(Mean))[1,"ID"]
MaxX## [1] 8.3mic.d$Mean = mic.d$Mean + LiveWeight.mod$summary.fixed[1,1]
mic.d$Q025 = mic.d$Q025 + LiveWeight.mod$summary.fixed[1,1]
mic.d$Q975 = mic.d$Q975 + LiveWeight.mod$summary.fixed[1,1]
Raw.plt2 = ggplot(mic.d, aes(ID, Mean), lab) +
geom_smooth(col = "black",
linetype= "solid",
method = "loess",
span = 0.9,
se = FALSE,
lwd = 1) +
geom_smooth(data = mic.d, aes(ID, Q025),
col = "grey",
method = "loess",
span = 0.9,
se = FALSE,
linetype= "dashed") +
geom_smooth(data = mic.d, aes(ID, Q975),
col = "grey",
method = "loess",
span = 0.9,
se = FALSE,
linetype= "dashed") +
geom_hline(yintercept = LiveWeight.mod$summary.fixed[1,1],
linetype = "solid",
col = "darkgray",
size = 0.5) +
ylab("Liveweight (Kg)") +
xlab(expression(Deer~km^2)) +
#ggtitle("Aerial Counts (COMPRECNTDEN)") +
theme_classic() +
scale_y_continuous(breaks = seq(50, 58, 2),
limits = c(50, 58)) +
scale_x_continuous(breaks = seq(4, 32, 8), limits =c(4,32)) +
theme(plot.title = element_text(hjust = 0.5),
title = element_text(face="bold", size=18, hjust=0.5),
axis.text=element_text(size=16),
strip.text = element_text(face="bold", size = 20),
axis.title.y = element_text(face="bold", size = 20),
axis.text.x = element_text(face="bold", size=14, vjust=0.5,
hjust=0.5, angle=0),
axis.title.x = element_text(face="bold", size = 20))
#Raw.plt2
grid.arrange(Raw.plt, Raw.plt2, ncol=1)
Sex-Age Class Effect
Complx.plt = as.data.frame(LiveWeight.mod$summary.random$SexAgeGroup.f)
names(Complx.plt) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975", "mode", "kld")
ggplot(Complx.plt, aes(factor(ID), y=Mean)) +
geom_point(size=2, pch=19, col = "red") +
geom_linerange(aes(ymin=Q025, ymax=Q975), colour="black") +
geom_hline(yintercept = 0,
linetype = "dotted",
colour = "darkgray",
size = 0.5) +
theme_classic() +
xlab(" ") +
ylab("LiveWeight") +
ggtitle("Sex-Age Group Effect") +
scale_x_discrete(labels = c("Fawn (m)", "Yearling (m)", "Adult (m)",
"Fawn (f)", "Yearling (f)", "Adult (f)")) +
theme(plot.title = element_text(hjust = 0.5),
axis.title.y = element_text(face="bold", size=18),
axis.title.x = element_text(face="bold", size=18),
title = element_text(face="bold", size=18, hjust=0.5),
strip.text.x = element_text(face="bold", size = 14, colour = "black"),
axis.text.y = element_text(face="bold", size=14),
axis.text.x = element_text(face="bold", size=14, vjust=0.5,
hjust=0.3, angle=90))
Month Effect
Month of harvest.
Complx.plt = as.data.frame(LiveWeight.mod$summary.random$Month)
names(Complx.plt) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975", "mode", "kld")
ggplot(Complx.plt, aes(factor(ID), y=Mean)) +
geom_point(size=2, pch=19, col = "red") +
geom_linerange(aes(ymin=Q025, ymax=Q975), colour="black") +
geom_hline(yintercept = 0,
linetype = "dotted",
colour = "darkgray",
size = 0.5) +
theme_classic() +
xlab(" ") +
ylab("LiveWeight") +
ggtitle("Month Group Effect") +
scale_x_discrete(labels = c("Jan", "Feb", "March",
"Oct", "Nov", "Dec")) +
theme(plot.title = element_text(hjust = 0.5),
axis.title.y = element_text(face="bold", size=18),
axis.title.x = element_text(face="bold", size=18),
title = element_text(face="bold", size=18, hjust=0.5),
strip.text.x = element_text(face="bold", size = 14, colour = "black"),
axis.text.y = element_text(face="bold", size=14),
axis.text.x = element_text(face="bold", size=14, vjust=0.5,
hjust=0.3, angle=90))
Year Effect
mic.d = as.data.frame(LiveWeight.mod$summary.random$Int.DATE_YR[,c(1:6)])
names(mic.d) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975")
mic.d$YLab = paste(levels(factor(FE.df1$DATE_YR)))
YLab = paste(mic.d$YLab)
ggplot(mic.d, aes(ID, Mean), lab) +
geom_smooth(col = "black",
linetype= "solid",
method = "loess",
span = 0.9,
se = FALSE,
lwd = 1) +
geom_smooth(data = mic.d, aes(ID, Q025),
col = "grey",
method = "loess",
span = 0.9,
se = FALSE,
linetype= "dashed") +
geom_smooth(data = mic.d, aes(ID, Q975),
col = "grey",
method = "loess",
span = 0.9,
se = FALSE,
linetype= "dashed") +
geom_hline(yintercept = 0,
linetype = "solid",
col = "darkgray",
size = 0.5) +
ylab("LiveWeight") +
xlab(" ") +
scale_x_continuous(breaks = 4:20,labels=paste0(mic.d$YLab)) +
ggtitle("Year Temporal Effect") +
theme_classic() +
theme(plot.title = element_text(hjust = 0.5),
title = element_text(face="bold", size=18, hjust=0.5),
axis.text=element_text(size=16),
strip.text = element_text(face="bold", size = 20),
axis.title.y = element_text(face="bold", size = 20),
axis.text.x = element_text(face="bold", size=14, vjust=0.5, angle=90),
axis.title.x = element_text(face="bold", size = 20))
Doe Pregnancy
Month of harvest.
Complx.plt = as.data.frame(LiveWeight.mod$summary.random$PREG)
names(Complx.plt) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975", "mode", "kld")
ggplot(Complx.plt, aes(factor(ID), y=Mean)) +
geom_point(size=2, pch=19, col = "red") +
geom_linerange(aes(ymin=Q025, ymax=Q975), colour="black") +
geom_hline(yintercept = 0,
linetype = "dotted",
colour = "darkgray",
size = 0.5) +
theme_classic() +
xlab(" ") +
ylab("LiveWeight") +
ggtitle("Doe Pregnancy Effect") +
scale_x_discrete(labels = c("Yes", "No", "March")) +
theme(plot.title = element_text(hjust = 0.5),
axis.title.y = element_text(face="bold", size=18),
axis.title.x = element_text(face="bold", size=18),
title = element_text(face="bold", size=18, hjust=0.5),
strip.text.x = element_text(face="bold", size = 14, colour = "black"),
axis.text.y = element_text(face="bold", size=14),
axis.text.x = element_text(face="bold", size=14, vjust=0.5))
Growth Season Effect
mic.d = as.data.frame(LiveWeight.mod$summary.random$GrowSeas[,c(1:6)])
names(mic.d) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975")
mic.d$IDsc = with(Climate.df,
GrowSeas[match(
mic.d$ID,
rLs_GrowSeas)])
ggplot(mic.d, aes(IDsc, Mean), lab) +
geom_smooth(col = "black",
linetype= "solid",
method = "loess",
span = 0.9,
se = FALSE,
lwd = 1) +
geom_smooth(data = mic.d, aes(IDsc, Q025),
col = "grey",
method = "loess",
span = 0.9,
se = FALSE,
linetype= "dashed") +
geom_smooth(data = mic.d, aes(IDsc, Q975),
col = "grey",
method = "loess",
span = 0.9,
se = FALSE,
linetype= "dashed") +
geom_hline(yintercept = 0,
linetype = "solid",
col = "darkgray",
size = 0.5) +
ylab("LiveWeight") +
xlab("Length of Growth Season (Days) ") +
ggtitle("Growth Season Effect") +
theme_classic() +
scale_x_continuous(breaks = seq(0, 265, 20), limits =c(150,265)) +
theme(plot.title = element_text(hjust = 0.5),
title = element_text(face="bold", size=18, hjust=0.5),
axis.text=element_text(size=16),
strip.text = element_text(face="bold", size = 20),
axis.title.y = element_text(face="bold", size = 20),
axis.text.x = element_text(face="bold", size=14, vjust=0.5,
hjust=0.5, angle=0),
axis.title.x = element_text(face="bold", size = 20))
Body Fat (Kistner - TOTAL)
BFat.frm = Kistner ~ -1 + intercept1 +
f(COMPRECNTDEN,
model="rw1",
constr = TRUE,
scale.model = TRUE,
hyper = nPrior) +
f(SexAgeGroup.f,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(Month,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(Int.DATE_YR,
model="rw1",
constr = TRUE,
scale.model = TRUE,
hyper = nPrior ) +
f(GrowSeas,
model="rw1",
constr = TRUE,
scale.model = TRUE,
hyper = nPrior ) +
f(SumDur,
model="rw1",
constr = TRUE,
scale.model = TRUE,
hyper = nPrior ) +
f(PREG,
model="iid",
constr = TRUE,
hyper = nPrior) +
DSubZero
pcprior = list(prec = list(prior="pc.prec", param = c(1, 0.001)))
pcprior2 = list(prec = list(prior="pc.prec", param = c(3, 0.001)))
nPrior = list(theta=list(prior = "normal", param=c(0, 10)))
theta1 = c(-3.99299354, -0.55776308, -2.91058370, -1.14850504, -0.15903006, -0.04047979, 0.08543686, -0.02441974)
BodyFat.mod = inla(BFat.frm,
data = inla.stack.data(SingF.stk),
family = "gaussian",
verbose = TRUE,
control.fixed = list(prec = 0.001,
prec.intercept = 0.0001),
control.predictor = list(
A = inla.stack.A(SingF.stk),
compute = TRUE,
link = 1),
control.mode = list(restart = TRUE, theta = theta1),
control.inla = list(strategy="gaussian",
int.strategy = "eb"),
control.results = list(return.marginals.random = TRUE,
return.marginals.predictor = TRUE),
control.compute=list(dic = TRUE, cpo = TRUE, waic = TRUE))
summary(BodyFat.mod)##
## Call:
## c("inla(formula = BFat.frm, family = \"gaussian\", data = inla.stack.data(SingF.stk), ", " verbose = TRUE, control.compute = list(dic = TRUE, cpo = TRUE, ", " waic = TRUE), control.predictor = list(A = inla.stack.A(SingF.stk), ", " compute = TRUE, link = 1), control.inla = list(strategy = \"gaussian\", ", " int.strategy = \"eb\"), control.results = list(return.marginals.random = TRUE, ", " return.marginals.predictor = TRUE), control.fixed = list(prec = 0.001, ", " prec.intercept = 1e-04), control.mode = list(restart = TRUE, ", " theta = theta1))")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 2.1300 9.1466 0.2892 11.5658
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept1 66.5913 0.7818 65.0563 66.5913 68.1251 66.5913 0
## DSubZero -0.7819 0.5621 -1.8854 -0.7819 0.3207 -0.7819 0
##
## Random effects:
## Name Model
## COMPRECNTDEN RW1 model
## SexAgeGroup.f IID model
## Month IID model
## Int.DATE_YR RW1 model
## GrowSeas RW1 model
## SumDur RW1 model
## PREG IID model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant
## Precision for the Gaussian observations 0.0044 0.0001 0.0042 0.0044
## Precision for COMPRECNTDEN 0.3155 0.1269 0.1350 0.2933
## Precision for SexAgeGroup.f 0.1350 0.0289 0.0926 0.1299
## Precision for Month 0.0675 0.0142 0.0402 0.0678
## Precision for Int.DATE_YR 0.3064 0.1282 0.1293 0.2822
## Precision for GrowSeas 0.7812 0.2499 0.4014 0.7446
## Precision for SumDur 0.7003 0.2493 0.3215 0.6648
## Precision for PREG 0.6617 0.1786 0.3917 0.6337
## 0.975quant mode
## Precision for the Gaussian observations 0.0046 0.0044
## Precision for COMPRECNTDEN 0.6252 0.2530
## Precision for SexAgeGroup.f 0.2046 0.1190
## Precision for Month 0.0944 0.0693
## Precision for Int.DATE_YR 0.6245 0.2400
## Precision for GrowSeas 1.3723 0.6764
## Precision for SumDur 1.2902 0.5959
## Precision for PREG 1.0884 0.5805
##
## Expected number of effective parameters(std dev): 22.36(0.00)
## Number of equivalent replicates : 132.76
##
## Deviance Information Criterion (DIC) ...............: 24450.62
## Deviance Information Criterion (DIC, saturated) ....: 2868.54
## Effective number of parameters .....................: 22.36
##
## Watanabe-Akaike information criterion (WAIC) ...: 24450.47
## Effective number of parameters .................: 22.03
##
## Marginal log-Likelihood: -12479.81
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computedBody Fat Curve
mic.d = as.data.frame(BodyFat.mod$summary.random$COMPRECNTDEN[,c(1:6)])
names(mic.d) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975")
Raw.plt = ggplot(mic.d, aes(ID, Mean), lab) +
geom_point() +
geom_line(col = "black",
linetype= "solid",
lwd = 1) +
geom_ribbon(aes(ymin=Q025, ymax=Q975),
alpha=0.1,
linetype=1,
colour="transparent",
size=1,
fill="gray") +
geom_hline(yintercept = 0,
linetype = "solid",
col = "darkgray",
size = 0.5) +
#geom_point(data=FE.df1,
# aes(COMPREHARKM, LiveWtKg)) +
ylab("Kistner TOTAL (log)") +
xlab(expression(Deer~km^2)) +
ggtitle("Kistner Vs Density (COMPRECNTDEN)") +
theme_classic() +
scale_y_continuous(breaks = c(-6, -5, -4, -3, -1, -2, 0, 1, 2, 3, 4), limits = c(-6, 4)) +
scale_x_continuous(breaks = seq(4, 32, 8), limits =c(4,32)) +
theme(plot.title = element_text(hjust = 0.5),
title = element_text(face="bold", size=18, hjust=0.5),
axis.text=element_text(size=16),
strip.text = element_text(face="bold", size = 20),
axis.title.y = element_text(face="bold", size = 20),
axis.text.x = element_text(face="bold", size=14, vjust=0.5,
hjust=0.5, angle=0),
axis.title.x = element_text(face="bold", size = 20))
#Raw.plt
FindZero(mic.d$Mean, mic.d$ID)## [1] 7.8 22.0MaxX = arrange(mic.d, desc(Mean))[1,"ID"]
MaxX## [1] 8.3mic.d$Mean = mic.d$Mean + BodyFat.mod$summary.fixed[1,1]
mic.d$Q025 = mic.d$Q025 + BodyFat.mod$summary.fixed[1,1]
mic.d$Q975 = mic.d$Q975 + BodyFat.mod$summary.fixed[1,1]
Raw.plt2 = ggplot(mic.d, aes(ID, Mean), lab) +
geom_smooth(col = "black",
linetype= "solid",
method = "loess",
span = 0.9,
se = FALSE,
lwd = 1) +
geom_smooth(data = mic.d, aes(ID, Q025),
col = "grey",
method = "loess",
span = 0.9,
se = FALSE,
linetype= "dashed") +
geom_smooth(data = mic.d, aes(ID, Q975),
col = "grey",
method = "loess",
span = 0.9,
se = FALSE,
linetype= "dashed") +
geom_hline(yintercept = BodyFat.mod$summary.fixed[1,1],
linetype = "solid",
col = "darkgray",
size = 0.5) +
ylab("Kistner (TOTAL)") +
xlab(expression(Deer~km^2)) +
#ggtitle("Aerial Counts (COMPRECNTDEN)") +
theme_classic() +
scale_y_continuous(breaks = seq(60, 71, 2),
limits = c(60, 71)) +
scale_x_continuous(breaks = seq(4, 32, 8), limits =c(4,32)) +
theme(plot.title = element_text(hjust = 0.5),
title = element_text(face="bold", size=18, hjust=0.5),
axis.text=element_text(size=16),
strip.text = element_text(face="bold", size = 20),
axis.title.y = element_text(face="bold", size = 20),
axis.text.x = element_text(face="bold", size=14, vjust=0.5,
hjust=0.5, angle=0),
axis.title.x = element_text(face="bold", size = 20))
#Raw.plt2
grid.arrange(Raw.plt, Raw.plt2, ncol=1)
## [1] 2726## [1] 116## [1] 1848Fecundity Model
Fecun.frm = FETUS ~ -1 + intercept1 +
f(COMPRECNTDEN,
model="rw1",
constr = TRUE,
scale.model = TRUE,
hyper = nPrior) +
f(SexAgeGroup.f,
model="iid",
constr = TRUE,
hyper = nPrior)
theta1 = c(0.49173512, -0.06059399)
Fecund.mod = inla(Fecun.frm,
data = inla.stack.data(FETUS1.stk),
family = "poisson",
verbose = TRUE,
control.fixed = list(prec = 0.001,
prec.intercept = 0.0001),
control.predictor = list(
A = inla.stack.A(FETUS1.stk),
compute = TRUE,
link = 1),
control.mode = list(restart = TRUE, theta = theta1),
control.inla = list(strategy="gaussian",
int.strategy = "eb"),
control.results = list(return.marginals.random = TRUE,
return.marginals.predictor = TRUE),
control.compute=list(dic = TRUE, cpo = TRUE, waic = TRUE))
summary(Fecund.mod)##
## Call:
## c("inla(formula = Fecun.frm, family = \"poisson\", data = inla.stack.data(FETUS1.stk), ", " verbose = TRUE, control.compute = list(dic = TRUE, cpo = TRUE, ", " waic = TRUE), control.predictor = list(A = inla.stack.A(FETUS1.stk), ", " compute = TRUE, link = 1), control.inla = list(strategy = \"gaussian\", ", " int.strategy = \"eb\"), control.results = list(return.marginals.random = TRUE, ", " return.marginals.predictor = TRUE), control.fixed = list(prec = 0.001, ", " prec.intercept = 1e-04), control.mode = list(restart = TRUE, ", " theta = theta1))")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 1.2617 2.0775 0.1935 3.5327
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept1 -0.174 0.061 -0.2937 -0.174 -0.0543 -0.174 0
##
## Random effects:
## Name Model
## COMPRECNTDEN RW1 model
## SexAgeGroup.f IID model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant 0.975quant
## Precision for COMPRECNTDEN 1.711 0.5249 0.9066 1.6357 2.951
## Precision for SexAgeGroup.f 0.984 0.2954 0.5259 0.9435 1.678
## mode
## Precision for COMPRECNTDEN 1.4956
## Precision for SexAgeGroup.f 0.8675
##
## Expected number of effective parameters(std dev): 13.87(0.00)
## Number of equivalent replicates : 75.97
##
## Deviance Information Criterion (DIC) ...............: 2448.21
## Deviance Information Criterion (DIC, saturated) ....: 344.31
## Effective number of parameters .....................: 13.92
##
## Watanabe-Akaike information criterion (WAIC) ...: 2437.80
## Effective number of parameters .................: 3.46
##
## Marginal log-Likelihood: -1257.15
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computedFecundity Curve
mic.d = as.data.frame(Fecund.mod$summary.random$COMPRECNTDEN[,c(1:6)])
names(mic.d) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975")
Raw.plt = ggplot(mic.d, aes(ID, Mean), lab) +
geom_point() +
geom_line(col = "black",
linetype= "solid",
lwd = 1) +
geom_ribbon(aes(ymin=Q025, ymax=Q975),
alpha=0.1,
linetype=1,
colour="transparent",
size=1,
fill="gray") +
geom_hline(yintercept = 0,
linetype = "solid",
col = "darkgray",
size = 0.5) +
ylab("Fetus Count (log)") +
xlab(expression(Deer~km^2)) +
ggtitle("Fecundity Vs Density (COMPRECNTDEN)") +
theme_classic() +
scale_y_continuous(breaks = c(-0.3, -0.2, -0.1, 0, 0.1, 0.2, 0.3), limits = c(-0.35, 0.36)) +
scale_x_continuous(breaks = seq(4, 32, 8), limits =c(4,32)) +
theme(plot.title = element_text(hjust = 0.5),
title = element_text(face="bold", size=18, hjust=0.5),
axis.text=element_text(size=16),
strip.text = element_text(face="bold", size = 20),
axis.title.y = element_text(face="bold", size = 20),
axis.text.x = element_text(face="bold", size=14, vjust=0.5,
hjust=0.5, angle=0),
axis.title.x = element_text(face="bold", size = 20))
#Raw.plt
FindZero(mic.d$Mean, mic.d$ID)## [1] 11.05MaxX = arrange(mic.d, desc(Mean))[1,"ID"]
MaxX## [1] 4.3mic.d$Mean = exp(mic.d$Mean)
mic.d$Q025 = exp(mic.d$Q025)
mic.d$Q975 = exp(mic.d$Q975)
Raw.plt2 = ggplot(mic.d, aes(ID, Mean), lab) +
geom_smooth(col = "black",
linetype= "solid",
method = "loess",
span = 0.9,
se = FALSE,
lwd = 1) +
geom_smooth(data = mic.d, aes(ID, Q025),
col = "grey",
method = "loess",
span = 0.9,
se = FALSE,
linetype= "dashed") +
geom_smooth(data = mic.d, aes(ID, Q975),
col = "grey",
method = "loess",
span = 0.9,
se = FALSE,
linetype= "dashed") +
geom_hline(yintercept = 1,
linetype = "solid",
col = "darkgray",
size = 0.5) +
ylab("Fetus Count") +
xlab(expression(Deer~km^2)) +
#ggtitle("Aerial Counts (COMPRECNTDEN)") +
theme_classic() +
scale_y_continuous(breaks = seq(0.5, 1.5, 0.25),
limits = c(0.5, 1.5)) +
scale_x_continuous(breaks = seq(4, 32, 8), limits =c(4,32)) +
theme(plot.title = element_text(hjust = 0.5),
title = element_text(face="bold", size=18, hjust=0.5),
axis.text=element_text(size=16),
strip.text = element_text(face="bold", size = 20),
axis.title.y = element_text(face="bold", size = 20),
axis.text.x = element_text(face="bold", size=14, vjust=0.5,
hjust=0.5, angle=0),
axis.title.x = element_text(face="bold", size = 20))
#Raw.plt2
grid.arrange(Raw.plt, Raw.plt2, ncol=1)