April 2021
library(shiny)
knitr::opts_chunk$set(echo = FALSE, warning = FALSE, message = FALSE)
enableBookmarking(store = "url")
varnames <- rbind(c("Transition ages","age_trans"),
c("Natural survival", "p_s"),
c("Harvest rate","h_s"),
c("Fecundities", "births"))
library(ggplot2)
use_plotly <- require(plotly)
survparams <- function(nstages, age_transition, p_stage, h_stage){
x <- (nstages != length(age_transition)) | (nstages != length(p_stage)) |
(nstages != length(h_stage))
if(any(x)) stop("Number of stages and lengths of demographic parameters don't match")
#Checks that age_transition, p_stage, and h_stage have the same length as nstages.
If the lengths mismatch, the function stops and throws an error.
P <- rep(NA, nstages)
G <- rep(NA, nstages)
for (i in 1:nstages){
p_natural_instant <- log(p_stage[i]) + log(1-h_stage[i])
# Computes the instantaneous survival probability, considering natural survival and harvest rates.
p_natural <- exp(p_natural_instant)
# Converts the instantaneous survival back to a regular probability using the exponential function.
if (i>1){
years_in_stage <- age_transition[i] - age_transition[i-1]
} else {
years_in_stage <- age_transition[i]
}
# For stages after the first (i > 1), duration is the difference between the current and previous transition ages.
For the first stage, duration is simply the first transition age.
P[i] <- ((1-p_natural^(years_in_stage-1))/(1-p_natural^(years_in_stage)))*p_natural
# P[i]: Proportion of individuals surviving within the current stage.
Accounts for survival across multiple years, adjusting for the probability of dying during the stage.
G[i] <- ((p_natural^(years_in_stage)) * (1-p_natural))/(1-p_natural^(years_in_stage))
# G[i]: Proportion od individuals transitionating to the next stage.
Accounts for the probability of surviving all years in the stage but not staying in the stage indefinitely, transitionating to the next one.
}
return(data.frame(P = P, G = G))
# Outputs a data frame containing the survival (P) and growth (G) values for all stages.
}
createA <- function(nstages, births, age_trans, p_s, h_s){
x <- survparams(nstages, age_trans, p_s, h_s)
# Calls survparams to compute survival (P) and growth (G) values for each stage.
A <- diag(nstages)*0
# Creates a square matrix (nstages x nstages) filled with zeros.
A[1,] <- births
diag(A) <- x$P
# A[1, ]: Assigns birth rates to the first row.
diag(A): Fills the diagonal with survival probabilities (P).
irow <- 2:nstages
icol <- 1:(nstages-1)
rowcolind <- cbind(irow, icol)
A[rowcolind] <- x$G[1:(nstages - 1)]
# Assigns growth probabilities (G) to the subdiagonal of the matrix.
colnames(A) <- paste("Age ", 1:nstages)
A
}
createX <- function(input, A){
xtemp <- matrix(0,
nrow = input$nstages, ncol = 3)
xtemp[,1] <- asnum(input$stage_sizes)
xtemp
# Creates a matrix xtemp with rows equal to nstages and 3 columns.
Populates the first column with initial population sizes (stage_sizes), parsed using asnum.
}
asnum <- function(x){
as.numeric(unlist(strsplit(x, ",")))
# Converts a comma-separated string into a numeric vector.
}
leslie <- function(tmax, Ninit, A){
if (length(Ninit) != nrow(A)) stop("Dimension mismatch between initial abundances and Leslie matrix")
N <- matrix(0, nrow = length(Ninit), ncol = tmax)
N[,1] <- Ninit
for (t in 1:(tmax-1)){
N[,t+1] <- A %*% N[,t]
}
N
# Multiplies the Leslie matrix A by the current population vector to calculate the next population vector.
}
inputPanel(
numericInput("nstages", label = "Number of stages",
value = 3,
min = 1, max = 20, step = 1),
textInput("age_trans", label = "Transition ages",
value = "1, 15, 50"),
textInput("p_s", label = "Natural survival probablities",
value = "0.5, 0.73, 0.97"),
textInput("h_s", label = "Harvest mortality rates",
value = c("0,0,0")),
textInput("births", label = "Fecundities",
value = "0,0,2"),
textInput("stage_sizes", label = "Stage abundances",
value = "1500, 4890, 750"),
numericInput("tmax", label = "Number of years",
value = 50,
min = 1, max = 200, step = 1),
numericInput("yr_init", label = "Initial year",
value = 2017,
min = 0, max = 2100, step = 1),
checkboxInput("logged", label = "Log y-axis?",
value = FALSE),
actionButton("create_leslie", label = "Run simulation"),
bookmarkButton()
)
ptable <- eventReactive(input$create_leslie, {
A <- createA(input$nstages, asnum(input$births), asnum(input$age_trans),
asnum(input$p_s), asnum(input$h_s))
A
})
xtable <- eventReactive(input$create_leslie, {
createX(input, ptable())
})
renderTable({
ptable()
}, digits = 3)
renderText({
paste0("The finite rate of growth for your Leslie matrix is ",
signif(eigen(ptable())$values[1],4), ". \n",
"Estimated as the leading eigenvalue of the
Leslie matrix. ")
})
renderTable({
xtemp <- xtable()
xtemp[,2] <- ptable() %*% xtemp[,1]
xtemp[,3] <- ptable() %*% xtemp[,2]
colnames(xtemp) <- paste("Year", c(1,2,3))
xtemp
})
N1 <- eventReactive(input$create_leslie, {
leslie(input$tmax, xtable()[,1], ptable())
})
p1 <- eventReactive(input$create_leslie, {
dplot <- data.frame(years = (1:input$tmax) + input$yr_init - 1,
Ntot = colSums(N1())
)
#Make long format data
Ndat <- data.frame(years = rep(dplot$years, each = input$nstages),
stage = factor(rep(1:input$nstages, input$tmax)),
N = as.numeric(N1()))
g1 <- ggplot(dplot, aes(x = years, y = Ntot)) +
geom_line() +
geom_point() +
ylab("Abundance") +
theme_classic()
g2 <- ggplot(Ndat, aes(x = years, y = N, color = stage)) +
geom_line() +
geom_point() +
ylab("Abundance") +
theme_classic() +
scale_color_brewer(palette = "Dark2")
# cat(file=stderr(), "\n",which(Ndat$years %in% baryears), "stopped", "\n")
if (input$logged){
g1 <- g1 + scale_y_log10()
g2 <- g2 + scale_y_log10()
}
return(list(g1 = g1, g2= g2))
})
if (use_plotly){
renderPlotly({
ggplotly(p1()[[1]])
})
} else {
renderPlot({
p1()[[1]]
})
}
if (use_plotly){
renderPlotly({
ggplotly(p1()[[2]])
})
} else {
renderPlot({
p1()[[2]]
})
}
Select one parameters to vary and which lifestage the parameter should be varied for. For the first parameter enter minimum and maximum values, separated by a comma. The app will plot the total population size, with the different parameter values as different colours.
inputPanel(
selectInput("param1", label = "First parameter to vary",
choices = c("Transition ages" = "age_trans",
"Natural survival" = "p_s",
"Harvest rate" = "h_s",
"Fecundities" = "births"),
selected = "p_s"),
numericInput("param1_stage", label = "Parameter 1 life stage",
value = 1),
numericInput("param1_n", label = "Number of values to simulate",
value = 5),
textInput("param1_values", label = "Parameter 1 values",
value ="0.5, 0.9"),
actionButton("run_model_sens1", label = "Run simulations")
)
psens <- eventReactive(input$run_model_sens1, {
nstages <- input$nstages
births <- asnum(input$births)
age_trans <- asnum(input$age_trans)
p_s <- asnum(input$p_s)
h_s <- asnum(input$h_s)
param1range <- asnum(input$param1_values)
param1vals <- seq(param1range[1], param1range[2],
length.out = input$param1_n)
Ndat <- NULL
for (i in 1:input$param1_n){
if (input$param1 == "p_s") p_s[input$param1_stage] <- param1vals[i]
if (input$param1 == "h_s") h_s[input$param1_stage] <- param1vals[i]
if (input$param1 == "age_trans") age_trans[input$param1_stage] <- round(param1vals[i])
if (input$param1 == "births") births[input$param1_stage] <- param1vals[i]
Atemp <- createA(nstages, births, age_trans,
p_s, h_s)
Xtemp <- createX(input, ptable())
Ntemp <- leslie(input$tmax, Xtemp[,1], Atemp)
dplot <- data.frame(years = (1:input$tmax) + input$yr_init - 1,
Ntot = colSums(Ntemp),
Parameter_value = param1vals[i])
Ndat <- c(Ndat, list(dplot))
}
Nout <- do.call("rbind", Ndat)
var1 <- varnames[match(input$param1, varnames[,2]),1]
g1 <- ggplot(Nout, aes(x = years, y = Ntot, color = Parameter_value)) +
geom_line() +
geom_point() +
ylab("Abundance") +
xlab("Year") +
theme_classic() +
labs(color = var1)
return(g1)
})
if (use_plotly){
renderPlotly({
ggplotly(psens())
})
} else {
renderPlot({
psens()
})
}
Select two parameters to vary and which lifestage the parameters should be varied for. For the first parameter enter minimum and maximum values, separated by a comma. For the second parameter, enter as many discrete values as you wish, seperated by commas. The app will plot the finite growth rate, with the first parameter on the x-axis and the second parameter as different colours.
inputPanel(
selectInput("param2", label = "First parameter to vary",
choices = c("Transition ages" = "age_trans",
"Natural survival" = "p_s",
"Harvest rate" = "h_s",
"Fecundities" = "births"),
selected = "p_s"),
numericInput("param2_stage", label = "Parameter 1 life stage",
value = 1),
numericInput("param2_n", label = "Number of values to simulate",
value = 30),
textInput("param2_values", label = "Parameter 1 values",
value ="0.5, 0.9"),
selectInput("param3", label = "Second parameter to vary",
choices = c("Transition ages" = "age_trans",
"Natural survival" = "p_s",
"Harvest rate" = "h_s",
"Fecundities" = "births"),
selected = "births"),
numericInput("param3_stage", label = "Parameter 2 life stage",
value = 3),
textInput("param3_values", label = "Parameter 2 values",
value = "2, 4, 6"),
actionButton("run_model_sens2", label = "Run simulations")
)
psens2 <- eventReactive(input$run_model_sens2, {
nstages <- input$nstages
births <- asnum(input$births)
age_trans <- asnum(input$age_trans)
p_s <- asnum(input$p_s)
h_s <- asnum(input$h_s)
param2range <- asnum(input$param2_values)
param2vals <- seq(param2range[1], param2range[2],
length.out = input$param2_n)
# cat(file=stderr(), "x", param1vals, "\n")
param3vals <- asnum(input$param3_values)
nparam3 <- length(param3vals)
datparams <- NULL
for (i in 1:input$param2_n){
if (input$param2 == "p_s") p_s[input$param2_stage] <- param2vals[i]
if (input$param2 == "h_s") h_s[input$param2_stage] <- param2vals[i]
if (input$param2 == "age_trans") age_trans[input$param2_stage] <- round(param2vals[i])
if (input$param2 == "births") births[input$param2_stage] <- param2vals[i]
for (j in 1:nparam3){
if (input$param3 == "p_s") p_s[input$param3_stage] <- param3vals[j]
if (input$param3 == "h_s") h_s[input$param3_stage] <- param3vals[j]
if (input$param3 == "age_trans") age_trans[input$param3_stage] <- round(param3vals[j])
if (input$param3 == "births") births[input$param3_stage] <- param3vals[j]
A <- createA(nstages, births, age_trans,
p_s, h_s)
dattemp <- data.frame(param2 = param2vals[i],
param3 = param3vals[j],
R = Re(eigen(A)$values[1]))
datparams <- c(datparams, list(dattemp))
}
}
dout <- do.call("rbind", datparams)
var2 <- varnames[match(input$param2, varnames[,2]),1]
var3 <- varnames[match(input$param3, varnames[,2]),1]
g1 <- ggplot(dout, aes(x = param2, y = R, color = param3)) +
geom_point() +
ylab("Finite growth rate per year") +
xlab(var2) +
theme_classic() +
labs(color = var3) +
geom_hline(yintercept = 1)
return(g1)
})
if (use_plotly){
renderPlotly({
ggplotly(psens2())
})
} else {
renderPlot({
psens2()
})
}