This supplementary is available under https://github.com/marruett/Supplementary_Ruett_2020.git and outlines the Decision Analysis approaches used in Model-based evaluation of management options in ornamental plant nurseries(Ruett et al., 2020). Here we provide our input table, code and results to allow for transparent reproduction of the model. To create and organize this supplementary file we used the packages devtools (Wickham, Hester, et al., 2019), knitr (Xie, 2019a), kableExtra (Zhu, 2019) and bookdown (Xie, 2019b). All other applied R packages are cited and listed in the References.

Input table for model-based evaluation of management options in ornamental plant nurseries

The input table Input_file.csv (attached) contains: Variable names, unit, distributions (posnorm = positive normal distribution, const = constant, tnorm_0_1= truncated normal distribution, norm = normal distribution), lower bound, mean, upper bound and definition of variables.

Implementing the model in R

Here we provide our mathematical model that enables simulation of decision outcomes. The calculations in this document are based on random draws from the distributions in the input table Input_file.csv. Explanations for code chunks are not part of the actual simulation. We used the following code descriptions to guide readers through our process of code development.

Load the input table Input_file.csv and make variables.

input_table <- "Input_file.csv"

make_variables <- function(est,n=1)
{ x <- random(rho=est, n=n)
for(i in colnames(x)) assign(i, as.numeric(x[1,i]),envir = .GlobalEnv)}

make_variables(estimate_read_csv(input_table))

We generated a function called Calluna_Simulation to evaluate the decision options in decisionSupport (Luedeling et al., 2019). The following calculations are part of this function (see the entire function below).

