reshape_array.Rdpredictnarwsim <- function(obj, yrs = 35, n = 100) # if(is.null(obj$gam)) stop("Insufficient data available. Cannot proceed with population projections.") # if(!identical(cohortID, 1:6)) stop("Missing cohorts in input <obj>. Cannot proceed with population projections.") # if(length(obj$gam$fit$surv$xlevels[[1]]) < 6 | length(obj$gam$fit$bc$xlevels[[1]]) < 6) stop("Missing factor levels in input <obj>. Cannot proceed with population projections.") # plogis("link" predictions + error) # Adapted from original code by Scott Creel # https://www.montana.edu/screel/teaching/bioe-440r-521/course-outline/stoch_projection_new.R # # Prediction intervals # https://stat.ethz.ch/pipermail/r-help/2011-April/275632.html # Population estimate as per 2022 NARW report card is 340 (+/- 7). # test <- mgcv::predict.gam(mod$gam[[1]]$surv, newdata = data.frame(start_bc = seq(0,find_maxBC(), 0.01)), type = "link", se.fit = TRUE) # plot(seq(0,find_maxBC(), 0.01), plogis(test$fit), type = "l") # lines(seq(0,find_maxBC(), 0.01), plogis(Reduce("+", test)), lty = 2) # lines(seq(0,find_maxBC(), 0.01), plogis(Reduce("-", test)), lty = 2) # test2 <- mgcv::predict.gam(mod$gam[[1]]$surv, newdata = data.frame(start_bc = 0.2), type = "response") # abline(v = 0.2) # abline(h = tesct2) cohortID <- obj$param$cohortID cohorts <- obj$param$cohorts |> dplyr::slice(1) |> dplyr::mutate(name = gsub(", female", "", name), abb = gsub(",f", "", abb)) |> dplyr::bind_rows(obj$param$cohorts) |> dplyr::mutate(name = gsub("male, female", "female", name), abb = gsub("m,f", "f", abb)) cohorts <- cohorts[c(1,3,5,2,4,6,7,8)] # Attributes to monitor during projection mat.attribs <- c("alive", "cohort", "female", "age", "length", "length_a", "length_b", "length_c", "tot_mass", "lean_mass", "bc", "mass_a", "mass_b", "p_surv", "min_bc", "calving_int", "t2calf", "birth") # Current year current.yr <- lubridate::year(lubridate::now()) # Extract terminal functions mod <- obj$gam$fit mod[["gest"]] <- gam_gest #'------------------------------------------------------ # INITIALIZATION #'------------------------------------------------------ # Define initial population vector N0 <- c(2, 5, 212, 2, 69, 1, 7, 61) names(N0) <- cohorts[, name] # Create matrices and initialize them # rows: years <yrs> # columns: <attributes> # layers: individuals <n> # 4th dimension: replicate projection -> later converted to list narw.indiv <- array(data = NA, c(yrs + 1, length(mat.attribs), sum(N0), n), dimnames = list(paste0("yr ", 0:yrs), mat.attribs, paste0("whale ", 1:sum(N0)), paste0("prj ", 1:n))) cohort.vec <- do.call(c, sapply(X = 1:nrow(cohorts), FUN = function(x) rep(cohorts$id[x], each = N0[x]))) # Alive and population cohort narw.indiv[1,"alive",,] <- 1 narw.indiv[1,"cohort",,] <- rep(cohort.vec, n) # Sex # -- Calves (male) narw.indiv[,"female", 1:N0[1], ] <- 0 # -- Calves (female) fem <- which(cohort.vec == 0) narw.indiv[,"female", fem[fem > N0[1]], ] <- 1 # -- Juveniles and adults narw.indiv[,"female", which(cohort.vec narw.indiv[,"female", which(cohort.vec # Age ages <- start_age_vec(rep(cohort.vec, n)) narw.indiv[1,"age", , ] <- ages # Total body length l.params <- agL_vec(ages) lengths <- age2length_vec(ages, l.params) narw.indiv[1,"length", , ] <- lengths narw.indiv[,"length_a",,] <- rep(l.params[,1], each = yrs + 1) narw.indiv[,"length_b",,] <- rep(l.params[,2], each = yrs + 1) narw.indiv[,"length_c",,] <- rep(l.params[,3], each = yrs + 1) # Total mass m.params <- mL(n*sum(N0)) mass <- length2mass_vec(lengths, m.params, FALSE) narw.indiv[,"mass_a",,] <- rep(m.params[,1], each = yrs + 1) narw.indiv[,"mass_b",,] <- rep(m.params[,2], each = yrs + 1) narw.indiv[1,"tot_mass", , ] <- mass # Body conditon narw.indiv[1,"bc",,] <- start_bcondition_vec(rep(cohort.vec, n)) # lean mass narw.indiv[1,"lean_mass", , ] <- mass - (narw.indiv[1,"bc",,] * mass) # Calving interval - mean of 5.3 years, up to apparent gaps of 13 years (Kraus et al. 2001) # NARWC Report Card 2022: Average inter-birth interval of 7.7 yrs, median of 8, min/max (2/13) # This corresponds to a Normal (7.7, 1.45) # Mean age at first calving = 9.53 +/- 2.32 (Kraus et al. 2001) # # Stewart et al. 2022 -- # The degree to which the energetic reserves of females are depleted during lactation may govern # the length of the resting period between successful pregnancies (Miller et al. 2011, Marón et al. 2015). narw.indiv[1,"calving_int",,] <- (rep(cohort.vec, n) == 6) * rnorm(sum(N0)*n, mean = 7.7, sd = 1.45) narw.indiv[1,"t2calf",,] <- round(narw.indiv[1,"calving_int",,],0) narw.indiv[1,"min_bc",,] <- predict_m(model = mod, values = as.vector(narw.indiv[1,"tot_mass",,]), prediction = "gest") * as.vector(narw.indiv[1, "cohort",, ] == 6) narw.indiv[1,"birth",,] <- 0 narw.indiv[1,"p_surv",,] <- 1 #' --------------------------- # IMPORTANT #' --------------------------- # Turn array into a list - otherwise cannot add new calves within projections narw.indiv <- purrr::array_branch(narw.indiv, 4) # Number of individuals in each cohort narw.pop <- array(data = NA, c(n, yrs + 1, length(N0)), dimnames = list(paste0("prj ", 1:n), paste0("yr ", 0:yrs), cohorts$name)) narw.pop[,1,] <- rep(N0, each = n) # Total population size tot.pop <- matrix(0, n, yrs + 1, dimnames = list(paste0("prj",1:n), paste0("yr ", 0:yrs))) tot.pop[,1] <- sum(N0) #'------------------------------------------------------ # RUN PROJECTIONS #'------------------------------------------------------ # This uses nested loops. # The prj loop (outermost loop) replicates the projection <n> times. # The i loop is next, and steps across all years of projection from an initial population vector. # Set up progress bar pb <- progress::progress_bar$new( format = "[:bar] :percent eta: :eta", total = n, clear = FALSE, width = 80 ) for(prj in 1:n) pb$tick() # Update progress bar #'------------------------------------------------------ # Loop over years #'------------------------------------------------------ for(i in 2:(yrs+1)) current.dat <- as.matrix(narw.indiv[[prj]][i-1, , ]) alive <- narw.indiv[[prj]][i-1, "alive", ] #' ---------------------------- # SURVIVAL #' ---------------------------- # Predict survival probability based on body condition ps <- alive * predict_m(model = mod, cohort = current.dat["cohort",], values = current.dat["bc",], prediction = "surv") narw.indiv[[prj]][i, "p_surv", ] <- ps # Determine whether the animal survived # Maximum longevity of 69 years alive <- rbinom(n = dim(narw.indiv[[prj]])[3], size = 1, prob = ps) * (narw.indiv[[prj]][i-1, "age", ] <=69) narw.indiv[[prj]][i, "alive", ] <- alive # Sex remains the same narw.indiv[[prj]][i, "female", ] <- narw.indiv[[prj]][i-1, "female", ] #' ---------------------------- # GROWTH #' ---------------------------- # Increment age narw.indiv[[prj]][i, "age", ] <- alive * narw.indiv[[prj]][i-1, "age", ] + 1 # Increment length narw.indiv[[prj]][i,"length", ] <- alive * age2length_vec(narw.indiv[[prj]][i, "age", ], t(narw.indiv[[prj]][i, c("length_a", "length_b", "length_c"),])) # Increment lean mass narw.indiv[[prj]][i,"lean_mass", ] <- alive * length2mass_vec(narw.indiv[[prj]][i, "length", ], t(narw.indiv[[prj]][i, c("mass_a", "mass_b"),]), lean = TRUE) # Predict new body condition from current body condition narw.indiv[[prj]][i, "bc", ] <- alive * predict_m(model = mod, cohort = current.dat["cohort",], values = current.dat["bc",], prediction = "bc") # Increment total mass narw.indiv[[prj]][i,"tot_mass", ] <- alive * narw.indiv[[prj]][i,"lean_mass", ] / (1 - narw.indiv[[prj]][i,"bc", ]) #' ---------------------------- # REPRODUCTION #' ---------------------------- narw.indiv[[prj]][i, "calving_int", ] <- alive * narw.indiv[[prj]][i-1, "calving_int", ] narw.indiv[[prj]][i, "t2calf", ] <- ifelse(narw.indiv[[prj]][i-1, "t2calf", ] - 1 < 0, 0, narw.indiv[[prj]][i-1, "t2calf", ] - 1) # Minimum body condition needed to successfully bring fetus to term without starving # No evidence of reproductive senescence in right whales - Hamilton et al. (1998) narw.indiv[[prj]][i, "min_bc", ] <- alive * predict_m(model = mod, values = narw.indiv[[prj]][i, "tot_mass", ], prediction = "gest") * (narw.indiv[[prj]][i-1, "cohort", ] == 6) # Birth of new calf narw.indiv[[prj]][i, "birth", ] <- alive * (narw.indiv[[prj]][i, "bc", ] >= narw.indiv[[prj]][i, "min_bc", ]) * (narw.indiv[[prj]][i-1, "cohort", ] == 6) * (narw.indiv[[prj]][i-1, "birth", ] == 0) # Maturity - transitions between cohorts narw.indiv[[prj]][i, "cohort", ] <- alive * increment_cohort(cohort = narw.indiv[[prj]][i-1, "cohort", ], age = narw.indiv[[prj]][i, "age", ], female = narw.indiv[[prj]][i, "female", ], t2calf = narw.indiv[[prj]][i, "t2calf", ]) # New births # narw.indiv[[prj]] <- new_calf(narw.indiv[[prj]][,,], year = i, attrib = mat.attribs) # Number of animals in each cohort # Calves (male) narw.pop[prj, i, 1] <- sum((narw.indiv[[prj]][i, "cohort", ] == 0) * (narw.indiv[[prj]][i, "female", ] == 0) * (narw.indiv[[prj]][i, "alive", ] == 1)) # Juveniles and adults (male) narw.pop[prj, i, 2] <- sum((narw.indiv[[prj]][i, "cohort", ] == 1) * (narw.indiv[[prj]][i, "alive", ] == 1)) narw.pop[prj, i, 3] <- sum((narw.indiv[[prj]][i, "cohort", ] == 3) * (narw.indiv[[prj]][i, "alive", ] == 1)) # Calves (female) narw.pop[prj, i, 4] <- sum((narw.indiv[[prj]][i, "cohort", ] == 0) * (narw.indiv[[prj]][i, "female", ] == 1) * (narw.indiv[[prj]][i, "alive", ] == 1)) # Juvenile and reproductive adults (female) narw.pop[prj, i, 5] <- sum((narw.indiv[[prj]][i, "cohort", ] == 2) * (narw.indiv[[prj]][i, "alive", ] == 1)) narw.pop[prj, i, 6] <- sum((narw.indiv[[prj]][i, "cohort", ] == 4) * (narw.indiv[[prj]][i, "alive", ] == 1)) narw.pop[prj, i, 7] <- sum((narw.indiv[[prj]][i, "cohort", ] == 5) * (narw.indiv[[prj]][i, "alive", ] == 1)) narw.pop[prj, i, 8] <- sum((narw.indiv[[prj]][i, "cohort", ] == 6) * (narw.indiv[[prj]][i, "alive", ] == 1)) # Total population size tot.pop[prj, i] <- sum(narw.indiv[[prj]][i, "alive", ], na.rm = TRUE) # End years # End projections # Compile outputs narw.df <- purrr::map(.x = cohorts$name, .f = ~ narw.pop[,,.x] |> tibble::as_tibble() |> tibble::rownames_to_column(var = "prj") |> tidyr::pivot_longer(!prj, names_to = "year", values_to = "N") |> dplyr::mutate(year = current.yr + as.numeric(gsub("yr ", "", year))) |> dplyr::mutate(cohort = stringr::str_to_sentence(.x)) ) |> do.call(what = rbind) |> data.table::data.table() tot.df <- tibble::as_tibble(tot.pop) |> tibble::rownames_to_column(var = "prj") |> tidyr::pivot_longer(!prj, names_to = "year", values_to = "N") |> dplyr::mutate(year = current.yr + as.numeric(gsub("yr ", "", year))) |> dplyr::mutate(cohort = "North Atlantic right whales") |> data.table::data.table() narw.conf <- narw.df[, list(mean = mean(N), lwr = quantile(N, 0.025), uppr = quantile(N, 0.975)), list(year, cohort)] |> dplyr::mutate(cohort = factor(cohort, levels = c("Calves (male)", "Calves (female)", "Juveniles (male)", "Juveniles (female)", "Adults (male)", "Adults (female, pregnant)", "Adults (female, lactating)", "Adults (female, resting)"))) tot.conf <- tot.df[, list(mean = mean(N), lwr = quantile(N, 0.025), uppr = quantile(N, 0.975)), list(year, cohort)]|> dplyr::mutate(cohort = factor(cohort, levels = c("Calves (male)", "Calves (female)", "Juveniles (male)", "Juveniles (female)", "Adults (male)", "Adults (female, pregnant)", "Adults (female, lactating)", "Adults (female, resting)", "North Atlantic right whales"))) p1 <- plot_projection(narw.df, narw.conf) p2 <- plot_projection(tot.df, tot.conf) print(p1) print(p2) # Find 95 cat("Final population size:") final.pop <- unname(tot.df[year == current.yr + yrs, quantile(N, probs = c(0.5, 0.025, 0.975))]) cat("N = ", final.pop[1], " (95 return(list(ind = narw.indiv, df = rbind(narw.df, tot.df), ci = rbind(narw.conf, tot.conf))) predict_leslie <- function(obj, n = 100, yrs = 35) gg.opts <- ggplot2::theme(axis.text = element_text(size = 10, color = "black"), axis.text.y = element_text(margin = margin(t = 0, r = 5, b = 0, l = 0)), axis.text.x = element_text(margin = margin(t = 10, r = 0, b = 0, l = 0)), axis.title = element_text(size = 12), axis.title.y = element_text(margin = margin(t = 0, r = 20, b = 0, l = 0)), axis.title.x = element_text(margin = margin(t = 20, r = 0, b = 0, l = 0)), strip.background = element_rect(fill = "grey20"), strip.text = element_text(colour = 'white', size = 12)) # Adapted from original code by Scott Creel # https://www.montana.edu/screel/teaching/bioe-440r-521/course-outline/stoch_projection_new.R # Population estimate as per 2022 NARW report card is 340 (+/- 7). # The percentage of the population estimated # cohort.id <- obj$param$cohort.id # cohort.ab <- obj$param$cohort.ab current.yr <- lubridate::year(lubridate::now()) cohorts <- c("Calves (male)", "Juveniles (male)", "Adults (male)", "Calves (female)", "Juveniles (female)", "Adults (female, pregnant)", "Adults (female, lactating)", "Adults (female, resting)") |> tolower() |> abbreviate(minlength = 6) |> tibble::enframe() |> dplyr::rename(cohort = name, abbr = value) # Vector to store total population size for N replicate projections narwpop <- purrr::set_names(x = cohorts$cohort) |> purrr::map(.f = ~matrix(0, n, yrs, dimnames = list(paste0("proj",1:n), paste0("yr ", 1:yrs)))) totpop <- matrix(0, n, yrs, dimnames = list(paste0("proj",1:n), paste0("yr ", 1:yrs))) # 1: Calves (male) # 2: Juveniles (male) # 3: Adults (male) # 4: Calves (female) # 5: Juveniles (female) # 6: Adults (female, pregnant) # 7: Adults (female, lactating) # 8: Adults (female, resting) p.fem <- 0.5 # Fraction of females at birth p.birth <- 0.5 # Probability that a pregnant female will give birth to a calf (fecundity) p.survival <- rep(0.8, nrow(cohorts)) # Probability of survival names(p.survival) <- cohorts p.recruit <- c(0.7, 0.7) # Male, female tau_rp <- 0.7 # Transition prob from resting to pregnant tau_pl <- 0.7 # Transition prob from pregnant to lactating tau_lr <- 0.8 # Transition prob from lactating to resting popmat <- matrix(0, nrow(cohorts), nrow(cohorts)) popmat[1, 6] <- p.survival[6] * p.birth * (1 - p.fem) popmat[2, 1] <- p.survival[1] popmat[2, 2] <- p.survival[2] popmat[3, 2] <- p.survival[2] * p.recruit[1] popmat[3, 3] <- p.survival[3] popmat[4, 4] <- p.survival[6] * p.birth * p.fem popmat[5, 5] <- p.survival[4] popmat[5, 6] <- p.survival[5] * (1 - p.recruit[2]) popmat[6, 5] <- p.survival[5] * p.recruit[2] popmat[6, 6] <- p.survival[6] * (1 - tau_pl) popmat[6, 8] <- p.survival[8] * tau_rp popmat[7, 6] <- p.survival[6] * tau_pl popmat[7, 7] <- p.survival[7] * (1 - tau_lr) popmat[8, 7] <- p.survival[7] * tau_lr popmat[8, 8] <- p.survival[8] * (1 - tau_rp) popmat <- matrix(popmat, nrow = nrow(cohorts), byrow = FALSE, dimnames = list(cohorts$abbr, cohorts$abbr)) popmat <- round(popmat, 5) # STOCHASTIC LESLIE PROJECTION #' ---------------------------- # This uses 4 nested loops. # The p loop (outermost loop) replicates the projection <n> times. # The y loop is next, and steps across all years of projection from an initial population vector. # The i and j loops are innermost and draw stochastic parameters (body condition and survival) # Set up progress bar pb <- progress::progress_bar$new( format = "[:bar] :percent eta: :eta", total = n, clear = FALSE, width = 80 ) for(p in 1:n) pb$tick() # Update progress bar # Matrix of age-class values, with one row per time step. # The columns represent the numbers of individuals in each age class. N.mat <- matrix(NA, nrow = yrs, ncol = nrow(cohorts), dimnames = list(paste("yr", 1:yrs), cohorts$abbr)) # Initial population vector # Sums to 340 (abundance estimate as of 2022). N.mat[1,] <- c(7, 10, 163, 8, 10, 50, 15, 77) # Begin loop for projections for(i in 2:yrs) # Set up a temporary Leslie matrix for each iteration # with the correct mean fecundities and survival rates leslie.mat <- popmat # Randomly draw fecundities for each class # Mean of Poisson is fecundity value # for(j in 1:n.stages) # leslie.mat[1,j] <- ifelse(popmat[1,j] == 0, 0, popmat[1,j] + rnorm(1, 0, 0.01)) # for(j in 1:nrow(cohorts)) for(k in 1:nrow(cohorts)) leslie.mat[k,j] <- ifelse(popmat[k,j] == 0, 0, popmat[k,j] + rnorm(1, 0, 0.01)) # Randomly draw survival probabilities for each class. # n.ind is number of individuals to calculate live/die for # for(k in 1:(n.stages-1)) # n.ind <- N.mat[(i-1),k] # # Need ifelse statement to deal with the possibility that there are no individuals in that age class # leslie.mat[(k+1),k] <- ifelse(n.ind > 1, rbinom(1,size = round(n.ind), p = les.mat[(k+1),k])/n.ind,0) # # Matrix multiplication for next time-step. N.mat[i,] <- leslie.mat # End i loop over time for(k in seq_along(narwpop)) narwpop[[k]][p,] <- N.mat[,k] totpop[p,] <- rowSums(N.mat) # Ends p loop over replicate projections # Compile outputs narw.df <- purrr::imap(.x = narwpop, .f = ~ tibble::as_tibble(.x) |> tibble::rownames_to_column(var = "proj") |> tidyr::pivot_longer(!proj, names_to = "year", values_to = "N") |> dplyr::mutate(year = current.yr + as.numeric(gsub("yr ", "", year))) |> dplyr::mutate(cohort = stringr::str_to_sentence(.y)) ) |> do.call(what = rbind) |> data.table::data.table() tot.df <- tibble::as_tibble(totpop) |> tibble::rownames_to_column(var = "proj") |> tidyr::pivot_longer(!proj, names_to = "year", values_to = "N") |> dplyr::mutate(year = current.yr + as.numeric(gsub("yr ", "", year))) |> dplyr::mutate(cohort = "North Atlantic right whales") |> data.table::data.table() narw.conf <- narw.df[, list(mean = mean(N), lwr = quantile(N, 0.025), uppr = quantile(N, 0.975)), list(year, cohort)] tot.conf <- tot.df[, list(mean = mean(N), lwr = quantile(N, 0.025), uppr = quantile(N, 0.975)), list(year, cohort)] # Generate plots make_plot <- function(df, conf) ggplot2::ggplot() + # ggplot2::geom_path(data = df, aes(x = year, y = N, group = factor(proj)), colour = "lightgrey", linewidth = 0.2) + ggplot2::geom_path(data = conf, aes(x = year, y = mean), colour = "#1565C0") + ggplot2::geom_ribbon(data = conf, aes(x = year, y = mean, ymin = lwr, ymax = uppr), alpha = 0.25, fill = "#1565C0") + ggplot2::facet_wrap(~cohort) + ggplot2::scale_x_continuous(breaks = pretty(seq(current.yr, current.yr + yrs))) + ggplot2::scale_y_continuous(breaks = pretty(df$N), labels = scales::comma) + xlab("") + ylab("Abundance") + gg.opts p1 <- make_plot(tot.df, tot.conf) p2 <- make_plot(narw.df, narw.conf) print(p1) print(p2) # Find 95 ci <- tot.df[year == current.yr + yrs, quantile(N, probs = c(0.025, 0.975))] return(ci)
reshape_array(a, value, ..., .id = NULL)