Working with a single population and transition
-## Warning: package 'tibble' was built under R version 3.5.2## Warning: package 'tidyr' was built under R version 3.5.2## Warning: package 'purrr' was built under R version 3.5.2## Warning: package 'dplyr' was built under R version 3.5.2## Warning: package 'stringr' was built under R version 3.5.2library(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"))The goal of this vignette is to demonstrate the use of package raretrans on a single population and transition period.
Part I: Obtaining the projection matrix
-raretrans assumes the projection matrix is a list of two matrices, a transition matrix and a fertility matrix. This is the output format of popbio::projection.matrix. If we have individual transitions in a dataframe we can use popbio::projection.matrix to obtain the required data. We demonstrate with the transition and fertility data for the epiphytic orchid Lepanthes elto POPNUM 250 in period 5.
raretrans assumes the projection matrix is a list of two matrices, a transition matrix and a fertility matrix. This is the output format of popbio::projection.matrix. If we have individual transitions in a dataframe we can use popbio::projection.matrix to obtain the required data. We demonstrate with the transition and fertility data for the epiphytic orchid Lepanthes elto POPNUM 250 in period 5.
Step 1: Load and munge the single population data for L. elto
- -## # A tibble: 6 x 15
-## POPNUM year seedlings adults fertility OCCUORNOT REMOVEDORNOT IND_NUM
-## <int> <dbl> <dbl> <dbl> <dbl> <chr> <chr> <int>
-## 1 209 1 1 6 0 Ocupado No removido 67
-## 2 209 2 0 9 0 Ocupado No removido 67
-## 3 209 3 0 9 0 Ocupado No removido 67
-## 4 209 4 0 9 0 Ocupado No removido 67
-## 5 209 1 1 6 0 Ocupado No removido 68
-## 6 209 2 0 9 0 Ocupado No removido 68
-## # … with 7 more variables: stage <chr>, next_stage <chr>,
-## # first_year <dbl>, last_year <dbl>, recruited <lgl>, died <dbl>,
-## # lifespan <int>Each row of this dataframe has columns for current stage, the next stage, and average per individual fertility.
-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## # A tibble: 6 x 13
+## POPNUM year seedlings adults fertility IND_NUM stage next_stage
+## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <chr> <chr>
+## 1 209 1 1 6 0 67 j j
+## 2 209 1 1 6 0 68 a a
+## 3 209 1 1 6 0 69 a a
+## 4 209 1 1 6 0 70 a a
+## 5 209 1 1 6 0 71 j a
+## 6 209 1 1 6 0 72 a a
+## # … with 5 more variables: first_year <dbl>, last_year <dbl>,
+## # recruited <lgl>, died <dbl>, lifespan <int>Each row of this dataframe has columns for current stage, the next stage, and average per individual fertility. Note that “p” stands from “plantula” which means “seedlings” in Spanish. The first set of lines below change the name of the life history stage from “p” to “s” after selecting the population and time period.
+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## $T
##
## s j a
@@ -158,11 +153,11 @@
Step 2: Obtain starting number of individuals per stage
Since our priors are based on counts (number of individuals, N) and the prior equivalent sampling size is expressed as a multiple of the number of individuals observed, we need to obtain the number of individuals in each stage (\(N\)) at the first time period.
We use function raretrans::get_state_vector() to obtain the starting individual count, N.
-N <- get_state_vector(onepop, stage = stage, sort = c("s", "j", "a"))
-N # A vector of # of starting individuals for each stage
+N <- get_state_vector(onepop, stage = stage, sort = c("s", "j", "a"))
+N # A vector of # of starting individuals for each stage
## [1] 11 47 34
-The list of matrices and vector of individual counts do not have to come from a dataframe as we’ve done here. As long as they have the expected format they can be created by hand. We use population 231 in period 2 as an example, splitting the matrix into transition T and fecundity F matrices.
-
+The list of matrices and vector of individual counts do not have to come from a dataframe as we’ve done here. As long as they have the expected format they can be created by hand. We use population 231 in period 2 as an example, splitting the matrix into transition T and fecundity F matrices. Below, “m” stands for “muerte” that is plants that are dead.
+
## $Tmat
## stage
## fate p j a
@@ -175,7 +170,7 @@
## [1,] 0 0 0.125
## [2,] 0 0 0.000
## [3,] 0 0 0.000
-
+
## p j a
## 2 6 16
This matrix is missing the transition from seedling to juvenile, and none of the 6 juveniles died leading to an overestimate of survival. The observed asymptotic population growth rate is \(\lambda =\) 0.88. The matrix is not Ergodic (you can’t get to every other state from one or more states), and reducible, meaning one or more columns and rows can be discarded and have the same eigen properties.
@@ -192,46 +187,46 @@
Transition matrix
So, let’s add a uniform dirichlet prior with weight = \(1\) to the transition matrix, \(T\). Here, we have 4 fates (3 + death), so each fate adds 0.25 to the matrix of observed fates (not the transition matrix!). When we specify a prior matrix for the transitions, there is one more row than columns. This extra row represents death.
-Tprior <- matrix(0.25, byrow = TRUE, ncol = 3, nrow = 4)
-fill_transitions(TF, N, P = Tprior) # returns transition matrix by default
+Tprior <- matrix(0.25, byrow = TRUE, ncol = 3, nrow = 4)
+fill_transitions(TF, N, P = Tprior) # returns transition matrix by default
## [,1] [,2] [,3]
## [1,] 0.10416667 0.005208333 0.007142857
## [2,] 0.60416667 0.567708333 0.007142857
## [3,] 0.02083333 0.296875000 0.835714286
We can get the same result ‘by hand’ - we need the vector of observations because the posterior is calculated from the observed transitions, not the transition matrix.
-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
+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
## s j a
## s 0.10416667 0.005208333 0.007142857
## j 0.60416667 0.567708333 0.007142857
## a 0.02083333 0.296875000 0.835714286
-The uniform prior fills in the missing transitions, but also creates problems because it provides transition values that are biologically impossible. For example, it provides a transition for adult->seedling, when this transition is only possible in the fecundity matrix \(F\).
+The uniform prior fills in the missing transitions, but also creates problems because it provides transition values that are biologically impossible. For example, it provides a transition for adult->seedling, when this transition is only possible in the fecundity matrix \(F\). For this reason we do not recommend using uniform priors.
-Fecundity matrix
-We must specify the parameters for the fecundity prior on the non-reproducing stages as missing, else the result will fill in the entire matrix.
- +Fertility matrix +We must specify the parameters for the fertility prior on the non-reproducing stages as missing, else the result will fill in the entire matrix.
+## [,1] [,2] [,3]
## [1,] 0 0 0.1176473
## [2,] 0 0 0.0000000
## [3,] 0 0 0.0000000The change in the fecundity is < 0.0001 compared to the observed value. Calculating by hand, prior alpha is the number of observed offspring, and prior beta is the number of observed adults.
-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 valueThe change in the fertility is < 0.0001 compared to the observed value. Calculating by hand, prior alpha is the number of observed offspring, and prior beta is the number of observed adults.
+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## [1] 0.1176473This demonstrates why the posterior point estimate of fecundity does not change much; the non-informative values for \(\alpha\) and \(\beta\) barely change the observed values.
Now we can put them together.
-unif <- list(T = fill_transitions(TF, N), F = fill_fecundity(TF, N, alpha = c(NA_real_,
- NA_real_, 1e-05), beta = c(NA_real_, NA_real_, 1e-05)))
-unifunif <- list(T = fill_transitions(TF, N), F = fill_fertility(TF, N, alpha = c(NA_real_,
+ NA_real_, 1e-05), beta = c(NA_real_, NA_real_, 1e-05)))
+unif## $T
## [,1] [,2] [,3]
## [1,] 0.10416667 0.005208333 0.007142857
@@ -245,16 +240,16 @@
## [3,] 0 0 0.0000000
The asymptotic population growth rate is now \(\lambda =\) 0.92. The growth rate shrinks slightly because applying the uniform prior to the transition probabilities causes the observed growth and survival transitions to shrink slightly relative to the unobserved transitions.
-Other options for returnType
+Other options for argument returnType
-There are three other returnTypes:
By default fill_transitions() returns the transition matrix \(T\), and fill_fertility() returns the fertility matrix \(F\). There are three other values that the argument returnType can take:
-
-
fill_transitions()returns an augmented matrix of fates, which is useful for simulation. The fourth row of this result (see below) is the mortality state.
+
fill_transitions(... returnType = "TN") can return an augmented matrix of fates, which is useful for simulation. The fourth row of this result (see below) is the mortality state.
## [,1] [,2] [,3]
## [1,] 1.25 0.25 0.25
## [2,] 7.25 27.25 0.25
@@ -262,10 +257,10 @@
## [4,] 3.25 6.25 5.25
-
-
fill_fecundity()returns the alpha and beta vectors of the posterior.
+
fill_fertility(... returnType = "ab") returns the alpha and beta vectors of the posterior.
fill_fecundity(TF, N, alpha = c(NA_real_, NA_real_, 1e-05), beta = c(NA_real_,
- NA_real_, 1e-05), returnType = "ab")fill_fertility(TF, N, alpha = c(NA_real_, NA_real_, 1e-05), beta = c(NA_real_,
+ NA_real_, 1e-05), returnType = "ab")## $alpha
## s j a
## NA NA 4.00001
@@ -273,9 +268,9 @@
## $beta
## [1] NA NA 34.00001
-
-
- Both functions can return the full matrix (not very useful) +
- Both functions can also return the full matrix, the sum of \(T\) and \(F\).
##
## s j a
## s 0.104166667 0.005208333 0.124789916
@@ -287,16 +282,16 @@
Incorporate informative priors
To fix the problem of creating impossible transitions, we specify a more informative prior elicited from an expert on epiphytic orchids (RLT). It still has to have the same shape, one more row than columns.
-RLT_Tprior <- matrix(c(0.25, 0.025, 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)
+RLT_Tprior <- matrix(c(0.25, 0.025, 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)
Note that the 1st row, 3rd column is 0.0, because this transition is impossible. This prior is constructed so the columns sum to 1, which creates the greatest flexibility for weighting the prior. By default the weight is 1, interpreted as a prior sample size of 1.
-
+
## [,1] [,2] [,3]
## [1,] 0.1041666667 0.0005208333 0.0000000000
## [2,] 0.5875000000 0.5812500000 0.0007142857
## [3,] 0.0008333333 0.2921875000 0.8557142857
We can specify the weight as a multiple of the sample size for each stage.
-
+
## [,1] [,2] [,3]
## [1,] 0.143939394 0.008333333 0.000000000
## [2,] 0.440909091 0.682978723 0.008333333
@@ -304,36 +299,36 @@
Here the prior is weighted half as much as the observed number of transitions. In this case, with only 2 transitions the prior effective sample size is still 1. If the number of observed transitions was larger, a prior weight of 0.5N would be larger than 1, but still allow the data to dominate.
-Part III: Obtain Credible Intervals
+Part 3: Obtain Credible IntervalsObtain CI for indvidual matrix entries
The marginal posterior distribution of one element in a multinomial is a beta distribution, and we use this to get credible intervals on our transition rates. We can use the TN return type to get the parameters of the desired multinomial.
-TN <- fill_transitions(TF, N, P = RLT_Tprior, priorweight = 0.5, returnType = "TN")
-alpha <- TN[, 1]
-beta <- sum(TN[, 1]) - TN[, 1]
-p <- alpha/(alpha + beta)
-lcl <- qbeta(0.025, alpha, beta)
-ucl <- qbeta(0.975, alpha, beta)
-knitr::kable(sprintf("%.2f (%.2f, %.2f)", p, lcl, ucl))TN <- fill_transitions(TF, N, P = RLT_Tprior, priorweight = 0.5, returnType = "TN")
+alpha <- TN[, 1] # change 1 to 2, 3 etc to get marginal distribution of other columns
+beta <- sum(TN[, 1]) - TN[, 1] # change 1 to 2, 3 etc to get marginal distribution of other columns
+p <- alpha/(alpha + beta)
+lcl <- qbeta(0.025, alpha, beta)
+ucl <- qbeta(0.975, alpha, beta)
+knitr::kable(sprintf("%.3f (%.3f, %.3f)", p, lcl, ucl))| x | |
|---|---|
| 0.14 (0.02, 0.34) | +0.144 (0.025, 0.343) |
| 0.44 (0.22, 0.68) | +0.441 (0.218, 0.677) |
| 0.00 (0.00, 0.04) | +0.003 (0.000, 0.038) |
| 0.41 (0.19, 0.65) | +0.412 (0.195, 0.649) |
Credible intervals on \(\lambda\)
Obtaining credible intervals on the asymptotic growth rate, \(\lambda\), requires simulating matrices from the posterior distributions. This is somewhat complicated to do correctly, and we have written a function raretrans::sim_transitions() to generate a list of simulated matrices given the observed matrix and prior specifications.
sim_transitions(TF, N, P = RLT_Tprior, alpha = c(NA, NA, 0.025), beta = c(NA,
- NA, 1), priorweight = 0.5)sim_transitions(TF, N, P = RLT_Tprior, alpha = c(NA, NA, 0.025), beta = c(NA,
+ NA, 1), priorweight = 0.5)## [[1]]
-## [,1] [,2] [,3]
-## [1,] 4.998590e-02 0.01503749 0.02148950
-## [2,] 5.646228e-01 0.68159789 0.01823104
-## [3,] 3.546484e-17 0.20384823 0.89981659Now we simulate 1000 times, calculate the leading eigenvalue of each matrix, and draw a histogram.
-set.seed(8390278) # make this part reproducible
-RLT_0.5 <- sim_transitions(TF, N, P = RLT_Tprior, alpha = c(NA, NA, 0.025),
- beta = c(NA, NA, 1), priorweight = 0.5, samples = 1000)
-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)set.seed(8390278) # make this part reproducible
+# generate 1000 matrices based on the prior transitions and fertilities plus
+# the data
+RLT_0.5 <- sim_transitions(TF, N, P = RLT_Tprior, alpha = c(NA, NA, 0.025),
+ beta = c(NA, NA, 1), 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
We can calculate some summary statistics too.
-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)We can calculate some summary statistics too. pincrease is the probability that \(\lambda > 1\).
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)| medianL | @@ -404,7 +403,7 @@
|---|
Data was collected from 1850 individuals which were permanently identified and tracked over a 6 year period in intervals of 6 months. Thus transitions are for six month intervals and not yearly transitions as is more frequently published. The life history stages included in this analysis are seedlings (individuals without a lepanthiform sheet on any of the leaves), juveniles (individuals with no evidence of present or past reproductive effort, the base of the inflorescence are persistent), and adult (with presence of active or inactive inflorescences). The average population growth rate over all years and populations is 0.86, suggesting a rapid decline.
-plotdata <- L_elto %>% group_by(POPNUM, year, stage) %>%
- tally() %>%
- group_by(POPNUM, year) %>%
- summarise(N = sum(n))
-
-gg1 <- ggplot(plotdata, aes(x=year, y=N, group=POPNUM)) +
- geom_line() + rlt_theme
-plotsumm <- plotdata %>% group_by(year) %>%
- summarize(N=sum(N))
-gg2 <- ggplot(plotsumm, aes(x=year, y=N)) +
- geom_line() + rlt_theme
-ggboth <- gridExtra::arrangeGrob(gg1, gg2)
-grid::grid.draw(ggboth)plotdata <- L_elto %>% group_by(POPNUM, year, stage) %>%
+ tally() %>%
+ group_by(POPNUM, year) %>%
+ summarise(N = sum(n))
+
+gg1 <- ggplot(plotdata, aes(x=year, y=N, group=POPNUM)) +
+ geom_line() + rlt_theme
+plotsumm <- plotdata %>% group_by(year) %>%
+ summarize(N=sum(N))
+gg2 <- ggplot(plotsumm, aes(x=year, y=N)) +
+ geom_line() + rlt_theme
+ggboth <- gridExtra::arrangeGrob(gg1, gg2)
+grid::grid.draw(ggboth)
Population size over time.
OK, now we want to generate matrices for each year and popnum. We have to change the variables into factors and control the levels so that the resulting matrices end up with the right shape with zeros for the missing transitions. Also have to coerce each grouped tbl_df to a regular data.frame using as.data.frame() before passing to projection.matrix().
allA <- L_elto %>%
- mutate(stage = factor(stage, levels=c("p","j","a","m")),
- fate = factor(next_stage, levels=c("p","j","a","m"))) %>%
- group_by(POPNUM, year) %>%
- do(A = projection.matrix(as.data.frame(.),
- stage="stage", fate="fate",
- fertility="fertility", sort=c("p","j","a")),
- TF = projection.matrix(as.data.frame(.),
- stage="stage", fate="fate",
- fertility="fertility", sort=c("p","j","a"), TF = TRUE),
- N = get_state_vector(as.data.frame(.),
- stage="stage", sort=c("p","j","a"))) %>%
-
- filter(year < 12) # last time period all zero for obvious reasonOK, now we want to generate matrices for each year and popnum. We have to change the variables into factors and control the levels so that the resulting matrices end up with the right shape with zeros for the missing transitions. Also have to coerce each grouped tbl_df to a regular data.frame using as.data.frame() before passing to projection.matrix().
allA <- L_elto %>%
+ mutate(stage = factor(stage, levels=c("p","j","a","m")),
+ fate = factor(next_stage, levels=c("p","j","a","m"))) %>%
+ group_by(POPNUM, year) %>%
+ do(A = projection.matrix(as.data.frame(.),
+ stage="stage", fate="fate",
+ fertility="fertility", sort=c("p","j","a")),
+ TF = projection.matrix(as.data.frame(.),
+ stage="stage", fate="fate",
+ fertility="fertility", sort=c("p","j","a"), TF = TRUE),
+ N = get_state_vector(as.data.frame(.),
+ stage="stage", sort=c("p","j","a"))) %>%
+
+ filter(year < 12) # last time period all zero for obvious reasonWe end up with all the matrices for each population and year in a list-column of the result. Working with this thing is a bit awkward. The best strategy is from Jenny Bryans’ talk on list columns and uses mutate with map_*(listcolumn, function).
library(popdemo)
-allA <- ungroup(allA) %>%
- mutate(A = map(A, matrix, nrow=3, ncol=3,
- dimnames=list(c("p","j","a"),c("p","j","a")))) %>%
- mutate(lmbda = map_dbl(A, lambda),
- ergodic = map_dbl(A, isErgodic),
- irreduc = map_dbl(A, isIrreducible),
- primitv = map_dbl(A, isPrimitive))
-ggplot(allA, aes(x=lmbda)) + geom_histogram(bins = 25) + facet_wrap(~year)+
- ylab("Frequency")+
- xlab(expression(paste(lambda, ", asymptotic population growth"))) + rlt_themelibrary(popdemo)
+allA <- ungroup(allA) %>%
+ mutate(A = map(A, matrix, nrow=3, ncol=3,
+ dimnames=list(c("p","j","a"),c("p","j","a")))) %>%
+ mutate(lmbda = map_dbl(A, lambda),
+ ergodic = map_dbl(A, isErgodic),
+ irreduc = map_dbl(A, isIrreducible),
+ primitv = map_dbl(A, isPrimitive))
+ggplot(allA, aes(x=lmbda)) + geom_histogram(bins = 25) + facet_wrap(~year)+
+ ylab("Frequency")+
+ xlab(expression(paste(lambda, ", asymptotic population growth"))) + rlt_theme
histograms of asymptotic population growth for each year.
Now, we need to figure out how many matrices are missing transitions. The full matrix has observed transitions in every cell.
-allA <- ungroup(allA) %>% mutate(missing = map(A, ~which(.x==0)),
- n_missing = map_dbl(missing, length))
-missing_summary <- summary(allA$n_missing)
-ergo_irr <- as.data.frame(with(allA, table(ergodic, irreduc)))
-knitr::kable(ergo_irr, caption = "Ergodicity and irreducibility of individual transition matrices.")allA <- ungroup(allA) %>% mutate(missing = map(A, ~which(.x==0)),
+ n_missing = map_dbl(missing, length))
+missing_summary <- summary(allA$n_missing)
+ergo_irr <- as.data.frame(with(allA, table(ergodic, irreduc)))
+knitr::kable(ergo_irr, caption = "Ergodicity and irreducibility of individual transition matrices.")
| 238 |
sgr_unif <- testing %>% split(.$POPNUM) %>%
- map(~stoch.growth.rate(.x$Adirch, verbose = FALSE))
-
-sgr_unif <- bind_rows(map(sgr_unif, xtrc), .id="POPNUM")
-compare_approx <- tibble(unif = sgr_unif$approx,
- ignored = sgr$approx)
-ggplot(compare_approx, aes(x=ignored, y = unif)) +
- geom_point() +
- geom_abline(slope = 1, intercept=0) +
- xlab(expression(paste("log(",lambda,") ignoring zeros"))) +
- ylab(expression(paste("log(",lambda,") with uniform prior"))) +
- rlt_themesgr_unif <- testing %>% split(.$POPNUM) %>%
+ map(~stoch.growth.rate(.x$Adirch, verbose = FALSE))
+
+sgr_unif <- bind_rows(map(sgr_unif, xtrc), .id="POPNUM")
+compare_approx <- tibble(unif = sgr_unif$approx,
+ ignored = sgr$approx)
+ggplot(compare_approx, aes(x=ignored, y = unif)) +
+ geom_point() +
+ geom_abline(slope = 1, intercept=0) +
+ xlab(expression(paste("log(",lambda,") ignoring zeros"))) +
+ ylab(expression(paste("log(",lambda,") with uniform prior"))) +
+ rlt_theme
So using a Dirichlet prior to fill in the gaps works, and doesn’t appear to distort the stochastic population growth rate; all populations still have negative growth. Most populations end up with stochastic growth rates that are lower than with the problematic matrices. All resulting matrices are now ergodic and irreducible.
Fill in using uninformative priors on fertility
-The other component of the matrix is the fertility. We do this by using a Gamma prior with the function fill_fecundity(). This function takes the matrix as a list of 2 components, a matrix of transition probabilities \(T\) and a matrix of fertility contributions \(F\). The Gamma prior only applies to the \(F\) matrix. If not specified, the function assumes an uniformative prior with \(\alpha = 0.00001, \beta = 0.00001\) across all stages. Stages that do not reproduce are marked with NA_real_ in the prior probability vectors. Again, demonstrate with a single example. The posterior \(\alpha\) is given by the sum of the prior \(\alpha = 0.00001\) and the number of offspring next year (2 in this case), while the posterior \(\beta\) is simply the sum of the prior and the number of adults contributing to the observed reproduction.
The other component of the matrix is the fertility. We do this by using a Gamma prior with the function fill_fertility(). This function takes the matrix as a list of 2 components, a matrix of transition probabilities \(T\) and a matrix of fertility contributions \(F\). The Gamma prior only applies to the \(F\) matrix. If not specified, the function assumes an uniformative prior with \(\alpha = 0.00001, \beta = 0.00001\) across all stages. Stages that do not reproduce are marked with NA_real_ in the prior probability vectors. Again, demonstrate with a single example. The posterior \(\alpha\) is given by the sum of the prior \(\alpha = 0.00001\) and the number of offspring next year (2 in this case), while the posterior \(\beta\) is simply the sum of the prior and the number of adults contributing to the observed reproduction.
The expected value for the posterior \(F\) is \(\alpha / \beta\).
- +##
## p j a
## p 0.0000000 0.0000000 0.1250005
## j 0.0000000 0.0000000 0.0000000
## a 0.0000000 0.0000000 0.0000000## [1] 1.055885## [1] FALSE## [1] FALSEThe uninformative prior on fecundity makes very little difference to lambda, and does not help with issues of ergodicity and irreducibility.
-testing <- allA %>% mutate(Agamma = map2(TF, N, fill_fecundity,
- alpha = c(NA_real_, NA_real_, 0.00001),
- beta = c(NA_real_, NA_real_, 0.00001),
- priorweight = 1, returnType = "A"),
- lgamma = map_dbl(Agamma, lambda),
- irrgamma = map_dbl(Agamma, isIrreducible),
- erggamma = map_dbl(Agamma, isErgodic))
-
-ggplot(testing, aes(x=lmbda, y=lgamma)) + geom_point() +
- geom_abline(slope = 1, intercept=0) +
- xlab(expression(paste(lambda, " ignoring zeros"))) +
- ylab(expression(paste(lambda, " using an uninformative prior"))) +
- rlt_themeThe uninformative prior on fertility makes very little difference to lambda, and does not help with issues of ergodicity and irreducibility.
+testing <- allA %>% mutate(Agamma = map2(TF, N, fill_fertility,
+ alpha = c(NA_real_, NA_real_, 0.00001),
+ beta = c(NA_real_, NA_real_, 0.00001),
+ priorweight = 1, returnType = "A"),
+ lgamma = map_dbl(Agamma, lambda),
+ irrgamma = map_dbl(Agamma, isIrreducible),
+ erggamma = map_dbl(Agamma, isErgodic))
+
+ggplot(testing, aes(x=lmbda, y=lgamma)) + geom_point() +
+ geom_abline(slope = 1, intercept=0) +
+ xlab(expression(paste(lambda, " ignoring zeros"))) +
+ ylab(expression(paste(lambda, " using an uninformative prior"))) +
+ rlt_theme
Asymptotic growth rate using a uniform prior with a total weight of 1 vs. the asymptotic growth rate for the observed data.
The matrices with \(\lambda = 0\) typically have very few individuals. Matrices with \(\lambda = 1\) usually have transitions only on the main diagonal. The uninformative Gamma prior does not change \(\lambda\), and does not improve the ergodicity and irreducibility issue. When there is at least 1 adult present, the expected value of Fecundity is very close to zero. When there are no adults present, all those matrices cannot reach the adult stage and so the expected value of Fecundity does not affect \(\lambda\). However, this is not a necessary condition.
-ergo_irr <- as.data.frame(with(testing, table(erggamma, irrgamma)))
-
-knitr::kable(ergo_irr, caption = "Ergodicity and irreducibility of individual transition matrices.")The matrices with \(\lambda = 0\) typically have very few individuals. Matrices with \(\lambda = 1\) usually have transitions only on the main diagonal. The uninformative Gamma prior does not change \(\lambda\), and does not improve the ergodicity and irreducibility issue. When there is at least 1 adult present, the expected value of fertility is very close to zero. When there are no adults present, all those matrices cannot reach the adult stage and so the expected value of fertility does not affect \(\lambda\). However, this is not a necessary condition.
+ergo_irr <- as.data.frame(with(testing, table(erggamma, irrgamma)))
+
+knitr::kable(ergo_irr, caption = "Ergodicity and irreducibility of individual transition matrices.")
| TF | +A list of two matrices, T and F, as ouput by |
+
|---|---|
| N | +A vector of observed stage distribution. |
+
| alpha | +A vector of the prior parameter for each stage. Stages that can't reproduce are NA_real_ |
+
| beta | +A vector of the prior parameter for each stage. Stages that can't reproduce are NA_real_ |
+
| priorweight | +total weight for each column of prior as a percentage of sample size or 1 if negative |
+
| returnType | +A character vector describing the desired return value. Defaults to "F" the fertility matrix |
+
Value
+ +The return value depends on parameter returnType.
-
+
A - the summed matrix "filled in" using a Gamma prior
+F - just the filled in fertility matrix
+ab - the posterior parameters alpha and beta as a list.
+
Details
+ +Assumes that only one stage reproduces ... needs generalizing.
+ + +