Calluna_Simulation <- function(x, varnames){

First we defined the number of plants in the whole production area.

  original_plant_number <- production_area * plants_per_ha

We made a vector called W to define the the risky months of the year (May, June, July and August).

  W <- weather_arguments_for_infection <- c(0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0)

Then we defined a vector called MPD to define the increased monitoring rate (May till June for the monitoring plan).

  MPD <- c(0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0)

We used the vv() (value varier function) function from decisionSupport (Luedeling et al., 2019) to account for the monthly variation in infection risks. Here it is used to generate values out of variables pre-defined upper and lower bounds over time. We use sapply (Team, 2019), which is available in base R to insert percentage values for risky months.

We also used Pipes %>% that allow faster calculation of connected calculations. Pipes are available in the package magrittr (Bache and Wickham, 2014).

In our simulation the general infection risk (risk_per_month) is the same for all management options.

  risk_per_month <- vv(infection_risk, var_CV, 12)* W %>%
  sapply(., function(x) max(c(min(c(1,x)),0)))

Impact of monitoring and direct pesticide applications may reduce the number of symptomatic plants, but risks like disease expansion can increase the number of symptomatic plants again. Therefore, our simulated management options (Normal, Watch, WatchReduce and Reduce) have different impacts on potential outcomes depending on options benefits and costs.

The Baseline scenario (Normal): We simulated fungus occurrence, detection rate, disease expansion and the effect of direct pesticide application over one year of heather production with normal preventive spraying.

  fungus_probability_N <- vv(fungus_probability_N, var_CV, 12)* W %>%
  sapply(., function(x) max(c(min(c(1,x)),0)))
  detection_factor_N <- vv(detection_factor_N, var_CV, 12)* W %>%
  sapply(., function(x) max(c(min(c(1,x)),0)))
  disease_expansion_factor_N <- vv(disease_expansion_factor_N, var_CV, 12)* W %>%
  sapply(., function(x) max(c(min(c(1,x)),0)))
  fungus_fight_effect_N <- vv(fungus_fight_effect_N, var_CV, 12)* W %>%
  sapply(., function(x) max(c(min(c(1,x)),0)))

The reduced preventive spraying scenario (Reduce): Then we simulated fungus occurrence, detection rate, disease expansion and the effect of direct pesticide application over one year of heather production with reduced preventive spraying.

  fungus_probability_R <- vv(fungus_probability_R, var_CV, 12)* W %>%
  sapply(., function(x) max(c(min(c(1,x)),0)))
  detection_factor_R <- vv(detection_factor_R, var_CV, 12)* W %>%
  sapply(., function(x) max(c(min(c(1,x)),0)))
  disease_expansion_factor_R <- vv(disease_expansion_factor_R, var_CV, 12)* W %>%
  sapply(., function(x) max(c(min(c(1,x)),0)))
  fungus_fight_effect_R <- vv(fungus_fight_effect_R, var_CV, 12)* W %>%
  sapply(., function(x) max(c(min(c(1,x)),0)))

The monitoring plan scenario (Watch): We simulated fungus occurrence, detection rate, disease expansion and the effect of direct pesticide application over one year of heather production with normal preventive spraying + monitoring plan.

  fungus_probability_MPN <- vv(fungus_probability_N, var_CV, 12)* W %>%
  sapply(., function(x) max(c(min(c(1,x)),0)))
  detection_factor_MPN <- vv(detection_factor_MP, var_CV, 12)* W %>%
  sapply(., function(x) max(c(min(c(1,x)),0)))
  disease_expansion_factor_MPN <- vv(disease_expansion_factor_N, var_CV, 12)* W %>%
  sapply(., function(x) max(c(min(c(1,x)),0)))
  fungus_fight_effect_MPN <- vv(fungus_fight_effect_MP, var_CV, 12)* W %>%
  sapply(., function(x) max(c(min(c(1,x)),0)))

The monitoring plan with reduced preventive spraying scenario (WatchReduce): Then we simulated fungus occurrence, detection rate, disease expansion and the effect of direct pesticide application over one year of heather production with reduced preventive spraying + monitoring plan.

  fungus_probability_MPR <- vv(fungus_probability_R, var_CV, 12)* W %>%
  sapply(., function(x) max(c(min(c(1,x)),0)))
  detection_factor_MPR <- vv(detection_factor_MP, var_CV, 12)* W %>%
  sapply(., function(x) max(c(min(c(1,x)),0)))
  disease_expansion_factor_MPR <- vv(disease_expansion_factor_R, var_CV, 12)* W %>%
  sapply(., function(x) max(c(min(c(1,x)),0)))
  fungus_fight_effect_MPR <- vv(fungus_fight_effect_MP, var_CV, 12)* W %>%
  sapply(., function(x) max(c(min(c(1,x)),0)))

Here we defined Normal and Reduce preventive spraying costs.

  costs_yearly_prophy_application_N <- round((number_yearly_prophy_application_N),digits=0) * 
  cost_one_prophy_application * production_area
  costs_yearly_prophy_application_R <- round((number_yearly_prophy_application_R),digits=0) * 
  cost_one_prophy_application * production_area

Additionally, we defined preventive spraying costs for Watch and WatchReduce.

  costs_yearly_prophy_application_MPN <- round((number_yearly_prophy_application_N),digits=0) * 
  cost_one_prophy_application * production_area
  costs_yearly_prophy_application_MPR <- round((number_yearly_prophy_application_R),digits=0) * 
  cost_one_prophy_application * production_area

We created the function need.random.integers to split up number of preventive applications per year in preventive applications per risky month.

need.random.integers <- function(a,b,n,k){
  if (n*b<k) stop("k too large")
  if (n*a>k) stop("k too small")
  while(TRUE){
      x <- sample(1:(k - n*a),n-1, replace = TRUE)
      x <- sort(x)
      x <- c(x,k-n*a) - c(0,x)
      if(max(x) <= b-a) return(a+x)}}

Then we defined the number of preventive applications for Normal.

  effect_application_N <- need.random.integers(0, round((number_yearly_prophy_application_N), digits = 0),
                                               4, round((number_yearly_prophy_application_N), digits = 0))* W

We defined the potential of preventive spraying to reduce fungus onset for Normal.

  effect_application_N <- effect_application_N %>%
    replace(., effect_application_N==0, effect_no_prophy_application) %>%
    replace(., effect_application_N==1, effect_one_prophy_application) %>%
    replace(., effect_application_N==2, effect_two_prophy_application) %>%
    replace(., effect_application_N==3, effect_three_prophy_application) %>%
    replace(., effect_application_N==4, effect_four_prophy_application) %>%
    replace(., effect_application_N==5, effect_five_prophy_application) %>%
    replace(., effect_application_N==6, effect_six_prophy_application) %>%
    replace(., effect_application_N==7, effect_seven_prophy_application) %>%
    replace(., effect_application_N>=8, effect_eight_prophy_application) * W

Then we defined the number of preventive applications for Reduce.

  effect_application_R <- need.random.integers(0, round((number_yearly_prophy_application_R), digits = 0),
                                               4, round((number_yearly_prophy_application_R), digits = 0))* W

We defined the potential of preventive spraying to reduce fungus onset for Reduce.

  effect_application_R <- effect_application_R %>%
    replace(., effect_application_R==0, effect_no_prophy_application) %>%
    replace(., effect_application_R==1, effect_one_prophy_application) %>%
    replace(., effect_application_R==2, effect_two_prophy_application) %>%
    replace(., effect_application_R==3, effect_three_prophy_application) %>%
    replace(., effect_application_R==4, effect_four_prophy_application) %>%
    replace(., effect_application_R==5, effect_five_prophy_application) %>%
    replace(., effect_application_R==6, effect_six_prophy_application) %>%
    replace(., effect_application_R==7, effect_seven_prophy_application) %>%
    replace(., effect_application_R>=8, effect_eight_prophy_application) * W

Because the same number of preventive applications were used for Watch and Normal we defined equal preventive spraying effects for these management strategies.

  effect_application_MPN <- effect_application_N

We also defined the same preventive spraying effects for WatchReduce and Reduce.

  effect_application_MPR <- effect_application_R

Here we simulated all management interactions and effects that influence the outcome of Normal heather production. For all management strategies, we adjusted all used vectors so that they had the same length to facilitate management simulation along an one year timescale.

  not_infected_plants_after_prophy_N <- W
  symptomatic_plants_N <- W
  detected_plants_after_monitoring_N <- W
  healthy_plants_after_fungus_fight_N <- W
  getting_again_sick_plants_N <- W

We applied probabilities for infection, effect of prophylactic applications, detection of symptomatic plants and reinfection to simulate impacts of management interventions on the number of detected, not marketable and healthy plants for Normal.

For loops are available in base R but can also be found in purrr (Henry and Wickham, 2019) as part of the tidyverse (Wickham, Averick, et al., 2019) package.

  infected_plant_number_N <- original_plant_number * risk_per_month
  
  for (i in 2:length(infected_plant_number_N)){ 
    infected_plant_number_N[i] <- infected_plant_number_N[i-1] + risk_per_month[i]*
      (original_plant_number - infected_plant_number_N[i-1])}
  for (j in 2:length(not_infected_plants_after_prophy_N)) {
    not_infected_plants_after_prophy_N[j] <- not_infected_plants_after_prophy_N[j-1] + effect_application_N[j]*
      (infected_plant_number_N[j] - not_infected_plants_after_prophy_N[j-1])}
  still_infected_plants_N <- infected_plant_number_N - not_infected_plants_after_prophy_N
  for (k in 2:length(symptomatic_plants_N)) {
    symptomatic_plants_N[k] <- symptomatic_plants_N[k-1] + fungus_probability_N[k]*
      (still_infected_plants_N[k] - symptomatic_plants_N[k-1])}
  for (l in 2:length(detected_plants_after_monitoring_N)) {
    detected_plants_after_monitoring_N[l] <- detected_plants_after_monitoring_N[l-1] + detection_factor_N[l]*
      (symptomatic_plants_N[l] - detected_plants_after_monitoring_N[l-1])}
  symptomatic_plants_after_monitoring_N <- symptomatic_plants_N - detected_plants_after_monitoring_N
  for (m in 2:length(getting_again_sick_plants_N)) {
    getting_again_sick_plants_N[m] <- getting_again_sick_plants_N[m-1] + disease_expansion_factor_N[m]*
      (symptomatic_plants_after_monitoring_N[m] - getting_again_sick_plants_N[m-1])}
  all_symptomatic_plants_N <- symptomatic_plants_after_monitoring_N + getting_again_sick_plants_N
  for (n in 2:length(healthy_plants_after_fungus_fight_N)) {
    healthy_plants_after_fungus_fight_N[n] <- healthy_plants_after_fungus_fight_N[n-1] + fungus_fight_effect_N[n]*
      (all_symptomatic_plants_N[n] - healthy_plants_after_fungus_fight_N[n-1])}

We ran the simulation to generate the number of infected plants, direct plant losses and actually marketable plants for Normal.

  final_fungus_infected_plants_N <- all_symptomatic_plants_N[12] - healthy_plants_after_fungus_fight_N[12]
  direct_plant_losses_N <- detected_plants_after_monitoring_N[12]
  actual_saleable_Callunas_N <- original_plant_number - (final_fungus_infected_plants_N + direct_plant_losses_N)

Here we simulated impacts of management interventions on the number of detected, not marketable and healthy plants for Reduce.

  not_infected_plants_after_prophy_R <- W
  symptomatic_plants_R <- W
  detected_plants_after_monitoring_R <- W
  healthy_plants_after_fungus_fight_R <- W
  getting_again_sick_plants_R <- W
  infected_plant_number_R <- original_plant_number * risk_per_month
  for (i in 2:length(infected_plant_number_R)){
    infected_plant_number_R[i] <- infected_plant_number_R[i-1] + risk_per_month[i]*
      (original_plant_number - infected_plant_number_R[i-1])}
  for (j in 2:length(not_infected_plants_after_prophy_R)) {
    not_infected_plants_after_prophy_R[j] <- not_infected_plants_after_prophy_R[j-1] + effect_application_R[j]*
      (infected_plant_number_R[j] - not_infected_plants_after_prophy_R[j-1])}
  still_infected_plants_R <- infected_plant_number_R - not_infected_plants_after_prophy_R
  for (k in 2:length(symptomatic_plants_R)) {
    symptomatic_plants_R[k] <- symptomatic_plants_R[k-1] + fungus_probability_R[k]*
      (still_infected_plants_R[k] - symptomatic_plants_R[k-1])}
  for (l in 2:length(detected_plants_after_monitoring_R)) {
    detected_plants_after_monitoring_R[l] <- detected_plants_after_monitoring_R[l-1] + detection_factor_R[l]*
      (symptomatic_plants_R[l] - detected_plants_after_monitoring_R[l-1])}
  symptomatic_plants_after_monitoring_R <- symptomatic_plants_R - detected_plants_after_monitoring_R
  for (m in 2:length(getting_again_sick_plants_R)) {
    getting_again_sick_plants_R[m] <- getting_again_sick_plants_R[m-1] + disease_expansion_factor_R[m]*
      (symptomatic_plants_after_monitoring_R[m] - getting_again_sick_plants_R[m-1])}
  all_symptomatic_plants_R <- symptomatic_plants_after_monitoring_R + getting_again_sick_plants_R
  for (n in 2:length(healthy_plants_after_fungus_fight_R)) {
    healthy_plants_after_fungus_fight_R[n] <- healthy_plants_after_fungus_fight_R[n-1] + fungus_fight_effect_R[n]*
      (all_symptomatic_plants_R[n] - healthy_plants_after_fungus_fight_R[n-1])}

Then we ran the simulation to generate the number of infected plants, direct plant losses and actually marketable plants for Reduce.

  final_fungus_infected_plants_R <- all_symptomatic_plants_R[12] - healthy_plants_after_fungus_fight_R[12]
  direct_plant_losses_R <- detected_plants_after_monitoring_R[12]
  actual_saleable_Callunas_R <- original_plant_number - (final_fungus_infected_plants_R + direct_plant_losses_R)

Here we simulated impacts of management interventions on the number of detected, not marketable and healthy plants for Watch.

  not_infected_plants_after_prophy_MPN <- W
  symptomatic_plants_MPN <- W
  detected_plants_after_monitoring_MPN <- MPD
  healthy_plants_after_fungus_fight_MPN <- W
  getting_again_sick_plants_MPN <- W
  infected_plant_number_MPN <- original_plant_number * risk_per_month
  for (i in 2:length(infected_plant_number_MPN)){ 
    infected_plant_number_MPN[i] <- infected_plant_number_MPN[i-1] + risk_per_month[i]*
      (original_plant_number - infected_plant_number_MPN[i-1])}
  for (j in 2:length(not_infected_plants_after_prophy_MPN)) {
    not_infected_plants_after_prophy_MPN[j] <- not_infected_plants_after_prophy_MPN[j-1] + effect_application_MPN[j]*
      (infected_plant_number_MPN[j] - not_infected_plants_after_prophy_MPN[j-1])}
   still_infected_plants_MPN <- infected_plant_number_MPN - not_infected_plants_after_prophy_MPN
  for (k in 2:length(symptomatic_plants_MPN)) {
    symptomatic_plants_MPN[k] <- symptomatic_plants_MPN[k-1] + fungus_probability_MPN[k]*
      (still_infected_plants_MPN[k] - symptomatic_plants_MPN[k-1])}
  for (l in 2:length(detected_plants_after_monitoring_MPN)) {
    detected_plants_after_monitoring_MPN[l] <- detected_plants_after_monitoring_MPN[l-1] + detection_factor_MPN[l]*
      (symptomatic_plants_MPN[l] - detected_plants_after_monitoring_MPN[l-1])}
  symptomatic_plants_after_monitoring_MPN <- symptomatic_plants_MPN - detected_plants_after_monitoring_MPN
  for (m in 2:length(getting_again_sick_plants_MPN)) {
    getting_again_sick_plants_MPN[m] <- getting_again_sick_plants_MPN[m-1] + disease_expansion_factor_MPN[m]*
      (symptomatic_plants_after_monitoring_MPN[m] - getting_again_sick_plants_MPN[m-1])}
  all_symptomatic_plants_MPN <- symptomatic_plants_after_monitoring_MPN + getting_again_sick_plants_MPN
  for (n in 2:length(healthy_plants_after_fungus_fight_MPN)) {
    healthy_plants_after_fungus_fight_MPN[n] <- healthy_plants_after_fungus_fight_MPN[n-1] + fungus_fight_effect_MPN[n]*
      (all_symptomatic_plants_MPN[n] - healthy_plants_after_fungus_fight_MPN[n-1])}

We ran the simulation to generate the final number of infected plants, direct plant losses and actually marketable plants for Watch.

  final_fungus_infected_plants_MPN <- all_symptomatic_plants_MPN[12] - healthy_plants_after_fungus_fight_MPN[12]
  direct_plant_losses_MPN <- detected_plants_after_monitoring_MPN[12]
  actual_saleable_Callunas_MPN <- original_plant_number - (final_fungus_infected_plants_MPN + direct_plant_losses_MPN)

Here we simulated impacts of management interventions on the number of detected, not marketable and healthy plants for WatchReduce.

  not_infected_plants_after_prophy_MPR <- W
  symptomatic_plants_MPR <- W
  detected_plants_after_monitoring_MPR <- MPD
  healthy_plants_after_fungus_fight_MPR <- W
  getting_again_sick_plants_MPR <- W
  infected_plant_number_MPR <- original_plant_number * risk_per_month
  for (i in 2:length(infected_plant_number_MPR)){ 
    infected_plant_number_MPR[i] <- infected_plant_number_MPR[i-1] + risk_per_month[i]*
      (original_plant_number - infected_plant_number_MPR[i-1])}
  for (j in 2:length(not_infected_plants_after_prophy_MPR)) {
    not_infected_plants_after_prophy_MPR[j] <- not_infected_plants_after_prophy_MPR[j-1] + effect_application_MPR[j]*
      (infected_plant_number_MPR[j] - not_infected_plants_after_prophy_MPR[j-1])}
  still_infected_plants_MPR <- infected_plant_number_MPR - not_infected_plants_after_prophy_MPR
  for (k in 2:length(symptomatic_plants_MPR)) {
    symptomatic_plants_MPR[k] <- symptomatic_plants_MPR[k-1] + fungus_probability_MPR[k]*
      (still_infected_plants_MPR[k] - symptomatic_plants_MPR[k-1])}
  for (l in 2:length(detected_plants_after_monitoring_MPR)) {
    detected_plants_after_monitoring_MPR[l] <- detected_plants_after_monitoring_MPR[l-1] + detection_factor_MPR[l]*
      (symptomatic_plants_MPR[l] - detected_plants_after_monitoring_MPR[l-1])}
  symptomatic_plants_after_monitoring_MPR <- symptomatic_plants_MPR - detected_plants_after_monitoring_MPR
  for (m in 2:length(getting_again_sick_plants_MPR)) {
    getting_again_sick_plants_MPR[m] <- getting_again_sick_plants_MPR[m-1] + disease_expansion_factor_MPR[m]*
      (symptomatic_plants_after_monitoring_MPR[m] - getting_again_sick_plants_MPR[m-1])}
  all_symptomatic_plants_MPR <- symptomatic_plants_after_monitoring_MPR + getting_again_sick_plants_MPR
  for (n in 2:length(healthy_plants_after_fungus_fight_MPR)) {
    healthy_plants_after_fungus_fight_MPR[n] <- healthy_plants_after_fungus_fight_MPR[n-1] + fungus_fight_effect_MPR[n]*
      (all_symptomatic_plants_MPR[n] - healthy_plants_after_fungus_fight_MPR[n-1])}

We ran the simulation to generate the final number of infected plants, direct plant losses and actually marketable plants for WatchReduce.

  final_fungus_infected_plants_MPR <- all_symptomatic_plants_MPR[12] - healthy_plants_after_fungus_fight_MPR[12]
  direct_plant_losses_MPR <- detected_plants_after_monitoring_MPR[12]
  actual_saleable_Callunas_MPR <- original_plant_number - (final_fungus_infected_plants_MPR + direct_plant_losses_MPR)

Then we defined the expected value for premium heather plants. The chance_event in the decisionSupport package simulates whether a premium price is achieved based on a certain probability. If arguments indicate that the chance of reaching a higher price is present, a premium price is selected.

value_of_saleable_more_sustainable_Calluna <- value_of_saleable_Calluna +
    chance_event(chance_higher_price_sustainable, price_premium_sustainable, 0)

We computed the value of marketable heather for Normal.

  value_saleable_plants_N <- actual_saleable_Callunas_N * value_of_saleable_Calluna
  value_of_sorted_out_plants_N <- direct_plant_losses_N * value_sorted_out_Calluna
  value_of_not_saleable_plants_N <- final_fungus_infected_plants_N * value_not_saleable_Calluna

We computed the value of marketable heather for Reduce.

  value_saleable_plants_R <- actual_saleable_Callunas_R * value_of_saleable_more_sustainable_Calluna
  value_of_sorted_out_plants_R <- direct_plant_losses_R * value_sorted_out_Calluna
  value_of_not_saleable_plants_R <- final_fungus_infected_plants_R * value_not_saleable_Calluna

We computed the value of marketable heather for Watch.

  value_saleable_plants_MPN <- actual_saleable_Callunas_MPN * value_of_saleable_Calluna
  value_of_sorted_out_plants_MPN <- direct_plant_losses_MPN * value_sorted_out_Calluna
  value_of_not_saleable_plants_MPN <- final_fungus_infected_plants_MPN * value_not_saleable_Calluna

We computed the value of marketable heather for WatchReduce.

  value_saleable_plants_MPR <- actual_saleable_Callunas_MPR * value_of_saleable_more_sustainable_Calluna
  value_of_sorted_out_plants_MPR <- direct_plant_losses_MPR * value_sorted_out_Calluna
  value_of_not_saleable_plants_MPR <- final_fungus_infected_plants_MPR * value_not_saleable_Calluna

Then we added up the yearly costs for Normal.

  costs_per_year_N <- (value_of_one_new_plant * original_plant_number) +
    value_of_sorted_out_plants_N +
    value_of_not_saleable_plants_N +
    costs_yearly_prophy_application_N +
    costs_normal_fertilizer * production_area +
   (costs_monitoring_per_ha_month * sum(W) * production_area) +
    costs_staff +
    yearly_costs_of_direct_fungus_fight_N * production_area

Here we added up the benefits for Normal.

  benefits_per_year_N <- value_saleable_plants_N

We simulated the general staff costs.

  costs_staff <- ifelse(production_area > threshold_big_area_more_staff, costs_staff + 
                          costs_more_staff, costs_staff)

And we added up yearly costs for Reduce.

  costs_per_year_R <- (value_of_one_new_plant * original_plant_number) +
    value_of_sorted_out_plants_R +
    value_of_not_saleable_plants_R +
    costs_yearly_prophy_application_R +
   (costs_normal_fertilizer + fertilizer_adjustment) * production_area +
   (costs_monitoring_per_ha_month * sum(W) * production_area) +
    additional_costs_more_monitoring_per_ha * production_area +
    costs_staff +
    yearly_costs_of_direct_fungus_fight_R * production_area

We added up the benefits for Reduce.

  benefits_per_year_R <- value_saleable_plants_R

Here we added up yearly costs for Watch.

  costs_per_year_MPN <- (value_of_one_new_plant * original_plant_number) +
    value_of_sorted_out_plants_MPN +
    value_of_not_saleable_plants_MPN +
   (direct_plant_losses_MPN * amount_of_samples_MP * costs_per_sample_MP) +
    costs_yearly_prophy_application_MPN +
    costs_normal_fertilizer * production_area +
   (costs_monitoring_plan_per_ha * production_area) +
    costs_staff +
    yearly_costs_of_direct_fungus_fight_N * production_area

Here we added up the benefits for Watch.

  benefits_per_year_MPN <- value_saleable_plants_MPN + (savings_due_to_MP * production_area)

And here we added up yearly costs for WatchReduce.

  costs_per_year_MPR <- (value_of_one_new_plant * original_plant_number) +
    value_of_sorted_out_plants_MPR +
    value_of_not_saleable_plants_MPR +
   (direct_plant_losses_MPR * amount_of_samples_MP * costs_per_sample_MP) +
    costs_yearly_prophy_application_MPR +
   (costs_normal_fertilizer + fertilizer_adjustment) * production_area +
   (costs_monitoring_plan_per_ha * production_area) +
    additional_costs_more_monitoring_per_ha * production_area +
    costs_staff +
    yearly_costs_of_direct_fungus_fight_R * production_area

Then we added up the benefits for WatchReduce.

  benefits_per_year_MPR <- value_saleable_plants_MPR + (savings_due_to_MP * production_area)

For all management strategies we subtracted the costs from the benefits.

  cashflow_N <- benefits_per_year_N - costs_per_year_N
  cashflow_R <- benefits_per_year_R - costs_per_year_R
  cashflow_MPN <- benefits_per_year_MPN - costs_per_year_MPN
  cashflow_MPR <- benefits_per_year_MPR - costs_per_year_MPR

To consider the relative merits of our stated management strategies, we compared each scenario and its outcomes with the Normal scenario.

Our model outcomes are normalized by production area to ensure comparability of Monte Carlo simulation runs. We did not analyze comp_MPR_MPN, because we chose to compare all management options with the current established production practices. Growers usually decide on adopting innovative practices based on their relative merits compared to the practices that are currently used, not relative to a hypothetical baseline situation that does not exist at present. Furthermore, results of the comp_MPR_MPN comparison indicated a very low likelihood of positive outcomes. We still show the code for this comparison here, since we did calculate it, but we do not report on the results.

  comp_R_N<-(cashflow_R-cashflow_N)/production_area
  comp_MPR_N<-(cashflow_MPR-cashflow_N)/production_area
  comp_MPN_N<-(cashflow_MPN-cashflow_N)/production_area
  comp_MPR_MPN<-(cashflow_MPR-cashflow_MPN)/production_area

At the end of the Calluna_Simulation function we listed all calculated outcomes.

  return(list(comp_R_N=comp_R_N,
              comp_MPR_N=comp_MPR_N,
              comp_MPN_N=comp_MPN_N,
              comp_MPR_MPN=comp_MPR_MPN,
              cashflow_N=cashflow_N,
              cashflow_R=cashflow_R,
              cashflow_MPR=cashflow_MPR,
              cashflow_MPN=cashflow_MPN))}

