When forecasting high-volatility financial assets like Bitcoin (BTC-USD), standard point-prediction models—whether classical ARIMA or deep neural networks—often fall short. A point estimate tells you where the model expects the price to go, but it completely ignores market regime shifts and changing risk levels.

In this article, we construct an end-to-end forecasting pipeline in R that combines Machine Learning (Multi-Layer Perceptron via brulee and tidymodels) for directional trend prediction with Financial Econometrics (GJR-GARCH via rugarch) to generate dynamic, volatility-adjusted confidence bounds.

The Architecture Overview

Instead of attempting to forecast level prices directly—which introduces non-stationarity issues—our pipeline operates on log returns and reconstructs price trajectories dynamically.

Key Architectural Decisions & Methodological Rationale

1. Modeling Log Returns Instead of Level Prices

Directly predicting raw prices often creates spurious regression results because price series are non-stationary. By transforming the price series to log returns:

This transformation achieves stationarity, making the target variable suitable for neural network training. Factual price predictions are reconstructed at test time via cumulative exponentiation:

2. Deep Learning for Trend: Multi-Layer Perceptron (brulee)

We use a 2-layer MLP architecture (32, 16) initialized via the tidymodels framework.

• Lag Features: We supply 20 days of lagged log returns to capture short- to medium-term autocorrelation.
• Normalizations & Dropping Zero Variance: Handled within a clean recipes pipeline to prevent data leakage

3. Econometric Volatility Modeling: GJR-GARCH(1,1) on Neural Network Residuals

Machine learning models capture the conditional mean, but their residuals, $\epsilon_t = r_t – \hat{r}_t$ , still contain volatility clustering—periods of high turbulence followed by high turbulence, and calm periods followed by calm.

To model this conditional variance ($\sigma_t^2$), we fit a GJR-GARCH(1,1) model on the MLP training residuals:

• Why GJR-GARCH over Standard GARCH? Crypto markets exhibit leverage effects: price drops often spark higher panic volatility than equivalent price rallies. The indicator term $I_{t-1}$ captures this asymmetry ($\gamma$).
• Why garchOrder = c(1, 1)? The order (1,1) specifies a single lag for the ARCH component ($\epsilon_{t-1}^2$) and a single lag for the GARCH component ($\sigma_{t-1}^2$). Empirical literature demonstrates that a single lag captures the vast majority of volatility persistence without overfitting.
• Why Skewed Student-t (sstd) Distribution? Bitcoin returns violate Gaussian normal distribution assumptions. They display fat tails (extreme events occur frequently) and skewness (asymmetric distribution). The sstd distribution explicitly models both heavy tails and skewness.
• Solver Selection: We use solver = "hybrid" in rugarch to combine optimization algorithms, ensuring convergence even under complex distribution parameters.

When the market is quiet, these confidence bounds contract automatically. During volatile regimes, the bounds expand, giving a realistic probabilistic envelope rather than static, uninformative confidence intervals.

Complete R Implementation

Below is the complete, self-contained R script implementing the full pipeline using tidymodels, brulee, rugarch, and plotly.

# Integrated Forecasting Pipeline with tidymodels, yardstick, and plotly
library(tidymodels)
library(timetk)
library(brulee)
library(modeltime)
library(tidyquant)
library(plotly)
library(scales)
library(glue)
library(ggtext)
library(rugarch) # Added for GJR-GARCH volatility modeling

# 1. Fetch Data & Calculate Log Returns
df_raw <- tq_get("BTC-USD") %>%
  select(date, close) %>%
  filter(date >= last(date) - months(24)) %>%
  drop_na() %>%
  mutate(log_ret = log(close / lag(close))) %>%
  drop_na()

# 2. Add Lag Features on Log Returns
df_prepared <- df_raw %>%
  tk_augment_lags(log_ret, .lags = 1:20) %>%
  drop_na()

# 3. Time Series Split (10 Assessment Window)
split <- time_series_split(
  df_prepared,
  assess = 10,
  cumulative = TRUE
)

df_train <- training(split)
df_test  <- testing(split)

# 4. Recipe Setup (Predict Log Return, Not Level Price)
recipe_mlp <- recipe(log_ret ~ ., data = df_train) %>%
  step_rm(date, close) %>%
  step_dummy(all_nominal_predictors()) %>%
  step_zv(all_predictors()) %>%
  step_normalize(all_numeric_predictors())

# modeltime::naive_reg() requires the date column as a predictor to parse the time index
recipe_naive <- recipe(log_ret ~ date, data = df_train)

# 5. Model Specifications
model_spec_mlp <- mlp(
  hidden_units = c(32, 16),
  epochs = 100,
  dropout = 0.05,
  learn_rate = 0.005
) %>%
  set_engine("brulee") %>%
  set_mode("regression")

model_spec_naive <- naive_reg() %>%
  set_engine("naive")

# 6. Workflow Construction & Model Fitting
workflow_mlp <- workflow() %>%
  add_recipe(recipe_mlp) %>%
  add_model(model_spec_mlp) %>%
  fit(df_train)

workflow_naive <- workflow() %>%
  add_recipe(recipe_naive) %>%
  add_model(model_spec_naive) %>%
  fit(df_train)

# 7. Predictions & Price Reconstruction (Reference Code Logic)
preds_mlp_ret   <- predict(workflow_mlp, df_test)
preds_naive_ret <- predict(workflow_naive, df_test)

df_results <- df_test %>%
  mutate(
    yesterday_close = lag(close, default = first(close)),
    actual      = close,
    pred_brulee = yesterday_close * exp(preds_mlp_ret$.pred),
    pred_naive  = yesterday_close * exp(preds_naive_ret$.pred)
  )

# ==============================================================================
# 7.5. GJR-GARCH (Skewed Student-t) Dynamic Volatility & Confidence Bands
# ==============================================================================
# Calculate MLP residuals on the training set
train_preds_mlp <- predict(workflow_mlp, df_train)$.pred
train_residuals <- df_train$log_ret - train_preds_mlp

# GJR-GARCH(1,1) specification with Skewed Student-t distribution for fat-tails and asymmetry
garch_spec <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = c(1, 1)),
  mean.model = list(armaOrder = c(0, 0), include.mean = FALSE),
  distribution.model = "sstd"
)

# Fit GJR-GARCH model on training residuals
garch_fit <- ugarchfit(spec = garch_spec, data = train_residuals, solver = "hybrid")

# Forecast conditional volatility for the 10-day test horizon
garch_fore <- ugarchforecast(garch_fit, n.ahead = 10)
dynamic_sigmas <- sigma(garch_fore) # 10-day dynamic conditional standard deviations

# Attach dynamic standard deviations to test results and compute dynamic confidence bounds
df_results <- df_results %>%
  mutate(
    conditional_sigma = as.vector(dynamic_sigmas),
    # Dynamically scaled confidence bounds in price space (± 1.96 * conditional_sigma)
    conf_hi = pred_brulee * exp(1.96 * conditional_sigma),
    conf_lo = pred_brulee * exp(-1.96 * conditional_sigma)
  )

# 8. Accuracy Table Construction (Yardstick on Original Price Scale)
accuracy_table <- tibble(
  .model_id   = c(1, 2),
  .model_desc = c("BRULEE MLP (32, 16) + GJR-GARCH", "NAIVE"),
  rmse        = c(rmse_vec(df_results$actual, df_results$pred_brulee),
                  rmse_vec(df_results$actual, df_results$pred_naive)),
  rsq         = c(rsq_vec(df_results$actual, df_results$pred_brulee),
                  rsq_vec(df_results$actual, df_results$pred_naive)),
  mape        = c(mape_vec(df_results$actual, df_results$pred_brulee),
                  mape_vec(df_results$actual, df_results$pred_naive)),
  mae         = c(mae_vec(df_results$actual, df_results$pred_brulee),
                  mae_vec(df_results$actual, df_results$pred_naive))
)

# Print Model Comparison Table
cat("=== MODEL ACCURACY EVALUATION TABLE (PRICE SCALE) ===\n")
cat(knitr::kable(accuracy_table, digits = 4, format = "simple"), sep = "\n")

# 9. Modern Interactive Plotly Visualization (Reference Style)
brulee_mape <- accuracy_table$mape[1]
df_eval <- df_results

# Custom Hover Text Data Frames
df_plot_actual <- df_eval %>% select(date, actual)  %>% mutate(text_actual = glue("<b>Actual Price:</b> {scales::dollar(actual, accuracy = 0.01)}<br><b>Date:</b> {format(date, '%b %d, %Y')}"))
df_plot_pred   <- df_eval %>% select(date, pred_brulee) %>% mutate(text_pred   = glue("<b>AI Prediction Price:</b> {scales::dollar(pred_brulee, accuracy = 0.01)}<br><b>Date:</b> {format(date, '%b %d, %Y')}"))
df_plot_hi     <- df_eval %>% select(date, conf_hi) %>% mutate(text_hi     = glue("<b>GJR-GARCH Upper Bound:</b> {scales::dollar(conf_hi, accuracy = 0.01)}<br><b>Date:</b> {format(date, '%b %d, %Y')}"))
df_plot_lo     <- df_eval %>% select(date, conf_lo) %>% mutate(text_lo     = glue("<b>GJR-GARCH Lower Bound:</b> {scales::dollar(conf_lo, accuracy = 0.01)}<br><b>Date:</b> {format(date, '%b %d, %Y')}"))

# Build Base GGPlot Layer
p <- ggplot() +
  geom_ribbon(data = df_eval, aes(x = date, ymin = conf_lo, ymax = conf_hi), fill = "#808080", alpha = 0.18) +
  geom_point(data = df_plot_hi, aes(x = date, y = conf_hi, text = text_hi), color = "transparent", alpha = 0, size = 3) +
  geom_point(data = df_plot_lo, aes(x = date, y = conf_lo, text = text_lo), color = "transparent", alpha = 0, size = 3) +
  geom_line(data = df_plot_actual, aes(x = date, y = actual), color = "#2c3e50", linewidth = 1.2) +
  geom_point(data = df_plot_actual, aes(x = date, y = actual, text = text_actual), color = "#2c3e50", size = 2) +
  geom_line(data = df_plot_pred, aes(x = date, y = pred_brulee), color = "#e74c3c", linetype = "dashed", linewidth = 1.2) +
  geom_point(data = df_plot_pred, aes(x = date, y = pred_brulee, text = text_pred), color = "#e74c3c", size = 2) +
  scale_y_continuous(labels = dollar_format(accuracy = 1)) +
  labs(
    x = "", y = "",
    title = paste0(
      "<span style='color:#2c3e50'>BTC-USD Prices</span> vs ",
      "<span style='color:#e74c3c'>MLP + GJR-GARCH Dynamic Forecast</span><br>",
      "<span style='font-size:12px; color:#555555;'>10-Day Horizon | Out-of-Sample MAPE: ", round(brulee_mape, 2), "</span>"
    )
  ) +
  theme_minimal() +
  theme(
    plot.title = element_markdown(hjust = 0.5, face = "bold"),
    plot.background = element_rect(fill = "#ffffff", color = NA),
    panel.background = element_rect(fill = "#ffffff", color = NA),
    panel.grid.minor = element_blank()
  )

# Typography Settings & Plotly Output
font_family <- list(family = "Roboto Slab, Sans-Serif", size = 16)
label_font  <- list(font = list(family = "Roboto Slab, Sans-Serif", size = 13))

ggplotly(p, tooltip = "text") %>%
  style(hoverlabel = label_font) %>%
  layout(font = font_family)

Conclusion

Combining machine learning with traditional econometrics produces a hybrid system that overcomes the limitations of both fields:

1. Machine Learning captures non-linear relationships and trend movements.
2. Econometrics (GARCH) supplies explicit statistical parameters for uncertainty and volatility clustering.

The resulting output provides realistic, dynamic risk envelopes that evolve with changing market regimes—giving quantitative traders and analysts a robust framework for volatility-aware decision-making.

Mark Padgham is leaving rOpenSci at the end of October, and we want to take this moment to recognize his impact and say thank you for helping to shape rOpenSci.

Mark joined rOpenSci as a staff member in December 2019 to help us expand software peer review to statistical topics. However, Mark’s impact on rOpenSci has gone far beyond this one project. Mark has not only shaped how rOpenSci reviews statistical software, but has also built much of the tooling that supports the peer review processes, and R packages which not only support but also strengthen our community. If you have submitted a package to rOpenSci, reviewed or edited one, you have almost certainly used something Mark made!

Making a Mark

Mark’s most visible contribution is to statistical software peer review. Reviewing a package that implements statistical methods is a different job from reviewing a package that accesses an API. Mark tackled the task of translating what “good” means in a statistical package and how a reviewer can assess it. He recruited a new set of editors with the expertise to tackel this work and turned insights into a concrete set of statistical software standards and the guide that accompanies them (rOpenSci Statistical Software Testing and Peer Review Community Call). Mark also built the machinery to go with this new type review: pkgstats gathers information about a package’s structure, helps authors and editors see what a submission contains and how it is organised; srr allows package authors to document their adherence to statistical software standards.

The result of all this work was our first peer-reviewed statistical packages, published in November 2022, and to date we’ve had over 26 statistical packages go through this review! Anyone who has gone through statistical software review has benefited from the groundwork Mark laid.

The engine of software review

Peer review at rOpenSci is a human process, but it runs on a lot of automation, and a good part of that automation is Mark’s work.

Mark built and maintained pkgcheck, which runs automated checks on submitted packages, and the GitHub Action that allows anyone to run the same checks on their own repository. pkgcheck uses goodpractice, another package maintained by Mark, that recently rolled out a host of updates and extension. Mark maintained Buffy, the bot behind our chat review system on GitHub, which lets editors assign reviewers, track deadlines and move submissions along with simple comments, and many other pieces of infrastructure that make the whole system work.

Review is a human process

In early 2025, Mark became our new Software Review Lead. In this role he was in charge of editorial management, curating the current editorial board and planning its future. Mark helped guide peer-review through new challenges, such as how review should handle AI-assisted contributions, resulting in rOpenSci’s first AI policy; or how to deal with large software submissions, resulting in our first group peer review.

Tools for everyone

Mark’s care for rOpenSci goes beyond software. He researched alternatives to GitHub, sharing the results in Beyond GitHub, a thoughtful look at what it would mean for a community like ours to depend less on a single platform.

He has also served on our Code of Conduct committee, work that helps ensure rOpenSci remains a welcoming and safe community.

Mark also builds tools that help the wider R community. The allcontributors package, inspired by the allcontributors project, makes it easy to recognise everyone who contributes to a project, not just those who write code. His work on repometrics explores how to understand the health and activity of software repositories, and his work on deposits looks at how research data and software can be deposited in trusted repositories.

Thank you, Mark

Mark’s work has shaped how we review software, support contributors, and build a more reliable and sustainable open science ecosystem. Thank you, Mark, for building so much of the infrastructure and process that we rely on, and for leaving it in better shape than you found it.

If you have benefited from Mark’s work, reach out to Mark and let him know! You can keep following his work through his website and social media.

Je dirige seul Qognito, une entreprise de conseil et de formation. Pour réaliser certaines activités de mon cœur de métier, je coordonne plusieurs équipes d’agents IA, mes « squads ». C’est ma manière de faire fonctionner une firme minimale : déléguer une partie du travail tout en restant responsable des orientations, des livrables et des engagements pris auprès de mes clients.

La veille sur l’IA alimente ce travail. Une note peut m’aider à actualiser une formation, à examiner le besoin d’un prospect ou à préparer une proposition. Sa valeur tient à ce que je peux en faire ensuite.

Je conserve ces notes dans Obsidian, une application de gestion de connaissances. Des liens internes, les wikilinks, les relient entre elles : je peux retrouver une analyse et la rapprocher d’autres sujets. Ces liens facilitent l’exploration ; ils ne garantissent pas que plusieurs notes apportent des preuves indépendantes. Une erreur peut ainsi devenir un appui pour un raisonnement ultérieur.

J’avais organisé la production avec soin. Un agent prépare la source, un autre rédige l’analyse, un troisième évalue la note. Je décide de sa publication dans ma base. Voici le circuit attendu :

Ce schéma décrit les responsabilités prévues. Il reste à vérifier que les outils respectent ce circuit.

Dans ma formation SAGA-IA, je demande aux décideurs d’examiner ce qu’ils délèguent, les preuves disponibles et le coût du contrôle. J’ai soumis mon propre dispositif à ces exigences. L’audit a commencé par une découverte très simple.

Les règles pouvaient changer avec celui qui les appliquait

Les fichiers qui définissaient les consignes, les critères du juge et les contrôles étaient modifiables par les agents. Celui qui produisait une note avait donc techniquement la possibilité de changer les critères utilisés pour l’évaluer.

Je n’avais pas constaté qu’un agent l’avait fait. Mais la séparation des rôles reposait sur une consigne que les droits d’accès ne faisaient pas respecter.

J’ai fait protéger ces fichiers. Un agent peut proposer une modification ; sa mise en place exige une intervention humaine avec les droits d’administration. Dix tentatives de modification, par des voies différentes, ont ensuite été refusées. Les règles sont restées intactes.

Ces premières corrections précédaient la campagne consacrée aux protections, appelée « Phase A » dans mon audit. Elles avaient déplacé la question : les règles tenaient mieux, mais que pouvait-il encore arriver aux documents sur lesquels je m’appuyais ?

Publier une note ne suffisait pas à la préserver

Au départ, une note rejoignait la base dès sa rédaction, avant le verdict du juge. J’avais ensuite instauré une étape d’attente, puis une autorisation humaine de publication. Le premier cycle réel avait suivi ce circuit : neuf sources traitées, neuf notes examinées et publiées par mes soins.

Pourtant, une note déjà publiée pouvait encore être modifiée. Une source aussi. Le contrôle existant ne repérait pas ces changements. Une approbation donnée hier ne disait donc pas ce que je retrouverais demain.

La Phase A a traité ce problème. Sur une copie de travail, trois altérations volontaires — une note modifiée, une note ajoutée hors circuit et une transcription modifiée — ont été détectées par un nouveau contrôle comparant les documents à une référence conservée. Des verrous ont également été posés sur les notes publiées et les sources préparées.

J’ai accepté une discipline pour moi aussi : plus de retouche directe dans Obsidian. Les changements doivent passer par le circuit prévu et laisser une trace. Sinon, comment distinguer ma correction d’une modification anormale ?

Préserver le contenu ne suffit toutefois pas à le rendre juste. Lors des premiers essais, le juge avait relevé « usually » traduit par « toujours », alors que les contrôles automatiques n’avaient rien signalé. Une nuance devenue certitude peut changer une décision.

Le programme compare les documents ; le juge examine le sens ; je conserve la décision de publication. Cette répartition rend les responsabilités identifiables. Elle ne dispense pas de vérifier la qualité du jugement.

J’ai dû examiner mes propres examens

Un banc d’essai est une répétition sur des cas choisis, dont certains contiennent volontairement des pièges. J’en avais utilisé pour éprouver la veille. Deux constats m’ont obligé à revoir la confiance accordée à leurs résultats.

Le premier concernait le corrigé : les réponses attendues restaient accessibles aux agents. Imaginez un examen dont les candidats pourraient ouvrir les solutions. Cela fragilise l’évaluation, même sans preuve qu’ils les aient consultées. L’accès au corrigé a été protégé, puis son refus éprouvé sur le poste.

Le second concernait les permissions. La copie utilisée pour certains essais était plus restrictive que l’environnement réel. Une action pouvait y demander un accord alors qu’elle passait librement au quotidien. Les résultats décrivaient bien les essais réalisés ; ils ne suffisaient pas à établir la même protection en usage réel.

La dernière étape de la Phase A devait donc réunir les protections dans une épreuve avec les permissions de production. Il fallait observer les tentatives, identifier ce qui les arrêtait et vérifier leurs effets sur les documents.

Vingt tentatives, et une limite qui reste visible

Le 19 octobre, vingt tentatives ont été exécutées et tracées dans les deux environnements d’agents utilisés. Les essais de modification de documents protégés et d’accès au corrigé ont rencontré des refus.

Deux créations de fichiers hors du circuit de publication ont toutefois réussi dans l’un des environnements. Le contrôle les a détectées ; elles ont ensuite été supprimées. La vérification finale a confirmé le retour de la base à son état initial.

