Deer Density Dependence
1 Data Preprocessing
date() ## [1] "Thu May 16 10:14:36 2019"1.1 Load Data
Load LiveWeightOutput Data
Deer.Models.Data = read.csv("LiveWeightOutput_110518.csv", sep = ",", header = TRUE)Convert LiveWtKg
Deer.Models.Data$LiveWtKg = Deer.Models.Data$EstLiveWWT * 0.453592Get state boundaries (for figures)
LL84 = "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
States = map("state",
fill = TRUE,
plot = FALSE)
IDs = sapply(strsplit(States$names, ":"),
function(x) x[1])
States = map2SpatialPolygons(States, IDs = IDs,
proj4string = CRS(LL84))
#Add a dataframe
pid = sapply(slot(States, "polygons"),
function(x) slot(x, "ID"))
p.df = data.frame(ID=1:length(States),
row.names = pid)
States = SpatialPolygonsDataFrame(States, p.df)
Ill = subset(States, ID == 12 | ID == 48 | ID == 13 | ID == 21)1.2 Climate Analysis
Locate Stations.
Stations = isd_stations_search(lat = 41.7508, lon = -88.1535,
radius = 100, bbox = NULL) #100km radius of Naperville, Dupage County, Illinois
Station.pnts = as.data.frame(Stations)
Station.pnts = SpatialPointsDataFrame(Station.pnts[,c("lon", "lat")], Station.pnts)
proj4string(Station.pnts) = proj4string(Ill)View Weather Stations
USs_df = fortify(Ill)## Regions defined for each PolygonsBr.plt = ggplot(USs_df,
aes(long,lat, group=group)) +
geom_polygon(fill="tan", col="black") +
geom_point(data=Station.pnts@data,
aes(lon, lat,
group=NULL),
col = "red", size=1, shape=2) +
xlab("Longitude") +
ylab("Latitude") +
coord_equal() +
ggtitle("Weather Stations") +
theme(panel.grid.minor = element_blank(),
panel.grid.major = element_blank(),
panel.background = element_blank(),
plot.background = element_blank(),
panel.border = element_blank(),
legend.title = element_text(size=12, face="bold"),
axis.title.x = element_text(size=16, face="bold"),
axis.title.y = element_text(size=16, face="bold"),
plot.title = element_text(size=22, face="bold", hjust = 0.5))
Br.plt
Download climate data by station.
plot(Ill, main="Weather Stations")
plot(Station.pnts, add=T, col="red")
Naperville = data.frame(id = "Naperville", latitude = 41.7508, longitude = -88.1535)
station_data = ghcnd_stations()
Stat.df = meteo_nearby_stations(lat_lon_df = Naperville, station_data = station_data,
radius = 250, limit=250, var = c("PRCP", "TMAX"),
year_min = 1990, year_max = 2010)
Clim.df = meteo_pull_monitors(Stat.df$Naperville$id,
date_min = "1990-01-01",
date_max = "2010-12-31")
Clim.df = as.data.frame(Clim.df)
#write.csv(Clim.df, "./Climate_102118.csv") #Save copy
#token "VmLVwVYuKJixCseSiebNvcoObqeMirtI" (needed for database acsess)Climate Analysis
ExtClim.df = read.csv("./Climate_102118.csv", stringsAsFactors=FALSE) %>% #preprocessed above
mutate(date = as.Date(date, "%Y-%m-%d"),
Date = as.PCICt(as.character(date),
cal = "gregorian",
format = "%Y-%m-%d"),
Year = year(date),
Month = month(date),
Day = day(date),
DOY = yday(date),
precip = prcp/10, #to mm
t_max = tmax/10, #to C
t_min = tmin/10) %>% #to C %>%
select(id, date, Date, Year, Month, Day, DOY, precip, snow, snwd, t_max, t_min)
ExtClim.df[is.na(ExtClim.df)] = 0
ExtClim.df$t_avg = (ExtClim.df$t_max + ExtClim.df$t_min)/2
Temp.df = as.data.frame(
ExtClim.df %>%
group_by(id) %>%
summarize(Tmax = mean(t_max, na.rm=T),
Tmin = mean(t_min, na.rm=T)))
loc.lvs = levels(factor(unique(ExtClim.df$id)))
for(i in 1:length(loc.lvs)){
tmp.f = ExtClim.df %>% filter(id == loc.lvs[i])
tmp.f$Tmx.chk = mean(tmp.f$t_max, na.rm=T)
tmp.f$Tmn.chk = mean(tmp.f$t_min, na.rm=T)
tmp.f$Filt = tmp.f$Tmx.chk+tmp.f$Tmn.chk
if(i == 1){Tmp.chk = tmp.f}
else{Tmp.chk = rbind(Tmp.chk, tmp.f)}
}
ExtClim.df = Tmp.chk %>% filter(Filt > 0) %>% select(-c(Tmx.chk, Tmn.chk, Filt))
Yr.lvs = levels(factor(ExtClim.df$Year))
loc.lvs = levels(factor(unique(ExtClim.df$id)))
for(i in 1:length(loc.lvs)){
snw2.tmp = ExtClim.df %>% filter(id == loc.lvs[i])
Srt.yr = as.numeric(min(levels(factor(snw2.tmp$Year))))
End.yr = as.numeric(max(levels(factor(snw2.tmp$Year))))
#ClimDex Extremes
ci = climdexInput.raw(snw2.tmp$t_max,
snw2.tmp$t_min,
snw2.tmp$precip,
snw2.tmp$Date,
snw2.tmp$Date,
snw2.tmp$Date,
base.range=c(Srt.yr, End.yr))
Frm1.df = as.data.frame(
cbind(
GrowSeas = climdex.gsl(ci, gsl.mode = "GSL"), #Grow season length
PrecEx10mm = climdex.r10mm(ci), #precip exceeding 10mm per day
PrecEx20mm = climdex.r20mm(ci), #precip exceeding 20mm per day
DSubZero = climdex.id(ci), #count of days maximum temp below 0
Prec95pct = climdex.r95ptot(ci), #precip over 95th %
MonMxDT = climdex.txx(ci, freq = "annual"), #Annual Max Daily temp
MMx1DPrec = climdex.rx1day(ci, freq = "annual"), #Ann Maximum 1-day Precipitation
AMinT = climdex.tnn(ci, freq = "annual"), #Annual Minimum Of Daily Minimum Temperature
AMxT = climdex.txn(ci, freq = "annual"), #Annual Minimum Of Daily Maximum Temperature
SumDur = climdex.su(ci), #Summer duration by year
FrstDays = climdex.fd(ci))) #Frost days
Frm1.df$Year = row.names(Frm1.df)
Frm1.df$Station = paste(loc.lvs[i])
if(i == 1){ClimChn = Frm1.df}
else{ClimChn = rbind(ClimChn, Frm1.df)}
}Match and Organize Climate Info by Year
length(ClimChn$Station) #Number of stations## [1] 279Climate.df = as.data.frame(
ClimChn %>%
group_by(Year) %>%
summarize(GrowSeas = mean(GrowSeas, na.rm=T),
PrecEx10mm = mean(PrecEx10mm, na.rm=T),
PrecEx20mm = mean(PrecEx20mm , na.rm=T),
DSubZero = mean(DSubZero, na.rm=T),
Prec95pct = mean(Prec95pct, na.rm=T),
MonMxDT = mean(MonMxDT, na.rm=T),
MMx1DPrec = mean(MMx1DPrec, na.rm=T),
AMinT = mean(AMinT, na.rm=T),
AMxT = mean(AMxT, na.rm=T),
SumDur = mean(SumDur, na.rm=T),
FrstDays = mean(FrstDays, na.rm=T)))
###Scale and Center
Cov.scale = Climate.df[,2:12]
for(i in 1:ncol(Cov.scale)){
Cov.scale[,i] = as.numeric(scale(Cov.scale[,i], scale = T, center = T))
}
names(Cov.scale) = paste("s_", names(Climate.df[,2:12]), sep="")
#Lag 1 year
lag.scale = Cov.scale
for(i in 1:ncol(lag.scale)){
lag.scale[,i] = lag(lag.scale[,i], 1)
lag.scale[,i][is.na(lag.scale[,i] )] = mean(lag.scale[,i], na.rm=T)
}
names(lag.scale) = paste("L", names(Cov.scale), sep="")
Climate.df = cbind(Climate.df, Cov.scale, lag.scale)
##Check for colinearity
corrplot(cor(Cov.scale))
library(perturb)
CI = colldiag(cor(Cov.scale))
CI## Condition
## Index Variance Decomposition Proportions
## intercept s_GrowSeas s_PrecEx10mm s_PrecEx20mm s_DSubZero
## 1 1.000 0.000 0.001 0.000 0.000 0.000
## 2 1.206 0.000 0.000 0.000 0.000 0.000
## 3 1.252 0.003 0.001 0.001 0.001 0.000
## 4 1.791 0.004 0.001 0.000 0.000 0.000
## 5 3.508 0.002 0.001 0.001 0.000 0.007
## 6 6.252 0.000 0.000 0.001 0.000 0.002
## 7 7.653 0.087 0.083 0.000 0.011 0.002
## 8 12.894 0.229 0.095 0.026 0.003 0.001
## 9 18.347 0.171 0.291 0.013 0.001 0.000
## 10 38.715 0.005 0.338 0.941 0.487 0.040
## 11 53.072 0.499 0.189 0.018 0.496 0.948
## s_Prec95pct s_MonMxDT s_MMx1DPrec s_AMinT s_AMxT s_SumDur s_FrstDays
## 1 0.000 0.000 0.001 0.000 0.000 0.001 0.000
## 2 0.000 0.000 0.001 0.002 0.001 0.000 0.000
## 3 0.000 0.000 0.000 0.000 0.000 0.000 0.000
## 4 0.000 0.000 0.007 0.002 0.000 0.001 0.001
## 5 0.004 0.002 0.013 0.006 0.000 0.003 0.000
## 6 0.001 0.003 0.002 0.001 0.010 0.089 0.020
## 7 0.007 0.006 0.038 0.003 0.002 0.023 0.007
## 8 0.001 0.038 0.118 0.008 0.065 0.041 0.027
## 9 0.065 0.000 0.001 0.685 0.090 0.036 0.019
## 10 0.150 0.216 0.413 0.247 0.301 0.027 0.066
## 11 0.771 0.735 0.407 0.047 0.530 0.779 0.860Cov.scale2 = Cov.scale %>% select(-c(s_Prec95pct, s_MonMxDT))
CI = colldiag(cor(Cov.scale2))
CI## Condition
## Index Variance Decomposition Proportions
## intercept s_GrowSeas s_PrecEx10mm s_PrecEx20mm s_DSubZero
## 1 1.000 0.000 0.001 0.000 0.000 0.005
## 2 1.158 0.003 0.006 0.001 0.001 0.000
## 3 1.279 0.000 0.000 0.002 0.005 0.000
## 4 1.598 0.010 0.005 0.000 0.000 0.002
## 5 3.940 0.001 0.033 0.000 0.001 0.121
## 6 5.634 0.000 0.048 0.004 0.002 0.035
## 7 7.429 0.088 0.199 0.004 0.076 0.025
## 8 18.905 0.066 0.002 0.374 0.632 0.810
## 9 22.542 0.832 0.706 0.615 0.283 0.000
## s_MMx1DPrec s_AMinT s_AMxT s_SumDur s_FrstDays
## 1 0.000 0.001 0.001 0.003 0.003
## 2 0.005 0.000 0.001 0.007 0.000
## 3 0.001 0.001 0.001 0.014 0.000
## 4 0.019 0.001 0.000 0.000 0.002
## 5 0.083 0.011 0.001 0.063 0.001
## 6 0.003 0.001 0.036 0.529 0.116
## 7 0.004 0.006 0.030 0.375 0.163
## 8 0.144 0.292 0.719 0.002 0.179
## 9 0.743 0.688 0.210 0.007 0.535###Match to deer data
Deer.Models.Data$GrowSeas = with(Climate.df,
Ls_GrowSeas[match(
Deer.Models.Data$DATE_YR,
Year)])
Deer.Models.Data$PrecEx10mm = with(Climate.df,
Ls_PrecEx10mm[match(
Deer.Models.Data$DATE_YR,
Year)])
Deer.Models.Data$PrecEx20mm = with(Climate.df,
Ls_PrecEx20mm[match(
Deer.Models.Data$DATE_YR,
Year)])
Deer.Models.Data$DSubZero = with(Climate.df,
Ls_DSubZero[match(
Deer.Models.Data$DATE_YR,
Year)])
Deer.Models.Data$Prec95pct = with(Climate.df,
Ls_Prec95pct[match(
Deer.Models.Data$DATE_YR,
Year)])
Deer.Models.Data$MonMxDT = with(Climate.df,
Ls_MonMxDT[match(
Deer.Models.Data$DATE_YR,
Year)])
Deer.Models.Data$MMx1DPrec = with(Climate.df,
Ls_MMx1DPrec[match(
Deer.Models.Data$DATE_YR,
Year)])
Deer.Models.Data$AMinT = with(Climate.df,
Ls_AMinT[match(
Deer.Models.Data$DATE_YR,
Year)])
Deer.Models.Data$AMxT = with(Climate.df,
Ls_AMxT[match(
Deer.Models.Data$DATE_YR,
Year)])
Deer.Models.Data$SumDur = with(Climate.df,
Ls_SumDur[match(
Deer.Models.Data$DATE_YR,
Year)])
Deer.Models.Data$FrstDays = with(Climate.df,
Ls_FrstDays[match(
Deer.Models.Data$DATE_YR,
Year)])
Deer.Models.Data$SumDur.sc = Scale(Deer.Models.Data$SumDur)
Deer.Models.Data$GrowSeas.sc = Scale(Deer.Models.Data$GrowSeas)1.3 Age Groups and Covariates
Age and Sex Groups
11 = male fawn, 12 = male yearling, 13 = male adult
21 = female fawn, 22 = female yearling, 23 = female adult
Deer.Models.Data$SexAgeGeoup = ((Deer.Models.Data$SEX * 10) + Deer.Models.Data$AGEGROUP)
Deer.Models.Data$SexAgeGroup.f = factor(Deer.Models.Data$SexAgeGeoup)
is.factor(Deer.Models.Data$SexAgeGroup.f)## [1] TRUEYears, Complex, and Data Correction
Deer.Models.Data$Int.DATE_YR = as.integer(factor(Deer.Models.Data$RMVYEAR1))
Deer.Models.Data$f_COMPLEX = Deer.Models.Data$COMPLEX
#Dwayne Email: IDNRTAG: A063321, Change COMPRECNTDEN from 11.5 to 11.8
Deer.Models.Data$COMPRECNTDEN = ifelse(Deer.Models.Data$IDNRTAG == "A063321", 11.8, Deer.Models.Data$COMPRECNTDEN)Match Sample Location Info to Data Frame
Goog.df = read.csv("./LatLongArea_042219.csv", stringsAsFactors=FALSE)
Lat_tab = as.data.frame(str_split(Goog.df$Location, ", ", simplify = TRUE))
names(Lat_tab) = c("Lat", "Long")
Lat_tab$Lat = as.numeric(as.character(Lat_tab$Lat))
Lat_tab$Long = as.numeric(as.character(Lat_tab$Long))
Goog.df = as.data.frame(cbind(Goog.df, Lat_tab))
Deer.Models.Data$Lat = with(Goog.df,
Lat[match(
Deer.Models.Data$COMPLEX,
Complex_Code)])
Deer.Models.Data$Long = with(Goog.df,
Long[match(
Deer.Models.Data$COMPLEX,
Complex_Code)])
Deer.Models.Data$CmpxArea = with(Goog.df,
ComAreaKm2[match(
Deer.Models.Data$COMPLEX,
Complex_Code)])Nearest Neighbor Distance
Between Complexes
LL84 = "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
Complx.pnts = Deer.Models.Data %>% select(COMPLEX, Lat, Long, CmpxArea)
Complx.pnts$Dups = duplicated(Complx.pnts[,c("COMPLEX", "Lat", "Long")])
Complx.pnts = Complx.pnts %>% filter(Dups == FALSE)
Complx.pnts = SpatialPointsDataFrame(Complx.pnts[,c("Long","Lat")], Complx.pnts)
proj4string(Complx.pnts) = LL84
suppressMessages(library(spatstat))
nProj = "+proj=utm +zone=16 +datum=NAD83 +units=km +no_defs +ellps=GRS80 +towgs84=0,0,0"
P.pp = spTransform(Complx.pnts, nProj)
P.pp = as.ppp(P.pp)
Complx.pnts$NN = (nndist(P.pp, k = 1))
range(Complx.pnts$NN) ## [1] 3.976131 26.198469Deer.Models.Data$NN = with(Complx.pnts@data,
NN[match(
Deer.Models.Data$COMPLEX,
COMPLEX)])
range(Deer.Models.Data$NN) ## [1] 3.976131 26.1984692 Deer Models (Complex 1)
2.1 Whole Weight
Organize Data as Lists Object
FE.df1 =Deer.Models.Data %>% filter(COMPLEX == 1) #Complex 1 only
#COMCOMCNT = "Survey completed" If not, set record to NA
FE.df1$COMPRECNTDEN2 = ifelse(FE.df1$COMCOMCNT == 2, NA, FE.df1$COMPRECNTDEN) #Complete aerial survey
FE.df1 = FE.df1 %>%
filter(is.na(AGEGROUP) == FALSE) #Remove records w/ missing age
#Integer version of Sex-Age groups
unique(FE.df1$SexAgeGroup.f)## [1] 11 12 13 21 22 23
## Levels: 11 12 13 21 22 23FE.df1$RepGroups = as.integer(FE.df1$SexAgeGroup.f)
unique(FE.df1$RepGroups)## [1] 1 2 3 4 5 6#Unique ID by Deer
FE.df1$myID = 1:nrow(FE.df1)
#set up model data as a list object
## Not all covariates are used.
FE.lst1 = list(list(intercept1 = rep(1, dim(FE.df1)[1])),
list(COMPREHARKM = FE.df1[,"COMPREHARKM"],
COMRECPOP = FE.df1[,"COMRECPOP"],
COMRECDOEPOP = FE.df1[,"COMRECDOEPOP"],
COMCUMPRE2HARKM = FE.df1[,"COMCUMPRE2HARKM"],
TOTAL = FE.df1[,"TOTAL"],
LiveWtKg = FE.df1[,"LiveWtKg"],
Lat = FE.df1[,"Lat"],
NN = FE.df1[,"NN"],
Area = FE.df1[,"CmpxArea"],
HINDFOOT = FE.df1[,"HINDFOOT"]/100,
COMPREDOEHARKM = FE.df1[,"COMPREDOEHARKM"],
COMCUMPRE2DOEHARKM = FE.df1[,"COMCUMPRE2DOEHARKM"],
COMPRECNTDEN = FE.df1[,"COMPRECNTDEN"],
GrowSeas = round(FE.df1[,"GrowSeas"],3),
GrowSeas.sc = FE.df1[,"GrowSeas.sc"],
PrecEx10mm = FE.df1[,"PrecEx10mm"],
PrecEx20mm = FE.df1[,"PrecEx20mm"],
DSubZero = FE.df1[,"DSubZero"],
Prec95pct = FE.df1[,"Prec95pct"],
MonMxDT = FE.df1[,"MonMxDT"],
MMx1DPrec = FE.df1[,"MMx1DPrec"],
AMinT = FE.df1[,"AMinT"],
AMxT = FE.df1[,"AMxT"],
SumDur = round(FE.df1[,"SumDur"],3),
SumDur.sc = FE.df1[,"SumDur.sc"],
FrstDays = FE.df1[,"FrstDays"],
Int.DATE_YR = FE.df1[,"Int.DATE_YR"],
RepGroups = FE.df1[,"RepGroups"],
SexAgeGroup.f = FE.df1[,"SexAgeGroup.f"],
AGEGROUP = FE.df1[,"AGEGROUP"],
PREG = as.factor(FE.df1[,"PREG"]),
f_COMPLEX = FE.df1[,"f_COMPLEX"],
myID = FE.df1[,"myID"],
Month = FE.df1[,"DATE_MO"]))
#Data Stack
WholeWeight.stk = inla.stack(data = list(Weight = FE.df1$LiveWtKg), #Weight is dependent
A = list(1,1),
effects = FE.lst1,
tag = "base.0")Model for Whole Weight (Reconstructed Population COMRECPOP)
nPrior = list(theta=list(prior = "normal", param=c(0, 10)))
LiveWeight.frm = Weight ~ -1 + intercept1 +
f(COMRECPOP, #Recon Pop (total)
model="rw1",
constr = TRUE,
scale.model = TRUE,
replicate = RepGroups, #Independent repeat for sex-age groups
hyper = nPrior) +
f(myID,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(SexAgeGroup.f,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(PREG,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(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 ) +
DSubZero + TOTAL
#theta.m1 = LiveWeight.mod$internal.summary.hyperpar$mean
theta.m1 = c(-3.8573024980, -0.6627697675, -0.3888027095, -2.9775732650, -0.0308187818, 0.0267963821,
-1.1014780788, -0.4258434063, 0.0302421314)
LiveWeight.mod = inla(LiveWeight.frm,
data = inla.stack.data(WholeWeight.stk),
family = "gaussian",
verbose = TRUE,
control.fixed = list(prec = 0.001,
prec.intercept = 0.0001),
control.predictor = list(
A = inla.stack.A(WholeWeight.stk),
compute = TRUE,
link = 1),
control.mode = list(restart = TRUE, theta = theta.m1),
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(WholeWeight.stk), ", " verbose = TRUE, control.compute = list(dic = TRUE, cpo = TRUE, ", " waic = TRUE), control.predictor = list(A = inla.stack.A(WholeWeight.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 = theta.m1))")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 2.3925 19.7116 0.6437 22.7479
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept1 54.2813 0.2965 53.6992 54.2813 54.8630 54.2813 0
## DSubZero -0.8112 0.3203 -1.4401 -0.8112 -0.1828 -0.8112 0
## TOTAL 0.0351 0.0042 0.0268 0.0351 0.0434 0.0351 0
##
## Random effects:
## Name Model
## COMRECPOP RW1 model
## myID IID model
## SexAgeGroup.f IID model
## PREG IID model
## Month IID model
## Int.DATE_YR RW1 model
## GrowSeas RW1 model
## SumDur RW1 model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant
## Precision for the Gaussian observations 0.0209 0.0004 0.0201 0.0210
## Precision for COMRECPOP 0.5416 0.1390 0.3237 0.5227
## Precision for myID 1.4429 0.9797 0.5608 1.1474
## Precision for SexAgeGroup.f 0.0523 0.0080 0.0388 0.0515
## Precision for PREG 1.0085 0.3048 0.5376 0.9662
## Precision for Month 0.3427 0.0745 0.2253 0.3324
## Precision for Int.DATE_YR 0.6710 0.1987 0.3623 0.6442
## Precision for GrowSeas 1.1421 0.3489 0.5838 1.1019
## Precision for SumDur 1.1099 0.3485 0.5925 1.0536
## 0.975quant mode
## Precision for the Gaussian observations 0.0217 0.0211
## Precision for COMRECPOP 0.8666 0.4869
## Precision for myID 4.0456 0.7811
## Precision for SexAgeGroup.f 0.0699 0.0499
## Precision for PREG 1.7263 0.8868
## Precision for Month 0.5168 0.3117
## Precision for Int.DATE_YR 1.1354 0.5938
## Precision for GrowSeas 1.9434 1.0231
## Precision for SumDur 1.9468 0.9508
##
## Expected number of effective parameters(std dev): 138.59(0.00)
## Number of equivalent replicates : 27.40
##
## Deviance Information Criterion (DIC) ...............: 25769.61
## Deviance Information Criterion (DIC, saturated) ....: 4139.32
## Effective number of parameters .....................: 138.59
##
## Watanabe-Akaike information criterion (WAIC) ...: 25782.36
## Effective number of parameters .................: 146.06
##
## Marginal log-Likelihood: -13167.52
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computedModel for Whole Weight (Reconstructed Doe Population COMRECDOEPOP)
nPrior = list(theta=list(prior = "normal", param=c(0, 10)))
LiveWeight.doe = Weight ~ -1 + intercept1 +
f(COMRECDOEPOP, #Recon Doe Pop
model="rw1",
constr = TRUE,
scale.model = TRUE,
replicate = RepGroups,
hyper = nPrior) +
f(myID,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(SexAgeGroup.f, #Does only = age groups
model="iid",
constr = TRUE,
hyper = nPrior) +
f(PREG,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(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 ) +
DSubZero + TOTAL
#theta.m2 = LiveWeight.mod.doe$internal.summary.hyperpar$mean
theta.m2 = c(-3.943684548, -0.657787303, 0.044425219, -2.970071171, -0.006181066, 0.010046463, -1.060661541,
-0.327049237, 0.044311737)
LiveWeight.mod.doe = inla(LiveWeight.doe,
data = inla.stack.data(WholeWeight.stk),
family = "gaussian",
verbose = TRUE,
control.fixed = list(prec = 0.001,
prec.intercept = 0.0001),
control.predictor = list(
A = inla.stack.A(WholeWeight.stk),
compute = TRUE,
link = 1),
control.mode = list(restart = TRUE, theta = theta.m2),
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.doe)##
## Call:
## c("inla(formula = LiveWeight.doe, family = \"gaussian\", data = inla.stack.data(WholeWeight.stk), ", " verbose = TRUE, control.compute = list(dic = TRUE, cpo = TRUE, ", " waic = TRUE), control.predictor = list(A = inla.stack.A(WholeWeight.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 = theta.m2))")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 2.2271 20.5259 0.6542 23.4072
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept1 54.2053 0.3059 53.6047 54.2053 54.8054 54.2053 0
## DSubZero -0.7988 0.3254 -1.4376 -0.7988 -0.1605 -0.7988 0
## TOTAL 0.0346 0.0044 0.0260 0.0346 0.0432 0.0346 0
##
## Random effects:
## Name Model
## COMRECDOEPOP RW1 model
## myID IID model
## SexAgeGroup.f IID model
## PREG IID model
## Month IID model
## Int.DATE_YR RW1 model
## GrowSeas RW1 model
## SumDur RW1 model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant
## Precision for the Gaussian observations 0.0205 0.0012 0.0189 0.0203
## Precision for COMRECDOEPOP 0.5458 0.1568 0.3043 0.5233
## Precision for myID 1.0658 0.3254 0.5547 1.0242
## Precision for SexAgeGroup.f 0.0523 0.0083 0.0381 0.0516
## Precision for PREG 1.0438 0.3227 0.5524 0.9965
## Precision for Month 0.3548 0.0805 0.2094 0.3519
## Precision for Int.DATE_YR 0.7838 0.2435 0.3974 0.7546
## Precision for GrowSeas 1.1152 0.3518 0.5945 1.0577
## Precision for SumDur 1.1164 0.3459 0.5982 1.0624
## 0.975quant mode
## Precision for the Gaussian observations 0.0233 0.0194
## Precision for COMRECDOEPOP 0.9160 0.4813
## Precision for myID 1.8232 0.9450
## Precision for SexAgeGroup.f 0.0706 0.0501
## Precision for PREG 1.8097 0.9088
## Precision for Month 0.5223 0.3479
## Precision for Int.DATE_YR 1.3449 0.6978
## Precision for GrowSeas 1.9626 0.9531
## Precision for SumDur 1.9450 0.9633
##
## Expected number of effective parameters(std dev): 112.23(0.00)
## Number of equivalent replicates : 33.84
##
## Deviance Information Criterion (DIC) ...............: 25786.41
## Deviance Information Criterion (DIC, saturated) ....: 3828.04
## Effective number of parameters .....................: 112.23
##
## Watanabe-Akaike information criterion (WAIC) ...: 25787.84
## Effective number of parameters .................: 110.39
##
## Marginal log-Likelihood: -13172.04
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computedModel for Whole Weight (Aerial Counts COMPRECNTDEN)
nPrior = list(theta=list(prior = "normal", param=c(0, 10)))
LiveWeight.aerial = Weight ~ -1 + intercept1 +
f(COMPRECNTDEN,
model="rw1",
constr = TRUE,
scale.model = TRUE,
replicate = RepGroups,
hyper = nPrior) +
f(myID,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(SexAgeGroup.f, #Does only = age groups
model="iid",
constr = TRUE,
hyper = nPrior) +
f(PREG,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(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 ) +
DSubZero + TOTAL
#theta.m3 = LiveWeight.mod.aerial$internal.summary.hyperpar$mean
theta.m3 = c(-3.949758515, -0.668447978, 0.023023690, -2.971971613, -0.026000289, 0.008109743,
-1.073993338, -0.302061119, 0.039326732)
LiveWeight.mod.aerial = inla(LiveWeight.aerial,
data = inla.stack.data(WholeWeight.stk),
family = "gaussian",
verbose = TRUE,
control.fixed = list(prec = 0.001,
prec.intercept = 0.0001),
control.predictor = list(
A = inla.stack.A(WholeWeight.stk),
compute = TRUE,
link = 1),
control.mode = list(restart = TRUE, theta = theta.m3),
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.aerial)##
## Call:
## c("inla(formula = LiveWeight.aerial, family = \"gaussian\", data = inla.stack.data(WholeWeight.stk), ", " verbose = TRUE, control.compute = list(dic = TRUE, cpo = TRUE, ", " waic = TRUE), control.predictor = list(A = inla.stack.A(WholeWeight.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 = theta.m3))")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 2.7528 20.5339 0.7900 24.0767
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept1 54.2246 0.3091 53.6177 54.2246 54.8310 54.2246 0
## DSubZero -0.8990 0.3383 -1.5631 -0.8990 -0.2354 -0.8990 0
## TOTAL 0.0343 0.0044 0.0256 0.0343 0.0430 0.0343 0
##
## Random effects:
## Name Model
## COMPRECNTDEN RW1 model
## myID IID model
## SexAgeGroup.f IID model
## PREG IID model
## Month IID model
## Int.DATE_YR RW1 model
## GrowSeas RW1 model
## SumDur RW1 model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant
## Precision for the Gaussian observations 0.0197 0.0007 0.0186 0.0196
## Precision for COMPRECNTDEN 0.7481 0.2383 0.3773 0.7169
## Precision for myID 1.0409 0.3190 0.5445 0.9985
## Precision for SexAgeGroup.f 0.0524 0.0083 0.0383 0.0517
## Precision for PREG 1.0266 0.3304 0.5282 0.9766
## Precision for Month 0.3461 0.0808 0.2220 0.3339
## Precision for Int.DATE_YR 0.4780 0.1248 0.2737 0.4647
## Precision for GrowSeas 1.1037 0.3542 0.5807 1.0459
## Precision for SumDur 1.1082 0.3499 0.5843 1.0536
## 0.975quant mode
## Precision for the Gaussian observations 0.0212 0.0193
## Precision for COMPRECNTDEN 1.3035 0.6573
## Precision for myID 1.7877 0.9183
## Precision for SexAgeGroup.f 0.0708 0.0501
## Precision for PREG 1.8126 0.8845
## Precision for Month 0.5374 0.3096
## Precision for Int.DATE_YR 0.7621 0.4394
## Precision for GrowSeas 1.9553 0.9409
## Precision for SumDur 1.9470 0.9536
##
## Expected number of effective parameters(std dev): 112.39(0.00)
## Number of equivalent replicates : 33.79
##
## Deviance Information Criterion (DIC) ...............: 25823.70
## Deviance Information Criterion (DIC, saturated) ....: 3842.26
## Effective number of parameters .....................: 112.39
##
## Watanabe-Akaike information criterion (WAIC) ...: 25825.45
## Effective number of parameters .................: 110.84
##
## Marginal log-Likelihood: -13225.09
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computedCompare DIC and WAIC
#DIC
LiveWeight.mod$dic$dic## [1] 25769.61LiveWeight.mod.doe$dic$dic## [1] 25786.41LiveWeight.mod.aerial$dic$dic## [1] 25823.7#WAIC
LiveWeight.mod$waic$waic## [1] 25782.36LiveWeight.mod.doe$waic$waic## [1] 25787.84LiveWeight.mod.aerial$waic$waic## [1] 25825.45Whole Weight Density Dependence (Reconstructed Population)
Horizontal line is global population mean whole weight.
mic.d = as.data.frame(LiveWeight.mod$summary.random$COMRECPOP[,c(1:6)])
names(mic.d) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975")
length(unique(mic.d$ID))## [1] 13length(mic.d$ID)## [1] 78mic.d$Sex = rep(1:2, each = 39)
mic.d$ID2 = as.factor(
rep(c("Fawn (m)", "Yearling (m)", "Adult (m)",
"Fawn (f)", "Yearling (f)", "Adult (f)"),
each = length(unique(mic.d$ID))))
mic.d$Meane = mic.d$Mean + LiveWeight.mod$summary.fixed[1,1]
mic.d$Q025e = mic.d$Q025 + LiveWeight.mod$summary.fixed[1,1]
mic.d$Q975e = mic.d$Q975 + LiveWeight.mod$summary.fixed[1,1]
mic.d$ID = mic.d$ID/30.63 #Area of Complex 1
range(mic.d$ID)## [1] 8.651649 38.132550male.set = mic.d %>% filter(Sex == 1)
range(male.set$Meane)## [1] 49.11197 57.14511female.set = mic.d %>% filter(Sex == 2)
range(female.set$Meane)## [1] 52.06804 56.20191ggplot(male.set, aes(ID, Meane, group = ID2)) +
geom_point(shape=1) +
geom_smooth(data= male.set,
aes(ID, Meane,
group = male.set$ID2,
col = factor(male.set$ID2),
linetype= factor(male.set$ID2)),
se=FALSE,
span = 0.9,
method="loess",
lwd = 1) +
geom_hline(yintercept = LiveWeight.mod$summary.fixed[1,1],
linetype = "solid",
col = "darkgray",
size = 0.5) +
scale_linetype_manual(name = "Age Group",
values=c("solid", "longdash", "dotted")) +
scale_colour_manual(name = "Age Group",
values=c("red", "brown", "blue")) +
ylab("Liveweight (kg)") +
xlab(expression(Deer~km^2)) +
ggtitle("Male (COMRECPOP)") +
theme_classic() +
scale_y_continuous(breaks = seq(45, 60, 5), limits = c(45, 60)) +
scale_x_continuous(breaks = seq(8, 40, 5), limits =c(8,40)) +
theme(axis.text=element_text(size=16),
legend.title = element_text(size=16, face="bold"),
legend.position = c(0.3,0.2),
legend.key = element_rect(fill = "white", linetype=0),
legend.background = element_rect(fill = "white", linetype=0),
plot.title = element_text(face="bold", size = 20, hjust=0.5),
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))
ggplot(female.set, aes(ID, Meane, group = ID2)) +
geom_point(shape=1) +
geom_smooth(data= female.set,
aes(ID, Meane,
group = female.set$ID2,
col = factor(female.set$ID2),
linetype= factor(female.set$ID2)),
se=FALSE,
span = 0.9,
method="loess",
lwd = 1) +
geom_hline(yintercept = LiveWeight.mod$summary.fixed[1,1],
linetype = "solid",
col = "darkgray",
size = 0.5) +
scale_linetype_manual(name = "Age Group",
values=c("solid", "longdash", "dotted")) +
scale_colour_manual(name = "Age Group",
values=c("red", "brown", "blue")) +
ylab("Liveweight (kg)") +
xlab(expression(Deer~km^2)) +
ggtitle("Female (COMRECPOP)") +
theme_classic() +
scale_y_continuous(breaks = seq(50, 60, 5), limits = c(50, 60)) +
scale_x_continuous(breaks = seq(8, 40, 5), limits =c(8,40)) +
theme(axis.text=element_text(size=16),
legend.title = element_text(size=16, face="bold"),
legend.position = c(0.3,0.2),
legend.key = element_rect(fill = "white", linetype=0),
legend.background = element_rect(fill = "white", linetype=0),
plot.title = element_text(face="bold", size = 20, hjust=0.5),
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))
Whole Weight Density Dependence (Aerial Counts, COMPRECNTDEN)
Horizontal line is global population mean whole weight.
mic.d = as.data.frame(LiveWeight.mod.aerial$summary.random$COMPRECNTDEN[,c(1:6)])
names(mic.d) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975")
length(unique(mic.d$ID))## [1] 17length(mic.d$ID)## [1] 102mic.d$Sex = rep(1:2, each = 51)
mic.d$ID2 = as.factor(
rep(c("Fawn (m)", "Yearling (m)", "Adult (m)",
"Fawn (f)", "Yearling (f)", "Adult (f)"),
each = length(unique(mic.d$ID))))
mic.d$Meane = mic.d$Mean + LiveWeight.mod$summary.fixed[1,1]
mic.d$Q025e = mic.d$Q025 + LiveWeight.mod$summary.fixed[1,1]
mic.d$Q975e = mic.d$Q975 + LiveWeight.mod$summary.fixed[1,1]
male.set = mic.d %>% filter(Sex == 1)
range(male.set$Meane)## [1] 50.52641 55.48707female.set = mic.d %>% filter(Sex == 2)
range(female.set$Meane)## [1] 52.45358 55.30290ggplot(male.set, aes(ID, Meane, group = ID2)) +
geom_point(shape=1) +
geom_smooth(data= male.set,
aes(ID, Meane,
group = male.set$ID2,
col = factor(male.set$ID2),
linetype= factor(male.set$ID2)),
se=FALSE,
span = 0.9,
method="loess",
lwd = 1) +
geom_hline(yintercept = LiveWeight.mod$summary.fixed[1,1],
linetype = "solid",
col = "darkgray",
size = 0.5) +
scale_linetype_manual(name = "Age Group",
values=c("solid", "longdash", "dotted")) +
scale_colour_manual(name = "Age Group",
values=c("red", "brown", "blue")) +
ylab("Liveweight (kg)") +
xlab(expression(Deer~km^2)) +
ggtitle("Male (COMPRECNTDEN)") +
theme_classic() +
scale_y_continuous(breaks = seq(45, 60, 5), limits = c(45, 60)) +
scale_x_continuous(breaks = seq(4, 40, 5), limits =c(4,40)) +
theme(axis.text=element_text(size=16),
legend.title = element_text(size=16, face="bold"),
legend.position = c(0.3,0.2),
legend.key = element_rect(fill = "white", linetype=0),
legend.background = element_rect(fill = "white", linetype=0),
plot.title = element_text(face="bold", size = 20, hjust=0.5),
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))
ggplot(female.set, aes(ID, Meane, group = ID2)) +
geom_point(shape=1) +
geom_smooth(data= female.set,
aes(ID, Meane,
group = female.set$ID2,
col = factor(female.set$ID2),
linetype= factor(female.set$ID2)),
se=FALSE,
span = 0.9,
method="loess",
lwd = 1) +
geom_hline(yintercept = LiveWeight.mod$summary.fixed[1,1],
linetype = "solid",
col = "darkgray",
size = 0.5) +
scale_linetype_manual(name = "Age Group",
values=c("solid", "longdash", "dotted")) +
scale_colour_manual(name = "Age Group",
values=c("red", "brown", "blue")) +
ylab("Liveweight (kg)") +
xlab(expression(Deer~km^2)) +
ggtitle("Female (COMPRECNTDEN)") +
theme_classic() +
scale_y_continuous(breaks = seq(50, 60, 5), limits = c(50, 60)) +
scale_x_continuous(breaks = seq(4, 40, 5), limits =c(4,40)) +
theme(axis.text=element_text(size=16),
legend.title = element_text(size=16, face="bold"),
legend.position = c(0.3,0.2),
legend.key = element_rect(fill = "white", linetype=0),
legend.background = element_rect(fill = "white", linetype=0),
plot.title = element_text(face="bold", size = 20, hjust=0.5),
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))
Model Fitted Estimates.
Temp.data = FE.df1
idat = inla.stack.index(WholeWeight.stk, "base.0")$data #Indexregion data from model
Temp.data$Pred.Weight = LiveWeight.mod$summary.fitted.values$mean[idat] #Fitted values for region
Temp.data$Low = LiveWeight.mod$summary.fitted.values$`0.025quant`[idat]
Temp.data$High = LiveWeight.mod$summary.fitted.values$`0.975quant`[idat]
LW.fit = as.data.frame(
Temp.data %>%
group_by(SexAgeGroup.f) %>%
summarise(Mean = mean(Pred.Weight, na.rm=T),
Low = mean(Low, na.rm=T),
High = mean(High, na.rm=T)))
LW.fit$Age = c("Fawn (m)", "Yearling (m)", "Adult (m)",
"Fawn (f)", "Yearling (f)", "Adult (f)")
LW.fit$Age = ordered(factor(LW.fit$Age),
levels = c("Fawn (m)", "Yearling (m)", "Adult (m)",
"Fawn (f)", "Yearling (f)", "Adult (f)"))
ggplot(LW.fit, aes(factor(Age), Mean)) +
geom_bar(stat = "identity", fill = "lightgray") +
geom_point(size=2, pch=19, col = "red") +
geom_linerange(aes(ymin=Low, ymax=High), colour="black") +
theme_classic() +
ylab("LiveWeight (kg)") +
xlab( "") +
ggtitle("Estimated LiveWeight \n(Complex 1, COMRECPOP)") +
theme(axis.text=element_text(size=16),
legend.title = element_text(size=16, face="bold"),
legend.position = c(0.7,0.7),
legend.key = element_rect(fill = "white", linetype=0),
legend.background = element_rect(fill = "white", linetype=0),
plot.title = element_text(face="bold", size = 20, hjust=0.5),
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))
Temp.data$Class = ifelse(Temp.data$SexAgeGeoup == "21", "Fawn (f)",
ifelse(Temp.data$SexAgeGeoup == "22", "Yearling (f)",
ifelse(Temp.data$SexAgeGeoup == "23", "Adult (f)",
ifelse(Temp.data$SexAgeGeoup == "11", "Fawn (m)",
ifelse(Temp.data$SexAgeGeoup == "12", "Yearling (m)",
ifelse(Temp.data$SexAgeGeoup == "13", "Adult (m)", NA))))))
Temp.data$Class = ordered(factor(Temp.data$Class,
levels = c("Fawn (m)", "Yearling (m)", "Adult (m)",
"Fawn (f)", "Yearling (f)", "Adult (f)")))
Temp.data$Density = Temp.data$COMRECPOP/30.63
COMRECPOP.plt = Temp.data %>% select(Class, Pred.Weight)
hline_dat1 = as.data.frame(
COMRECPOP.plt %>%
group_by(Class) %>%
summarise(Med = mean(Pred.Weight)))
ggplot(Temp.data, aes(Density, Pred.Weight)) +
geom_point(shape=1) +
geom_smooth(method="loess",
se=FALSE,
col="black",
lwd = 1) +
facet_wrap(~Class) +
geom_hline(data=hline_dat1,
aes(yintercept=Med),
col = "darkgray",
size = 0.5) +
ylab("Live Weight (kg)") +
xlab(expression(Deer~km^2)) +
ggtitle("Estimated Live Weight \n(Complex 1, COMRECPOP)") +
theme_classic() +
theme(axis.text=element_text(size=16),
plot.title = element_text(face="bold", size = 20, hjust=0.5),
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))
Sex-Age Class Effect
Horizontal line is population average.
mic.d = as.data.frame(LiveWeight.mod$summary.random$SexAgeGroup.f)
names(mic.d) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975", "mode", "kld")
mic.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]
range(mic.d$Mean)## [1] 35.88074 78.75895ggplot(mic.d, aes(factor(ID), y=Mean)) +
geom_point(size=2, pch=19, col = "red") +
geom_linerange(aes(ymin=Q025, ymax=Q975), colour="black") +
geom_hline(yintercept = LiveWeight.mod$summary.fixed[1,1],
linetype = "solid",
col = "darkgray",
size = 0.5) +
theme_classic() +
xlab(" ") +
ylab("LiveWeight (kg)") +
ggtitle("Sex-Age Group Effect") +
scale_x_discrete(labels = c("Fawn (m)", "Yearling (m)", "Adult (m)",
"Fawn (f)", "Yearling (f)", "Adult (f)")) +
scale_y_continuous(breaks = seq(35, 80, 5), limits = c(35, 80)) +
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. Horizontal line is population mean.
mic.d = as.data.frame(LiveWeight.mod$summary.random$Month)
names(mic.d) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975", "mode", "kld")
mic.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]
range(mic.d$Mean)## [1] 48.22542 58.79737ggplot(mic.d, aes(factor(ID), y=Mean)) +
geom_point(size=2, pch=19, col = "red") +
geom_linerange(aes(ymin=Q025, ymax=Q975), colour="black") +
geom_hline(yintercept = LiveWeight.mod$summary.fixed[1,1],
linetype = "solid",
col = "darkgray",
size = 0.5) +
theme_classic() +
xlab(" ") +
ylab("LiveWeight") +
ggtitle("Month Effect") +
scale_x_discrete(labels = c("Jan", "Feb", "March",
"Oct", "Nov", "Dec")) +
scale_y_continuous(breaks = seq(45, 60, 5), limits = c(45, 60)) +
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$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]
range(mic.d$Mean)## [1] 51.25916 55.74334mic.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 = LiveWeight.mod$summary.fixed[1,1],
linetype = "solid",
col = "darkgray",
size = 0.5) +
scale_y_continuous(breaks = seq(45, 60, 5), limits = c(45, 60)) +
ylab("LiveWeight (kg)") +
xlab(" ") +
scale_x_continuous(breaks = 4:20,labels=paste0(mic.d$YLab)) +
ggtitle("Year 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 Horizontal line is population mean.
mic.d = as.data.frame(LiveWeight.mod$summary.random$PREG)
names(mic.d) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975", "mode", "kld")
mic.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]
range(mic.d$Mean)## [1] 53.51592 55.04675ggplot(mic.d, aes(factor(ID), y=Mean)) +
geom_point(size=2, pch=19, col = "red") +
geom_linerange(aes(ymin=Q025, ymax=Q975), colour="black") +
geom_hline(yintercept = LiveWeight.mod$summary.fixed[1,1],
linetype = "solid",
col = "darkgray",
size = 0.5) +
scale_y_continuous(breaks = seq(50, 60, 5), limits = c(50, 60)) +
theme_classic() +
xlab(" ") +
ylab("LiveWeight (kg)") +
ggtitle("Doe Pregnancy Effect") +
scale_x_discrete(labels = c("Yes", "No")) +
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$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]
range(mic.d$Mean)## [1] 53.24481 55.03222Climate.df$rLs_GrowSeas = round(Climate.df[,"Ls_GrowSeas"],3)
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 = LiveWeight.mod$summary.fixed[1,1],
linetype = "solid",
col = "darkgray",
size = 0.5) +
scale_y_continuous(breaks = seq(50, 60, 5), limits = c(50, 60)) +
ylab("LiveWeight (kg)") +
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))
2.2 TOTAL Kistner
Organize Data as List Object
#Data Stack
Kistner.stk = inla.stack(data = list(Kistner = FE.df1$TOTAL), #Weight is dependent
A = list(1,1),
effects = FE.lst1,
tag = "Kist.0")Model for TOTAL KISTNER (Reconstructed Population COMRECPOP)
Kistner.frm = Kistner ~ -1 + intercept1 +
f(COMRECPOP, #Recon Pop (total)
model="rw1",
constr = TRUE,
scale.model = TRUE,
replicate = RepGroups,
hyper = nPrior) +
f(myID,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(SexAgeGroup.f,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(PREG,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(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 ) +
DSubZero + LiveWtKg
#theta.m1 = Kistner.mod$internal.summary.hyperpar$mean
theta.m1 = c(-5.31381402, -0.92406946, -0.02804894, -2.36447538, -0.04321129, -0.27580659, -2.51593932,
-1.77086861, -0.39001213)
Kistner.mod = inla(Kistner.frm,
data = inla.stack.data(Kistner.stk),
family = "gaussian",
verbose = TRUE,
control.fixed = list(prec = 0.001,
prec.intercept = 0.0001),
control.predictor = list(
A = inla.stack.A(Kistner.stk),
compute = TRUE,
link = 1),
control.mode = list(restart = TRUE, theta = theta.m1),
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(Kistner.mod)##
## Call:
## c("inla(formula = Kistner.frm, family = \"gaussian\", data = inla.stack.data(Kistner.stk), ", " verbose = TRUE, control.compute = list(dic = TRUE, cpo = TRUE, ", " waic = TRUE), control.predictor = list(A = inla.stack.A(Kistner.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 = theta.m1))")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 2.4406 23.2579 0.9655 26.6640
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept1 54.5162 1.4034 51.7610 54.5162 57.2692 54.5162 0
## DSubZero -1.6098 0.5116 -2.6142 -1.6098 -0.6061 -1.6098 0
## LiveWtKg 0.2258 0.0227 0.1811 0.2258 0.2704 0.2258 0
##
## Random effects:
## Name Model
## COMRECPOP RW1 model
## myID IID model
## SexAgeGroup.f IID model
## PREG IID model
## Month IID model
## Int.DATE_YR RW1 model
## GrowSeas RW1 model
## SumDur RW1 model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant
## Precision for the Gaussian observations 0.0048 0.0003 0.0042 0.0048
## Precision for COMRECPOP 0.4182 0.1119 0.2382 0.4050
## Precision for myID 1.0285 0.3279 0.5455 0.9743
## Precision for SexAgeGroup.f 0.0941 0.0177 0.0652 0.0921
## Precision for PREG 0.7792 0.2457 0.4201 0.7378
## Precision for Month 0.0820 0.0144 0.0558 0.0815
## Precision for Int.DATE_YR 0.1731 0.0553 0.0894 0.1648
## Precision for GrowSeas 0.7331 0.2638 0.3501 0.6890
## Precision for SumDur 0.6710 0.2563 0.2806 0.6360
## 0.975quant mode
## Precision for the Gaussian observations 0.0052 0.0050
## Precision for COMRECPOP 0.6755 0.3801
## Precision for myID 1.8194 0.8755
## Precision for SexAgeGroup.f 0.1343 0.0878
## Precision for PREG 1.3747 0.6624
## Precision for Month 0.1120 0.0809
## Precision for Int.DATE_YR 0.3046 0.1495
## Precision for GrowSeas 1.3718 0.6094
## Precision for SumDur 1.2729 0.5630
##
## Expected number of effective parameters(std dev): 48.75(0.00)
## Number of equivalent replicates : 60.86
##
## Deviance Information Criterion (DIC) ...............: 24234.02
## Deviance Information Criterion (DIC, saturated) ....: 3082.93
## Effective number of parameters .....................: 48.75
##
## Watanabe-Akaike information criterion (WAIC) ...: 24236.36
## Effective number of parameters .................: 50.24
##
## Marginal log-Likelihood: -12445.70
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computedModel for TOTAL KISTNER (Reconstructed Doe Population COMRECDOEPOP)
Kistner.doe = Kistner ~ -1 + intercept1 +
f(COMRECDOEPOP, #Recon Doe Pop
model="rw1",
constr = TRUE,
scale.model = TRUE,
replicate = RepGroups,
hyper = nPrior) +
f(myID,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(SexAgeGroup.f, #Does only = age groups
model="iid",
constr = TRUE,
hyper = nPrior) +
f(PREG,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(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 ) +
DSubZero + LiveWtKg
#theta.m2 = Kistner.mod.doe$internal.summary.hyperpar$mean
theta.m2 = c(-5.31381402, -0.92406946, -0.02804894, -2.36447538, -0.04321129, -0.27580659, -2.51593932,
-1.77086861, -0.39001213)
Kistner.mod.doe = inla(Kistner.doe,
data = inla.stack.data(Kistner.stk),
family = "gaussian",
verbose = TRUE,
control.fixed = list(prec = 0.001,
prec.intercept = 0.0001),
control.predictor = list(
A = inla.stack.A(Kistner.stk),
compute = TRUE,
link = 1),
control.mode = list(restart = TRUE, theta = theta.m2),
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(Kistner.mod.doe)##
## Call:
## c("inla(formula = Kistner.doe, family = \"gaussian\", data = inla.stack.data(Kistner.stk), ", " verbose = TRUE, control.compute = list(dic = TRUE, cpo = TRUE, ", " waic = TRUE), control.predictor = list(A = inla.stack.A(Kistner.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 = theta.m2))")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 2.2879 28.7818 0.5682 31.6379
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept1 54.5883 1.4148 51.8105 54.5882 57.3638 54.5883 0
## DSubZero -1.6008 0.5213 -2.6243 -1.6008 -0.5782 -1.6008 0
## LiveWtKg 0.2244 0.0230 0.1792 0.2244 0.2695 0.2244 0
##
## Random effects:
## Name Model
## COMRECDOEPOP RW1 model
## myID IID model
## SexAgeGroup.f IID model
## PREG IID model
## Month IID model
## Int.DATE_YR RW1 model
## GrowSeas RW1 model
## SumDur RW1 model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant
## Precision for the Gaussian observations 0.0050 0.0001 0.0047 0.0050
## Precision for COMRECDOEPOP 0.3991 0.1054 0.2327 0.3853
## Precision for myID 1.0214 0.3372 0.5224 0.9662
## Precision for SexAgeGroup.f 0.0935 0.0172 0.0628 0.0926
## Precision for PREG 0.7783 0.2339 0.4093 0.7489
## Precision for Month 0.0821 0.0146 0.0554 0.0816
## Precision for Int.DATE_YR 0.1638 0.0531 0.0868 0.1546
## Precision for GrowSeas 0.6878 0.2330 0.3404 0.6513
## Precision for SumDur 0.6607 0.2404 0.2997 0.6250
## 0.975quant mode
## Precision for the Gaussian observations 0.0052 0.0050
## Precision for COMRECDOEPOP 0.6449 0.3591
## Precision for myID 1.8433 0.8664
## Precision for SexAgeGroup.f 0.1300 0.0913
## Precision for PREG 1.3213 0.6928
## Precision for Month 0.1122 0.0812
## Precision for Int.DATE_YR 0.2928 0.1380
## Precision for GrowSeas 1.2434 0.5843
## Precision for SumDur 1.2341 0.5565
##
## Expected number of effective parameters(std dev): 49.05(0.00)
## Number of equivalent replicates : 60.50
##
## Deviance Information Criterion (DIC) ...............: 24221.37
## Deviance Information Criterion (DIC, saturated) ....: 3030.98
## Effective number of parameters .....................: 49.05
##
## Watanabe-Akaike information criterion (WAIC) ...: 24223.00
## Effective number of parameters .................: 49.82
##
## Marginal log-Likelihood: -12446.02
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computedModel for TOTAL KISTNER (Aerial Counts COMPRECNTDEN)
Kistner.aerial = Kistner ~ -1 + intercept1 +
f(COMPRECNTDEN, #Recon Doe Pop
model="rw1",
constr = TRUE,
scale.model = TRUE,
replicate = RepGroups,
hyper = nPrior) +
f(myID,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(SexAgeGroup.f, #Does only = age groups
model="iid",
constr = TRUE,
hyper = nPrior) +
f(PREG,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(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 ) +
DSubZero + LiveWtKg
#theta.m3 = Kistner.mod.aerial$internal.summary.hyperpar$mean
theta.m3 = c(-3.949758515, -0.668447978, 0.023023690, -2.971971613, -0.026000289, 0.008109743,
-1.073993338, -0.302061119, 0.039326732)
Kistner.mod.aerial = inla(Kistner.aerial,
data = inla.stack.data(Kistner.stk),
family = "gaussian",
verbose = TRUE,
control.fixed = list(prec = 0.001,
prec.intercept = 0.0001),
control.predictor = list(
A = inla.stack.A(Kistner.stk),
compute = TRUE,
link = 1),
control.mode = list(restart = TRUE, theta = theta.m3),
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(Kistner.mod.aerial)##
## Call:
## c("inla(formula = Kistner.aerial, family = \"gaussian\", data = inla.stack.data(Kistner.stk), ", " verbose = TRUE, control.compute = list(dic = TRUE, cpo = TRUE, ", " waic = TRUE), control.predictor = list(A = inla.stack.A(Kistner.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 = theta.m3))")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 2.1695 24.4960 0.8238 27.4894
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept1 55.5460 1.4472 52.7046 55.5459 58.3849 55.5460 0
## DSubZero -1.5927 0.5199 -2.6134 -1.5927 -0.5729 -1.5927 0
## LiveWtKg 0.2096 0.0233 0.1638 0.2096 0.2553 0.2096 0
##
## Random effects:
## Name Model
## COMPRECNTDEN RW1 model
## myID IID model
## SexAgeGroup.f IID model
## PREG IID model
## Month IID model
## Int.DATE_YR RW1 model
## GrowSeas RW1 model
## SumDur RW1 model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant
## Precision for the Gaussian observations 0.0048 0.0001 0.0046 0.0048
## Precision for COMPRECNTDEN 0.3268 0.0914 0.1864 0.3134
## Precision for myID 1.0630 0.3277 0.5559 1.0180
## Precision for SexAgeGroup.f 0.0987 0.0179 0.0667 0.0977
## Precision for PREG 0.6820 0.1921 0.3828 0.6558
## Precision for Month 0.0798 0.0149 0.0565 0.0777
## Precision for Int.DATE_YR 0.1637 0.0497 0.0846 0.1578
## Precision for GrowSeas 0.7552 0.2633 0.3743 0.7099
## Precision for SumDur 0.7180 0.2580 0.3305 0.6796
## 0.975quant mode
## Precision for the Gaussian observations 0.0051 0.0048
## Precision for COMPRECNTDEN 0.5422 0.2884
## Precision for myID 1.8331 0.9337
## Precision for SexAgeGroup.f 0.1367 0.0963
## Precision for PREG 1.1329 0.6066
## Precision for Month 0.1147 0.0730
## Precision for Int.DATE_YR 0.2783 0.1463
## Precision for GrowSeas 1.3963 0.6289
## Precision for SumDur 1.3336 0.6066
##
## Expected number of effective parameters(std dev): 47.40(0.00)
## Number of equivalent replicates : 62.60
##
## Deviance Information Criterion (DIC) ...............: 24228.84
## Deviance Information Criterion (DIC, saturated) ....: 2931.30
## Effective number of parameters .....................: 47.40
##
## Watanabe-Akaike information criterion (WAIC) ...: 24229.03
## Effective number of parameters .................: 46.79
##
## Marginal log-Likelihood: -12484.17
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computedCompare DIC and WAIC
#DIC
Kistner.mod$dic$dic## [1] 24234.02Kistner.mod.doe$dic$dic## [1] 24221.37Kistner.mod.aerial$dic$dic## [1] 24228.84#WAIC
Kistner.mod$waic$waic## [1] 24236.36Kistner.mod.doe$waic$waic## [1] 24223Kistner.mod.aerial$waic$waic## [1] 24229.03Total Kistner Density Dependence (Reconstructed Population)
Horizontal line is global population mean Total Kistner.
mic.d = as.data.frame(Kistner.mod$summary.random$COMRECPOP[,c(1:6)])
names(mic.d) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975")
length(unique(mic.d$ID))## [1] 13length(mic.d$ID)## [1] 78mic.d$Sex = rep(1:2, each = 39)
mic.d$ID2 = as.factor(
rep(c("Fawn (m)", "Yearling (m)", "Adult (m)",
"Fawn (f)", "Yearling (f)", "Adult (f)"),
each = length(unique(mic.d$ID))))
mic.d$Meane = mic.d$Mean + Kistner.mod.aerial$summary.fixed[1,1]
mic.d$Q025e = mic.d$Q025 + Kistner.mod.aerial$summary.fixed[1,1]
mic.d$Q975e = mic.d$Q975 + Kistner.mod.aerial$summary.fixed[1,1]
mic.d$ID = mic.d$ID/30.63 #Area of Complex 1
range(mic.d$ID)## [1] 8.651649 38.132550male.set = mic.d %>% filter(Sex == 1)
range(male.set$Meane)## [1] 50.43272 61.75611female.set = mic.d %>% filter(Sex == 2)
range(female.set$Meane)## [1] 50.58203 57.83789ggplot(male.set, aes(ID, Meane, group = ID2)) +
geom_point(shape=1) +
geom_smooth(data= male.set,
aes(ID, Meane,
group = male.set$ID2,
col = factor(male.set$ID2),
linetype= factor(male.set$ID2)),
se=FALSE,
span = 0.9,
method="loess",
lwd = 1) +
geom_hline(yintercept = Kistner.mod.aerial$summary.fixed[1,1],
linetype = "solid",
col = "darkgray",
size = 0.5) +
scale_linetype_manual(name = "Age Group",
values=c("solid", "longdash", "dotted")) +
scale_colour_manual(name = "Age Group",
values=c("red", "brown", "blue")) +
ylab("Kistner (kg)") +
xlab(expression(Deer~km^2)) +
ggtitle("Male (COMRECPOP)") +
theme_classic() +
scale_y_continuous(breaks = seq(50, 65, 5), limits = c(50, 65)) +
scale_x_continuous(breaks = seq(8, 40, 5), limits =c(8,40)) +
theme(axis.text=element_text(size=16),
legend.title = element_text(size=16, face="bold"),
legend.position = c(0.2,0.8),
legend.key = element_rect(fill = "white", linetype=0),
legend.background = element_rect(fill = "white", linetype=0),
plot.title = element_text(face="bold", size = 20, hjust=0.5),
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))
ggplot(female.set, aes(ID, Meane, group = ID2)) +
geom_point(shape=1) +
geom_smooth(data= female.set,
aes(ID, Meane,
group = female.set$ID2,
col = factor(female.set$ID2),
linetype= factor(female.set$ID2)),
se=FALSE,
span = 0.9,
method="loess",
lwd = 1) +
geom_hline(yintercept = Kistner.mod.aerial$summary.fixed[1,1],
linetype = "solid",
col = "darkgray",
size = 0.5) +
scale_linetype_manual(name = "Age Group",
values=c("solid", "longdash", "dotted")) +
scale_colour_manual(name = "Age Group",
values=c("red", "brown", "blue")) +
ylab("Kistner") +
xlab(expression(Deer~km^2)) +
ggtitle("Female (COMRECPOP)") +
theme_classic() +
scale_y_continuous(breaks = seq(50, 60, 5), limits = c(50, 60)) +
scale_x_continuous(breaks = seq(8, 40, 5), limits =c(8,40)) +
theme(axis.text=element_text(size=16),
legend.title = element_text(size=16, face="bold"),
legend.position = c(0.7,0.8),
legend.key = element_rect(fill = "white", linetype=0),
legend.background = element_rect(fill = "white", linetype=0),
plot.title = element_text(face="bold", size = 20, hjust=0.5),
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))
Total Kistner Density Dependence (Aerial Counts)
Horizontal line is global population mean Total Kistner.
mic.d = as.data.frame(Kistner.mod.aerial$summary.random$COMPRECNTDEN[,c(1:6)])
names(mic.d) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975")
length(unique(mic.d$ID))## [1] 17length(mic.d$ID)## [1] 102mic.d$Sex = rep(1:2, each = 51)
mic.d$ID2 = as.factor(
rep(c("Fawn (m)", "Yearling (m)", "Adult (m)",
"Fawn (f)", "Yearling (f)", "Adult (f)"),
each = length(unique(mic.d$ID))))
mic.d$Meane = mic.d$Mean + Kistner.mod.aerial$summary.fixed[1,1]
mic.d$Q025e = mic.d$Q025 + Kistner.mod.aerial$summary.fixed[1,1]
mic.d$Q975e = mic.d$Q975 + Kistner.mod.aerial$summary.fixed[1,1]
#mic.d$ID = mic.d$ID/30.63 #Area of Complex 1
range(mic.d$ID)## [1] 4.3 31.0male.set = mic.d %>% filter(Sex == 1)
range(male.set$Meane)## [1] 50.37919 63.93442female.set = mic.d %>% filter(Sex == 2)
range(female.set$Meane)## [1] 50.69919 57.30351ggplot(male.set, aes(ID, Meane, group = ID2)) +
geom_point(shape=1) +
geom_smooth(data= male.set,
aes(ID, Meane,
group = male.set$ID2,
col = factor(male.set$ID2),
linetype= factor(male.set$ID2)),
se=FALSE,
span = 0.9,
method="loess",
lwd = 1) +
geom_hline(yintercept = Kistner.mod.aerial$summary.fixed[1,1],
linetype = "solid",
col = "darkgray",
size = 0.5) +
scale_linetype_manual(name = "Age Group",
values=c("solid", "longdash", "dotted")) +
scale_colour_manual(name = "Age Group",
values=c("red", "brown", "blue")) +
ylab("Kistner (kg)") +
xlab(expression(Deer~km^2)) +
ggtitle("Male (COMPRECNTDEN)") +
theme_classic() +
scale_y_continuous(breaks = seq(50, 65, 5), limits = c(50, 65)) +
scale_x_continuous(breaks = seq(4, 35, 5), limits =c(4,35)) +
theme(axis.text=element_text(size=16),
legend.title = element_text(size=16, face="bold"),
legend.position = c(0.2,0.8),
legend.key = element_rect(fill = "white", linetype=0),
legend.background = element_rect(fill = "white", linetype=0),
plot.title = element_text(face="bold", size = 20, hjust=0.5),
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))
ggplot(female.set, aes(ID, Meane, group = ID2)) +
geom_point(shape=1) +
geom_smooth(data= female.set,
aes(ID, Meane,
group = female.set$ID2,
col = factor(female.set$ID2),
linetype= factor(female.set$ID2)),
se=FALSE,
span = 0.9,
method="loess",
lwd = 1) +
geom_hline(yintercept = Kistner.mod.aerial$summary.fixed[1,1],
linetype = "solid",
col = "darkgray",
size = 0.5) +
scale_linetype_manual(name = "Age Group",
values=c("solid", "longdash", "dotted")) +
scale_colour_manual(name = "Age Group",
values=c("red", "brown", "blue")) +
ylab("Kistner") +
xlab(expression(Deer~km^2)) +
ggtitle("Female (COMPRECNTDEN)") +
theme_classic() +
scale_y_continuous(breaks = seq(50, 60, 5), limits = c(50, 60)) +
scale_x_continuous(breaks = seq(4, 35, 5), limits =c(4,35)) +
theme(axis.text=element_text(size=16),
legend.title = element_text(size=16, face="bold"),
legend.position = c(0.7,0.8),
legend.key = element_rect(fill = "white", linetype=0),
legend.background = element_rect(fill = "white", linetype=0),
plot.title = element_text(face="bold", size = 20, hjust=0.5),
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))
Model Fitted Estimates.
Temp.data = FE.df1
idat = inla.stack.index(Kistner.stk, "Kist.0")$data #Indexregion data from model
Temp.data$Pred.Kist = Kistner.mod$summary.fitted.values$mean[idat] #Fitted values for region
Temp.data$Low = Kistner.mod$summary.fitted.values$`0.025quant`[idat]
Temp.data$High = Kistner.mod$summary.fitted.values$`0.975quant`[idat]
LW.fit = as.data.frame(
Temp.data %>%
group_by(SexAgeGroup.f) %>%
summarise(Mean = mean(Pred.Kist, na.rm=T),
Low = mean(Low, na.rm=T),
High = mean(High, na.rm=T)))
LW.fit$Age = c("Fawn (m)", "Yearling (m)", "Adult (m)",
"Fawn (f)", "Yearling (f)", "Adult (f)")
LW.fit$Age = ordered(factor(LW.fit$Age),
levels = c("Fawn (m)", "Yearling (m)", "Adult (m)",
"Fawn (f)", "Yearling (f)", "Adult (f)"))
ggplot(LW.fit, aes(factor(Age), Mean)) +
geom_bar(stat = "identity", fill = "lightgray") +
geom_point(size=2, pch=19, col = "red") +
geom_linerange(aes(ymin=Low, ymax=High), colour="black") +
theme_classic() +
ylab("Kistner (TOTAL)") +
xlab( "") +
ggtitle("Estimated Total Kistner \n(Complex 1, COMRECPOP)") +
theme(axis.text=element_text(size=16),
legend.title = element_text(size=16, face="bold"),
legend.position = c(0.7,0.7),
legend.key = element_rect(fill = "white", linetype=0),
legend.background = element_rect(fill = "white", linetype=0),
plot.title = element_text(face="bold", size = 20, hjust=0.5),
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))
Temp.data$Class = ifelse(Temp.data$SexAgeGeoup == "21", "Fawn (f)",
ifelse(Temp.data$SexAgeGeoup == "22", "Yearling (f)",
ifelse(Temp.data$SexAgeGeoup == "23", "Adult (f)",
ifelse(Temp.data$SexAgeGeoup == "11", "Fawn (m)",
ifelse(Temp.data$SexAgeGeoup == "12", "Yearling (m)",
ifelse(Temp.data$SexAgeGeoup == "13", "Adult (m)", NA))))))
Temp.data$Class = ordered(factor(Temp.data$Class,
levels = c("Fawn (m)", "Yearling (m)", "Adult (m)",
"Fawn (f)", "Yearling (f)", "Adult (f)")))
Temp.data$Density = Temp.data$COMRECPOP/30.63
COMRECPOP.plt = Temp.data %>% select(Class, Pred.Kist)
hline_dat1 = as.data.frame(
COMRECPOP.plt %>%
group_by(Class) %>%
summarise(Med = mean(Pred.Kist)))
ggplot(Temp.data, aes(Density, Pred.Kist)) +
geom_point(shape=1) +
geom_smooth(method="loess",
se=FALSE,
col="black",
lwd = 1) +
facet_wrap(~Class) +
geom_hline(data=hline_dat1,
aes(yintercept=Med),
col = "darkgray",
size = 0.5) +
ylab("Kistner (TOTAL)") +
xlab(expression(Deer~km^2)) +
ggtitle("Estimated Total Kistner \n(Complex 1, COMRECPOP)") +
theme_classic() +
theme(axis.text=element_text(size=16),
plot.title = element_text(face="bold", size = 20, hjust=0.5),
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))
Sex-Age Class Effect
Horizontal line is population average.
mic.d = as.data.frame(Kistner.mod.aerial$summary.random$SexAgeGroup.f)
names(mic.d) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975", "mode", "kld")
mic.d$Mean = mic.d$Mean + Kistner.mod$summary.fixed[1,1]
mic.d$Q025 = mic.d$Q025 + Kistner.mod$summary.fixed[1,1]
mic.d$Q975 = mic.d$Q975 + Kistner.mod$summary.fixed[1,1]
range(mic.d$Mean)## [1] 39.03594 65.45139ggplot(mic.d, aes(factor(ID), y=Mean)) +
geom_point(size=2, pch=19, col = "red") +
geom_linerange(aes(ymin=Q025, ymax=Q975), colour="black") +
geom_hline(yintercept = Kistner.mod.aerial$summary.fixed[1,1],
linetype = "solid",
col = "darkgray",
size = 0.5) +
theme_classic() +
xlab(" ") +
ylab("Kistner") +
ggtitle("Sex-Age Group Effect") +
scale_x_discrete(labels = c("Fawn (m)", "Yearling (m)", "Adult (m)",
"Fawn (f)", "Yearling (f)", "Adult (f)")) +
scale_y_continuous(breaks = seq(35, 70, 5), limits = c(35, 70)) +
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. Horizontal line is population mean.
mic.d = as.data.frame(Kistner.mod.aerial$summary.random$Month)
names(mic.d) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975", "mode", "kld")
mic.d$Mean = mic.d$Mean + Kistner.mod$summary.fixed[1,1]
mic.d$Q025 = mic.d$Q025 + Kistner.mod$summary.fixed[1,1]
mic.d$Q975 = mic.d$Q975 + Kistner.mod$summary.fixed[1,1]
range(mic.d$Mean)## [1] 34.43344 66.72099ggplot(mic.d, aes(factor(ID), y=Mean)) +
geom_point(size=2, pch=19, col = "red") +
geom_linerange(aes(ymin=Q025, ymax=Q975), colour="black") +
geom_hline(yintercept = Kistner.mod.aerial$summary.fixed[1,1],
linetype = "solid",
col = "darkgray",
size = 0.5) +
theme_classic() +
xlab(" ") +
ylab("Kistner") +
ggtitle("Month Effect") +
scale_x_discrete(labels = c("Jan", "Feb", "March",
"Oct", "Nov", "Dec")) +
scale_y_continuous(breaks = seq(30, 70, 5), limits = c(30, 70)) +
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(Kistner.mod.aerial$summary.random$Int.DATE_YR[,c(1:6)])
names(mic.d) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975")
mic.d$Mean = mic.d$Mean + Kistner.mod.aerial$summary.fixed[1,1]
mic.d$Q025 = mic.d$Q025 + Kistner.mod.aerial$summary.fixed[1,1]
mic.d$Q975 = mic.d$Q975 + Kistner.mod.aerial$summary.fixed[1,1]
range(mic.d$Mean)## [1] 50.45915 66.31360range(mic.d$Q975)## [1] 53.28662 69.16879range(mic.d$Q025)## [1] 47.62932 63.45603mic.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 = Kistner.mod.aerial$summary.fixed[1,1],
linetype = "solid",
col = "darkgray",
size = 0.5) +
scale_y_continuous(breaks = seq(45, 70, 5), limits = c(45, 70)) +
ylab("Kistner") +
xlab(" ") +
scale_x_continuous(breaks = 4:20,labels=paste0(mic.d$YLab)) +
ggtitle("Year 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 Horizontal line is population mean.
mic.d = as.data.frame(Kistner.mod.aerial$summary.random$PREG)
names(mic.d) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975", "mode", "kld")
mic.d$Mean = mic.d$Mean + Kistner.mod.aerial$summary.fixed[1,1]
mic.d$Q025 = mic.d$Q025 + Kistner.mod.aerial$summary.fixed[1,1]
mic.d$Q975 = mic.d$Q975 + Kistner.mod.aerial$summary.fixed[1,1]
range(mic.d$Mean)## [1] 53.02354 58.06838range(mic.d$Q975)## [1] 53.92174 58.96658range(mic.d$Q025)## [1] 52.12459 57.16943ggplot(mic.d, aes(factor(ID), y=Mean)) +
geom_point(size=2, pch=19, col = "red") +
geom_linerange(aes(ymin=Q025, ymax=Q975), colour="black") +
geom_hline(yintercept = Kistner.mod.aerial$summary.fixed[1,1],
linetype = "solid",
col = "darkgray",
size = 0.5) +
scale_y_continuous(breaks = seq(50, 60, 5), limits = c(50, 60)) +
theme_classic() +
xlab(" ") +
ylab("Kistner") +
ggtitle("Doe Pregnancy Effect") +
scale_x_discrete(labels = c("Yes", "No")) +
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))
2.3 Individual Kistners
Prepare Data for Models
Formulas for Indivdual KISTNER Loop (Changing density measure)
Formula.set = list()
Formula.set[["Kistner.A"]] = Kistner ~ -1 + intercept1 +
f(COMRECPOP, #Recon Pop (total)
model="rw1",
constr = TRUE,
scale.model = TRUE,
replicate = RepGroups,
hyper = nPrior) +
f(myID,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(SexAgeGroup.f,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(PREG,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(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 ) +
DSubZero + LiveWtKg
Formula.set[["Kistner.B"]] = Kistner ~ -1 + intercept1 +
f(COMRECDOEPOP, #Recon doe pop
model="rw1",
constr = TRUE,
scale.model = TRUE,
replicate = RepGroups,
hyper = nPrior) +
f(myID,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(SexAgeGroup.f,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(PREG,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(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 ) +
DSubZero + LiveWtKg
Formula.set[["Kistner.C"]] = Kistner ~ -1 + intercept1 +
f(COMPRECNTDEN, #Aerial
model="rw1",
constr = TRUE,
scale.model = TRUE,
replicate = RepGroups,
hyper = nPrior) +
f(myID,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(SexAgeGroup.f,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(PREG,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(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 ) +
DSubZero + LiveWtKgModel Loop for Individual Kistner Scores
theta1 = c(-3.949758515, -0.668447978, 0.023023690, -2.971971613, -0.026000289, 0.008109743,
-1.073993338, -0.302061119, 0.039326732)
Ind.Kistner = c("BRISK","OMEN", "KDNY", "PERCDM", "HEART", "RUMP", "MUSCLE")
Density.metrics = c("COMRECPOP", "COMRECDOEPOP", "COMPRECNTDEN")
for(i in 1:length(Density.metrics)){
Loop.metric = Density.metrics[i]
for(j in 1:length(Ind.Kistner)){
Loop.score = Ind.Kistner[j]
FE.df1$Ind.Kist = FE.df1[,Loop.score]
IndKist.stk = inla.stack(data = list(Kistner = FE.df1$Ind.Kist),
A = list(1,1),
effects = FE.lst1,
tag = "iK.0")
IndKist.mod = inla(Formula.set[[i]],
data = inla.stack.data(IndKist.stk),
family = "gaussian",
verbose = TRUE,
control.fixed = list(prec = 0.001,
prec.intercept = 0.0001),
control.predictor = list(
A = inla.stack.A(IndKist.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))
mic.d = as.data.frame(IndKist.mod$summary.random[[1]][,c(1:6)])
names(mic.d) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975")
mic.d$Score = Loop.score
mic.d$Metric = Loop.metric
mic.d$Int = IndKist.mod$summary.fixed[1,1]
mic.d$DIC = IndKist.mod$dic$dic
mic.d$WAIC = IndKist.mod$waic$waic
if(j == 1){Ind.Scores.df = mic.d}
else{Ind.Scores.df = rbind(Ind.Scores.df, mic.d)}
}
if(i == 1){Ind.Scores.m.df = Ind.Scores.df}
else{Ind.Scores.m.df = rbind(Ind.Scores.m.df, Ind.Scores.df)}
}Individual Kistner Scores
IndKist.plt = Ind.Scores.m.df
IndKist.plt$ID = ifelse(IndKist.plt$Metric == "COMRECPOP" | IndKist.plt$Metric == "COMRECDOEPOP", IndKist.plt$ID/30.63, IndKist.plt$ID)
IndKist.plt$Mean = IndKist.plt$Mean + IndKist.plt$Int
IndKist.plt$Q025 = IndKist.plt$Q025 + IndKist.plt$Int
IndKist.plt$Q975 = IndKist.plt$Q975 + IndKist.plt$Int
range(IndKist.plt$Mean)## [1] 4.394922 10.678574range(IndKist.plt$Q975)## [1] 5.268521 11.878574range(IndKist.plt$Q025)## [1] 3.114421 9.547753COMRECPOP.plt = IndKist.plt %>% filter(Metric == "COMRECPOP")
COMRECDOEPOP.plt = IndKist.plt %>% filter(Metric == "COMRECDOEPOP")
COMPRECNTDEN.plt = IndKist.plt %>% filter(Metric == "COMPRECNTDEN")
hline_dat1 = as.data.frame(
COMRECPOP.plt %>%
group_by(Score) %>%
summarise(Med = median(Int)))
hline_dat2 = as.data.frame(
COMRECDOEPOP.plt %>%
group_by(Score) %>%
summarise(Med = median(Int)))
hline_dat3 = as.data.frame(
COMPRECNTDEN.plt %>%
group_by(Score) %>%
summarise(Med = median(Int)))
ggplot(COMRECPOP.plt, aes(ID, Mean)) +
geom_smooth(col = "black",
linetype= "solid",
method = "loess",
span = 0.5,
se = FALSE,
lwd = 1) +
geom_smooth(data = COMRECPOP.plt, aes(ID, Q975),
col = "grey",
method = "loess",
span = 0.5,
se = FALSE,
linetype= "dashed") +
geom_smooth(data = COMRECPOP.plt, aes(ID, Q025),
col = "grey",
method = "loess",
span = 0.5,
se = FALSE,
linetype= "dashed") +
xlab(expression(Deer~km^2)) +
ylab("Kistner Score") +
facet_wrap(~Score, ncol = 2, scales="free_y") +
geom_hline(data=hline_dat1,
aes(yintercept=Med),
col = "darkgray",
size = 0.5) +
ggtitle("Individual Kistner Scores (COMRECPOP)") +
theme_classic() +
scale_x_continuous(breaks = seq(4, 40, 10), limits =c(4,40)) +
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 = 14),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
strip.background = element_blank(),
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))
ggplot(COMRECDOEPOP.plt, aes(ID, Mean)) +
geom_smooth(col = "black",
linetype= "solid",
method = "loess",
span = 0.5,
se = FALSE,
lwd = 1) +
geom_smooth(data = COMRECDOEPOP.plt, aes(ID, Q975),
col = "grey",
method = "loess",
span = 0.5,
se = FALSE,
linetype= "dashed") +
geom_smooth(data = COMRECDOEPOP.plt, aes(ID, Q025),
col = "grey",
method = "loess",
span = 0.5,
se = FALSE,
linetype= "dashed") +
xlab(expression(Deer~km^2)) +
ylab("Kistner Score") +
facet_wrap(~Score, ncol = 2, scales="free_y") +
geom_hline(data=hline_dat2,
aes(yintercept=Med),
col = "darkgray",
size = 0.5) +
ggtitle("Individual Kistner Scores (COMRECDOEPOP)") +
theme_classic() +
scale_x_continuous(breaks = seq(4, 40, 10), limits =c(4,40)) +
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 = 14),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
strip.background = element_blank(),
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))
ggplot(COMPRECNTDEN.plt, aes(ID, Mean)) +
geom_smooth(col = "black",
linetype= "solid",
method = "loess",
span = 0.5,
se = FALSE,
lwd = 1) +
geom_smooth(data = COMPRECNTDEN.plt, aes(ID, Q975),
col = "grey",
method = "loess",
span = 0.5,
se = FALSE,
linetype= "dashed") +
geom_smooth(data = COMPRECNTDEN.plt, aes(ID, Q025),
col = "grey",
method = "loess",
span = 0.5,
se = FALSE,
linetype= "dashed") +
xlab(expression(Deer~km^2)) +
ylab("Kistner Score") +
facet_wrap(~Score, ncol = 2, scales="free_y") +
geom_hline(data=hline_dat3,
aes(yintercept=Med),
col = "darkgray",
size = 0.5) +
ggtitle("Individual Kistner Scores (COMPRECNTDEN)") +
theme_classic() +
scale_x_continuous(breaks = seq(4, 40, 10), limits =c(4,40)) +
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 = 14),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
strip.background = element_blank(),
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))
2.4 Pregnancy
2.4.1 All Age Groups
Organize Doe Data as List Object
Doe.df = Deer.Models.Data %>% filter(COMPLEX == 1 & SEX == 2) #Complex 1 and females only
#Yes = 1, No = 2
Doe.df$PREG2 = ifelse(Doe.df$PREG == 1, 1,
ifelse(Doe.df$PREG == 2, 0, NA))
#Missing value: Was pregnancy detectable, if yes then "0", else NA
Doe.df$PREG2 = ifelse(Doe.df$USEREPRO == 1 & is.na(Doe.df$PREG2) == TRUE, 0, Doe.df$PREG2)
#Check totals
length(which(is.na(Doe.df$PREG2)))## [1] 495length(which(Doe.df$PREG2 == 0))## [1] 734length(which(Doe.df$PREG2 == 1))## [1] 1013#Missing value: Was pregnancy detectable, if yes then "0", else NA
Doe.df$FETUS2 = ifelse(Doe.df$USEREPRO == 1 & is.na(Doe.df$FETUS) == TRUE, 0, Doe.df$FETUS)
length(which(is.na(Doe.df$FETUS)))## [1] 1384#COMCOMCNT = "Survey completed" If not, set record to NA
Doe.df$COMPRECNTDEN2 = ifelse(Doe.df$COMCOMCNT == 2, NA, Doe.df$COMPRECNTDEN) #Complete aerial survey
Doe.df = Doe.df %>%
filter(is.na(AGEGROUP) == FALSE) #Remove records w/ missing age
#Integer version of Sex-Age groups
unique(Doe.df$SexAgeGroup.f)## [1] 21 22 23
## Levels: 11 12 13 21 22 23Doe.df$RepGroups = as.integer(Doe.df$SexAgeGroup.f)
unique(Doe.df$RepGroups)## [1] 4 5 6Doe.df$RepGroups[Doe.df$RepGroups == 4] = 1
Doe.df$RepGroups[Doe.df$RepGroups == 5] = 2
Doe.df$RepGroups[Doe.df$RepGroups == 6] = 3
unique(Doe.df$RepGroups)## [1] 1 2 3Doe.df$SexAgeGroup.f = factor(as.character(Doe.df$SexAgeGroup.f))
XXX = Doe.df %>% select(SexAgeGroup.f, RepGroups)
#Unique ID by Deer
Doe.df$myID = 1:nrow(Doe.df)
#set up model data as a list object
## Not all covariates are used.
FE.lst2 = list(list(intercept1 = rep(1, dim(Doe.df)[1])),
list(COMPREHARKM = Doe.df[,"COMPREHARKM"],
COMRECPOP = Doe.df[,"COMRECPOP"],
COMRECDOEPOP = Doe.df[,"COMRECDOEPOP"],
COMCUMPRE2HARKM = Doe.df[,"COMCUMPRE2HARKM"],
TOTAL = Doe.df[,"TOTAL"],
LiveWtKg = Doe.df[,"LiveWtKg"],
Lat = Doe.df[,"Lat"],
NN = Doe.df[,"NN"],
Area = Doe.df[,"CmpxArea"],
HINDFOOT = Doe.df[,"HINDFOOT"]/100,
COMPREDOEHARKM = Doe.df[,"COMPREDOEHARKM"],
COMCUMPRE2DOEHARKM = Doe.df[,"COMCUMPRE2DOEHARKM"],
COMPRECNTDEN = Doe.df[,"COMPRECNTDEN"],
GrowSeas = round(Doe.df[,"GrowSeas"],3),
GrowSeas.sc = Doe.df[,"GrowSeas.sc"],
PrecEx10mm = Doe.df[,"PrecEx10mm"],
PrecEx20mm = Doe.df[,"PrecEx20mm"],
DSubZero = Doe.df[,"DSubZero"],
Prec95pct = Doe.df[,"Prec95pct"],
MonMxDT = Doe.df[,"MonMxDT"],
MMx1DPrec = Doe.df[,"MMx1DPrec"],
AMinT = Doe.df[,"AMinT"],
AMxT = Doe.df[,"AMxT"],
SumDur = round(Doe.df[,"SumDur"],3),
SumDur.sc = Doe.df[,"SumDur.sc"],
FrstDays = Doe.df[,"FrstDays"],
Int.DATE_YR = Doe.df[,"Int.DATE_YR"],
RepGroups = Doe.df[,"RepGroups"],
SexAgeGroup.f = Doe.df[,"SexAgeGroup.f"],
AGEGROUP = Doe.df[,"AGEGROUP"],
PREG = as.factor(Doe.df[,"PREG"]),
f_COMPLEX = Doe.df[,"f_COMPLEX"],
myID = Doe.df[,"myID"],
Month = Doe.df[,"DATE_MO"]))
#Data for Combined Adult + Yearling Model
Doe.ya.df = Doe.df %>% filter(SexAgeGroup.f != "21")
Doe.ya.df$SexAgeGroup.f = factor(as.character(Doe.ya.df$SexAgeGroup.f))
#Unique ID by Deer
Doe.ya.df$myID = 1:nrow(Doe.ya.df)
FE.lst2B = list(list(intercept1 = rep(1, dim(Doe.ya.df)[1])),
list(COMPREHARKM = Doe.ya.df[,"COMPREHARKM"],
COMRECPOP = Doe.ya.df[,"COMRECPOP"],
COMRECDOEPOP = Doe.ya.df[,"COMRECDOEPOP"],
COMCUMPRE2HARKM = Doe.ya.df[,"COMCUMPRE2HARKM"],
TOTAL = Doe.ya.df[,"TOTAL"],
LiveWtKg = Doe.ya.df[,"LiveWtKg"],
Lat = Doe.ya.df[,"Lat"],
NN = Doe.ya.df[,"NN"],
Area = Doe.ya.df[,"CmpxArea"],
HINDFOOT = Doe.ya.df[,"HINDFOOT"]/100,
COMPREDOEHARKM = Doe.ya.df[,"COMPREDOEHARKM"],
COMCUMPRE2DOEHARKM = Doe.ya.df[,"COMCUMPRE2DOEHARKM"],
COMPRECNTDEN = Doe.ya.df[,"COMPRECNTDEN"],
GrowSeas = round(Doe.ya.df[,"GrowSeas"],3),
GrowSeas.sc = Doe.ya.df[,"GrowSeas.sc"],
PrecEx10mm = Doe.ya.df[,"PrecEx10mm"],
PrecEx20mm = Doe.ya.df[,"PrecEx20mm"],
DSubZero = Doe.ya.df[,"DSubZero"],
Prec95pct = Doe.ya.df[,"Prec95pct"],
MonMxDT = Doe.ya.df[,"MonMxDT"],
MMx1DPrec = Doe.ya.df[,"MMx1DPrec"],
AMinT = Doe.ya.df[,"AMinT"],
AMxT = Doe.ya.df[,"AMxT"],
SumDur = round(Doe.ya.df[,"SumDur"],3),
SumDur.sc = Doe.ya.df[,"SumDur.sc"],
FrstDays = Doe.ya.df[,"FrstDays"],
Int.DATE_YR = Doe.ya.df[,"Int.DATE_YR"],
RepGroups = Doe.ya.df[,"RepGroups"],
SexAgeGroup.f = Doe.ya.df[,"SexAgeGroup.f"],
AGEGROUP = Doe.ya.df[,"AGEGROUP"],
PREG = as.factor(Doe.ya.df[,"PREG"]),
f_COMPLEX = Doe.ya.df[,"f_COMPLEX"],
myID = Doe.ya.df[,"myID"],
Month = Doe.ya.df[,"DATE_MO"]))
#All Ages
Preg.stk = inla.stack(data = list(Preg = Doe.df$PREG2),
A = list(1,1),
effects = FE.lst2,
tag = "base.0")
#Yearling + Adult Combination
Preg2.stk = inla.stack(data = list(Preg = Doe.ya.df$PREG2),
A = list(1,1),
effects = FE.lst2B,
tag = "ya.0")Model for Pregnancy (Reconstructed Population COMRECPOP)
Pregnancy.frm = Preg ~ -1 + intercept1 +
f(COMRECPOP,
model="rw1",
constr = TRUE,
scale.model = TRUE,
replicate = RepGroups,
hyper = nPrior) +
f(myID,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(SexAgeGroup.f,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(Month,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(Int.DATE_YR, #year
model="rw1",
constr = TRUE,
scale.model = TRUE,
hyper = nPrior ) +
f(SumDur,#summer duration
model="rw1",
constr = TRUE,
scale.model = TRUE,
hyper = nPrior ) +
DSubZero + TOTAL
#theta.m1 = Pregnancy.mod$internal.summary.hyperpar$mean
theta.m1 = c(-0.2969193, 0.6452064, -0.5640170, -0.6115045, -0.2670323, -0.0485563)
Pregnancy.mod = inla(Pregnancy.frm,
data = inla.stack.data(Preg.stk),
family = "binomial",
verbose = TRUE,
control.fixed = list(prec = 0.001,
prec.intercept = 0.0001),
control.predictor = list(
A = inla.stack.A(Preg.stk),
compute = TRUE,
link = 1),
control.mode = list(restart = TRUE, theta = theta.m1),
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(Pregnancy.mod)##
## Call:
## c("inla(formula = Pregnancy.frm, family = \"binomial\", data = inla.stack.data(Preg.stk), ", " verbose = TRUE, control.compute = list(dic = TRUE, cpo = TRUE, ", " waic = TRUE), control.predictor = list(A = inla.stack.A(Preg.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 = theta.m1))")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 1.8640 6.5815 0.2872 8.7327
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept1 -2.7313 0.3754 -3.4683 -2.7313 -1.9949 -2.7313 0
## DSubZero 0.0731 0.2462 -0.4102 0.0731 0.5561 0.0731 0
## TOTAL 0.0359 0.0038 0.0285 0.0359 0.0434 0.0359 0
##
## Random effects:
## Name Model
## COMRECPOP RW1 model
## myID IID model
## SexAgeGroup.f IID model
## Month IID model
## Int.DATE_YR RW1 model
## SumDur RW1 model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant 0.975quant
## Precision for COMRECPOP 0.8039 0.2481 0.4240 0.7683 1.3909
## Precision for myID 1.9499 0.4843 1.1736 1.8911 3.0645
## Precision for SexAgeGroup.f 0.5950 0.1518 0.3514 0.5768 0.9448
## Precision for Month 0.5530 0.1376 0.3306 0.5369 0.8688
## Precision for Int.DATE_YR 0.8303 0.2645 0.4299 0.7909 1.4592
## Precision for SumDur 0.7754 0.2280 0.4223 0.7439 1.3112
## mode
## Precision for COMRECPOP 0.7020
## Precision for myID 1.7793
## Precision for SexAgeGroup.f 0.5422
## Precision for Month 0.5064
## Precision for Int.DATE_YR 0.7181
## Precision for SumDur 0.6850
##
## Expected number of effective parameters(std dev): 91.88(0.00)
## Number of equivalent replicates : 19.01
##
## Deviance Information Criterion (DIC) ...............: 853.69
## Deviance Information Criterion (DIC, saturated) ....: 853.69
## Effective number of parameters .....................: 94.84
##
## Watanabe-Akaike information criterion (WAIC) ...: 837.85
## Effective number of parameters .................: 71.98
##
## Marginal log-Likelihood: -536.79
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computedModel for Pregnancy (Reconstructed Doe Population COMRECDOEPOP)
Pregnancy.rd.frm = Preg ~ -1 + intercept1 +
f(COMRECDOEPOP,
model="rw1",
constr = TRUE,
scale.model = TRUE,
replicate = RepGroups,
hyper = nPrior) +
f(myID,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(SexAgeGroup.f,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(Month,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(Int.DATE_YR, #year
model="rw1",
constr = TRUE,
scale.model = TRUE,
hyper = nPrior ) +
f(SumDur,#summer duration
model="rw1",
constr = TRUE,
scale.model = TRUE,
hyper = nPrior ) +
DSubZero + TOTAL
#theta.m1 = Pregnancy.rd.mod$internal.summary.hyperpar$mean
theta.m1 = c(-0.2969193, 0.6452064, -0.5640170, -0.6115045, -0.2670323, -0.0485563)
Pregnancy.rd.mod = inla(Pregnancy.rd.frm,
data = inla.stack.data(Preg.stk),
family = "binomial",
verbose = TRUE,
control.fixed = list(prec = 0.001,
prec.intercept = 0.0001),
control.predictor = list(
A = inla.stack.A(Preg.stk),
compute = TRUE,
link = 1),
control.mode = list(restart = TRUE, theta = theta.m1),
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(Pregnancy.rd.mod)##
## Call:
## c("inla(formula = Pregnancy.rd.frm, family = \"binomial\", data = inla.stack.data(Preg.stk), ", " verbose = TRUE, control.compute = list(dic = TRUE, cpo = TRUE, ", " waic = TRUE), control.predictor = list(A = inla.stack.A(Preg.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 = theta.m1))")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 2.0364 6.5582 0.2768 8.8714
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept1 -2.6972 0.3731 -3.4297 -2.6972 -1.9653 -2.6972 0
## DSubZero 0.0643 0.2428 -0.4125 0.0643 0.5407 0.0643 0
## TOTAL 0.0354 0.0038 0.0280 0.0354 0.0428 0.0354 0
##
## Random effects:
## Name Model
## COMRECDOEPOP RW1 model
## myID IID model
## SexAgeGroup.f IID model
## Month IID model
## Int.DATE_YR RW1 model
## SumDur RW1 model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant 0.975quant
## Precision for COMRECDOEPOP 0.8267 0.2574 0.4352 0.7890 1.4370
## Precision for myID 1.9287 0.4801 1.1597 1.8704 3.0342
## Precision for SexAgeGroup.f 0.6089 0.1565 0.3578 0.5900 0.9697
## Precision for Month 0.5547 0.1379 0.3316 0.5386 0.8712
## Precision for Int.DATE_YR 0.8290 0.2615 0.4301 0.7911 1.4490
## Precision for SumDur 0.8230 0.2449 0.4447 0.7889 1.3996
## mode
## Precision for COMRECDOEPOP 0.7190
## Precision for myID 1.7595
## Precision for SexAgeGroup.f 0.5543
## Precision for Month 0.5082
## Precision for Int.DATE_YR 0.7205
## Precision for SumDur 0.7250
##
## Expected number of effective parameters(std dev): 92.58(0.00)
## Number of equivalent replicates : 18.87
##
## Deviance Information Criterion (DIC) ...............: 858.97
## Deviance Information Criterion (DIC, saturated) ....: 858.97
## Effective number of parameters .....................: 95.67
##
## Watanabe-Akaike information criterion (WAIC) ...: 843.23
## Effective number of parameters .................: 72.86
##
## Marginal log-Likelihood: -540.43
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computedModel for Pregnancy (Aerial Counts COMPRECNTDEN)
Pregnancy.ac.frm = Preg ~ -1 + intercept1 +
f(COMPRECNTDEN,
model="rw1",
constr = TRUE,
scale.model = TRUE,
replicate = RepGroups,
hyper = nPrior) +
f(myID,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(SexAgeGroup.f,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(Month,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(Int.DATE_YR, #year
model="rw1",
constr = TRUE,
scale.model = TRUE,
hyper = nPrior ) +
f(SumDur,#summer duration
model="rw1",
constr = TRUE,
scale.model = TRUE,
hyper = nPrior ) +
DSubZero + TOTAL
#theta.m1 = Pregnancy.ac.mod$internal.summary.hyperpar$mean
theta.m1 = c(-0.06524442, 0.62028369, -0.52218859, -0.58439412, -0.35774291, -0.05702578)
Pregnancy.ac.mod = inla(Pregnancy.ac.frm,
data = inla.stack.data(Preg.stk),
family = "binomial",
verbose = TRUE,
control.fixed = list(prec = 0.001,
prec.intercept = 0.0001),
control.predictor = list(
A = inla.stack.A(Preg.stk),
compute = TRUE,
link = 1),
control.mode = list(restart = TRUE, theta = theta.m1),
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(Pregnancy.ac.mod)##
## Call:
## c("inla(formula = Pregnancy.ac.frm, family = \"binomial\", data = inla.stack.data(Preg.stk), ", " verbose = TRUE, control.compute = list(dic = TRUE, cpo = TRUE, ", " waic = TRUE), control.predictor = list(A = inla.stack.A(Preg.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 = theta.m1))")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 2.1124 6.1961 0.3092 8.6176
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept1 -2.6820 0.3705 -3.4094 -2.6820 -1.9552 -2.6820 0
## DSubZero 0.0906 0.2446 -0.3896 0.0906 0.5704 0.0906 0
## TOTAL 0.0351 0.0037 0.0278 0.0351 0.0425 0.0351 0
##
## Random effects:
## Name Model
## COMPRECNTDEN RW1 model
## myID IID model
## SexAgeGroup.f IID model
## Month IID model
## Int.DATE_YR RW1 model
## SumDur RW1 model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant 0.975quant
## Precision for COMPRECNTDEN 1.0439 0.3324 0.5412 0.9945 1.8347
## Precision for myID 1.8961 0.4749 1.1369 1.8381 2.9911
## Precision for SexAgeGroup.f 0.6201 0.1600 0.3636 0.6008 0.9893
## Precision for Month 0.5630 0.1402 0.3361 0.5467 0.8847
## Precision for Int.DATE_YR 0.7291 0.2195 0.3914 0.6981 1.2465
## Precision for SumDur 0.8600 0.2562 0.4648 0.8241 1.4630
## mode
## Precision for COMPRECNTDEN 0.9031
## Precision for myID 1.7279
## Precision for SexAgeGroup.f 0.5641
## Precision for Month 0.5158
## Precision for Int.DATE_YR 0.6402
## Precision for SumDur 0.7570
##
## Expected number of effective parameters(std dev): 93.82(0.00)
## Number of equivalent replicates : 18.62
##
## Deviance Information Criterion (DIC) ...............: 865.29
## Deviance Information Criterion (DIC, saturated) ....: 865.29
## Effective number of parameters .....................: 96.76
##
## Watanabe-Akaike information criterion (WAIC) ...: 850.47
## Effective number of parameters .................: 74.62
##
## Marginal log-Likelihood: -558.41
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computedCompare DIC and WAIC
#DIC
Pregnancy.mod$dic$dic## [1] 853.692Pregnancy.rd.mod$dic$dic## [1] 858.9729Pregnancy.ac.mod$dic$dic## [1] 865.2875#WAIC
Pregnancy.mod$waic$waic## [1] 837.8519Pregnancy.rd.mod$waic$waic## [1] 843.2301Pregnancy.ac.mod$waic$waic## [1] 850.4666Pregnancy Density Dependence (Reconstructed Population, COMRECPOP)
mic.d = as.data.frame(Pregnancy.mod$summary.random$COMRECPOP[,c(1:6)])
names(mic.d) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975")
length(unique(mic.d$ID))## [1] 13length(mic.d$ID)## [1] 39mic.d$ID2 = as.factor(
rep(c("Fawn (f)", "Yearling (f)", "Adult (f)"),
each = length(unique(mic.d$ID))))
mic.d$Meane = plogis(mic.d$Mean + Pregnancy.mod$summary.fixed[1,1])
mic.d$Q025e = plogis(mic.d$Q025 + Pregnancy.mod$summary.fixed[1,1])
mic.d$Q975e = plogis(mic.d$Q975 + Pregnancy.mod$summary.fixed[1,1])
mic.d$ID = mic.d$ID/30.63 #Area of Complex 1
range(mic.d$ID)## [1] 8.651649 38.132550ggplot(mic.d, aes(ID, Meane, group = ID2)) +
geom_point(shape=1) +
geom_smooth(data= mic.d,
aes(ID, Meane,
group = mic.d$ID2,
col = factor(mic.d$ID2),
linetype= factor(mic.d$ID2)),
se=FALSE,
span = 0.8,
method="loess",
lwd = 1) +
geom_hline(yintercept = plogis(Pregnancy.mod$summary.fixed[1,1]),
linetype = "solid",
col = "darkgray",
size = 0.5) +
scale_linetype_manual(name = "Age Group",
values=c("solid", "longdash", "dotted")) +
scale_colour_manual(name = "Age Group",
values=c("red", "brown", "blue")) +
ylab("Pregnancy (probability)") +
xlab(expression(Deer~km^2)) +
ggtitle("Pregnancy (COMRECPOP)") +
theme_classic() +
scale_y_continuous(breaks = seq(0, 0.25, 0.05), limits = c(-0.1, 0.25)) +
scale_x_continuous(breaks = seq(8, 40, 5), limits =c(8,40)) +
theme(axis.text=element_text(size=16),
legend.title = element_text(size=16, face="bold"),
legend.position = c(0.8,0.8),
legend.key = element_rect(fill = "white", linetype=0),
legend.background = element_rect(fill = "white", linetype=0),
plot.title = element_text(face="bold", size = 20, hjust=0.5),
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))
Pregnancy Density Dependence (Aerial Counts, COMPRECNTDEN)
mic.d = as.data.frame(Pregnancy.ac.mod$summary.random$COMPRECNTDEN[,c(1:6)])
names(mic.d) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975")
length(unique(mic.d$ID))## [1] 17length(mic.d$ID)## [1] 51range(mic.d$ID)## [1] 4.3 31.0mic.d$ID2 = as.factor(
rep(c("Fawn (f)", "Yearling (f)", "Adult (f)"),
each = length(unique(mic.d$ID))))
mic.d$Meane = plogis(mic.d$Mean + Pregnancy.ac.mod$summary.fixed[1,1])
mic.d$Q025e = plogis(mic.d$Q025 + Pregnancy.ac.mod$summary.fixed[1,1])
mic.d$Q975e = plogis(mic.d$Q975 + Pregnancy.ac.mod$summary.fixed[1,1])
ggplot(mic.d, aes(ID, Meane, group = ID2)) +
geom_point(shape=1) +
geom_smooth(data= mic.d,
aes(ID, Meane,
group = mic.d$ID2,
col = factor(mic.d$ID2),
linetype= factor(mic.d$ID2)),
se=FALSE,
span = 0.8,
method="loess",
lwd = 1) +
geom_hline(yintercept = plogis(Pregnancy.mod$summary.fixed[1,1]),
linetype = "solid",
col = "darkgray",
size = 0.5) +
scale_linetype_manual(name = "Age Group",
values=c("solid", "longdash", "dotted")) +
scale_colour_manual(name = "Age Group",
values=c("red", "brown", "blue")) +
ylab("Pregnancy (probability)") +
xlab(expression(Deer~km^2)) +
ggtitle("Pregnancy (COMPRECNTDEN)") +
theme_classic() +
scale_x_continuous(breaks = seq(4, 35, 5), limits =c(4,35)) +
theme(axis.text=element_text(size=16),
legend.title = element_text(size=16, face="bold"),
legend.position = c(0.3,0.8),
legend.key = element_rect(fill = "white", linetype=0),
legend.background = element_rect(fill = "white", linetype=0),
plot.title = element_text(face="bold", size = 20, hjust=0.5),
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))
Model Fitted Estimates.
Temp.data = Doe.df
idat = inla.stack.index(Preg.stk, "base.0")$data #Indexregion data from model
Temp.data$Pred.Preg = Pregnancy.mod$summary.fitted.values$mean[idat] #Fitted values for region
Temp.data$Low = Pregnancy.mod$summary.fitted.values$`0.025quant`[idat]
Temp.data$High = Pregnancy.mod$summary.fitted.values$`0.975quant`[idat]
Preg.fit = as.data.frame(
Temp.data %>%
group_by(SexAgeGroup.f) %>%
summarise(Mean = mean(Pred.Preg, na.rm=T),
Low = mean(Low, na.rm=T),
High = mean(High, na.rm=T)))
Preg.fit$Age = c("Fawn (f)", "Yearling (f)", "Adult (f)")
Preg.fit$Age = ordered(factor(Preg.fit$Age), levels = c("Fawn (f)", "Yearling (f)", "Adult (f)"))
ggplot(Preg.fit, aes(factor(Age),Mean)) +
geom_bar(stat = "identity", fill = "lightgray") +
geom_point(size=2, pch=19, col = "red") +
geom_linerange(aes(ymin=Low, ymax=High), colour="black") +
theme_classic() +
ylab("Mean Pregnancy Probabilty") +
xlab( "") +
ggtitle("Population Mean Pregnancy (COMRECPOP)") +
scale_x_discrete(labels = c("Fawn (f)", "Yearling (f)", "Adult (f)")) +
scale_y_continuous(breaks = seq(0, 1, 0.25), limits = c(0, 1)) +
theme(axis.text=element_text(size=16),
legend.title = element_text(size=16, face="bold"),
legend.position = c(0.7,0.7),
legend.key = element_rect(fill = "white", linetype=0),
legend.background = element_rect(fill = "white", linetype=0),
plot.title = element_text(face="bold", size = 20, hjust=0.5),
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))
Temp.data$Class = ifelse(Temp.data$SexAgeGeoup == "21", "Fawn (f)",
ifelse(Temp.data$SexAgeGeoup == "22", "Yearling (f)",
ifelse(Temp.data$SexAgeGeoup == "23", "Adult (f)", NA)))
Temp.data$Class = ordered(factor(Temp.data$Class, levels = c("Fawn (f)", "Yearling (f)", "Adult (f)")))
Temp.data$Density = Temp.data$COMRECPOP/30.63
COMRECPOP.plt = Temp.data %>% select(Class, Pred.Preg)
hline_dat1 = as.data.frame(
COMRECPOP.plt %>%
group_by(Class) %>%
summarise(Med = mean(Pred.Preg)))
ggplot(Temp.data, aes(Density, Pred.Preg)) +
geom_point(shape=1) +
geom_smooth(method="loess",
se=FALSE,
col="black",
lwd = 1) +
facet_wrap(~Class) +
geom_hline(data=hline_dat1,
aes(yintercept=Med),
col = "darkgray",
size = 0.5) +
ylab("Pregnancy Probability") +
xlab(expression(Deer~km^2)) +
ggtitle("Estimated Pregnancy (Complex 1) \n COMRECPOP") +
theme_classic() +
theme(axis.text=element_text(size=16),
plot.title = element_text(face="bold", size = 20, hjust=0.5),
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))
Sex-Age Class Effect
Horizontal line is doe population average.
mic.d = as.data.frame(Pregnancy.mod$summary.random$SexAgeGroup.f)
names(mic.d) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975", "mode", "kld")
mic.d$Mean = plogis(mic.d$Mean + Pregnancy.mod$summary.fixed[1,1])
mic.d$Q025 = plogis(mic.d$Q025 + Pregnancy.mod$summary.fixed[1,1])
mic.d$Q975 = plogis(mic.d$Q975 + Pregnancy.mod$summary.fixed[1,1])
range(mic.d$Mean)## [1] 0.002109522 0.343072101ggplot(mic.d, aes(factor(ID), y=Mean)) +
geom_point(size=2, pch=19, col = "red") +
geom_linerange(aes(ymin=Q025, ymax=Q975), colour="black") +
geom_hline(yintercept = plogis(Pregnancy.mod$summary.fixed[1,1]),
linetype = "solid",
col = "darkgray",
size = 0.5) +
theme_classic() +
xlab(" ") +
ylab("Pregnancy (probability") +
ggtitle("Age Effect") +
scale_x_discrete(labels = c("Fawn (f)", "Yearling (f)", "Adult (f)")) +
scale_y_continuous(breaks = seq(0, 0.5, 0.05), limits = c(-0.05, 0.5)) +
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. Horizontal line is doe population mean.
mic.d = as.data.frame(Pregnancy.mod$summary.random$Month)
names(mic.d) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975", "mode", "kld")
mic.d$Mean = plogis(mic.d$Mean + Pregnancy.mod$summary.fixed[1,1])
mic.d$Q025 = plogis(mic.d$Q025 + Pregnancy.mod$summary.fixed[1,1])
mic.d$Q975 = plogis(mic.d$Q975 + Pregnancy.mod$summary.fixed[1,1])
range(mic.d$Mean)## [1] 0.002024723 0.362773247ggplot(mic.d, aes(factor(ID), y=Mean)) +
geom_point(size=2, pch=19, col = "red") +
geom_linerange(aes(ymin=Q025, ymax=Q975), colour="black") +
geom_hline(yintercept = plogis(Pregnancy.mod$summary.fixed[1,1]),
linetype = "solid",
col = "darkgray",
size = 0.5) +
theme_classic() +
xlab(" ") +
ylab("Pregnancy (probability)") +
ggtitle("Month Effect") +
scale_x_discrete(labels = c("Jan", "Feb", "March",
"Oct", "Nov", "Dec")) +
scale_y_continuous(breaks = seq(0, 0.6, 0.05), limits = c(-0.05, 0.6)) +
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))
2.4.2 Combined Yearling and Adult
Model for Yearling + Adult Pregnancy (Reconstructed Population COMRECPOP)
YA.Pregnancy.frm = Preg ~ -1 + intercept1 +
f(COMRECPOP,
model="rw1",
constr = TRUE,
scale.model = TRUE,
hyper = nPrior) +
f(myID,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(SexAgeGroup.f,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(Month,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(Int.DATE_YR, #year
model="rw1",
constr = TRUE,
scale.model = TRUE,
hyper = nPrior ) +
f(SumDur,#summer duration
model="rw1",
constr = TRUE,
scale.model = TRUE,
hyper = nPrior ) +
TOTAL
#theta.m1 = YA.Pregnancy.mod$internal.summary.hyperpar$mean
theta.m1 = c(-0.2969193, 0.6452064, -0.5640170, -0.6115045, -0.2670323, -0.0485563)
YA.Pregnancy.mod = inla(YA.Pregnancy.frm,
data = inla.stack.data(Preg2.stk),
family = "binomial",
verbose = TRUE,
control.fixed = list(prec = 0.001,
prec.intercept = 0.0001),
control.predictor = list(
A = inla.stack.A(Preg2.stk),
compute = TRUE,
link = 1),
control.mode = list(restart = TRUE, theta = theta.m1),
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(YA.Pregnancy.mod)##
## Call:
## c("inla(formula = YA.Pregnancy.frm, family = \"binomial\", data = inla.stack.data(Preg2.stk), ", " verbose = TRUE, control.compute = list(dic = TRUE, cpo = TRUE, ", " waic = TRUE), control.predictor = list(A = inla.stack.A(Preg2.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 = theta.m1))")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 1.9870 5.6568 0.2064 7.8502
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept1 -1.2476 0.4280 -2.0880 -1.2477 -0.4079 -1.2476 0
## TOTAL 0.0421 0.0046 0.0331 0.0421 0.0511 0.0421 0
##
## Random effects:
## Name Model
## COMRECPOP RW1 model
## myID IID model
## SexAgeGroup.f IID model
## Month IID model
## Int.DATE_YR RW1 model
## SumDur RW1 model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant 0.975quant
## Precision for COMRECPOP 0.8574 0.2788 0.4391 0.8148 1.5227
## Precision for myID 1.7157 0.4523 1.0045 1.6564 2.7640
## Precision for SexAgeGroup.f 1.0868 0.3494 0.5585 1.0349 1.9185
## Precision for Month 0.4674 0.1122 0.2839 0.4550 0.7233
## Precision for Int.DATE_YR 0.9660 0.3088 0.4986 0.9201 1.7004
## Precision for SumDur 0.7110 0.2020 0.3959 0.6838 1.1837
## mode
## Precision for COMRECPOP 0.7363
## Precision for myID 1.5442
## Precision for SexAgeGroup.f 0.9387
## Precision for Month 0.4315
## Precision for Int.DATE_YR 0.8350
## Precision for SumDur 0.6326
##
## Expected number of effective parameters(std dev): 65.16(0.00)
## Number of equivalent replicates : 18.85
##
## Deviance Information Criterion (DIC) ...............: 556.18
## Deviance Information Criterion (DIC, saturated) ....: 556.18
## Effective number of parameters .....................: 67.90
##
## Watanabe-Akaike information criterion (WAIC) ...: 544.83
## Effective number of parameters .................: 51.15
##
## Marginal log-Likelihood: -348.07
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computedModel for Pregnancy Yearling + Adult (Reconstructed Doe Population COMRECDOEPOP)
YA.Pregnancy.rd.frm = Preg ~ -1 + intercept1 +
f(COMRECDOEPOP,
model="rw1",
constr = TRUE,
scale.model = TRUE,
hyper = nPrior) +
f(myID,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(SexAgeGroup.f,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(Month,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(Int.DATE_YR, #year
model="rw1",
constr = TRUE,
scale.model = TRUE,
hyper = nPrior ) +
f(SumDur,#summer duration
model="rw1",
constr = TRUE,
scale.model = TRUE,
hyper = nPrior ) +
TOTAL
#theta.m1 = YA.Pregnancy.rd.mod$internal.summary.hyperpar$mean
theta.m1 = c(-0.2969193, 0.6452064, -0.5640170, -0.6115045, -0.2670323, -0.0485563)
YA.Pregnancy.rd.mod = inla(YA.Pregnancy.rd.frm,
data = inla.stack.data(Preg2.stk),
family = "binomial",
verbose = TRUE,
control.fixed = list(prec = 0.001,
prec.intercept = 0.0001),
control.predictor = list(
A = inla.stack.A(Preg2.stk),
compute = TRUE,
link = 1),
control.mode = list(restart = TRUE, theta = theta.m1),
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(YA.Pregnancy.rd.mod)##
## Call:
## c("inla(formula = YA.Pregnancy.rd.frm, family = \"binomial\", data = inla.stack.data(Preg2.stk), ", " verbose = TRUE, control.compute = list(dic = TRUE, cpo = TRUE, ", " waic = TRUE), control.predictor = list(A = inla.stack.A(Preg2.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 = theta.m1))")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 1.9420 5.1857 0.2385 7.3662
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept1 -1.2053 0.4221 -2.0341 -1.2053 -0.3772 -1.2053 0
## TOTAL 0.0415 0.0045 0.0326 0.0415 0.0504 0.0415 0
##
## Random effects:
## Name Model
## COMRECDOEPOP RW1 model
## myID IID model
## SexAgeGroup.f IID model
## Month IID model
## Int.DATE_YR RW1 model
## SumDur RW1 model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant 0.975quant
## Precision for COMRECDOEPOP 0.9113 0.2927 0.4699 0.8673 1.6082
## Precision for myID 1.7117 0.4526 1.0002 1.6522 2.7616
## Precision for SexAgeGroup.f 1.0875 0.3496 0.5592 1.0354 1.9197
## Precision for Month 0.4731 0.1137 0.2873 0.4605 0.7323
## Precision for Int.DATE_YR 0.9075 0.2847 0.4752 0.8655 1.5830
## Precision for SumDur 0.7510 0.2159 0.4154 0.7215 1.2569
## mode
## Precision for COMRECDOEPOP 0.7860
## Precision for myID 1.5398
## Precision for SexAgeGroup.f 0.9391
## Precision for Month 0.4366
## Precision for Int.DATE_YR 0.7877
## Precision for SumDur 0.6661
##
## Expected number of effective parameters(std dev): 65.29(0.00)
## Number of equivalent replicates : 18.81
##
## Deviance Information Criterion (DIC) ...............: 558.35
## Deviance Information Criterion (DIC, saturated) ....: 558.35
## Effective number of parameters .....................: 68.00
##
## Watanabe-Akaike information criterion (WAIC) ...: 547.15
## Effective number of parameters .................: 51.39
##
## Marginal log-Likelihood: -349.02
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computedModel for Pregnancy Yearling + Adult (Aerial Counts COMPRECNTDEN)
YA.Pregnancy.ac.frm = Preg ~ -1 + intercept1 +
f(COMPRECNTDEN,
model="rw1",
constr = TRUE,
scale.model = TRUE,
hyper = nPrior) +
f(myID,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(SexAgeGroup.f,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(Month,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(Int.DATE_YR, #year
model="rw1",
constr = TRUE,
scale.model = TRUE,
hyper = nPrior ) +
f(SumDur,#summer duration
model="rw1",
constr = TRUE,
scale.model = TRUE,
hyper = nPrior ) +
TOTAL
#theta.m1 = YA.Pregnancy.ac.mod$internal.summary.hyperpar$mean
theta.m1 = c(-0.06524442, 0.62028369, -0.52218859, -0.58439412, -0.35774291, -0.05702578)
YA.Pregnancy.ac.mod = inla(YA.Pregnancy.ac.frm,
data = inla.stack.data(Preg2.stk),
family = "binomial",
verbose = TRUE,
control.fixed = list(prec = 0.001,
prec.intercept = 0.0001),
control.predictor = list(
A = inla.stack.A(Preg2.stk),
compute = TRUE,
link = 1),
control.mode = list(restart = TRUE, theta = theta.m1),
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(YA.Pregnancy.ac.mod)##
## Call:
## c("inla(formula = YA.Pregnancy.ac.frm, family = \"binomial\", data = inla.stack.data(Preg2.stk), ", " verbose = TRUE, control.compute = list(dic = TRUE, cpo = TRUE, ", " waic = TRUE), control.predictor = list(A = inla.stack.A(Preg2.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 = theta.m1))")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 1.8897 5.4433 0.2424 7.5754
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept1 -1.2439 0.4119 -2.0526 -1.2439 -0.4358 -1.2439 0
## TOTAL 0.0419 0.0046 0.0329 0.0419 0.0508 0.0419 0
##
## Random effects:
## Name Model
## COMPRECNTDEN RW1 model
## myID IID model
## SexAgeGroup.f IID model
## Month IID model
## Int.DATE_YR RW1 model
## SumDur RW1 model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant 0.975quant
## Precision for COMPRECNTDEN 0.9747 0.3269 0.4857 0.9243 1.7570
## Precision for myID 1.7261 0.4556 1.0093 1.6663 2.7831
## Precision for SexAgeGroup.f 1.0873 0.3496 0.5587 1.0354 1.9193
## Precision for Month 0.4861 0.1176 0.2943 0.4730 0.7544
## Precision for Int.DATE_YR 0.8601 0.2822 0.4376 0.8167 1.5343
## Precision for SumDur 0.7499 0.2107 0.4198 0.7219 1.2414
## mode
## Precision for COMPRECNTDEN 0.8317
## Precision for myID 1.5534
## Precision for SexAgeGroup.f 0.9392
## Precision for Month 0.4482
## Precision for Int.DATE_YR 0.7369
## Precision for SumDur 0.6692
##
## Expected number of effective parameters(std dev): 65.69(0.00)
## Number of equivalent replicates : 18.69
##
## Deviance Information Criterion (DIC) ...............: 565.11
## Deviance Information Criterion (DIC, saturated) ....: 565.11
## Effective number of parameters .....................: 68.17
##
## Watanabe-Akaike information criterion (WAIC) ...: 554.07
## Effective number of parameters .................: 51.79
##
## Marginal log-Likelihood: -357.70
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computedCompare DIC and WAIC
#DIC
YA.Pregnancy.mod$dic$dic## [1] 556.1793YA.Pregnancy.rd.mod$dic$dic## [1] 558.3519YA.Pregnancy.ac.mod$dic$dic## [1] 565.1072#WAIC
YA.Pregnancy.mod$waic$waic## [1] 544.8303YA.Pregnancy.rd.mod$waic$waic## [1] 547.1549YA.Pregnancy.ac.mod$waic$waic## [1] 554.0672Yearling + Adult Pregnancy Density Dependence (Reconstructed Population, COMRECPOP)
mic.d = as.data.frame(YA.Pregnancy.mod$summary.random$COMRECPOP[,c(1:6)])
names(mic.d) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975")
mic.d$Meane = plogis(mic.d$Mean + YA.Pregnancy.mod$summary.fixed[1,1])
mic.d$Q025e = plogis(mic.d$Q025 + YA.Pregnancy.mod$summary.fixed[1,1])
mic.d$Q975e = plogis(mic.d$Q975 + YA.Pregnancy.mod$summary.fixed[1,1])
mic.d$ID = mic.d$ID/30.63 #Area of Complex 1
range(mic.d$ID)## [1] 8.651649 38.132550ggplot(mic.d, aes(ID, Meane)) +
geom_point(shape=1) +
geom_smooth(data= mic.d,
aes(ID, Meane),
col = "black",
se=FALSE,
span = 0.9,
method="loess",
lwd = 1) +
geom_smooth(data = mic.d,
aes(ID, Q025e),
col = "grey",
method = "loess",
span = 0.9,
se = FALSE,
linetype= "dashed") +
geom_smooth(data = mic.d, aes(ID, Q975e),
col = "grey",
method = "loess",
span = 0.9,
se = FALSE,
linetype= "dashed") +
geom_hline(yintercept = plogis(YA.Pregnancy.mod$summary.fixed[1,1]),
linetype = "solid",
col = "darkgray",
size = 0.5) +
ylab("Pregnancy (probability)") +
xlab(expression(Deer~km^2)) +
ggtitle("Pregnancy (Yearling + Adult)") +
theme_classic() +
#scale_y_continuous(breaks = seq(0, 0.30, 0.05), limits = c(-0.1, 0.30)) +
scale_x_continuous(breaks = seq(8, 40, 5), limits =c(8,40)) +
theme(axis.text=element_text(size=16),
legend.title = element_text(size=16, face="bold"),
legend.position = "none",
legend.key = element_rect(fill = "white", linetype=0),
legend.background = element_rect(fill = "white", linetype=0),
plot.title = element_text(face="bold", size = 20, hjust=0.5),
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))
Sex-Age Class Effect
Horizontal line is doe population average.
mic.d = as.data.frame(YA.Pregnancy.mod$summary.random$SexAgeGroup.f)
names(mic.d) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975", "mode", "kld")
mic.d$Mean = plogis(mic.d$Mean + YA.Pregnancy.mod$summary.fixed[1,1])
mic.d$Q025 = plogis(mic.d$Q025 + YA.Pregnancy.mod$summary.fixed[1,1])
mic.d$Q975 = plogis(mic.d$Q975 + YA.Pregnancy.mod$summary.fixed[1,1])
range(mic.d$Mean)## [1] 0.1702061 0.2867718ggplot(mic.d, aes(factor(ID), y=Mean)) +
geom_point(size=2, pch=19, col = "red") +
geom_linerange(aes(ymin=Q025, ymax=Q975), colour="black") +
geom_hline(yintercept = plogis(YA.Pregnancy.mod$summary.fixed[1,1]),
linetype = "solid",
col = "darkgray",
size = 0.5) +
theme_classic() +
xlab(" ") +
ylab("Pregnancy (probability)") +
ggtitle("Age Effect (Yearling + Adult)") +
scale_x_discrete(labels = c("Yearling (f)", "Adult (f)")) +
scale_y_continuous(breaks = seq(0, 0.5, 0.05), limits = c(-0.05, 0.5)) +
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. Horizontal line is doe population mean.
mic.d = as.data.frame(YA.Pregnancy.mod$summary.random$Month)
names(mic.d) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975", "mode", "kld")
mic.d$Mean = plogis(mic.d$Mean + YA.Pregnancy.mod$summary.fixed[1,1])
mic.d$Q025 = plogis(mic.d$Q025 + YA.Pregnancy.mod$summary.fixed[1,1])
mic.d$Q975 = plogis(mic.d$Q975 + YA.Pregnancy.mod$summary.fixed[1,1])
range(mic.d$Mean)## [1] 0.005494404 0.802064554ggplot(mic.d, aes(factor(ID), y=Mean)) +
geom_point(size=2, pch=19, col = "red") +
geom_linerange(aes(ymin=Q025, ymax=Q975), colour="black") +
geom_hline(yintercept = plogis(YA.Pregnancy.mod$summary.fixed[1,1]),
linetype = "solid",
col = "darkgray",
size = 0.5) +
theme_classic() +
xlab(" ") +
ylab("Pregnancy (probability)") +
ggtitle("Month Effect (Yearling + Adult)") +
scale_x_discrete(labels = c("Jan", "Feb", "March",
"Oct", "Nov", "Dec")) +
#scale_y_continuous(breaks = seq(0, 0.6, 0.05), limits = c(-0.05, 0.6)) +
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))
2.5 Sex Ratio
2.5.1 All Age Groups
Organize Data.
Doe.df$M0 = Doe.df$M
length(Doe.df$M0)## [1] 2239length(which(is.na(Doe.df$M0)))## [1] 1483length(which(Doe.df$M0 >= 1))## [1] 547Doe.df$F0 = Doe.df$F
length(Doe.df$F0)## [1] 2239length(which(is.na(Doe.df$F0)))## [1] 1485length(which(Doe.df$F0 >= 1))## [1] 485Doe.df$M0z = Doe.df$M0
Doe.df$M0z[is.na(Doe.df$M0z)] = 0
Doe.df$F0z = Doe.df$F0
Doe.df$F0z[is.na(Doe.df$F0z)] = 0
Doe.df$FETUS3 = Doe.df$M0z + Doe.df$F0z
length(Doe.df$FETUS3)## [1] 2239length(which(is.na(Doe.df$FETUS3)))## [1] 0length(which(Doe.df$FETUS3 >= 1))## [1] 756Doe2.df = Doe.df %>% filter(is.na(M0) == FALSE)
dim(Doe.df)## [1] 2239 113dim(Doe2.df)## [1] 756 113#Unique ID by Deer
Doe2.df$myID = 1:nrow(Doe2.df)
#set up model data as a list object
## Not all covariates are used.
FE.lst3 = list(list(intercept1 = rep(1, dim(Doe2.df)[1])),
list(COMPREHARKM = Doe2.df[,"COMPREHARKM"],
COMRECPOP = Doe2.df[,"COMRECPOP"],
COMRECDOEPOP = Doe2.df[,"COMRECDOEPOP"],
COMCUMPRE2HARKM = Doe2.df[,"COMCUMPRE2HARKM"],
TOTAL = Doe2.df[,"TOTAL"],
LiveWtKg = Doe2.df[,"LiveWtKg"],
Lat = Doe2.df[,"Lat"],
NN = Doe2.df[,"NN"],
Area = Doe2.df[,"CmpxArea"],
HINDFOOT = Doe2.df[,"HINDFOOT"]/100,
COMPREDOEHARKM = Doe2.df[,"COMPREDOEHARKM"],
COMCUMPRE2DOEHARKM = Doe2.df[,"COMCUMPRE2DOEHARKM"],
COMPRECNTDEN = Doe2.df[,"COMPRECNTDEN"],
GrowSeas = round(Doe2.df[,"GrowSeas"],3),
GrowSeas.sc = Doe2.df[,"GrowSeas.sc"],
PrecEx10mm = Doe2.df[,"PrecEx10mm"],
PrecEx20mm = Doe2.df[,"PrecEx20mm"],
DSubZero = Doe2.df[,"DSubZero"],
Prec95pct = Doe2.df[,"Prec95pct"],
MonMxDT = Doe2.df[,"MonMxDT"],
MMx1DPrec = Doe2.df[,"MMx1DPrec"],
AMinT = Doe2.df[,"AMinT"],
AMxT = Doe2.df[,"AMxT"],
SumDur = round(Doe2.df[,"SumDur"],3),
SumDur.sc = Doe2.df[,"SumDur.sc"],
FrstDays = Doe2.df[,"FrstDays"],
Int.DATE_YR = Doe2.df[,"Int.DATE_YR"],
RepGroups = Doe2.df[,"RepGroups"],
SexAgeGroup.f = Doe2.df[,"SexAgeGroup.f"],
AGEGROUP = Doe2.df[,"AGEGROUP"],
PREG = as.factor(Doe2.df[,"PREG"]),
f_COMPLEX = Doe2.df[,"f_COMPLEX"],
myID = Doe2.df[,"myID"],
Month = Doe2.df[,"DATE_MO"]))
Doe2B.df = Doe2.df %>% filter(SexAgeGroup.f != "21")
dim(Doe2.df)## [1] 756 113dim(Doe2B.df)## [1] 733 113Doe2B.df$SexAgeGroup.f = factor(as.character(Doe2B.df$SexAgeGroup.f))
#Unique ID by Deer
Doe2B.df$myID = 1:nrow(Doe2B.df)
#Yearling + Adult Set
FE.lst3B = list(list(intercept1 = rep(1, dim(Doe2B.df)[1])),
list(COMPREHARKM = Doe2B.df[,"COMPREHARKM"],
COMRECPOP = Doe2B.df[,"COMRECPOP"],
COMRECDOEPOP = Doe2B.df[,"COMRECDOEPOP"],
COMCUMPRE2HARKM = Doe2B.df[,"COMCUMPRE2HARKM"],
TOTAL = Doe2B.df[,"TOTAL"],
LiveWtKg = Doe2B.df[,"LiveWtKg"],
Lat = Doe2B.df[,"Lat"],
NN = Doe2B.df[,"NN"],
Area = Doe2B.df[,"CmpxArea"],
HINDFOOT = Doe2B.df[,"HINDFOOT"]/100,
COMPREDOEHARKM = Doe2B.df[,"COMPREDOEHARKM"],
COMCUMPRE2DOEHARKM = Doe2B.df[,"COMCUMPRE2DOEHARKM"],
COMPRECNTDEN = Doe2B.df[,"COMPRECNTDEN"],
GrowSeas = round(Doe2B.df[,"GrowSeas"],3),
GrowSeas.sc = Doe2B.df[,"GrowSeas.sc"],
PrecEx10mm = Doe2B.df[,"PrecEx10mm"],
PrecEx20mm = Doe2B.df[,"PrecEx20mm"],
DSubZero = Doe2B.df[,"DSubZero"],
Prec95pct = Doe2B.df[,"Prec95pct"],
MonMxDT = Doe2B.df[,"MonMxDT"],
MMx1DPrec = Doe2B.df[,"MMx1DPrec"],
AMinT = Doe2B.df[,"AMinT"],
AMxT = Doe2B.df[,"AMxT"],
SumDur = round(Doe2B.df[,"SumDur"],3),
SumDur.sc = Doe2B.df[,"SumDur.sc"],
FrstDays = Doe2B.df[,"FrstDays"],
Int.DATE_YR = Doe2B.df[,"Int.DATE_YR"],
RepGroups = Doe2B.df[,"RepGroups"],
SexAgeGroup.f = Doe2B.df[,"SexAgeGroup.f"],
AGEGROUP = Doe2B.df[,"AGEGROUP"],
PREG = as.factor(Doe2B.df[,"PREG"]),
f_COMPLEX = Doe2B.df[,"f_COMPLEX"],
myID = Doe2B.df[,"myID"],
Month = Doe2B.df[,"DATE_MO"]))
Ratio.stk = inla.stack(data = list(Ratio = Doe2.df$M0, Trials = Doe2.df$FETUS3), #successes and trials (male)
A = list(1,1),
effects = FE.lst3,
tag = "sr.0")
Ratio2.stk = inla.stack(data = list(Ratio = Doe2B.df$M0, Trials = Doe2B.df$FETUS3), #Yearling + Adult Set
A = list(1,1),
effects = FE.lst3B,
tag = "ya.0")Model for Sex Ratio (proportion male) (Reconstructed Population COMRECPOP)
SexRatio.frm = Ratio ~ -1 + intercept1 +
f(COMRECPOP,
model="rw1",
constr = TRUE,
scale.model = TRUE,
replicate = RepGroups,
hyper = nPrior) +
f(myID,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(Int.DATE_YR,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(SexAgeGroup.f,
model="iid",
constr = TRUE,
hyper = nPrior)
#theta.m1 = SexRatio.mod$internal.summary.hyperpar$mean
theta.m1 = c(0.36618147, 0.78887100, 0.48440643, 0.08022418)
SexRatio.mod = inla(SexRatio.frm,
data = inla.stack.data(Ratio.stk),
family = "binomial",
verbose = TRUE,
Ntrials = inla.stack.data(Ratio.stk)$Trials, #Trials
control.fixed = list(prec = 0.001,
prec.intercept = 0.0001),
control.predictor = list(
A = inla.stack.A(Ratio.stk),
compute = TRUE,
link = 1),
control.mode = list(restart = TRUE, theta = theta.m1),
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(SexRatio.mod)##
## Call:
## c("inla(formula = SexRatio.frm, family = \"binomial\", data = inla.stack.data(Ratio.stk), ", " Ntrials = inla.stack.data(Ratio.stk)$Trials, verbose = TRUE, ", " control.compute = list(dic = TRUE, cpo = TRUE, waic = TRUE), ", " control.predictor = list(A = inla.stack.A(Ratio.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 = theta.m1))")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 1.5684 1.4780 0.1416 3.1881
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept1 0.2714 0.1608 -0.0442 0.2714 0.5868 0.2714 0
##
## Random effects:
## Name Model
## COMRECPOP RW1 model
## myID IID model
## Int.DATE_YR IID model
## SexAgeGroup.f IID model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant 0.975quant
## Precision for COMRECPOP 1.507 0.4607 0.8015 1.442 2.596
## Precision for myID 2.253 0.5019 1.4425 2.192 3.405
## Precision for Int.DATE_YR 1.689 0.5037 0.9096 1.619 2.874
## Precision for SexAgeGroup.f 1.139 0.3660 0.5850 1.085 2.010
## mode
## Precision for COMRECPOP 1.3192
## Precision for myID 2.0729
## Precision for Int.DATE_YR 1.4889
## Precision for SexAgeGroup.f 0.9847
##
## Expected number of effective parameters(std dev): 138.67(0.00)
## Number of equivalent replicates : 5.452
##
## Deviance Information Criterion (DIC) ...............: 1430.10
## Deviance Information Criterion (DIC, saturated) ....: 1044.21
## Effective number of parameters .....................: 132.58
##
## Watanabe-Akaike information criterion (WAIC) ...: 1421.11
## Effective number of parameters .................: 107.89
##
## Marginal log-Likelihood: -770.92
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computedModel for Sex Ratio (proportion male) (Reconstructed Doe Population COMRECDOEPOP)
SexRatio.rd.frm = Ratio ~ -1 + intercept1 +
f(COMRECDOEPOP,
model="rw1",
constr = TRUE,
scale.model = TRUE,
replicate = RepGroups,
hyper = nPrior) +
f(myID,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(Int.DATE_YR,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(SexAgeGroup.f,
model="iid",
constr = TRUE,
hyper = nPrior)
#theta.m1 = SexRatio.rd.mod$internal.summary.hyperpar$mean
theta.m1 = c(0.3715700, 0.7865159, 0.4925481, 0.0806010)
SexRatio.rd.mod = inla(SexRatio.rd.frm,
data = inla.stack.data(Ratio.stk),
family = "binomial",
verbose = TRUE,
Ntrials = inla.stack.data(Ratio.stk)$Trials, #Trials
control.fixed = list(prec = 0.001,
prec.intercept = 0.0001),
control.predictor = list(
A = inla.stack.A(Ratio.stk),
compute = TRUE,
link = 1),
control.mode = list(restart = TRUE, theta = theta.m1),
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(SexRatio.rd.mod)##
## Call:
## c("inla(formula = SexRatio.rd.frm, family = \"binomial\", data = inla.stack.data(Ratio.stk), ", " Ntrials = inla.stack.data(Ratio.stk)$Trials, verbose = TRUE, ", " control.compute = list(dic = TRUE, cpo = TRUE, waic = TRUE), ", " control.predictor = list(A = inla.stack.A(Ratio.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 = theta.m1))")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 1.5369 1.4631 0.1287 3.1287
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept1 0.2684 0.1617 -0.0491 0.2684 0.5856 0.2684 0
##
## Random effects:
## Name Model
## COMRECDOEPOP RW1 model
## myID IID model
## Int.DATE_YR IID model
## SexAgeGroup.f IID model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant 0.975quant
## Precision for COMRECDOEPOP 1.516 0.4629 0.8066 1.450 2.610
## Precision for myID 2.248 0.5008 1.4392 2.187 3.397
## Precision for Int.DATE_YR 1.703 0.5082 0.9171 1.633 2.899
## Precision for SexAgeGroup.f 1.140 0.3661 0.5851 1.085 2.010
## mode
## Precision for COMRECDOEPOP 1.3265
## Precision for myID 2.0685
## Precision for Int.DATE_YR 1.5015
## Precision for SexAgeGroup.f 0.9849
##
## Expected number of effective parameters(std dev): 138.86(0.00)
## Number of equivalent replicates : 5.444
##
## Deviance Information Criterion (DIC) ...............: 1431.12
## Deviance Information Criterion (DIC, saturated) ....: 1045.23
## Effective number of parameters .....................: 132.75
##
## Watanabe-Akaike information criterion (WAIC) ...: 1422.23
## Effective number of parameters .................: 108.12
##
## Marginal log-Likelihood: -773.82
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computedModel for Sex Ratio (proportion male) (Aerial Counts COMPRECNTDEN)
SexRatio.ac.frm = Ratio ~ -1 + intercept1 +
f(COMPRECNTDEN,
model="rw1",
constr = TRUE,
scale.model = TRUE,
replicate = RepGroups,
hyper = nPrior) +
f(myID,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(Int.DATE_YR,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(SexAgeGroup.f,
model="iid",
constr = TRUE,
hyper = nPrior)
#theta.m1 = SexRatio.ac.mod$internal.summary.hyperpar$mean
theta.m1 = c(0.31736554, 0.79131261, 0.48778867, 0.07937508)
SexRatio.ac.mod = inla(SexRatio.ac.frm,
data = inla.stack.data(Ratio.stk),
family = "binomial",
verbose = TRUE,
Ntrials = inla.stack.data(Ratio.stk)$Trials, #Trials
control.fixed = list(prec = 0.001,
prec.intercept = 0.0001),
control.predictor = list(
A = inla.stack.A(Ratio.stk),
compute = TRUE,
link = 1),
control.mode = list(restart = TRUE, theta = theta.m1),
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(SexRatio.ac.mod)##
## Call:
## c("inla(formula = SexRatio.ac.frm, family = \"binomial\", data = inla.stack.data(Ratio.stk), ", " Ntrials = inla.stack.data(Ratio.stk)$Trials, verbose = TRUE, ", " control.compute = list(dic = TRUE, cpo = TRUE, waic = TRUE), ", " control.predictor = list(A = inla.stack.A(Ratio.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 = theta.m1))")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 1.4813 1.4612 0.1287 3.0711
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept1 0.3001 0.1625 -0.019 0.3001 0.619 0.3001 0
##
## Random effects:
## Name Model
## COMPRECNTDEN RW1 model
## myID IID model
## Int.DATE_YR IID model
## SexAgeGroup.f IID model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant 0.975quant
## Precision for COMPRECNTDEN 1.436 0.4390 0.7636 1.373 2.474
## Precision for myID 2.258 0.5032 1.4455 2.197 3.413
## Precision for Int.DATE_YR 1.698 0.5081 0.9137 1.628 2.895
## Precision for SexAgeGroup.f 1.138 0.3654 0.5842 1.084 2.007
## mode
## Precision for COMPRECNTDEN 1.2568
## Precision for myID 2.0771
## Precision for Int.DATE_YR 1.4952
## Precision for SexAgeGroup.f 0.9839
##
## Expected number of effective parameters(std dev): 138.73(0.00)
## Number of equivalent replicates : 5.45
##
## Deviance Information Criterion (DIC) ...............: 1429.09
## Deviance Information Criterion (DIC, saturated) ....: 1043.20
## Effective number of parameters .....................: 132.64
##
## Watanabe-Akaike information criterion (WAIC) ...: 1420.06
## Effective number of parameters .................: 107.87
##
## Marginal log-Likelihood: -783.86
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computedCompare DIC and WAIC
#DIC
SexRatio.mod$dic$dic## [1] 1430.098SexRatio.rd.mod$dic$dic## [1] 1431.119SexRatio.ac.mod$dic$dic## [1] 1429.087#WAIC
SexRatio.mod$waic$waic## [1] 1421.109SexRatio.rd.mod$waic$waic## [1] 1422.233SexRatio.ac.mod$waic$waic## [1] 1420.062Sex Ratio (Reconstructed Population, COMRECPOP)
Zero is excluded from credible interval.
mic.d = as.data.frame(SexRatio.mod$summary.random$COMRECPOP[,c(1:6)])
names(mic.d) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975")
length(unique(mic.d$ID))## [1] 13length(mic.d$ID)## [1] 39mic.d$ID2 = as.factor(
rep(c("Fawn (f)", "Yearling (f)", "Adult (f)"),
each = length(unique(mic.d$ID))))
mic.d$Meane = plogis(mic.d$Mean + SexRatio.mod$summary.fixed[1,1])
mic.d$Q025e = plogis(mic.d$Q025 + SexRatio.mod$summary.fixed[1,1])
mic.d$Q975e = plogis(mic.d$Q975 + SexRatio.mod$summary.fixed[1,1])
mic.d$ID = mic.d$ID/30.63 #Area of Complex 1
range(mic.d$ID)## [1] 8.651649 38.132550ggplot(mic.d, aes(ID, Meane, group = ID2)) +
geom_point(shape=1) +
geom_smooth(data= mic.d,
aes(ID, Meane,
group = mic.d$ID2,
col = factor(mic.d$ID2),
linetype= factor(mic.d$ID2)),
se=FALSE,
span = 0.9,
method="loess",
lwd = 1) +
geom_hline(yintercept = plogis(SexRatio.mod$summary.fixed[1,1]),
linetype = "solid",
col = "darkgray",
size = 0.5) +
scale_linetype_manual(name = "Age Group",
values=c("solid", "longdash", "dotted")) +
scale_colour_manual(name = "Age Group",
values=c("red", "brown", "blue")) +
ylab("Male (probability)") +
xlab(expression(Deer~km^2)) +
ggtitle("Proportion Male (COMRECPOP)") +
theme_classic() +
scale_x_continuous(breaks = seq(8, 40, 5), limits =c(8,40)) +
theme(axis.text=element_text(size=16),
legend.title = element_text(size=16, face="bold"),
legend.position = c(0.2,0.2),
legend.key = element_rect(fill = "white", linetype=0),
legend.background = element_rect(fill = "white", linetype=0),
plot.title = element_text(face="bold", size = 20, hjust=0.5),
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))
Sex Ratio (Aerial Counts, COMPRECNTDEN)
mic.d = as.data.frame(SexRatio.ac.mod$summary.random$COMPRECNTDEN[,c(1:6)])
names(mic.d) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975")
length(unique(mic.d$ID))## [1] 16length(mic.d$ID)## [1] 48mic.d$ID2 = as.factor(
rep(c("Fawn (f)", "Yearling (f)", "Adult (f)"),
each = length(unique(mic.d$ID))))
mic.d$Meane = plogis(mic.d$Mean + SexRatio.mod$summary.fixed[1,1])
mic.d$Q025e = plogis(mic.d$Q025 + SexRatio.mod$summary.fixed[1,1])
mic.d$Q975e = plogis(mic.d$Q975 + SexRatio.mod$summary.fixed[1,1])
range(mic.d$ID)## [1] 4.3 31.0ggplot(mic.d, aes(ID, Meane, group = ID2)) +
geom_point(shape=1) +
geom_smooth(data= mic.d,
aes(ID, Meane,
group = mic.d$ID2,
col = factor(mic.d$ID2),
linetype= factor(mic.d$ID2)),
se=FALSE,
span = 0.9,
method="loess",
lwd = 1) +
geom_hline(yintercept = plogis(SexRatio.mod$summary.fixed[1,1]),
linetype = "solid",
col = "darkgray",
size = 0.5) +
scale_linetype_manual(name = "Age Group",
values=c("solid", "longdash", "dotted")) +
scale_colour_manual(name = "Age Group",
values=c("red", "brown", "blue")) +
ylab("Male (probability)") +
xlab(expression(Deer~km^2)) +
ggtitle("Proportion Male (COMPRECNTDEN)") +
theme_classic() +
scale_x_continuous(breaks = seq(4, 35, 5), limits =c(5,35)) +
theme(axis.text=element_text(size=16),
legend.title = element_text(size=16, face="bold"),
legend.position = c(0.8,0.8),
legend.key = element_rect(fill = "white", linetype=0),
legend.background = element_rect(fill = "white", linetype=0),
plot.title = element_text(face="bold", size = 20, hjust=0.5),
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))
Model Fitted Estimates.
Temp.data = Doe2.df
idat = inla.stack.index(Ratio.stk, "sr.0")$data #Indexregion data from model
Temp.data$Pred.Preg = SexRatio.mod$summary.fitted.values$mean[idat] #Fitted values for region
Temp.data$Low = SexRatio.mod$summary.fitted.values$`0.025quant`[idat]
Temp.data$High = SexRatio.mod$summary.fitted.values$`0.975quant`[idat]
Preg.fit = as.data.frame(
Temp.data %>%
group_by(SexAgeGroup.f) %>%
summarise(Mean = mean(Pred.Preg, na.rm=T),
Low = mean(Low, na.rm=T),
High = mean(High, na.rm=T)))
Preg.fit$Age = c("Fawn (f)", "Yearling (f)", "Adult (f)")
Preg.fit$Age = ordered(factor(Preg.fit$Age), levels = c("Fawn (f)", "Yearling (f)", "Adult (f)"))
ggplot(Preg.fit, aes(factor(Age),Mean)) +
geom_bar(stat = "identity", fill = "lightgray") +
geom_point(size=2, pch=19, col = "red") +
geom_linerange(aes(ymin=Low, ymax=High), colour="black") +
theme_classic() +
ylab("Mean Probabilty of Male Fetus") +
xlab( "") +
ggtitle("Population Mean Probabilty of Male Fetus \n (COMRECPOP)") +
scale_x_discrete(labels = c("Fawn (f)", "Yearling (f)", "Adult (f)")) +
scale_y_continuous(breaks = seq(0, 1, 0.25), limits = c(0, 1)) +
theme(axis.text=element_text(size=16),
legend.title = element_text(size=16, face="bold"),
legend.position = c(0.7,0.7),
legend.key = element_rect(fill = "white", linetype=0),
legend.background = element_rect(fill = "white", linetype=0),
plot.title = element_text(face="bold", size = 20, hjust=0.5),
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))
Temp.data$Class = ifelse(Temp.data$SexAgeGeoup == "21", "Fawn (f)",
ifelse(Temp.data$SexAgeGeoup == "22", "Yearling (f)",
ifelse(Temp.data$SexAgeGeoup == "23", "Adult (f)", NA)))
Temp.data$Class = ordered(factor(Temp.data$Class, levels = c("Fawn (f)", "Yearling (f)", "Adult (f)")))
Temp.data$Density = Temp.data$COMRECPOP/30.63
COMRECPOP.plt = Temp.data %>% select(Class, Pred.Preg)
hline_dat1 = as.data.frame(
COMRECPOP.plt %>%
group_by(Class) %>%
summarise(Med = mean(Pred.Preg)))
ggplot(Temp.data, aes(Density, Pred.Preg)) +
geom_point(shape=1) +
geom_smooth(method="loess",
se=FALSE,
col="black",
lwd = 1) +
facet_wrap(~Class) +
geom_hline(data=hline_dat1,
aes(yintercept=Med),
col = "darkgray",
size = 0.5) +
ylab("Probability of Male Fetus") +
xlab(expression(Deer~km^2)) +
ggtitle("Estimated Probabilty of Male Fetus \n (Complex 1, COMRECPOP)") +
theme_classic() +
theme(axis.text=element_text(size=16),
plot.title = element_text(face="bold", size = 20, hjust=0.5),
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))
Sex-Age Class Effect
Horizontal line is doe population average.
Zero is excluded from credible interval.
mic.d = as.data.frame(SexRatio.mod$summary.random$SexAgeGroup.f)
names(mic.d) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975", "mode", "kld")
mic.d$Mean = plogis(mic.d$Mean + SexRatio.mod$summary.fixed[1,1])
mic.d$Q025 = plogis(mic.d$Q025 + SexRatio.mod$summary.fixed[1,1])
mic.d$Q975 = plogis(mic.d$Q975 + SexRatio.mod$summary.fixed[1,1])
range(mic.d$Mean)## [1] 0.4889161 0.6134600ggplot(mic.d, aes(factor(ID), y=Mean)) +
geom_point(size=2, pch=19, col = "red") +
geom_linerange(aes(ymin=Q025, ymax=Q975), colour="black") +
geom_hline(yintercept = plogis(SexRatio.mod$summary.fixed[1,1]),
linetype = "solid",
col = "darkgray",
size = 0.5) +
theme_classic() +
xlab(" ") +
ylab("Male (probability)") +
ggtitle("Age Effect") +
scale_x_discrete(labels = c("Fawn (f)", "Yearling (f)", "Adult (f)")) +
#scale_y_continuous(breaks = seq(0, 0.7, 0.05), limits = c(-0.05, 0.7)) +
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))
2.5.2 Yearling and Adult
Model for Yearling and Adult Sex Ratio (proportion male) (Reconstructed Population COMRECPOP)
YA.SexRatio.frm = Ratio ~ -1 + intercept1 +
f(COMRECPOP,
model="rw1",
constr = TRUE,
scale.model = TRUE,
hyper = nPrior) +
f(myID,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(Int.DATE_YR,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(SexAgeGroup.f,
model="iid",
constr = TRUE,
hyper = nPrior)
#theta.m1 = YA.SexRatio.mod$internal.summary.hyperpar$mean
theta.m1 = c(0.36618147, 0.78887100, 0.48440643, 0.08022418)
YA.SexRatio.mod = inla(YA.SexRatio.frm,
data = inla.stack.data(Ratio2.stk),
family = "binomial",
verbose = TRUE,
Ntrials = inla.stack.data(Ratio2.stk)$Trials, #Trials
control.fixed = list(prec = 0.001,
prec.intercept = 0.0001),
control.predictor = list(
A = inla.stack.A(Ratio2.stk),
compute = TRUE,
link = 1),
control.mode = list(restart = TRUE, theta = theta.m1),
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(YA.SexRatio.mod)##
## Call:
## c("inla(formula = YA.SexRatio.frm, family = \"binomial\", data = inla.stack.data(Ratio2.stk), ", " Ntrials = inla.stack.data(Ratio2.stk)$Trials, verbose = TRUE, ", " control.compute = list(dic = TRUE, cpo = TRUE, waic = TRUE), ", " control.predictor = list(A = inla.stack.A(Ratio2.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 = theta.m1))")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 1.6735 1.2631 0.1355 3.0720
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept1 0.1781 0.081 0.019 0.1781 0.3371 0.1781 0
##
## Random effects:
## Name Model
## COMRECPOP RW1 model
## myID IID model
## Int.DATE_YR IID model
## SexAgeGroup.f IID model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant 0.975quant
## Precision for COMRECPOP 1.221 0.3845 0.6373 1.165 2.134
## Precision for myID 2.275 0.5056 1.4553 2.214 3.430
## Precision for Int.DATE_YR 1.643 0.4921 0.8797 1.575 2.799
## Precision for SexAgeGroup.f 1.099 0.3551 0.5631 1.046 1.945
## mode
## Precision for COMRECPOP 1.060
## Precision for myID 2.096
## Precision for Int.DATE_YR 1.449
## Precision for SexAgeGroup.f 0.948
##
## Expected number of effective parameters(std dev): 129.26(0.00)
## Number of equivalent replicates : 5.671
##
## Deviance Information Criterion (DIC) ...............: 1391.53
## Deviance Information Criterion (DIC, saturated) ....: 1008.41
## Effective number of parameters .....................: 123.86
##
## Watanabe-Akaike information criterion (WAIC) ...: 1383.04
## Effective number of parameters .................: 101.03
##
## Marginal log-Likelihood: -727.31
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computedModel for Yearling and Adult Sex Ratio (proportion male) (Reconstructed Doe Population COMRECDOEPOP)
YA.SexRatio.rd.frm = Ratio ~ -1 + intercept1 +
f(COMRECDOEPOP,
model="rw1",
constr = TRUE,
scale.model = TRUE,
hyper = nPrior) +
f(myID,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(Int.DATE_YR,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(SexAgeGroup.f,
model="iid",
constr = TRUE,
hyper = nPrior)
#theta.m1 = YA.SexRatio.rd.mod$internal.summary.hyperpar$mean
theta.m1 = c(0.3715700, 0.7865159, 0.4925481, 0.0806010)
YA.SexRatio.rd.mod = inla(YA.SexRatio.rd.frm,
data = inla.stack.data(Ratio2.stk),
family = "binomial",
verbose = TRUE,
Ntrials = inla.stack.data(Ratio2.stk)$Trials, #Trials
control.fixed = list(prec = 0.001,
prec.intercept = 0.0001),
control.predictor = list(
A = inla.stack.A(Ratio2.stk),
compute = TRUE,
link = 1),
control.mode = list(restart = TRUE, theta = theta.m1),
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(YA.SexRatio.rd.mod)##
## Call:
## c("inla(formula = YA.SexRatio.rd.frm, family = \"binomial\", data = inla.stack.data(Ratio2.stk), ", " Ntrials = inla.stack.data(Ratio2.stk)$Trials, verbose = TRUE, ", " control.compute = list(dic = TRUE, cpo = TRUE, waic = TRUE), ", " control.predictor = list(A = inla.stack.A(Ratio2.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 = theta.m1))")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 1.7053 1.4323 0.1446 3.2822
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept1 0.1785 0.0809 0.0197 0.1785 0.3373 0.1785 0
##
## Random effects:
## Name Model
## COMRECDOEPOP RW1 model
## myID IID model
## Int.DATE_YR IID model
## SexAgeGroup.f IID model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant 0.975quant
## Precision for COMRECDOEPOP 1.212 0.3819 0.6325 1.157 2.119
## Precision for myID 2.275 0.5057 1.4557 2.215 3.431
## Precision for Int.DATE_YR 1.649 0.4934 0.8837 1.581 2.808
## Precision for SexAgeGroup.f 1.099 0.3552 0.5635 1.046 1.945
## mode
## Precision for COMRECDOEPOP 1.0528
## Precision for myID 2.0973
## Precision for Int.DATE_YR 1.4545
## Precision for SexAgeGroup.f 0.9469
##
## Expected number of effective parameters(std dev): 129.17(0.00)
## Number of equivalent replicates : 5.675
##
## Deviance Information Criterion (DIC) ...............: 1391.42
## Deviance Information Criterion (DIC, saturated) ....: 1008.30
## Effective number of parameters .....................: 123.78
##
## Watanabe-Akaike information criterion (WAIC) ...: 1382.94
## Effective number of parameters .................: 100.97
##
## Marginal log-Likelihood: -727.99
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computedModel for Yearling and Adult Sex Ratio (proportion male) (Aerial Counts COMPRECNTDEN)
YA.SexRatio.ac.frm = Ratio ~ -1 + intercept1 +
f(COMPRECNTDEN,
model="rw1",
constr = TRUE,
scale.model = TRUE,
hyper = nPrior) +
f(myID,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(Int.DATE_YR,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(SexAgeGroup.f,
model="iid",
constr = TRUE,
hyper = nPrior)
#theta.m1 = YA.SexRatio.ac.mod$internal.summary.hyperpar$mean
theta.m1 = c(0.31736554, 0.79131261, 0.48778867, 0.07937508)
YA.SexRatio.ac.mod = inla(YA.SexRatio.ac.frm,
data = inla.stack.data(Ratio2.stk),
family = "binomial",
verbose = TRUE,
Ntrials = inla.stack.data(Ratio2.stk)$Trials, #Trials
control.fixed = list(prec = 0.001,
prec.intercept = 0.0001),
control.predictor = list(
A = inla.stack.A(Ratio2.stk),
compute = TRUE,
link = 1),
control.mode = list(restart = TRUE, theta = theta.m1),
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(YA.SexRatio.ac.mod)##
## Call:
## c("inla(formula = YA.SexRatio.ac.frm, family = \"binomial\", data = inla.stack.data(Ratio2.stk), ", " Ntrials = inla.stack.data(Ratio2.stk)$Trials, verbose = TRUE, ", " control.compute = list(dic = TRUE, cpo = TRUE, waic = TRUE), ", " control.predictor = list(A = inla.stack.A(Ratio2.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 = theta.m1))")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 1.4980 1.3738 0.1247 2.9965
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept1 0.1768 0.0808 0.0182 0.1768 0.3354 0.1768 0
##
## Random effects:
## Name Model
## COMPRECNTDEN RW1 model
## myID IID model
## Int.DATE_YR IID model
## SexAgeGroup.f IID model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant 0.975quant
## Precision for COMPRECNTDEN 1.216 0.3830 0.6338 1.160 2.125
## Precision for myID 2.274 0.5057 1.4553 2.213 3.431
## Precision for Int.DATE_YR 1.646 0.4950 0.8804 1.578 2.811
## Precision for SexAgeGroup.f 1.099 0.3547 0.5637 1.046 1.944
## mode
## Precision for COMPRECNTDEN 1.0554
## Precision for myID 2.0948
## Precision for Int.DATE_YR 1.4499
## Precision for SexAgeGroup.f 0.9474
##
## Expected number of effective parameters(std dev): 129.40(0.00)
## Number of equivalent replicates : 5.665
##
## Deviance Information Criterion (DIC) ...............: 1391.59
## Deviance Information Criterion (DIC, saturated) ....: 1008.47
## Effective number of parameters .....................: 123.99
##
## Watanabe-Akaike information criterion (WAIC) ...: 1383.10
## Effective number of parameters .................: 101.13
##
## Marginal log-Likelihood: -731.71
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computedCompare DIC and WAIC
#DIC
YA.SexRatio.mod$dic$dic## [1] 1391.527YA.SexRatio.rd.mod$dic$dic## [1] 1391.421YA.SexRatio.ac.mod$dic$dic## [1] 1391.586#WAIC
YA.SexRatio.mod$waic$waic## [1] 1383.043YA.SexRatio.rd.mod$waic$waic## [1] 1382.937YA.SexRatio.ac.mod$waic$waic## [1] 1383.098Adult and Yearling
Sex Ratio (Reconstructed Population, COMRECPOP)
Zero is excluded from credible interval.
mic.d = as.data.frame(YA.SexRatio.mod$summary.random$COMRECPOP[,c(1:6)])
names(mic.d) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975")
length(unique(mic.d$ID))## [1] 13length(mic.d$ID)## [1] 13mic.d$Meane = plogis(mic.d$Mean + YA.SexRatio.mod$summary.fixed[1,1])
mic.d$Q025e = plogis(mic.d$Q025 + YA.SexRatio.mod$summary.fixed[1,1])
mic.d$Q975e = plogis(mic.d$Q975 + YA.SexRatio.mod$summary.fixed[1,1])
mic.d$ID = mic.d$ID/30.63 #Area of Complex 1
range(mic.d$ID)## [1] 8.651649 38.132550ggplot(mic.d, aes(ID, Meane)) +
geom_point(shape=1) +
geom_smooth(data= mic.d,
aes(ID, Meane),
col = "black",
se=FALSE,
span = 0.9,
method="loess",
lwd = 1) +
geom_smooth(data = mic.d,
aes(ID, Q025e),
col = "grey",
method = "loess",
span = 0.9,
se = FALSE,
linetype= "dashed") +
geom_smooth(data = mic.d, aes(ID, Q975e),
col = "grey",
method = "loess",
span = 0.9,
se = FALSE,
linetype= "dashed") +
geom_hline(yintercept = plogis(YA.SexRatio.mod$summary.fixed[1,1]),
linetype = "solid",
col = "darkgray",
size = 0.5) +
ylab("Male (probability)") +
xlab(expression(Deer~km^2)) +
ggtitle("Proportion Male (Yearling & Adult)") +
theme_classic() +
scale_x_continuous(breaks = seq(8, 40, 5), limits =c(8,40)) +
theme(axis.text=element_text(size=16),
legend.title = element_text(size=16, face="bold"),
legend.position = "none",
legend.key = element_rect(fill = "white", linetype=0),
legend.background = element_rect(fill = "white", linetype=0),
plot.title = element_text(face="bold", size = 20, hjust=0.5),
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))
Sex-Age Class Effect
Horizontal line is doe population average.
Zero is excluded from credible interval.
mic.d = as.data.frame(YA.SexRatio.mod$summary.random$SexAgeGroup.f)
names(mic.d) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975", "mode", "kld")
mic.d$Mean = plogis(mic.d$Mean + YA.SexRatio.mod$summary.fixed[1,1])
mic.d$Q025 = plogis(mic.d$Q025 + YA.SexRatio.mod$summary.fixed[1,1])
mic.d$Q975 = plogis(mic.d$Q975 + YA.SexRatio.mod$summary.fixed[1,1])
range(mic.d$Mean)## [1] 0.4985906 0.5894887ggplot(mic.d, aes(factor(ID), y=Mean)) +
geom_point(size=2, pch=19, col = "red") +
geom_linerange(aes(ymin=Q025, ymax=Q975), colour="black") +
geom_hline(yintercept = plogis(YA.SexRatio.mod$summary.fixed[1,1]),
linetype = "solid",
col = "darkgray",
size = 0.5) +
theme_classic() +
xlab(" ") +
ylab("Male (probability)") +
ggtitle("Proportion Male (Yearling & Adult)") +
scale_x_discrete(labels = c("Yearling (f)", "Adult (f)")) +
#scale_y_continuous(breaks = seq(0, 0.7, 0.05), limits = c(-0.05, 0.7)) +
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))
2.6 Fecundity
2.6.1 All Age Groups
Organize data.
length(Doe.df$FETUS2)## [1] 2239length(which(is.na(Doe.df$FETUS2)))## [1] 1069length(which(Doe.df$FETUS2 == 0))## [1] 312length(which(Doe.df$FETUS2 >= 0))## [1] 1170Doe3.df = Doe.df %>% filter(is.na(Doe.df$FETUS2) == FALSE)
hist(Doe3.df$FETUS2, col="lightgray")
#Unique ID by Deer
Doe3.df$myID = 1:nrow(Doe3.df)
#set up model data as a list object
## Not all covariates are used.
FE.lst4 = list(list(intercept1 = rep(1, dim(Doe3.df)[1])),
list(COMPREHARKM = Doe3.df[,"COMPREHARKM"],
COMRECPOP = Doe3.df[,"COMRECPOP"],
COMRECDOEPOP = Doe3.df[,"COMRECDOEPOP"],
COMCUMPRE2HARKM = Doe3.df[,"COMCUMPRE2HARKM"],
TOTAL = Doe3.df[,"TOTAL"],
LiveWtKg = Doe3.df[,"LiveWtKg"],
Lat = Doe3.df[,"Lat"],
NN = Doe3.df[,"NN"],
Area = Doe3.df[,"CmpxArea"],
HINDFOOT = Doe3.df[,"HINDFOOT"]/100,
COMPREDOEHARKM = Doe3.df[,"COMPREDOEHARKM"],
COMCUMPRE2DOEHARKM = Doe3.df[,"COMCUMPRE2DOEHARKM"],
COMPRECNTDEN = Doe3.df[,"COMPRECNTDEN"],
GrowSeas = round(Doe3.df[,"GrowSeas"],3),
GrowSeas.sc = Doe3.df[,"GrowSeas.sc"],
PrecEx10mm = Doe3.df[,"PrecEx10mm"],
PrecEx20mm = Doe3.df[,"PrecEx20mm"],
DSubZero = Doe3.df[,"DSubZero"],
Prec95pct = Doe3.df[,"Prec95pct"],
MonMxDT = Doe3.df[,"MonMxDT"],
MMx1DPrec = Doe3.df[,"MMx1DPrec"],
AMinT = Doe3.df[,"AMinT"],
AMxT = Doe3.df[,"AMxT"],
SumDur = round(Doe3.df[,"SumDur"],3),
SumDur.sc = Doe3.df[,"SumDur.sc"],
FrstDays = Doe3.df[,"FrstDays"],
Int.DATE_YR = Doe3.df[,"Int.DATE_YR"],
RepGroups = Doe3.df[,"RepGroups"],
SexAgeGroup.f = Doe3.df[,"SexAgeGroup.f"],
AGEGROUP = Doe3.df[,"AGEGROUP"],
PREG = as.factor(Doe3.df[,"PREG"]),
f_COMPLEX = Doe3.df[,"f_COMPLEX"],
myID = Doe3.df[,"myID"],
Month = Doe3.df[,"DATE_MO"]))
##Yearling and Adult COmbined Set
Doe3B.df = Doe3.df %>% filter(SexAgeGroup.f !="21")
dim(Doe3.df)## [1] 1170 113dim(Doe3B.df)## [1] 932 113Doe3B.df$SexAgeGroup.f = factor(as.character(Doe3B.df$SexAgeGroup.f))
#Unique ID by Deer
Doe3B.df$myID = 1:nrow(Doe3B.df)
#set up model data as a list object
## Not all covariates are used.
FE.lst4B = list(list(intercept1 = rep(1, dim(Doe3B.df)[1])),
list(COMPREHARKM = Doe3B.df[,"COMPREHARKM"],
COMRECPOP = Doe3B.df[,"COMRECPOP"],
COMRECDOEPOP = Doe3B.df[,"COMRECDOEPOP"],
COMCUMPRE2HARKM = Doe3B.df[,"COMCUMPRE2HARKM"],
TOTAL = Doe3B.df[,"TOTAL"],
LiveWtKg = Doe3B.df[,"LiveWtKg"],
Lat = Doe3B.df[,"Lat"],
NN = Doe3B.df[,"NN"],
Area = Doe3B.df[,"CmpxArea"],
HINDFOOT = Doe3B.df[,"HINDFOOT"]/100,
COMPREDOEHARKM = Doe3B.df[,"COMPREDOEHARKM"],
COMCUMPRE2DOEHARKM = Doe3B.df[,"COMCUMPRE2DOEHARKM"],
COMPRECNTDEN = Doe3B.df[,"COMPRECNTDEN"],
GrowSeas = round(Doe3B.df[,"GrowSeas"],3),
GrowSeas.sc = Doe3B.df[,"GrowSeas.sc"],
PrecEx10mm = Doe3B.df[,"PrecEx10mm"],
PrecEx20mm = Doe3B.df[,"PrecEx20mm"],
DSubZero = Doe3B.df[,"DSubZero"],
Prec95pct = Doe3B.df[,"Prec95pct"],
MonMxDT = Doe3B.df[,"MonMxDT"],
MMx1DPrec = Doe3B.df[,"MMx1DPrec"],
AMinT = Doe3B.df[,"AMinT"],
AMxT = Doe3B.df[,"AMxT"],
SumDur = round(Doe3B.df[,"SumDur"],3),
SumDur.sc = Doe3B.df[,"SumDur.sc"],
FrstDays = Doe3B.df[,"FrstDays"],
Int.DATE_YR = Doe3B.df[,"Int.DATE_YR"],
RepGroups = Doe3B.df[,"RepGroups"],
SexAgeGroup.f = Doe3B.df[,"SexAgeGroup.f"],
AGEGROUP = Doe3B.df[,"AGEGROUP"],
PREG = as.factor(Doe3B.df[,"PREG"]),
f_COMPLEX = Doe3B.df[,"f_COMPLEX"],
myID = Doe3B.df[,"myID"],
Month = Doe3B.df[,"DATE_MO"]))
Fecund.stk = inla.stack(data = list(Fecund = Doe3.df$FETUS2),
A = list(1,1),
effects = FE.lst4,
tag = "base.0")
Fecund2.stk = inla.stack(data = list(Fecund = Doe3B.df$FETUS2),
A = list(1,1),
effects = FE.lst4B,
tag = "ya.0")Model for Fecundity (Reconstructed Population COMRECPOP)
Fecund.frm = Fecund ~ -1 + intercept1 +
f(COMRECPOP,
model="rw1",
constr = TRUE,
scale.model = TRUE,
replicate = RepGroups,
hyper = nPrior) +
f(myID,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(Int.DATE_YR,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(SexAgeGroup.f,
model="iid",
constr = TRUE,
hyper = nPrior)
#theta.m1 = Fecund.mod$internal.summary.hyperpar$mean
theta.m1 = c(0.5070721, 2.8727497, 0.5030679, -0.0437353)
Fecund.mod = inla(Fecund.frm,
data = inla.stack.data(Fecund.stk),
family = "poisson",
verbose = TRUE,
control.fixed = list(prec = 0.001,
prec.intercept = 0.0001),
control.predictor = list(
A = inla.stack.A(Fecund.stk),
compute = TRUE,
link = 1),
control.mode = list(restart = TRUE, theta = theta.m1),
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 = Fecund.frm, family = \"poisson\", data = inla.stack.data(Fecund.stk), ", " verbose = TRUE, control.compute = list(dic = TRUE, cpo = TRUE, ", " waic = TRUE), control.predictor = list(A = inla.stack.A(Fecund.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 = theta.m1))")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 1.5239 1.6726 0.1552 3.3518
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept1 -0.2408 0.0724 -0.3829 -0.2408 -0.0988 -0.2408 0
##
## Random effects:
## Name Model
## COMRECPOP RW1 model
## myID IID model
## Int.DATE_YR IID model
## SexAgeGroup.f IID model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant 0.975quant
## Precision for COMRECPOP 1.7321 0.5120 0.9410 1.6608 2.936
## Precision for myID 17.9312 2.9963 12.8452 17.6498 24.586
## Precision for Int.DATE_YR 1.7291 0.5334 0.9124 1.6525 2.990
## Precision for SexAgeGroup.f 0.9992 0.3028 0.5328 0.9565 1.713
## mode
## Precision for COMRECPOP 1.5273
## Precision for myID 17.0779
## Precision for Int.DATE_YR 1.5096
## Precision for SexAgeGroup.f 0.8768
##
## Expected number of effective parameters(std dev): 105.32(0.00)
## Number of equivalent replicates : 11.11
##
## Deviance Information Criterion (DIC) ...............: 2768.54
## Deviance Information Criterion (DIC, saturated) ....: 664.64
## Effective number of parameters .....................: 107.45
##
## Watanabe-Akaike information criterion (WAIC) ...: 2697.05
## Effective number of parameters .................: 33.27
##
## Marginal log-Likelihood: -1465.08
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computedModel for Fecundity (Reconstructed Doe Population COMRECDOEPOP)
Fecund.rd.frm = Fecund ~ -1 + intercept1 +
f(COMRECDOEPOP,
model="rw1",
constr = TRUE,
scale.model = TRUE,
replicate = RepGroups,
hyper = nPrior) +
f(myID,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(Int.DATE_YR,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(SexAgeGroup.f,
model="iid",
constr = TRUE,
hyper = nPrior)
#theta.m1 = Fecund.rd.mod$internal.summary.hyperpar$mean
theta.m1 = c(0.49631206, 2.87371629, 0.50501719, -0.04383419)
Fecund.rd.mod = inla(Fecund.rd.frm,
data = inla.stack.data(Fecund.stk),
family = "poisson",
verbose = TRUE,
control.fixed = list(prec = 0.001,
prec.intercept = 0.0001),
control.predictor = list(
A = inla.stack.A(Fecund.stk),
compute = TRUE,
link = 1),
control.mode = list(restart = TRUE, theta = theta.m1),
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.rd.mod)##
## Call:
## c("inla(formula = Fecund.rd.frm, family = \"poisson\", data = inla.stack.data(Fecund.stk), ", " verbose = TRUE, control.compute = list(dic = TRUE, cpo = TRUE, ", " waic = TRUE), control.predictor = list(A = inla.stack.A(Fecund.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 = theta.m1))")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 1.4790 2.0143 0.1812 3.6746
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept1 -0.2428 0.074 -0.3881 -0.2428 -0.0975 -0.2428 0
##
## Random effects:
## Name Model
## COMRECDOEPOP RW1 model
## myID IID model
## Int.DATE_YR IID model
## SexAgeGroup.f IID model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant 0.975quant
## Precision for COMRECDOEPOP 1.7079 0.5036 0.9275 1.6385 2.891
## Precision for myID 17.9391 2.9977 12.8522 17.6568 24.599
## Precision for Int.DATE_YR 1.7312 0.5345 0.9125 1.6546 2.995
## Precision for SexAgeGroup.f 0.9988 0.3026 0.5325 0.9562 1.712
## mode
## Precision for COMRECDOEPOP 1.5083
## Precision for myID 17.0828
## Precision for Int.DATE_YR 1.5117
## Precision for SexAgeGroup.f 0.8766
##
## Expected number of effective parameters(std dev): 105.28(0.00)
## Number of equivalent replicates : 11.11
##
## Deviance Information Criterion (DIC) ...............: 2767.69
## Deviance Information Criterion (DIC, saturated) ....: 663.79
## Effective number of parameters .....................: 107.41
##
## Watanabe-Akaike information criterion (WAIC) ...: 2696.34
## Effective number of parameters .................: 33.38
##
## Marginal log-Likelihood: -1467.14
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computedModel for Fecundity (Aerial Counts COMPRECNTDEN)
Fecund.ac.frm = Fecund ~ -1 + intercept1 +
f(COMPRECNTDEN,
model="rw1",
constr = TRUE,
scale.model = TRUE,
replicate = RepGroups,
hyper = nPrior) +
f(myID,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(Int.DATE_YR,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(SexAgeGroup.f,
model="iid",
constr = TRUE,
hyper = nPrior)
#theta.m1 = Fecund.ac.mod$internal.summary.hyperpar$mean
theta.m1 = c(0.52886462, 2.87382500, 0.47945051, -0.04602817)
Fecund.ac.mod = inla(Fecund.ac.frm,
data = inla.stack.data(Fecund.stk),
family = "poisson",
verbose = TRUE,
control.fixed = list(prec = 0.001,
prec.intercept = 0.0001),
control.predictor = list(
A = inla.stack.A(Fecund.stk),
compute = TRUE,
link = 1),
control.mode = list(restart = TRUE, theta = theta.m1),
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.ac.mod)##
## Call:
## c("inla(formula = Fecund.ac.frm, family = \"poisson\", data = inla.stack.data(Fecund.stk), ", " verbose = TRUE, control.compute = list(dic = TRUE, cpo = TRUE, ", " waic = TRUE), control.predictor = list(A = inla.stack.A(Fecund.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 = theta.m1))")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 1.5100 2.0455 0.1596 3.7151
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept1 -0.2531 0.0748 -0.3999 -0.2531 -0.1064 -0.2531 0
##
## Random effects:
## Name Model
## COMPRECNTDEN RW1 model
## myID IID model
## Int.DATE_YR IID model
## SexAgeGroup.f IID model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant 0.975quant
## Precision for COMPRECNTDEN 1.7631 0.5111 0.9682 1.694 2.962
## Precision for myID 17.9405 2.9981 12.8519 17.659 24.600
## Precision for Int.DATE_YR 1.6838 0.5122 0.8954 1.611 2.892
## Precision for SexAgeGroup.f 0.9964 0.3017 0.5316 0.954 1.707
## mode
## Precision for COMPRECNTDEN 1.5629
## Precision for myID 17.0860
## Precision for Int.DATE_YR 1.4765
## Precision for SexAgeGroup.f 0.8747
##
## Expected number of effective parameters(std dev): 106.15(0.00)
## Number of equivalent replicates : 11.02
##
## Deviance Information Criterion (DIC) ...............: 2769.82
## Deviance Information Criterion (DIC, saturated) ....: 665.92
## Effective number of parameters .....................: 108.29
##
## Watanabe-Akaike information criterion (WAIC) ...: 2697.30
## Effective number of parameters .................: 33.09
##
## Marginal log-Likelihood: -1484.86
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computedCompare DIC and WAIC
#DIC
Fecund.mod$dic$dic## [1] 2768.544Fecund.rd.mod$dic$dic## [1] 2767.689Fecund.ac.mod$dic$dic## [1] 2769.825#WAIC
Fecund.mod$waic$waic## [1] 2697.055Fecund.rd.mod$waic$waic## [1] 2696.34Fecund.ac.mod$waic$waic## [1] 2697.299Fecundity (Reconstructed Population, COMRECPOP)
Zero is excluded from credible interval.
mic.d = as.data.frame(Fecund.mod$summary.random$COMRECPOP[,c(1:6)])
names(mic.d) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975")
length(unique(mic.d$ID))## [1] 13length(mic.d$ID)## [1] 39mic.d$ID2 = as.factor(
rep(c("Fawn (f)", "Yearling (f)", "Adult (f)"),
each = length(unique(mic.d$ID))))
mic.d$Meane = exp(mic.d$Mean + Fecund.mod$summary.fixed[1,1])
mic.d$Q025e = exp(mic.d$Q025 + Fecund.mod$summary.fixed[1,1])
mic.d$Q975e = exp(mic.d$Q975 + Fecund.mod$summary.fixed[1,1])
mic.d$ID = mic.d$ID/30.63 #Area of Complex 1
ggplot(mic.d, aes(ID, Meane, group = ID2)) +
geom_point(shape=1) +
geom_smooth(data= mic.d,
aes(ID, Meane,
group = mic.d$ID2,
col = factor(mic.d$ID2),
linetype= factor(mic.d$ID2)),
se=FALSE,
span = 0.9,
method="loess",
lwd = 1) +
geom_hline(yintercept = exp(Fecund.mod$summary.fixed[1,1]),
linetype = "solid",
col = "darkgray",
size = 0.5) +
scale_linetype_manual(name = "Age Group",
values=c("solid", "longdash", "dotted")) +
scale_colour_manual(name = "Age Group",
values=c("red", "brown", "blue")) +
ylab("Fecundity (Fetus Count)") +
xlab(expression(Deer~km^2)) +
ggtitle("Fecundity (COMRECPOP)") +
theme_classic() +
scale_y_continuous(breaks = seq(0.25, 2, 0.1), limits = c(0.2, 2)) +
scale_x_continuous(breaks = seq(8, 40, 5), limits =c(8,40)) +
theme(axis.text=element_text(size=16),
legend.title = element_text(size=16, face="bold"),
legend.position = c(0.7,0.7),
legend.key = element_rect(fill = "white", linetype=0),
legend.background = element_rect(fill = "white", linetype=0),
plot.title = element_text(face="bold", size = 20, hjust=0.5),
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))
Fecundity (Aerial Counts, COMPRECNTDEN)
mic.d = as.data.frame(Fecund.ac.mod$summary.random$COMPRECNTDEN[,c(1:6)])
names(mic.d) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975")
length(unique(mic.d$ID))## [1] 17length(mic.d$ID)## [1] 51mic.d$ID2 = as.factor(
rep(c("Fawn (f)", "Yearling (f)", "Adult (f)"),
each = length(unique(mic.d$ID))))
mic.d$Meane = exp(mic.d$Mean + Fecund.mod$summary.fixed[1,1])
mic.d$Q025e = exp(mic.d$Q025 + Fecund.mod$summary.fixed[1,1])
mic.d$Q975e = exp(mic.d$Q975 + Fecund.mod$summary.fixed[1,1])
ggplot(mic.d, aes(ID, Meane, group = ID2)) +
geom_point(shape=1) +
geom_smooth(data= mic.d,
aes(ID, Meane,
group = mic.d$ID2,
col = factor(mic.d$ID2),
linetype= factor(mic.d$ID2)),
se=FALSE,
span = 0.9,
method="loess",
lwd = 1) +
geom_hline(yintercept = exp(Fecund.mod$summary.fixed[1,1]),
linetype = "solid",
col = "darkgray",
size = 0.5) +
scale_linetype_manual(name = "Age Group",
values=c("solid", "longdash", "dotted")) +
scale_colour_manual(name = "Age Group",
values=c("red", "brown", "blue")) +
ylab("Fecundity (Fetus Count)") +
xlab(expression(Deer~km^2)) +
ggtitle("Fecundity (COMPRECNTDEN)") +
theme_classic() +
#scale_y_continuous(breaks = seq(0.25, 2, 0.1), limits = c(0.2, 2)) +
scale_x_continuous(breaks = seq(4, 40, 5), limits =c(4,40)) +
theme(axis.text=element_text(size=16),
legend.title = element_text(size=16, face="bold"),
legend.position = c(0.85,0.2),
legend.key = element_rect(fill = "white", linetype=0),
legend.background = element_rect(fill = "white", linetype=0),
plot.title = element_text(face="bold", size = 20, hjust=0.5),
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))
Model Fitted Estimates. Note that plots above are the density effect only, WHILE holding all other covariates constant. In actuality, the other covariates (Year, SexAgeGroups, and individuals (MyID) vary). To understand estimates consoidering all covariates, we need to look at the fitted values.
Temp.data = Doe3.df
idat = inla.stack.index(Fecund.stk, "base.0")$data #Indexregion data from model
Temp.data$Pred.Fecund = Fecund.mod$summary.fitted.values$mean[idat] #Fitted values for region
Temp.data$Low = Fecund.mod$summary.fitted.values$`0.025quant`[idat]
Temp.data$High = Fecund.mod$summary.fitted.values$`0.975quant`[idat]
Temp.data = Temp.data %>% filter(is.na(FETUS) == FALSE)
Fecund.fit = as.data.frame(
Temp.data %>%
group_by(SexAgeGroup.f) %>%
summarise(Mean = mean(Pred.Fecund, na.rm=T),
Low = mean(Low, na.rm=T),
High = mean(High, na.rm=T)))
Fecund.fit$Age = c("Fawn (f)", "Yearling (f)", "Adult (f)")
Fecund.fit$Age = ordered(factor(Fecund.fit$Age), levels = c("Fawn (f)", "Yearling (f)", "Adult (f)"))
ggplot(Fecund.fit, aes(factor(Age),Mean)) +
geom_bar(stat = "identity", fill = "lightgray") +
geom_point(size=2, pch=19, col = "red") +
geom_linerange(aes(ymin=Low, ymax=High), colour="black") +
theme_classic() +
ylab("Mean Fecundity (Fetus Count)") +
xlab( "") +
ggtitle("Population Mean Fecundity (COMRECPOP)") +
scale_x_discrete(labels = c("Fawn (f)", "Yearling (f)", "Adult (f)")) +
scale_y_continuous(breaks = seq(0, 3, 0.5), limits = c(0, 3)) +
theme(axis.text=element_text(size=16),
legend.title = element_text(size=16, face="bold"),
legend.position = c(0.7,0.7),
legend.key = element_rect(fill = "white", linetype=0),
legend.background = element_rect(fill = "white", linetype=0),
plot.title = element_text(face="bold", size = 20, hjust=0.5),
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))
Temp.data$Class = ifelse(Temp.data$SexAgeGeoup == "21", "Fawn (f)",
ifelse(Temp.data$SexAgeGeoup == "22", "Yearling (f)",
ifelse(Temp.data$SexAgeGeoup == "23", "Adult (f)", NA)))
Temp.data$Class = ordered(factor(Temp.data$Class, levels = c("Fawn (f)", "Yearling (f)", "Adult (f)")))
Temp.data$Density = Temp.data$COMRECPOP/30.63
COMRECPOP.plt = Temp.data %>% select(Class, Pred.Fecund)
hline_dat1 = as.data.frame(
COMRECPOP.plt %>%
group_by(Class) %>%
summarise(Med = mean(Pred.Fecund)))
ggplot(Temp.data, aes(Density, Pred.Fecund)) +
geom_point(shape=1) +
geom_smooth(method="loess",
se=FALSE,
col="black",
lwd = 1) +
facet_wrap(~Class) +
geom_hline(data=hline_dat1,
aes(yintercept=Med),
col = "darkgray",
size = 0.5) +
ylab("Fecundity (Fetus Count)") +
xlab(expression(Deer~km^2)) +
ggtitle("Estimated Fecundity (Complex 1) \n COMRECPOP") +
theme_classic() +
theme(axis.text=element_text(size=16),
plot.title = element_text(face="bold", size = 20, hjust=0.5),
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))
Sex-Age Class Effect
Horizontal line is doe population average.
Zero is excluded from credible interval.
mic.d = as.data.frame(Fecund.mod$summary.random$SexAgeGroup.f)
names(mic.d) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975", "mode", "kld")
mic.d$Mean = exp(mic.d$Mean + Fecund.mod$summary.fixed[1,1])
mic.d$Q025 = exp(mic.d$Q025 + Fecund.mod$summary.fixed[1,1])
mic.d$Q975 = exp(mic.d$Q975 + Fecund.mod$summary.fixed[1,1])
range(mic.d$Mean)## [1] 0.2057161 1.6412071ggplot(mic.d, aes(factor(ID), y=Mean)) +
geom_point(size=2, pch=19, col = "red") +
geom_linerange(aes(ymin=Q025, ymax=Q975), colour="black") +
geom_hline(yintercept = plogis(Fecund.mod$summary.fixed[1,1]),
linetype = "solid",
col = "darkgray",
size = 0.5) +
theme_classic() +
xlab(" ") +
ylab("Fecundity (Fetus Count)") +
ggtitle("Age Effect") +
scale_x_discrete(labels = c("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))
2.6.2 Yearling and Adult Combined
Model for Yearling and Adult Fecundity (Reconstructed Population COMRECPOP)
YA.Fecund.frm = Fecund ~ -1 + intercept1 +
f(COMRECPOP,
model="rw1",
constr = TRUE,
scale.model = TRUE,
hyper = nPrior) +
f(myID,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(Int.DATE_YR,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(SexAgeGroup.f,
model="iid",
constr = TRUE,
hyper = nPrior)
#theta.m1 = YA.Fecund.mod$internal.summary.hyperpar$mean
theta.m1 = c(0.5070721, 2.8727497, 0.5030679, -0.0437353)
YA.Fecund.mod = inla(YA.Fecund.frm,
data = inla.stack.data(Fecund2.stk),
family = "poisson",
verbose = TRUE,
control.fixed = list(prec = 0.001,
prec.intercept = 0.0001),
control.predictor = list(
A = inla.stack.A(Fecund2.stk),
compute = TRUE,
link = 1),
control.mode = list(restart = TRUE, theta = theta.m1),
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(YA.Fecund.mod)##
## Call:
## c("inla(formula = YA.Fecund.frm, family = \"poisson\", data = inla.stack.data(Fecund2.stk), ", " verbose = TRUE, control.compute = list(dic = TRUE, cpo = TRUE, ", " waic = TRUE), control.predictor = list(A = inla.stack.A(Fecund2.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 = theta.m1))")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 1.5059 1.6097 0.1317 3.2473
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept1 0.4345 0.048 0.3401 0.4345 0.5287 0.4345 0
##
## Random effects:
## Name Model
## COMRECPOP RW1 model
## myID IID model
## Int.DATE_YR IID model
## SexAgeGroup.f IID model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant 0.975quant
## Precision for COMRECPOP 1.234 0.3873 0.6433 1.178 2.152
## Precision for myID 18.008 3.0100 12.9098 17.720 24.713
## Precision for Int.DATE_YR 1.634 0.4999 0.8623 1.564 2.812
## Precision for SexAgeGroup.f 1.105 0.3577 0.5635 1.052 1.956
## mode
## Precision for COMRECPOP 1.0735
## Precision for myID 17.1314
## Precision for Int.DATE_YR 1.4341
## Precision for SexAgeGroup.f 0.9535
##
## Expected number of effective parameters(std dev): 91.83(0.00)
## Number of equivalent replicates : 10.15
##
## Deviance Information Criterion (DIC) ...............: 2556.83
## Deviance Information Criterion (DIC, saturated) ....: 518.00
## Effective number of parameters .....................: 93.35
##
## Watanabe-Akaike information criterion (WAIC) ...: 2490.11
## Effective number of parameters .................: 24.58
##
## Marginal log-Likelihood: -1330.54
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computedModel for Yearling and Adult Fecundity (Reconstructed Doe Population COMRECDOEPOP)
YA.Fecund.rd.frm = Fecund ~ -1 + intercept1 +
f(COMRECDOEPOP,
model="rw1",
constr = TRUE,
scale.model = TRUE,
hyper = nPrior) +
f(myID,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(Int.DATE_YR,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(SexAgeGroup.f,
model="iid",
constr = TRUE,
hyper = nPrior)
#theta.m1 = YA.Fecund.rd.mod$internal.summary.hyperpar$mean
theta.m1 = c(0.49631206, 2.87371629, 0.50501719, -0.04383419)
YA.Fecund.rd.mod = inla(YA.Fecund.rd.frm,
data = inla.stack.data(Fecund2.stk),
family = "poisson",
verbose = TRUE,
control.fixed = list(prec = 0.001,
prec.intercept = 0.0001),
control.predictor = list(
A = inla.stack.A(Fecund2.stk),
compute = TRUE,
link = 1),
control.mode = list(restart = TRUE, theta = theta.m1),
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(YA.Fecund.rd.mod)##
## Call:
## c("inla(formula = YA.Fecund.rd.frm, family = \"poisson\", data = inla.stack.data(Fecund2.stk), ", " verbose = TRUE, control.compute = list(dic = TRUE, cpo = TRUE, ", " waic = TRUE), control.predictor = list(A = inla.stack.A(Fecund2.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 = theta.m1))")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 1.5359 1.5458 0.1466 3.2284
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept1 0.4348 0.0479 0.3407 0.4348 0.5288 0.4348 0
##
## Random effects:
## Name Model
## COMRECDOEPOP RW1 model
## myID IID model
## Int.DATE_YR IID model
## SexAgeGroup.f IID model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant 0.975quant
## Precision for COMRECDOEPOP 1.221 0.3837 0.6362 1.165 2.130
## Precision for myID 18.012 3.0091 12.9091 17.727 24.704
## Precision for Int.DATE_YR 1.649 0.5054 0.8690 1.579 2.840
## Precision for SexAgeGroup.f 1.104 0.3580 0.5644 1.050 1.957
## mode
## Precision for COMRECDOEPOP 1.0615
## Precision for myID 17.1456
## Precision for Int.DATE_YR 1.4474
## Precision for SexAgeGroup.f 0.9505
##
## Expected number of effective parameters(std dev): 91.74(0.00)
## Number of equivalent replicates : 10.16
##
## Deviance Information Criterion (DIC) ...............: 2556.70
## Deviance Information Criterion (DIC, saturated) ....: 517.87
## Effective number of parameters .....................: 93.26
##
## Watanabe-Akaike information criterion (WAIC) ...: 2490.04
## Effective number of parameters .................: 24.56
##
## Marginal log-Likelihood: -1331.18
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computedModel for Yearling and Adult Fecundity (Aerial Counts COMPRECNTDEN)
YA.Fecund.ac.frm = Fecund ~ -1 + intercept1 +
f(COMPRECNTDEN,
model="rw1",
constr = TRUE,
scale.model = TRUE,
hyper = nPrior) +
f(myID,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(Int.DATE_YR,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(SexAgeGroup.f,
model="iid",
constr = TRUE,
hyper = nPrior)
#theta.m1 = YA.Fecund.ac.mod$internal.summary.hyperpar$mean
theta.m1 = c(0.52886462, 2.87382500, 0.47945051, -0.04602817)
YA.Fecund.ac.mod = inla(YA.Fecund.ac.frm,
data = inla.stack.data(Fecund2.stk),
family = "poisson",
verbose = TRUE,
control.fixed = list(prec = 0.001,
prec.intercept = 0.0001),
control.predictor = list(
A = inla.stack.A(Fecund2.stk),
compute = TRUE,
link = 1),
control.mode = list(restart = TRUE, theta = theta.m1),
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(YA.Fecund.ac.mod)##
## Call:
## c("inla(formula = YA.Fecund.ac.frm, family = \"poisson\", data = inla.stack.data(Fecund2.stk), ", " verbose = TRUE, control.compute = list(dic = TRUE, cpo = TRUE, ", " waic = TRUE), control.predictor = list(A = inla.stack.A(Fecund2.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 = theta.m1))")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 1.6675 1.5519 0.1376 3.3570
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept1 0.4302 0.0492 0.3336 0.4302 0.5267 0.4302 0
##
## Random effects:
## Name Model
## COMPRECNTDEN RW1 model
## myID IID model
## Int.DATE_YR IID model
## SexAgeGroup.f IID model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant 0.975quant
## Precision for COMPRECNTDEN 1.227 0.3851 0.6394 1.171 2.139
## Precision for myID 18.014 3.0094 12.9107 17.728 24.707
## Precision for Int.DATE_YR 1.616 0.4913 0.8558 1.548 2.772
## Precision for SexAgeGroup.f 1.103 0.3577 0.5642 1.050 1.956
## mode
## Precision for COMPRECNTDEN 1.067
## Precision for myID 17.147
## Precision for Int.DATE_YR 1.421
## Precision for SexAgeGroup.f 0.950
##
## Expected number of effective parameters(std dev): 91.79(0.00)
## Number of equivalent replicates : 10.15
##
## Deviance Information Criterion (DIC) ...............: 2556.60
## Deviance Information Criterion (DIC, saturated) ....: 517.77
## Effective number of parameters .....................: 93.32
##
## Watanabe-Akaike information criterion (WAIC) ...: 2489.89
## Effective number of parameters .................: 24.56
##
## Marginal log-Likelihood: -1337.00
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computedCompare DIC and WAIC
#DIC
YA.Fecund.mod$dic$dic## [1] 2556.83YA.Fecund.rd.mod$dic$dic## [1] 2556.701YA.Fecund.ac.mod$dic$dic## [1] 2556.6#WAIC
YA.Fecund.mod$waic$waic## [1] 2490.109YA.Fecund.rd.mod$waic$waic## [1] 2490.04YA.Fecund.ac.mod$waic$waic## [1] 2489.894Yearling and Adult
Fecundity (Reconstructed Population, COMRECPOP)
Zero is excluded from credible interval.
mic.d = as.data.frame(YA.Fecund.mod$summary.random$COMRECPOP[,c(1:6)])
names(mic.d) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975")
length(unique(mic.d$ID))## [1] 13length(mic.d$ID)## [1] 13mic.d$Meane = exp(mic.d$Mean + YA.Fecund.mod$summary.fixed[1,1])
mic.d$Q025e = exp(mic.d$Q025 + YA.Fecund.mod$summary.fixed[1,1])
mic.d$Q975e = exp(mic.d$Q975 + YA.Fecund.mod$summary.fixed[1,1])
mic.d$ID = mic.d$ID/30.63 #Area of Complex 1
ggplot(mic.d, aes(ID, Meane)) +
geom_point(shape=1) +
geom_smooth(data= mic.d,
aes(ID, Meane),
col = "black",
se=FALSE,
span = 0.9,
method="loess",
lwd = 1) +
geom_smooth(data = mic.d,
aes(ID, Q025e),
col = "grey",
method = "loess",
span = 0.9,
se = FALSE,
linetype= "dashed") +
geom_smooth(data = mic.d, aes(ID, Q975e),
col = "grey",
method = "loess",
span = 0.9,
se = FALSE,
linetype= "dashed") +
geom_hline(yintercept = exp(YA.Fecund.mod$summary.fixed[1,1]),
linetype = "solid",
col = "darkgray",
size = 0.5) +
ylab("Fecundity (Fetus count)") +
xlab(expression(Deer~km^2)) +
ggtitle("Fecundity (Yearling & Adult)") +
theme_classic() +
scale_x_continuous(breaks = seq(8, 40, 5), limits =c(8,40)) +
theme(axis.text=element_text(size=16),
legend.title = element_text(size=16, face="bold"),
legend.position = "none",
legend.key = element_rect(fill = "white", linetype=0),
legend.background = element_rect(fill = "white", linetype=0),
plot.title = element_text(face="bold", size = 20, hjust=0.5),
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))
Yearling and Adult
Sex-Age Class Effect
Horizontal line is doe population average.
Zero is excluded from credible interval.
mic.d = as.data.frame(YA.Fecund.mod$summary.random$SexAgeGroup.f)
names(mic.d) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975", "mode", "kld")
mic.d$Mean = exp(mic.d$Mean + YA.Fecund.mod$summary.fixed[1,1])
mic.d$Q025 = exp(mic.d$Q025 + YA.Fecund.mod$summary.fixed[1,1])
mic.d$Q975 = exp(mic.d$Q975 + YA.Fecund.mod$summary.fixed[1,1])
range(mic.d$Mean)## [1] 1.426238 1.671794ggplot(mic.d, aes(factor(ID), y=Mean)) +
geom_point(size=2, pch=19, col = "red") +
geom_linerange(aes(ymin=Q025, ymax=Q975), colour="black") +
geom_hline(yintercept = exp(YA.Fecund.mod$summary.fixed[1,1]),
linetype = "solid",
col = "darkgray",
size = 0.5) +
theme_classic() +
xlab(" ") +
ylab("Fecundity (Fetus Count)") +
ggtitle("Age Effect (Yearling and Adult") +
scale_x_discrete(labels = c("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))
3 Prediction (validation)
Out of Sample Prediction.
Models Fit from Complex 1 data are used to predict rates from the regional species pool (all deer records except those from Complex 1). Note that density was “significant” for all models above as judged by credible intervals, EXCEPT in the case of Sex Ratio and Fecundity models combining Yearling and Adult Does.
3.1 Predicted Whole Weight
Orginize Regional Pool Data (Exclude Complex 1)
Pool.df = Deer.Models.Data %>% filter(COMPLEX != 1)# & COMPLEX != 1 = Exclude Complex 1
Pool.df$COMPRECNTDEN2 = ifelse(Pool.df$COMCOMCNT == 2, NA, Pool.df$COMPRECNTDEN) #Complete aerial survey
Pool.df = Pool.df %>%
filter(is.na(AGEGROUP) == FALSE) #Remove records w/ missing age
unique(Pool.df$SexAgeGroup.f)## [1] 11 12 13 21 22 23
## Levels: 11 12 13 21 22 23Pool.df$RepGroups = as.integer(Pool.df$SexAgeGroup.f)
unique(Pool.df$RepGroups)## [1] 1 2 3 4 5 6#Unique ID
Pool.df$myID = paste(seq(1:nrow(Pool.df)),"P", sep="_")
#set up model data as a list object
Pool.lst1 = list(list(intercept1 = rep(1, dim(Pool.df)[1])),
list(COMPREHARKM = Pool.df[,"COMPREHARKM"],
COMRECPOP = Pool.df[,"COMRECPOP"],
COMRECDOEPOP = Pool.df[,"COMRECDOEPOP"],
COMCUMPRE2HARKM = Pool.df[,"COMCUMPRE2HARKM"],
TOTAL = Pool.df[,"TOTAL"],
Lat = Pool.df[,"Lat"],
Area = Pool.df[,"CmpxArea"],
NN = Pool.df[,"NN"],
HINDFOOT = Pool.df[,"HINDFOOT"]/100,
COMPREDOEHARKM = Pool.df[,"COMPREDOEHARKM"],
COMCUMPRE2DOEHARKM = Pool.df[,"COMCUMPRE2DOEHARKM"],
COMPRECNTDEN = Pool.df[,"COMPRECNTDEN"],
GrowSeas = round(Pool.df[,"GrowSeas"],3),
GrowSeas.sc = Pool.df[,"GrowSeas.sc"],
PrecEx10mm = Pool.df[,"PrecEx10mm"],
PrecEx20mm = Pool.df[,"PrecEx20mm"],
DSubZero = Pool.df[,"DSubZero"],
Prec95pct = Pool.df[,"Prec95pct"],
MonMxDT = Pool.df[,"MonMxDT"],
MMx1DPrec = Pool.df[,"MMx1DPrec"],
AMinT = Pool.df[,"AMinT"],
AMxT = Pool.df[,"AMxT"],
SumDur = round(Pool.df[,"SumDur"],3),
SumDur.sc = Pool.df[,"SumDur.sc"],
FrstDays = Pool.df[,"FrstDays"],
Int.DATE_YR = Pool.df[,"Int.DATE_YR"],
RepGroups = Pool.df[,"RepGroups"],
SexAgeGroup.f = Pool.df[,"SexAgeGroup.f"],
AGEGROUP = Pool.df[,"AGEGROUP"],
PREG = as.factor(Pool.df[,"PREG"]),
f_COMPLEX = Pool.df[,"f_COMPLEX"],
myID = Pool.df[,"myID"],
Month = Pool.df[,"DATE_MO"]))
Pool.stk = inla.stack(data = list(Weight = rep(NA, dim(Pool.df)[1])), #Weight is set to NA
A = list(1,1),
effects = Pool.lst1,
tag = "pool.0")
Pool.kist.stk = inla.stack(data = list(Kistner = rep(NA, dim(Pool.df)[1])), #TOTAL Kistner Score is set to NA
A = list(1,1),
effects = Pool.lst1,
tag = "kpool.0")Combine with Complex 1 Model Data
Prediction.stk = inla.stack(WholeWeight.stk, Pool.stk) #For Whole Weight
Prediction.Kist.stk = inla.stack(Kistner.stk, Pool.kist.stk) #For Total KistnerPrediction for weight (Aerial Counts COMPRECNTDEN)
Pred.LiveWeight = Weight ~ -1 + intercept1 +
f(COMPRECNTDEN,
model="rw1",
constr = TRUE,
scale.model = TRUE,
replicate = RepGroups,
hyper = nPrior) +
f(myID,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(f_COMPLEX,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(SexAgeGroup.f,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(PREG,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(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 ) +
DSubZero + TOTAL
#theta.m3 = Pred.LiveWeight.mod$internal.summary.hyperpar$mean
theta.m3 = c(-3.864739128, -0.791401360, -0.573948340, 0.004445255, -2.868251163, 0.115819905,
-1.408524181, -1.086211376, 0.529612968, 0.232323111)
Pred.LiveWeight.mod = inla(Pred.LiveWeight,
data = inla.stack.data(Prediction.stk),
family = "gaussian",
verbose = TRUE,
control.fixed = list(prec = 0.001,
prec.intercept = 0.0001),
control.predictor = list(
A = inla.stack.A(Prediction.stk),
compute = TRUE,
link = 1),
control.mode = list(restart = TRUE, theta = theta.m3),
control.inla = list(strategy="gaussian",
int.strategy = "eb"),
control.results = list(return.marginals.random = TRUE,
return.marginals.predictor = TRUE),
control.compute=list(dic = FALSE, cpo = FALSE, waic = FALSE)) Extract Results for Regional Pool whole Weight and Compare.
idat = inla.stack.index(Prediction.stk, "pool.0")$data #Indexregion data from model
Pool.df$Pred.Weight = Pred.LiveWeight.mod$summary.fitted.values$mean[idat] #Fitted values for region
Pool.df$Low = Pred.LiveWeight.mod$summary.fitted.values$`0.025quant`[idat]
Pool.df$High = Pred.LiveWeight.mod$summary.fitted.values$`0.975quant`[idat]
Pool2.df = Pool.df %>% filter(is.na(LiveWtKg) == FALSE) #Can't compare to missing value, even though it was predicted
#Range
#Actual
range(Pool2.df$LiveWtKg)## [1] 22.95429 117.32484#Predicted
range(Pool2.df$Pred.Weight)## [1] 24.78030 87.69713#Correlation
cor(Pool2.df$Pred.Weight, Pool2.df$LiveWtKg)## [1] 0.8826331#Absolute Error by Sex Age Group
GrpMSE = as.data.frame(
Pool2.df %>%
mutate(Error = abs(Pred.Weight - LiveWtKg)))
GrpMSE2 = as.data.frame(
GrpMSE %>%
group_by(SexAgeGroup.f) %>%
summarise(Actual = mean(LiveWtKg),
Error = mean(Error)))
GrpMSE2$Percent = round(GrpMSE2$Error/GrpMSE2$Actual, 2)
GrpMSE2$Lable = c("Fawn (m)", "Yearling (m)", "Adult (m)",
"Fawn (f)", "Yearling (f)", "Adult (f)")
kable(GrpMSE2,
caption = "Whole Weight Prediction", font_size = 22, bold=T) %>%
kable_styling("striped", full_width = F) %>%
row_spec(0, font_size = 20) %>%
column_spec(1, bold = T)| SexAgeGroup.f | Actual | Error | Percent | Lable |
|---|---|---|---|---|
| 11 | 44.36459 | 5.702348 | 0.13 | Fawn (m) |
| 12 | 67.02380 | 7.057961 | 0.11 | Yearling (m) |
| 13 | 87.31139 | 10.828411 | 0.12 | Adult (m) |
| 21 | 39.97784 | 4.636660 | 0.12 | Fawn (f) |
| 22 | 59.11583 | 5.002457 | 0.08 | Yearling (f) |
| 23 | 65.29888 | 6.094993 | 0.09 | Adult (f) |
3.2 Predicted Kistner
Prediction for Total Kistner Score (Aerial Counts COMPRECNTDEN)
Kistner.pred = Kistner ~ -1 + intercept1 +
f(COMPRECNTDEN,
model="rw1",
constr = TRUE,
scale.model = TRUE,
replicate = RepGroups,
hyper = nPrior) +
f(myID,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(SexAgeGroup.f,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(PREG,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(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 ) +
DSubZero + LiveWtKg
#theta.m3 = Kistner.mod.aerial$internal.summary.hyperpar$mean
theta.m3 = c(-3.949758515, -0.668447978, 0.023023690, -2.971971613, -0.026000289, 0.008109743,
-1.073993338, -0.302061119, 0.039326732)
Kistner.mod.pred = inla(Kistner.pred,
data = inla.stack.data(Prediction.Kist.stk),
family = "gaussian",
verbose = TRUE,
control.fixed = list(prec = 0.001,
prec.intercept = 0.0001),
control.predictor = list(
A = inla.stack.A(Prediction.Kist.stk),
compute = TRUE,
link = 1),
control.mode = list(restart = TRUE, theta = theta.m3),
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)) Extract Results for Regional Pool Total Kistner and Compare.
idat = inla.stack.index(Prediction.Kist.stk, "kpool.0")$data #Indexregion data from model
Pool.df$Pred.kist = Kistner.mod.pred$summary.fitted.values$mean[idat] #Fitted values for region
Pool.df$Low = Kistner.mod.pred$summary.fitted.values$`0.025quant`[idat]
Pool.df$High = Kistner.mod.pred$summary.fitted.values$`0.975quant`[idat]
Pool2.df = Pool.df %>% filter(is.na(TOTAL) == FALSE) #Can't compare to missing value, even though it was predicted
#Range
#Actual
range(Pool2.df$TOTAL)## [1] 0 100#Predicted
range(Pool2.df$Pred.kist)## [1] 7.82504 80.98510#Correlation
cor(Pool2.df$Pred.kist, Pool2.df$TOTAL)## [1] 0.6298667#Absolute Error by Sex Age Group
GrpMSE = as.data.frame(
Pool2.df %>%
mutate(Error = abs(Pred.kist - TOTAL)))
GrpMSE2 = as.data.frame(
GrpMSE %>%
group_by(SexAgeGroup.f) %>%
summarise(Actual = mean(TOTAL),
Error = mean(Error)))
GrpMSE2$Percent = round(GrpMSE2$Error/GrpMSE2$Actual, 2)
GrpMSE2$Lable = c("Fawn (m)", "Yearling (m)", "Adult (m)",
"Fawn (f)", "Yearling (f)", "Adult (f)")
kable(GrpMSE2,
caption = "Total Kistner Prediction", font_size = 22, bold=T) %>%
kable_styling("striped", full_width = F) %>%
row_spec(0, font_size = 20) %>%
column_spec(1, bold = T)| SexAgeGroup.f | Actual | Error | Percent | Lable |
|---|---|---|---|---|
| 11 | 57.85408 | 16.11042 | 0.28 | Fawn (m) |
| 12 | 55.19517 | 18.21094 | 0.33 | Yearling (m) |
| 13 | 50.81854 | 21.22222 | 0.42 | Adult (m) |
| 21 | 62.74709 | 16.71469 | 0.27 | Fawn (f) |
| 22 | 79.38871 | 16.20263 | 0.20 | Yearling (f) |
| 23 | 79.35360 | 20.09587 | 0.25 | Adult (f) |
3.3 Predicted Pregnancy
Organize Regional Doe Data
Pool.Doe.df = Deer.Models.Data %>% filter(COMPLEX != 1 & SEX == 2) #Complex 1 and females only
#Yes = 1, No = 2
Pool.Doe.df$PREG2 = ifelse(Pool.Doe.df$PREG == 1, 1,
ifelse(Pool.Doe.df$PREG == 2, 0, NA))
#Missing value: Was pregnancy detectable, if yes then "0", else NA
Pool.Doe.df$PREG2 = ifelse(Pool.Doe.df$USEREPRO == 1 & is.na(Pool.Doe.df$PREG2) == TRUE, 0, Pool.Doe.df$PREG2)
#Missing value: Was pregnancy detectable, if yes then "0", else NA
Pool.Doe.df$FETUS2 = ifelse(Pool.Doe.df$USEREPRO == 1 & is.na(Pool.Doe.df$FETUS) == TRUE, 0, Pool.Doe.df$FETUS)
#COMCOMCNT = "Survey completed" If not, set record to NA
Pool.Doe.df$COMPRECNTDEN2 = ifelse(Pool.Doe.df$COMCOMCNT == 2, NA, Pool.Doe.df$COMPRECNTDEN) #Complete aerial survey
Pool.Doe.df = Pool.Doe.df %>%
filter(is.na(AGEGROUP) == FALSE) #Remove records w/ missing age
#Integer version of Sex-Age groups
unique(Pool.Doe.df$SexAgeGroup.f)## [1] 21 22 23
## Levels: 11 12 13 21 22 23Pool.Doe.df$RepGroups = as.integer(Pool.Doe.df$SexAgeGroup.f)
unique(Pool.Doe.df$RepGroups)## [1] 4 5 6Pool.Doe.df$RepGroups[Pool.Doe.df$RepGroups == 4] = 1
Pool.Doe.df$RepGroups[Pool.Doe.df$RepGroups == 5] = 2
Pool.Doe.df$RepGroups[Pool.Doe.df$RepGroups == 6] = 3
unique(Pool.Doe.df$RepGroups)## [1] 1 2 3Pool.Doe.df$SexAgeGroup.f = factor(as.character(Pool.Doe.df$SexAgeGroup.f))
#Unique ID by Deer
Pool.Doe.df$myID = paste("P", 1:nrow(Pool.Doe.df), sep=".")
#set up model data as a list object
## Not all covariates are used.
FE.lst2p = list(list(intercept1 = rep(1, dim(Pool.Doe.df)[1])),
list(COMPREHARKM = Pool.Doe.df[,"COMPREHARKM"],
COMRECPOP = Pool.Doe.df[,"COMRECPOP"],
COMRECDOEPOP = Pool.Doe.df[,"COMRECDOEPOP"],
COMCUMPRE2HARKM = Pool.Doe.df[,"COMCUMPRE2HARKM"],
TOTAL = Pool.Doe.df[,"TOTAL"],
LiveWtKg = Pool.Doe.df[,"LiveWtKg"],
Lat = Pool.Doe.df[,"Lat"],
NN = Pool.Doe.df[,"NN"],
Area = Pool.Doe.df[,"CmpxArea"],
HINDFOOT = Pool.Doe.df[,"HINDFOOT"]/100,
COMPREDOEHARKM = Pool.Doe.df[,"COMPREDOEHARKM"],
COMCUMPRE2DOEHARKM = Pool.Doe.df[,"COMCUMPRE2DOEHARKM"],
COMPRECNTDEN = Pool.Doe.df[,"COMPRECNTDEN"],
GrowSeas = round(Pool.Doe.df[,"GrowSeas"],3),
GrowSeas.sc = Pool.Doe.df[,"GrowSeas.sc"],
PrecEx10mm = Pool.Doe.df[,"PrecEx10mm"],
PrecEx20mm = Pool.Doe.df[,"PrecEx20mm"],
DSubZero = Pool.Doe.df[,"DSubZero"],
Prec95pct = Pool.Doe.df[,"Prec95pct"],
MonMxDT = Pool.Doe.df[,"MonMxDT"],
MMx1DPrec = Pool.Doe.df[,"MMx1DPrec"],
AMinT = Pool.Doe.df[,"AMinT"],
AMxT = Pool.Doe.df[,"AMxT"],
SumDur = round(Pool.Doe.df[,"SumDur"],3),
SumDur.sc = Pool.Doe.df[,"SumDur.sc"],
FrstDays = Pool.Doe.df[,"FrstDays"],
Int.DATE_YR = Pool.Doe.df[,"Int.DATE_YR"],
RepGroups = Pool.Doe.df[,"RepGroups"],
SexAgeGroup.f = Pool.Doe.df[,"SexAgeGroup.f"],
AGEGROUP = Pool.Doe.df[,"AGEGROUP"],
PREG = as.factor(Pool.Doe.df[,"PREG"]),
f_COMPLEX = Pool.Doe.df[,"f_COMPLEX"],
myID = Pool.Doe.df[,"myID"],
Month = Pool.Doe.df[,"DATE_MO"]))
#Regional Pool
Pool.Preg.stk = inla.stack(data = list(Preg = rep(NA, dim(Pool.Doe.df)[1])), #NA for dependent - to be predicted
A = list(1,1),
effects = FE.lst2p,
tag = "pool.0")Combine with COmplex 1 Model Data
Pred.Preg.stk = inla.stack(Preg.stk, Pool.Preg.stk)Model for Pregnancy (Aerial Counts COMPRECNTDEN)
Pregnancy.pred.frm = Preg ~ -1 + intercept1 +
f(COMPRECNTDEN,
model="rw1",
constr = TRUE,
scale.model = TRUE,
replicate = RepGroups,
hyper = nPrior) +
f(myID,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(SexAgeGroup.f,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(Month,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(Int.DATE_YR, #year
model="rw1",
constr = TRUE,
scale.model = TRUE,
hyper = nPrior ) +
f(SumDur,#summer duration
model="rw1",
constr = TRUE,
scale.model = TRUE,
hyper = nPrior ) +
DSubZero + TOTAL
#theta.m1 = Pregnancy.pred.mod $internal.summary.hyperpar$mean
theta.m1 = c(-0.005399492, 0.609777159, -0.510220899, -0.604797308, -0.359177299, -0.193281100)
Pregnancy.pred.mod = inla(Pregnancy.pred.frm,
data = inla.stack.data(Pred.Preg.stk),
family = "binomial",
verbose = TRUE,
control.fixed = list(prec = 0.001,
prec.intercept = 0.0001),
control.predictor = list(
A = inla.stack.A(Pred.Preg.stk),
compute = TRUE,
link = 1),
control.mode = list(restart = TRUE, theta = theta.m1),
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)) Extract Results for Regional Pool and Compare Pregnancy.
idat = inla.stack.index(Pred.Preg.stk, "pool.0")$data #Indexregion data from model
Pool.Doe.df$Pred.preg = Pregnancy.pred.mod$summary.fitted.values$mean[idat] #Fitted values for region
Pool.Doe.df$Low = Pregnancy.pred.mod$summary.fitted.values$`0.025quant`[idat]
Pool.Doe.df$High = Pregnancy.pred.mod$summary.fitted.values$`0.975quant`[idat]
Pool2.df = Pool.Doe.df %>% filter(is.na(PREG2) == FALSE) #Can't compare to missing value, even though it was predicted
#Range
#Actual
range(Pool2.df$PREG2)## [1] 0 1#Predicted
range(Pool2.df$Pred.preg)## [1] 0.0002311584 0.9985277031#Correlation
cor(Pool2.df$Pred.preg, Pool2.df$PREG2)## [1] 0.6047499Pool2.df$Brier1 = (Pool2.df$Pred.preg - Pool2.df$PREG2)^2
GrpMSE2 = as.data.frame(
Pool2.df %>%
group_by(SexAgeGroup.f) %>%
summarise(Brier = mean(Brier1)))
GrpMSE2$Lable = c("Fawn (f)", "Yearling (f)", "Adult (f)")
kable(GrpMSE2,
caption = "Pregnancy Prediction", font_size = 22, bold=T) %>%
kable_styling("striped", full_width = F) %>%
row_spec(0, font_size = 20) %>%
column_spec(1, bold = T)| SexAgeGroup.f | Brier | Lable |
|---|---|---|
| 21 | 0.1196339 | Fawn (f) |
| 22 | 0.1416930 | Yearling (f) |
| 23 | 0.1769187 | Adult (f) |
3.4 Predicted Sex Ratio
Organize Regional Data.
Pool.Doe.df$M0 = Pool.Doe.df$M
Pool.Doe.df$F0 = Pool.Doe.df$F
Pool.Doe.df$M0z = Pool.Doe.df$M0
Pool.Doe.df$M0z[is.na(Pool.Doe.df$M0z)] = 0
Pool.Doe.df$F0z = Pool.Doe.df$F0
Pool.Doe.df$F0z[is.na(Pool.Doe.df$F0z)] = 0
Pool.Doe.df$FETUS3 = Pool.Doe.df$M0z + Pool.Doe.df$F0z
Pool.Doe2.df = Pool.Doe.df
#Unique ID by Deer
Pool.Doe2.df$myID = paste("p", 1:nrow(Pool.Doe2.df), ".")
#set up model data as a list object
## Not all covariates are used.
FE.lst3p = list(list(intercept1 = rep(1, dim(Pool.Doe2.df)[1])),
list(COMPREHARKM = Pool.Doe2.df[,"COMPREHARKM"],
COMRECPOP = Pool.Doe2.df[,"COMRECPOP"],
COMRECDOEPOP = Pool.Doe2.df[,"COMRECDOEPOP"],
COMCUMPRE2HARKM = Pool.Doe2.df[,"COMCUMPRE2HARKM"],
TOTAL = Pool.Doe2.df[,"TOTAL"],
LiveWtKg = Pool.Doe2.df[,"LiveWtKg"],
Lat = Pool.Doe2.df[,"Lat"],
NN = Pool.Doe2.df[,"NN"],
Area = Pool.Doe2.df[,"CmpxArea"],
HINDFOOT = Pool.Doe2.df[,"HINDFOOT"]/100,
COMPREDOEHARKM = Pool.Doe2.df[,"COMPREDOEHARKM"],
COMCUMPRE2DOEHARKM = Pool.Doe2.df[,"COMCUMPRE2DOEHARKM"],
COMPRECNTDEN = Pool.Doe2.df[,"COMPRECNTDEN"],
GrowSeas = round(Pool.Doe2.df[,"GrowSeas"],3),
GrowSeas.sc = Pool.Doe2.df[,"GrowSeas.sc"],
PrecEx10mm = Pool.Doe2.df[,"PrecEx10mm"],
PrecEx20mm = Pool.Doe2.df[,"PrecEx20mm"],
DSubZero = Pool.Doe2.df[,"DSubZero"],
Prec95pct = Pool.Doe2.df[,"Prec95pct"],
MonMxDT = Pool.Doe2.df[,"MonMxDT"],
MMx1DPrec = Pool.Doe2.df[,"MMx1DPrec"],
AMinT = Pool.Doe2.df[,"AMinT"],
AMxT = Pool.Doe2.df[,"AMxT"],
SumDur = round(Pool.Doe2.df[,"SumDur"],3),
SumDur.sc = Pool.Doe2.df[,"SumDur.sc"],
FrstDays = Pool.Doe2.df[,"FrstDays"],
Int.DATE_YR = Pool.Doe2.df[,"Int.DATE_YR"],
RepGroups = Pool.Doe2.df[,"RepGroups"],
SexAgeGroup.f = Pool.Doe2.df[,"SexAgeGroup.f"],
AGEGROUP = Pool.Doe2.df[,"AGEGROUP"],
PREG = as.factor(Pool.Doe2.df[,"PREG"]),
f_COMPLEX = Pool.Doe2.df[,"f_COMPLEX"],
myID = Pool.Doe2.df[,"myID"],
Month = Pool.Doe2.df[,"DATE_MO"]))
pRatio.stk = inla.stack(data = list(Ratio = rep(NA, dim(Pool.Doe2.df)[1]),
Trials = dim(Pool.Doe2.df)[1]), #successes and trials set to NA
A = list(1,1),
effects = FE.lst3p,
tag = "pool.0")Combine with Complex 1 Data
Pred.Ratio.stk = inla.stack(Ratio.stk, pRatio.stk)Prediction for Sex Ratio (proportion male) (Aerial Counts COMPRECNTDEN)
SexRatio.pred.frm = Ratio ~ -1 + intercept1 +
f(COMPRECNTDEN,
model="rw1",
constr = TRUE,
scale.model = TRUE,
replicate = RepGroups,
hyper = nPrior) +
f(myID,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(Int.DATE_YR,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(SexAgeGroup.f,
model="iid",
constr = TRUE,
hyper = nPrior)
#theta.m1 = SexRatio.ac.mod$internal.summary.hyperpar$mean
theta.m1 = c(0.31736554, 0.79131261, 0.48778867, 0.07937508)
SexRatio.pred.mod = inla(SexRatio.pred.frm,
data = inla.stack.data(Pred.Ratio.stk),
family = "binomial",
verbose = TRUE,
Ntrials = inla.stack.data(Pred.Ratio.stk)$Trials, #Trials
control.fixed = list(prec = 0.001,
prec.intercept = 0.0001),
control.predictor = list(
A = inla.stack.A(Pred.Ratio.stk),
compute = TRUE,
link = 1),
control.mode = list(restart = TRUE, theta = theta.m1),
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)) Extract Results for Regional Pool and Compare Sex Ratio.
idat = inla.stack.index(Pred.Ratio.stk, "pool.0")$data #Indexregion data from model
Pool.Doe.df$Pred.ratio = SexRatio.pred.mod$summary.fitted.values$mean[idat] #Fitted values for region
Pool.Doe.df$Low = SexRatio.pred.mod$summary.fitted.values$`0.025quant`[idat]
Pool.Doe.df$High = SexRatio.pred.mod$summary.fitted.values$`0.975quant`[idat]
Pool2.df = Pool.Doe.df %>% filter(is.na(M) == FALSE) #Can't compare to missing value, even though it was predicted
Pool2.df$EstRatio = Pool2.df$M0z/Pool2.df$FETUS3
#Range
#Actual
range(Pool2.df$EstRatio)## [1] 0 1#Predicted
range(Pool2.df$Pred.ratio)## [1] 0.4267888 0.6840289#Correlation
cor(Pool2.df$Pred.ratio, Pool2.df$EstRatio)## [1] 0.009605882Pool2.df$Brier1 = (Pool2.df$Pred.ratio - Pool2.df$EstRatio)^2
GrpMSE2 = as.data.frame(
Pool2.df %>%
group_by(SexAgeGroup.f) %>%
summarise(Brier = mean(Brier1)))
GrpMSE2$Lable = c("Fawn (f)", "Yearling (f)", "Adult (f)")
kable(GrpMSE2,
caption = "Sex Ratio Prediction", font_size = 22, bold=T) %>%
kable_styling("striped", full_width = F) %>%
row_spec(0, font_size = 20) %>%
column_spec(1, bold = T)| SexAgeGroup.f | Brier | Lable |
|---|---|---|
| 21 | 0.2179128 | Fawn (f) |
| 22 | 0.1747227 | Yearling (f) |
| 23 | 0.1289388 | Adult (f) |
3.5 Predicted Fecundity
Organize data.
Doe5.df = Pool.Doe.df %>% filter(is.na(FETUS) == FALSE)
#Unique ID by Deer
Doe5.df$myID = paste("P", 1:nrow(Doe5.df), sep=".")
#set up model data as a list object
## Not all covariates are used.
FE.lst4p = list(list(intercept1 = rep(1, dim(Doe5.df)[1])),
list(COMPREHARKM = Doe5.df[,"COMPREHARKM"],
COMRECPOP = Doe5.df[,"COMRECPOP"],
COMRECDOEPOP = Doe5.df[,"COMRECDOEPOP"],
COMCUMPRE2HARKM = Doe5.df[,"COMCUMPRE2HARKM"],
TOTAL = Doe5.df[,"TOTAL"],
LiveWtKg = Doe5.df[,"LiveWtKg"],
Lat = Doe5.df[,"Lat"],
NN = Doe5.df[,"NN"],
Area = Doe5.df[,"CmpxArea"],
HINDFOOT = Doe5.df[,"HINDFOOT"]/100,
COMPREDOEHARKM = Doe5.df[,"COMPREDOEHARKM"],
COMCUMPRE2DOEHARKM = Doe5.df[,"COMCUMPRE2DOEHARKM"],
COMPRECNTDEN = Doe5.df[,"COMPRECNTDEN"],
GrowSeas = round(Doe5.df[,"GrowSeas"],3),
GrowSeas.sc = Doe5.df[,"GrowSeas.sc"],
PrecEx10mm = Doe5.df[,"PrecEx10mm"],
PrecEx20mm = Doe5.df[,"PrecEx20mm"],
DSubZero = Doe5.df[,"DSubZero"],
Prec95pct = Doe5.df[,"Prec95pct"],
MonMxDT = Doe5.df[,"MonMxDT"],
MMx1DPrec = Doe5.df[,"MMx1DPrec"],
AMinT = Doe5.df[,"AMinT"],
AMxT = Doe5.df[,"AMxT"],
SumDur = round(Doe5.df[,"SumDur"],3),
SumDur.sc = Doe5.df[,"SumDur.sc"],
FrstDays = Doe5.df[,"FrstDays"],
Int.DATE_YR = Doe5.df[,"Int.DATE_YR"],
RepGroups = Doe5.df[,"RepGroups"],
SexAgeGroup.f = Doe5.df[,"SexAgeGroup.f"],
AGEGROUP = Doe5.df[,"AGEGROUP"],
PREG = as.factor(Doe5.df[,"PREG"]),
f_COMPLEX = Doe5.df[,"f_COMPLEX"],
myID = Doe5.df[,"myID"],
Month = Doe5.df[,"DATE_MO"]))
Fecund.pred.stk = inla.stack(data = list(Fecund = rep(NA, dim(Doe5.df)[1])),
A = list(1,1),
effects = FE.lst4p,
tag = "pool.0")Combine with Complex 1 Data from Fecundity Model
Pred.Fecund.stk = inla.stack(Fecund.stk, Fecund.pred.stk)Model for Predicted Fecundity (Aerial Counts COMPRECNTDEN)
Fecund.pred.frm = Fecund ~ -1 + intercept1 +
f(COMPRECNTDEN,
model="rw1",
constr = TRUE,
scale.model = TRUE,
replicate = RepGroups,
hyper = nPrior) +
f(myID,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(Int.DATE_YR,
model="iid",
constr = TRUE,
hyper = nPrior) +
f(SexAgeGroup.f,
model="iid",
constr = TRUE,
hyper = nPrior)
#theta.m1 = Fecund.pred.mod$internal.summary.hyperpar$mean
theta.m1 = c(0.45849459, 2.87420421, 0.48666681, -0.06799408)
Fecund.pred.mod = inla(Fecund.pred.frm,
data = inla.stack.data(Pred.Fecund.stk),
family = "poisson",
verbose = TRUE,
control.fixed = list(prec = 0.001,
prec.intercept = 0.0001),
control.predictor = list(
A = inla.stack.A(Pred.Fecund.stk),
compute = TRUE,
link = 1),
control.mode = list(restart = TRUE, theta = theta.m1),
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)) Extract Results for Regional Pool Fecundity
idat = inla.stack.index(Pred.Fecund.stk, "pool.0")$data #Indexregion data from model
Doe5.df$Pred.Fecund = Fecund.pred.mod$summary.fitted.values$mean[idat] #Fitted values for region
Doe5.df$Low = Fecund.pred.mod$summary.fitted.values$`0.025quant`[idat]
Doe5.df$High = Fecund.pred.mod$summary.fitted.values$`0.975quant`[idat]
Pool2.df = Doe5.df %>% filter(is.na(FETUS) == FALSE) #Can't compare to missing value, even though it was predicted
dim(Pool2.df)## [1] 2385 118#Range
#Actual
range(Pool2.df$FETUS)## [1] 0 4#Predicted
range(Pool2.df$Pred.Fecund)## [1] 0.07076114 3.63492488#Correlation
cor(Pool2.df$Pred.Fecund, Pool2.df$FETUS)## [1] 0.2650681#Absolute Error by Sex Age Group
GrpMSE = as.data.frame(
Pool2.df %>%
mutate(Error = abs(Pred.Fecund - FETUS)))
GrpMSE2 = as.data.frame(
GrpMSE %>%
group_by(SexAgeGroup.f) %>%
summarise(Actual = mean(FETUS),
Error = mean(Error)))
GrpMSE2$Percent = round(GrpMSE2$Error/GrpMSE2$Actual, 2)
GrpMSE2$Lable = c("Fawn (f)", "Yearling (f)", "Adult (f)")
kable(GrpMSE2,
caption = "Fecundity Prediction", font_size = 22, bold=T) %>%
kable_styling("striped", full_width = F) %>%
row_spec(0, font_size = 20) %>%
column_spec(1, bold = T)| SexAgeGroup.f | Actual | Error | Percent | Lable |
|---|---|---|---|---|
| 21 | 0.9017857 | 0.7126208 | 0.79 | Fawn (f) |
| 22 | 1.6672474 | 0.5230239 | 0.31 | Yearling (f) |
| 23 | 1.8940553 | 0.5996212 | 0.32 | Adult (f) |
4 Data Glossary
MetaData = read.csv("METADATA.csv", sep = ",", header = TRUE)[1:80,] %>% select(Attribute, Description)
kable(MetaData,
caption = " ") %>%
kable_styling("striped", full_width = F) %>%
row_spec(0, font_size = 20) %>%
column_spec(1, bold = T)| Attribute | Description |
|---|---|
| ID | Unique ID assigned to each record. Deer culled from Argonne National Lab begin with the year “1995” of harvest. X10 added to all Cook County Deer so they’re unique from DuPage. |
| Deer | Count of deer |
| IDNR TAG | IL DNR leg tag number. Corresponds to an individual deer, but it was not entered into database every year |
| COOKDUPAGE | County ID (DuPage=1, Cook=2) |
| RMVYEAR1 | Removal year in which deer was harvested (e.g., Fall 2004-Spring 2005=2004; Fall 2005-Spring 2006=2005). |
| REDUCMANAG | For Complex 1, 1=1992-95 and 2=1996-2007 |
| DATE | Date deer was harvested (xx/xx/xxxx) |
| DATE_YR | Year deer was harvested |
| DATE_MO | Month deer was harvested |
| DATE_DAY | Day deer was harvested |
| JULIAN_DATE | Julian date (year, 1-366) |
| PRESERVE | Preserve location coded by number. See attached Preserve Codes tab. |
| PREAREAKM2 | Preserve area in square KM |
| PRE_AG | Was there row cropland (corn or beans) directly adjacent to the preserve (1=yes, 2=no). |
| COMPLEX | Preserve complexes that were either directly connected by deer habitat and within 5 Km from each other; or where collared deer were known to move between preserves. I have a more detailed writeup. |
| COMAREAKM2 | Complex area in square KM |
| COMP_AG | Was there row cropland (corn or beans) directly adjacent to the majority of preserves in the complex (1=yes, 2=no). |
| RMVYEAR2 | Series of years culling has occurred in a given preserve. Years are not always consecutive. |
| WFG | Was deer harvested from Waterfall Glen Forest Preserve (1=yes, 2=no). Used to distinguish deer that were harvested inside Argonne National Lab from deer harvested in Waterfall Glen which are both in LOCAT=1. |
| PREVCOUNT | Number of deer counted in preserve during aerial flights in winter previous to RMVYEAR1. Any deer culled after count were subtracted to produce a minimum population beginning in April. |
| PREVCNTDEN | PREVCOUNT/AREAKM2 |
| COMPRECNT | Number of deer counted in complex during aerial flights in winter previous to RMVYEAR1. Any deer culled after count were subtracted to produce a minimum population beginning in April. |
| COMPRECNTDEN | COMPRECNT/COMAREAKM2 |
| COMCOMCNT | A complete aerial count of all preserves in a complex was conducted for using COMPRECNT and COMPRECNTDEN (1=yes, 2=no). |
| SPOTCOUNT | Spotlight counts conducted inside Argonne National Laboratoy immediately prior to RMVYEAR1 harvest. |
| HARVEST | Total # of deer harvested from preserve in current year. |
| DOEHARV | Total # of does harvested from preserve in current year. |
| PREHARKM | Total deer/km2 culled from preserve in year previous to RMVYEAR1. |
| CUMPRE2HARKM | Total deer/km2 culled from preserve in 2-years previous to RMVYEAR1. |
| PREDOEHARKM | Total does/km2 culled from preserve in year previous to RMVYEAR1. |
| CUMPRE2DOEHARKM | Total does/km2 culled from preserve in 2-years previous to RMVYEAR1. |
| COMHARVEST | Total # of deer harvested from complex in RMVYEAR1. |
| COMDOEHARV | Total # of does harvested from complex in RMVYEAR1. |
| COMPREHARKM | Total deer/km2 culled from complex in year previous to RMVYEAR1. |
| COMCUMPRE2HARKM | Total deer/km2 culled from complex in 2-years previous to RMVYEAR1. |
| COMPREDOEHARKM | Total does/km2 culled from complex in year previous to RMVYEAR1. |
| COMCUMPRE2DOEHARKM | Total does/km2 culled from complex in 2-years previous to RMVYEAR1. |
| RECPOP | Preserve pre-cull total reconstructed population. Only exists for a few preserves. Unaged deer were added to total. |
| RECDOEPOP | Preserve pre-cull total reconstructed doe population. Only exists for a few preserves. Unaged does were added to total. |
| PRERECPOP | Preserve reconstructed population in year previous to RMVYEAR1. Only exists for a few preserves. Unaged deer were added to total. |
| PRERECDOEPOP | Preserve reconstructed doe population in year previous to RMVYEAR1. Only exists for a few preserves. Unaged does were added to total. |
| COMRECPOP | Complex pre-cull total reconstructed population. Only exists for a few preserves. Unaged deer were added to total. |
| COMDOEPOP | Complex pre-cull total reconstructed doe population. Only exists for a few preserves. Unaged does were added to total. |
| COMRECPOP | Complex reconstructed population in year previous to RMVYEAR1. Only exists for a few preserves. Unaged deer were added to total. |
| COMRECDOEPOP | Complex reconstructed doe population in year previous to RMVYEAR1. Only exists for a few preserves. Unaged does were added to total. |
| SEX | 1=male, 2=female |
| AGEWEAR | Age of deer determined by tooth replacement and wear. |
| AGECEM | Age of deer determined by cementum aging. |
| AGE | Age of deer in years (e.g., 1.5). All fawns and most yearlings were determined by tooth replacement and wear. Used AGECEM over toothwear age when it was available. |
| AGEGROUP | Fawn=1, Yearling=2, Adult=3 |
| AGECAT | Fawn = 1, Yearling = 2, 2.5-4.5 = 3, 5.5-8.5 = 4, 9.5+ = 5 |
| USEAGE | Can AGE (i.e., 5.5) be used (1=yes, 2=no). Includes all cementum aged deer, and all fawns and yearlings aged by tooth replacement and wear. |
| DRWTLBS | Weight of deer in pounds after removal of entrails. Original data. |
| DRWTKG | Weight of deer in kilograms after removal of entrails. All were converted from DRWTLBS (LBS/2.2046). |
| WWTLBS | Whole weight of deer in pounds. Original data. |
| WWTKG | Whole weight of deer in kilograms. All were converted from WWTLBS (LBS/2.2046). |
| RECWWT | Indicates a whole weight and dressed weight was recorded for the same deer. 1=YES, 2=NO |
| USEWT | For MALES and females, indicates that the weights correspond with dates used for USEREPRO (after Jan 10 for yearlings and adults; after January 31 for fawns). 1=YES, 2=NO |
| PREG | 1=yes, 2=no, blank=unknown |
| PREGUKN | Culled during USEREPRO, but pregancy status unknown (1=yes, 2=no) |
| FETUS | Total number of fetuses. Blank if unknown. |
| M | Number of male fetuses. |
| F | Number of female fetuses. |
| FU | Number of unknown fetuses |
| CLP | Number of corpus luteum of pregnancy. |
| GF | Number of Graafion follicles; indicating cycling. |
| CLE | Number of corpus luteum of estrus. |
| BRISK | Kistner fat indice score for brisket (0, 5, 10, 15). Blank if unknown. |
| OMEN | Kistner fat indice score for omentum (0, 5, 10, 15) Blank if unknown. |
| KDNY | Kistner fat indice score for kidneys (0, 5, 10, 15). Blank if unknown. |
| PERCDM | Kistner fat indice score for pericardium (0, 5, 10, 15). Blank if unknown. |
| HEART | Kistner fat indice score for heart (0, 5, 10, 15). Blank if unknown. |
| RUMP | Kistner fat indice score for rump (0, 5, 10, 15). Blank if unknown. |
| MUSCLE | Kistner fat indice score for muscle (0, 5, 10). Blank if unknown. |
| TOTAL | Sum of all Kistner fat indice scores (0-100). When scores were not available for all organs TOTAL was left blank. |
| HINDFOOT | Right hind foot length (mm) |
| USEREPRO | Defines whether doe was harvested late enough in the year to determine pregnancy (AFTER Jan 10 for yearlings and adults; AFTER January 31 for fawns). 1=YES, 2=NO. |
| USETOTAL | Total Kistner scores for males and females which include all indice scores; also corresponds with dates used for USEREPRO (after Jan 10 for yearlings and adults; after January 31 for fawns). 1=YES, 2=NO |
| RECHF | Indicates recorded hind foot length (mm); 1=YES, 2=NO |