Here we defined input table, legend file, results folder and figures folder.

input_table <- "Input_file.csv"
legend_file <- "Legend_file.csv"
results_folder <- "Results_low_prophy"
figures_folder <- "EVPI_tables_low_prophy"

The decisionSupport function facilitated the start of the Monte Carlo simulation.

decisionSupport(input_table,
                results_folder,
                write_table = TRUE, Calluna_Simulation, 10000,
                functionSyntax = "plainNames")

We produced all outcome variables with the decisionSupport function. Files which start with cashflow_ refer to management options while comp_ files refer to simulated decisions in comparison to current established production practices.

outvars <- list("comp_R_N","comp_MPR_N","comp_MPN_N","comp_MPR_MPN","cashflow_N","cashflow_R","cashflow_MPR","cashflow_MPN")
legend_table <- read.csv(legend_file)

labels <- list("comp_R_N","comp_MPR_N","comp_MPN_N","comp_MPR_MPN","cashflow_N","cashflow_R","cashflow_MPR","cashflow_MPN")

The multi_EVPI function is also part of the decisionSupport package and enables EVPI analysis by pointing out variables uncertainties.

MC_file <- read.csv(paste(results_folder,"/mcSimulationResults.csv",sep = ""))


mc_EVPI <- MC_file[,-grep("cashflow",colnames(MC_file))]
mc_EVPI <- MC_file
multi_EVPI(mc = mc_EVPI,"comp_R_N",write_table = TRUE,outfolder = figures_folder)