Vingt tentatives documentées ne signifient donc pas vingt actions empêchées. Dans le fonctionnement courant, un ajout hors circuit peut rester présent jusqu’au contrôle déclenché par la prochaine publication. Cette limite reste ouverte.

La Phase A est close, avec une réserve précise : le prochain cycle réel doit encore montrer qu’une publication autorisée aboutit à une note effectivement verrouillée. Les refus constatés pendant les essais ne remplacent pas cette vérification de bout en bout.

L’évaluation de qualité menée en parallèle a donné sept cas réussis sur dix, avec une fausse alerte de fidélité et deux cas dont la portée restait à trancher. Protéger les documents n’établit pas que leur contenu s’est amélioré. Un examen rétrospectif des notes publiées reste prévu.

Le temps que le chronomètre ne voyait pas

Pour comprendre le coût de cette maîtrise, il faut revenir au premier cycle réel. L’outil avait enregistré une à trois secondes par publication. J’avais pourtant consacré neuf minutes à relire les neuf notes avant de l’ouvrir. Ce temps a dû être déclaré séparément.

Ces neuf minutes suffisaient à ma décision de publier ; elles ne constituaient pas une comparaison complète de chaque note avec sa source. Une validation humaine doit préciser ce qui a réellement été vérifié et avec quel temps disponible.

La même exigence vaut pour les corrections. Pendant les premiers travaux, j’ai accepté un défaut de recoupement compté deux fois : une seule occurrence avait été constatée, et le script envisagé pour la traiter aurait été coûteux et lui-même faillible. Chaque anomalie ne justifie pas une règle supplémentaire.

J’ai aussi fixé des critères d’arrêt, des responsables et des conditions de réexamen. Ces dispositions devront vivre dans l’usage. À ce stade, je n’ai démontré aucun gain global de rentabilité : la relecture, les reprises, la maintenance et l’audit font partie du coût à apprécier.

Ce que je peux désormais justifier

Le chemin parcouru tient en trois étapes :

Les protections testées apportent des preuves circonscrites. La dernière étape reste à observer.

Je peux mieux expliquer ce que je confie aux agents, ce qui résiste aux tentatives testées et ce qui demande encore un contrôle. Ma responsabilité de dirigeant consiste aussi à décider quelles incertitudes sont acceptables et lesquelles doivent limiter la délégation.

Si vous coordonnez vous aussi une activité avec des agents, mon livre blanc gratuit La facture fantôme de l’IA propose une première lecture sur la supervision, les coûts complets et l’usage du temps gagné.

Pour travailler ces décisions, SAGA-IA — Volume 1 réunit le diagnostic de délégation et l’évaluation de la viabilité économique, avec des cas et des exercices. C’est le prolongement de ce REX : apprendre à justifier ce que l’on engage, les preuves que l’on demande et les conditions auxquelles on poursuivra.

I often need to get mainland from country / administrative region polygons, so I finally decided to write a function to help automate it. It selects the N largest feature(s) in a polygons or lines map. Here’s the code:

