Deer Density Dependence
1 Data Preprocessing
date() ## [1] "Mon Feb 17 12:01:26 2020"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 s_MMx1DPrec
## 1 1.000 0.000 0.001 0.000 0.000 0.005 0.000
## 2 1.158 0.003 0.006 0.001 0.001 0.000 0.005
## 3 1.279 0.000 0.000 0.002 0.005 0.000 0.001
## 4 1.598 0.010 0.005 0.000 0.000 0.002 0.019
## 5 3.940 0.001 0.033 0.000 0.001 0.121 0.083
## 6 5.634 0.000 0.048 0.004 0.002 0.035 0.003
## 7 7.429 0.088 0.199 0.004 0.076 0.025 0.004
## 8 18.905 0.066 0.002 0.374 0.632 0.810 0.144
## 9 22.542 0.832 0.706 0.615 0.283 0.000 0.743
## s_AMinT s_AMxT s_SumDur s_FrstDays
## 1 0.001 0.001 0.003 0.003
## 2 0.000 0.001 0.007 0.000
## 3 0.001 0.001 0.014 0.000
## 4 0.001 0.000 0.000 0.002
## 5 0.011 0.001 0.063 0.001
## 6 0.001 0.036 0.529 0.116
## 7 0.006 0.030 0.375 0.163
## 8 0.292 0.719 0.002 0.179
## 9 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 Pregnancy
2.1.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"],
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"],
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)
nPrior = list(theta=list(prior = "normal", param=c(0, 10)))
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 = FALSE,
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 = FALSE, 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.9298 6.0889 0.3421 8.3608
##
## 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 mode
## Precision for COMRECPOP 0.8039 0.2481 0.4240 0.7683 1.3909 0.7020
## Precision for myID 1.9499 0.4843 1.1736 1.8911 3.0645 1.7793
## Precision for SexAgeGroup.f 0.5950 0.1518 0.3514 0.5768 0.9448 0.5422
## Precision for Month 0.5530 0.1376 0.3306 0.5369 0.8688 0.5064
## Precision for Int.DATE_YR 0.8303 0.2645 0.4299 0.7909 1.4592 0.7181
## Precision for SumDur 0.7754 0.2280 0.4223 0.7439 1.3112 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 = FALSE,
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 = FALSE, 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.7802 6.1020 0.3042 8.1865
##
## 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 mode
## Precision for COMRECDOEPOP 0.8267 0.2574 0.4352 0.7890 1.4370 0.7190
## Precision for myID 1.9287 0.4801 1.1597 1.8704 3.0342 1.7595
## Precision for SexAgeGroup.f 0.6089 0.1565 0.3578 0.5900 0.9697 0.5543
## Precision for Month 0.5547 0.1379 0.3316 0.5386 0.8712 0.5082
## Precision for Int.DATE_YR 0.8290 0.2615 0.4301 0.7911 1.4490 0.7205
## Precision for SumDur 0.8230 0.2449 0.4447 0.7889 1.3996 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 = FALSE,
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 = FALSE, 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.7753 5.6112 0.3410 7.7275
##
## 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 mode
## Precision for COMPRECNTDEN 1.0439 0.3324 0.5412 0.9945 1.8347 0.9031
## Precision for myID 1.8961 0.4749 1.1369 1.8381 2.9911 1.7279
## Precision for SexAgeGroup.f 0.6201 0.1600 0.3636 0.6008 0.9893 0.5641
## Precision for Month 0.5630 0.1402 0.3361 0.5467 0.8847 0.5158
## Precision for Int.DATE_YR 0.7291 0.2195 0.3914 0.6981 1.2465 0.6402
## Precision for SumDur 0.8600 0.2562 0.4648 0.8241 1.4630 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). One estimate for each density value.
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, col="darkgray") +
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),
plot.margin = unit(c(1,3,1,1),"cm"),
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))
Linear version of above
ggplot(mic.d, aes(ID, Meane, group = ID2)) +
geom_point(shape=1, col="darkgray") +
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="lm",
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),
plot.margin = unit(c(1,1,1,1),"cm"),
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))
Above, but as panels (Model estimate probability for each density value)
my.formula = y ~ x
ggplot(mic.d, aes(ID, Meane)) +
geom_point(shape=1, cex=3) +
geom_smooth(data= mic.d,
aes(ID, Meane,
#group = mic.d$ID2,
linetype= "solid"),
se=FALSE,
#span = 0.8,
method="lm",
col = "black",
formula = my.formula,
lwd = 1) +
stat_poly_eq(formula = my.formula,
aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")),
parse = TRUE, label.x = 35, label.y = 0.75, cex=5) +
geom_hline(yintercept = plogis(Pregnancy.mod$summary.fixed[1,1]),
linetype = "solid",
col = "darkgray",
size = 0.5) +
facet_wrap(~ID2, ncol = 1, scales = "free_y") +
ylab("Pregnancy (probability)") +
xlab(expression(Deer~km^2)) +
ggtitle("Pregnancy (COMRECPOP)") +
theme_classic() +
scale_y_continuous(breaks = scales::pretty_breaks(5)) +
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",
plot.margin = unit(c(1,3,1,1),"cm"),
strip.background = element_blank(),
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 = 18),
axis.title.y = element_text(face="bold", size = 23),
axis.title.x = element_text(face="bold", size = 22, vjust=-2))
Model Fitted Estimates.
Pullling estimates for individual deer.
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)))Bar graph version. Mean values across individual deer.
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),
plot.margin = unit(c(1,1,1,1),"cm"),
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))
Smooth line values.
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)))
hline_dat1## Class Med
## 1 Fawn (f) 0.07913857
## 2 Yearling (f) 0.61005621
## 3 Adult (f) 0.70878201ggplot(Temp.data, aes(Density, Pred.Preg)) +
geom_point(shape=1, col="darkgray") +
geom_smooth(method="loess",
se=FALSE,
span=1,
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 (COMRECPOP)") +
theme_classic() +
theme(axis.text=element_text(size=16),
plot.margin = unit(c(1,3,1,1),"cm"),
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))
Linear version of above (Model estimate probability for individual deer)
ggplot(Temp.data, aes(Density, Pred.Preg)) +
geom_point(shape=1, col="darkgray") +
geom_smooth(method="lm",
se=FALSE,
col="black",
lwd = 1) +
facet_wrap(~Class) +
ylab("Pregnancy Probability") +
xlab(expression(Deer~km^2)) +
ggtitle("Estimated Pregnancy (COMRECPOP)") +
theme_classic() +
theme(axis.text=element_text(size=16),
plot.margin = unit(c(1,1,1,1),"cm"),
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))
my.formula = y ~ x
ggplot(Temp.data, aes(Density, Pred.Preg)) +
geom_point(shape=1, cex=3, col ="darkgray") +
geom_smooth(linetype= "solid",
se=FALSE,
method="lm",
col = "black",
formula = my.formula,
lwd = 1) +
stat_poly_eq(formula = my.formula,
aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")),
parse = TRUE, label.x = 35, label.y = 0.25, cex=5) +
facet_wrap(~Class, ncol = 1, scales = "free_y") +
ylab("Pregnancy (probability)") +
xlab(expression(Deer~km^2)) +
ggtitle("Pregnancy (COMRECPOP)") +
theme_classic() +
#scale_y_continuous(breaks = scales::pretty_breaks(5)) +
#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",
plot.margin = unit(c(1,3,1,1),"cm"),
strip.background = element_blank(),
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 = 18),
axis.title.y = element_text(face="bold", size = 23),
axis.title.x = element_text(face="bold", size = 22, vjust=-2))