Outcome distributions

Outcome distributions for management strategies

The following plots are generated from the mcSimulationResults.csv from the decisionSupport function. We plotted probability distributions and merged them to illustrate the frequency of Monte Carlo simulation outcomes. Additionally, we integrated boxplots to compare Partial Farm Budget of Watch, WatchReduce and Reduce and Normal.

For data visualization we applied the ggplot2 (Wickham, 2016), ggstance (Henry et al., 2019) and dplyr (Wickham, François, et al., 2019) packages. We used the function ggarrange from the ggpubr package (Kassambara, 2019) to merge figures.

Probability density distributions (scaled density among 10,000 runs of a Monte Carlo simulation) of Partial Farm Budgets (Euros) for four heather production management strategies on an area of one hectare. Watch = Monitoring Plan and normal prophylactic pesticide application. WatchReduce = Monitoring Plan and reduced prophylactic pesticide application. Reduce = Reduced prophylactic pesticide application. Normal = Normal prophylactic pesticide application.

Outcome distributions for management decisions

Comparison of Partial Farm Budget of WatchMore, WatchMoreSprayLess and SprayLess.

Probability density distributions (scaled density among 10,000 runs of a Monte Carlo simulation) of Partial Farm Budgets (Euros) for a heather production area of one hectare for three decision options for heather production. WatchMore = Improved monitoring combined with normal prophylactic application compared to current standard practices. WatchMoreSprayLess = Monitoring combined with reduced pesticide application compared to standard practices. SprayLess = Reduced prophylactic pesticide application compared to standard practice, in a system without improved monitoring.

