Newer
Older
param$fatalities_se <- c(0, input$fatalities_mat_cumulated[,2])
param$onset_year <- c(min(input$fatalities_mat_cumulated[,3]), input$fatalities_mat_cumulated[,3])
param$onset_time <- param$onset_year - min(param$onset_year) + 1
} # end (if1)
}) # end observe
###~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~###
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
# Make sure fatalities are expressed as "number" (not rate) for the run_simul function
se_prod2 <- function(mu1, se1, mu2, se2) sqrt((se1^2 * se2^2) + (se1^2 * mu2^2) + (mu1^2 * se2^2))
observeEvent({
input$run
},{
if(input$fatalities_unit == "h"){
pop_size_tot <- sum(pop_vector(pop_size = param$pop_size_mean, pop_size_type = param$pop_size_type, s = param$survivals, f = param$fecundities)[-1])
param$fatalities_mean_nb <- (param$fatalities_mean/100) * pop_size_tot
param$fatalities_se_nb <- se_prod2(mu1 = param$fatalities_mean/100,
se1 = param$fatalities_se/100,
mu2 = pop_size_tot,
se2 = (pop_size_tot/param$pop_size_mean) * param$pop_size_se)
}else{
param$fatalities_mean_nb <- param$fatalities_mean
param$fatalities_se_nb <- param$fatalities_se
}
})
#################################
## Population size
##-------------------------------
observe({
# Case 1 : Values from expert elicitation
if(input$pop_size_input_type == "eli_exp"){
if(!(is.null(param$pop_size_eli_result))){
param$pop_size_mean <- round(param$pop_size_eli_result$mean)
param$pop_size_se <- round(param$pop_size_eli_result$SE)
ready$pop_size <- TRUE
} else {
ready$pop_size <- FALSE
}
} else {
if(input$pop_size_input_type == "val"){
# Case 2 : Values directly provided as mean & SE
ready$pop_size <- TRUE
param$pop_size_mean <- input$pop_size_mean
param$pop_size_se <- input$pop_size_se
}else{
# Case 3 : Values directly provided as lower/upper interval
ready$pop_size <- TRUE
param$pop_size_mean <- round(get_mu(lower = input$pop_size_lower, upper = input$pop_size_upper), 2)
param$pop_size_se <- round(get_sd(lower = input$pop_size_lower, upper = input$pop_size_upper, coverage = CP), 3)
} # end (if3)
}
param$pop_size_unit <- input$pop_size_unit
})
#################################
## Population growth
##-------------------------------
observe({
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
# Case 1 : Values from expert elicitation
if(input$pop_growth_input_type == "eli_exp"){
if(!(is.null(param$pop_growth_eli_result))){
param$pop_growth_mean <- round(min(1 + param$rMAX_species, round(param$pop_growth_eli_result$mean, 2)), 2)
param$pop_growth_se <- round(param$pop_growth_eli_result$SE, 2)
ready$pop_growth <- TRUE
} else {
ready$pop_growth <- FALSE
}
} else {
# Case 2 : Trend information
if(input$pop_growth_input_type == "trend"){
ready$pop_growth <- TRUE
if(input$pop_trend == "growth") {
if(input$pop_trend_strength == "weak") {
param$pop_growth_mean <- 1.01
} else if(input$pop_trend_strength == "average"){
param$pop_growth_mean <- 1.03
} else {
param$pop_growth_mean <- 1.06
}
} else if(input$pop_trend == "decline"){
if(input$pop_trend_strength == "weak") {
param$pop_growth_mean <- 0.99
} else if(input$pop_trend_strength == "average"){
param$pop_growth_mean <- 0.97
} else {
param$pop_growth_mean <- 0.94
}
} else {
param$pop_growth_mean <- 1
}
param$pop_growth_se <- 0
# Case 3 : Values directly provided (i.e., not from expert elicitation)
} else {
if(input$pop_growth_input_type == "val"){
# Case 2 : Values directly provided as mean & SE
ready$pop_growth <- TRUE
param$pop_growth_mean <- round(min(1 + param$rMAX_species, make_lambda(input$pop_growth_mean)), 3)
param$pop_growth_se <- input$pop_growth_se/100
}else{
# Case 3 : Values directly provided as lower/upper interval
ready$pop_growth <- TRUE
param$pop_growth_mean <- round(min(1 + param$rMAX_species,
round(get_mu(lower = make_lambda(input$pop_growth_lower),
upper = make_lambda(input$pop_growth_upper)), 2)
param$pop_growth_se <- round(get_sd(lower = make_lambda(input$pop_growth_lower),
upper = make_lambda(input$pop_growth_upper), coverage = CP), 3)
}
}
})
#################################
## Carrying capacity
##------------------------------
observe({
if(input$carrying_cap_input_type == "eli_exp"){
if(!(is.null(param$carrying_cap_eli_result))){
param$carrying_capacity <- round(param$carrying_cap_eli_result$mean)
ready$carrying_capacity <- TRUE
} else {
ready$carrying_capacity <- FALSE
}
} else {
if(input$carrying_cap_input_type == "unknown"){
ready$carrying_capacity <- TRUE
param$carrying_capacity <- max(param$pop_size_mean*100, 1e8) # use a very large K
}else{
ready$carrying_capacity <- TRUE
param$carrying_capacity <- input$carrying_capacity
}
}
})
#############################################
##-------------------------------------------
param$survivals <- input$mat_fill_vr[,1]
param$fecundities <- input$mat_fill_vr[,2]
}) # end observeEvent
#####
#############################################
## Calibration of survivals & fecundities
##-------------------------------------------
observeEvent({
input$run
},{
# we also define rMAX here
param$rMAX_species <- rMAX_spp(surv = tail(param$survivals,1), afr = min(which(param$fecundities != 0)))
param$vr_calibrated <- calibrate_params(
inits = init_calib(s = param$survivals, f = param$fecundities, lam0 = param$pop_growth_mean),
f = param$fecundities, s = param$survivals, lam0 = param$pop_growth_mean
)
param$s_calibrated <- head(param$vr_calibrated, length(param$survivals))
param$f_calibrated <- tail(param$vr_calibrated, length(param$fecundities))
})
#####
############################################################
## Observe parameter values to be used in simulations run
##----------------------------------------------------------
observe({
param # required to ensure up-to-date values are run
# simple inputs
param$nsim <- input$nsim
param$fatal_constant <- "h" # input$fatalities_unit
# fixed in global environment (for now)
param$theta = theta
param$time_horzion = time_horzion
param$coeff_var_environ = coeff_var_environ
}) # end observe
#####
#####
##-----------------------------------------------------------------------------------
## RUN SIMULATIONS
##-----------------------------------------------------------------------------------
observeEvent({
input$run
}, {
thierrychambert
committed
if(ready$fatalities & ready$pop_size & ready$pop_growth & ready$carrying_capacity){
withProgress(message = 'Simulation progress', value = 0, {
out$run <- run_simul_shiny(nsim = param$nsim,
cumulated_impacts = param$cumulated_impacts,
fatalities_mean = param$fatalities_mean_nb,
fatalities_se = param$fatalities_se_nb,
thierrychambert
committed
onset_time = param$onset_time,
thierrychambert
committed
pop_size_mean = param$pop_size_mean,
pop_size_se = param$pop_size_se,
pop_size_type = param$pop_size_unit,
thierrychambert
committed
pop_growth_mean = param$pop_growth_mean,
pop_growth_se = param$pop_growth_se,
thierrychambert
committed
survivals = param$s_calibrated,
fecundities = param$f_calibrated,
thierrychambert
committed
carrying_capacity = param$carrying_capacity,
theta = param$theta,
rMAX_species = param$rMAX_species,
thierrychambert
committed
model_demo = NULL,
time_horzion = param$time_horzion,
coeff_var_environ = param$coeff_var_environ,
fatal_constant = param$fatal_constant)
}) # Close withProgress
}else{
out$run <- NULL
out$msg <- "error_not_ready"
}
}) # Close observEvent
#####
#####
##-----------------------------------------------------------------------------------
## OUTPUTS
##-----------------------------------------------------------------------------------
##-------------------------------------------
## Impact text
##-------------------------------------------
## Functions to print the output as text (non cumulated impacts)
print_impact_text <- function(impact, lci, uci){
paste0("Impact : ", round(impact, 2)*100, "%",
"[", round(lci, 2)*100, "% ; ", round(uci, 2)*100, "%]")
} # end function print_impact_text
## Functions to print the output as text (non cumulated impacts)
print_impact_table <- function(res){
n_farm <- (dim(res$indiv_farm$impact)[3]-1)
fil <- paste0(round(t(res$indiv_farm$impact[time_horzion, -2, -1]),2)*100, "%")
matrix(fil,
nrow = n_farm,
dimnames = list(paste("Parc",1:n_farm), c("Impact", "IC (min)", "IC (max)"))
)
} # end function print_impact_table
thierrychambert
committed
print_out <- function(){
thierrychambert
committed
if(!is.null(out$run)) {
# Print the result
if(param$cumulated_impacts){
# cumulated impact ==> Table
print_impact_table(res = get_metrics(N = out$run$N, cumulated_impacts = TRUE))
}else{
# non cumulated impact ==> Text
print_impact_text(impact = get_metrics(N = out$run$N)$scenario$impact[time_horzion, "avg",-1],
lci = get_metrics(N = out$run$N)$scenario$impact[time_horzion, "lci",-1],
uci = get_metrics(N = out$run$N)$scenario$impact[time_horzion, "uci",-1])
}
thierrychambert
committed
} else {
# When run is NULL
if(!is.null(out$msg)){
# Print the error msg, if there is one
if(out$msg == "error_not_ready"){
paste0("Erreur: Vous n'avez pas lancer l'analyse 'valeurs experts'")
}else{
paste0("Some other error occurred")
}
}else{
# When no error msg : nothing happens
} # if "msg"
} # if "run
} # end function print_out
# Display title
output$title_impact_result <- renderText({
if(input$run > 0){
"Rsultat : Impact estim au bout de 30 ans"
}
})
# Display result (text for non cumulated impacts)
output$impact_text <- renderText({
if(input$run == 0){
}else{
if(!param$cumulated_impacts){
print_out()
} else{
NULL
}
# Display result (table for cumulated impacts)
output$impact_table <- renderTable({
if(input$run == 0){
NULL
}else{
if(param$cumulated_impacts){
print_out()
} else{
NULL
}
}
}, rownames = TRUE)
##-------------------------------------------
## Plot Impacts
##-------------------------------------------
## Function to plot the impact
plot_out_impact <- function(){
if(is.null(out$run)) {} else {
plot_impact(N = out$run$N, onset_year = param$onset_year, percent = TRUE,
xlab = "\nAnne", ylab = "Impact relatif (%)\n")
}
output$title_impact_plot <- renderText({
if(input$run > 0){
"Rsultat : Impact relatif au cours du temps"
}
})
output$impact_plot <- renderPlot({
plot_out_impact()
})
##-------------------------------------------
## Plot Demographic Trajectories
##-------------------------------------------
# Function to plot trajectories
plot_out_traj <- function(){
if(is.null(out$run)) {
} else {
plot_traj(N = out$run$N, xlab = "Anne", ylab = "Taille de population (toutes classes d'ges)")}
output$title_traj_plot <- renderText({
if(input$run > 0){
"Graphique : Trajectoire dmographique"
}
})
output$traj_plot <- renderPlot({
plot_out_traj()
})
#####
###################################################################################