LandSciTech / caribouMetrics

Models and metrics of boreal caribou responses to forest landscapes
https://landscitech.github.io/caribouMetrics
Other
3 stars 3 forks source link

Simulating probability of sustaining stable growth #109

Closed northernbio closed 2 months ago

northernbio commented 2 months ago

About a month ago I was testing the caribouMetrics for estimating recruitment and survivorship, relative to disturbance, and it is working well for me. I've also tried simulating population growth to replicate Figure 11. "Probability of observing stable or positive growth (λ ≥ stable)", in the ECCC 2011 update report (a sigmoidal graph). Given a set level of disturbance, I would simulate population growth 1000 times, then count the number of time lambda was the equal to or greater than 1. I expected, based on Fig. 11, that 35% disturbance should result in about 65% of the runs having lambda >= 1, but instead found only a tiny number of runs to have lambda > 1. I have never found a link to the actual model on the web or the appendix it is documented in, and it is referenced as an ensemble model in the report, so perhaps that means it takes the average of multiple models. In any case, I'm guessing that the variance in your caribouPopGrowth function results is much tighter than the ECCC 2011 model. My question: Have I interpreted what is going on correctly, and do you have any suggestions in how I can use caribouMetrics to simulate the model used in ECCC Fig. 11? Also, if you have documentation on that model, would you mind sending it to me? Finally, I see that new updates have been posted on the caribouPopGrowth code, but I have not tested the new code.

josie-hughes commented 2 months ago

The disturbance management threshold was informed by analysis and modelling that differs in a number of ways from the Johnson et al 2020 model that is included in caribouMetrics. However, the most important differences are the survival and recruitment regression model parameters. Essentially, the more recent Johnson et al. 2020 analysis of more recent data from more herds gives somewhat different demographic disturbance relationships (i.e. compare Fig 8 of the ECCC 2011 science assessment to Fig 3 of Johnson et al 2020). This code shows how to use older parameters (from table 56 of ECCC 2011 science assessment) with the newer model to get an answer similar to the ECCC 2011 result.

library(caribouMetrics)
library(tidyverse)

#Use ECCC 2011 model parameters - from table 56
coefTab2011 = read.table(text="modelVersion responseVariable    ModelNumber Type    Coefficient Value   StdErr  lowerCI upperCI
ECCC    femaleSurvival  M2  National    Intercept   0.87    0   NA  NA
ECCC    recruitment M3  National    Intercept   44.265  2.942   NA  NA
ECCC    recruitment M3  National    Total_dist  -0.429  0.061   NA  NA",header=T)

covTableSim <-data.frame(Total_dist=seq(0,100,by=5))
popGrowthPars <- demographicCoefficients(
  1000,
  modelVersion = "ECCC",
  survivalModelNumber = "M2",
  recruitmentModelNumber = "M3",
  populationGrowthTable = coefTab2011,
  useQuantiles=F
)

rateSamples <- demographicRates(
  covTable = covTableSim,
  popGrowthPars = popGrowthPars,
  ignorePrecision = T,
  returnSample = T,
  useQuantiles = F
)

pars <- cbind(rateSamples,
              caribouPopGrowth(rep(300,nrow(rateSamples)),
                               R_bar = rateSamples$R_bar, S_bar = rateSamples$S_bar,
                               numSteps = 20))

propViable = pars %>% group_by(Total_dist) %>%
  summarise(propViable = sum((lambda>0.99)/length(lambda)))

base = ggplot(propViable,aes(x=Total_dist,y=propViable*100))+geom_line()+
  scale_x_continuous(breaks = seq(0, 100, by = 10)) +
  scale_y_continuous(breaks = seq(0, 100, by = 10)) +
  geom_hline(yintercept=60)+geom_vline(xintercept=35)+
  ylab("probability of lambda > 0.99")+xlab("% total disturbance")
print(base)