EVPI and VIP analysis of management decisions

Comparison of Information Value and Variable Importance for WatchMore, WatchMoreSprayLess and SprayLess.

Here we present Variable Importance in the Projection (VIP) and Information Value (Expected Value of Perfect Information; EVPI) for the implementation of WatchMore, WatchMoreSprayLess and SprayLess in comparison to established heather production. VIP and EVPI specify the influence of each variable on the projected economic value of adopting each strategy (Variable Importance) as well as on whether the simulation makes adoption or non-adoption appear economically preferable. For Variable Importance, we show the top eight variables with the respective EVPI. In the Variable Importance plot, variables that positively affect projected Partial Farm Budget are shown in green, those with negative effects in red.

Acknowledgements

This research was financially supported by Stiftung Zukunft NRW within the research project inruga (Innovationen fuer NRW zur Steigerung der Ressourceneffizienz und Umweltvertraeglichkeit im Gartenbau Entscheidungshilfen im Zierpflanzenbau). We thank Uwe Rascher and Laura Junker (Forschungszentrum Juelich, IBG-II), Andrew Gallik, Elisabeth Goette, Monika Heupel, Rainer Peters, Michael Stuch, Peter Tiede-Arlt and Rainer Wilke (Landwirtschaftskammer Nordrhein-Westfahlen), Hannah Jaenicke and Miriam Robertz (University of Bonn) and Martin Balmer (DLR) for their workshop participation, commitment and advice in this study. We would like to acknowledge the heather farmers Gerd Canders, Tom Canders, Matthias Kueppers and Verena Kueppers for their participation and contributions throughout the research process.

