Newer
Older
thierrychambert
committed
param$fatalities_se <- c(0, round(get_sd(lower = input$fatalities_lower, upper = input$fatalities_upper, coverage = CP), 3))
} # end (if2)
# Case 2 : Cumulated effects
} else {
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
if(input$analysis_choice == "cumulated"){
ready$fatalities <- TRUE
param$fatalities_mean <- c(0, input$fatalities_mat_cumulated[,1])
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
# Case 3 : Scenarios
}else{
req(input$fatalities_vec_scenario)
vec01 <- as.numeric(unlist(strsplit(input$fatalities_vec_scenario, " ")))
param$fatalities_mean <- c(0, vec01)
param$fatalities_se <- rep(0, length(vec01)+1)
param$onset_time <- NULL
ready$fatalities <- TRUE
}
} # end (if1)
}) # end observe
###~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~###
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
# 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({
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
# 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){
out$analysis_choice <- input$analysis_choice
thierrychambert
committed
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
##-----------------------------------------------------------------------------------
##-------------------------------------------
##-------------------------------------------
print_impact <- function(){
req(out$run)
# cumulated impact
if(param$cumulated_impacts){
res = get_metrics(N = out$run$N, cumulated_impacts = TRUE)
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)"))
)
# Not cumulated impacts
}else{
res = get_metrics(N = out$run$N, cumulated_impacts = FALSE)
n_scen <- (dim(res$scenario$impact)[3]-1)
fil <- paste0(round(t(res$scenario$impact[time_horzion, -2, -1]),2)*100, "%")
matrix(fil,
nrow = n_scen,
dimnames = list(paste("Scenario",1:n_scen), c("Impact", "IC (min)", "IC (max)"))
)
# Display title
output$title_impact_result <- renderText({
if(input$run > 0){
"Rsultat : Impact estim au bout de 30 ans"
}
})
1313
1314
1315
1316
1317
1318
1319
1320
1321
1322
1323
1324
1325
1326
1327
1328
1329
1330
1331
1332
1333
1334
1335
1336
1337
1338
1339
1340
1341
1342
1343
1344
# Display impact result (table)
output$impact_table <- renderTable({
req(input$run)
print_impact()
}, rownames = TRUE)
##-------------------------------------------
## Probability of extinction
##-------------------------------------------
print_PrExt <- function(){
req(out$run)
res = get_metrics(N = out$run$N, cumulated_impacts = FALSE)
n_scen <- dim(res$scenario$impact)[3]
RowNam <- NULL
if(out$analysis_choice == "single_farm") RowNam <- c("Sans parc", "Avec parc")
if(out$analysis_choice == "cumulated") RowNam <- c("Sans parc", "+ Parc 1", paste("... + Parc", (3:n_scen)-1))
if(out$analysis_choice == "multi_scenario") RowNam <- paste("Scenario", (1:n_scen)-1)
fil <- paste0(round(t(res$scenario$Pext),2)*100, "%")
matrix(fil,
nrow = n_scen,
dimnames = list(RowNam, c("Probabilit d'extinction"))
)
} # end function print_PrExt
# Display title
output$title_PrExt_result <- renderText({
if(input$run > 0){
"Rsultat : Probabilit d'extinction 30 ans"
# Display impact result (table)
output$PrExt_table <- renderTable({
req(input$run)
print_PrExt()
}, 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()
})
#####
###################################################################################