selectLargest <- function(x,  # terra SpatVector polygons or lines                          n = 1,  # number of largest features to retain                          dissolve = FALSE,  # whether to group contiguous features                          singleparts = TRUE  # whether to disaggregate multi-part features) {  stopifnot(inherits(x, "SpatVector"),            terra::geomtype(x) %in% c("polygons", "lines"))  if (dissolve) {    x <- terra::aggregate(x)  }  if (singleparts) {    x <- terra::disagg(x)  }  if (terra::is.polygons(x)) {    size_rank <- order(terra::expanse(x), decreasing = TRUE)  } else if (terra::is.lines(x)) {    size_rank <- order(terra::perim(x), decreasing = TRUE)  }  x[size_rank[1:n], ]}

The dissolving of internal borders is optional. Here are some usage examples:

library(terra)library(geodata)# import a USA map and get the conterminous and continental parts:usa <- geodata::gadm("United States", level = 1, path = tempdir())mainland_usa <- selectLargest(usa, dissolve = TRUE)contin_usa <- selectLargest(usa, dissolve = TRUE, n = 2)  # two largest polygons, to include Alaskaplot(usa, main = "USA")plot(mainland_usa, main = "Conterminous USA")plot(contin_usa, main = "Continental USA")# import a UK map and get Great Britain:uk <- geodata::gadm("United Kingdom", path = tempdir())great_britain <- selectLargest(uk, dissolve = TRUE)plot(uk, main = "UK")plot(great_britain, main = "Great Britain")# import a Portugal + Spain map and get the Iberian Peninsula:ptes <- geodata::gadm(c("Portugal", "Spain"), path = tempdir())iberia <- selectLargest(ptes, dissolve = TRUE)plot(ptes, main = "Portugal and Spain")plot(iberia, main = "Iberian Peninsula")

As usual, feedback is welcome!

TLDR: For the life of me, I couldn’t get Ctrl+Shift+M to work for inserting pipes (|>) when editing a Quarto document in Positron visual mode. After countless hours fiddling with keybindings.json, I developed a functional, but obtuse workaround for Windows. It’s not pretty, but it works. This post outlines the approach in case it helps others facing the same problem.

Note: this problem is already on the radar of the development team, so will no doubt be fixed soon. I’ve also only outlined a workaround for Windows. Feel free to reach out if you’ve found a better solution (or a solution that works on other platforms) and I’ll add this to the post.

The problem

I’m a big fan of Quarto as it provides a seamless way to mix data, code and narrative text without ever leaving my IDE. I’ve also been trying to move from RStudio to Positron as I find it snappier, works better with AI tools and works well with Python. Unfortunately, Positron is still experiencing growing pains as a newish IDE, including some shortcuts not working when editing Quarto documents in visual mode.

This isn’t a new issue (e.g: 1, 2 and 3), so is sure to be fixed in the not-too-distant future (likely just after I post this). But, after countless hours messing around with keybindings.json and hassling people from Posit (apologies, Isabella) I asked Claude if it had a solution.

It did.

And it works.

So, I thought I’d share it here, despite the fact it doesn’t spark joy with the minimalist inside of me:

Duct tape and thumbtacks: intercept the shortcut before it reaches Positron

As far as I can tell, the problem comes from Positron / VSCode intercepting some shortcuts when Quarto is in visual mode, which is probably why adjusting the settings did nothing for me. As a (temporary) solution to this, Claude suggested installing a program called AutoHotkey v2 and creating a script operates when Positron is active.

Obviously, as an independent-minded person I did what I’m told and here’s how you can do the same:

1. Install AutoHotkey v2

Get Version 2 from autohotkey.com (the code below will not run under v1).

2. Create the script

Open Autohotkey2 (or AutoHotKey Dash) and select New Script to create a new blank .ahk script in a secluded folder of your choice. Give it a name like ‘positron-keys.ahk’ so future you knows what it’s doing there and use a text editor to replace the contents of the file with the code below (or the equivalent depending on your shortcut preferences):

(SendText sends the literal characters to the Quarto document, so nothing gets re-interpreted by Positron, while #HotIf restricts the behavior to Positron, so it doesn’t interfere with other applications.)

#Requires AutoHotkey v2.0
#SingleInstance Force

#HotIf WinActive("ahk_exe Positron.exe")

^+m::SendText(" |> ")    ;Ctrl+Shift+M -> native pipe
!-::SendText(" <- ")     ;Alt+-        -> assignment arrow

#HotIf

3. Run it

To run the code and launch AutoHotkey, open the ‘positron-keys.ahk’ file. If it worked, a green H should appear in the system tray.

4. Test it

To test it you can now open Positron and try using Ctrl + Shift + M to insert a pipe in a Quarto when in visual mode. Hopefully it works, as this is the only idea covered in this post.

5. Make it persistent

To avoid needing to load the code each time you want to use Positron, you can then add a shortcut to the script that loads at startup. For Windows, you can open the startup folder by typing Win+R and executing the command ‘shell:startup‘. Right-drag positron-keys.ahk into the folder and choose Create shortcut here.

Note: if you have other solutions, please let me know. I’ll happily update this post if others find it useful.

AI use: Claude suggested this solution and wrote the first draft of this post, which I heavily edited.

The post Getting the pipe shortcut to work in Positron appeared first on Giles.

In June 2026, rOpenSci marked its 15th anniversary with three community celebration events: a themed coworking event and two extra events. This post captures the reflections and ideas shared by participants across the two special gatherings named “rOpenSci Community Celebration: Then, Now & Next” that were held to accommodate different time zones.

How the events worked

Our goal was to create a space for people to have conversations, get to know each other and also reflect on rOpenSci’s past, present and future.

After a welcome and a brief introduction, participants joined three successive breakout rounds of 15 minutes each, one for each theme: the past, the present, and the future. Between rounds, attendees took part in a Scavenger Hunt (a regular feature of our co-working sessions), finding something in their physical surroundings that matched a given prompt.

The Scavenger Hunt is a fun moment together, and we also used it as a natural prompt for introductions at the start of the next breakout room — “what did you find?” becomes a nice and easy opening for conversations.

As in our community calls, workshops and co-working sessions, a shared document served as the place where everyone could take notes in real time, add thoughts, and keep contributing for a couple of days after the event ended. Noam Ross named this “process is product” approach, caring about how things happen, not just what happens, as one of rOpenSci’s defining traits.

Round 1 – The Past: What made rOpenSci special?

The first breakout asked participants how they first connected with rOpenSci and what their favourite memory was.

Several people traced their connection back to a talk about rOpenSci or to a personal introduction from a community member who made them feel welcome from the start.

The unconferences were remembered with particular warmth and excitement: the chance to sit next to people, see how they work, talk about R, and be part of something as it was taking shape. “Meeting my heroes”, as one participant put it.

Other things people named as favourite:

• The peer-review process: praised by editors, reviewers and authors alike, both for the quality of the process and for the connections it created.

Being invited to write a blog post after a first package review … made me feel very cool and a valued part of the community.” – Alex Koiter

• The dev guide: a detailed and evolving documentation with knowledge that isn’t easily found elsewhere.

• The Champions Program and the mini-hackathons, spaces that provide friendly and guided entry points for open source contributors and R package developers.

• Co-working sessions are a place were you can find “your people,” a place where you don’t need to explain what R is.

It’s a satisfying community ‘distraction/procrastination’ from our normal jobs, that keeps us from being lonely and brings us together from different worlds to work on something that’s benefitting people.” – Ronny Hernandez Mora

Round 2 – The Present: What does rOpenSci do well?

The second round invited reflection on what rOpenSci does well today and what participants want to get more involved in.

The most common answer was creating a welcoming space. Participants then commented on the different places where they experience the friendly and welcoming community:

• Co-working came up again here, noted as something other communities have adopted from rOpenSci.

• The Slack workspace and the atmosphere it creates: a common space for learning, getting advice, sharing ideas, and feeling part of something larger.

• Good documentation for package submission, which lowers the barrier to participating.

• Some participants expressed interest in getting more involved in the Champions Program, along with honest questions about whether the entry points are accessible enough.

• The sense that rOpenSci’s quality in different activities comes from its community participation.

You can jump into a role and the community will support you. If there is anything you want to do, start doing it and you will get the support you need.” – Adam Sparks

Round 3 – The Future: What comes next?

The third round asked what participants would like to see rOpenSci do next, and what ideas from other communities could be worth adopting or partnering on.

In-person and hybrid gatherings were a recurring theme: local meetups, local chapters, satellite events at larger R conferences like useR!, and the return of unconferences, potentially co-organised with partner communities.

Other ideas from across the rooms:

• First-time contributor hackathons, possibly with a GenAI focus

• A 24-hour unconference modeled on formats used by other communities like SORTEE

• Career panels showing the range of roles in research software engineering and data-centric work

• An rOpenSci-led package for R package development, like usethis, that supports GitLab, Codeberg, and other platforms beyond GitHub

• Deeper coalitions with open science organizations

• Connecting package authors with graphic designers to help them with their image and visual identity

• Exploring Zarr in R as an area where rOpenSci could help shape the ecosystem

• Better serving very beginners, the people who are still asking “what is R, anyway?”

• Merchandise and swag, we all want our rOpenSci hex sticker

One thread in the conversation pointed to a real tension: rOpenSci has spent 15 years building carefully and steadily. GenAI is moving much faster. How do you maintain that care and those foundations when the pace of change accelerates? Even when we didn’t have any easy answers, this question is the right one to be discussing with the community.

Closing thoughts

Reading through the notes from both events, a few things are consistent.

People talk about rOpenSci not just as a set of tools or infrastructure, but as a place where they found colleagues, had their work taken seriously, and felt like valued contributors.

The software peer-review process, the dev guide, the co-working sessions, the Champions Program, these work because the community behind them is invested in supporting them and find them useful.

At 15 years, the ideas that surfaced – local chapters, hackathons, coalitions, better onboarding for newcomers – suggest a community that still wants to try new things.

It works because everyone is welcome.” – several attendees at the events

Thank you to everyone who joined and shared their thoughts.

These events were part of rOpenSci’s 15th anniversary celebration series. Keep an eye on rOpenSci’s events calendar, our blog and our newsletter for more content around our anniversary.

On 22 July 2026 at 15:05, I made the first commit to a new repository. It held an empty R package and one Excel workbook with field data from a household waste study in Kampala, Uganda. Three minutes later I sent the first prompt to Claude Code:

I am building a new data package and received a dataset in @data-raw/ to process for clean up. […] strictly start with cleaning the data, writing a dictionary. […] make a plan first.

On 24 July at 00:18, the repository held version 1.0.0 of solidwastekampala with a Zenodo DOI: a tidy and validated dataset of 103 households, a data dictionary for 37 variables, documentation, a website, machine-readable metadata, a citation with the ORCID iDs of the authors, and automated checks on five platforms.

My own time on it was no more than about three hours, spread over four sessions. The repository records every step in its commits, issues and pull requests, and in an archive of every prompt I sent.

In version 0.39.0 of the ahead package, the conformalize() function has been updated to include a new parameter called split_prob. This parameter allows users to specify the proportion of the data that should be used for calibration when constructing conformal prediction intervals.

0 – Install packages

install.packages("pak")

pak::pak("Techtonique/ahead")

library(ahead)

1 – USAccDeaths

par(mfrow=c(2, 2))
obj <- ahead::conformalize(FUN = forecast::ets, y = USAccDeaths,
                           h = 20L, split_prob = 0.4)
plot(obj)
obj <- ahead::conformalize(FUN = forecast::ets, y = USAccDeaths,
                           h = 20L, split_prob = 0.3)
plot(obj)
obj <- ahead::conformalize(FUN = forecast::ets, y = USAccDeaths,
                           h = 20L, split_prob = 0.6)
plot(obj)
obj <- ahead::conformalize(FUN = forecast::ets, y = USAccDeaths,
                           h = 20L, split_prob = 0.7)
plot(obj)

2 - AirPassengers

par(mfrow=c(2, 2))
obj <- ahead::conformalize(FUN = forecast::ets, y = AirPassengers,
                           h = 20L, split_prob = 0.4)
plot(obj)
obj <- ahead::conformalize(FUN = forecast::ets, y = AirPassengers,
                           h = 20L, split_prob = 0.3)
plot(obj)
obj <- ahead::conformalize(FUN = forecast::ets, y = AirPassengers,
                           h = 20L, split_prob = 0.6)
plot(obj)
obj <- ahead::conformalize(FUN = forecast::ets, y = AirPassengers,
                           h = 20L, split_prob = 0.7)
plot(obj)

3 - fdeaths

par(mfrow=c(2, 2))
obj <- ahead::conformalize(FUN = ahead::dynrmf, y = fdeaths,
                           h = 20L, split_prob = 0.4)
plot(obj)
obj <- ahead::conformalize(FUN = ahead::dynrmf, y = fdeaths,
                           h = 20L, split_prob = 0.3)
plot(obj)
obj <- ahead::conformalize(FUN = ahead::dynrmf, y = fdeaths,
                           h = 20L, split_prob = 0.6)
plot(obj)
obj <- ahead::conformalize(FUN = ahead::dynrmf, y = fdeaths,
                           h = 20L, split_prob = 0.7)
plot(obj)

4 - UKDriverDeaths

par(mfrow=c(2, 2))
obj <- ahead::conformalize(FUN = ahead::dynrmf, y = UKDriverDeaths,
                           h = 20L, split_prob = 0.4)
plot(obj)
obj <- ahead::conformalize(FUN = ahead::dynrmf, y = UKDriverDeaths,
                           h = 20L, split_prob = 0.3)
plot(obj)
obj <- ahead::conformalize(FUN = ahead::dynrmf, y = UKDriverDeaths,
                           h = 20L, split_prob = 0.6)
plot(obj)
obj <- ahead::conformalize(FUN = ahead::dynrmf, y = UKDriverDeaths,
                           h = 20L, split_prob = 0.7)
plot(obj)

5 - WWWusage

par(mfrow=c(2, 2))
obj <- ahead::conformalize(FUN = ahead::dynrmf, y = WWWusage,
                           h = 20L, split_prob = 0.4)
plot(obj)
obj <- ahead::conformalize(FUN = ahead::dynrmf, y = WWWusage,
                           h = 20L, split_prob = 0.3)
plot(obj)
obj <- ahead::conformalize(FUN = ahead::dynrmf, y = WWWusage,
                           h = 20L, split_prob = 0.6)
plot(obj)
obj <- ahead::conformalize(FUN = ahead::dynrmf, y = WWWusage,
                           h = 20L, split_prob = 0.5)
plot(obj)

Thanks to the awesomeness of webR and Shinylive, the tool I made in back in 2014 to estimate your blood alcohol concentration (BAC) now runs 100% in the browser, no server needed!

Try the browser version of drinkR here!

library(knitr)
library(ggplot2)
library(dplyr)
library(tidyr)
library(broom)

set.seed(123)
n <- 100
X1 <- rnorm(n, mean = 0, sd = 1)
X2 <- rnorm(n, mean = 0, sd = 1)
Z <- 5 * X2 + rnorm(n, mean = 0, sd = 1)
Y <- 3 * X1 + 2 * X2 + rnorm(n, mean = 0, sd = 1)

kable(cor(cbind(X1, X2, Z)), caption="Correlation of the predictors")

fit1 <- lm(Y ~ X1 + X2)
fit2 <- lm(Y ~ X1 + X2 + Z)
coef_results <- rbind(
  tidy(fit1) %>% mutate(model = "fit1"),
  tidy(fit2) %>% mutate(model = "fit2")
)
# move model to the first column
coef_results <- coef_results %>% select(model, everything())
# now R-squared and adjusted R-squared
rsq_results <- rbind(
  glance(fit1) %>% mutate(model = "fit1"),
  glance(fit2) %>% mutate(model = "fit2")
) %>%
  select(model, everything())
kable(coef_results, caption = "The coefficients and statistics")

kable(rsq_results, caption = "The goodness of fit")

simulate_correlated_data<-function(rho, n=100, seed=NULL) {
  if(is.null(seed)) {
    seed < Sys.time()
  }
  set.seed(seed)
  x <- rnorm(n)
  z <- rnorm(n)
  y <- rho*x + sqrt(1-rho^2)*z
  return(cbind(x,y))
}

n<-100
rhos<-c(0.1, 0.3, 0.5, 0.8, 0.95)
simulated_data <- lapply(rhos, simulate_correlated_data, n=n)
# let's calculate the correlation coefficients for the simulated data
simulated_rhos<-sapply(simulated_data, function(d) cor(d[,"x"],d[, "y"]))
kable(cbind(expected_rho=rhos, actual_rho=simulated_rhos))

 fits <- lapply(simulated_data, function(d) {
          x1<-d[,"x"]
          x2<-d[, "y"]
          y <- 10*x1 + 3*x2 + rnorm(n)
          fit<-lm(y~x1 + x2)
          return(fit)
          } 
          )
coef_results<-lapply(seq_along(rhos),
                     function(i) {
                       tidy(fits[[i]]) %>%
                         mutate(rho = rhos[i]) %>%
                         select(rho, everything())
                     }
                     )
coef_results <- do.call(rbind, coef_results)
kable(coef_results, caption = "Regression statistics with different correlation strengths among predictors")

coef_results <- coef_results %>% 
    mutate(rho = factor(rho)) %>%
    pivot_longer(
      cols = estimate:p.value,
      names_to = "metric",
      values_to = "value"
    )
plt <- ggplot(coef_results, aes(x=rho, y=value)) +
      geom_point(aes(color=term)) +
      facet_wrap(~ metric, ncol=2, scales="free_y") +
      theme(legend.position = "top")
plt

If you plan to learn R or Python in the new year, book before the end of November. Our spring 2027 data science training courses run from January to June, and early bird pricing on the first one ends on 29 November.

There are 22 courses on the schedule. All of them run live online, and the trainers are data scientists and engineers who spend the rest of their week on client projects. Each course runs over two afternoons, 13:30 to 17:00 UK time. The Bayesian inference course is the exception: it runs over four.

Three routes through our spring 2027 training courses

The courses are scheduled so you can take them in order. There are three obvious routes, each ending with a set of skills you can use at work the following week.

Learning R from scratch

Start with Introduction to R on 11–12 January. Data Wrangling in the Tidyverse follows on 25–26 January, then Data Visualisation with ggplot2 on 1–2 February and Programming with R on 8–9 February. That gets you from your first script to writing your own functions in a month.

Why take this route: Clean, chart and automate your data in code. Monthly reporting becomes a script you rerun in minutes, with fewer manual errors.

January too soon? The same four courses run again between 12 April and 11 May. If you are still deciding whether R is worth the time, we wrote about why learning R is a good career move.

Learning Python

Introduction to Python runs on 15–16 February. Programming with Python follows on 22–23 February, and Data Visualisation with Python on 1–2 March. The same three run again in May and June.

Why take this route: Analyse data in pandas, automate repeat tasks and produce report-ready charts, in the same language as your engineering team. It is also the groundwork for machine learning in Python.

Statistics and machine learning in R

This route is for people who already use R. Statistical Modelling with R on 8–9 March covers hypothesis testing, regression, clustering and principal components analysis. Machine Learning with Tidymodels on 15–16 March and Advanced Machine Learning with Tidymodels on 22–23 March take you through pre-processing, model fitting and tuning.

Why take this route: Build, validate and tune predictive models in tidymodels, and explain which variables drive them. Your team moves from reporting what happened to forecasting what happens next, with evidence for how far to trust it.

Two courses that run once this spring

Introduction to Bayesian Inference using RStan runs over four afternoons, 18–21 January. You write Stan programs for a range of statistical models and learn to read the MCMC diagnostics in R. It is the only date for it in the first half of 2027.

Reporting with Quarto runs on 21–22 June. You build reports, slides and documents straight from R or Python code, so the numbers in the report change when the data does. It pairs well with the ideas in our post on reproducible analytical pipelines.

Early bird deadlines

Early bird pricing ends about six weeks before each course. Booking closes a week before the start date. The first three deadlines are:

• 29 November for Introduction to R
• 6 December for Bayesian inference
• 13 December for Data Wrangling in the Tidyverse

Dates for the second half of 2027 will follow.

Using this year’s training budget

If your team has training budget left this December, booking a January course now spends it before the year ends. It also gets you the early bird price, because every January deadline falls before Christmas.

What you get on every course

• PDF notes and scripts to keep
• Live demonstrations and hands-on exercises
• Time with the trainer to ask questions about your own work

We have run more than 1,000 courses, for organisations including NHS Scotland, Shell, Wessex Water and the Royal Statistical Society.

Training a whole team

Every course in our full catalogue can also run privately for one team, using your data and your tools. Group bookings and returning clients get a discount. Email training@jumpingrivers.com with the size of the team and what they need to learn.

See the full public schedule and book a place.

For updates and revisions to this article, see the original post

The Model Context Protocol (MCP) is quietly becoming one of the most useful ideas in AI tooling. Here, we describe what MCP is, why it matters, and how to build an MCP server yourself.

MCP is a standard way to link a Large Language Model (LLM) to the outside world. An MCP server is a small piece of software that exposes some data or capability to an AI model. Everyday examples include connecting your AI agent to your calendar or your file system, or allowing it to search the web. These are things an LLM cannot do (or would struggle to do) by itself. On its own, a model only knows what it learned during training, plus whatever you provide it, and it has no way to check, fetch or act on anything beyond that. Even then, it can only give out text responses, which you have to be the one to act on. An MCP server can allow an LLM to carry out actions itself. Less talking, more doing.

Why not just use an API?

There are many ways to connect external capabilities to an AI model, usually via an API (Application Programming Interface), and this was the industry standard. However, APIs are all configured slightly differently, and once you start to scale up the number of connections you need between multiple pieces of software, you hit a combinatorial problem: the number of custom connections a product needs multiplies by the number of tools that use it.

Anthropic introduced MCP in late 2024 as an open standard for connecting tools and data to LLM applications. It works like a USB port on a computer: any MCP server can be used with any agent/model, whether you are using Claude, Gemini, ChatGPT or one you’ve made yourself.

Building your own MCP server

Full instructions on how to build an MCP server are available in the MCP documentation, but here is a very quick overview using Claude Code and the Sonnet 5 model.

The good news is that a working MCP server can be as short as a few lines of code. Here’s the shape of it, step by step, using Python.

Step 1: Install the dependencies

The MCP project maintains official SDKs for Python, TypeScript and other languages. The Python ‘mcp’ package embodies this SDK. However, here we are going to use the ‘fastmcp’ package. This is built upon ‘mcp’, but has some benefits (in particular, ‘fastmcp’-based tools can readily be deployed on Posit Connect, see later).

mkdir -p ~/projects/
cd ~/projects
uv init --python=3.15 dice-roll
# Initialized project `dice-roll`
# at `/home/<user>/projects/dice-roll`
cd dice-roll
uv add fastmcp

This gives us the following project structure, which is typical for a ‘uv’ project. We will define our MCP server in ./src/dice_roll/__init__.py, although if the project got much larger, we would split it into separate modules or sub-packages.

tree -L 2
.
├── pyproject.toml
├── README.md
├── src
│   └── dice_roll
└── uv.lock

Step 2: Define a tool

An MCP server exposes “tools”: functions the model can call. Here is a mini example: a server with a single tool that rolls a dice.

# src/dice_roll/__init__.py
from typing import Annotated
import random

from fastmcp import FastMCP
from pydantic import Field

# Give your server a name: this is what shows up when a model connects
mcp = FastMCP("dice-pick")

@mcp.tool
def roll_dice(
 sides: Annotated[int, Field(description="Number of sides on the dice", ge=2)] = 6,
) -> int:
 """Roll a single fair dice with the given number of sides and return the number rolled."""
 return random.randint(1, sides)

def main() -> None:
 mcp.run()

The @mcp.tool docstring (“Roll a single …”) is important. This becomes the description of the tool that the LLM reads. It is like the blurb of a book: the LLM will look at the description and from it decide whether it wishes to use the tool or not. If this description isn’t clear, it may lead to the tool not being used.

The type annotations for both the parameter sides and the value returned by roll_dice() also provide useful metadata for the LLM. Our sides parameter has a constraint that ensures at least 2 sides are present on the dice (ge=2; this is enforced within Python by Pydantic). When describing the tools provided by this MCP server, an LLM will be able to explain that the minimum number of sides is 2 (and more advanced models may not call the tool otherwise).

Step 3: Connect it

Save the code in the project, and we will test whether an MCP-compatible agent can use it. There are a few ways to register our locally-defined MCP server with Claude Code.

Since we have used the fastmcp framework, there is a command line tool that installs the MCP server for use from the project directory:

# From the ~/projects/dice-roll/ directory
uv run \
 fastmcp install claude-code src/dice_roll/__init__.py:mcp --with-editable .

The uv run part is necessary to ensure the correct Python environment is used. The fastmcp command needs to find the mcp object that we have defined in our source code. If we open claude and ask “What dice-pick tools do you know about?”, it replies:

The dice-pick MCP server gives me one tool:
roll_dice: rolls one fair dice and returns the number it lands on.
┌───────────┬─────────┬─────────┬───────────────────┐
│ Parameter | Type │ Default │ Constraint │
├───────────┼─────────┼─────────┼───────────────────┤
│ sides │ integer │ 6 │ must be 2 or more │
└───────────┴─────────┴─────────┴───────────────────┘
So I can roll a d6, d20, d100 or any other single dice. It can't roll several dice in one call, add
modifiers or pick from a list. For something like 3d6, I'd call it three times.

Then we can use the tool:

❯ Roll a dice
Called dice-pick
I rolled a 6 on a standard six-sided dice.
...
❯ Roll a 20-side dice
Called dice-pick
I rolled a 20 on the d20.
That's the highest number on both rolls so far. [... blah blah blah ...]
...
❯ Roll a d12
Called dice-pick
I rolled a 4 on the d12, so the run of top scores is over.

We have installed our locally in-development MCP server on our machine, but it’s only usable from the ~/projects/dice_roll directory. This is fine while you are testing, but you may want to use your MCP server in other locations once it works well.

Our project is uv-based. If we tell it which module and function should be called to start the project running, we can make the project more easily usable from other locations. But, the entry-point function used to start our MCP server running is dice_roll.main(), which uv has already added as the entry-point to the project in the pyproject.toml:

[project.scripts]
dice-roll = "dice_roll:main"

On running uv run python or uv sync, a dice-roll entry-point will be added to the virtual environment for the project.

For Claude Code, we can then register our MCP server with a single command:

claude mcp add --scope user dice-pick \
 -- /abs/path/to/projects/dice-roll/.venv/bin/dice-roll

Now, you can ask claude to “roll a 20-sided dice” from any directory on your computer and it will call your tool directly rather than guessing at an answer. This is possible because claude mcp add ... appends to the mcpServers entry of a file (~/.claude.json).

A few tips before you build something bigger:

• Keep tool descriptions (the docstrings for your @mcp.tool-decorated functions) short and specific. The model relies on the docstring to decide when to use a tool, so vague descriptions lead to it being ignored or misused.
• Think about what should require confirmation. Read-only tools (like looking something up) are usually safe to let a model call freely; tools that change or delete something are worth gating behind an explicit user confirmation step.

You can remove the MCP server from your claude setup by using

claude mcp remove dice-pick

Or just ask Claude to remove it!

Deploying to Posit Connect

A server running on your own machine is fine for testing, but it is only useful to you. MCP servers can be deployed to a range of different cloud providers. Posit Connect, the platform we normally reach for when publishing Shiny apps, Plumber APIs and Quarto reports, added first-class support for hosting MCP servers in 2026, complete with a dedicated icon and its own tab in the content list.

Posit Connect does expect you to write your Python-based MCP servers using the fastmcp framework that we used above. Connect also needs to import a runnable application object that can handle HTTP-based requests. When deployed, ‘fastmcp’-based MCP servers can expose an ASGI app object that handles these requests. Locally run servers use ‘stdio’ routing, but we only have to modify our current source code slightly to deploy to Posit Connect.

# imports ...

mcp = ...

def roll_dice():
 ...

app = mcp.http_app() # For deployment on Connect

def main():
 mcp.run() # For running using stdio locally

For deployment of the MCP server to Posit Connect, we also need to create a requirements.txt file. These are used to specify the run-time (and sometimes development) dependencies of a project. Since we used uv to manage the project, this is easy to create:

uv export --no-dev --no-hashes -o requirements.txt

The requirements.txt file contains a line -e . that dynamically includes any source code from src/. The Python environment that is built for your MCP server depends upon requirements.txt and is rebuilt when that file is updated.

We use the rsconnect-python package (and the rsconnect command line tool that it installs) to interact with Posit Connect. You may first have to register with your Connect system using rsconnect add ... or rsconnect login ....

Then you can deploy with the rsconnect tool, using the FastAPI deployment mode and pointing the entrypoint at the app object that we have just added to the source code. This is done using a --entrypoint argument of the form module:object_name. Here this is dice_roll:app, since our app object is in the dice_roll module:

rsconnect deploy fastapi -n <nickname-of-posit-server> --entrypoint dice_roll:app ./

The Posit Publisher extension for VS Code and Positron works too, if you’d rather deploy from your editor. Posit has some good documentation on publishing MCP servers.

On first deployment, Connect scans requirements.txt, recognises fastmcp, and automatically sets the content type to MCP. Once deployed, it keeps at least one process warm so there is no cold-start delay, and gives it its own icon in the content list, with no extra configuration needed.

Using your deployed MCP server

To use an MCP server that is deployed in the cloud, the process is very similar to how we added our locally-defined MCP server. On Posit Connect (or wherever you have published), you will have a URL for the MCP server and you can configure claude to use this URL as follows:

claude mcp add --transport http <mcp-name> https://<connect-server>/content/<content-hash>/mcp

Here we used dice-pick as the name of the MCP server, consistent with the name in the source code. When you next open claude you will be prompted to authenticate before you can use the tool.

Want to see this in action? Watch this walkthrough, recorded live at our AI in Production 2026 conference, of an MCP server being built and used from scratch.

Key takeaways

• An MCP server is the standard way to give an AI model access to things it can’t do alone: checking a calendar, reading a file, searching the web, or looking up live data.
• It works like a USB port: build the connection once, and any compatible model can use it.
• A working MCP server can be as short as a few lines of code. The model calls your tools directly instead of guessing at an answer.
• Keep tool descriptions short and specific, and gate anything destructive behind explicit confirmation.
• Posit Connect can host your MCP server directly, so your whole team can use it rather than just your own machine.

For updates and revisions to this article, see the original post

One of the fundamental aims of field experiments is to explorare the relationships between several experimental traits, which have been recorded/measured altogheter for all the experimental units. For example, think about a genotype experiment with 27 winter wheat genotypes in a Randomised Complete Block Design with three replicates. If we aim at assessing yield and yield quality, we have to measure, in each plot, several traits, such as grain yield and the weight of thousand kernels (TKW). Apart from separately analysing these two variables, we might be interested in describing their reciprocal relationships, by using the Pearson’s correlation coefficient. Let’s do this, by using the dataset ‘WheatQuality’, which is available in the ‘statforbiology’ package; we have 81 plots and just as many couples of measures in all.

library(Hmisc)
library(dplyr)
library(statforbiology)
dataset <- getAgroData("WheatQuality")
head(dataset)
##     Genotype Block Height  TKW Whectol Yield
## 1 arcobaleno     1     90 44.5    83.2 64.40
## 2 arcobaleno     2     90 42.8    82.2 60.58
## 3 arcobaleno     3     88 42.7    83.1 59.42
## 4       baio     1     80 40.6    81.8 51.93
## 5       baio     2     75 42.7    81.3 51.34
## 6       baio     3     76 41.1    81.1 47.78

Now, we consider the whole lot of measures and calculate the correlation coefficient between yield and TKW; we also use a t-test for the hypothesis that the correlation coefficient is equal to 0. This is the result:

ra <- with(dataset, rcorr(Yield, TKW) )
# 
# Correlation coefficient
ra$r[1,2] 
## [1] 0.540957
#
# P-level
ra$P[1,2] 
## [1] 1.850931e-07

The conclusion is: the correlation coefficient is positive and significantly different from 0, which supports the idea that plots with a high value on yield tend to have a high value on TKW, as well.

Such an approach is very often seen also in high ranking Journals, but: is it correct? Can we trust the conclusions? The answer is no, not completely. And here is why. Please, note that in this post I am following the point made by Bland and Altman (1994), by using an example that is more closely related to plant/crop science.

# Correlations within genotypes
wgCors <- by(dataset, list(dataset$Genotype), function(x) cor(x$Yield, x$TKW))
wgCors
## : arcobaleno
## [1] 0.9847228
## ------------------------------------------------------------ 
## : baio
## [1] 0.1611952
## ------------------------------------------------------------ 
## : claudio
## [1] -0.9993872
## ------------------------------------------------------------ 
## : colorado
## [1] 0.9837293
## ------------------------------------------------------------ 
## : colosseo
## [1] 0.4564855
## ------------------------------------------------------------ 
## : creso
## [1] -0.5910193
## ------------------------------------------------------------ 
## : duilio
## [1] -0.988233
## ------------------------------------------------------------ 
## : gianni
## [1] -0.7603802
## ------------------------------------------------------------ 
## : giotto
## [1] 0.9520211
## ------------------------------------------------------------ 
## : grazia
## [1] 0.4980828
## ------------------------------------------------------------ 
## : iride
## [1] 0.7563338
## ------------------------------------------------------------ 
## : meridano
## [1] 0.1174342
## ------------------------------------------------------------ 
## : neodur
## [1] 0.4805871
## ------------------------------------------------------------ 
## : orobel
## [1] 0.8907754
## ------------------------------------------------------------ 
## : pietrafitta
## [1] -0.9633891
## ------------------------------------------------------------ 
## : portobello
## [1] 0.9210135
## ------------------------------------------------------------ 
## : portofino
## [1] -0.9900764
## ------------------------------------------------------------ 
## : portorico
## [1] 0.1394211
## ------------------------------------------------------------ 
## : preco
## [1] 0.9007067
## ------------------------------------------------------------ 
## : quadrato
## [1] -0.5840238
## ------------------------------------------------------------ 
## : sancarlo
## [1] -0.646067
## ------------------------------------------------------------ 
## : simeto
## [1] -0.4051779
## ------------------------------------------------------------ 
## : solex
## [1] -0.6066363
## ------------------------------------------------------------ 
## : terrabianca
## [1] -0.4076416
## ------------------------------------------------------------ 
## : verdi
## [1] 0.5801404
## ------------------------------------------------------------ 
## : vesuvio
## [1] -0.7797493
## ------------------------------------------------------------ 
## : vitromax
## [1] -0.8056514

# Correlations within blocks
wbCors <- by(dataset, list(dataset$Block), function(x) cor(x$Yield, x$TKW))
wbCors
## : 1
## [1] 0.5998137
## ------------------------------------------------------------ 
## : 2
## [1] 0.539999
## ------------------------------------------------------------ 
## : 3
## [1] 0.5370398

mod1 <- lm(Yield ~ factor(Block) + Genotype, data = dataset)
mod2 <- lm(TKW ~ factor(Block) + Genotype, data = dataset)
wCor <- rcorr(residuals(mod1), residuals(mod2))
wCor$r
##            x          y
## x 1.00000000 0.03133109
## y 0.03133109 1.00000000
wCor$P
##           x         y
## x        NA 0.7812693
## y 0.7812693        NA

means <- aggregate(dataset[,c(4,6)], list(dataset$Genotype), mean)
rgPear <- rcorr( as.matrix(means[,2:3]) )
rgPear$r
##             TKW     Yield
## TKW   1.0000000 0.5855966
## Yield 0.5855966 1.0000000
rgPear$P
##              TKW      Yield
## TKW           NA 0.00133149
## Yield 0.00133149         NA

means <- aggregate(dataset[,c(4,6)], list(dataset$Block), mean)
rbPear <- cor( as.matrix(means[,2:3]) )
rbPear
##               TKW       Yield
## TKW    1.00000000 -0.08812544
## Yield -0.08812544  1.00000000

You can read the original post in its original format on Rtask website by ThinkR here: {brochure}: multipage {shiny} apps, now on CRAN

A {shiny} app is a single page: one url, one UI, one server function. That is comfortable while the app stays small. As soon as you want an actual site — a home page, a contact page, one page per user, urls people can bookmark or send by email — you end up faking pages with tabs and conditionalPanel().

{brochure} takes the problem the other way around: every page is a real page, with its own url, its own UI, its own server function and its own Shiny session.

install.packages("brochure")
library(shiny)
library(brochure)

brochureApp(
  page(
    href = "/",
    ui = fluidPage(
      h1("Home"),
      tags$a(href = "/contact", "Contact us")
    ),
    server = function(input, output, session) {
      # ...
    }
  ),
  page(
    href = "/contact",
    ui = fluidPage(h1("Contact"))
  )
)

You run it as you would any Shiny app, and moving from one page to another is a plain link. A url can also carry a parameter — /user/:id, which the page reads back with get_keys().

Around that, {brochure} gives you what you expect from a web server: request and response handlers to log, serve a healthcheck or answer before Shiny ever runs, helpers to set and read cookies (and so carry state from one page to the next), redirections, and what it takes to serve the app behind a reverse proxy or on Posit Connect.

The trade-off: you leave the single page paradigm behind, and nothing computed on one page survives a navigation to another. That is what vignette("design") is about — worth reading before you start.

The full post, with examples and what this release brings:

{brochure} 1.0.0 is now on CRAN

Documentation: colinfay.me/brochure.

This post is better presented on its original ThinkR website here: {brochure}: multipage {shiny} apps, now on CRAN

The 1.4.0 release of maestro is out, and with it comes the ability to record custom metrics from your pipelines. Execution time and success/failure were already tracked, but there was no good way to answer questions like “how many rows did this pipeline process?” or “how many API calls did it take?” Now there is.

If you haven’t heard of maestro, it’s a pipeline orchestration package. You can learn more about it here.

Get it from CRAN:

install.packages("maestro")

#' ./pipelines/metric_example.R
#' @maestroFrequency 1 hour
extract_flights <- function() {
  
  # Pretends to fetch data from an API
  df <- data.frame(
    flight_id = 1:5,
    airline_code = c("UA", "AC", "AC", "AA", "DE"),
    departure_time = as.POSIXct("2026-10-14 12:00:00") + c(100, 450, 750, 1450, 1750)
  )
  
  update_metric("rows_read", nrow(df))
  update_metric("bytes_read", as.numeric(object.size(df)))
  
  df
}

#' @maestroInputs extract_flights
load_flights <- function(.input) {
  
  write.csv(.input, tempfile(fileext = ".csv"))
  
  update_metric("rows_written", nrow(.input))
}

library(maestro)

schedule <- build_schedule(quiet = TRUE)

status <- run_schedule(
  schedule,
  orch_frequency = "1 hour",
  quiet = TRUE
)

get_metrics(schedule)

# A tibble: 3 × 4
  pipe_name       run_id metric       value
  <chr>           <chr>  <chr>        <dbl>
1 extract_flights vPnJoC rows_read        5
2 extract_flights vPnJoC bytes_read    1648
3 load_flights    ulLl5m rows_written     5

status <- run_schedule(
  schedule,
  orch_frequency = "1 hour",
  quiet = TRUE
)

get_metrics(schedule)

# A tibble: 6 × 4
  pipe_name       run_id metric       value
  <chr>           <chr>  <chr>        <dbl>
1 extract_flights vPnJoC rows_read        5
2 extract_flights vPnJoC bytes_read    1648
3 extract_flights bilN7s rows_read        5
4 extract_flights bilN7s bytes_read    1648
5 load_flights    ulLl5m rows_written     5
6 load_flights    R9iMP0 rows_written     5

#' ./pipelines/metric_increment.R
#' @maestroFrequency 1 day
fetch_paginated <- function() {
  
  results <- list()
  
  for (page in 1:3) {
    attempts <- 0
    
    repeat {
      attempts <- attempts + 1
      update_metric("attempts", 1, increment = TRUE)
      
      # Pretends to call an API that sometimes fails
      ok <- runif(1) > 0.5
      if (ok || attempts >= 3) break
    }
    
    update_metric("api_calls", 1, increment = TRUE)
    results[[page]] <- data.frame(page = page, value = rnorm(10))
  }
  
  df <- do.call(rbind, results)
  
  # Overwriting (the default) is fine for a one-off value
  update_metric("rows_read", nrow(df))
  
  df
}

Flow cytometry researchers have a growing collection of R and Bioconductor tools available to them. But for someone who is comfortable working with cytometry data and has little programming experience, knowing where to begin can be difficult. That is the gap the free Cytometry in R course aims to fill.

© 2026 Bioconductor. Content is published under Creative Commons CC-BY-4.0 License for the text and BSD 3-Clause License for any code. | R-Bloggers

When we have a multi-trait table of data, we might be interested in investigating variability structure and the reciprocal relationships (whether causal or non-causal) among variables. For example, let’s consider the ‘mtcars’ dataset, shown below:

data(mtcars)
head(mtcars)
##                    mpg cyl disp  hp drat    wt  qsec vs am gear carb
## Mazda RX4         21.0   6  160 110 3.90 2.620 16.46  0  1    4    4
## Mazda RX4 Wag     21.0   6  160 110 3.90 2.875 17.02  0  1    4    4
## Datsun 710        22.8   4  108  93 3.85 2.320 18.61  1  1    4    1
## Hornet 4 Drive    21.4   6  258 110 3.08 3.215 19.44  1  0    3    1
## Hornet Sportabout 18.7   8  360 175 3.15 3.440 17.02  0  0    3    2
## Valiant           18.1   6  225 105 2.76 3.460 20.22  1  0    3    1

Considering the different cars (along the rows), we might be interested in investigating the structure of variability for fuel consumption (mpg), number of cylinders (cyl), displacement (disp) and horsepower (hp). How much does each variable change across the different cars? How much does each variable change when another variable changes? Is there any kind of joint variation, either positive or negative?

One suitable statistic to answer these question is the covariance, that is:

where and are the means for each variable and is the number of individuals. In words, covariance is the sum of products of the i^th deviations from the mean for both variables. In this setting, the variance is just the covariance of a variable with itself.

For a data table, variances and covariances can be calculated by using the cov() function, that returns a matrix (the variance-covariance matrix, indeed):

matr <- mtcars[,1:4]

# Covariances
Sigma <- cov(matr)
Sigma
##              mpg        cyl       disp        hp
## mpg    36.324103  -9.172379  -633.0972 -320.7321
## cyl    -9.172379   3.189516   199.6603  101.9315
## disp -633.097208 199.660282 15360.7998 6721.1587
## hp   -320.732056 101.931452  6721.1587 4700.8669

While covariances indicate direction (positive or negative association), their values depend on the scale of measurement, making it difficult to gauge the strength of the relationship at a glance. For this reason, covariances are usually ‘standardized’ by dividing by the product of the standard deviations to obtain Pearson’s correlation coefficient:

where and are the standard deviations for and . Unlike covariances, correlations go from -1 to 1 and the correlation matrix has 1s along its main diagonal (since the correlation of any variable with itself is 1). In R, the matrix of correlation is obtained by using the cor() function:

# Correlations
R <- cor(matr)
R
##             mpg        cyl       disp         hp
## mpg   1.0000000 -0.8521620 -0.8475514 -0.7761684
## cyl  -0.8521620  1.0000000  0.9020329  0.8324475
## disp -0.8475514  0.9020329  1.0000000  0.7909486
## hp   -0.7761684  0.8324475  0.7909486  1.0000000

This brings us to the main question of this post: Given a variance-covariance matrix, can we compute the correlation matrix? And can we do the reverse—moving from a correlation matrix back to a variance-covariance matrix? With reference to the matrices above, can we go from ‘Sigma’ to ‘R’ and from ‘R’ to ‘Sigma’? I often find myself needing to make these conversions, especially when working with mixed models or time series analysis.

-633.097208 / (sqrt(36.324103) * sqrt(15360.7998))
## [1] -0.8475514

# from covariance to correlation - 1
StDev <- sqrt(diag(Sigma))
StDevMat <- diag(StDev)
InvStDev <- solve(StDevMat)
InvStDev %*% Sigma %*% InvStDev
##            [,1]       [,2]       [,3]       [,4]
## [1,]  1.0000000 -0.8521620 -0.8475514 -0.7761684
## [2,] -0.8521620  1.0000000  0.9020329  0.8324475
## [3,] -0.8475514  0.9020329  1.0000000  0.7909486
## [4,] -0.7761684  0.8324475  0.7909486  1.0000000

# From covariance to correlation - 2
cov2cor(Sigma)
##             mpg        cyl       disp         hp
## mpg   1.0000000 -0.8521620 -0.8475514 -0.7761684
## cyl  -0.8521620  1.0000000  0.9020329  0.8324475
## disp -0.8475514  0.9020329  1.0000000  0.7909486
## hp   -0.7761684  0.8324475  0.7909486  1.0000000

# From covariance to correlation - 3
sweep(sweep(Sigma, 1, StDev, FUN = "/"), 2, StDev, FUN = "/")
##             mpg        cyl       disp         hp
## mpg   1.0000000 -0.8521620 -0.8475514 -0.7761684
## cyl  -0.8521620  1.0000000  0.9020329  0.8324475
## disp -0.8475514  0.9020329  1.0000000  0.7909486
## hp   -0.7761684  0.8324475  0.7909486  1.0000000

# From covariance to correlation - 4
scale(t(scale(t(Sigma), center = F, scale = StDev)), 
      center = F, scale = StDev)
##             mpg        cyl       disp         hp
## mpg   1.0000000 -0.8521620 -0.8475514 -0.7761684
## cyl  -0.8521620  1.0000000  0.9020329  0.8324475
## disp -0.8475514  0.9020329  1.0000000  0.7909486
## hp   -0.7761684  0.8324475  0.7909486  1.0000000
## attr(,"scaled:scale")
##        mpg        cyl       disp         hp 
##   6.026948   1.785922 123.938694  68.562868

# From covariance to correlation - 5
Sigma/outer(StDev, StDev)
##             mpg        cyl       disp         hp
## mpg   1.0000000 -0.8521620 -0.8475514 -0.7761684
## cyl  -0.8521620  1.0000000  0.9020329  0.8324475
## disp -0.8475514  0.9020329  1.0000000  0.7909486
## hp   -0.7761684  0.8324475  0.7909486  1.0000000

-0.8475514 * sqrt(36.324103) * sqrt(15360.7998)
## [1] -633.0972

# from correlation to covariance - 1
StDev <- sqrt(diag(Sigma))
StDevMat <- diag(StDev)
StDevMat %*% R %*% StDevMat
##             [,1]       [,2]       [,3]      [,4]
## [1,]   36.324103  -9.172379  -633.0972 -320.7321
## [2,]   -9.172379   3.189516   199.6603  101.9315
## [3,] -633.097208 199.660282 15360.7998 6721.1587
## [4,] -320.732056 101.931452  6721.1587 4700.8669

# From correlation to covariance - 2
lavaan::cor2cov(R, StDev)
##              mpg        cyl       disp        hp
## mpg    36.324103  -9.172379  -633.0972 -320.7321
## cyl    -9.172379   3.189516   199.6603  101.9315
## disp -633.097208 199.660282 15360.7998 6721.1587
## hp   -320.732056 101.931452  6721.1587 4700.8669

# From correlation to covariance - 3
sweep(sweep(R, 1, StDev, FUN = "*"), 2, StDev, FUN = "*")
##              mpg        cyl       disp        hp
## mpg    36.324103  -9.172379  -633.0972 -320.7321
## cyl    -9.172379   3.189516   199.6603  101.9315
## disp -633.097208 199.660282 15360.7998 6721.1587
## hp   -320.732056 101.931452  6721.1587 4700.8669

# From correlation to covariance - 5
scale(t(scale(t(R), center = F, scale = 1/StDev)), 
      center = F, scale = 1/StDev)
##              mpg        cyl       disp        hp
## mpg    36.324103  -9.172379  -633.0972 -320.7321
## cyl    -9.172379   3.189516   199.6603  101.9315
## disp -633.097208 199.660282 15360.7998 6721.1587
## hp   -320.732056 101.931452  6721.1587 4700.8669
## attr(,"scaled:scale")
##         mpg         cyl        disp          hp 
## 0.165921457 0.559934979 0.008068505 0.014585154

# From correlation to covariance - 6
R * outer(StDev, StDev)
##              mpg        cyl       disp        hp
## mpg    36.324103  -9.172379  -633.0972 -320.7321
## cyl    -9.172379   3.189516   199.6603  101.9315
## disp -633.097208 199.660282 15360.7998 6721.1587
## hp   -320.732056 101.931452  6721.1587 4700.8669

I recently updated nycdogs, a data package for R that I originally made for teaching purposes. It bundles up New York City’s dog licensing and dog bite datasets into an R package. The data is already provided by the city in a directly usable way (as CSV files), but packaging it up makes it easier to integrate with other packages and get things moving a little more quickly in class.

The nyc_license table looks like this:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15

nyc_license
#> # A tibble: 819,323 × 12
#>    animal_name animal_gender animal_birth_year breed_name       breed_rc zip_code zip   license_issued_date license_expired_date extract_year borough city 
#>    <chr>       <chr>                     <int> <chr>            <chr>    <chr>    <chr> <date>              <date>                      <int> <chr>   <chr>
#>  1 Paige       F                          2014 American Pit Bu… Pit Bul… 10035    10035 2014-09-12          2017-09-12                   2016 Manhat… New …
#>  2 Yogi        M                          2010 Boxer            Boxer    10465    10465 2014-09-12          2017-10-02                   2016 Bronx   Bronx
#>  3 Ali         M                          2014 Basenji          Basenji  10013    10013 2014-09-12          2019-09-12                   2016 Manhat… New …
#>  4 Queen       F                          2013 Akita Crossbreed Akita C… 10013    10013 2014-09-12          2017-09-12                   2016 Manhat… New …
#>  5 Lola        F                          2009 Maltese          Maltese  10028    10028 2014-09-12          2017-10-09                   2016 Manhat… New …
#>  6 Ian         M                          2006 Unknown          Unknown  10013    10013 2014-09-12          2019-10-30                   2016 Manhat… New …
#>  7 Buddy       M                          2008 Unknown          Unknown  10025    10025 2014-09-12          2017-10-20                   2016 Manhat… New …
#>  8 Chewbacca   F                          2012 Labrador Retrie… Labrado… 10013    10013 2014-09-12          2019-10-01                   2016 Manhat… New …
#>  9 Heidi-bo    F                          2007 Dachshund Smoot… Dachshu… 11215    11215 2014-09-13          2017-04-16                   2016 Brookl… Broo…
#> 10 Massimo     M                          2009 Bull Dog, French French … 11201    11201 2014-09-13          2017-09-17                   2016 Brookl… Broo…
#> # ℹ 819,313 more rows

The dataset is good for teaching with for several reasons. First, people like dogs. Second, it has about 800,000 records, which is a lot in many ways but also trivial for a modern laptop to work with. Third, the data are pretty clean, but they’re not that clean. There are various things going on in there that invite us to explore it but also to be careful about what we can and can’t be confident about.

The first thing we ask about any data set like this is: what is each row in this table? In the case of nyc_license, each row is a record of one time-limited license for a dog. A row is not necessarily the record of a unique dog. When an owner renews a dog’s license, the dog will appear one time for each license period. In New York City, dog licenses can be issued for one year or five years. The table does not provide any dog_id column that ties sequences of licenses to individual dogs.

Second, the table contains several “extracts” of the license data. It seems like the way the provider updates the NYC Open Data record is to periodically pull from its internal database every year or so and append it to the existing, previously-released version. The extract_year column identifies the extract, so at least we know that. But it’s also the case that a license that was active in more than one extract year appears in each extract. As a result, the table has a large number of duplicate rows.

Third, license information is provided by dog owners. This means that there’s a fair amount of noise in the records, because people make typos and other errors. For instance, the ZIP code data has incomplete ZIPs. Some people provide invalid ZIP codes (e.g. ones that are fewer than five digits long). Others (about 3,300 records) have ZIP codes that are not in New York City. Eighty percent of these (about 2,700) are ZIP codes located elsewhere in New York state. Another 77 are located in Florida, which I kind of like. Five of them live in Boca Raton, for example. Their kids never call. You’d think it would be easy for them to just pick up the bone, I mean phone, but you never hear from them.

About 650 rows have valid ZIP codes that are in Manhattan, but these ZIPs are not ZCTAs. The trouble with ZIP codes is that people think of them as areas, as polygons, but in fact they are lines. They are defined by the Post Office as postal routes, i.e. as journeys a mail carrier takes. While they tend to form closed polygons, they are absolutely not guaranteed to. But people know their ZIP code and all kinds of analysts (from marketers to political operatives) like thinking in terms of ZIP codes. So the Census defines ZIP Code Tabulation Areas that are a kind of best-effort at answering the question “What if this ZIP code were a closed polygon?” ZCTAs exclude every P.O. box and single-address ZIP nationwide. Manhattan has several buildings that have their own ZIP code. New York City also defines MODZCTAs that merge some low- or no-population ZCTAs with their neighbors. I use these below.

Fourth, beyond the ZIP codes, owners also provide information about the age and sex of their dogs, the dog’s name, and the dog’s breed. Again, this kind of information tends to be a little noisy. For example, people vary in how they spell names. They also vary in how they classify their dog by breed. Beyond the list of “recognized” breeds (itself something that dog owners variably conform to), some people might describe their dog as, e.g., “Labrador-Sheltie mix” while others might just say “Unknown”. The dataset includes the breed_name column from the original dataset and also breed_rc, which is a recoding by me to clean up equivalent terms (“Yorkshire Terrier” and “Terrier, Yorkshire”, for example) along with a few other coding decisions.

Finally, and only because I’m choosing to stop here, the table just shows the number of licensed dogs in New York City. That is definitely not the same as the number of dogs there. Moreover, license compliance varies quite a bit in ways that are related to characteristics of the owner that we can’t observe, and also ones that we partially can, like where they live.

All of which is related to the fact that one of the lolsob-mantras of courses I teach is “Now you know the true meaning of data”. In the graphs and maps below, I’ve done the best I can, where appropriate, to uniquely identify individual dogs in the dataset (even though ultimately we cannot).

Popular Breeds

With that out of the way, it’s time to indulge your prejudices about Dogs, New Yorkers, and New York Dogs. The most common “breed” of dog in New York is “Unknown”, as befits a city of immigrants, hustlers, scrappy up-and-comers, and dogs seeking to leave their small-town past behind and remake their identities. Somewhere between eight and ten percent of records each year have “Unknown” recorded as their breed. There’s some change over time (a declining share), but I think this is an artifact of how breeds are getting recorded rather than a change in the composition of the dog population. If we leave the Unknowns out, here are the trends in annual licenses issued to the twenty most popular breeds.

Some breeds are pretty steadily popular. Smaller breeds do well in a city full of apartments. We can also see some evidence here of trends in popularity. However, these annual numbers are what we might call the prevalence of dogs, rather than the incidence of new dogs. If we want to get a better sense of fashion trends in dog breeds, we should look at puppies. (Who doesn’t want to look at puppies?) Because of the structure of the data we can’t quite do this. The best we can do, based on the dog’s age, is to examine licenses issued to dogs younger than a year old. Let’s pick out some breeds that seem to be rising, a few that seem relatively steady, and a few that might be declining. I’ve indexed these changes relative to each breed in 2016 by setting the value for that year to 100. Even though our range is narrow, we use a log scale on the y-axis because, e.g., a rise in a breed’s index from 100 to 200 reflects the same relative change as a fall from 100 to 50.

Maltipoos and Goldendoodles seem to be doing pretty well for themselves. Golden Retrievers remain steadily popular. Meanwhile, Shiba Inus, beloved fad-breed of the 2010s, seem to have been affected by their association with crypto and are strongly in decline. Much drop. Very decline. Wow. Beagles are also doing poorly. Good grief.

Where Breeds Live

We can look at the relative prevalence of breeds by their ZIP code (strictly speaking, by their MODZCTA). I hesitate a little on this just because, even with 800,000 licenses, we’ll run out of degrees of freedom fairly quickly if we break things down by breed and ZIP. Relatedly, we should ask what sort of metric we’re trying to show here. What’s the denominator? (A surprising amount of Statistics, I claim, is essentially exploring the question “What happens if I divide by this number instead of that number?”) Popular dogs are, to a first approximation, pretty popular all over the city. Most ZIP codes have a solid contingent of Yorkies. If, for each popular breed, we calculate the share of that breed as a percentage of all dogs in the ZIP code, we’ll see some ZIP-by-ZIP variation come out.

It would also be nice to adjust for the fact that different ZIP codes have quite different shares of the total number of dogs licensed in the city. (Which, again, is not the same as share of the total number of dogs in the city.) So we calculate the within-ZIP share and make a ratio of it by dividing that number by the breed’s citywide share. The ratio is 1 when the breed is just as common in that ZIP code as it is across the whole city. If it’s 2, the breed is twice as common in that ZIP as it is across the city. If it’s 0.5, the breed is half as common in that ZIP as it is in the city as a whole. The main downside of ratios like this is that they can exaggerate differences a little, especially if the base rates are small. We have enough breeds in each zip code to make me feel OK about doing this. But if you refer to the bar charts above you can see that even for disproportionately popular breeds in zips, you’re looking at changes like “twice as common with respect to a baseline of five percent of all dogs”. Here’s the map:

We see some plausible evidence that breeds are distributed across ZIP codes in a manner that tracks the demography of owners. French Bulldogs in Chelsea, Soho, and TriBeCa. Havanese on the Upper East Side. Pit Bulls in the Bronx. German Shepherds on Staten Island. Goldendoodles and King Charles Spaniels in Manhattan south of 96th.

The Diffusion of Goldendoodles

Let’s take a closer look at the Goldendoodles. In 2015, Goldendoodles weren’t that common. We can ask: where were they initially popular? Did they diffuse across New York City from focal points where people were more likely to set trends in breed fashion? Restricting ourselves to young dogs, we can try to look at each ZIP’s share of all Goldendoodle puppies licensed each year from 2015 to 2025. Here’s what that looks like:

In 2015, Goldendoodle puppies are licensed in a relatively small number of ZIP codes, mostly in Manhattan. Most ZIP codes are doodle-free. Over the next five years, they gradually spread across the city. By 2023 or 2024, most—though still not all—ZIP codes are licensing young Goldendoodles.

Again, the underlying numbers aren’t huge. So take all of this with a delicious treat. I mean, grain of salt. I also wonder what the next few years are going to look like in terms of dog popularity. Based on the maps above, I think Australian Shepherds are likely to continue to rise in popularity amongst owners looking to acquire a fashionable breed. And while they’re not on any of the maps above, the data also suggest that the Jindo, a Korean breed that is a relatively recent arrival in the US, has been sharply on the rise. If you ask me, though, if you’re in New York City and planning to get a dog, you should throw in your lot with Unknown and get a good, solid, dog-shaped mutt. Preferably one with strong views on the state of the Subway these days, or who knows a guy who can get you a deal on some shirts, or the special ramen, or tickets. Not tickets for anything in particular. Just, like, tickets. Works out of an upstairs bar downtown on Water St. Trust me.

Data centers required to fuel the ambitions of Artificial Intelligence companies have sparked a nationwide conversation. AI companies are moving quickly to build out infrastructure and local communities are pushing back to resist what they see as the appropriation of public resources such as electricity and water to further enrich wealthy corporations. Moreover, the sentiment shared by both the political left and right is that data centers will significantly degrade the quality of life for those living adjacent to a data center and to the extended community as well. This post implements a non-trivial Bayesian Network to study the risks involved. I hope it is of some interest to a broad audience including, people interested in clarifying their thoughts about the risks that data centers may pose to their particular community and analysts curious about Bayesian Networks. However, I think the sample model built out in this post will be especially interesting to decision analysts who know quite a bit about Bayesian Networks but may be unfamiliar with the R programming language and unaware that there is capable opensource software which can provide capabilities beyond those offered by proprietary software systems.

Other than addressing the topic of data center risks the analysis contains no original research except perhaps the contribution of ChatGPT 5.6 Luna which conjured up the innovative parts of the synthetic data generation algorithm from a series of requirements.

%%| label: fig-risk-network
%%| fig-cap: "Data Center Risk Network"
%%| fig-alt: "Directed network showing data-center impacts and risks leading to overall community risk."



flowchart LR

  Scale --> LandFootprint
  Scale --> ElectricityDemand
  Scale --> WaterDemand
  Scale --> ConstructionIntensity
  Scale --> BackupGeneration
  Scale --> WorkforceDemand
  Scale --> NoiseExposure
  Scale --> AirQualityExposure
  Scale --> AccessToBenefits
  Scale --> EconomicBenefit

  Cooling --> WaterDemand
  Cooling --> NoiseExposure

  Distance --> NoiseExposure
  Distance --> AirQualityExposure
  Distance --> TrafficExposure
  Distance --> VisualLandUseExposure
  Distance --> AccessToBenefits
  Distance --> ElectricityCostRisk

  LandFootprint --> VisualLandUseExposure
  ConstructionIntensity --> TrafficExposure
  BackupGeneration --> AirQualityExposure

  ElectricityDemand --> GridStress
  GridCapacity --> GridStress
  GridStress --> ElectricityCostRisk
  GridStress --> ReliabilityRisk

  WaterDemand --> WaterCompetition
  WaterAvailability --> WaterCompetition
  WaterCompetition --> WaterSecurityRisk
  Vulnerability --> WaterSecurityRisk

  WorkforceDemand --> HousingPressure
  HousingAvailability --> HousingPressure
  WorkforceDemand --> ServicePressure
  HousingAvailability --> ServicePressure

  LandFootprint --> EnvironmentalDegradation
  WaterDemand --> EnvironmentalDegradation
  WaterAvailability --> EnvironmentalDegradation

  NoiseExposure --> QualityOfLifeRisk
  AirQualityExposure --> QualityOfLifeRisk
  TrafficExposure --> QualityOfLifeRisk
  VisualLandUseExposure --> QualityOfLifeRisk
  HousingPressure --> QualityOfLifeRisk
  Vulnerability --> QualityOfLifeRisk

  AccessToBenefits --> EconomicBenefit

  QualityOfLifeRisk --> OverallCommunityRisk
  WaterSecurityRisk --> OverallCommunityRisk
  ElectricityCostRisk --> OverallCommunityRisk
  ReliabilityRisk --> OverallCommunityRisk
  EnvironmentalDegradation --> OverallCommunityRisk
  HousingPressure --> OverallCommunityRisk
  ServicePressure --> OverallCommunityRisk
  EconomicBenefit --> OverallCommunityRisk

# ============================================================
# Bayesian Network for Data-Center Community Risk
# ============================================================

# Install packages once if necessary:
# install.packages(c("bnlearn", "ggplot2"))

library(bnlearn)
library(ggplot2)
library(Rgraphviz)
library(igraph)
library(dplyr)
library(tidyr)
library(ggridges)
library(ggtern)

# ------------------------------------------------------------
# 2. Helper functions for generating categorical variables
# ------------------------------------------------------------

# sample_category <- function(probabilities, labels = names(probabilities)) {
#   sample(labels, size = 1, prob = probabilities)
# }

# Convert an ordinal category to a numeric score.
# Useful for simple synthetic-data rules.
ordinal_score <- function(x) {
  match(x, c("Low", "Moderate", "High", "VeryHigh"))
}

set.seed(20260914)

# ------------------------------------------------------------
# 1. Define the expert network structure
# ------------------------------------------------------------

# Short node names are used because they are easier to reference
# in bnlearn expressions and plots.

network_string <- paste0(
  "[Distance]",
  "[Scale]",
  "[Cooling]",
  "[GridCapacity]",
  "[WaterAvailability]",
  "[HousingAvailability]",
  "[Vulnerability]",
  "[LandFootprint|Scale]",
  "[ElectricityDemand|Scale]",
  "[WaterDemand|Scale:Cooling]",
  "[ConstructionIntensity|Scale]",
  "[BackupGeneration|Scale]",
  "[WorkforceDemand|Scale]",
  "[NoiseExposure|Distance:Cooling:Scale]",
  "[AirQualityExposure|Distance:BackupGeneration:Scale]",
  "[TrafficExposure|Distance:ConstructionIntensity]",
  "[VisualLandUseExposure|Distance:LandFootprint]",
  "[AccessToBenefits|Distance:Scale]",
  "[GridStress|ElectricityDemand:GridCapacity]",
  "[ElectricityCostRisk|GridStress:Distance]",
  "[ReliabilityRisk|GridStress]",
  "[WaterCompetition|WaterDemand:WaterAvailability]",
  "[WaterSecurityRisk|WaterCompetition:Vulnerability]",
  "[HousingPressure|WorkforceDemand:HousingAvailability]",
  "[ServicePressure|WorkforceDemand:HousingAvailability]",
  "[EnvironmentalDegradation|LandFootprint:WaterDemand:WaterAvailability]",
  "[QualityOfLifeRisk|NoiseExposure:AirQualityExposure:TrafficExposure:",
  "VisualLandUseExposure:HousingPressure:Vulnerability]",
  "[EconomicBenefit|Scale:AccessToBenefits]",
  "[OverallCommunityRisk|QualityOfLifeRisk:WaterSecurityRisk:",
  "ElectricityCostRisk:ReliabilityRisk:EnvironmentalDegradation:",
  "HousingPressure:ServicePressure:EconomicBenefit]"
)

dag <- model2network(network_string)

# Optional inspection:
#print(dag)
#print(modelstring(dag))

#list(nodes(dag))
#The `arcs()` function explicitly lists the from/to structure of the network.

arcs_dag <- as.data.frame(arcs(dag))
#rownames(arcs_dag) <- NULL
#print(arcs_dag[order(arcs_dag$from),],row.names=FALSE)

# Replace each original node name with its alias.
# The lookaround expressions ensure that only complete node names are replaced.

# Alias mapping: original node name = abbreviated node name
node_alias <- c(
  Distance                  = "dst",
  Scale                     = "scl",
  Cooling                   = "clg",
  GridCapacity              = "grd",
  WaterAvailability         = "wat",
  HousingAvailability       = "hou",
  Vulnerability              = "vul",
  LandFootprint              = "lnd",
  ElectricityDemand         = "eld",
  WaterDemand               = "wdr",
  ConstructionIntensity     = "con",
  BackupGeneration          = "gen",
  WorkforceDemand           = "wrk",
  NoiseExposure             = "noi",
  AirQualityExposure        = "air",
  TrafficExposure           = "trf",
  VisualLandUseExposure     = "vis",
  AccessToBenefits          = "ben",
  GridStress                = "gst",
  ElectricityCostRisk       = "ecr",
  ReliabilityRisk           = "rel",
  WaterCompetition          = "wcp",
  WaterSecurityRisk         = "wsr",
  HousingPressure           = "hpr",
  ServicePressure           = "spr",
  EnvironmentalDegradation  = "env",
  QualityOfLifeRisk         = "qol",
  EconomicBenefit            = "eco",
  OverallCommunityRisk      = "ocr"
)
network_string_alias <- network_string

for (node_name in names(node_alias)) {
  alias_name <- unname(node_alias[node_name])

  pattern <- paste0(
    "(?<![A-Za-z0-9_.])",
    node_name,
    "(?![A-Za-z0-9_.])"
  )

  network_string_alias <- gsub(
    pattern,
    alias_name,
    network_string_alias,
    perl = TRUE
  )
}

# Create the original DAG
dag <- model2network(network_string)

# Create the DAG using aliases
dag_alias <- model2network(network_string_alias)

# View the aliased network string
#network_string_alias

# View the aliased DAG
#dag_alias

graphviz.plot(dag_alias,
              shape = "ellipse",
              main = "Data-Center Community-Risk Bayesian Network")

#library(bnlearn)

# dag should be a bnlearn network created with model2network(),
# empty.graph(), hc(), tabu(), gs(), etc.

node_names <- nodes(dag)

# Markov blanket size for each node
markov_blanket_sizes <- sapply(
  node_names,
  function(node) {
    length(mb(dag, node))
  }
)

# Direct neighbours = parents + children
neighbourhood_sizes <- sapply(
  node_names,
  function(node) {
    direct_neighbours <- unique(c(
      parents(dag, node),
      children(dag, node)
    ))
    
    length(direct_neighbours)
  }
)

# Branching factor = number of children
branching_factors <- sapply(
  node_names,
  function(node) {
    length(children(dag, node))
  }
)

# Network-level averages
network_statistics <- data.frame(
  Statistic = c(
    "Average Markov blanket size",
    "Average neighbourhood size",
    "Average branching factor"
  ),
  Value = c(
    mean(markov_blanket_sizes),
    mean(neighbourhood_sizes),
    mean(branching_factors)
  )
)

network_statistics$Value <- round(
  network_statistics$Value,
  digits = 2
)

network_statistics

                    Statistic Value
1 Average Markov blanket size  7.31
2  Average neighbourhood size  3.52
3    Average branching factor  1.76

node_statistics <- data.frame(
  Node = node_names,
  MB_Size = as.numeric(markov_blanket_sizes),
  Neigh_Size = as.numeric(neighbourhood_sizes),
  BF_Size = as.numeric(branching_factors),
  row.names = NULL
)

#node_statistics

# ================================================================
# Synthetic data generation from a bnlearn DAG
# ================================================================

#library(bnlearn)

# ----------------------------------------------------------------
# Main function
# ----------------------------------------------------------------
generate_synthetic_data <- function(
  n = 10000, 
  global_risk_shift = 0.00, 
  parent_effect_strength = 1.00, 
  conditional_noise = 0.10
) {
  # ARCHITECTURAL NOTE: Dynamic environment lookup lets the function bypass bnlearn's 
  # restrictive function signatures while maintaining access to the global network structure.
  if (!exists("dag", envir = parent.frame(), inherits = TRUE)) {
    stop("An object named 'dag' must exist before calling generate_synthetic_data().")
  }
  dag <- get("dag", envir = parent.frame(), inherits = TRUE)
  
  # Structural defensive validation steps
  validate_dag(dag)
  validate_generator_arguments(n, global_risk_shift, parent_effect_strength, conditional_noise)
  
  # Bundle configurations into a centralized configuration object
  config <- create_network_config(
    n = n, 
    global_risk_shift = global_risk_shift, 
    parent_effect_strength = parent_effect_strength, 
    conditional_noise = conditional_noise
  )
  
  # Execute the core simulation pipeline
  generate_rows_from_dag(dag = dag, config = config)
}

# ----------------------------------------------------------------
# Configuration
# ----------------------------------------------------------------
create_network_config <- function(n, global_risk_shift, parent_effect_strength, conditional_noise, seed = 20260921) {
  list(
    n = n,
    seed = seed,
    global_risk_shift = global_risk_shift,
    parent_effect_strength = parent_effect_strength,
    conditional_noise = conditional_noise,
    states = create_node_states(),
    root_probabilities = create_root_probabilities(),
    
    # PARENT WEIGHTING MATRIX: Allows specific upstream nodes to exert a stronger 
    # or weaker relative pull on their children during the weighted average calculation.
    parent_weights = c(
      Distance = 1.20, Scale = 1.00, Cooling = 1.00, GridCapacity = 1.00, 
      HousingAvailability = 1.00, Vulnerability = 1.10, WaterAvailability = 1.00
    ),
    risk_mappings = create_risk_mappings()
  )
}

create_node_states <- function() {
  # Definition of ordinal categorical states for every node in the entire network universe
  list(
    Distance = c("Adjacent", "Neighborhood", "Municipality", "Watershed", "Regional", "Global"),
    Scale = c("Small", "Medium", "Large", "Hyperscale"),
    Cooling = c("AirCooled", "Hybrid", "Evaporative", "WaterIntensive"),
    GridCapacity = c("Adequate", "Constrained", "SeverelyConstrained"),
    HousingAvailability = c("Adequate", "Tight", "SeverelyConstrained"),
    Vulnerability = c("Low", "Moderate", "High"),
    WaterAvailability = c("Abundant", "SeasonalStress", "Scarce"),
    AccessToBenefits = c("Low", "Moderate", "High"),
    BackupGeneration = c("Low", "Moderate", "Extensive"),
    ConstructionIntensity = c("Low", "Moderate", "High"),
    ElectricityDemand = c("Low", "Moderate", "High", "Extreme"),
    LandFootprint = c("Low", "Moderate", "High"),
    NoiseExposure = c("Low", "Moderate", "High"),
    WaterDemand = c("Low", "Moderate", "High"),
    WorkforceDemand = c("Low", "Moderate", "High"),
    AirQualityExposure = c("Low", "Moderate", "High"),
    EconomicBenefit = c("Low", "Moderate", "High"),
    EnvironmentalDegradation = c("Low", "Moderate", "High"),
    GridStress = c("Low", "Moderate", "High"),
    HousingPressure = c("Low", "Moderate", "High"),
    ServicePressure = c("Low", "Moderate", "High"),
    TrafficExposure = c("Low", "Moderate", "High"),
    VisualLandUseExposure = c("Low", "Moderate", "High"),
    WaterCompetition = c("Low", "Moderate", "High"),
    ElectricityCostRisk = c("Low", "Moderate", "High"),
    QualityOfLifeRisk = c("Low", "Moderate", "High"),
    ReliabilityRisk = c("Low", "Moderate", "High"),
    WaterSecurityRisk = c("Low", "Moderate", "High"),
    OverallCommunityRisk = c("Low", "Moderate", "High")
  )
}

# create_root_probabilities sets the probability distributions for how frequently these sates occur.
create_root_probabilities <- function() {
  # Baseline categorical probability distributions applied strictly to exogenous Root Nodes (no parents)
  list(
    Distance = c(Adjacent = 0.10, Neighborhood = 0.18, 
                 Municipality = 0.25, Watershed = 0.22, Regional = 0.18, Global = 0.07),
    Scale = c(Small = 0.15, Medium = 0.35, Large = 0.35, Hyperscale = 0.15),
    Cooling = c(AirCooled = 0.30, Hybrid = 0.30, Evaporative = 0.25, WaterIntensive = 0.15),
    GridCapacity = c(Adequate = 0.50, Constrained = 0.35, SeverelyConstrained = 0.15),
    HousingAvailability = c(Adequate = 0.45, Tight = 0.35, SeverelyConstrained = 0.20),
    Vulnerability = c(Low = 0.45, Moderate = 0.35, High = 0.20),
    WaterAvailability = c(Abundant = 0.45, SeasonalStress = 0.35, Scarce = 0.20)
  )
}

# risk mappings are not probabilities. They do not sum to 1.0. They quantify the downstream impact or "hazard level" of a state once it has been selected. 
create_risk_mappings <- function() {
  # DIRECT NODES: Hand-calibrated risk translations mapping structural categories to numeric impacts.
  mappings <- list(
    Distance = c(Adjacent = 1.00, Neighborhood = 0.80, Municipality = 0.45, Watershed = 0.25, Regional = 0.10, Global = 0.02),
    Scale = c(Small = 0.25, Medium = 0.50, Large = 0.75, Hyperscale = 1.00),
    Cooling = c(AirCooled = 0.25, Hybrid = 0.55, Evaporative = 0.75, WaterIntensive = 0.95),
    GridCapacity = c(Adequate = 0.10, Constrained = 0.55, SeverelyConstrained = 0.95),
    HousingAvailability = c(Adequate = 0.10, Tight = 0.60, SeverelyConstrained = 0.95),
    Vulnerability = c(Low = 0.10, Moderate = 0.55, High = 0.95),
    WaterAvailability = c(Abundant = 0.10, SeasonalStress = 0.60, Scarce = 0.95),
    LandFootprint = c(Low = 0.25, Moderate = 0.60, High = 1.00),
    ElectricityDemand = c(Low = 0.15, Moderate = 0.40, High = 0.70, Extreme = 1.00),
    WaterDemand = c(Low = 0.15, Moderate = 0.55, High = 1.00),
    ConstructionIntensity = c(Low = 0.25, Moderate = 0.60, High = 1.00),
    BackupGeneration = c(Low = 0.25, Moderate = 0.60, Extensive = 1.00),
    WorkforceDemand = c(Low = 0.15, Moderate = 0.55, High = 1.00),
    
    # PROTECTIVE VARIABLE CONVERSION: Inverts risk mapping so high benefits lower the downstream risk score
    AccessToBenefits = c(Low = 0.85, Moderate = 0.50, High = 0.15),
    EconomicBenefit = c(Low = 0.85, Moderate = 0.50, High = 0.15)
  )
  
  # DERIVED NODES: Programmatic batch generation of standardized downstream risk states
  standard_risk_nodes <- c(
    "NoiseExposure", "AirQualityExposure", "TrafficExposure", "VisualLandUseExposure", 
    "EnvironmentalDegradation", "GridStress", "HousingPressure", "ServicePressure", 
    "WaterCompetition", "ElectricityCostRisk", "QualityOfLifeRisk", "ReliabilityRisk", 
    "WaterSecurityRisk", "OverallCommunityRisk"
  )
  standard_mapping <- c(Low = 0.15, Moderate = 0.50, High = 0.85)
  for (node in standard_risk_nodes) {
    mappings[[node]] <- standard_mapping
  }
  mappings
}

# ----------------------------------------------------------------
# Generate all observations
# ----------------------------------------------------------------
generate_rows_from_dag <- function(dag, config) {
  if (!is.null(config$seed)) { set.seed(config$seed) }
  
  dag_nodes <- bnlearn::nodes(dag)
  
  # CRITICAL ARCHITECTURAL STEP: Nodes must be sorted topologically. 
  # This guarantees a parent node has simulated data before a child tries to read it.
  generation_order <- get_topological_order(dag)
  
  validate_configuration(dag = dag, dag_nodes = dag_nodes, config = config)
  
  # OPTIMIZATION: Memory pre-allocation prevents R from copying data frames repeatedly in memory loops.
  result <- setNames(vector("list", length(dag_nodes)), dag_nodes)
  for (node in dag_nodes) {
    result[[node]] <- character(config$n)
  }
  
  # Row-by-row simulation loop
  for (i in seq_len(config$n)) {
    one_row <- generate_one_observation(dag = dag, generation_order = generation_order, config = config)
    for (node in dag_nodes) {
      result[[node]][i] <- one_row[[node]]
    }
  }
  
  # Formatting raw arrays back into structured, clean R Factor columns
  result <- as.data.frame(result, stringsAsFactors = FALSE, check.names = FALSE)
  for (node in dag_nodes) {
    result[[node]] <- factor(result[[node]], levels = config$states[[node]])
  }
  result
}

generate_one_observation <- function(dag, generation_order, config) {
  row <- setNames(vector("list", length(generation_order)), generation_order)
  
  # Sequentially process each node down the causal chain
  for (node in generation_order) {
    node_states <- config$states[[node]]
    node_parents <- bnlearn::parents(dag, node)
    
    # ROOT NODE INTERCEPT: If a node has no parents, sample from its static baseline distribution.
    if (length(node_parents) == 0) {
      probabilities <- config$root_probabilities[[node]]
      probabilities <- probabilities[node_states]
      probabilities <- normalize_probabilities(probabilities)
      names(probabilities) <- node_states
      row[[node]] <- sample_state(probabilities)
      next
    }
    
    # CHILD NODE COMPUTATION: Compress multi-parent states into a singular continuous score
    parent_risk <- calculate_parent_risk(parents = node_parents, row = row, config = config)
    
    # LINEAR SENSITIVITY TRANSFORM: Shifts values relative to a 0.5 mid-point based on network effects strength
    conditional_risk <- 0.5 + config$parent_effect_strength * (parent_risk - 0.5)
    
    # Map the scalar risk target back into a probability distribution across discrete states
    probabilities <- calculate_state_probabilities(
      node = node, states = node_states, 
      conditional_risk = conditional_risk, config = config
    )
    row[[node]] <- sample_state(probabilities)
  }
  row
}

# ----------------------------------------------------------------
# Parent-risk calculations

## ----------------------------------------------------------------
calculate_parent_risk <- function(parents, row, config) {
## Extract configured or default (1.00) relative weights for each parent
parent_weights <- vapply(parents, function(parent) {
get_parent_weight(parent = parent, config = config)
}, numeric(1))
## Convert each parent's current simulated text state into its corresponding numeric risk score
weighted_scores <- vapply(seq_along(parents), function(i) {
parent <- parents[[i]]
parent_score <- state_risk_score(node = parent, state = row[[parent]], config = config)
parent_score * parent_weights[[i]]
}, numeric(1))
## COMPRESSION METRIC: Calculated as a weighted average to safely scale multi-parent risks between 0 and 1
sum(weighted_scores) / sum(parent_weights)
}
get_parent_weight <- function(parent, config) {
weights <- config$parent_weights
## Defensive safety check guarding against structural mismatches between the input DAG and configuration vectors
if (is.null(weights) || is.null(names(weights)) || !parent %in% names(weights)) {
return(1.00)
}
weight <- unname(weights[[parent]])
if (!is.finite(weight) || weight < 0) {
stop("The weight for parent '", parent, "' must be a finite, non-negative number.")
}
weight
}
## ----------------------------------------------------------------## Conditional probability calculations## ----------------------------------------------------------------
calculate_state_probabilities <- function(node, states, conditional_risk, config) {
## STOCHASTIC PERTURBATION: Incorporates systemic shifts (scenarios) and random environmental noise via Gaussian distributions.
target_risk <- conditional_risk + config$global_risk_shift + rnorm(n = 1, mean = 0, sd = config$conditional_noise)
## Hard boundaries preventing mathematically invalid probabilities outside (0, 1)
target_risk <- max(0.01, min(target_risk, 0.99))
state_scores <- vapply(states, function(state) {
state_risk_score(node = node, state = state, config = config)
}, numeric(1) )
## EXPONENTIAL DISTANCE KERNEL (RBF-style): States whose static risk mappings sit closest## to the dynamic target_risk calculation receive exponentially higher generation probabilities.
concentration <- 6
probabilities <- exp(-concentration * abs(state_scores - target_risk))
## Convert distances to a standard, valid probability array summing to 1.00
probabilities <- normalize_probabilities(probabilities)
names(probabilities) <- states
probabilities
}
state_risk_score <- function(node, state, config) {
mapping <- config$risk_mappings[[node]]
if (is.null(mapping)) { stop("No risk mapping is defined for node '", node, "'.") }
if (!state %in% names(mapping)) { stop("No risk score is defined for node '", node, "' and state '", state, "'.") }
unname(mapping[[state]])
}
normalize_probabilities <- function(probabilities) {
probabilities <- pmax(probabilities, 0)
if (length(probabilities) == 0) { stop("Cannot normalize an empty probability vector.") }
## Uniform fallback handling to preserve math stability if all states evaluate to a zero score
if (sum(probabilities) == 0) {
probabilities <- rep(1 / length(probabilities), length(probabilities))
} else {
probabilities <- probabilities / sum(probabilities)
}
probabilities
}
sample_state <- function(probabilities) {
## Executes a standard weighted random draw based on the dynamically formulated distribution
sample(x = names(probabilities), size = 1, prob = probabilities)
}
## ----------------------------------------------------------------## Validation & Topological ordering## ----------------------------------------------------------------
validate_dag <- function(dag) {
if (!inherits(dag, "bn")) { stop("dag must be a bnlearn network object with class 'bn'.") }
}
validate_generator_arguments <- function(n, global_risk_shift, parent_effect_strength, conditional_noise) {
if (length(n) != 1 || !is.numeric(n) || is.na(n) || n < 1 || n != as.integer(n)) { stop("n must be one positive integer.") }
numeric_values <- c(global_risk_shift, parent_effect_strength, conditional_noise)
if (any(!is.numeric(numeric_values)) || any(!is.finite(numeric_values))) { stop("global_risk_shift, parent_effect_strength, and conditional_noise must be finite numeric values.") }
if (conditional_noise < 0) { stop("conditional_noise must be non-negative.") }
}
validate_configuration <- function(dag, dag_nodes, config) {
missing_states <- setdiff(dag_nodes, names(config$states))
if (length(missing_states) > 0) { stop("No state definitions supplied for: ", paste(missing_states, collapse = ", ")) }
root_nodes <- dag_nodes[vapply(dag_nodes, function(node) { length(bnlearn::parents(dag, node)) == 0 }, logical(1))]
missing_root_probabilities <- setdiff(root_nodes, names(config$root_probabilities))
if (length(missing_root_probabilities) > 0) { stop("No root probabilities supplied for: ", paste(missing_root_probabilities, collapse = ", ")) }
missing_mappings <- setdiff(dag_nodes, names(config$risk_mappings))
if (length(missing_mappings) > 0) { stop("No risk mappings supplied for: ", paste(missing_mappings, collapse = ", ")) }
}
get_topological_order <- function(dag) {
## KAHN'S ALGORITHM EXTENSION: Resolves dependence hierarchies.## Iteratively extracts nodes whose parents have already crossed the resolution boundary.
dag_nodes <- bnlearn::nodes(dag)
remaining_nodes <- dag_nodes
ordering <- character(0)
while (length(remaining_nodes) > 0) {
available_nodes <- remaining_nodes[vapply(remaining_nodes, function(node) {
node_parents <- bnlearn::parents(dag, node)
all(node_parents %in% ordering)
}, logical(1))]
# If no nodes are free but un-generated items remain, a cyclic loop exists.
if (length(available_nodes) == 0) {
stop("The graph does not appear to be acyclic, or its parent relationships could not be resolved.")
}
available_nodes <- sort(available_nodes)
ordering <- c(ordering, available_nodes)
remaining_nodes <- setdiff(remaining_nodes, available_nodes)
}
ordering
}

# ================================================================
# Usage
# ================================================================
#
# A DAG named `dag` must exist first.

# Global scenario controls.
  # global_risk_shift
  # 0.00 = baseline moderate-risk scenario
  # Positive values increase risk throughout the network
  # Negative values decrease risk throughout the network
  
  #  parent_effect_strength
  # Controls how strongly parent risk affects child risk.
  # 1.00 is the default.
 
  # conditional_noise 
  # Controls random variation in conditional probabilities.
  # Larger values produce less deterministic relationships.

#
synthetic <- generate_synthetic_data(
  n = 5000,
  global_risk_shift = 0.00,
  parent_effect_strength = 1.0,
  conditional_noise = 0.10
)
#
# str(synthetic)
# head(synthetic)
# ================================================================

state_proportions <- function(data) {
  lapply(data, function(x) {
    round(prop.table(table(x)), 3)
  })
}

proportions <- state_proportions(synthetic)
proportions$Distance

x
    Adjacent Neighborhood Municipality    Watershed     Regional       Global 
       0.100        0.166        0.253        0.231        0.179        0.071 

data_summary_table <- do.call(
  rbind,
  lapply(names(synthetic), function(node) {
    counts <- table(synthetic[[node]])
    proportions <- prop.table(counts)

    data.frame(
      Node = node,
      State = names(counts),
      Count = as.integer(counts),
      Proportion = round(as.numeric(proportions), 3),
      row.names = NULL
    )
  })
)

#head(data_summary_table,10)

risk_nodes <- c(
  "GridStress",
  "HousingPressure",
  "ServicePressure",
  "TrafficExposure",
  "VisualLandUseExposure",
  "WaterCompetition",
  "ElectricityCostRisk",
  "QualityOfLifeRisk",
  "ReliabilityRisk",
  "WaterSecurityRisk",
  "OverallCommunityRisk"
)

risk_summary <- do.call(
  rbind,
  lapply(risk_nodes, function(node) {
    p <- prop.table(table(synthetic[[node]]))

    data.frame(
      Node = node,
      Low = unname(p["Low"]),
      Moderate = unname(p["Moderate"]),
      High = unname(p["High"])
    )
  })
)

risk_summary

                    Node    Low Moderate   High
1             GridStress 0.2788   0.4614 0.2598
2        HousingPressure 0.2562   0.4352 0.3086
3        ServicePressure 0.2438   0.4474 0.3088
4        TrafficExposure 0.2432   0.4764 0.2804
5  VisualLandUseExposure 0.2632   0.4674 0.2694
6       WaterCompetition 0.2602   0.4496 0.2902
7    ElectricityCostRisk 0.3258   0.4554 0.2188
8      QualityOfLifeRisk 0.2260   0.5550 0.2190
9        ReliabilityRisk 0.3070   0.4132 0.2798
10     WaterSecurityRisk 0.3138   0.4476 0.2386
11  OverallCommunityRisk 0.2068   0.5744 0.2188

round(
  with(
    synthetic,
    prop.table(table(Scale, OverallCommunityRisk), margin = 1)
  ),
  3
)

            OverallCommunityRisk
Scale          Low Moderate  High
  Small      0.350    0.549 0.100
  Medium     0.229    0.598 0.173
  Large      0.165    0.579 0.256
  Hyperscale 0.113    0.535 0.352

# library(dplyr)
# library(tidyr)
# library(ggplot2)

risk_nodes <- c(
  "GridStress",
  "HousingPressure",
  "ServicePressure",
  "TrafficExposure",
  "VisualLandUseExposure",
  "WaterCompetition",
  "ElectricityCostRisk",
  "QualityOfLifeRisk",
  "ReliabilityRisk",
  "WaterSecurityRisk",
  "OverallCommunityRisk"
)

risk_long <- synthetic %>%
  select(all_of(risk_nodes)) %>%
  pivot_longer(
    cols = everything(),
    names_to = "Node",
    values_to = "Risk"
  ) %>%
  count(Node, Risk) %>%
  group_by(Node) %>%
  mutate(Proportion = n / sum(n)) %>%
  ungroup()

risk_long$Risk <- factor(
  risk_long$Risk,
  levels = c("Low", "Moderate", "High")
)

risk_long$Node <- factor(
  risk_long$Node,
  levels = rev(risk_nodes)
)

ggplot(risk_long, aes(x = Proportion, y = Node, fill = Risk)) +
  geom_col(width = 0.75) +
  scale_x_continuous(
    labels = scales::percent_format(),
    limits = c(0, 1)
  ) +
  scale_fill_manual(
    values = c(
      Low = "#2CA25F",
      Moderate = "#FDD835",
      High = "#DE2D26"
    )
  ) +
  labs(
    title = "Distribution of Risk States Across Synthetic Data",
    x = "Proportion of observations",
    y = NULL,
    fill = "Risk state"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    legend.position = "top",
    panel.grid.major.y = element_blank()
  )

high_risk <- synthetic %>%
  select(all_of(risk_nodes)) %>%
  summarise(
    across(
      everything(),
      ~ mean(.x == "High")
    )
  ) %>%
  pivot_longer(
    cols = everything(),
    names_to = "Node",
    values_to = "HighRiskProportion"
  ) %>%
  arrange(HighRiskProportion)

ggplot(
  high_risk,
  aes(
    x = HighRiskProportion,
    y = reorder(Node, HighRiskProportion)
  )
) +
  geom_col(fill = "#DE2D26") +
  scale_x_continuous(
    labels = scales::percent_format(),
    limits = c(0, 1)
  ) +
  labs(
    title = "Prevalence of High-Risk States",
    x = "Observations classified as High",
    y = NULL
  ) +
  theme_minimal(base_size = 12)

# ------------------------------------------------------------
# 5. Fit the Bayesian Network
# ------------------------------------------------------------

fitted_bn <- bn.fit(
  dag,
  data = synthetic,
  method = "bayes",
  iss = 10
)

## Inspect Overall Community Risk, or any other node directly from the object
query <- cpquery(
  fitted_bn,
  event = (OverallCommunityRisk == "High"),
  evidence = (Scale == "Hyperscale")
)

query

[1] 0.3472178

cpquery(
  fitted_bn,
  event = (QualityOfLifeRisk == "High"),
  evidence = (Distance == "Adjacent")
)

[1] 0.4408929

predicted <- predict(
  fitted_bn,
  node = "OverallCommunityRisk",
  data = synthetic,
  prob = TRUE,
  method = "parents"
)

# 1. Extract and transpose the probability matrix
probs <- t(attr(predicted, "prob"))
# 2. Assign column names for clarity
colnames(probs) <- c("Low", "Moderate", "High")
# 3. Generate summary statistics (Min, 1st Qu., Median, Mean, 3rd Qu., Max)
#summary(probs)
rounded_summary <- sapply(as.data.frame(probs), summary)
rounded_summary <- round(rounded_summary, digits = 4)
rounded_summary

           Low Moderate   High
Min.    0.0000   0.0001 0.0001
1st Qu. 0.0003   0.3333 0.0003
Median  0.0005   0.5000 0.0005
Mean    0.2069   0.5742 0.2189
3rd Qu. 0.3333   0.9990 0.3750
Max.    0.9998   0.9998 0.9998

# 1. Extract the probability matrix and convert to a data frame
prob_matrix <- t(attr(predicted, "prob"))
colnames(prob_matrix) <- c("Low", "Moderate", "High")
df_wide <- as.data.frame(prob_matrix)

# 2. Reshape data from WIDE to LONG format for ggplot2
df_long <- df_wide %>%
  pivot_longer(
    cols = everything(), 
    names_to = "Risk_Category", 
    values_to = "Probability"
  ) %>%
  # Ensure the categories stay in a logical order (Low -> Moderate -> High)
  mutate(Risk_Category = factor(Risk_Category, levels = c("Low", "Moderate", "High")))


# 3. Create the ggplot Boxplot
ggplot(df_long, aes(x = Risk_Category, y = Probability, fill = Risk_Category)) +
  # Draw the boxplots with custom outlier styling
  geom_boxplot(
    alpha = 0.7, 
    outlier.color = "grey60", 
    outlier.size = 1.5, 
    outlier.shape = 16
  ) +
  # Overlay a diamond indicating the MEAN probability for each category
  stat_summary(fun = mean, geom = "point", shape = 18, size = 4, color = "black") +
  # Clean, minimal visual theme
  #theme_minimal(base_size = 13) +
  # Custom color palette (Green = Low, Yellow/Orange = Moderate, Red = High)
  scale_fill_manual(values = c("Low" = "#2ecc71", "Moderate" = "#f39c12", "High" = "#e74c3c")) +
  # Descriptive labels
  labs(
    title = "Comparison of Risk Category Probabilities",
    subtitle = "Boxplots showing medians and quartiles; black diamonds represent the mean",
    x = "Risk Category",
    y = "Predicted Probability (0 to 1)"
  ) +
  # Remove the legend since the x-axis already labels the categories
  theme(legend.position = "none")

library(ggtern)

# ----------------------------------------------------------------
# 1. Ingest Your bnlearn Bootstrap Matrix
# ----------------------------------------------------------------
prob_matrix <- t(attr(predicted, "prob"))
colnames(prob_matrix) <- c("Low", "Moderate", "High")
df_wide <- as.data.frame(prob_matrix)

# ----------------------------------------------------------------
# 2. Modernized Ternary Contoured Density Heatmap (No Warnings)
# ----------------------------------------------------------------
plot_density <- ggtern(data = df_wide, aes(x = Low, y = High, z = Moderate)) +
  # Computes and fills a 2D density surface inside the triangular space
  stat_density_tern(
    geom  = "polygon", 
    # FIXED: Replaced deprecated ..level.. with the modern after_stat(level)
    aes(fill = after_stat(level), alpha = after_stat(level)), 
    n     = 200, 
    bdl   = 0.01,
    # FIXED: Mutes the out-of-bounds scale range warning safely
    na.rm = TRUE 
  ) +
  scale_fill_gradient(low = "#ffeda0", high = "#f03b20") +
  scale_alpha_continuous(range = c(0.2, 0.85)) + 
  geom_point(alpha = 0.2, size = 0.8, color = "black") +
  theme_rgbw() +
  labs(
    title    = "Density Contour of Predicted Community Risk Profile",
    subtitle = "Heatmap showing where bootstrap predictions cluster",
    x        = "Low Risk",
    y        = "High Risk",
    z        = "Moderate Risk",
    # FIX 1: Customizing the arrow labels around the triangle (Counterclockwise)
    xarrow   = "Decreasing Risk",
    yarrow   = "Increasing Risk",
    zarrow   = "Increasing Risk"
  ) +
  theme(
    plot.title      = element_text(face = "bold", size = 14),
    legend.position = "none",
    tern.axis.arrow.show = TRUE
    #plot.margin = margin(t = 10, r = 5, b = 10, l = 5, unit = "pt") 
  )

print(plot_density)

# library(ggplot2)
# library(ggridges)
# library(tidyr)

# ----------------------------------------------------------------
# 1. Process Your Exact bnlearn Bootstrap Matrix
# ----------------------------------------------------------------
# Extract the probability matrix and convert to a data frame as specified
prob_matrix <- t(attr(predicted, "prob"))
colnames(prob_matrix) <- c("Low", "Moderate", "High")
df_wide <- as.data.frame(prob_matrix)

# Add a distinct identifier for each bootstrap replication row
df_wide$Bootstrap_ID <- seq_len(nrow(df_wide))

# ----------------------------------------------------------------
# 2. Reshape Data for Ridge Profiling
# ----------------------------------------------------------------
# Pivot data into a long format optimized for ggplot2 tracking
df_long <- pivot_longer(
  df_wide,
  cols      = c("Low", "Moderate", "High"),
  names_to  = "Risk_State",
  values_to = "Predicted_Probability"
)

# Enforce vertical layout hierarchy (High on top, Low on bottom)
df_long$Risk_State <- factor(df_long$Risk_State, levels = c("High", "Moderate", "Low"))

# ----------------------------------------------------------------
# 3. Render Ridgeline Distribution Profiles
# ----------------------------------------------------------------
ggplot(df_long, aes(x = Predicted_Probability, y = Risk_State, fill = Risk_State)) +
  # Adds smooth density waves with the actual bootstrap data points marked underneath as tick marks
  geom_density_ridges(
    alpha = 0.75, scale = 1.4, color = "white", 
    jittered_points = TRUE, point_shape = "|", point_size = 2, point_alpha = 0.3
  ) +
  # Consistent color schema matching the network architecture
  scale_fill_manual(values = c("Low" = "#2b8cbe", "Moderate" = "#feb24c", "High" = "#f03b20")) +
  scale_x_continuous(labels = scales::percent_format(), limits = c(0, 1)) +
  labs(
    title    = "Bootstrapped Target Predictions: Overall Community Risk",
    subtitle = paste("Density distribution profiles across", nrow(df_wide), "bnlearn bootstrap iterations"),
    x        = "Predicted Probability Assigned to State",
    y        = "Risk Classification Spectrum"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    legend.position  = "none",
    plot.title       = element_text(face = "bold"),
    panel.grid.minor = element_blank()
  )

# Extract the conditional-probability table

fitted_bn$HousingPressure$prob

, , WorkforceDemand = Low

               HousingAvailability
HousingPressure   Adequate      Tight SeverelyConstrained
       Low      0.82534532 0.34708510          0.11661093
       Moderate 0.15809619 0.58509638          0.67045708
       High     0.01655850 0.06781852          0.21293200

, , WorkforceDemand = Moderate

               HousingAvailability
HousingPressure   Adequate      Tight SeverelyConstrained
       Low      0.49983482 0.09918699          0.05124471
       Moderate 0.44720846 0.61065041          0.32834061
       High     0.05295672 0.29016260          0.62041467

, , WorkforceDemand = High

               HousingAvailability
HousingPressure   Adequate      Tight SeverelyConstrained
       Low      0.13366947 0.02910443          0.01941528
       Moderate 0.64492997 0.22603268          0.10475223
       High     0.22140056 0.74486289          0.87583249

# Extract the conditional-probability table
prob_table <- fitted_bn$OverallCommunityRisk$prob

# Convert the multidimensional table to a data frame
high_risk_probabilities <- as.data.frame(prob_table)

#dfh <- high_risk_probabilities
# the entire flat table for Overall Community Risk has 19,683 rows
# Rename the probability column
names(high_risk_probabilities)[ncol(high_risk_probabilities)] <- "Probability"

# Keep only cases where OverallCommunityRisk is High
high_risk_probabilities <- subset(
  high_risk_probabilities,
  OverallCommunityRisk == "High"
)

# high_risk_probabilities
# sort the combinations from the highest to the lowest high-risk probability:

high_risk_probabilities <- high_risk_probabilities[
  order(high_risk_probabilities$Probability,decreasing = TRUE),]

# change names to print table
#original_names <- names(high_risk_probabilities)

#Table with short column names and no row names for printing

HRP <- high_risk_probabilities
names(HRP) <- c("ocr", "eco", "ecr", "env", "hpr", "qol", "rel", "spr", "wsr", "Prob")
rownames(HRP) <- NULL
HRP$Prob <- round(HRP$Prob, digits = 4)
head(HRP, 5)

   ocr      eco      ecr  env  hpr      qol      rel  spr  wsr   Prob
1 High Moderate Moderate High High     High     High High High 0.9998
2 High      Low     High High High Moderate Moderate High  Low 0.9997
3 High      Low Moderate High High     High Moderate High  Low 0.9997
4 High      Low Moderate High High     High     High High  Low 0.9997
5 High      Low Moderate High  Low Moderate Moderate  Low High 0.9997

# Change names back
#names(high_risk_probabilities) <- original_names

#some useful redundant code
# risk_flat <- flatten_bnlearn_cpt(
#   fitted_bn = fitted_bn,
#   child_name = "OverallCommunityRisk"
# )
# 
# head(risk_flat)
# 
# 
# # Descending
# sorted_risk_flat <- risk_flat %>%
#   arrange(desc("Probability"))
# 
# head(sorted_risk_flat)

HRP[333,]

     ocr eco  ecr      env      hpr      qol      rel      spr      wsr   Prob
333 High Low High Moderate Moderate Moderate Moderate Moderate Moderate 0.7999

### First or random
plot_radar <- function(data, n_highlight = 5, selection_method = c("first", "random"), level_order = c("Low", "Moderate", "High"), tleg = 45) {
  selection_method <- match.arg(selection_method)
  
  categorical_vars <- c(
    "OverallCommunityRisk", "EconomicBenefit", "ElectricityCostRisk",
    "EnvironmentalDegradation", "HousingPressure", "QualityOfLifeRisk",
    "ReliabilityRisk", "ServicePressure", "WaterSecurityRisk"
  )
  
  required_vars <- c(categorical_vars, "Probability")
  missing_vars <- setdiff(required_vars, names(data))
  if (length(missing_vars) > 0) {
    stop(
      "Missing columns: ", paste(missing_vars, collapse = ", ")
    )
  }
  
  data$Probability <- as.numeric(data$Probability)
  
  # Capture original dataframe row numbers before filtering and sorting
  plot_data <- data %>% 
    dplyr::mutate(orig_row = dplyr::row_number()) %>%
    dplyr::select(dplyr::all_of(required_vars), orig_row) %>% 
    dplyr::filter(dplyr::if_all(dplyr::all_of(required_vars), ~ !is.na(.x)))
  
  if (nrow(plot_data) == 0) {
    stop("There are no complete rows available for plotting.")
  }
  
  # Sort to isolate the highest-probability rows
  top_pool <- plot_data %>% 
    dplyr::arrange(dplyr::desc(Probability))
  
  n_top <- min(n_highlight, nrow(top_pool))
  
  # Select paths based on user choice
  if (selection_method == "first") {
    plot_data <- top_pool %>% dplyr::slice_head(n = n_top)
  } else {
    plot_data <- top_pool %>% dplyr::slice_sample(n = n_top)
  }
  
  plot_data <- plot_data %>% 
    dplyr::mutate(
      path = dplyr::row_number(),
      path_label = paste0("Top path ", path)
    )
  
  # Convert Low, Moderate, High to 1, 2, 3
  for (v in categorical_vars) {
    unexpected <- setdiff(
      unique(plot_data[[v]]), level_order
    )
    if (length(unexpected) > 0) {
      stop(
        "Unexpected level(s) in ", v, ": ", paste(unexpected, collapse = ", ")
      )
    }
    plot_data[[v]] <- as.numeric(
      factor(
        plot_data[[v]], levels = level_order, ordered = TRUE
      )
    )
  }
  
  n_axes <- length(categorical_vars)
  angles <- seq(
    from = pi / 2, to = pi / 2 - 2 * pi, length.out = n_axes + 1
  )[-(n_axes + 1)]
  
  axis_positions <- data.frame(
    variable = categorical_vars,
    angle = angles,
    stringsAsFactors = FALSE
  )
  
  # Convert radar values to Cartesian coordinates
  radar_data <- plot_data %>% 
    tidyr::pivot_longer(
      cols = dplyr::all_of(categorical_vars),
      names_to = "variable",
      values_to = "value"
    ) %>% 
    dplyr::left_join(axis_positions, by = "variable") %>% 
    dplyr::mutate(
      x = value * cos(angle),
      y = value * sin(angle),
      path_label = factor(path_label, levels = paste0("Top path ", seq_len(n_top)))
    )
  
  # Close each profile polygon by repeating its first point
  radar_data_closed <- radar_data %>% 
    dplyr::group_by(path, path_label) %>% 
    dplyr::arrange(angle, .by_group = TRUE) %>% 
    dplyr::slice(c(dplyr::row_number(), 1)) %>% 
    dplyr::ungroup()
  
  # Polygonal grid and outer boundary
  grid_data <- tidyr::crossing(
    radius = 1:3,
    axis_positions
  ) %>% 
    dplyr::mutate(
      x = radius * cos(angle),
      y = radius * sin(angle)
    ) %>% 
    dplyr::group_by(radius) %>% 
    dplyr::arrange(angle, .by_group = TRUE) %>% 
    dplyr::slice(c(dplyr::row_number(), 1)) %>% 
    dplyr::ungroup()
  
  # Axis spokes
  spoke_data <- axis_positions %>% 
    dplyr::mutate(
      x = 3 * cos(angle),
      y = 3 * sin(angle)
    )
  
  highlight_colours <- c(
    "#D73027", "#4575B4", "#1A9850", "#984EA3", "#E66101"
  )[seq_len(n_top)]
  
  colour_values <- stats::setNames(highlight_colours, paste0("Top path ", seq_len(n_top)))
  
  # Legend labels matching selection strategy with 5 significant digits & row numbers
  legend_labels <- plot_data %>% 
    dplyr::arrange(path) %>% 
    dplyr::mutate(
      label = paste0(
        "[Row ", orig_row, "] (risk = ", format(signif(Probability, 5), scientific = FALSE), ")"
      )
    ) %>% 
    dplyr::pull(label)
  
  names(legend_labels) <- paste0("Top path ", seq_len(n_top))
  
  # Axis-label positions
  label_data <- axis_positions %>% 
    dplyr::mutate(
      x = 3.25 * cos(angle),
      y = 3.25 * sin(angle),
      hjust = dplyr::case_when(
        cos(angle) > 0.1 ~ 0,
        cos(angle) < -0.1 ~ 1,
        TRUE ~ 0.5
      ),
      vjust = dplyr::case_when(
        sin(angle) > 0.1 ~ 0,
        sin(angle) < -0.1 ~ 1,
        TRUE ~ 0.5
      )
    )
  
  ggplot2::ggplot() +
    ggplot2::geom_path(
      data = grid_data,
      ggplot2::aes(x = x, y = y, group = radius),
      colour = "grey75", linewidth = 0.5
    ) +
    ggplot2::geom_segment(
      data = spoke_data,
      ggplot2::aes(x = 0, y = 0, xend = x, yend = y),
      colour = "grey80", linewidth = 0.5
    ) +
    ggplot2::geom_polygon(
      data = radar_data_closed,
      ggplot2::aes(x = x, y = y, group = path, colour = path_label, fill = path_label),
      alpha = 0.08, linewidth = 0.5
    ) +
    ggplot2::geom_path(
      data = radar_data_closed,
      ggplot2::aes(x = x, y = y, group = path, colour = path_label),
      linewidth = 0.9, alpha = 0.85
    ) +
    ggplot2::geom_point(
      data = radar_data,
      ggplot2::aes(x = x, y = y, colour = path_label),
      size = 1.5, alpha = 0.9
    ) +
    ggplot2::geom_text(
      data = label_data,
      ggplot2::aes(x = x, y = y, label = variable, hjust = hjust, vjust = vjust),
      size = 3.5, fontface = "bold"
    ) +
    ggplot2::annotate(
      "text", x = 0, y = 1, label = level_order[1], colour = "grey45", size = 3, vjust = 1.5
    ) +
    ggplot2::annotate(
      "text", x = 0, y = 2, label = level_order[2], colour = "grey45", size = 3, vjust = 1.5
    ) +
    ggplot2::annotate(
      "text", x = 0, y = 3, label = level_order[3], colour = "grey45", size = 3, vjust = 1.5
    ) +
    ggplot2::scale_colour_manual(
      values = colour_values,
      breaks = names(legend_labels),
      labels = unname(legend_labels),
      name = NULL
    ) +
    ggplot2::scale_fill_manual(
      values = colour_values, guide = "none"
    ) +
    ggplot2::coord_equal(
      xlim = c(-4.5, 4.5), ylim = c(-3.5, 3.8), clip = "off"
    ) +
    ggplot2::labs(
      title = paste(
        "Overall Community Risk radar profiles:", n_top, "paths selected via", selection_method
      ),
      x = NULL, y = NULL
    ) +
    ggplot2::theme(
      legend.margin = ggplot2::margin(t = tleg, r = 0, b = 0, l = 0)
    )
}

#df <- sorted_risk_flat
set.seed(456)
df <- high_risk_probabilities

plot_radar(
  data=df, n_highlight = 3, selection_method = "random", level_order = c("Low", "Moderate", "High"), tleg = 195
)

summary(high_risk_probabilities$Probability)

     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
5.644e-05 5.073e-04 3.333e-01 2.859e-01 3.333e-01 9.998e-01 

# ================================================================
# Usage
# ================================================================
#
# A DAG named `dag` must exist first.

# Global scenario controls.
  # global_risk_shift
  # 0.00 = baseline moderate-risk scenario
  # Positive values increase risk throughout the network
  # Negative values decrease risk throughout the network
  
  #  parent_effect_strength
  # Controls how strongly parent risk affects child risk.
  # 1.00 is the default.
 
  # conditional_noise 
  # Controls random variation in conditional probabilities.
  # Larger values produce less deterministic relationships.

#
synthetic_2 <- generate_synthetic_data(
  n = 5000,
  global_risk_shift = 0.20,
  parent_effect_strength = 1.00,
  conditional_noise = 0.10
)
#
# str(synthetic)
# head(synthetic)
# ================================================================

risk_nodes <- c(
  "GridStress",
  "HousingPressure",
  "ServicePressure",
  "TrafficExposure",
  "VisualLandUseExposure",
  "WaterCompetition",
  "ElectricityCostRisk",
  "QualityOfLifeRisk",
  "ReliabilityRisk",
  "WaterSecurityRisk",
  "OverallCommunityRisk"
)

risk_summary <- do.call(
  rbind,
  lapply(risk_nodes, function(node) {
    p <- prop.table(table(synthetic_2[[node]]))

    data.frame(
      Node = node,
      Low = unname(p["Low"]),
      Moderate = unname(p["Moderate"]),
      High = unname(p["High"])
    )
  })
)

risk_summary

                    Node    Low Moderate   High
1             GridStress 0.0698   0.3416 0.5886
2        HousingPressure 0.0606   0.2964 0.6430
3        ServicePressure 0.0596   0.3144 0.6260
4        TrafficExposure 0.0608   0.3348 0.6044
5  VisualLandUseExposure 0.0610   0.3366 0.6024
6       WaterCompetition 0.0598   0.3144 0.6258
7    ElectricityCostRisk 0.0732   0.3844 0.5424
8      QualityOfLifeRisk 0.0262   0.2422 0.7316
9        ReliabilityRisk 0.0556   0.2348 0.7096
10     WaterSecurityRisk 0.0748   0.3656 0.5596
11  OverallCommunityRisk 0.0194   0.1638 0.8168

# ------------------------------------------------------------
# 5. Fit the Bayesian Network
# ------------------------------------------------------------

fitted_bn_2 <- bn.fit(
  dag,
  data = synthetic_2,
  method = "bayes",
  iss = 10
)

## Inspect Overall Community Risk, or any other node directly from the object

predicted_2 <- predict(
  fitted_bn_2,
  node = "OverallCommunityRisk",
  data = synthetic_2,
  prob = TRUE
)

# 1. Extract and transpose the probability matrix
probs <- t(attr(predicted_2, "prob"))
# 2. Assign column names for clarity
colnames(probs) <- c("Low", "Moderate", "High")
# 3. Generate summary statistics (Min, 1st Qu., Median, Mean, 3rd Qu., Max)
#summary(probs)
rounded_summary <- sapply(as.data.frame(probs), summary)
rounded_summary <- round(rounded_summary, digits = 4)
rounded_summary

           Low Moderate   High
Min.    0.0000   0.0000 0.0003
1st Qu. 0.0000   0.0005 0.7777
Median  0.0002   0.1183 0.8767
Mean    0.0195   0.1638 0.8167
3rd Qu. 0.0068   0.2076 0.9434
Max.    0.9990   0.9995 0.9999

#library(ggtern)

# ----------------------------------------------------------------
# 1. Ingest Your bnlearn Bootstrap Matrix
# ----------------------------------------------------------------
prob_matrix <- t(attr(predicted_2, "prob"))
colnames(prob_matrix) <- c("Low", "Moderate", "High")
df_wide <- as.data.frame(prob_matrix)

# ----------------------------------------------------------------
# 2. Modernized Ternary Contoured Density Heatmap (No Warnings)
# ----------------------------------------------------------------
plot_density <- ggtern(data = df_wide, aes(x = Low, y = High, z = Moderate)) +
  # Computes and fills a 2D density surface inside the triangular space
  stat_density_tern(
    geom  = "polygon",
    # FIXED: Replaced deprecated ..level.. with the modern after_stat(level)
    aes(fill = after_stat(level), alpha = after_stat(level)),
    n     = 200,
    bdl   = 0.01,
    # FIXED: Mutes the out-of-bounds scale range warning safely
    na.rm = TRUE
  ) +
  scale_fill_gradient(low = "#ffeda0", high = "#f03b20") +
  scale_alpha_continuous(range = c(0.4, 0.99)) +
  geom_point(alpha = 0.2, size = 0.8, color = "black") +
  theme_rgbw() +
  labs(
    title    = "Density Contour of Predicted Community Risk Profile",
    subtitle = "Heatmap showing where bootstrap predictions cluster",
    x        = "Low Risk",
    y        = "High Risk",
    z        = "Mod Risk",
    xarrow   = "Decreasing Risk",
    yarrow   = "Increasing Risk",
    zarrow   = "Increasing Risk"
  ) +
  theme(
    plot.title      = element_text(face = "bold", size = 14),
    legend.position = "none",
    tern.axis.arrow.show = TRUE
  )

plot_density

# plot_density <- ggtern(
#   data = df_wide,
#   aes(x = Low, y = High, z = Moderate)
# ) +
#   stat_density_tern(
#     geom = "polygon",
#     aes(
#       fill = after_stat(level),
#       group = after_stat(level)
#     ),
#     n = 200,
#     bdl = 0.01,
#     na.rm = TRUE
#   ) +
#   scale_fill_gradient(
#     low = "#ffeda0",
#     high = "#f03b20"
#   ) +
#   geom_point(
#     alpha = 0.2,
#     size = 0.8,
#     color = "black"
#   ) +
#   theme_rgbw() +
#   theme(
#     plot.title = element_text(face = "bold", size = 14),
#     legend.position = "none",
#     tern.axis.arrow.show = TRUE
#   )
# 
# plot_density

count_displayed_cpt_tables <- function(fitted_bn) {
  
  result <- lapply(
    names(fitted_bn),
    function(node_name) {
      
      prob <- fitted_bn[[node_name]]$prob
      
      if (is.null(dim(prob))) {
        
        # Root node: a probability vector, not a 2D table
        data.frame(
          Node = node_name,
          CPTDimensions = length(prob),
          NumberOfTables = 1L,
          stringsAsFactors = FALSE
        )
        
      } else {
        
        dimensions <- dim(prob)
        
        # For a 3D CPT, the printed tables are slices of the third dimension.
        # For higher-dimensional CPTs, count all slices after dimensions 1 and 2.
        number_of_tables <- if (length(dimensions) <= 2L) {
          1L
        } else {
          prod(dimensions[-c(1, 2)])
        }
        
        data.frame(
          Node = node_name,
          CPTDimensions = paste(dimensions, collapse = " x "),
          NumberOfTables = number_of_tables,
          stringsAsFactors = FALSE
        )
      }
    }
  )
  
  result <- do.call(rbind, result)
  rownames(result) <- NULL
  
  result
}

cpt_table_counts <- count_displayed_cpt_tables(fitted_bn)

cpt_table_counts

                       Node                     CPTDimensions NumberOfTables
1          AccessToBenefits                         3 x 6 x 4              4
2        AirQualityExposure                     3 x 3 x 6 x 4             24
3          BackupGeneration                             3 x 4              1
4     ConstructionIntensity                             3 x 4              1
5                   Cooling                                 4              1
6                  Distance                                 6              1
7           EconomicBenefit                         3 x 3 x 4              4
8       ElectricityCostRisk                         3 x 6 x 3              3
9         ElectricityDemand                             4 x 4              1
10 EnvironmentalDegradation                     3 x 3 x 3 x 3              9
11             GridCapacity                                 3              1
12               GridStress                         3 x 4 x 3              3
13      HousingAvailability                                 3              1
14          HousingPressure                         3 x 3 x 3              3
15            LandFootprint                             3 x 4              1
16            NoiseExposure                     3 x 4 x 6 x 4             24
17     OverallCommunityRisk 3 x 3 x 3 x 3 x 3 x 3 x 3 x 3 x 3           2187
18        QualityOfLifeRisk         3 x 3 x 3 x 3 x 3 x 3 x 3            243
19          ReliabilityRisk                             3 x 3              1
20                    Scale                                 4              1
21          ServicePressure                         3 x 3 x 3              3
22          TrafficExposure                         3 x 3 x 6              6
23    VisualLandUseExposure                         3 x 6 x 3              3
24            Vulnerability                                 3              1
25        WaterAvailability                                 3              1
26         WaterCompetition                         3 x 3 x 3              3
27              WaterDemand                         3 x 4 x 4              4
28        WaterSecurityRisk                         3 x 3 x 3              3
29          WorkforceDemand                             3 x 4              1

tot_tables <- sum(cpt_table_counts$NumberOfTables)

# ----------------------------------------------------------------
# 1. Define Core Simulation Parameters
# ----------------------------------------------------------------
concentration <- 6
state_scores  <- c(Low = 0.15, Moderate = 0.50, High = 0.85)
target_range  <- seq(0.00, 1.00, by = 0.01)

# ----------------------------------------------------------------
# 2. Compute Normalized Probabilities over Target Risk Continuum
# ----------------------------------------------------------------
plot_data <- data.frame()

for (target_risk in target_range) {
  # Step 3: Measure raw exponential distances
  raw_probs <- exp(-concentration * abs(state_scores - target_risk))
  
  # Step 4: Normalize scores to sum to 1.00
  final_probs <- raw_probs / sum(raw_probs)
  
  # Bind results for plotting
  row_df <- data.frame(
    TargetRisk = target_risk,
    Low        = final_probs["Low"],
    Moderate   = final_probs["Moderate"],
    High       = final_probs["High"]
  )
  plot_data <- rbind(plot_data, row_df)
}

# Pivot data into a long format optimized for ggplot2 tracking
plot_data_long <- pivot_longer(
  plot_data, 
  cols      = c("Low", "Moderate", "High"), 
  names_to  = "State", 
  values_to = "Probability"
)

# Enforce clean factor ordering matching the network architecture
plot_data_long$State <- factor(plot_data_long$State, levels = c("Low", "Moderate", "High"))

# ----------------------------------------------------------------
# 3. Generate the RMarkdown Visualization
# ----------------------------------------------------------------
ggplot(plot_data_long, aes(x = TargetRisk, y = Probability, color = State)) +
  geom_line(size = 1.2) +
  # Visual reference lines matching the static state risk mappings
  geom_vline(xintercept = c(0.15, 0.50, 0.85), linetype = "dashed", color = "gray60", alpha = 0.7) +
  scale_color_manual(values = c("Low" = "#2b8cbe", "Moderate" = "#feb24c", "High" = "#f03b20")) +
  labs(
    title    = "Dynamic Probability Shifts for Derived Non-Root Nodes",
    subtitle = paste0("Exponential RBF Kernel Behavior (Concentration = ", concentration, ")"),
    x        = "Calculated Target Risk Score (From Parents + Noise)",
    y        = "Final State Generation Probability",
    color    = "Assigned State"
  ) +
  scale_x_continuous(breaks = seq(0, 1, 0.1)) +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1)) 

  #theme_minimal(base_size = 12) +
  # theme(
  #   plot.title         = element_text(face = "bold", linewidth = 1),
  #   panel.grid.minor   = element_blank(),
  #   legend.position    = "bottom"
  #)