References

Bache, S.M., Wickham, H., 2014. Magrittr: A Forward-Pipe Operator for R.

Henry, L., Wickham, H., 2019. Purrr: Functional Programming Tools.

Henry, L., Wickham, H., Chang, W., 2019. Ggstance: Horizontal ’ggplot2’ Components.

Kassambara, A., 2019. Ggpubr: ’Ggplot2’ Based Publication Ready Plots.

Luedeling, E., Goehring, L., Schiffers, K., 2019. decisionSupport: Quantitative Support of Decision Making under Uncertainty, decisionSupport Version 1.105.2.

Ruett, M., Whitney, C., Luedeling, E., 2020. Model-based evaluation of management options in ornamental plant nurseries. Journal of Cleaner Production 271, 122653. https://doi.org/10.1016/j.jclepro.2020.122653

Team, R.D.C., 2019. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria.

Wickham, H., 2016. Ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York.

Wickham, H., Averick, M., Bryan, J., Chang, W., D’Agostino McGowan, L., François, R., Grolemund, G.G., Hayes, A., Henry, L., Hester, J., Kuhn, M., Lin Pedersen, T., Miller, E., Bache, S.M., Müller, K., Ooms, J., Robinson, D., Paige Seidel, D., Spinu, V., Takahashi, K., Vaughan, D., Wilke, K., Woo, K., Yutani, H., 2019. Welcome to the {tidyverse}. J. Open Source Softw. 4, 1686. https://doi.org/10.21105/joss.01686

Wickham, H., François, R., Henry, L., Müller, K., 2019. Dplyr: A Grammar of Data Manipulation.

Wickham, H., Hester, J., Chang, W., 2019. Devtools: Tools to Make Developing R Packages Easier.

Xie, Y., 2019a. Knitr: A General-Purpose Package for Dynamic Report Generation in R.

Xie, Y., 2019b. Bookdown: Authoring Books and Technical Documents with R Markdown.

Zhu, H., 2019. kableExtra: Construct Complex Table with ’kable’ and Pipe Syntax.