Overview
This vignette demonstrates how Traditional Ecological Knowledge (TEK) can be incorporated into decision models using Bayesian Networks (BNs) and Monte Carlo sampling. It is designed as a short tutorial showing practical, reproducible code snippets and guidance on eliciting and encoding TEK.
Note: this vignette is both pedagogical and practical. It gives worked examples you can run locally. The goal is to show how TEK — whether collected as counts (k of n informants), multinomial categories, or qualitative rankings (high/med/low) — can be translated into probability distributions that feed Bayesian Networks. Monte Carlo sampling is then used to propagate TEK-derived uncertainty through the BN and produce decision-relevant outcome distributions.
The examples use real survey data, homegardens: plants
and their uses in 102 homegardens in southwestern Uganda, one garden per
informant (Whitney et al. 2017; Whitney et al.
2018). Please cite both papers when you use these data. Items the
survey did not collect (informant confidence, management judgements) are
simulated or illustrative and labelled as such. For heavier examples you
may conditionally skip chunks if packages are not installed (see notes
in chunks). By default the main code chunks here are enabled so the
vignette will run in an environment with the required packages
installed.
Required packages
We will show examples that use the following packages (install if you want to run the code):
- bnlearn
- gRain
- dplyr
- ggplot2
- purrr
You can install them with
install.packages(c("bnlearn","gRain","dplyr","ggplot2","purrr"))
1. Short motivation and data
We often collect TEK as counts (k informants out of n) or qualitative categories (high/medium/low). A convenient approach is to convert counts into Beta (binary) or Dirichlet (categorical) priors and then sample from these priors inside a Monte Carlo loop to create Conditional Probability Tables (CPTs) for a Bayesian Network (BN).
Below we illustrate the workflow with a tiny, didactic model: a
management Decision (Harvest vs Conserve), a latent
Abundance node (High/Low), and an Impact node
(Recovery/Decline). The goal is to compare management options under TEK
uncertainty.
2. Example: encode TEK as Beta priors and sample CPTs
library(ethnobotanyR)
# Informants are gardens. k of n gardens growing avocado (Persea americana) sell it.
avocado <- homegardens[homegardens$sp_name == "Persea americana", ]
tek_k <- sum(avocado$sale) # gardens that sell avocado
tek_n <- nrow(avocado) # gardens that grow avocado
# Beta posterior with a uniform prior: Beta(k+1, n-k+1)
alpha <- tek_k + 1
beta <- tek_n - tek_k + 1
# draw Monte Carlo samples for that probability
set.seed(123)
p_samples <- rbeta(1000, alpha, beta)
# Inspect distribution
library(ggplot2)
ggplot(data.frame(p = p_samples), aes(p)) +
geom_density() +
labs(x = "P(a garden growing avocado sells it)", title = "Beta posterior from survey counts")
Informant dataset from the survey
The avocado gardens give a small informant-level data set: - a binary response, whether the garden sells avocado (from the survey), - a categorical response, the breadth of avocado uses reported (1, 2 or 3+ uses; from the survey), and - a confidence score (1-5). The survey did not collect confidence, so this column is simulated to illustrate weighting.
set.seed(42)
uses <- setdiff(names(homegardens), c("informant", "sp_name"))
informants <- data.frame(
id = avocado$informant,
sells = avocado$sale,
breadth = cut(rowSums(avocado[uses]), c(0, 1, 2, Inf), labels = c("1", "2", "3+")),
confidence = sample(3:5, nrow(avocado), replace = TRUE, prob = c(0.2, 0.6, 0.2)) # simulated
)
knitr::kable(head(informants, 8))| id | sells | breadth | confidence |
|---|---|---|---|
| BUR1 | 1 | 2 | 3 |
| BUR10 | 1 | 2 | 3 |
| BUR11 | 0 | 1 | 4 |
| BUR2 | 0 | 1 | 3 |
| BUR3 | 0 | 1 | 5 |
| BUR4 | 0 | 1 | 4 |
| BUR5 | 0 | NA | 5 |
| BUR6 | 0 | 1 | 4 |
We use these data to illustrate aggregation approaches below: simple counts to Beta/Dirichlet, and weighted pooling by confidence.
3. Build a small BN and run Monte Carlo parameter sampling
High-level algorithm
- Define BN structure (nodes and conditional relationships).
- For each Monte Carlo iteration:
- Sample each uncertain probability from its Beta/Dirichlet prior.
- Populate CPTs with the sampled probabilities.
- Compile the BN and query target node probabilities or expected utilities.
- Aggregate the distribution of outcomes across iterations and compare decisions.
if (requireNamespace("bnlearn", quietly = TRUE) && requireNamespace("gRain", quietly = TRUE)) {
library(bnlearn); library(gRain); library(purrr)
} else {
message("Skipping BN + gRain workflow: 'bnlearn' and/or 'gRain' not available")
}## Loading required package: gRbase##
## Attaching package: 'gRbase'## The following objects are masked from 'package:bnlearn':
##
## ancestors, children, nodes, parents# ILLUSTRATIVE elicited judgement, not in the survey: 7 of 10 informants say harvest lowers abundance
tek_k <- 7 # number of informants saying 'decline'
tek_n <- 10 # total informants
alpha <- tek_k + 1
beta <- tek_n - tek_k + 1
# Structure: Decision -> Abundance -> Impact (Decision also acts on Impact)
# p_low_harvest is the probability that abundance is Low under harvest. The Beta prior
# below is uncertainty about that probability, taken here from informants' judgements;
# the share of informants who agree is a stand-in, not an observed frequency.
yn <- c("Yes", "No")
sample_one <- function(){
p_low_harvest <- rbeta(1, alpha, beta)
p_low_conserve <- 0.1 # toy value, fixed
cpt_decision <- cptable(~Decision, values = c(0.5, 0.5), levels = c("Harvest", "Conserve"))
cpt_abundance <- cptable(~Abundance | Decision,
values = c(1 - p_low_harvest, p_low_harvest, 1 - p_low_conserve, p_low_conserve),
levels = c("High", "Low"))
# Impact | Abundance, Decision: 8 values, ordered Recovery/Decline for each
# combination, Abundance varying fastest, then Decision
cpt_impact <- cptable(~Impact | Abundance + Decision,
values = c(0.95, 0.05, 0.40, 0.60, # Harvest: High, Low
0.98, 0.02, 0.70, 0.30), # Conserve: High, Low
levels = c("Recovery", "Decline"))
grain_net <- grain(compileCPT(list(cpt_decision, cpt_abundance, cpt_impact)))
# set each decision as evidence and read P(Recovery)
sapply(c("Harvest", "Conserve"), function(d) {
querygrain(setEvidence(grain_net, nodes = "Decision", states = d),
nodes = "Impact", type = "marginal")$Impact["Recovery"]
})
}
set.seed(1)
results <- as.data.frame(t(replicate(500, sample_one())))
colnames(results) <- c("Harvest", "Conserve")
# posterior of P(Recovery) per decision, and of the difference
apply(results, 2, quantile, probs = c(0.1, 0.5, 0.9))## Harvest Conserve
## 10% 0.4967453 0.952
## 50% 0.5782777 0.952
## 90% 0.6872939 0.952## 10% 50% 90%
## 0.2647061 0.3737223 0.45525473.1 Learn a network from the survey
When you have data rather than judgements, learn the structure and
fit the parameters from it. Here each row is a species in a garden;
nodes are life form, native status and three uses (Whitney et al. 2018).
boot.strength() shows how stable each arc is across
resamples; arcs below about 0.5 should not be interpreted.
library(bnlearn)
hg <- merge(homegardens, homegardens_species, by = "sp_name")
hg <- hg[hg$type != "Palm", ]
bn_data <- data.frame(
type = factor(hg$type), native = factor(hg$native, labels = c("introduced", "native")),
food = factor(hg$food), sale = factor(hg$sale), medicine = factor(hg$medicine))
set.seed(1)
strength <- boot.strength(bn_data, R = 200, algorithm = "hc")
strength[strength$strength > 0.5 & strength$direction > 0.5, ]## from to strength direction
## 5 native type 1.000 0.6350000
## 6 native food 1.000 0.5050000
## 8 native medicine 0.840 0.7410714
## 9 food type 1.000 0.6000000
## 12 food medicine 1.000 0.6300000
## 13 sale type 0.940 0.8430851
## 14 sale native 0.585 0.7820513
## 15 sale food 0.720 0.8055556
## 16 sale medicine 0.995 0.8341709net <- averaged.network(strength, threshold = 0.5)
fitted_net <- bn.fit(net, bn_data, method = "bayes", iss = 10) # Dirichlet priorQuery the fitted network, for example how likely a native tree is to have a medicinal use:
g <- bnlearn::as.grain(fitted_net)
q <- gRain::querygrain(gRain::setEvidence(g, nodes = c("type", "native"),
states = c("Tree", "native")),
nodes = "medicine", type = "marginal")
q$medicine## medicine
## 0 1
## 0.5524984 0.4475016Rows are species-in-gardens, so the same garden or species appears many times and these intervals are too narrow. Treat this as exploration and use a mixed model (see the modeling vignette) for inference.
4. Decision rules and visualization
- Compute expected utility per decision for each Monte Carlo iteration, then compare mean, median and quantiles.
- Use conservative decision rules if stakeholders prefer lower risk (e.g., maximize the 10th percentile of utility).
Visualization ideas:
- Density plots of expected utility for each option.
- Fan chart showing quantiles across iterations.
5. Practical notes on elicitation and aggregation
- For binary outcomes from TEK use Beta priors (k + 1, n - k + 1).
- For multinomial TEK responses use Dirichlet with counts + 1 as pseudo-counts.
- Weight informants by confidence or expertise when aggregating (weighted counts or hierarchical models).
- Record provenance and do sensitivity by omitting individual informants.
5.1 Realistic aggregation examples
We now show practical aggregation methods using the
informants data above.
- Simple pooling (counts -> Beta/Dirichlet)
# Binary: counts for 'sells'
k <- sum(informants$sells)
n <- nrow(informants)
alpha <- k + 1; beta <- n - k + 1
cat(sprintf("P(sells) prior ~ Beta(%d, %d) -> mean = %.2f\n", alpha, beta, alpha/(alpha+beta)))## P(sells) prior ~ Beta(22, 71) -> mean = 0.24# Multinomial: counts for breadth categories -> Dirichlet
br_counts <- table(informants$breadth)
br_dirichlet_alpha <- as.numeric(br_counts) + 1
cat("Breadth counts:\n")## Breadth counts:##
## 1 2 3+
## 57 22 6- Weighted pooling by informant confidence
Confidence is simulated here; use it only if you collect it.
# Use confidence as weight
weights <- informants$confidence / sum(informants$confidence)
# Weighted proportion, with the effective sample size (Kish) so that
# weighting does not make the prior look more certain than n informants allow
p_w <- sum(informants$sells * weights)
n_eff <- 1 / sum(weights^2)
walpha <- p_w * n_eff + 1; wbeta <- (1 - p_w) * n_eff + 1
cat(sprintf("Weighted Beta prior approx: Beta(%.1f, %.1f), n_eff = %.1f -> mean = %.2f\n", walpha, wbeta, n_eff, walpha/(walpha+wbeta)))## Weighted Beta prior approx: Beta(20.7, 69.7), n_eff = 88.5 -> mean = 0.23# Weighted multinomial: use weighted frequencies
weighted_br <- tapply(weights, informants$breadth, sum)
weighted_br_counts <- weighted_br * n_eff
cat("Weighted breadth pseudo-counts:\n")## Weighted breadth pseudo-counts:## 1 2 3+
## 56.314976 20.478173 5.363331- Qualitative judgements (e.g. a stated confidence) as informant strength
Do not map response categories to pseudo-counts (e.g. High = 5, Low = 1): that reweights the categories, not the informants, and inflates the sample size. Weight informants as in method 2, and use the Dirichlet with the effective sample size.
w <- informants$confidence / sum(informants$confidence)
p_br <- tapply(w, informants$breadth, sum)
qual_alpha <- p_br * n_eff + 1
print(round(qual_alpha, 2))## 1 2 3+
## 57.31 21.48 6.36All three approaches give you numeric (pseudo-)counts that can be used as Beta or Dirichlet parameters. Which is appropriate depends on your elicitation protocol and whether you want to weight informants by confidence or expertise.