In my last post on permute, I looked at ordered set partitions and the new formula interface for plots and blocks. This time I want to show two helpers for the next part of writing a permutation test: taking those permutations and doing something with them.

permutations() gives you the shuffled values directly. permutator() gives you a function that hands out the next set of observation indices each time you call it. Both use the same permutation machinery underneath, so we can use them to build the same test in two slightly different ways.

These functions are in the development version of permute. To follow along, install it from my R-universe:

install.packages("permute", repos = c(
  "https://gavinsimpson.r-universe.dev",
  "https://cloud.r-project.org"
))

library("permute")
data(PlantGrowth, package = "datasets")

aggregate(weight ~ group, data = PlantGrowth, FUN = mean)

  group weight
1  ctrl  5.032
2  trt1  4.661
3  trt2  5.526

f_stat <- function(weight, group) {
  fit <- lm(weight ~ group)
  unname(summary(fit)$fstatistic[["value"]])
}

observed <- with(PlantGrowth, f_stat(weight, group))
observed

[1] 4.846088

ctrl <- how(nperm = 9999)

set.seed(42)
shuffled <- permutations(PlantGrowth$weight, control = ctrl)

dim(shuffled)

[1] 9999   30

as.matrix(shuffled)[seq_len(3), seq_len(5)]

     [,1] [,2] [,3] [,4] [,5]
