install.packages("mcmodule", repos = "https://nataliaciria.com/r", type = "source")mcmodule: An R Package for Multi-Pathway Monte Carlo Risk Assessment
R!SK 2026 — 19 February 2026, 10:15 AM Eastern
Package Website
https://nataliaciria.com/mcmodule
CRAN
https://cran.r-project.org/package=mcmodule
Presentation materials
Code and data: https://nataliaciria.com/talks/risk2026/risk2026-mcmodule-materials.zip
⚠️ Install mcmodule (for R!SK 2026)
Normally, mcmodule is installed from CRAN. This tutorial, however, was built using the development version, which was due to be released this week as mcmodule v1.2.0.
At the moment I cannot submit to CRAN because my GitHub account was suspended on 16 February 2026. GitHub has not provided a clear reason, but it appears to be related to academic affiliation re-verification, similar to this case: My Account Was Suspended Unexpectedly · Discussion #171069.
To run this example on 19 February 2026, you will need to install a development build of mcmodule from my personal repository.
If installing a package from a personal repository feels risky, that is completely fair. If you would rather wait for an official release, you can follow me on LinkedIn or subscribe to my blog. I will share an update as soon as my GitHub account is restored and v1.2.0 is available on CRAN.
You can also directly download the tar.gz file.
Intro
Quantifying risk through probability steps is a clean and simple way to estimate the likelihood of an unwanted event. With mc2d, adding a stochastic dimension to probability steps is straightforward and allows us to account for uncertainty.
library(mc2d)
# Probability an animal is infected
inf_animal <- mcstoc(runif, min = 0.002, max = 0.08)
# Probability an infected animal is tested
test_animal <- mcstoc(rpert, min = 0.4, mode = 0.5, max = 0.6)
test_inf_animal <- inf_animal * test_animal
# Probability an animal infected animal is not detected
test_se <- mcstoc(runif, min = 0.8, max = 0.95)
not_detected <- test_inf_animal * (1 - test_se)However, while mc2d is excellent for single-pathway models, it becomes very complex to implement with multiple pathways, granularities, and scenarios. Real-life risk pathways are often messy, interconnected, and long, involving many variables and mismatched pathway components. The mcmodule package provides a framework to handle complex risk pathways in a structured way, and to compare scenarios and perform model diagnostics.
Core components
This package is structured around an S3 object class called mcmodule. These modules are created based on four core components:
| Component | Purpose |
|---|---|
| Data | Table(s) with one row per variate and columns for input distribution parameters |
| Data keys | Unique identifiers of rows/variates |
| mctable | Table with node specifications (distributions, transformations) |
| Expressions | Mathematical risk calculations |
This allows to:
- Automatically generate nodes from data
- Identify variates by keys (country, species, scenario, etc.)
- Align operations automatically
- Combine pathways safely
Case study
This case study illustrates how mcmodule supports a multi-pathway import risk question where inputs differ by commodity and origin.
Question: What is the probability of introducing sheep/goat pox to Spain through imports?
We treat the event “introduction” as occurring if at least one imported unit (animal or product lot) is infected and results in effective exposure in Spain.
Pathways:
Live animals (sheep and goat)
Animal products with contamination risk (raw wool, hides)
Each pathway is represented as a sequence of probability steps. It will be computed with one variate per origin country and commodity, then aggregated across origins and products.
Data:
Trade volumes by commodity and country (COMTRADE + Spanish Cámara de Comercio)
Disease status by country (WAHIS)
Animal holdings and heads census by country (Eurostat + FAO)
Set-up
# Load packages
library(mcmodule) # version 1.2.0
library(dplyr)
library(tidyr)
library(ggplot2)
ndvar(1001) # Uncertainty simulations[1] 1001Data
imports_by_country_commodity <- read.csv("data/imports_by_country_commodity.csv", stringsAsFactors = TRUE)
disease_by_country <- read.csv("data/disease_by_country.csv", stringsAsFactors = TRUE)
census_by_country <- read.csv("data/census_by_country.csv", stringsAsFactors = TRUE)When input parameters exist at different aggregation levels, it’s best to load them in separate tables and aggregate them as needed for each pathway. Input data will be defined by different keys and related through one or more shared keys (see Figure 1). This approach helps trace input data resolution and saves memory as pathways that don’t require highly detailed parameters will have fewer variables.
Data keys
These data keys identified in Figure 1 must be declared to be used in mcmodule, using a list with the name of each data set, its columns, and which of them are keys. Use set_data_keys() to set the keys in the environment, avoiding the need to specify them in mcmodule functions.
sgp_data_keys <- list(
imports_by_country_commodity = list(
cols = names(imports_by_country_commodity),
keys = c("country_code", "species", "commodity")
),
disease_by_country = list(
cols = names(disease_by_country),
keys = "country_code"
),
census_by_country = list(
cols = names(census_by_country),
keys = c("country_code", "species")
)
)
set_data_keys(sgp_data_keys) # Set in envir Nodes definition: mctable
To avoid manually building hundreds of nodes, mcmodule creates them automatically based on the information in the mctable. It creates the nodes listed in mcnode using the information in the other columns.
If a probability distribution is specified in mc_func, it will look for the columns with the parameters needed to build those distributions with mcstoc(). For example, if import_qty is a pert distribution, it will look for imports_qty_min, imports_qty_mode, and imports_qty_max columns in the data and build the node with the number of variates the data has. It no probability distribution is specified, it will create mcdata() nodes with the values of columns that don’t have a probability function specified.
The column from_variable allows to use alternative column names present in the data as inputs for that mcnode.
The transformation column lets you convert non-numeric columns to mcnode, for example, our input data has a column that specifies whether an import is live animals or products. We’ll transform this to a boolean so we can use this categorical information in the model.
We usually create this table outside R, in a CSV or similar, so it is easy to modify and can be used as documentation.
| mcnode | description | mc_func | from_variable | transformation |
|---|---|---|---|---|
| import_qty | Annual imports quantity (units for live animals, kilograms for animal products) | rpert | NA | NA |
| import_operations | Annual imports operations | rpert | NA | NA |
| n_years_outbreaks | Number of years with new outbreaks | NA | NA | NA |
| n_years_info | Number of years with reported health status | NA | NA | NA |
| census_heads | Number of animal heads in the country | NA | NA | NA |
| census_holdings | Number of animal holdings in the country | NA | NA | NA |
| present | Disease was present in the country the last year reported | NA | NA | NA |
| border_present | TRUE if disease is present in any neighbouring country | NA | NA | NA |
| n_border_present | Number of neighbouring countries with disease present | NA | NA | NA |
| n_border_years_info | Total years of health status information for neighbouring countries | NA | NA | NA |
| n_border_years_outbreaks | Total years of outbreaks in neighbouring countries | NA | NA | NA |
| live_animals | Live animals imports (TRUE if live animals, FALSE if animal products) | NA | pathway | value==‘live_animals’ |
| p_herds_surveillance | Proportion of herds visited for surveillance at origin | runif | NA | NA |
| p_animals_inspected | Proportion of animals inspected per herd at origin | runif | NA | NA |
| clinical_se_origin | Sensitivity of clinical inspection at origin | runif | NA | NA |
| clinical_se_dest | Sensitivity of clinical inspection at destination | runif | NA | NA |
sgp_mctable <- read.csv(
"data/sgp_mctable.csv",
stringsAsFactors = TRUE
)
set_mctable(sgp_mctable) # Set in envir Expressions
Now we finally get to risk calculations! Here we have a model with some mathematical expressions. Note that it is quoted, which means it is unevaluated R code until we call for its evaluation. Input mcnodes will be created from data and mctable specifications. When they are evaluated, they will keep their metadata with them (their keys, the table they come from, etc.). Output mcnodes metadata will record their inputs and the combined keys of their inputs. Nodes can also be created inline, but in that case no variates must be specified, as they will be automatically assigned by data rows.
# Probability that the origin country is infected with sheep/goat pox
inf_origin_exp <- quote({
# Proportion of years with documented outbreaks
p_years_outbreaks <- n_years_outbreaks / n_years_info
# Number of herds visited for surveillance
n_herd_surveillance <- p_herds_surveillance * census_holdings
# Average herd size
herd_size <- census_heads / census_holdings
# Number of animals inspected
n_animals_inspected <- p_animals_inspected * herd_size
# Surveillance programme design prevalences
herd_design_prev <- 0.1 # Herd-level target
system_design_prev <- 0.001 # System-level target
# Herd-level surveillance sensitivity
herd_se <- 1 -
(1 -
(n_animals_inspected *
clinical_se_origin /
herd_size))^(herd_design_prev *
herd_size)
# System-level surveillance sensitivity
system_se <- 1 - (1 - system_design_prev * herd_se)^n_herd_surveillance
# Probability of disease outbreak occurring in origin country
# Only applies to animal products as live animals are only imported from disease-free regions
p_outbreak <- (1 - live_animals) * present * p_years_outbreaks
# Outbreak in infected countries that escape protection zones
# Time window from outbreak occurrence to detection, days (sources: ANSES, EFSA 2013)
time_det <- mcstoc(rpert, min = 15, mode = 28, max = 60)
# Probability that at least one new outbreak occurs AND an animal is moved
# during the detection window before outbreak is detected
p_new_outbreak <- 1 - exp(-p_years_outbreaks / 365 * time_det)
# Transmission beyond protection zone (source: EFSA 2022)
escape <- mcstoc(rpert, min = 0.01, mode = 0.023, max = 0.055)
# Probability of escaped outbreak in disease-free area within infected country
p_scaped_new_outbreak <- present * p_new_outbreak * escape
# Introduction to disease-free country from infected neighbors
# Proportion of years with documented new outbreaks in neighboring countries
p_border_new_outbreak <- mcnode_na_rm(
n_border_years_outbreaks / n_border_years_info
)
# Probability of undetected introduction escaping surveillance
p_undetected_introduction <- border_present *
p_border_new_outbreak *
(1 - system_se)
# Combined probability that origin is infected
# - p_outbreak: region currently infected
# - p_scaped_new_outbreak: new outbreak escaped protection area in infected countries
# - p_undetected_introduction: introduction from border that is not detected by surveillance
p_inf_origin <- 1 -
(1 - p_outbreak) *
(1 - p_scaped_new_outbreak) *
(1 - p_undetected_introduction)
})Combine expressions
One mcmodule can evaluate multiple expressions, and expressions can use mcnodes from previous expressions. If expressions are organized in lists (recommended), the name of each list element will be used in the mcnode metadata.
# Probability of disease introduction through live animals
# Given that the origin is infected
animal_entry_exp <- quote({
# Within-herd prevalence in infected herds (source: OIE Greece data 2005-2024)
w_prev <- mcstoc(rpert, min = 0.002, mode = 0.08, max = 0.32)
# Number of infected herds before detection
n_herds_pos <- mcstoc(rpert, min = 1, mode = 3, max = 9)
# Number of infected animals before detection
n_animals_pos <- n_herds_pos * w_prev * (census_heads / census_holdings)
# Average number of live animals per import operation
import_size <- import_qty / import_operations
# Probability an animal is infected
p_inf_animal <- n_animals_pos / census_heads
# Probability of survival of infected animals
p_animal_survival <- mcstoc(rpert, min = 0, mode = 0.05, max = 0.80)
# Probability that an infected animal is contagious
p_contag_animal <- mcstoc(rpert, min = 0.002, mode = 0.08, max = 0.32)
# Probability of successful disease introduction through live animals
p_entry_animal <- p_inf_animal *
(1 - p_animal_survival) *
(1 - clinical_se_dest) *
p_contag_animal
})
live_animals_exp <- list(
animal_inf_origin = inf_origin_exp,
animal_entry = animal_entry_exp
)Scenarios
Scenarios are a special type of key that lets you create what-if scenarios. The baseline scenario joins all “real” input datasets and must always be labeled “0” (this is crucial because mcmodule treats the baseline scenario differently in certain operations). You can copy this baseline and modify it to create alternative scenarios. For instance: “What would happen to import risk if clinical sensitivity was higher at destination?”. You can give these alternative scenarios any name you like. Then, bind the rows for all current and alternative scenarios and use it as input data to create the mcmodule.
# Create baseline scenario with current surveillance parameters
baseline_data <- imports_by_country_commodity %>%
left_join(disease_by_country, by = "country_code") %>%
left_join(census_by_country, by = c("country_code", "species")) %>%
filter(!is.na(census_heads)) %>%
arrange(n_years_outbreaks) %>%
mutate(
scenario_id = "0",
p_herds_surveillance_min = 0.3,
p_herds_surveillance_max = 0.5,
p_animals_inspected_min = 0.25,
p_animals_inspected_max = 0.35,
clinical_se_origin_min = 0.55,
clinical_se_origin_max = 0.65,
clinical_se_dest_min = 0.55,
clinical_se_dest_max = 0.70
)
# Enhanced surveillance in destination scenario
enhanced_dest_data <- baseline_data %>%
mutate(
scenario_id = "Destination surveillance",
clinical_se_dest_min = 0.85,
clinical_se_dest_max = 0.95
)
# Enhanced surveillance in origin scenario
enhanced_origin_data <- baseline_data %>%
mutate(
scenario_id = "Origin surveillance",
p_herds_surveillance_min = 0.6,
p_herds_surveillance_max = 0.8,
p_animals_inspected_min = 0.60,
p_animals_inspected_max = 0.70,
clinical_se_origin_min = 0.75,
clinical_se_origin_max = 0.85
)
# Combine scenarios
imports_data <- bind_rows(
baseline_data,
enhanced_dest_data,
enhanced_origin_data
)Create an mcmodule
We’ve filtered the data just the inputs for this pathway, in this case, live animals. We already have our list of expressions, and now all we need to do is evaluate them eval_module(). The mctable and data_keys are already in the environment, so we don’t need to specify them.
live_animals_data <- imports_data %>%
filter(pathway == "live_animals")
# Evaluate live animals model
live_animals <- eval_module(
exp = live_animals_exp,
data = live_animals_data
)This creates an mcmodule class object, a special type of list that contains the input data (which serves as an “index” of keys and variates), the expressions, and a node list containing all the evaluated nodes and their metadata.
For each node, we can access its type, the description from the mctable or the expression that created it, its input columns and dataset, its keys, and the mcnode object itself. This metadata is used in many ways throughout mcmodule, but the most crucial function is associating keys per row to each variate, allowing modules to be manipulated in very flexible ways later on.
Trial totals
mcmodule includes functions to calculate the probability of something happening in at least one occasion after a number of repetitions. We call each of these repetitions a “trial”. These calculations can sometimes be simple and built directly into expressions. However, if, for example, we calculate all expressions at the animal level, it’s very useful to have a function that determines the probability at animal batch level, at any step of the expression (whether at least one animal is infected, not detected, contagious, etc.).
Both the input and ouput of this function are mcmodules. The function automatically creates an “total” mcnode with metadata and inserts it into the mcmodule.
# Calculate pathway totals: all animals in one import operation
live_animals <- trial_totals(
mcmodule = live_animals,
mc_names = "p_entry_animal",
trials_n = "import_size"
)We can also calculate hierarchical trial totals (see Figure 3). These work in two steps: first, select subsets from a larger set (for example, herds from a country’s total sheep population). Then, if the subset is positive, calculate the risk that individual trials within it (in this case, individual sheep) are positive.
# Calculate pathway totals: all animals and import operations
live_animals <- trial_totals(
mcmodule = live_animals,
mc_names = "p_entry_animal",
trials_n = "import_size",
subsets_p = "p_inf_origin",
subsets_n = "import_operations"
)For more details, see the multivariate operations vignette, which explains how to produce aggregated hierarchical trials and other types of totals.
Combine modules
Let’s create and evaluate another mcmodule for animal products imports. We add a new expression to calculate the probability of disease introduction through animal products, and reuse the expression for the probability that the origin is infected. Then we will evaluate it and calculate the totals.
# Prepare data for animal products pathway with scenario variation
animal_products_data <- imports_data %>%
filter(pathway == "products")
# Probability of disease introduction through animal products
# Given that the origin is infected
products_entry_exp <- quote({
# Within-herd prevalence in infected herds
w_prev <- mcstoc(rpert, min = 0.002, mode = 0.08, max = 0.32)
# Number of infected herds and animals in origin
n_herds_pos <- mcstoc(rpert, min = 1, mode = 3, max = 9)
herd_size <- census_heads / census_holdings
n_animals_pos <- n_herds_pos * w_prev * herd_size
# Probability an animal is infected
p_inf_animal <- n_animals_pos / census_heads
# Average quantity per animal product (e.g., kg of wool per animal)
qty_per_animal <- mcstoc(func = rnorm, mean = 5, sd = 1)
# Average number of animals whose products are imported per operation
import_size <- import_qty / qty_per_animal / import_operations
# Probability that product is contaminated (from infected animal)
p_contam_product <- mcstoc(rpert, min = 0.002, mode = 0.08, max = 0.32)
# Probability that contaminated product remains contagious (infectious)
p_contag_product <- mcstoc(rpert, min = 0.00032, mode = 0.0016, max = 0.008)
# Probability that infected product is eliminated by visual screening
p_elimin <- mcstoc(rpert, min = 0, mode = 0.5, max = 0.975)
# Probability that virus survives transport and storage
p_virus_survival <- mcstoc(rpert, min = 0.04, mode = 0.5, max = 1)
# Probability of animal survival (same as live animals)
p_animal_survival <- mcstoc(rpert, min = 0, mode = 0.05, max = 0.80)
# Probability of product introducing contamination at destination
p_entry_product <- p_inf_animal *
p_contam_product *
(1 - p_animal_survival) *
(1 - clinical_se_origin) *
p_contag_product *
(1 - p_elimin) *
p_virus_survival
# Probability of exposure to susceptible animals after entry
p_exposure <- mcstoc(
rpert,
min = 1e-08,
mode = 1e-05,
max = 0.0001
)
# Probability that product is still contaminated at exposure
p_contam_exposure <- mcstoc(runif, min = 0.01, max = 0.1)
# Probability of indirect transmission (e.g., through contaminated environment)
p_indir_trans <- mcstoc(runif, min = 0.01, max = 0.1)
# Probability of infecting animals through contaminated products
p_product_inf_animal <- p_entry_product *
p_exposure *
p_contam_exposure *
p_indir_trans
})
animal_products_exp <- list(
product_inf_origin = inf_origin_exp,
products_entry = products_entry_exp
)
# Evaluate animal products model
animal_products <- eval_module(
exp = animal_products_exp,
data = animal_products_data
)
# Calculate pathway totals: aggregate across all products and import operations
animal_products <- trial_totals(
mcmodule = animal_products,
mc_names = "p_product_inf_animal",
trials_n = "import_size",
subsets_p = "p_inf_origin",
subsets_n = "import_operations"
)We can combine live animals and animal product pathways in one mcmodule using combine_modules(). Note that we used the expression “inf_origin” in both modules, so they’ll have some nodes with the same names. To avoid conflicts regarding name duplication, we can add a prefix to all the nodes in the pathway using add_prefix().
# Add prefixes to distinguish nodes by pathway
live_animals <- add_prefix(live_animals, prefix = "live")
animal_products <- add_prefix(animal_products, prefix = "products")
# Combine modules
animal_imports <- combine_modules(live_animals, animal_products)Now that we have both pathways in the same mcmodule, we can combine the probability of their totals to get the probability of entry from either pathway.
# Calculate combined probability of entry from either pathway
animal_imports <- at_least_one(
mcmodule = animal_imports,
mc_names = c("products_p_product_inf_animal_set", "live_p_entry_animal_set"),
name = "p_entry_total"
)Some groups in the products pathway were not found in the live animals pathway, as not all countries that export wool or hides export live animals (see Figure 4)
However at_least_one() can match the dimensions between both nodes, account for different scenarios independently, and assume that the probability of the pathway is 0 when the match is null. With this, we’ll have a complete mcnode combining nodes with very different dimensions while still propagating uncertainty (see Figure 5)
Aggregation totals
If we have an analysis by species and country, we might want to see the results by country alone. To do this, we calculate the combined probability across all subsets of variates for each country, following the element-wise logic of mcnode to maintain uncertainty propagation (see Figure 6). By default, it calculates the combined probability of independent events, but it can also take functions such as sum, product, or average.
Using this approach we can calculate the results by different aggregation levels for reporting.
# 1. Total probability by country and scenario
animal_imports <- agg_totals(
mcmodule = animal_imports,
mc_name = "p_entry_total",
agg_keys = c("country_code", "scenario_id"),
name = "p_entry_total_by_country"
)
# 2. Total probability by commodity and scenario
animal_imports <- agg_totals(
mcmodule = animal_imports,
mc_name = "p_entry_total",
agg_keys = c("commodity", "scenario_id"),
name = "p_entry_total_by_commodity"
)
# 3. Total probability by scenario only (overall risk aggregation)
animal_imports <- agg_totals(
mcmodule = animal_imports,
mc_name = "p_entry_total",
agg_keys = "scenario_id",
name = "p_entry_total_by_scenario"
)Module results and visualization
Node summary
We can show the results in different ways, the simplest is using an mc_summay(). This provides the summary statistics of an mc2d node adding the mcmodule keys.
# Summary statistics by scenario (overall risk)
mc_summary(animal_imports, "p_entry_total_by_scenario") mc_name scenario_id mean sd
1 p_entry_total_by_scenario 0 1.523511e-06 1.451889e-06
41 p_entry_total_by_scenario Destination surveillance 4.063556e-07 4.250882e-07
81 p_entry_total_by_scenario Origin surveillance 1.431641e-06 1.292199e-06
Min 2.5% 25% 50% 75%
1 2.768819e-08 1.453322e-07 5.433060e-07 1.028348e-06 2.030297e-06
41 1.118113e-08 3.643254e-08 1.482629e-07 2.824225e-07 5.058123e-07
81 3.434629e-08 1.677904e-07 5.666681e-07 1.071101e-06 1.855722e-06
97.5% Max nsv Na's
1 5.110509e-06 1.094106e-05 1001 0
41 1.598643e-06 3.777209e-06 1001 0
81 4.819469e-06 1.106696e-05 1001 0Network visualization
We can create a network representation of an mcmodule that plots every node in the model and lets you explore the result for the first variable it finds. This visualization is still experimental, but it has great potential to help with communication and debugging.
mc_network(animal_imports, legend = TRUE)Probability distribution plot
The mc_plot() function creates a boxplot combined with scattered points to represent the uncertainty dimension and show the probability distribution of each variate in an mcnode. It returns a ggplot object that can be further customised.
# Plot total probability of entry by scenario
mc_plot(animal_imports, "p_entry_total_by_scenario", order_by = "median")
# Plot total probability of entry by commodity
mc_plot(animal_imports, "p_entry_total_by_commodity", order_by = "median")
# Plot with log scale to better visualize differences
mc_plot(
animal_imports,
"p_entry_total_by_commodity",
scale = "log10",
order_by = "median"
)
Module diagnostics
Correlation analysis
The function mcmodule_corr() evaluates how variability and uncertainty in the inputs correlate with variability and uncertainty in the outputs. This helps identify the most influential inputs. It is based on mc2d::tornado() but applied to many distinct variates.
# Correlation betwen inputs and total output
corr_analysis <- mcmodule_corr(
mcmodule = animal_imports,
output = "p_entry_total",
progress = TRUE,
print_summary = TRUE,
method = "spearman"
)
[Correlation analysis] Expression animal_inf_origin, animal_entry (1/2)
[Correlation analysis] Expression product_inf_origin, products_entry (2/2)
=== Correlation Analysis Summary ===
Analysis Parameters:
- Analysis type: Global output
- Output node: p_entry_total
- Correlation method(s): spearman
- Missing value handling: all.obs
Expression Information:
- Module: live_animals
- Expression: animal_inf_origin, animal_entry
- Input nodes: live_p_herds_surveillance, live_p_animals_inspected, live_clinical_se_origin, live_import_qty, live_import_operations, live_clinical_se_dest
- Variates analyzed: 153
- Module: animal_products
- Expression: product_inf_origin, products_entry
- Input nodes: products_p_herds_surveillance, products_p_animals_inspected, products_clinical_se_origin, products_import_qty, products_import_operations
- Variates analyzed: 153
Results Summary:
- Total correlations calculated: 1836
- Top 5 most influential inputs (by absolute mean correlation):
1. live_import_qty: 0.0598
2. products_import_qty: 0.0454
3. products_import_operations: -0.0403
4. live_p_animals_inspected: -0.0353
5. live_clinical_se_dest: 0.0227
Input Correlation Strength Distribution:
- None: 64 (3.5%)
Inputs by Correlation Strength:
- None: live_p_herds_surveillance, live_p_animals_inspected, live_clinical_se_origin, live_clinical_se_dest, live_import_qty, live_import_operations, products_p_herds_surveillance, products_p_animals_inspected, products_clinical_se_origin, products_import_qty, products_import_operations
Warnings and Errors:
- Number of correlations with warnings: 1788 (97.39%)
- Input nodes with warnings: live_p_herds_surveillance, live_p_animals_inspected, live_clinical_se_origin, live_import_qty, live_import_operations, live_clinical_se_dest, products_p_herds_surveillance, products_p_animals_inspected, products_clinical_se_origin, products_import_qty, products_import_operations
- Unique warning/error types:
1. the standard deviation is zero; the standard deviation is zero
2. the standard deviation is zero; the standard deviation is zero; the standard dev...Convergence analysis
The mcmodule_converg() function assesses whether your Monte Carlo simulations have run for enough iterations to produce stable results. It tracks how statistics (mean, median, quantiles) change as iterations accumulate.
# Analyze convergence of live animals pathway
converg_live <- mcmodule_converg(
mcmodule = live_animals,
from_quantile = 0.95,
to_quantile = 1.0,
print_summary = TRUE
)
=== Convergence Analysis Summary ===
Analysis Parameters:
- Number of simulations: 1001
- Simulation quantile range: 0.95 to 1
- Simulations range: 950 to 1001 (51 simulations)
Convergence Results:
- Total nodes analyzed: 711
- Nodes with divergence below 0.001: 437 (61.46%)
- Nodes with divergence below 1% of their mean: 539 (75.81%)
- Nodes with divergence below 0.001 or 1% of their mean: 626 (88.05%)
- Nodes with divergence below 0.001 or 2.5% of their mean: 687 (96.62%)
- Nodes with divergence below 0.001 or 5% of their mean: 708 (99.58%)
Stochastic Distributions Stability:
- Maximum deviation of mean: 44.010307 (standardized: 0.023160)
- Maximum deviation of median: 73.284349 (standardized: 0.019482)
- Maximum deviation of 2.5% quantile: 53.112381 (standardized: 0.107945)
- Maximum deviation of 97.5% quantile: 505.442493 (standardized: 0.033336)
3 (0.42%) nodes did not converge at 5% threshold :(# Analyze convergence of animal products pathway
converg_products <- mcmodule_converg(
mcmodule = animal_products,
from_quantile = 0.99,
to_quantile = 1.0,
print_summary = TRUE
)
=== Convergence Analysis Summary ===
Analysis Parameters:
- Number of simulations: 1001
- Simulation quantile range: 0.99 to 1
- Simulations range: 990 to 1001 (11 simulations)
Convergence Results:
- Total nodes analyzed: 3377
- Nodes with divergence below 0.001: 2404 (71.19%)
- Nodes with divergence below 1% of their mean: 3146 (93.16%)
- Nodes with divergence below 0.001 or 1% of their mean: 3272 (96.89%)
- Nodes with divergence below 0.001 or 2.5% of their mean: 3342 (98.96%)
- Nodes with divergence below 0.001 or 5% of their mean: 3367 (99.70%)
Stochastic Distributions Stability:
- Maximum deviation of mean: 18054.706658 (standardized: 0.033319)
- Maximum deviation of median: 24791.939809 (standardized: 0.026880)
- Maximum deviation of 2.5% quantile: 1842.220445 (standardized: 0.163391)
- Maximum deviation of 97.5% quantile: 3448.129501 (standardized: 0.078528)
10 (0.30%) nodes did not converge at 5% threshold :(Tricks and tweaks
Replace missing or infinite values to prevent calculation errors using
mcnode_na_rm().Use the name parameter in functions like
at_least_one()andagg_totals()to customize output node nameSome functions can work outside modules:
create_mcnodes()andmc_summary()(useful for prototyping)Convert to other formats: Transform mcmodules into matrices or mc objects with
mcmodule_to_matrices()andmcmodule_to_mc()
Next directions
We are currently planning to add global sensitivity analysis based on Morris screening, using the existing local sensitivity (OAT) support
We are working on building compatibility with other R packages, to facilitate integration with broader risk assessment workflows
To stay updated watch our repository: https://github.com/NataliaCiria/mcmodule
We encourage you to explore the package further, adapt it to your own use cases, report bugs, and contribute feedback or improvements.
Acknowledgments
Thanks to my thesis supervisors Alberto Allepuz and Giovanna Ciaravino, and to all my colleages and friends for listening to me complain about my bugs.
mcmodule was developed with support from:
BIOSECURE, a European Union’s HORIZON Europe FARM2FORK project
Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or REA. Neither the European Union nor the granting authority can be held responsible for them.