## ----setup, message = FALSE--------------------------------------------------- knitr::opts_chunk$set(echo = TRUE, tidy = TRUE) library(dplyr) library(purrr) library(tibble) library(ggplot2) library(popbio) # for projection.matrix() library(raretrans) # Raymond's theme modifications rlt_theme <- theme(axis.title.y = element_text(colour="grey20",size=15,face="bold"), axis.text.x = element_text(colour="grey20",size=15, face="bold"), axis.text.y = element_text(colour="grey20",size=15,face="bold"), axis.title.x = element_text(colour="grey20",size=15,face="bold")) ## ----viewData----------------------------------------------------------------- data("L_elto") # load the dataset `L_elto` into memory (from package `raretrans`) head(L_elto) ## ----mungeData---------------------------------------------------------------- onepop <- L_elto %>% # Filter out population # 250, period (year) 5 filter(POPNUM == 250, year == 5) %>% # redefine p for el plantón to s for seedling mutate(stage = case_when(stage == "p" ~ "s", TRUE ~ stage), next_stage = case_when(next_stage == "p"~ "s", TRUE ~ next_stage)) # popbio::projection.matrix doesn't like tibbles # set TF = TRUE to get the right format. TF <- popbio::projection.matrix(as.data.frame(onepop), stage = stage, fate = next_stage, fertility="fertility", sort=c("s","j","a"), TF = TRUE) TF # This is the stage-structured population model ## ----------------------------------------------------------------------------- N <- get_state_vector(onepop, stage = stage, sort=c("s","j","a")) N # A vector of # of starting individuals for each stage ## ----echo=FALSE--------------------------------------------------------------- Tmat <- L_elto %>% mutate(stage = factor(stage, levels=c("p","j","a","m")), fate = factor(next_stage, levels=c("p","j","a","m"))) %>% filter(POPNUM == 231, year == 2) %>% as.data.frame() %>% xtabs(~fate + stage, data = .) %>% as.matrix() Tmat <- Tmat[,-4] #throw away the column for transitions to death N2 <- colSums(Tmat) #get the total number ... CHECK should be before 86 Tmat <- sweep(Tmat[-4,], 1, N2, "/") #normalize to 1, discarding transitions FROM death # figure out how much reproduction happened in year 2 by looking at year 3 # L_elto %>% # mutate(stage = factor(stage, levels=c("p","j","a","m")), # fate = factor(next_stage, levels=c("p","j","a","m"))) %>% # filter(POPNUM == 231, year == 3, stage == "p") %>% # count(stage) # 2 offspring from N[3] == 16 adults Fmat <- matrix(0, nrow=3, ncol=3) # create a matrix full of zeros Fmat[1,3] <- 2/16 # counted 16 adults in this stage, and 2 seedlings in next year TF2 <- list(T = Tmat, F = Fmat) ## ----------------------------------------------------------------------------- TF2 N2 ## ----posterior1--------------------------------------------------------------- Tprior <- matrix(0.25, byrow = TRUE, ncol = 3, nrow=4) fill_transitions(TF, N, P = Tprior) # returns transition matrix by default ## ----manualTransitionPost----------------------------------------------------- Tobs <- sweep(TF$T, 2, N, "*") # get the observed transitions Tobs <- rbind(Tobs, N - colSums(Tobs)) # add the mortality row Tobs <- Tobs + 0.25 # add the prior sweep(Tobs, 2, colSums(Tobs), "/")[-4,] # divide by the column sum, and discard mortalityrow ## ----------------------------------------------------------------------------- alpha <- matrix(c(NA_real_, NA_real_, 1e-5, NA_real_, NA_real_, NA_real_, NA_real_, NA_real_, NA_real_), nrow=3, ncol = 3, byrow = TRUE) beta <- matrix(c(NA_real_, NA_real_, 1e-5, NA_real_, NA_real_, NA_real_, NA_real_, NA_real_, NA_real_), nrow=3, ncol = 3, byrow = TRUE) fill_fertility(TF, N, alpha = alpha, beta = beta) ## ----------------------------------------------------------------------------- obs_offspring <- N[3]*TF$F[1,3] prior_alpha <- 1e-05 prior_beta <- 1e-05 posterior_alpha <- obs_offspring + prior_alpha posterior_beta <- N[3] + prior_beta posterior_alpha / posterior_beta # expected value ## ----uninformativePrior------------------------------------------------------- unif <- list(T = fill_transitions(TF, N), F = fill_fertility(TF, N, alpha = alpha, beta = beta)) unif ## ----------------------------------------------------------------------------- fill_transitions(TF, N, returnType = "TN") ## ----------------------------------------------------------------------------- fill_fertility(TF, N, alpha = alpha, beta = beta, returnType = "ab") ## ----------------------------------------------------------------------------- fill_transitions(TF, N, returnType = "A") ## ----------------------------------------------------------------------------- RLT_Tprior <- matrix(c(0.25, 0.025, 0.0, 0.05, 0.9, 0.025, 0.01, 0.025, 0.95, 0.69, 0.05, 0.025), byrow = TRUE, nrow = 4, ncol = 3) ## ----------------------------------------------------------------------------- fill_transitions(TF, N, P = RLT_Tprior) ## ----------------------------------------------------------------------------- fill_transitions(TF, N, P = RLT_Tprior, priorweight = 0.5) ## ----------------------------------------------------------------------------- TN <- fill_transitions(TF, N, P = RLT_Tprior, priorweight = 0.5, returnType = "TN") a <- TN[,1] # change 1 to 2, 3 etc to get marginal distribution of other columns b <- sum(TN[,1]) - TN[,1]# change 1 to 2, 3 etc to get marginal distribution of other columns p <- a / (a + b) lcl <- qbeta(0.025, a, b) ucl <- qbeta(0.975, a, b) knitr::kable(sprintf("%.3f (%.3f, %.3f)", p, lcl, ucl)) ## ----echo=FALSE--------------------------------------------------------------- TN <- fill_transitions(TF, N, P = RLT_Tprior, priorweight = 9, returnType = "TN") a <- TN[,1] b <- sum(TN[,1]) - TN[,1] # beta parameter is sum of all other dirichlet parameters. p2 <- a / (a + b) lcl2 <- qbeta(0.025, a, b) ucl2 <- qbeta(0.975, a, b) alltogethernow <- tibble(priorweight = factor(rep(c(1,9),each = 4)), fate = factor(rep(c("s","j","a","dead"), times = 2), levels = c("s","j","a","dead")), p = c(p, p2), lcl = c(lcl, lcl2), ucl = c(ucl, ucl2)) ggplot(data = alltogethernow, mapping = aes(x = fate, y = p, group= priorweight)) + geom_point(mapping = aes(shape = priorweight), position = position_dodge(width =0.5)) + geom_errorbar(mapping = aes(ymin = lcl, ymax = ucl), position = position_dodge(width = 0.5)) + rlt_theme ## ----------------------------------------------------------------------------- sim_transitions(TF, N, P = RLT_Tprior, alpha = alpha, beta = beta, priorweight = 0.5) ## ----------------------------------------------------------------------------- set.seed(8390278) # make this part reproducible alpha2 <- matrix(c(NA_real_, NA_real_, 0.025, NA_real_, NA_real_, NA_real_, NA_real_, NA_real_, NA_real_), nrow=3, ncol = 3, byrow = TRUE) beta2 <- matrix(c(NA_real_, NA_real_, 1, NA_real_, NA_real_, NA_real_, NA_real_, NA_real_, NA_real_), nrow=3, ncol = 3, byrow = TRUE) # generate 1000 matrices based on the prior transitions and fertilities plus the data RLT_0.5 <- sim_transitions(TF, N, P = RLT_Tprior, alpha = alpha2, beta = beta2, priorweight = 0.5, samples = 1000) # extract the lambdas for each matrix RLT_0.5 <- tibble(lposterior = map_dbl(RLT_0.5, lambda)) ggplot(data = RLT_0.5, mapping = aes(x = lposterior)) + geom_histogram(binwidth = 0.02) + rlt_theme ## ----------------------------------------------------------------------------- RLT_0.5_summary <- summarize(RLT_0.5, medianL = median(lposterior), meanL = mean(lposterior), lcl = quantile(lposterior, probs = 0.025), ucl = quantile(lposterior, probs = 0.975), pincrease = sum(lposterior > 1.)/n()) knitr::kable(RLT_0.5_summary, digits = 2) ## ----------------------------------------------------------------------------- cri <- transition_CrI(TF, N, stage_names = c("plantula", "juvenile", "adult")) plot_transition_CrI(cri)