[1,] 6.03 4.50 4.17 5.37 5.14
[2,] 4.81 5.87 5.12 5.29 4.53
[3,] 5.14 4.53 4.50 4.17 6.03

null_values <- apply(shuffled, 1, f_stat, group = PlantGrowth$group)

tol <- sqrt(.Machine$double.eps)
p_value <- (1 + sum(null_values >= observed - tol)) /
  (length(null_values) + 1)

c(F = observed, p_value = p_value)

       F  p_value 
4.846088 0.014200 

old_par <- par(mfrow = c(1, 2))
boxplot(weight ~ group, data = PlantGrowth,
        xlab = "Group", ylab = "Dried plant weight",
        col = "grey85")
hist(null_values, breaks = 40, main = "",
     xlab = "F statistic after shuffling weights",
     col = "grey85", border = "white")
abline(v = observed, col = "firebrick", lwd = 2)
par(old_par)

set.seed(42)
next_perm <- permutator(nrow(PlantGrowth), control = ctrl)

null_iter <- numeric(attr(next_perm, "nperm"))

for (b in seq_along(null_iter)) {
  i <- next_perm()
  null_iter[b] <- f_stat(PlantGrowth$weight[i], PlantGrowth$group)
}

next_perm()

NULL

all.equal(null_values, null_iter)

