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
##-----------------------------------------------------------------------------------
1279
1280
1281
1282
1283
1284
1285
1286
1287
1288
1289
1290
1291
1292
1293
1294
1295
1296
1297
1298
1299
1300
1301
1302
1303
1304
1305
1306
1307
1308
#######################################################################
## Impact : individual farms (for "cumulated impact" analysis only)
##---------------------------------------------------------------------
print_indiv_impact <- function(){
req(out$run)
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)"))
)
} # end function print_impact
# Display title
output$title_indiv_impact_result <- renderText({
req(input$run > 0, out$analysis_choice == "cumulated")
"Rsultat : Impact de chaque parc olien, estim au bout de 30 ans"
})
# Display impact result (table)
output$indiv_impact_table <- renderTable({
req(input$run & out$analysis_choice == "cumulated")
print_indiv_impact()
}, rownames = TRUE)
##################################################
## Impact : GLOBAL (for all types of analysis)
##------------------------------------------------
print_impact <- function(){
req(out$run)
res = get_metrics(N = out$run$N, cumulated_impacts = FALSE)
n_scen <- (dim(res$scenario$impact)[3]-1)
RowNam <- NULL
if(out$analysis_choice == "single_farm") RowNam <- c("Parc 1")
if(out$analysis_choice == "cumulated") RowNam <- c("Parc 1", paste("... + Parc", (2:n_scen)))
if(out$analysis_choice == "multi_scenario") RowNam <- paste("Scenario", (1:n_scen))
fil <- paste0(round(t(res$scenario$impact[time_horzion, -2, -1]),2)*100, "%")
matrix(fil,
nrow = n_scen,
dimnames = list(RowNam, c("Impact", "IC (min)", "IC (max)"))
)
# Display title
output$title_impact_result <- renderText({
req(input$run)
"Rsultat : Impact global estim au bout de 30 ans"
})
# Display impact result (table)
output$impact_table <- renderTable({
req(input$run)
print_impact()
}, rownames = TRUE)
#############################################
1340
1341
1342
1343
1344
1345
1346
1347
1348
1349
1350
1351
1352
1353
1354
1355
1356
1357
1358
1359
1360
## 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({
req(input$run)
"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(){
n_scen <- dim(out$run$N)[3]
Legend <- NULL
if(out$analysis_choice == "single_farm") Legend <- c("Sans parc", "Avec parc")
if(out$analysis_choice == "cumulated") Legend <- c("Sans parc", "+ Parc 1", paste("... + Parc", (3:n_scen)-1))
if(out$analysis_choice == "multi_scenario") Legend <- paste("Scenario", (1:n_scen)-1)
plot_impact(N = out$run$N, onset_year = param$onset_year, percent = TRUE,
xlab = "\nAnne", ylab = "Impact relatif (%)\n", Legend = Legend)
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 {
n_scen <- dim(out$run$N)[3]
Legend <- NULL
if(out$analysis_choice == "single_farm") Legend <- c("Sans parc", "Avec parc")
if(out$analysis_choice == "cumulated") Legend <- c("Sans parc", "+ Parc 1", paste("... + Parc", (3:n_scen)-1))
if(out$analysis_choice == "multi_scenario") Legend <- paste("Scenario", (1:n_scen)-1)
plot_traj(N = out$run$N, onset_year = param$onset_year,
xlab = "\nAnne", ylab = "Taille de population\n", Legend = Legend)}
output$title_traj_plot <- renderText({
if(input$run > 0){
"Graphique : Trajectoire dmographique"
}
})
output$traj_plot <- renderPlot({
plot_out_traj()
})
#####
###################################################################################