|
| 1 | +#' Function to project SRKW population in the future using individual based models |
| 2 | +#' |
| 3 | +#' @param whale_data the expanded data frame, returned by `expand()` |
| 4 | +#' @param seed the random seed used to initialize simulations, defaults to `123` |
| 5 | +#' @param n_years the length of the time series projections, defaults to 30 time steps |
| 6 | +#' @param n_iter the number of iterations per scenario, defaults to 200 |
| 7 | +#' @param scenarios a character string, or vector of strings describing which scenario |
| 8 | +#' to be run. Should be entered as a range of years separated by a colon, |
| 9 | +#' e.g. "1981:1985" to use rates from 1981 to 1985 or ("1981:1985","1981:2021") |
| 10 | +#' @param verbose whether to print status updates to screen, defaults to FALSE |
| 11 | +#' @param p_female The probability of a female birth, defaults to the empirical value since 1980 |
| 12 | +#' @param n_births Used for uncertainty in the sex ratio at birth. The distribution of female births is |
| 13 | +#' assumed to be beta binomial, with shape parameters `(n_births+1)*p_female` and `(n_births+1)*(1-p_female)`. |
| 14 | +#' There have been approximately 100 births since 1980, and this defaults to the empirical number |
| 15 | +#' of births |
| 16 | +#' formula above |
| 17 | +#' |
| 18 | +#' @import dplyr |
| 19 | +#' @import mgcv |
| 20 | +#' @importFrom stats rnorm runif |
| 21 | +#' @importFrom stringr str_locate |
| 22 | +#' |
| 23 | +#' @export |
| 24 | +project = function(whale_data, seed = 123, |
| 25 | + n_years = 30, |
| 26 | + scenarios = c("1981:2020"), |
| 27 | + n_iter = 200, |
| 28 | + verbose = FALSE, |
| 29 | + p_female = NA, |
| 30 | + sd_female = NA) { |
| 31 | + set.seed(seed) |
| 32 | + |
| 33 | + data = whale_data |
| 34 | + |
| 35 | + # First fit models with survival |
| 36 | + whale_data[which(whale_data$animal %in% c("L095","L098","L112") & whale_data$alive==0),] = NA |
| 37 | + # filter out animals born > 1960 |
| 38 | + whale_data = dplyr::filter(whale_data, birth > 1960) |
| 39 | + # fill in missing sexes randomly |
| 40 | + whale_data$sexF1M2[which(whale_data$sexF1M2==0)] = sample(1:2, size=length(which(whale_data$sexF1M2==0)), replace=T) |
| 41 | + whale_data = dplyr::rename(whale_data, sex = sexF1M2) |
| 42 | + surv_fit <- gam(alive ~ s(year) + age+I(age^2)+sex,data=dplyr::filter(whale_data,includeSurv==1), |
| 43 | + family="binomial") |
| 44 | + |
| 45 | + # Second fit models to fecundity |
| 46 | + whale_data = data |
| 47 | + whale_data = dplyr::filter(whale_data, sexF1M2==1, |
| 48 | + includeFec==1, |
| 49 | + age >= 10, |
| 50 | + age<=43, |
| 51 | + alive == 1, |
| 52 | + !is.na(gave_birth)) |
| 53 | + fec_fit <- gam(gave_birth ~ s(year) + age+I(age^2),data=whale_data, |
| 54 | + family="binomial") |
| 55 | + |
| 56 | + total_popsize = list() |
| 57 | + repro_females = list() |
| 58 | + |
| 59 | + # do the projections here |
| 60 | + for(s in 1:length(scenarios)) { |
| 61 | + indx = str_locate(scenarios[s],":") |
| 62 | + if(!is.na(indx[1])) { |
| 63 | + year_start = as.numeric(substr(scenarios[s],1,indx[1]-1)) |
| 64 | + year_end = as.numeric(substr(scenarios[s],indx[2]+1, nchar(scenarios[s]))) |
| 65 | + years_to_sample = seq(year_start, year_end) |
| 66 | + } |
| 67 | + |
| 68 | + # These were the sex ratios we'd used during the workshops -- this is |
| 69 | + # a random variable, rather than fixed, so the uncertainty is propagated |
| 70 | + srb = dplyr::filter(data,birth>1980) %>% |
| 71 | + dplyr::group_by(animal) %>% |
| 72 | + dplyr::summarize(sex=sexF1M2[1]) %>% |
| 73 | + dplyr::filter(sex!=0) |
| 74 | + if(is.na(p_female)) { |
| 75 | + p_female = table(srb$sex)[1]/nrow(srb) |
| 76 | + } |
| 77 | + if(is.na(n_births)) { |
| 78 | + n_births = nrow(srb) |
| 79 | + } |
| 80 | + beta_1 = p_female*(n_births+1) |
| 81 | + beta_2 = (1-p_female)*(n_births+1) |
| 82 | + |
| 83 | + popSize = array(0, dim=c(n_iter, n_years)) |
| 84 | + popFems = array(0, dim=c(n_iter, n_years)) |
| 85 | + currentPop = dplyr::filter(data, |
| 86 | + year == max(whale_data$year), |
| 87 | + alive==1) |
| 88 | + for(i in 1:n_iter) { |
| 89 | + # get the current age / sex / pod structure of the population |
| 90 | + # to make this more general, I called "populationSegment" = pod, and "breedingGroup" = matriline |
| 91 | + initPopCurrent = currentPop[,c("animal", "pod", "matriline", "sexF1M2", "age")] # animal, pod, matriline, sex, age |
| 92 | + initPopCurrent = dplyr::rename(initPopCurrent, sex=sexF1M2) |
| 93 | + initPopCurrent$animal = as.character(initPopCurrent$animal) |
| 94 | + initPopCurrent$matriline = as.numeric(as.factor(initPopCurrent$matriline)) # , -29 because of NRs |
| 95 | + initPopCurrent$dad = NA # keep track of dads |
| 96 | + |
| 97 | + numWhale = dim(initPopCurrent)[1] |
| 98 | + initPopCurrent$hasCalf = 0 |
| 99 | + initPopCurrent$hasCalf[which(as.character(currentPop$animal) %in% as.character(currentPop$mom[which(currentPop$age==1)]))] = 1 |
| 100 | + initPopCurrent$mom = as.character(currentPop$mom) |
| 101 | + |
| 102 | + initPopCurrent$sex = as.numeric(initPopCurrent$sex) |
| 103 | + # assign sexes to unsexed individuals [0s] |
| 104 | + initPopCurrent$sex[which(initPopCurrent$Sex==0)] = ifelse(runif(length(which(initPopCurrent$sex==0))) < rnorm(length(which(initPopCurrent$sex==0)),p_female, sd_female), 1, 2) |
| 105 | + newID = 9999 |
| 106 | + |
| 107 | + for(yrs in 1:dim(popSize)[2]) { |
| 108 | + # first step is find females available to give birth |
| 109 | + indx = which(initPopCurrent$sex == 1 & as.numeric(initPopCurrent$age) >= 10 & as.numeric(initPopCurrent$age) < 43 & initPopCurrent$hasCalf == 0) |
| 110 | + # second step is to calculate predicted fecundity rates and make fecundity stochastic - every individual's pregnancy is a coin flip |
| 111 | + |
| 112 | + # if a range of years is used to draw demographic rates from sample those randomly |
| 113 | + initPopCurrent$year = as.numeric(sample(c(as.character(years_to_sample)), size=1)) |
| 114 | + |
| 115 | + if(length(indx) > 0) { |
| 116 | + # bind together the current moms with the matching fecundity data |
| 117 | + p_fec = mgcv::predict.gam(fec_fit, initPopCurrent[indx,], type="response") |
| 118 | + |
| 119 | + pregMoms = indx[which(runif(length(p_fec)) < p_fec)] |
| 120 | + # third step is to make moms that aren't mate limited give birth to calves of known sexes |
| 121 | + if(length(pregMoms) > 0) { |
| 122 | + |
| 123 | + for(ll in 1:length(pregMoms)) { |
| 124 | + # loop over moms and only let them breed if there's a mature male in a different matriline |
| 125 | + dads = which(initPopCurrent$sex == 2 & as.numeric(initPopCurrent$age) > 12) |
| 126 | + #dads = which(initPopCurrent$Sex == 2 & initPopCurrent$breedingGroup != initPopCurrent$breedingGroup[pregMoms[ll]] & as.numeric(initPopCurrent$age1) > 12) |
| 127 | + if(length(dads) > 0) { |
| 128 | + # assign the pod / matriline to be the same of the mom |
| 129 | + newpod = initPopCurrent$pod[pregMoms[ll]] |
| 130 | + newmat = initPopCurrent$matriline[pregMoms[ll]] |
| 131 | + # sex is stochastic |
| 132 | + newsex = ifelse(runif(1) < rbeta(1, beta_1, beta_2), 1, 2) |
| 133 | + # bookkeeping |
| 134 | + newage = 0 |
| 135 | + newcalf = 0 |
| 136 | + newmom = initPopCurrent$animal[pregMoms[ll]] |
| 137 | + |
| 138 | + # sample from potential dads in proprtion to their estimated relative reproductive output |
| 139 | + # their ages are initPopCurrent$age1[dads], and the fecundity model defined above outside of loops |
| 140 | + # sample from potential dads in proprtion to their estimated relative reproductive output |
| 141 | + # their ages are initPopCurrent$age1[dads], and the fecundity model defined above outside of loops |
| 142 | + |
| 143 | + #inflation_factor = 5 |
| 144 | + #pred.male.fecundity[32:200] = pred.male.fecundity[31] |
| 145 | + #probs = pred.male.fecundity[as.numeric(initPopCurrent$age[dads])] |
| 146 | + #probs[which.max(initPopCurrent$age1[dads])] = probs[which.max(initPopCurrent$age[dads])] * inflation_factor |
| 147 | + #newdad = initPopCurrent$animal[sample(dads,1, prob=probs)] |
| 148 | + newdad = initPopCurrent$animal[sample(dads,1)] |
| 149 | + newID = newID + 1 |
| 150 | + # add calves to the population |
| 151 | + newdf = data.frame(animal=newID, |
| 152 | + pod=newpod,matriline=newmat,sex=newsex, |
| 153 | + age=newage,dad=newdad,hasCalf=newcalf, |
| 154 | + mom=newmom,year=initPopCurrent$year[1]) |
| 155 | + initPopCurrent = rbind(initPopCurrent, newdf) |
| 156 | + |
| 157 | + }# end if(length(dads) > 0) { |
| 158 | + } # end ll loop |
| 159 | + } # end if(length(pregMoms) |
| 160 | + |
| 161 | + } |
| 162 | + # bookkeeping: update whales that have calves |
| 163 | + initPopCurrent$hasCalf[which(initPopCurrent$hasCalf==1)] = 0 |
| 164 | + initPopCurrent$hasCalf[pregMoms] = 1 |
| 165 | + |
| 166 | + # step 4 is calculate predicted survival at age |
| 167 | + initPopCurrent$sexF1M2 = as.numeric(initPopCurrent$sex) |
| 168 | + p_surv = mgcv::predict.gam(surv_fit, initPopCurrent, type="response") |
| 169 | + |
| 170 | + initPopCurrent = dplyr::select(initPopCurrent, -sexF1M2) |
| 171 | + |
| 172 | + # step 5: stochastic survival to kill whales |
| 173 | + liveOrDie = rep(1,length(p_surv)) |
| 174 | + dead = which(runif(length(p_surv)) > p_surv) |
| 175 | + liveOrDie[dead] = 0 |
| 176 | + |
| 177 | + # step 6: see if any of these dead animals has a calf - if so, kill the calf too |
| 178 | + if(length(which(liveOrDie == 0)) > 0) { |
| 179 | + for(ll in 1:length(dead)) { |
| 180 | + # kill the calf |
| 181 | + if(is.na(initPopCurrent$hasCalf[dead[ll]]) == FALSE & initPopCurrent$hasCalf[dead[ll]] == 1) { |
| 182 | + liveOrDie[which(initPopCurrent$mom == dead[ll])] = 0 |
| 183 | + } |
| 184 | + } |
| 185 | + } |
| 186 | + |
| 187 | + # step 7: bookkeeping at the end of the time step |
| 188 | + # first remove dead animals from the population |
| 189 | + if(length(dead) > 0 ) deadWhales = initPopCurrent[which(liveOrDie==0),] ## MIKE ADDITION - list of who is dead |
| 190 | + if(length(dead) > 0) initPopCurrent = initPopCurrent[-which(liveOrDie==0),] |
| 191 | + # second age remaining whales |
| 192 | + initPopCurrent$age = as.numeric(initPopCurrent$age) + 1 |
| 193 | + # third record pop size to later calculate recovery targets |
| 194 | + popSize[i,yrs] = dim(initPopCurrent)[1] |
| 195 | + popFems[i,yrs] = dim(dplyr::filter(initPopCurrent, sex==1, |
| 196 | + age <= 42, age >= 10))[1] |
| 197 | + if(verbose==TRUE) print(paste0("Iteration ",i, " : year ", yrs)) |
| 198 | + } # this is end of yrs loop |
| 199 | + |
| 200 | + } # end i loop |
| 201 | + |
| 202 | + total_popsize[[s]] = popSize |
| 203 | + repro_females[[s]] = popFems |
| 204 | + } |
| 205 | + |
| 206 | + return(list(fec_fit = fec_fit, surv_fit = surv_fit, |
| 207 | + total_popsize = total_popsize, |
| 208 | + repro_females = repro_females)) |
| 209 | +} |
| 210 | + |
0 commit comments