[1] TRUE

p_iter <- (1 + sum(null_iter >= observed - tol)) /
  (length(null_iter) + 1)

c(values = p_value, iterator = p_iter)

  values iterator 
  0.0142   0.0142 

set.seed(42)
indices <- shuffleSet(nrow(PlantGrowth), control = ctrl)

next_perm <- permutator(perms = indices)
i <- next_perm()
head(PlantGrowth$weight[i])

[1] 6.03 4.50 4.17 5.37 5.14 6.11

replay <- permutator(perms = indices)
identical(i, replay())

[1] TRUE

What difference can a sponsored place on a bioinformatics course make? For some researchers, it means access to training they would otherwise struggle to afford. For others, it is an opportunity to explore a new field, work through unfamiliar analytical methods, or find new ways to approach an existing research question.

Through the Bioconductor–Physalia collaboration, Physalia offers two free places on eligible courses to researchers based in low- and middle-income countries (LMICs). The Bioconductor Community Manager coordinates the application process and works with the Training Committee co-chairs on selecting participants, while Bioconductor helps promote the courses through its community channels. Between November 2025 and October 2026, the programme provided 18 sponsored places across nine courses, covering topics ranging from RNA-seq and proteomics to microbiome analysis, epigenomics and network biology.

The aim is to make specialist training more accessible to researchers who might otherwise miss out on these opportunities. We asked participants what they found useful about their courses and how they planned to use what they had learned. Their responses offer a glimpse into the different ways the training connects with their work.

© 2026 Bioconductor. Content is published under Creative Commons CC-BY-4.0 License for the text and BSD 3-Clause License for any code. | R-Bloggers

I’ve been giving permute a bit of attention recently. Two additions should make it easier to set up permutation tests: an option for ordered set partitions, and formulas for specifying plots and blocks. That’s a rather grand name for the first addition, so let’s start with what it does.

These changes are in the development version. To follow along, install it from my R-universe:

install.packages("permute", repos = c(
  "https://gavinsimpson.r-universe.dev",
  "https://cloud.r-project.org"
))

library("permute")
data(jackal, package = "permute")

table(jackal$Sex)

  Male Female 
    10     10 

numPerms(nrow(jackal), control = how(within = Within(type = "free")))

[1] 2.432902e+18

ctrl <- how(
  plots = Plots(strata = ~ Sex, type = "partition"),
  data = jackal
)

numPerms(nrow(jackal), control = ctrl)

[1] 184756

set.seed(42)
p <- shuffle(nrow(jackal), control = ctrl)

table(jackal$Sex[p])

  Male Female 
    10     10 

fit <- lm(Length ~ Sex, data = jackal)
coef(fit)

(Intercept)   SexFemale 
      113.4        -4.8 

set.seed(42)
perms <- shuffleSet(nrow(jackal), nset = 999, control = ctrl)

observed <- unname(coef(fit)[2])
null <- apply(perms, 1, function(i) {
  permuted <- jackal
  permuted$Sex <- jackal$Sex[i]
  unname(coef(lm(Length ~ Sex, data = permuted))[2])
})

hist(null, breaks = "FD", main = "",
     xlab = "Female minus male mean mandible length (mm)",
     col = "grey85", border = "white")
abline(v = observed, col = "firebrick", lwd = 2)

tol <- sqrt(.Machine$double.eps)
p_value <- (1 + sum(abs(null) >= abs(observed) - tol)) /
  (nrow(perms) + 1)

c(difference_mm = observed, p_value = p_value)

difference_mm       p_value 
       -4.800         0.005 

library("vegan")
data(dune, package = "vegan")
data(dune.env, package = "vegan")

formula_ctrl <- how(
  plots = Plots(strata = ~ Management / Use),
  blocks = ~ Management,
  data = dune.env
)

set.seed(42)
i <- shuffle(nrow(dune.env), control = formula_ctrl)

all(getBlocks(formula_ctrl)[i] == getBlocks(formula_ctrl))

[1] TRUE

all(getStrata(formula_ctrl)[i] == getStrata(formula_ctrl))

[1] TRUE

env <- dune.env
env$group <- getStrata(formula_ctrl)
ord <- rda(dune ~ A1 + Condition(group), data = env)

set.seed(42)
anova(ord, permutations = formula_ctrl)

Permutation test for rda under reduced model
Blocks:  Management 
Plots: Management/Use, plot permutation: none
Permutation: free
Number of permutations: 199

Model: rda(formula = dune ~ A1 + Condition(group), data = env)
         Df Variance      F Pr(>F)
Model     1   4.3492 1.3265  0.225
Residual  8  26.2298              

plot_groups <- Plots(strata = ~ Management / Use, data = dune.env)

In statistics, dependence or association is any statistical relationship, whether causal or not, between two random variables or bivariate data. It is often measured by the Pearson’s correlation coefficient:

Other measures of correlation can be thought of, such as the Spearman rank correlation coefficient or Kendall rank correlation coefficient.

r <- cor(mtcars$cyl, mtcars$gear)
r
## [1] -0.4926866

T <- r * sqrt(32 - 2) / sqrt(1 - r^2)
T
## [1] -3.101051

2 * pt(T, 30)
## [1] 0.004173297

# Correlations with significance levels
library(Hmisc)
corr2 <- rcorr(as.matrix(mtcars), type="pearson")
print(corr2$r, digits = 2)
##        mpg   cyl  disp    hp   drat    wt   qsec    vs     am  gear   carb
## mpg   1.00 -0.85 -0.85 -0.78  0.681 -0.87  0.419  0.66  0.600  0.48 -0.551
## cyl  -0.85  1.00  0.90  0.83 -0.700  0.78 -0.591 -0.81 -0.523 -0.49  0.527
## disp -0.85  0.90  1.00  0.79 -0.710  0.89 -0.434 -0.71 -0.591 -0.56  0.395
## hp   -0.78  0.83  0.79  1.00 -0.449  0.66 -0.708 -0.72 -0.243 -0.13  0.750
## drat  0.68 -0.70 -0.71 -0.45  1.000 -0.71  0.091  0.44  0.713  0.70 -0.091
## wt   -0.87  0.78  0.89  0.66 -0.712  1.00 -0.175 -0.55 -0.692 -0.58  0.428
## qsec  0.42 -0.59 -0.43 -0.71  0.091 -0.17  1.000  0.74 -0.230 -0.21 -0.656
## vs    0.66 -0.81 -0.71 -0.72  0.440 -0.55  0.745  1.00  0.168  0.21 -0.570
## am    0.60 -0.52 -0.59 -0.24  0.713 -0.69 -0.230  0.17  1.000  0.79  0.058
## gear  0.48 -0.49 -0.56 -0.13  0.700 -0.58 -0.213  0.21  0.794  1.00  0.274
## carb -0.55  0.53  0.39  0.75 -0.091  0.43 -0.656 -0.57  0.058  0.27  1.000
print(corr2$P, digits = 2)
##          mpg     cyl    disp      hp    drat      wt    qsec      vs      am
## mpg       NA 6.1e-10 9.4e-10 1.8e-07 1.8e-05 1.3e-10 1.7e-02 3.4e-05 2.9e-04
## cyl  6.1e-10      NA 1.8e-12 3.5e-09 8.2e-06 1.2e-07 3.7e-04 1.8e-08 2.2e-03
## disp 9.4e-10 1.8e-12      NA 7.1e-08 5.3e-06 1.2e-11 1.3e-02 5.2e-06 3.7e-04
## hp   1.8e-07 3.5e-09 7.1e-08      NA 1.0e-02 4.1e-05 5.8e-06 2.9e-06 1.8e-01
## drat 1.8e-05 8.2e-06 5.3e-06 1.0e-02      NA 4.8e-06 6.2e-01 1.2e-02 4.7e-06
## wt   1.3e-10 1.2e-07 1.2e-11 4.1e-05 4.8e-06      NA 3.4e-01 9.8e-04 1.1e-05
## qsec 1.7e-02 3.7e-04 1.3e-02 5.8e-06 6.2e-01 3.4e-01      NA 1.0e-06 2.1e-01
## vs   3.4e-05 1.8e-08 5.2e-06 2.9e-06 1.2e-02 9.8e-04 1.0e-06      NA 3.6e-01
## am   2.9e-04 2.2e-03 3.7e-04 1.8e-01 4.7e-06 1.1e-05 2.1e-01 3.6e-01      NA
## gear 5.4e-03 4.2e-03 9.6e-04 4.9e-01 8.4e-06 4.6e-04 2.4e-01 2.6e-01 5.8e-08
## carb 1.1e-03 1.9e-03 2.5e-02 7.8e-07 6.2e-01 1.5e-02 4.5e-05 6.7e-04 7.5e-01
##         gear    carb
## mpg  5.4e-03 1.1e-03
## cyl  4.2e-03 1.9e-03
## disp 9.6e-04 2.5e-02
## hp   4.9e-01 7.8e-07
## drat 8.4e-06 6.2e-01
## wt   4.6e-04 1.5e-02
## qsec 2.4e-01 4.5e-05
## vs   2.6e-01 6.7e-04
## am   5.8e-08 7.5e-01
## gear      NA 1.3e-01
## carb 1.3e-01      NA

# Correlation matrix from mtcars
x <- mtcars[1:3]
y <- mtcars[4:6]
cor(x, y)
##              hp       drat         wt
## mpg  -0.7761684  0.6811719 -0.8676594
## cyl   0.8324475 -0.6999381  0.7824958
## disp  0.7909486 -0.7102139  0.8879799

library(MASS)
cov <- matrix(c(2.20, 0.48, 0.48, 0.20), 2, 2)
cov
##      [,1] [,2]
## [1,] 2.20 0.48
## [2,] 0.48 0.20

set.seed(1234)
dataset <- data.frame( mvrnorm(n=10, mu = c(10, 2), Sigma = cov) )
names(dataset) <- c("X", "Y")
dataset
##            X        Y
## 1  11.756647 2.547203
## 2   9.522180 2.199740
## 3   8.341254 1.862362
## 4  13.480005 2.772031
## 5   9.428296 1.573435
## 6   9.242788 1.861756
## 7  10.817449 2.343918
## 8  10.749047 2.451999
## 9  10.780400 2.436263
## 10 11.480301 1.590436

r <- with(dataset, cor(X, Y))
b1 <- coef( lm(Y ~ X, data=dataset) )[2]
r
## [1] 0.6372927
b1
##         X 
## 0.1785312

sigmaX <- with(dataset, sd(X) )
sigmaY <- with(dataset, sd(Y) )
b1 * sigmaX / sigmaY 
##         X 
## 0.6372927
r * sigmaY / sigmaX
## [1] 0.1785312

Yst <- with(dataset, scale(Y, scale=T) )
summary( lm(Yst ~ I(scale(X, scale = T) ), data = dataset) )
## 
## Call:
## lm(formula = Yst ~ I(scale(X, scale = T)), data = dataset)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.082006 -0.067143 -0.036850  0.009214  0.237923 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)  
## (Intercept)            1.065e-18  3.478e-02   0.000   1.0000  
## I(scale(X, scale = T)) 1.785e-01  7.633e-02   2.339   0.0475 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.11 on 8 degrees of freedom
## Multiple R-squared:  0.4061, Adjusted R-squared:  0.3319 
## F-statistic: 5.471 on 1 and 8 DF,  p-value: 0.04748

Intra-class correlation (ICC)

This is a related yet distinct concept. While Pearson’s correlation coefficient is used for bivariate data—where two different variables are observed on the same independent units—intraclass correlation arises when the same variable is repeatedly observed on the same units (e.g., across different times or spaces), so that each unit acts as a grouping factor that creates clusters of measurements.

Let’s illustrate this with an example. Consider the following data table:

#|echo: false
m1 <- c(10, 12, 13, 16)
m2 <- c(21, 20, 17, 19)
dtab <- cbind(m1, m2)
row.names(dtab) <- paste0("g", 1:4)
knitr::kable(dtab)

	m1	m2
g1	10	21
g2	12	20
g3	13	17
g4	16	19

When g1, g2, g3, and g4 represent four independent units (e.g., four maize plants) on which two variables, m1 and m2 (e.g., weight and height), were measured, we are in the realm of Pearson’s correlation. On the other hand, when g1, g2, g3, and g4 are entities on which multiple measurements of the same variable (e.g., yield) are taken, we enter the realm of intraclass correlation. In this second scenario, we are not assessing joint variability between two different variables; rather, we want to determine whether belonging to the same group (‘proximity’) makes observations more similar to each other than observations from different groups.

The most common intraclass correlation coefficient is:

where is the between-group variance component and is the within-group variance component (or, more clearly, the variance of observations within a group). For balanced data, those two components can be obtained from one-way ANOVA, considering that there are four ‘replicates’ per group:

y <- c(10, 12, 13, 16, 21, 20, 17, 19)
group <- factor( rep(c("A", "B"), each = 4) )
nmeasures <- 4
mod <- lm(y ~ group)
anova(mod)
## Analysis of Variance Table
## 
## Response: y
##           Df Sum Sq Mean Sq F value   Pr(>F)   
## group      1   84.5  84.500  18.436 0.005128 **
## Residuals  6   27.5   4.583                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
sigma2_alpha <- (84.5 - 4.583)/nmeasures
sigma2_epsilon <- 4.583
icc <- sigma2_alpha/(sigma2_alpha + sigma2_epsilon)
icc
## [1] 0.8134129

In general, if subjects in one group are very similar (small ), while different groups vary widely (high ), the ICC is high. Conversely, if measurements are just as similar across different groups as they are within the same group, the ICC is low. In theory, the ICC cannot be negative, as variance components are expected to be positive. In practice, however, the ICC calculated via the method of moments shown above can turn out negative when the group mean square is lower than the residual mean square from the ANOVA. Such negative values are usually constrained to zero or interpreted as indicating an absence of positive clustering; in both cases, intraclass correlation remains a meaningful tool.

ICC is crucial and is most often studied using linear mixed models. But… that is a story for another day…

Thanks for reading, and happy coding!

(P.S. If you enjoyed this post, check out my book “Field Research Methods in Agriculture: An Introduction with R” via the link below!)

Prof. Andrea Onofri
Department of Agricultural, Food and Environmental Sciences
University of Perugia (Italy)
Send comments to: andrea.onofri@unipg.it

---

This post was originally published on 2019-02-07