With BiocAsia 2026 coming up in Melbourne, community members have been organising seminars and exploring ways to take part locally.

More than 200 people joined the latest BiocAsia seminar on 13 October 2026. The session ran for two and a half hours, with discussion continuing well beyond the planned finish. It was the first seminar in the series hosted from China, with Mandarin as the main language. Shixiang Wang spoke about building an open-source ecosystem around UCSC Xena and visualisation tools, and Zhougeng Xu shared his experience developing open-source tools for AI in science.

The seminar was organised by Zuguang Gu, a member of the Bioconductor Community Advisory Board and an active contributor to BiocAsia based in China. The next seminar will take place in December.

Alongside the seminars, Jyoti Bhogal and Malvika Kharbanda have been discussing watch party options in India for BiocAsia 2026. The aim is to give people who cannot travel to Melbourne a chance to take part in the conference and meet others nearby.

Join us at BiocAsia 2026 on 19 and 20 November at RMIT City Campus, Melbourne, for workshops, talks and a chance to meet others working with R and Bioconductor. Registration is open, and we’d love to see you there.

To keep in touch between events, join the BiocAsia channel on Bioconductor Zulip. It’s a place to hear about future seminars, discuss watch parties and share ideas for other activities.

A big thank you to the speakers, organisers and everyone giving their time to make these activities happen.

© 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 volunteer at my local junior AFL club and my role as club registrar grants access to a sports management platform called PlayHQ. It’s a rich source of data so for the end of the 2026 season, I thought it would be fun to use some of that data to create a season summary in an informational video format. You can see the result in the video, above.

This post is not a how-to, just a brief summary of the tools that I used.

The input data: CSV files downloaded from the available reports in PlayHQ. These include data describing player registration, games played and venues.

A lot of R users seem to be moving towards the Positron IDE and Quarto to generate presentation material. I’m still using RStudio and RMarkdown, simply because that’s what I know best right now. The key packages for this project: ggplot2 of course; waffle to show player counts; sf, leaflet, geosphere, tidyterra, osrm and maptiles for geographical data; wordcloud for, well, word clouds.

I used the RMarkdown file to sketch out each slide and generate one high-quality graphic per slide in PNG format, or in some cases as an animated GIF. Next, I used Keynote to create the final slides from the imported PNG files. I’ve had Keynote sitting largely unopened on every Mac that I’ve owned for years – turns out that it’s pretty good – not least because you can export the final product straight to video. I used a custom resolution of 1170×2532 to match my (iPhone 13) screen, which seems to work OK on most current phone screens.

There’s plenty of free-to-use cheesy corporate music online to choose from, then it’s just a case of fiddling with the slide transitions to sync the music and any animations. Getting the timing right for slides that contain animated GIFs can be a bit tricky – quite honestly I don’t recall what worked best, other than clicking on the item, selecting the “Movie” tab in the right sidebar menu and “adjusting until it works”.

Once the run-through in Keynote looks and sounds good, export to video and it’s ready to share.

ahead::ctxthetaf has no more contextual information about the time series (ahead::contextridge2f still does), but it still works, and can still be used notably with caretS3 (a caret derivative introduced in https://thierrymoudiki.github.io/blog/2026/10/25/r/caretS3-ahead) to forecast time series data. The following code demonstrates how to use caretS3 models for time series forecasting with ahead::ctxthetaf.

list_models <- caret::modelLookup()[caret::modelLookup()[,"forReg"]==TRUE, 1]

n_models <- length(list_models)

pb <- utils::txtProgressBar(min = 0, max = n_models, style = 3)

for (i in seq_len(n_models))
{
  ts_model <- try(ahead::ctxthetaf(y=fdeaths, 
                                   theta = 1,
                                   h=25L, 
                                   fit_func = function (x, y) caretS3::caretS3(x, y, 
                                                                               method=list_models[i])), 
                  silent = TRUE)
  if (class(ts_model) != "try-error")
  {
    plot(ts_model) 
  } 
  utils::setTxtProgressBar(pb, i)
}

close(pb)

Free AI and Data Courses with 365 Data Science: Unlimited Access until Nov. 20

365 Data Science opens its doors this November 100% for free!

From November 5 to November 20, 2026, enjoy unlimited access to the entire platform, including a wide range of R programming courses and projects.

This offer features expert-led courses, practical projects, and interactive exercises covering diverse topics in AI and data science. With flexible 24/7 access, learners can study at their own pace, making it an ideal opportunity to grow and advance in these fields.

The Initiative’s Sixth and Biggest Year

Now in its sixth year, 365 Data Science reaffirms its dedication to delivering accessible, high-quality education through its annual Free Access Initiative, originally launched in 2020 during the global pandemic. CEO Ned Krastev emphasizes the accelerating demand for AI-driven skills, stating, “The AI and data landscape is evolving faster than ever, creating extraordinary opportunities for career growth for those ready to embrace new technologies.”

The initiative’s reach has grown remarkably: 2025 marked its most successful edition yet, with over 220,000 unique learners from 216 countries participating, accumulating 7.2 million minutes of learning and earning more than 36,000 certificates.

Krastev adds, “Artificial intelligence is reshaping industries at an unprecedented pace. Understanding how AI systems function and how to build, deploy, and integrate them has become vital for anyone pursuing a data-driven career. At 365 Data Science, our mission is to close that gap by equipping learners not only with data literacy but also with hands-on expertise in AI engineering and intelligent agents: the skills that will define the next generation of tech professionals.”

365 Data Science empowers learners to move beyond traditional analytics and step confidently into the era of AI engineering and intelligent agents. The platform equips users with the skills to design, deploy, and collaborate with AI systems capable of reasoning, planning, and acting autonomously.

What’s Included

During this special period, learners gain unlimited access to the entire 365 Data Science platform, a comprehensive destination for mastering data and AI. The platform offers 134 expert-led courses, covering everything from fundamental data skills to advanced topics in machine learning, AI, and AI engineering.

Participants can apply their knowledge through real-world projects that mirror practical work scenarios, gaining hands-on experience in solving meaningful problems. Newly added interactive exercises and guided challenges deepen understanding and strengthen key concepts.

To support career growth, 365 Data Science also provides structured, career-focused learning paths guiding learners step by step from beginner to job-ready professional with a clear roadmap to success in today’s AI-driven world.

Certifications that Open Doors

In today’s fast-changing job market, recognized certifications are essential to stand out. Through this free access initiative, 365 Data Science offers learners the chance to earn industry-recognized certificates at no cost.

These credentials validate practical expertise in data analytics, AI, and machine learning, enhancing employability and building credibility with employers worldwide. By offering verifiable, career-boosting certifications, the initiative bridges the gap between learning and professional achievement by showcasing real-world competence.

R Programming: A Core Focus

365 Data Science acknowledges the central role of R in modern data analysis and visualization, reflected in its carefully curated course selection:

• Introduction to R Programming: A beginner-friendly course that builds a strong foundation in R syntax, data structures, and essential operations—perfect for kick-starting your R journey.
• The Complete Data Visualization Course with Python, R, Tableau, and Excel: A comprehensive program that highlights R’s powerful visualization tools while exploring multiple platforms for data storytelling.
• Data-Driven Business Growth: Learn how to apply R to real-world business cases, adopting a growth mindset to generate meaningful insights and drive measurable impact.

Additionally, learners can explore a variety of complementary courses designed to deepen their analytical toolkit, covering Statistics, Probability, Statistical Tests in Sales and Marketing, Data Preprocessing, and more—helping you master both R and other key programming languages in data science.

Don’t Miss this Opportunity

Free access runs from November 5 to November 20, 2026.

In a world increasingly powered by data and artificial intelligence, staying ahead of the curve is essential. This three-week open-access period from 365 Data Science offers a rare opportunity to invest in your future: whether you’re just getting started, switching careers, or advancing your expertise in AI and data.

Don’t miss your chance to build in-demand skills, earn industry-recognized certificates, and take the next step toward a successful career in data science and AI engineering.

The future belongs to those who prepare for it today: start your journey for free with 365 Data Science.

In 2020 I trained a U-Net on the 2018 Data Science Bowl nuclei with the first version of this package, watched a Dice score climb during training, and finished the post with four predicted masks on the competition’s test images. Those images have no masks. The pictures looked good, and I had no way of knowing where the model was wrong – so I never asked.

This post asks. Nothing in it trains: the model is the published dsbowl-unet from the first post of the series, scored on the images it was held out from. Scoring all of them took about twelve seconds on a processor.

data_dir <- "stage1_train"   # BBBC038v1, unzipped

data_dir = "stage1_train"   # BBBC038v1, unzipped

library(rplatypus)

split <- split_dataset(data_dir, tempfile("splits"),
                       fractions = c(0.8, 0.2), seed = 1)
split
#> platypus split
#>   train        536 samples
#>   validation   134 samples
#>   every sample its own group
#>   seed 1 - the same inputs give this split again

import tempfile
import pyplatypus as pp

split = pp.split_dataset(data_dir, tempfile.mkdtemp(),
                         fractions=(0.8, 0.2), seed=1)
split["samples"]
#> {'train': 536, 'validation': 134, 'test': 0}

spec <- platypus_spec(
  data = segmentation_data(split, colormap = binary_colormap),
  models = list(
    u_net("nuclei", input_shape = c(256, 256), blocks = 4, filters = 16,
          weights = "dsbowl-unet", fit = FALSE,
          metrics = list(metric_dice(include_background = FALSE)))
  )
)
fit <- platypus_fit(spec, verbose = FALSE)

cases <- evaluate_cases(fit)
summary(cases)
#>   metric   n      mean         sd    median       min       max
#> 1   dice 134 0.9205236 0.05322748 0.9313503 0.7239583 0.9860661
#>                                                         worst_case
#> 1 7f38885521586fc6011bef1314a9fb2aa1e4935bd581b2991e1d963395eab770

engine = pp.build_engine(pp.from_dict({
    "task": "semantic_segmentation",
    "data": {
        "train_path": str(split["paths"]["train"]),
        "validation_path": str(split["paths"]["validation"]),
        "colormap": [[0, 0, 0], [255, 255, 255]],
        "mode": "config_file",
    },
    "models": [{
        "name": "nuclei",
        "architecture": "u_net",
        "input_shape": [256, 256],
        "blocks": 4,
        "filters": 16,
        "weights": "dsbowl-unet",
        "fit": False,
        "metrics": [{"name": "dice", "include_background": False}],
    }],
}))
engine.fit()

cases = engine.evaluate_cases("nuclei")
pp.summarise_cases(cases)
#> [{'metric': 'dice', 'n': 134, 'mean': 0.9205236150257623,
#>   'sd': 0.05322747958922138, 'median': 0.9313502907752991,
#>   'min': 0.7239583134651184, 'max': 0.9860661029815674,
#>   'worst_case':
#>     '7f38885521586fc6011bef1314a9fb2aa1e4935bd581b2991e1d963395eab770'}]

One row per image

evaluate() gives one row per model. evaluate_cases() gives one row per image, and that is the table worth reading:

• R
• Python

head(cases[order(cases$dice), ], 5)
#> platypus scores by case (5 rows)
#>                                                              case      dice
#>  7f38885521586fc6011bef1314a9fb2aa1e4935bd581b2991e1d963395eab770 0.7239583
#>  29ea4f6eb4545f43868a9b40a60000426bf8dfd9d062546656a37bd2a2aaf9ec 0.7531800
#>  8f94a80b95a881d0efdec36affc915dca9609f4cba8134c4a91b219d418778aa 0.7541562
#>  a6515d73077866808ad4cb837ecdac33612527b8a1041e82135e40fce2bb9380 0.7588775
#>  2a1a294e21d76efd0399e4eb321b45f44f7510911acd92c988480195c5b4c812 0.7702851

worst = sorted(cases, key=lambda row: row["dice"])[:5]
[(row["case"][:16], round(row["dice"], 4)) for row in worst]
#> [('7f38885521586fc6', 0.724), ('29ea4f6eb4545f43', 0.7532),
#>  ('8f94a80b95a881d0', 0.7542), ('a6515d7307786680', 0.7589),
#>  ('2a1a294e21d76efd', 0.7703)]

Those names say nothing, so give each image a label. The Data Science Bowl mixed three kinds of microscopy, and two numbers tell them apart well enough: whether the picture has colour, and whether its background is light or dark.

• R
• Python

image_kind <- function(case) {
  path <- list.files(file.path(data_dir, case, "images"), full.names = TRUE)
  pixels <- png::readPNG(path[1])[, , 1:3]
  colour <- mean(abs(pixels[, , 1] - pixels[, , 2])) +
    mean(abs(pixels[, , 2] - pixels[, , 3]))
  if (colour > 0.02) return("stained")
  if (mean(pixels) > 0.5) "brightfield" else "fluorescence"
}
cases$kind <- vapply(cases$case, image_kind, character(1))

aggregate(dice ~ kind, cases, function(x) {
  round(c(n = length(x), mean = mean(x), min = min(x)), 3)
})
#>           kind  dice.n dice.mean dice.min
#> 1  brightfield   3.000     0.749    0.724
#> 2 fluorescence 105.000     0.935    0.753
#> 3      stained  26.000     0.883    0.805

from pathlib import Path
import numpy as np
from PIL import Image

def image_kind(case):
    path = next((Path(data_dir) / case / "images").glob("*.png"))
    pixels = np.asarray(Image.open(path).convert("RGB")) / 255
    colour = (np.abs(pixels[..., 0] - pixels[..., 1]).mean()
              + np.abs(pixels[..., 1] - pixels[..., 2]).mean())
    if colour > 0.02:
        return "stained"
    return "brightfield" if pixels.mean() > 0.5 else "fluorescence"

for row in cases:
    row["kind"] = image_kind(row["case"])

kinds = sorted({row["kind"] for row in cases})
dice = {k: [r["dice"] for r in cases if r["kind"] == k] for k in kinds}
{k: (len(d), round(float(np.mean(d)), 3), round(min(d), 3))
 for k, d in dice.items()}
#> {'brightfield': (3, 0.749, 0.724), 'fluorescence': (105, 0.935, 0.753),
#>  'stained': (26, 0.883, 0.805)}

• R
• Python

library(ggplot2)
ggplot(cases, aes(dice, kind)) +
  geom_vline(xintercept = mean(cases$dice), linetype = "dashed",
             colour = "grey50") +
  geom_jitter(height = 0.15, width = 0, alpha = 0.7, size = 2) +
  labs(x = "Dice, per image", y = NULL) +
  theme_minimal(base_size = 13)

import matplotlib.pyplot as plt

rng = np.random.default_rng(1)
fig, ax = plt.subplots(figsize=(7, 3.2))
ax.axvline(np.mean([r["dice"] for r in cases]), linestyle="--", color="grey")
for i, k in enumerate(kinds):
    jitter = rng.uniform(-0.15, 0.15, len(dice[k]))
    ax.scatter(dice[k], i + jitter, alpha=0.7)
ax.set_yticks(range(len(kinds)), kinds)
ax.set_xlabel("Dice, per image")

There are three brightfield images among the 134, and they score 0.72, 0.75 and 0.77 - three of the five worst. Fluorescence averages 0.93 over 105 images, and it is most of what the mean measures. Stained tissue sits between the two, at 0.88.

• R
• Python

predictions <- predict(fit, split = "validation")

files_of <- function(case, what) {
  list.files(file.path(data_dir, case, what), full.names = TRUE)
}
truth_of <- function(case) {
  masks <- lapply(files_of(case, "masks"), function(f) {
    read_masks(f, binary_colormap, size = c(256, 256))[1, , ]
  })
  unite_masks(masks)
}

middle <- which.min(abs(cases$dice - median(cases$dice)))
shown <- c(which(cases$kind == "brightfield"), middle)
images <- read_images(
  vapply(cases$case[shown], function(k) files_of(k, "images")[1], ""),
  size = c(256, 256)
)
truth <- simplify2array(lapply(cases$case[shown], truth_of))

plot_masks(
  images,
  prediction = predictions[shown, , ],
  truth      = aperm(truth, c(3, 1, 2)),
  colormap   = binary_colormap,
  labels     = sprintf("%s, Dice %.2f", cases$kind[shown], cases$dice[shown])
)

predictions = engine.predict("nuclei", "validation")
data = engine.dataset(engine.runs["nuclei"].spec, "validation")

middle = np.median([r["dice"] for r in cases])
shown = [i for i, r in enumerate(cases) if r["kind"] == "brightfield"]
shown.append(min(range(len(cases)),
                 key=lambda i: abs(cases[i]["dice"] - middle)))

pp.plot_masks(
    np.stack([data[i][0] for i in shown]),
    prediction=predictions[shown],
    truth=np.stack([data[i][1] for i in shown]),
    labels=[f"{cases[i]['kind']}, Dice {cases[i]['dice']:.2f}"
            for i in shown],
)

nucleus_area <- function(case) {
  files <- list.files(file.path(data_dir, case, "masks"), full.names = TRUE)
  areas <- vapply(files, function(f) {
    sum(read_masks(f, binary_colormap) == 2)
  }, numeric(1))
  size <- dim(read_masks(files[1], binary_colormap))[2:3]
  mean(areas) * 256^2 / prod(size)
}
cases$area <- vapply(cases$case, nucleus_area, numeric(1))

aggregate(area ~ kind, cases, function(x) round(range(x)))
#>           kind area.1 area.2
#> 1  brightfield     11     17
#> 2 fluorescence     69    957
#> 3      stained    137   1893

def nucleus_area(case):
    files = sorted((Path(data_dir) / case / "masks").glob("*.png"))
    masks = [np.asarray(Image.open(f)) > 0 for f in files]
    height, width = masks[0].shape[:2]
    return np.mean([m.sum() for m in masks]) * 256**2 / (height * width)

for row in cases:
    row["area"] = nucleus_area(row["case"])

areas = {k: [r["area"] for r in cases if r["kind"] == k] for k in kinds}
{k: (round(min(a)), round(max(a))) for k, a in areas.items()}
#> {'brightfield': (11, 17), 'fluorescence': (69, 957),
#>  'stained': (137, 1893)}

summary(evaluate_cases(fit, at = "source"))
#>   metric   n     mean         sd    median       min       max
#> 1   dice 134 0.920808 0.05371718 0.9313503 0.7164674 0.9860661
#>                                                         worst_case
#> 1 7f38885521586fc6011bef1314a9fb2aa1e4935bd581b2991e1d963395eab770

pp.summarise_cases(engine.evaluate_cases("nuclei", at="source"))
#> [{'metric': 'dice', 'n': 134, 'mean': 0.920808045276955,
#>   'sd': 0.053717183966230486, 'median': 0.9313502907752991,
#>   'min': 0.7164673805236816, 'max': 0.9860661029815674,
#>   'worst_case':
#>     '7f38885521586fc6011bef1314a9fb2aa1e4935bd581b2991e1d963395eab770'}]

From postal points to areas using Voronoi polygons

13 min.

I recently published my project Postal codes of Spain, with estimated postal areas, an interactive map and GeoPackage downloads. But there is something we have not looked at yet: how do those address points become polygons? The answer is… Voronoi polygons.

What is a Voronoi polygon?

Imagine a few shops on a flat map. Which one is closest to your house? Ask that question for every location and give each shop its own color. The result is a Voronoi diagram, also known as a Voronoi tessellation, Voronoi decomposition, Voronoi partition or Dirichlet tessellation. This technique divides a plane into regions based on the distance to a given set of points.

In our example, each shop is a site, or generating point. Its Voronoi cell, also called a Thiessen polygon, contains the locations at least as close to that site as to any other. Along a shared edge, two sites are equally close. Where several cells meet, there may be a tie between more sites. With ordinary straight-line distance, the edges follow perpendicular bisectors between sites. Some cells extend indefinitely, so maps usually show them clipped to a study area.

The name honors Georgy Voronoi, who studied these partitions in the early twentieth century. The idea has earlier roots, including Dirichlet’s nineteenth-century work (hence the name Dirichlet tessellation). Pokojski and Pokojska (2018) review the history and geographical applications.

These partitions are useful for exploring proximity to shops or public services, assigning areas to weather stations and studying spatial patterns. For example, Thiessen polygons can supply area weights when estimating rainfall from rain gauges. The underlying question is always about the chosen distance: the nearest hospital in a straight line may not be the quickest one to reach.

A tiny example in R

Let’s start with five points in an imaginary square, including one at its center, using just sf and ggplot2. The coordinates have arbitrary units and no geographical coordinate reference system.

library(sf)
library(ggplot2)

# Define a common theme for the maps.
background <- "#faf8f2"
ink <- "#203b47"

theme_set(
  theme_void() +
    theme(
      plot.background = element_rect(fill = background, colour = NA),
      panel.background = element_rect(fill = background, colour = NA),
      plot.title = element_text(face = "bold", colour = ink),
      plot.subtitle = element_text(colour = ink),
      plot.caption = element_text(colour = "grey40", hjust = 0),
      legend.position = "bottom",
      plot.margin = margin(10, 10, 10, 10)
    )
)

sites <- data.frame(
  site = c("A", "B", "C", "D", "E"),
  x = c(2, 7, 3, 8, 5),
  y = c(2, 3, 8, 7, 5)
)

toy_points <- st_as_sf(sites, coords = c("x", "y"))

study_area <- st_as_sfc(st_bbox(c(xmin = 0, ymin = 0, xmax = 10, ymax = 10)))

toy_cells <- toy_points |>
  st_union() |>
  st_voronoi(envelope = study_area) |>
  st_collection_extract("POLYGON") |>
  st_intersection(study_area)

ggplot() +
  geom_sf(
    data = toy_cells,
    fill = c("#a8dadc", "#f4c89a", "#c6b8df", "#bfd59b", "#f2b5c4"),
    colour = "white", linewidth = 1
  ) +
  geom_sf(data = toy_points, colour = "#203b47", size = 3) +
  geom_text(
    data = sites, aes(x, y, label = site),
    nudge_y = 0.5, colour = "#203b47", fontface = "bold"
  ) +
  coord_sf(datum = NA) +
  labs(
    title = "Five points, five Voronoi cells",
    subtitle = "Each location belongs to the closest site"
  )

Five points, five cells. The point E at (5, 5) has a bounded cell surrounded by the other four. The outer cells extend beyond the square and are clipped to its edges. st_union() collects the points into one geometry, st_voronoi() builds the diagram and st_intersection() clips it to our square. The colors help us tell the cells apart. We have not assigned any postal codes yet.

The diagram knows nothing about roads, rivers or administrative boundaries. As we move to postal points, we will need to make our choices about those boundaries explicit.

Before we start: what are we estimating?

The Spanish project uses CartoCiudad address points with their postal codes. It follows Goerlich (2022), who develops an idea presented by Sánchez (2020). The updated online version describes the method in detail.

For Slovenia, we will use GISCO’s postal points, retrieved with giscoR. There is an important difference here: a representative point for a code tells us much less than thousands of addresses with that code. Detailed addresses help the partition follow their distribution more closely, although proximity alone still cannot reproduce postal delivery decisions.

We have not established that Slovenian postal areas should follow municipal boundaries. Imposing those boundaries is a modeling choice for this exercise, inspired by the Spanish workflow, rather than a verified property of the Slovenian postal system.

This post explains the workflow. The polygons we create are not official postal boundaries, and this exercise does not validate the GISCO postal data.

Getting the data

First, let’s get the municipal boundaries and postal points. We use the 2024 release of both datasets.

library(giscoR)
library(sf)
library(dplyr)

municipalities <- gisco_get_lau(
  year = 2024, country = "SI", epsg = 3035
) |>
  select(GISCO_ID, LAU_NAME) |>
  st_make_valid()

raw_points <- gisco_get_postalcodes(
  year = 2024, country = "SI", epsg = 3035
)

We work in ETRS89 / LAEA Europe (EPSG:3035), a projected coordinate system with coordinates in meters. This makes the 200 m and 50 m margins used later meaningful. Voronoi uses Euclidean distances in this plane, rather than road distances. The projection does not preserve every distance exactly either.

Before we start drawing, what have we actually downloaded?

postal_points <- raw_points |>
  mutate(POSTCODE = trimws(as.character(POSTCODE))) |>
  filter(!is.na(POSTCODE), POSTCODE != "") |>
  select(POSTCODE)

code_counts <- postal_points |>
  st_drop_geometry() |>
  count(POSTCODE, name = "points_per_code")

knitr::kable(data.frame(
  Indicator = c(
    "Municipalities", "Downloaded points", "Points with a postal code",
    "Distinct postal codes", "Codes represented by one point"
  ),
  Value = c(
    nrow(municipalities), nrow(raw_points), nrow(postal_points),
    nrow(code_counts), sum(code_counts$points_per_code == 1)
  )
))

These counts give us a first look at the file, but they are not a quality check. They do not establish whether every code is present or whether the coordinates are correct. With one point per code, we will be drawing areas around representatives rather than estimating them from individual addresses.

Two points at exactly the same coordinates cannot generate separate Voronoi cells. We sort by postal code and keep the first point at each location, as in the original workflow. This is a reproducible tie-breaking rule, not evidence that the retained code is correct. The counts below show how many records we remove and how many locations have conflicting codes.

postal_points <- postal_points |> arrange(POSTCODE)

conflicts <- postal_points |>
  group_by(geometry) |>
  summarise(codes = n_distinct(POSTCODE), .groups = "drop") |>
  filter(codes > 1)

duplicates <- duplicated(postal_points$geometry)

knitr::kable(data.frame(
  Indicator = c("Duplicate locations removed", "Locations with conflicting codes"),
  Value = c(sum(duplicates), nrow(conflicts))
))

postal_points <- postal_points[!duplicates, ]

Starting with Ljubljana

Let’s pick Ljubljana and keep the points inside it. We use st_within() so that a point exactly on a shared boundary is not assigned to both neighboring municipalities.

city <- municipalities |> filter(LAU_NAME == "Ljubljana")

city_points <- st_filter(postal_points, city, .predicate = st_within)

base_map <- ggplot() +
  geom_sf(data = city, fill = "#e8e4d8", colour = ink, linewidth = 0.4)

gisco_attr <- gisco_attributions()

base_map +
  geom_sf(data = city_points, colour = "#bd5137", size = 2) +
  labs(
    title = "Ljubljana: the starting point",
    subtitle = "A municipal boundary and points with postal codes",
    caption = gisco_attr
  )

Building the cells

Each cell contains locations at least as close to its generating point as to any other point. Cell boundaries occur where distances to two points are equal. For this first demonstration, we keep all interior points, without the margins used in the national workflow.

cells <- city_points |>
  st_union() |>
  st_voronoi(envelope = st_as_sfc(st_bbox(city))) |>
  st_collection_extract("POLYGON") |>
  st_sf(geometry = _)

base_map +
  geom_sf(data = cells, fill = NA, colour = "#bd5137", linewidth = 0.5) +
  geom_sf(data = city_points, colour = ink, size = 2) +
  labs(
    title = "An area of proximity for every point",
    subtitle = "Cells still extend beyond the municipal boundary",
    caption = "Own elaboration using data from Eurostat / GISCO, 2024"
  )

st_union() gathers the points so we build one diagram. One detail to keep in mind: envelope supplies a working extent, but does not clip the cells to the municipality. We need an intersection for that.

Clipping and recovering postal codes

It would be convenient if cells came out in the same order as the input points, but we should not rely on that. Instead, we recover each cell’s code with a spatial relationship before clipping. Each cell must contain exactly one generating point.

point_matches <- st_intersects(cells, city_points)

cells$POSTCODE <- city_points$POSTCODE[unlist(point_matches)]

clipped_geometry <- st_intersection(st_geometry(cells), st_geometry(city))

city_cells <- st_sf(
  POSTCODE = cells$POSTCODE[attr(clipped_geometry, "idx")[, 1]],
  geometry = clipped_geometry
)

city_areas <- city_cells |>
  group_by(POSTCODE) |>
  summarise(.groups = "drop") |>
  st_make_valid()

city_colours <- setNames(
  hcl.colors(nrow(city_areas), palette = "Dark 3"),
  sort(city_areas$POSTCODE)
)
city_labels <- st_sf(
  POSTCODE = city_areas$POSTCODE,
  geometry = st_point_on_surface(st_geometry(city_areas))
)

ggplot(city_areas) +
  geom_sf(aes(fill = POSTCODE), colour = background, linewidth = 0.6) +
  geom_sf(
    data = city_points, shape = 21, fill = ink,
    colour = background, size = 2.4
  ) +
  geom_sf_text(
    data = city_labels,
    aes(label = POSTCODE), colour = ink, fontface = "bold", size = 3.5
  ) +
  scale_fill_manual(values = city_colours, guide = "none") +
  coord_sf(datum = NA) +
  labs(
    title = "From points to estimated areas",
    subtitle = "Cells clipped and merged by postal code",
    caption = "Teaching example · Not official postal boundaries"
  )

We intersect the geometry columns, then use the intersection’s idx matrix to recover each cell’s postal code. This keeps the attribute assignment explicit. Label positions are also calculated from geometries alone.

Grouping by postal code with summarise() merges the geometries. With many addresses sharing a code, their cells would dissolve into one geometry, possibly with several separate parts.

From one municipality to the whole country

Ljubljana is done. Now let’s repeat the operation for each component of each municipality. Municipalities can have separate pieces: normalizing to MULTIPOLYGON and then casting to POLYGON lets us handle them independently, without automatically extending a code to a piece with no points.

The municipality ID and name apply to every component, so we mark those attributes as spatially constant before splitting the geometries.

components <- municipalities |>
  st_set_agr("constant") |>
  st_cast("MULTIPOLYGON") |>
  st_cast("POLYGON") |>
  mutate(component_id = row_number())

Why use an inward margin?

In the Spanish project, we first select addresses within an inward margin of 200 m. If none are found, we try 50 m. This reduces the influence of points close to a boundary, where small mismatches between layers may affect their assignment.

There is no universal rule behind these numbers. Goerlich describes a 20 m margin. The 200 m and 50 m choices come from my implementation. We keep them here to explain that workflow, although sparse postal points make it easier to exclude useful information. Comparing other margins would be a useful next step.

select_points <- function(component, postal_points) {
  for (margin in c(200, 50)) {
    interior <- st_buffer(st_geometry(component), -margin)

    if (all(st_is_empty(interior))) next

    selected_points <- st_filter(postal_points, interior,
      .predicate = st_within
    )

    if (nrow(selected_points) > 0) {
      return(selected_points)
    }
  }

  postal_points[0, ]
}

The margin only selects points: the resulting areas cover the original component, not the smaller buffered geometry. If both attempts fail, we leave the component unassigned. We do not borrow the nearest point from outside the municipality.

Three cases, one function

There are three cases to handle. With no points, we keep the component and assign NA. If all selected points share a code, we assign it to the whole component. We only need Voronoi when several codes are present: no need to split an area that already has a single assignment.

build_areas <- function(component, postal_points) {
  selected_points <- select_points(component, postal_points)

  codes <- unique(selected_points$POSTCODE)

  if (length(codes) <= 1L) {
    component$POSTCODE <- if (length(codes)) codes else NA_character_
    return(component)
  }

  cells <- selected_points |>
    st_union() |>
    st_voronoi(envelope = st_as_sfc(st_bbox(component))) |>
    st_collection_extract("POLYGON") |>
    st_sf(geometry = _)

  point_matches <- st_intersects(cells, selected_points)

  cells$POSTCODE <- selected_points$POSTCODE[unlist(point_matches)]

  clipped_geometry <- st_intersection(
    st_geometry(cells),
    st_geometry(component)
  )

  result <- st_sf(
    POSTCODE = cells$POSTCODE[attr(clipped_geometry, "idx")[, 1]],
    geometry = clipped_geometry
  ) |>
    group_by(POSTCODE) |>
    summarise(.groups = "drop") |>
    st_make_valid()

  result$GISCO_ID <- component$GISCO_ID

  result$LAU_NAME <- component$LAU_NAME

  result$component_id <- component$component_id

  result
}

result <- lapply(seq_len(nrow(components)), function(i) {
  build_areas(components[i, ], postal_points)
}) |>
  bind_rows()

Handling unassigned areas

Can we fill any of the gaps? A component without points may belong to a municipality whose other components have information. We recover its code only if the municipality has one distinct recovered code. This is an imputation from the available data, not proof that the municipality really has just one code.

unique_codes <- result |>
  st_drop_geometry() |>
  filter(!is.na(POSTCODE)) |>
  distinct(GISCO_ID, POSTCODE) |>
  add_count(GISCO_ID) |>
  filter(n == 1)

unassigned_before <- sum(is.na(result$POSTCODE))

index <- match(result$GISCO_ID, unique_codes$GISCO_ID)

recoverable <- is.na(result$POSTCODE) & !is.na(index)

result$POSTCODE[recoverable] <- unique_codes$POSTCODE[index[recoverable]]

Before drawing the final map, let’s tidy up the geometries and keep only surfaces. Repairs can produce collections containing residual lines or points, which do not belong in an area layer. Dissolving by postal code also removes municipal boundaries between areas sharing that code.

keep_surfaces <- function(x) {
  geometry <- st_make_valid(st_geometry(x))

  geometry <- lapply(geometry, function(part) {
    if (!inherits(part, c("POLYGON", "MULTIPOLYGON", "GEOMETRYCOLLECTION"))) {
      return(st_multipolygon())
    }

    surfaces <- suppressWarnings(st_collection_extract(st_sfc(part), "POLYGON"))

    if (length(surfaces) == 0L || all(st_is_empty(surfaces))) {
      return(st_multipolygon())
    }

    st_cast(st_union(surfaces), "MULTIPOLYGON", warn = FALSE)[[1]]
  })

  st_geometry(x) <- st_sfc(geometry, crs = st_crs(x))

  x[!st_is_empty(x), ]
}

unassigned <- result |>
  filter(is.na(POSTCODE)) |>
  keep_surfaces()

postal_areas <- result |>
  filter(!is.na(POSTCODE)) |>
  group_by(POSTCODE) |>
  summarise(.groups = "drop") |>
  keep_surfaces() |>
  st_cast("MULTIPOLYGON")

The complete map

For the national overview, let’s color the areas by the first digit of their code. This keeps the map readable, without claiming that we have validated a classification of postal regions. Gray is reserved for unassigned components, with their own legend entry.

postal_colours <- c(
  "1" = "#007F7F", "2" = "#E87500", "3" = "#3B6FB6", "4" = "#A64E9A",
  "5" = "#649A32", "6" = "#D9AF00", "8" = "#B94350", "9" = "#6255A4",
  "Unassigned" = "#727A80"
)

postal_groups <- sort(unique(substr(postal_areas$POSTCODE, 1, 1)))

ggplot() +
  geom_sf(data = municipalities, fill = "#d5d2ca", colour = NA) +
  geom_sf(
    data = postal_areas,
    aes(fill = substr(POSTCODE, 1, 1)), colour = background, linewidth = 0.15
  ) +
  geom_sf(
    data = unassigned, aes(fill = "Unassigned"),
    colour = background, linewidth = 0.15
  ) +
  scale_fill_manual(
    values = postal_colours, name = "Postal code group",
    breaks = c(postal_groups, "Unassigned"),
    labels = c(paste0(postal_groups, "xxx"), "Unassigned")
  ) +
  coord_sf(datum = NA) +
  labs(
    title = "Slovenia, point by point",
    subtitle = "Proximity areas clipped to municipal boundaries",
    caption = paste("Eurostat / GISCO, 2024 · Own elaboration",
      "Dark grey: unassigned areas · Estimated boundaries, not official",
      sep = "\n"
    )
  )

Checking the result

The map is ready, but let’s check the geometry before going further. Are the polygons valid? Do assigned and unassigned areas cover the territory? Are there any appreciable overlaps? We allow one square meter for numerical differences in the geometry operations.

territory <- st_union(st_geometry(municipalities))

coverage <- st_union(c(st_geometry(postal_areas), st_geometry(unassigned)))

difference_m2 <- as.numeric(sum(st_area(st_sym_difference(territory, coverage))))

overlap_m2 <- as.numeric(
  sum(st_area(postal_areas)) + sum(st_area(unassigned)) - sum(st_area(coverage))
)

knitr::kable(data.frame(
  Indicator = c(
    "Codes with an estimated area", "Components assigned by municipal recovery",
    "Unassigned records", "Unassigned area (km²)",
    "Coverage difference (m²)", "Overlap (m²)"
  ),
  Value = round(c(
    nrow(postal_areas), sum(recoverable), nrow(unassigned),
    as.numeric(sum(st_area(unassigned))) / 1e6, difference_m2,
    overlap_m2
  ), 4)
))

These checks verify the geometry, not postal accuracy. We could cover an entire country with valid polygons and still have a poor approximation of its delivery areas.

Back to Spain

We now have the full sequence: points, proximity areas, clipping and merging by postal code. The same steps work with much more detailed addresses, such as the CartoCiudad data in my Spanish postal codes project.

The input data matters as much as the algorithm. With one representative point per code, Voronoi mainly maps distances between representatives. With address points, it can better approximate where each code’s addresses are concentrated. It still knows nothing about delivery routes, physical barriers or official boundaries.

We have not compared these results with an independent postal source for Slovenia. Before using them for analysis, we would need to review coverage, coordinates, the relationship with municipalities and sensitivity to the margins. For now, we have a working example that shows how the polygons are built, one step at a time.

References

• Pokojski, W. & Pokojska, P. (2018). Voronoi diagrams – inventor, method, applications. Polish Cartographical Review, 50(3), 141–150. https://doi.org/10.2478/pcr-2018-0009.
• Goerlich, F. J. (2022). Elaboración de un mapa de códigos postales de España con recursos libres. Cómo evitar pagar por disponer de información de referencia. Working Papers Ivie, no. 2022-3. València: Ivie. https://doi.org/10.12842/WPIVIE_0322. Updated online version: https://www.uv.es/goerlich/Ivie/CodPost.html.
• Sánchez, J. I. (2020, November 13). “Así hice el mapa de los códigos postales de España con SIG y datos abiertos.” NoSoloSIG. https://www.nosolosig.com/articulos/asi-hice-el-mapa-de-los-codigos-postales-de-espana-con-sig-y-datos-abiertos.
• Hernangómez, D. (2026). Códigos postales de España [Dataset]. https://doi.org/10.5281/zenodo.23698418. https://dieghernan.github.io/codigospostales/.

I recall when the R shiny framework first came out. It suddenly got relatively easy for us ecological modellers and environmental statisticians to turn our R code into web apps. This improved accessibility of complex analyses to more general audiences.

You started seeing ‘web apps’ popping up as outcomes of all sorts of projects. (Here’s one great example, Biological Ocean Observer let’s you analyse timeseries of plankton data from Australia).

The barrier to doing these apps before was that most of the modelling community didn’t know how to code in the languages of the web. An R based framework for web apps made that possible.

But there was still a downside. With R shiny you still need to host the R code on a server somewhere. This requires either a free, but limited account, or paying for compute.

The ideal situation is to have your app written in Java Script and html. These languages run web pages. The key difference from R Shiny is the computations run in the web browser, on the viewer’s computer. This means you can host a web page cheaply or for free (e.g. as a github page). No server for computing R required.

Back in 2018, we paid a JS programmer to turn my R code into a web app for one of our projects (on projecting carbon emissions from mangrove deforestation). This app become very useful for people assessing potential outcomes of carbon projects.

Well AI coding agents have removed that barrier coding in JS.

I don’t think the huge potential of this change is widely appreciated in the environmental data analysis and modelling community yet. So I wanted to call it out here.

The tech community itself is only just starting to appreciate the impacts of this change. Here’s a nice explanation of how this new power of running apps in the browser is transformative for the tech industry. (Check out tinylittleoffice, a browser based office suite that saves the document as a URL you can copy and paste, or Battle Drones a multiplayer web game that connects two devices through sound and microphones!).

This change is also going to be transformative for modelling and communication of science.

I got inspired by examples like those in the post above, and set-off on my own AI coding experiment. The result is FANGS, a browser based Bayesian computation engine.

Brief background on my motivation here. I’m a huge fan of Bayesian statistics. More people should be using it. But a barrier to sharing that enthusiasm is the complexity of installing many Bayesian R packages on other people’s computers.

So I built (vibe coded) FANGS as an educational tool. It means my students can write and run simple Bayesian models, directly in their browsers. This means that in a 3 hour lesson I can: cover an intro to Bayesian concepts, have students write models in R like syntax, run the models and interpret the results.

The idea is that is enough to get them motivated to then install the proper R software. The use of model syntax that emulates the serious Bayesian R packages is also a stepping stone to the more powerful R software.

I’ll write up the development pipeline in a later post. But I just want to say it blew me away how well this experiment with vibe coding worked. FANGS gets you the right answers and runs faster on small teaching scale problems than the equivalent R packages.

I think this new power will revolutionize scientific communication and education, through improving accessibility of many types of data and science.

What does this new wave of scientific app development what do you need?

On a personal level, curiosity about what is possible in web browsers. Once you know something is possible its easy to get an AI agent to make it. Once I learned about saving data in URLs from the tinylittleoffice example, it was easy to ask Claude to implement that feature in FANGS (see the copy share link button).

On a community level, translating more of the most common R packages to Java Script libraries will help. I actually haven’t fully explored what is available in JS yet. But for example, efficient JS packages for mapping exist, extending those for complex GIS operations will extend the potential of our apps.

Standards for validation. When you vibe code a package, you need to know it is accurate. I’ll come to how to do this when I write up the FANGS development process.

Finally, a bottleneck with apps has always been user experience (UX). For any scientific application a key question for me has always been: how do we make sure users don’t misinterpret the information this app provides? I’m not sure AI can solve that problem yet. It requires thoughtful design and lots of testing.

Let me know cool ideas you have or have seen, I’m keen to see more examples.

This post is part of our SAS to R series:

• Part 1: Is SAS still used, and is it worth keeping?
• Part 2: Why R, and where your migration budget should go (this post)
• Part 3: coming soon

If you’re moving regulated analysis from SAS to R, two things are worth knowing up front. R is a choice you can defend: regulators don’t tell you which software to use, and R has already been through FDA review in public pilot submissions. And the language itself is the cheap part. Most of the cost and risk comes from six other places, and most of those are about planning, not programming.

Part one asked whether SAS is still worth keeping. This post covers why R, and where to put your budget and attention once you commit.

Why R is a defensible choice

Regulators ask for evidence, not a vendor

The FDA’s Statistical Software Clarifying Statement says the agency does not require any specific software for statistical analysis. It asks for the software to be documented in the submission, including version and build, and ICH E9, the international guideline on statistics in clinical trials, asks for documented testing.

The business implication is that the obligation stays the same whichever language you choose. What changes is who does the work. With SAS, much of the evidence arrives with the licence. With R, your organisation produces it.

There is public precedent to show your risk committee

The R Consortium’s R Submissions Working Group has taken R-based tables, figures, standard analysis datasets and an interactive app through FDA review in a series of public pilots, with the code and lessons published. That is not a guarantee for every submission, but it is the most useful thing you can put in front of a risk committee.

Hex stickers and pilot details: R Consortium R Submissions Working Group.

The clinical tooling already exists

Most of what a clinical SAS team does has a maintained R equivalent, much of it built by pharma companies working together in the pharmaverse. You are not funding the tooling, only the move.

One finding comes up repeatedly. A lot of “we have to keep SAS” turns out to mean “we have to produce SAS transport files”, the format regulators expect submitted data in. That format is a published specification, not a licence, and R writes it.

For your team: admiral and sdtm.oak build ADaM and SDTM datasets, mmrm reproduces SAS mixed models, and xportr and haven write and read v5 transport files.

What about Python?

For statistical analysis and regulated reporting, R is usually the closer match to the work SAS did. Python tends to be the better answer for machine learning engineering and production services. Many organisations run both, on shared infrastructure and review processes, so choosing R for statistics does not lock in the rest of the stack.

Where the cost and risk come from

Teams often budget a SAS to R migration as a code rewrite. The rewrite is the smallest of the seven layers below. The six above it are where timelines slip and budgets get cut.

1. The two systems will disagree, and both will be right

SAS and R make different default choices for the same statistical method. For example, they handle tied survival times differently. Run the same analysis in both and the numbers can differ in the last decimal places, with neither tool wrong.

The cost comes when nobody notices until validation, and nobody recorded which choice the analysis plan meant. The fix is cheap if you do it early: add the method choices to each analysis specification before conversion starts. That turns an open-ended rewrite into a bounded validation exercise.

For your team: the cross-industry CAMIS project documents these differences method by method, including survival analysis in R’s survival package and rounding.

2. Data does not mean quite the same thing

The two languages treat missing values and dates differently. A filter carried over from SAS can quietly drop rows in R, changing the denominator. The result is a table that is nearly right, which is harder to catch than one that is obviously wrong. Budget for checking the data, not just converting the code.

For your team: SAS sorts missing numerics below every value and has 28 kinds of missing, while R has one NA. Date origins also differ (1960 against 1970), which haven handles on import.

3. Your macro library holds rules nobody wrote down

As we covered in part one, a mature macro library holds business rules that exist nowhere else. A line-by-line conversion reproduces the code and loses the reasoning. Recovering that intent is analysis work, and it needs people who know the SAS estate. In mature estates it can be one of the longest parts of the migration, so plan for it explicitly.

4. Reproducibility becomes your responsibility

SAS gives you one vendor managing one product. R spreads that across R itself and the add-on packages (community-built libraries of functions) your analyses depend on. Being able to rerun an analysis years later is now something your team owns. The tooling is mature, and it should be mandatory from day one rather than added later.

For your team: pin dependencies with renv, serve packages from a dated snapshot with Posit Package Manager (which we wrote about alongside diffify), and containerise anything long-lived. In the UK public sector, the same idea goes by the name reproducible analytical pipelines.

5. Proving it works is a separate project

Getting the right answer and being able to prove you will keep getting it are different projects, and only the first one is fun. Budget them separately.

The R Validation Hub set out the widely used approach: assess each package by risk rather than validating everything to the same depth. That keeps the effort proportionate. We build on the same idea in Litmus, which scores packages at function level.

For your team: how the scoring works for validation and code quality, the shorter argument in Should I Use Your R Package?, unit testing with testthat, and validating Shiny apps in regulated environments.

6. Dual running is the cost that gets cut

For the whole parallel period you pay for both systems, both sets of people and both sets of sign-off. This is what de-risks the migration, and it is reliably the first thing trimmed when the budget tightens. Protect it.

Agree in advance how close the two sets of results need to be, and why. Then check automatically rather than by eye, so a difference shows up the day it appears.

For your team: diffdf was built for this comparison, and arsenal::comparedf and waldo cover adjacent cases. Run them in CI.

What to do before converting any code

1. Inventory the macro library and the business rules inside it.
2. Add method choices to every analysis specification.
3. Decide your package validation approach.
4. Budget dual running as a protected line.

Each of these is manageable once it has been named and written down. That is what an assessment produces before a single line of code is converted.

If you are planning a migration, get in touch. If your team is earlier than that, our Why Use R? course is written for people arriving from SAS and SPSS.

Frequently asked questions

Is R accepted by regulators?

The FDA does not require any specific software, provided the software used is documented with version and build and its testing is documented. The R Consortium’s pilots have taken R-based tables, figures, standard analysis datasets and an interactive app through FDA review.

Will R give exactly the same numbers as SAS?

Not by default, and that is normal. The two make different default choices for some methods, so results can differ slightly while both are correct. The fix is to state those choices in the analysis specification rather than hunt for a bug.

References

• US FDA, Statistical Software Clarifying Statement
• US FDA, Study Data Technical Conformance Guide
• R Consortium, R Submissions Working Group and 2026 plans and 2025 progress
• PHUSE, PSI, ASA and IASCT, CAMIS: Comparing Analysis Method Implementations in Software
• CAMIS, R versus SAS: Kaplan-Meier and Cox proportional hazards and R versus SAS rounding
• R Validation Hub, A Risk-based Approach for Assessing R Package Accuracy within a Validated Infrastructure
• The R Foundation, R: Regulatory Compliance and Validation Issues
• CRAN, mmrm for mixed models for repeated measures
• CRAN, admiral and sdtm.oak for ADaM and SDTM builds
• CRAN, xportr for writing v5 transport files

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

As part of rOpenSci’s 15th anniversary celebrations, Community Manager Yani sits down with community member Mo to talk about how a first software review turned into years of involvement across R-Universe, the Champions Program, the Slack community, and core package work. Mo reflects on what has kept her engaged and also offers advice for anyone just discovering this community.

• How did you first end up in rOpenSci’s orbit, what was the very first interaction or moment that pulled you in?

I actually don’t recall. I think it was a natural evolution of my general career progression using R as a scientist and slowly transitioning into a research software engineer. Still working in academia, and mostly on tooling, analysis, and data management using R, rOpenSci felt like a very natural community to get involved in.

I believe I was first introduced to it through Maëlle Salmon, who was/is a fellow RLadies+ global team member. My first interaction may have been as a software reviewer for a package, which Maëlle asked me if I was interested in doing.

• Over the years, you’ve shown up in a lot of different ways: R-Universe, the Champions Program, an active participant on the Slack, and a contributor to our core packages. Looking back, how would you describe the arc of your involvement?

It’s been such an easy community to be involved in. I’ve felt that every suggestion, discussion, and help request has been met with kindness, collegiality and with a true sense of people wanting to share their knowledge. Seeing this [atmosphere] in the Slack has lowered any barrier I had for asking “stupid” questions, and always felt like any contribution has been genuinely valued. I’ve also felt very humble in how some interactions have altered my own perception of myself as an Research Software Engineer (RSE). I remember in particular when you asked me if I would consider being a Champion’s mentor, and being so humbled. I had actually considered applying as a Champion, but didn’t have the time to set aside for it at that moment. Being asked to be a mentor instead made me shift my own understanding of my expertise level.

• Of all the ways you’ve participated, is there one that feels like “the” rOpenSci experience for you? The one you’d point to first?

I actually think the Slack is the pinnacle of how it is to be part of the community. Which is likely frustrating for anyone not in that Slack to read. While it’s not the most active space, the activity is generally always useful and/or enjoying. It’s like a friendly bubble of R-geeks just wanting to make cool things.

• Across all these different roles, has your sense of what rOpenSci is, or what it’s for, changed over time?

That’s a hard one. I feel like it’s been quite stable over the years I have participated, but there have likely been subtle changes over time that are hard to put a pin in.

• Was there ever a moment you felt surprised by how much you’d ended up doing, or by where it led?

All the time. It’s almost like having a second set of colleagues, that technically don’t work with me (as in at the same company), but are equally important to my work and its robustness. Between the software review, blog posts, collaborative tutorials, Champion mentoring, rOpenSci package tooling development, discussions in Slack (let us not forget the awesome cat pictures in the cat channel), I feel like I have quite immersed myself in the community. But I have gotten much more back that I have contributed, I believe. Not only in the obvious ways visible, like mentor training, learning more about collaborative development etc., but also in getting feedback on my funding applications and on my personal development projects. It’s all been invaluable for my own development.

• How does participating in rOpenSci activities compare to other ways you’ve been involved in open source more broadly?

I’m not someone who has contributed a lot to open source in general. I’m pretty much in my lane with rOpenSci, RLadies+ and The Carpentries, three organizations I feel all emulate pretty much the same core values in how to treat contributors. This has given me the confidence to participate openly and without fear, and has also given me strength when I’ve received some less than courteous responses in other places.

But comparing it to certain Stack overflow and social media interactions I have had outside these communities, it’s for sure a much more wholesome and pleasant experience. Not being afraid of being seen as incompetent for asking what someone else thinks is a basic question, is so important for how a community can keep growing and thriving.

• What does rOpenSci mean to you at this point, fifteen years in?

A lot. As I am sure you can tell from all that I write here. rOpenSci is really my second professional “home”, where I know I can get support and feedback from those willing to take the time (and even if no one has the time, one generally always gets a positive nudge to continue).

• What would you want someone just discovering rOpenSci today to know?

Don’t hesitate to contact and get involved. Start wherever you feel you can; doing a software review as a reviewer is a nice experience and a really good place to start. You get to learn a little about reviewing others’ work, how to work with git and github, and communicate with other professionals. It really is a very good starting point.

What a kind vision of rOpenSci! This post is part of our 15th anniversary celebration series, and we still have more to come!

Following the recent migration of the teal framework to pharmaverse, Roche and Appsilon are building on their existing collaboration to support its continued development. Appsilon will continue contributing developer expertise and will be recognized as a core contributor and copyright holder alongside Roche.

teal helps developers build Shiny applications for exploring clinical trial data. If you’ve worked with it, you’ll know its modular approach: combining data, filters, and analysis modules into interactive applications.

Appsilon has been working with Roche on {teal} development and has been a core contributor to the open-source library. For example, it has co-implemented the new teal design, teal.reporter, and teal.picks, changes that have improved the quality and capabilities of the framework. Formal recognition as a core contributor is a new step forward that formalizes Appsilon’s role in the project, alongside Roche and the wider community.

@online{unardi2026,
  author = {Unardi, Dony and Veríssimo, André and Senan, Oriol},
  title = {Roche and {Appsilon} {Build} on {Their} {Collaboration} to
    {Support} Teal’s {Continued} {Development}},
  date = {2026-11-04},
  url = {https://pharmaverse.github.io/blog/posts/2026-11-04-roche-appsilon-teal/roche-appsilon-teal.html},
  langid = {en}
}

En 1983, Lisanne Bainbridge, psychologue à l’University College London, publie un article de cinq pages devenu un classique : « Les ironies de l’automatisation ». Elle y observe que le concepteur qui cherche à éliminer l’opérateur lui laisse encore les tâches qu’il ne voit pas comment automatiser. Elle pensait aux procédés industriels. Sa phrase décrit tout aussi bien les équipes qui confient aujourd’hui une partie de leur travail à des agents d’IA.

Un dossier de projet sait chiffrer le temps que l’IA fera gagner. Il sait beaucoup moins bien chiffrer ce qui reste. Or ce qui reste n’est pas un résidu : c’est de la surveillance, de la vérification et du jugement, c’est-à-dire la part du travail où l’erreur de confiance est la plus probable. Automatiser ne retire pas le travail humain, cela en change la nature. Ce changement a un coût, mais ce coût n’est pas une fatalité : il dépend de ce qu’on retire à l’humain et de ce qu’on lui laisse.

Dans un précédent article, je montrais que le coût d’une tâche réussie dépend du travail humain qui rattrape les échecs. Celui-ci pose la question suivante : que fait l’humain quand l’IA réussit ?

Surveiller est un travail

Bainbridge en tirait une conséquence gênante : l’opérateur peut se retrouver avec « un assortiment arbitraire de tâches », celles que la machine ne sait pas faire. Il hérite aussi de la surveillance d’un système où il ne se passe presque rien.

On a longtemps tenu la surveillance pour une tâche peu exigeante. La recherche sur la vigilance dit le contraire. En 2008, une revue de Warm, Parasuraman et Matthews réunit des mesures comportementales, neurales et subjectives. Elles convergent : la vigilance « demande un travail mental intense et elle est stressante ». Sa charge de travail est élevée, et ces tâches réduisent l’engagement tout en augmentant la détresse. Surveiller n’est pas s’ennuyer : c’est un effort.

Ces résultats viennent d’études expérimentales. Ils posent pourtant une question simple à tout projet d’automatisation : qui fournit cet effort, et dans quelles conditions ?

Le travail glisse vers la vérification

L’IA générative rend cette question quotidienne. En 2025, des chercheurs de Microsoft Research et de l’université Carnegie Mellon interrogent 319 travailleurs du savoir. Selon leurs réponses, l’IA générative déplace la pensée critique « vers la vérification de l’information, l’intégration des réponses et la tenue de la tâche ». On rédige moins, on contrôle davantage. Ces résultats sont déclarés par les participants : ils disent ce que les personnes perçoivent, pas ce qu’elles font.

Ce qu’elles font, Xingqi Maggie Ye et Aruna Ranganathan, de Berkeley Haas, l’ont observé pendant huit mois dans une entreprise technologique de 200 personnes dont les salariés avaient largement accès aux outils d’IA générative. Leur étude ethnographique s’appuie aussi sur plus de 40 entretiens. Le temps gagné a été réinvesti dans le travail : les salariés « travaillaient plus vite, prenaient en charge un périmètre plus large et étendaient leur travail à davantage d’heures de la journée, souvent sans qu’on le leur demande ». Ye et Ranganathan notent que cette expansion semble souvent venir des salariés eux-mêmes. Une seule entreprise, des travaux encore en cours : c’est un signal, pas une loi.

Reste la perception. En 2025, l’organisme de recherche METR mène un essai randomisé auprès de 16 développeurs expérimentés de logiciels libres, sur 246 tâches de projets qu’ils connaissent depuis cinq ans en moyenne. Avec l’IA, ils mettent 19 % de temps en plus. Après l’étude, ils estiment pourtant que l’IA leur a fait gagner 20 % de temps. L’effectif est petit et les outils sont ceux de février à juin 2025. L’écart entre l’impression et la mesure devrait cependant alerter quiconque évalue un projet sur les seuls retours de ses utilisateurs.

L’effet « bûcheron »

Ce déplacement ne serait qu’une affaire d’emploi du temps s’il épargnait la qualité du jugement. Il ne l’épargne pas.

En 2014, une méta-analyse de Linda Onnasch et de ses collègues rassemble 18 expériences. Plus le degré d’automatisation augmente, meilleure est la performance de routine ; tant que l’automate fonctionne, la charge baisse aussi, mais plus faiblement. Quand il tombe en panne, en revanche, la performance et la conscience de la situation se dégradent. Les auteurs parlent d’effet « bûcheron » : plus les arbres sont hauts, plus ils tombent loin. Dans toutes les études sauf quatre, les participants étaient des étudiants.

L’expertise ne protège pas. En 2023, la revue Radiology publie une expérience menée auprès de 27 radiologues qui lisent des mammographies avec l’aide d’un système présenté comme une IA. En moyenne, les plus expérimentés classent correctement 82,3 % des clichés quand le système a raison, et 45,5 % quand il se trompe. Chez les moins expérimentés, la justesse passe de 79,7 % à 19,8 %.

Enfin, une revue systématique de Lyell et Coiera (2017) compare la littérature des facteurs humains et celle de la santé : sur 40 études retenues, 6 portent sur la santé. Les auteurs observent le biais d’automatisation (le fait de trop s’en remettre à l’aide à la décision) jusque dans des tâches uniques, lorsque la vérification est complexe, et l’associent à la charge cognitive de la décision. Pour le limiter, ils suggèrent de réduire cette charge.

Mettons ces résultats bout à bout. Quand on automatise les vérifications faciles, l’humain garde les difficiles : les exceptions, les cas ambigus, les erreurs que l’automate n’a pas vues. Ses décisions deviennent moins nombreuses, mais chacune porte sur une erreur plus rare et plus difficile à voir. C’est mon interprétation, pas un résultat mesuré : ce qui reste à vérifier après l’automatisation est ce qui expose le plus au biais.

Ce que l’automatisation retire, ce qu’elle laisse. Interprétation de l’auteur à partir des études citées ; ce schéma ne représente aucune mesure.

L’objection : l’IA allège aussi

L’objection est sérieuse, et elle s’appuie sur des études solides. Erik Brynjolfsson, Danielle Li et Lindsey Raymond ont étudié l’introduction échelonnée d’un assistant conversationnel d’IA générative auprès de 5 172 agents de service client. La productivité, mesurée en demandes résolues par heure, augmente de 15 % en moyenne, avec de fortes différences : les agents les moins expérimentés et les moins qualifiés gagnent en vitesse et en qualité, les plus expérimentés et les plus qualifiés gagnent un peu en vitesse et perdent un peu en qualité. Les auteurs observent aussi une amélioration de l’expérience du travail : les clients sont plus polis et demandent moins souvent à parler à un responsable.

En médecine, une étude publiée en 2025 dans JAMA Network Open interroge 263 cliniciens de six systèmes de santé américains, avant et après trente jours d’usage d’une même plateforme de scribe IA, qui documente les consultations. La part de ceux qui déclarent un burn-out passe de 51,9 % à 38,8 %, et la charge cognitive liée aux notes diminue. L’étude n’a pas de groupe témoin, et trois de ses auteurs déclarent des liens avec l’entreprise Abridge AI, dont deux comme salariés. Elle va néanmoins dans le même sens.

Ces résultats ne contredisent pas les précédents, ils les précisent. La méta-analyse d’Onnasch le montrait déjà : tant que l’automate fonctionne, la charge baisse. Le scribe allège une charge administrative. L’assistant propose des réponses dans une conversation. À mon sens, aucun des deux ne laisse à l’humain une vérification plus difficile que celle qu’il faisait avant. Le coût d’attention n’est donc pas une propriété de l’IA, mais de la répartition du travail : ce qu’on retire à l’humain, et ce qu’on lui laisse.

Ce qu’on ne sait pas encore

Trois questions restent ouvertes, et je préfère les nommer.

Je n’ai trouvé aucune étude qui mesure si l’effort par minute augmente quand l’humain ne garde que le jugement. L’intuition est forte, la preuve manque.

Sur la fatigue des personnes qui supervisent des outils d’IA au travail, la seule source directe que j’ai trouvée est une étude du BCG, menée auprès de 1 488 salariés américains de grandes entreprises et publiée dans la Harvard Business Review. Le BCG conseille lui-même les entreprises sur leur transformation par l’IA. L’étude conclut que l’usage excessif ou la surveillance constante des outils d’IA peut entraîner une fatigue mentale. Une seule étude, publiée par un acteur intéressé : c’est un indice, pas une preuve.

Enfin, je n’ai trouvé aucune étude relue par des pairs qui mesure comment la charge évolue avec le nombre d’agents qu’une même personne supervise.

Trois questions avant d’automatiser

On ne sait pas tout, mais on en sait assez pour poser de meilleures questions à un projet.

1. Que reste-t-il à vérifier une fois l’automatisation en place ? Pas seulement combien de temps : quelles décisions, et sont-elles plus difficiles qu’avant ?
2. Qui le vérifie, et avec quel temps réellement alloué ? Une vérification sans temps alloué devient une signature, pas un contrôle.
3. Que se passe-t-il le jour où l’automate se trompe ? Qui s’en aperçoit, et qui sait encore reprendre la main ?

Le temps gagné reste un bon argument. Il ne dit simplement pas tout ce que l’automatisation coûte, ni à qui.

Pour aller plus loin, le livre blanc gratuit La facture fantôme de l’IA part d’une question voisine, « qui vérifie le résultat ? », et montre comment faire apparaître la relecture, la correction et la supervision dans le coût complet d’un projet.

Les citations d’études publiées en anglais sont traduites par l’auteur.

Sources

• Lisanne Bainbridge, « Ironies of automation », Automatica, 19(6), 1983, p. 775-779 : https://ckrybus.com/static/papers/Bainbridge_1983_Automatica.pdf
• J. S. Warm, R. Parasuraman, G. Matthews, « Vigilance requires hard mental work and is stressful », Human Factors, 50(3), 2008, p. 433-441 : https://pubmed.ncbi.nlm.nih.gov/18689050/
• Hao-Ping (Hank) Lee et al., « The Impact of Generative AI on Critical Thinking: Self-Reported Reductions in Cognitive Effort and Confidence Effects From a Survey of Knowledge Workers », CHI 2025 : https://www.microsoft.com/en-us/research/publication/the-impact-of-generative-ai-on-critical-thinking-self-reported-reductions-in-cognitive-effort-and-confidence-effects-from-a-survey-of-knowledge-workers/
• Berkeley Haas, « AI promised to free up workers’ time. UC Berkeley Haas researchers found the opposite », 18 février 2026, sur les travaux en cours de Xingqi Maggie Ye et Aruna Ranganathan, présentés dans la Harvard Business Review : https://newsroom.haas.berkeley.edu/ai-promised-to-free-up-workers-time-uc-berkeley-haas-researchers-found-the-opposite/
• Joel Becker et al. (METR), « Measuring the Impact of Early-2025 AI on Experienced Open-Source Developer Productivity », 2025 : https://arxiv.org/abs/2507.09089
• Linda Onnasch, Christopher D. Wickens, Huiyang Li, Dietrich Manzey, « Human Performance Consequences of Stages and Levels of Automation: An Integrated Meta-Analysis », Human Factors, 56(3), 2014, p. 476-488 : https://depositonce.tu-berlin.de/items/7df9879f-d77b-4f2b-8760-f0e3d174ca67/full
• Thomas Dratsch et al., « Automation Bias in Mammography: The Impact of Artificial Intelligence BI-RADS Suggestions on Reader Performance », Radiology, 307(4), 2023, e222176 : https://pubmed.ncbi.nlm.nih.gov/37129490/
• D. Lyell, E. Coiera, « Automation bias and verification complexity: a systematic review », JAMIA, 24(2), 2017, p. 423-431 : https://pubmed.ncbi.nlm.nih.gov/27516495/
• Erik Brynjolfsson, Danielle Li, Lindsey Raymond, « Generative AI at Work », Quarterly Journal of Economics, 140(2), 2025, p. 889-942 : https://doi.org/10.1093/qje/qjae044
• Kristine D. Olson et al., « Use of Ambient AI Scribes to Reduce Administrative Burden and Professional Burnout », JAMA Network Open, 8(10), 2025, e2534976 : https://pubmed.ncbi.nlm.nih.gov/41037268/
• BCG, « When Using AI Leads to “Brain Fry” », 5 mars 2026, présentation de l’article de Julie Bedard et al. dans la Harvard Business Review : https://www.bcg.com/news/5march2026-when-using-ai-leads-brain-fry

NOTE: Any updates to this article will be maintained here: https://medley.bryer.org/articles/stratified_glm.html

results <- medley::stratified_glm(
    formula = Private ~ .,
    data = pisa,
    ratio = 1
)

pisa_summary <- summary(results)
pisa_summary |> dplyr::filter(p < 0.1)

                     coefficients   se     p
(Intercept)                 -5.21 2.63 0.051
ST04Q01Male                  0.84 0.28 0.004
ST08Q04TRUE                 -0.50 0.27 0.071
ST08Q05TRUE                  0.75 0.36 0.038
ST20Q06TRUE                  1.79 0.95 0.060
ST20Q07TRUE                  0.72 0.31 0.021
ST21Q05Three or more         0.93 0.45 0.042
ST22Q01201-500 books         1.18 0.59 0.050

plot(results)

pisa_predicted <- predict(results, newdata = pisa)
head(pisa_predicted)

   1    2    3    4    5    6 
0.79 0.16 0.30 0.24 0.22 0.46 

medley::stratified_glm(
    formula = Private ~ .,
    data = pisa,
    ratio = 3
) |> plot()

set.seed(2112); split <- medley::stratified_split(pisa, 'Private', 0.8)
pisa_train <- split$train
pisa_valid <- split$validation
pisa_train$Private |> table() |> prop.table()

FALSE  TRUE 
0.934 0.066 

pisa_valid$Private |> table() |> prop.table()

FALSE  TRUE 
0.934 0.066 

full_model <- glm(Private ~ ., data = pisa_train, family = binomial(link = 'logit'))

full_new_predictions <- predict(full_model, newdata = pisa_valid, type = 'response')

full_roc <- medley::calculate_roc(
    predictions = full_new_predictions,
    observed = pisa_valid$Private
)
summary(full_roc)

AUC = 0.747
Cost of false-positive = 1
Cost of false-negative = 1
Threshold with minimum cost = 0.465 (accuracy = 93.4%)
Threshold with maximum accuracy = 0.465 (accuracy = 93.4%)
Sensitivity and specificity intersection = 0.059 (accuracy = 66.4%)

plot(full_roc)

medley::confusion_matrix(
    observed = pisa_valid$Private,
    predicted = full_new_predictions > attr(full_roc, 'intersection')
)

              predicted             
  observed        FALSE         TRUE
     FALSE 649 (61.99%) 329 (31.42%)
      TRUE   23 (2.20%)   46 (4.39%)
Accuracy: 66.38%
Sensitivity: 66.36%
Specificity: 66.67%

stratified_model <- medley::stratified_glm(
    formula = Private ~ .,
    data = pisa_train,
    ratio = 1
)
stratified_new_predictions <- predict(stratified_model, newdata = pisa_valid)

stratified_roc <- medley::calculate_roc(
    predictions = stratified_new_predictions,
    observed = pisa_valid$Private
)
summary(stratified_roc)

AUC = 0.745
Cost of false-positive = 1
Cost of false-negative = 1
Threshold with minimum cost = 1 (accuracy = 93.4%)
Threshold with maximum accuracy = 1 (accuracy = 93.4%)
Sensitivity and specificity intersection = 0.509 (accuracy = 68.6%)

plot(stratified_roc)

medley::confusion_matrix(
    observed = pisa_valid$Private,
    predicted = full_new_predictions > attr(stratified_roc, 'intersection')
)

              predicted          
  observed        FALSE      TRUE
     FALSE 975 (93.12%) 3 (0.29%)
      TRUE   68 (6.49%) 1 (0.10%)
Accuracy: 93.22%
Sensitivity: 99.69%
Specificity: 1.45%

In October 2020 I wrote about platypus, an R package for object detection and image segmentation: YOLOv3 and U-Net in a few lines of R, on top of the keras and tensorflow packages. People used it. It collected 134 stars and a queue of issues – and then it quietly stopped working, and I stopped answering them.

This year I went back to it. Computer vision is the part of data science I enjoy most, and this package has always been my passion project – the place where I get to learn things my day job does not ask of me. What came out is not an update but a different package, with a new name, a different backend and a narrower purpose: rplatypus, deep learning for medical images in R. It has been submitted to CRAN and is waiting for review; until it is accepted, install it from GitHub:

# install.packages("remotes")
remotes::install_github("maju116/rplatypus")

This post is about why it had to be rewritten rather than fixed, what is in it now, and where the series of posts that follows is going.

pip install "pyplatypus[hub]"

A mask, without training anything

The quickest way to see whether a package does anything useful is to borrow a model someone already trained. rplatypus publishes weights by name – here a U-Net trained on the 2018 Data Science Bowl, the same nuclei dataset I used on this blog in 2020.

The folder nuclei/ holds four microscopy images in the dataset’s own layout – one directory per image, with images/ and masks/ inside.

library(rplatypus)

spec <- platypus_spec(
  data = segmentation_data("nuclei", "nuclei", colormap = binary_colormap),
  models = list(u_net("nuclei", input_shape = c(256, 256), blocks = 4, filters = 16,
                      weights = "dsbowl-unet", fit = FALSE))
)
fit <- platypus_fit(spec, verbose = FALSE)
masks <- predict(fit, split = "validation")
dim(masks)
#> [1]   4 256 256

Four images, 256 by 256 pixels, one class label per pixel. This ran on the processor, with no GPU; platypus_device() tells you whether yours will use one, and why not if it can’t.

images <- list.files("nuclei", pattern = "[.]png$", recursive = TRUE, full.names = TRUE)
images <- read_images(grep("/images/", images, value = TRUE), size = c(256, 256))
plot_masks(images, prediction = masks, colormap = binary_colormap)

The first three are fluorescence images, and the nuclei come out as clean blobs. The fourth is a brightfield image of stained cells, and the model finds specks - the small dark nuclei inside the cells. That is the hard kind of image in this dataset, and the next post looks at it closely.

There is no score here, and that is on purpose. These four images come from the same collection the weights were trained on, so any number measured on them would flatter the model. The numbers that mean something were measured on a held-out split when the weights were published, and available_weights() says what they are and what the model cannot do - for instance, that it is a semantic model, so touching nuclei come back as one region and a count taken from them would be wrong.

nuclei.yaml

task: semantic_segmentation
data:
  train_path: nuclei
  validation_path: nuclei
  colormap: [[0, 0, 0], [255, 255, 255]]
models:
- name: nuclei
  architecture: u_net
  input_shape: [256, 256]
  blocks: 4
  filters: 16
  weights: dsbowl-unet
  fit: false

identical(as.list(platypus_spec("nuclei.yaml")), as.list(spec))
#> [1] TRUE

from pyplatypus import build_engine, from_yaml

engine = build_engine(from_yaml("nuclei.yaml"))
engine.fit()  # with `fit: false` this loads the published weights instead of training
masks = engine.predict("nuclei", split="validation")

Written on behalf of the NEST team.

Over the past few months, we have reached an important milestone in our open-source journey: transitioning key clinical reporting and exploratory analysis R packages from the insightsengineering GitHub organization into the pharmaverse.

The recent migration of the {teal} framework marks another step in this transition, bringing these framework packages under pharmaverse governance.

pak::pkg_install("pharmaverse/teal")
pak::pkg_install("pharmaverse/rtables")

@online{zhu2026,
  author = {Zhu, Joe and Unardi, Dony},
  title = {Roche {Open-Source} {TLG} {Packages} {Move} to Pharmaverse},
  date = {2026-11-03},
  url = {https://pharmaverse.github.io/blog/posts/2026-11-03-roche-open-source-t/roche-open-source-tlg-packages-move-to-pharmaverse.html},
  langid = {en}
}

This is part one of a five-part series on putting AI agents into production with R. Each part builds on the same app, a clinical trial safety assistant. The parts are:

• From Chatbot to Agent (this post)
• Putting the Agent in a Shiny App
• Can You Trust It?
• Audit Trails, Costs and Models
• Sharing Your Tools with MCP

We asked Gemini how many people on the high dose of a drug had itching in a well-known clinical trial. It said 21, and sounded sure of it. We asked again and got 25, with different numbers for the other groups too. Both answers cited the trial’s data files, but the model hadn’t been given them. It was answering from memory, so neither answer could be checked.

In this series we build an assistant whose numbers you can check. The model picks an R function, and R does the counting. We tell the model never to do the maths itself, and part 3 checks that it doesn’t. It can still ask R the wrong question, so every result says what was counted.

No pharma background is needed, since the pattern works for any data your team asks about. Later parts add a Shiny app, tests, limits on who sees patient records, an audit trail, and a way to share the tools with other AI agents.

What is an AI agent?

A chatbot answers from what it learned in training. An AI agent can also ask your code to do things, which is called tool calling (or function calling).

You give the model a list of tools, each with a name, a description and its arguments. When a question needs a tool, the model replies with a request instead of an answer, such as “run count_subjects() with arm = "Placebo"”.

Your code runs the function and sends back the result. This repeats until the model has enough to answer, and that cycle is the agent loop.

The model never runs your code. It can only ask. Your R session decides what runs and what goes back.

{ellmer} runs this loop for you. You write ordinary R functions, describe them and register them with a chat. New to {ellmer}? Our LLMs in R & Python webinar is a good start. As the {ellmer} tool-calling guide recommends, all exact computation happens inside the tools.

What you need

Get this part’s folder from our blog repository, by cloning the repository or downloading it as a ZIP file, and open it as a project in R 4.5. The project uses {renv} to pin every package version, so run this once:

install.packages("renv")
renv::restore()

If you’d rather not use {renv}, install the packages yourself, though you’ll get newer versions than we tested:

install.packages(c("ellmer", "pharmaverseadam", "dplyr", "jsonlite"))

You’ll also need a Gemini API key, and our alt-text post shows how to get one. To store it, run file.edit(".Renviron") in the project, add the line GEMINI_API_KEY=your-key, save the file and restart R. The code’s README has more setup tips, for Windows in particular.

The free tier is fine here, because this trial’s data is public, but not for your own patient or confidential data. Google’s terms also require the paid service for any app offered to users in the UK, the EEA or Switzerland (more on this in part 4).

We tested everything with R 4.5.1, {ellmer} 0.5.0 and {pharmaverseadam} 1.3.0 in September 2026. {ellmer} changes quickly, so if something behaves differently, check your versions first.

The trial data

Our example is CDISC pilot study 01, a worked example published by CDISC, the body that sets clinical trial data standards. It tested a skin patch of an experimental drug, xanomeline, in people with mild to moderate Alzheimer’s disease. Participants, called subjects in trial data, wore a placebo, low-dose or high-dose patch for 26 weeks.

The {pharmaverseadam} package rebuilds the study’s analysis datasets in CDISC’s ADaM format, and we need two of them: adsl, with one row per subject, and adae, with one row per adverse event. An adverse event is anything medically unwanted that happened to a subject during the trial, whether or not the drug caused it.

Our tools count the way a safety table usually does, and two of their choices are easy to miss. Both change every count:

• Not everyone in adsl joined the trial. Of its 306 people, 52 are screen failures, so the safety population (SAFFL == "Y") is 254 subjects.
• Planned and received treatment can differ. Each subject has a planned arm (ARM) and the treatment they received (ACTARM). Twelve subjects randomised to the high dose stopped within two weeks, before their dose went up to 81 mg, so here they count as low dose.

If you work with trial data, the note for pharma readers near the end covers a few choices a formal analysis plan might make differently.

The tools also keep only the treatment-emergent events (TRTEMFL == "Y"), which leaves 1,122 of the 1,191 event rows. One function, in R/data.R, does all of this:

ARMS <- c("Placebo", "Xanomeline Low Dose", "Xanomeline High Dose")

load_trial_data <- function() {
 subjects <- pharmaverseadam::adsl |>
 dplyr::filter(SAFFL == "Y") |>
 dplyr::select(USUBJID, ACTARM, AGE, SEX, RACE)
 events <- pharmaverseadam::adae |>
 dplyr::filter(TRTEMFL %in% "Y", USUBJID %in% subjects$USUBJID) |>
 dplyr::select(USUBJID, ACTARM, AEDECOD, AETERM, AEBODSYS, AESER, AESEV, ASTDT)
 list(
 subjects = subjects,
 events = events,
 # Every term and body system recorded in the trial, including the few
 # never treatment-emergent, so asking about them gives a real zero.
 terms = sort(unique(pharmaverseadam::adae$AEDECOD)),
 body_systems = sort(unique(pharmaverseadam::adae$AEBODSYS))
 )
}

Asking without tools

We started by asking gemini-3.7-flash the opening question twice, with no tools:

# Without tools, the model can only guess.
bare <- ellmer::chat_google_gemini()
bare$chat("In CDISC pilot study 01, how many subjects on high-dose xanomeline had pruritus?", echo = "output")

Pruritus is itching. Both answers are pasted as they came back, Markdown and all, except that their headings are shown in bold.

Neither answer matches the data we hold, but each high-dose figure is a correct number for something else:

On placebo, the data says 8 of 86 (9.3%), and the other counts in both answers match nothing.

The second answer is the more dangerous one. Its count is right, and only the arm size differs from ours, so it survives a quick check. And nothing in either answer says which data, which population or which events it counted.

Writing the tools

The model never sees a tool’s code, only its description, so the description is what tells it when to use the tool. Most of the work happens in count_subjects(), in R/tools.R:

count_subjects <- function(data, arm = NULL, serious_only = FALSE, term = NULL, body_system = NULL) {
 arm <- check_arm(arm)
 serious_only <- check_flag(serious_only, "serious_only")
 term <- check_known(
 check_string(term, "term"), data$terms, "preferred term",
 "Call ae_summary to see the terms used in this trial."
 )
 body_system <- check_known(
 check_string(body_system, "body_system"), data$body_systems, "body system",
 "Call ae_summary with by = 'body_system' to see them."
 )
 data <- filter_arm(data, arm)
 events <- data$events
 if (serious_only) events <- dplyr::filter(events, AESER == "Y")
 if (!is.null(term)) events <- dplyr::filter(events, AEDECOD == term)
 if (!is.null(body_system)) events <- dplyr::filter(events, AEBODSYS == body_system)
 n <- dplyr::n_distinct(events$USUBJID)
 total <- nrow(data$subjects)
 jsonlite::toJSON(list(
 arm = arm %||% "All arms",
 # Say what was counted, so a wrong reading of the question shows.
 filters = list(
 serious_only = serious_only,
 term = term %||% "any",
 body_system = body_system %||% "any"
 ),
 subjects_with_event = n,
 subjects_in_arm = total,
 percent = round(100 * n / total, 1)
 ), auto_unbox = TRUE)
}

Four choices matter most:

• It returns JSON. The model reads text, and {ellmer} 0.5.0 deprecates returning data frames from tools.
• It says what it counted. The result repeats the filters it applied, so if the model misread the question, you can see it.
• Unknown input is an error, never a zero. Ask for a term the trial doesn’t have and the tool fails, with a message telling the model how to recover. A zero that means “I didn’t understand” is worse than no answer.
• Every argument is checked. Models don’t always send the types they’re asked for, and in part 5 other programs call these tools over MCP, a standard way for AI apps to share tools.

We added that last check after a problem in testing. Passing serious_only = "true", as text rather than TRUE, made the tool quietly count every subject with any adverse event. It answered 217 instead of 3.

Because it’s an ordinary R function, you can call it without a model. Source the files in R/ first:

for (file in list.files("R", pattern = "[.]R$", full.names = TRUE)) source(file)
count_subjects(load_trial_data(), arm = "Xanomeline High Dose", term = "PRURITUS")
#> {"arm":"Xanomeline High Dose","filters":{"serious_only":false,"term":"PRURITUS","body_system":"any"},"subjects_with_event":25,"subjects_in_arm":72,"percent":34.7}

The tool definition is what the model sees:

ellmer::tool(
 function(arm = NULL, serious_only = FALSE, term = NULL, body_system = NULL) {
 count_subjects(data, arm, serious_only, term, body_system)
 },
 name = "count_subjects",
 description = paste(
 "Count subjects in the safety population with at least one treatment-emergent",
 "adverse event matching the filters. Returns the count, the number of subjects",
 "in the arm and the percentage."
 ),
 arguments = list(
 arm = arm_type,
 serious_only = ellmer::type_boolean("Only count serious adverse events.", required = FALSE),
 term = ellmer::type_string(
 "MedDRA preferred term, e.g. 'PRURITUS'. Leave out to count any event.",
 required = FALSE
 ),
 body_system = ellmer::type_string(
 "MedDRA body system, e.g. 'NERVOUS SYSTEM DISORDERS'. Leave out for any.",
 required = FALSE
 )
 ),
 annotations = ellmer::tool_annotations(title = "Counting subjects", read_only_hint = TRUE)
)

The annotations give the tool a title, which part 2’s app uses to label each tool call. The arm is an enum, a fixed list of allowed values, so the model sees the three valid arms:

arm_type <- ellmer::type_enum(ARMS, "Treatment arm. Leave out to cover all arms.", required = FALSE)

check_arm() still rejects anything else, for callers that ignore the list. In all, trial_tools() builds four tools around one copy of the data: list_arms() for arm sizes, ae_summary() for the most common events, count_subjects(), and subject_profile() for one patient’s record.

Running the loop

One small function in R/agent.R puts the agent together, creating a Gemini chat with the system prompt (next section) and registering the tools:

make_agent <- function(tools = trial_tools(), prompt_path = "prompts/system.md",
 model = NULL, chat_fun = ellmer::chat_google_gemini) {
 chat <- chat_fun(system_prompt = read_prompt(prompt_path), model = model)
 chat$register_tools(tools)
 chat
}

With model = NULL you get {ellmer} 0.5.0’s default, gemini-3.7-flash. That default changes between releases, so part 4 pins the model by name.

Switching provider is usually one argument, such as make_agent(chat_fun = ellmer::chat_anthropic), though local models through chat_ollama() have a catch that part 4 covers.

We asked two questions with echo = "all", which prints the question, every tool request and every result, as well as the answer. Below are the tool requests and some of the results. We’ve kept the high dose’s ae_summary() result in full, so you can match it to the answer.

chat <- make_agent()
chat$chat("Which adverse events were most common on the high dose, compared with placebo?", echo = "all")
chat$chat("How many subjects had a serious adverse event, and on which arms?", echo = "all")

< [tool request (ae_summary)]: ae_summary(arm = "Xanomeline High Dose", by =
< "term", top_n = 10L)
< [tool request (ae_summary)]: ae_summary(arm = "Placebo", by = "term", top_n =
< 10L)
[tool result (ae_summary)]: {"arm":"Xanomeline High Dose","serious_only":false,"subjects_in_arm":72,"by":"term","rows":[{"event":"PRURITUS","subjects":25,"percent":34.7},{"event":"APPLICATION SITE PRURITUS","subjects":21,"percent":29.2},{"event":"APPLICATION SITE ERYTHEMA","subjects":14,"percent":19.4},{"event":"ERYTHEMA","subjects":14,"percent":19.4},{"event":"DIZZINESS","subjects":10,"percent":13.9},{"event":"APPLICATION SITE IRRITATION","subjects":9,"percent":12.5},{"event":"HYPERHIDROSIS","subjects":8,"percent":11.1},{"event":"RASH","subjects":8,"percent":11.1},{"event":"SINUS BRADYCARDIA","subjects":8,"percent":11.1},{"event":"APPLICATION SITE DERMATITIS","subjects":7,"percent":9.7}]}
< [tool request (count_subjects)]: count_subjects(arm = "Placebo", term =
< "DIZZINESS")
< [tool request (count_subjects)]: count_subjects(arm = "Placebo", term =
< "APPLICATION SITE IRRITATION")
< [tool request (count_subjects)]: count_subjects(arm = "Placebo", term =
< "HYPERHIDROSIS")
< [tool request (count_subjects)]: count_subjects(arm = "Placebo", term =
< "SINUS BRADYCARDIA")

< [tool request (count_subjects)]: count_subjects(serious_only = TRUE)
< [tool request (count_subjects)]: count_subjects(arm = "Placebo", serious_only
< = TRUE)
< [tool request (count_subjects)]: count_subjects(arm = "Xanomeline Low Dose",
< serious_only = TRUE)
< [tool request (count_subjects)]: count_subjects(arm = "Xanomeline High Dose",
< serious_only = TRUE)
[tool result (count_subjects)]: {"arm":"All arms","filters":{"serious_only":true,"term":"any","body_system":"any"},"subjects_with_event":3,"subjects_in_arm":254,"percent":1.2}
[tool result (count_subjects)]: {"arm":"Placebo","filters":{"serious_only":true,"term":"any","body_system":"any"},"subjects_with_event":0,"subjects_in_arm":86,"percent":0}
[tool result (count_subjects)]: {"arm":"Xanomeline Low Dose","filters":{"serious_only":true,"term":"any","body_system":"any"},"subjects_with_event":2,"subjects_in_arm":96,"percent":2.1}
[tool result (count_subjects)]: {"arm":"Xanomeline High Dose","filters":{"serious_only":true,"term":"any","body_system":"any"},"subjects_with_event":1,"subjects_in_arm":72,"percent":1.4}

We never told the model which tools to use. For the first question, it called ae_summary() once per arm for the top ten events. Four of the high dose’s ten weren’t in placebo’s, so it counted those on placebo with count_subjects() before it answered. For the second, it called count_subjects() with serious_only = TRUE for all arms together, then for each arm. Every number in both answers appears in a tool result.

To run all of this yourself, source("demo/part-1.R") asks the no-tools question and both agent questions, which takes about six requests. The offline tests run with source("tests/run.R") and make no API calls.

The system prompt

The system prompt sets the rules for every conversation:

You answer questions about the safety data of one clinical trial: CDISC pilot study 01, which compared a xanomeline skin patch at two doses with placebo.

- Use the tools for every fact about the trial. Every number in your answer must appear in a tool result from this conversation. Never calculate, estimate or remember numbers yourself.
- Quote numbers exactly as the tools give them.
- Say which arm each number is for. Counts are subjects in the safety population, by the treatment they actually received, with treatment-emergent adverse events.
- If the tools cannot answer a question, say so.
- Keep answers short.

Two lines do most of the work. The first, that every number must come from a tool, is the rule the rest of the series tests. The third makes the model name the arm and population, so a reader always has a denominator to check.

A prompt is still only a request, and a model can ignore it. That’s why part 3 adds a check on where each number came from.

Why narrow tools

You could give the model a single run_r_code() tool and let it answer anything. But then you’d be validating whatever code it happened to write that day.

A narrow tool can be tested against numbers checked by hand, logged, and allowed or refused on its own. It’s easier for the model too, since picking the right tool is a smaller job than writing correct {dplyr} code.

What could still go wrong

The agent works, but it isn’t ready for a team yet:

• The model invents or rounds a number anyway. The prompt only asks it not to. Part 3 checks every answer’s numbers against the tool results.
• Patient-level data. Right now the model can open any patient’s record. Part 3 limits subject_profile() to a reviewer group, whose members confirm each record before it opens, and part 5 leaves it out of the tools shared with other AI agents.
• Many users at once. A console session suits one person. In part 2’s Shiny app, each person gets their own conversation.
• Cost and outages. Within minutes, the free tier sent us “503 high demand” and “429 quota exceeded” errors. It’s small, too: when we wrote this, about 20 requests a day per model, and a question takes one request per round of tool calls. {ellmer} tries each request up to three times by default, and options(ellmer_max_tries = …) changes that. Parts 2 and 4 deal with the rest.

Summary

An agent is a model, some tools and a loop. The model decides what to count, R does the counting on carefully filtered data, and the prompt asks the model to quote only what the tools return.

Without tools, it gave us two confident answers that disagreed, and neither matched our data. With tools, every number came from R.

Next time, we put the agent in a Shiny app, so that several people can use it at once.

All the code for this part, including demo/part-1.R and the tests, is in our blog repository.

Building something like this? Tell us about it at AI in Production 2027, where the call for talks is open until 29 January 2027. For a guided tour of LLMs from R and Python, see our LLM-Driven Applications with R and Python course.

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

We started this article as an interview series with Champions from our program, and we decided to publish for our 15th anniversary series. This is the interview with Carolina Pradier, main maintainer of the eph package.

Carolina Pradier has long worked to make the world a more equitable place. Her doctoral research focused on inequalities in science with a focus on the mechanisms that structure gender inequality in the Latin American scientific community. ¹ She is also part of Ecofeminita (Asociación Civil Economía Feminista), an interdisciplinary organization created and led by women. Carolina has a degree in Economics from the University of Buenos Aires and a holds a Master’s degree in Social Sciences of Work from the Faculty of Social Sciences of the University of Buenos Aires. She is a PhD candidate in Information Science from the University of Montreal and teaches at post-secondary institutions.

As a researcher, Carolina was part of the Center for Studies on Population, Employment and Development and now works for the UNESCO Chair on Open Science. She also is the maintainer of the eph R package, whose objective is to facilitate working with the Permanent Household Survey data.

rOpenSci

Carolina’s encounter with the rOpenSci Champions Program seemed natural. Her goals and those of the organization were aligned. In her application to the Champions Program she expressed:

I would like to be part of the rOpenSci Champions Program to improve my skills as a developer, as well as to improve my abilities to teach and communicate these skills.”

Next followed over 12 months of training. Through the coaching, mentoring, and exposure to colleagues, Carolina received information on tools and processes, as well as advice and feedback to improve her package. She was also able to develop her communication skills and feel more confident in communicating more effectively with the people who use eph.

Exchanges and meeting with other people made her work easier, and helped her stay organized and anticipate future problems she might encounter:

I think I felt great when we finished the training and I started to work closely with my mentor, we had a good relationship and we could work very well and efficiently.”

In addition to being selected to participate in the rOpenSci Champions Program through the rOpenSci application, Carolina had already gone through a previous selection. The rOpenSci Champions Program accepts only individual applications. Therefore, when Carolina’s team wished to participate, her colleagues encouraged her to participate on behalf of the team.

During the program, Carolina developed different events to help disseminate her work, share it with others, and invite them to get involved in open science. This idea was also present in the objectives that Carolina expressed in her first contact with rOpenSci:

One of my main motivations is to bring new tools to the communities of social scientists in which I actively participate in order to strengthen our capacities to transform reality.”

It was with this goal that Carolina organized an event in Spanish with news and ideas for the use of the eph package. This achieved a double objective: sharing content about R packages in Spanish and bringing more people closer to the R and open science community.

Carolina also participated in a Community Call on the rOpenSci Champions Program where she shared her experiences as a champion, co-hosted a co-working session about Building your first R package, and joined as a mentor for the cohort 2025-2026 of the Champions Program.

Perfecting a job well done

The eph package processes data from the Permanent Household Survey (INDEC) in Argentina to replicate or conduct independent analyses from those published in the reports prepared by the government agency. It also allows the analysis of different categories and indicators not provided in the official reports. Carolina explained that throughout the program they managed to improve the package:

“We tried to improve the quality of the code (reliability and efficiency), we improved the documentation and created new functions.

They also worked hard to translate the documentation and make it more accessible to Spanish-speaking users. The rOpenSci peer review process showed how a rigorous approach can also be welcoming and kind.

Open Science

The eph package is a tool for academic research, journalism and evidence-based public policy formulation. The Champions Program allowed Carolina to improve the software using different points of view from peer and expert review.

Beyond improving the package itself, one of Carolina’s goals was to help increase the impact of her community’s work:

I am particularly excited about the prospect of participating in a program aimed at making the R community more inclusive. in addition, I feel that this experience would be a great opportunity as it would allow me to meet and work with highly skilled professionals whose work I deeply admire. Not only am I interested in learning from their experience, but I also believe that building these networks can help increase the impact of my community’s efforts and access new sources of support.”

Some existing uses of eph show the kind of impact Carolina hopes the package can have. At EcoFeminita, the civil association Carolina is part of, the Ecofemidata team uses eph to provide data on gender inequality through indicators, statistics, and reproducible thematic reports. Among other uses, the team maintains an app that analyzes the labor market, income, time use, and domestic work.

The package has also inspired programmers, data scientists, and people from other disciplines to develop the PyEPH library, which has a similar purpose to the eph package in R but is used and supported by people in the Python data community. As eph keeps improving and the network around it grows, Carolina hopes to see more experiences like these.

Beyond learning: building with a community

Carolina Pradier was able to improve her R skills, acquire new skills, and offer a multilingual package so that more people can access open science. Experiences with her mentor, her peers, the people who participated in the review process, and the people who use the package were a great contribution to learning and personal growth. Her participation in the Champions Program was the beginning of a network of links that promise new projects.

I always had someone to talk to to help me solve problems; it was never a frustrating process.”

Carolina is now handing over maintenance of eph to a new maintainer. The package is used by many people and has a community around it, which is what makes it sustainable.

Shared knowledge allows us to anticipate, help and include. The more experiences like Carolina Pradier’s there are, the more shared and open open science there will be.

Clinical medicine is by far the largest field in the 2026 Stanford–Elsevier top 2% list: 76,155 of the 249,310 listed authors, 31%, have it as their main field, more than the next three of the 20 fields combined. I wanted to know where those scientists work, and whether the country ranking still looks the same once you ask how many people each country has.

The data is John Ioannidis’s science-wide author database, built with Elsevier from Scopus (Ioannidis, J.P.A. (2026), “August 2026 data-update for ‘Updated science-wide author databases of standardized citation indicators’”, Mendeley Data, v9, doi:10.17632/btchxktzyw.9, licensed CC BY-NC 3.0). An author is listed if they rank in the top 100,000 by a composite citation score or in the top 2% of their main subfield. I use the single-year table, which scores citations received during 2025. Each row carries a main field from the Science-Metrix classification (sm-field) and the country of the author’s most recent institution (cntry).

Find a top-cited scientist or institution: If you would rather skip the tutorial and find your own name or institution, topcitedscientists.com has the whole 2026 list: search for a scientist by name or see institutions ranked by their number of top-cited scientists.

How do I count the Stanford–Elsevier top 2% scientists by country in R?

Download the single-year author table from Mendeley Data, read its Data sheet with readxl, fill the blank countries from the top_cntry_career column, then filter on sm-field and count() the countries. This block runs on its own with readxl 1.5.0 and tidyverse 2.0.0 on R 4.6.1. The file is about 100 MB, so it raises the download timeout first:

library(tidyverse)
library(readxl)

url <- "https://data.mendeley.com/public-files/datasets/btchxktzyw/files/a517a94c-3fb0-4cee-b97f-f2d476df66f0/file_downloaded"
file <- "Table_1_Authors_singleyr_2025_pubs_since_1788_wopp_extracted_202608.xlsx"

options(timeout = 600)
if (!file.exists(file)) {
  download.file(url, file, mode = "wb")
}

authors <- read_excel(file, sheet = "Data") |>
  mutate(
    filled = is.na(cntry),
    cntry = coalesce(cntry, str_to_lower(top_cntry_career))
  )

clinical <- authors |>
  filter(`sm-field` == "Clinical Medicine") |>
  count(cntry, sort = TRUE) |>
  mutate(share = n / sum(n))

clinical
## # A tibble: 147 × 3
##    cntry     n  share
##    <chr> <int>  <dbl>
##  1 usa   28458 0.374 
##  2 gbr    6303 0.0828
##  3 chn    3924 0.0515
##  4 deu    3770 0.0495
##  5 ita    3710 0.0487
##  6 can    3140 0.0412
##  7 aus    2943 0.0386
##  8 fra    2120 0.0278
##  9 jpn    2074 0.0272
## 10 nld    1955 0.0257
## # ℹ 137 more rows

Three details in that code matter. download.file() gives up after 60 seconds by default, which a slow connection hits long before 100 MB. The workbook’s first sheet is a column key, so read_excel(file) without sheet = quietly returns its 53 rows of definitions instead of the data. And 169 authors have no cntry, 77 of them in clinical medicine, the most of any field; the dataset’s top_cntry_career column (the country that appears most often on their papers) fills them, and the filled column keeps a record of which rows were changed. That moves the US from 28,424 to 28,458.

The United States has 28,458 of the field’s top-cited scientists, 37.4%. The United Kingdom follows with 6,303, then China (3,924), Germany (3,770) and Italy (3,710). Those five hold 61% of the field, and 146 countries have at least one listed author.

Related posts on DataScience+:

• Machine Learning Results in R: one plot to rule them all! (Part 1 – Classification Models)
• Visualize your Portfolio’s Performance and Generate a Nice Report with R
• Creating Slopegraphs with R

Per resident, Switzerland and Denmark lead

A raw count mostly measures how big a country’s research system is. Dividing by population asks a different question: how dense the country’s pool of highly cited clinicians and medical researchers is. I take 2024 population from the World Bank’s API (indicator SP.POP.TOTL, keyless JSON) and keep countries with at least a million residents:

wb <- jsonlite::fromJSON(paste0(
  "https://api.worldbank.org/v2/country/all/indicator/",
  "SP.POP.TOTL?date=2024&format=json&per_page=400"
))

population <- wb[[2]] |>
  as_tibble() |>
  transmute(
    cntry = str_to_lower(countryiso3code),
    country = country$value,
    pop = value
  ) |>
  filter(cntry != "")

per_million <- clinical |>
  inner_join(population, by = "cntry") |>
  filter(pop >= 1e6) |>
  mutate(per_million = n / pop * 1e6) |>
  arrange(desc(per_million))

per_million |> select(country, n, per_million) |> head(8)
## # A tibble: 8 × 3
##   country            n per_million
##   <chr>          <int>       <dbl>
## 1 Switzerland     1345       149. 
## 2 Denmark          820       137. 
## 3 Netherlands     1955       109. 
## 4 Australia       2943       108. 
## 5 Sweden          1028        97.3
## 6 United Kingdom  6303        91.0
## 7 United States  28458        83.7
## 8 Canada          3140        76.1

Switzerland has 149 top-cited clinical medicine scientists per million residents and Denmark 137, against 84 in the United States, which ranks 7th per resident behind the Netherlands, Australia, Sweden and the United Kingdom. China is 2.8 per million and India 0.9. The floor matters: without it, Monaco’s three listed authors would rank it among the leaders. Taiwan (428 authors) drops out of this table because the World Bank does not publish its population.

A country’s count also reflects its field mix

China is second in the whole list, with 37,204 listed authors across all fields, yet only third in clinical medicine. The reason is what its top-cited scientists work on:

field_mix <- authors |>
  summarise(
    listed = n(),
    clinical = mean(`sm-field` == "Clinical Medicine"),
    .by = cntry
  ) |>
  arrange(desc(listed))

field_mix |> head(5)
## # A tibble: 5 × 3
##   cntry listed clinical
##   <chr>  <int>    <dbl>
## 1 usa    73722    0.386
## 2 chn    37204    0.105
## 3 gbr    18320    0.344
## 4 deu    10942    0.345
## 5 aus     8974    0.328

Clinical medicine is the main field of 39% of listed US authors and 34% of British ones, against 31% across the whole list and 11% in China, where the largest field is Enabling & Strategic Technologies (8,807 authors). So a field-by-country ranking says as much about where a country’s research effort goes as about how strong it is in that field.

Two caveats travel with every number here. The list counts citations, not quality, and the dataset’s authors warn against reading it as absolute evidence of excellence. And cntry is where an author works now, not their nationality or where they trained.

This article was first published on DataScience+, a community of R and Python tutorial authors. Have a data-science technique worth sharing? Write for us — no pitch required.

Learn the R skills needed to analyse data confidently, build powerful statistical models, create publication-quality figures, and answer complex research questions with the Centre for Biodiversity and Conservation Science (CBCS) 2027 Summer R Workshops @ UQ.

These hands-on workshops are designed for both beginners and researchers looking to strengthen their analytical skills. You’ll learn how to wrangle data efficiently, develop reproducible workflows, apply appropriate statistical methods, and avoid common analytical pitfalls. By the end of the workshop, you’ll have practical skills that can be immediately applied to your own research projects. Although most of our example datasets are ecological in nature, the approaches and skills are equally applicable in many other research contexts.

MICE always felt like magic to me So I tried rebuilding it from scratch in R to understand it better. Lots of mistakes and help along the way, but my version ended up close to the mice package. Might be just luck. Still learning, but here are my messy notes!

Motivations

We used mulitple imputation by chained equations (MICE) quite often in imputing missing data in many occasion, but I never really understood the mechanism behind it except that it appears magical! I want to see if I can build this from scratch as a mean to understand the foundation of MICE so that at least it will give me the illusion of understanding. Join me in my messy journey of trying to code it from the ground up (FYI, I’m not coding OLS or MLE from the ground up…), and let’s see if we can create some “magic” here! This is going to be a bumpy one with lots of coding, and maybe less of mathematical formulae. I’ll try my best to provide a TL;DR for each section so that we can get a high level view as well. Buckle up! And oh, Happy Halloween! What a better time to look at magical recipes of statistics! Let’s go!

Objectives:

• What Is MICE
• Simulate Our DGP
• The Procedure To Impute And Create Multiple Datasets
  • Pool them!
• Opportunities For Improvement
• Lessons Learnt

What Is MICE

TL;DR: fill in each gap several times with plausible random guesses, analyse each completed dataset, then combine the results.

When data are missing, the tempting fix is to plug in one best guess, like the mean. The problem is that a single guess pretends we know the missing value exactly, so our standard errors come out too small. MICE takes a different approach: it fills in each missing value with a plausible random draw, and does this several times to create multiple complete datasets. The differences between those datasets capture our uncertainty about what the missing values really were.

The “chained equations” part is how those draws are made. Instead of one giant model for all variables at once, MICE fits a simple regression for each incomplete variable, using all the others as predictors, and cycles through them, each one feeding its freshly imputed values into the next. After a few rounds, the imputations settle down. Finally, we analyse each completed dataset and pool the results with Rubin’s rules.

Simulate Our DGP

Alright let’s simulate our data generating process (the truth) and add in missing data mechanism to create our dataset that we want to impute from (imagine this is the dataset that we have and we want to apply mice on it).

library(tidyverse)
library(broom)

set.seed(1)
n <- 1000
p1 <- rnorm(n)
p2 <- rnorm(n, 0.5*p1)
p3 <- rnorm(n, 0.4*p1 + 0.3*p2)
p4 <- rbinom(n, 1, 0.5)
p5 <- rbinom(n, 1, plogis(-0.5 + 0.3*p1 - 0.2*p2 + 0.1*p3))
outcome <- rbinom(n, 1, plogis(-0.5 + 0.3*p1 - 0.2*p2 + 0.1*p3 + 0.5*p4 + 0.4*p5))

df <- tibble(p1,p2,p3,p4,p5,outcome)

# missing mechanism 
## mcar
p1_miss <- ifelse(runif(n) < 0.1, NA, p1)
p3_miss <- ifelse(runif(n) < 0.2, NA, p3)

## mar 
p2_miss <- ifelse(runif(n) < plogis(-0.5 + 0.3*p1), NA, p2)
p4_miss <- ifelse(runif(n) < plogis(-0.5 + 0.3*p1 - 0.2*p2), NA, p4)

## mnar
p5_miss <- ifelse(runif(n) < plogis(-0.5 + 0.8*p5), NA, p5)

## recombine
df_miss <- tibble(p1 = p1_miss, p2 = p2_miss, p3 = p3_miss, p4 = p4_miss, p5 = p5_miss, outcome = outcome)

The Procedure To Impute And Create Multiple Datasets

Alright, let’s go through the procedure for MICE real quick.

1. Index all missing values and record them in a dataframe.
2. For each chain:

• Fill every missing value with a random starting value.
• For each column with missing data:
  • Bootstrap the rows where that column is observed, using the current filled-in values for the other columns.
  • Fit a model with that column as the response and all other columns, including the outcome, as predictors: linear regression for continuous columns, logistic for binary.
  • Predict the missing values, then add noise. For continuous columns, add a normal draw using the model’s residual SD. For binary columns, draw 0 or 1 from the predicted probability.
  • Write the new values back immediately, so the next column’s model uses them. This is the “chained” part.
• Repeat the column loop for 10 iterations.
• Keep the final imputed values from the 10th iteration.

3. Repeat step 2 for 10 chains.
4. Build 10 completed datasets, one per chain.
5. Fit the analysis model on each dataset.
6. Pool by averaging the coefficients, plus Rubin’s rules if you need standard errors.

Phew. What a procedure! At least we now know the machinery behind MICE!

Index all missing values and record them in a dataframe

## procedure of simple mice
### 1. index all missing values
df_miss_idx <- df_miss |> is.na()

#### create an empty dataframe based on number of missing data
empty_row <- df_miss_idx |> as.vector() |> sum()

empty_df <- which(df_miss_idx, arr.ind = T) |> as_tibble()

data_type <- df_miss |> 
  map_chr(~class(.x)[1]) 

data_type_colnum <- 1:length(data_type)
colname <- colnames(df_miss)

data_type_df <- tibble(colnum = data_type_colnum, dtype=data_type, colname = colname)

empty_df <- empty_df |>
  left_join(data_type_df, by = c("col" = "colnum")) |>
  mutate(value = NA)

Notice that we added datatype to help us identify continuous data vs binary data?

Fill in missing values with random starting values

Let’s create a an empty dataframe then write functions to do these so we can use them on our for loop

### 2. First iteration, first chain
### 2a. randomly impute all values, use conditions per column to assess types of data before random normal, binomial distribution
fill_empty_df <- function(empty_df) { 
  df <- empty_df |>
  mutate(value = case_when(
    dtype == "numeric" ~ rnorm(n()),
    dtype == "integer" ~ rbinom(n(),1,0.5),
    TRUE ~ NA_real_
  )) |>
  mutate(model = case_when(
    dtype == "numeric" ~ "linreg",
    dtype == "integer" ~ "logreg",
    TRUE ~ NA_character_
  ))
  
  return(df)
}

grab_xvec <- function(df, outcome) {
  df |> 
    select(-all_of(outcome)) |>
    colnames()
}

gen_formula <- function(df, outcome) {
  x_vec <- grab_xvec(df, outcome)
  x <- paste(x_vec,collapse = "+")
  formula <- paste0(outcome,"~",x)
  return(formula)
}


formula_df <- empty_df |>
  distinct(colname, .keep_all = T) |>
  select(col, colname) |>
  mutate(formula = map(.x=colname, .f=~gen_formula(df_miss,.x)) |> unlist()) 


missing_fill <- function(empty_df) {
  df_miss_idx_fill <- matrix(nrow = nrow(df_miss), ncol = ncol(df_miss)) 
  
  for (i_mis in 1:empty_row) {
    rownum <- empty_df[[i_mis, "row"]]
    colnum <- empty_df[[i_mis, "col"]]
    data <- empty_df[[i_mis, "value"]]
    df_miss_idx_fill[rownum,colnum] <- data
  }
  
  df_miss_idx_fill <- df_miss_idx_fill |> as_tibble() 
  colnames(df_miss_idx_fill) <- colnames(df_miss)
  df_miss_idx_fill <- df_miss_idx_fill |>
    mutate(p4 = as.integer(p4),
           p5 = as.integer(p5))
  
  df_miss_df_miss_idx_fill <<- map2_dfc(df_miss, df_miss_idx_fill, coalesce) #learn about map2_dfc
}

empty_df contains randomly pulled number from respective datatype. grab_xvec grabs the predictor columns then gen_formula turns that into model formula for fitting. And missing_fill starts filling out the missing data Sorry if my code is not efficient. I’m sure we can have LLM make it even more efficient

Next, we iterate!

# grab the colomuns that have missing values 
col_missing <- empty_df |> distinct(col) |> pull()

# create empty df so we can fill in
empty_df_mice <- tibble()
  
# set up chain
for (chain in 1:10) {
  
  # internally create another df so we can use and override data
  impute_df <- fill_empty_df(empty_df)
  
# for each chain iterate up to 10 times (this is where the chain equation starts)
for (iter in 1:10) {
  
# for each iteration look at each column with missing data
for (i in col_missing) {
  
  # refresh the filled dataset with the latest imputations
  missing_fill(impute_df)
  
  # set seed for each chain, iter, and column so that we can reproduce the same for testing purposes
  seed <- paste0(chain,iter,i) |> as.numeric()
  set.seed(seed)
  
  # look for datatype of the column we're evaluating now
  dtype <- impute_df |>
    filter(col == i) |>
    distinct(model)
  
  # a long winded way of sorting what columns to use for modeling and what to readd back
  empty_df_i <- impute_df |>
    filter(col == i)
  
  empty_df_i_not <- impute_df |>
    filter(col != i)
  
  # finding out what rows are missing and remove those for modelling
  missing_vec <- empty_df_i |>
    pull(row)
  
  df_miss_i_full <- df_miss_df_miss_idx_fill |>
    slice(-missing_vec)
  
  # resample (in place of posterior sampling)
  df_miss_i <- slice_sample(df_miss_i_full, replace = T, n = nrow(df_miss_i_full))
  
  
  # find out the formula for model
  formula <- formula_df |>
    filter(col == i) |>
    pull(formula)
  
  # algo for which estimator to use
  if (dtype == "linreg") {
  model_i <- lm(formula, data = df_miss_i)
  } 
  
  if (dtype == "logreg") {
   model_i <- glm(formula, data = df_miss_i, family = "binomial")   
  }

  newdata <- df_miss_df_miss_idx_fill |>
    slice(missing_vec)
  
  if (dtype == "logreg") {
  empty_df_i <- empty_df_i |>
    mutate(value = rbinom(nrow(empty_df_i), 1, predict(model_i, newdata = newdata, type = "response")))
  } else {
  empty_df_i <- empty_df_i |>
    mutate(value = predict(model_i, newdata = newdata) + rnorm(nrow(empty_df_i),0,sigma(model_i)))
  }
  
  impute_df <- empty_df_i_not |>
    rbind(empty_df_i)
  
}
  # only save the last iter data for each chain
  if (iter == 10) {
    empty_df_mice <- empty_df_mice |>
      rbind(impute_df |> mutate(chain = chain)) 
  }
  
}
}

# write a function to create dataframe based on chains
missing_fill_full <- function(x) {
  df_miss_idx_fill <- matrix(nrow = nrow(df_miss), ncol = ncol(df_miss)) 
  
  for (i_mis in 1:empty_row) {
    rownum <- x[[i_mis, "row"]]
    colnum <- x[[i_mis, "col"]]
    data <- x[[i_mis, "value"]]
    df_miss_idx_fill[rownum,colnum] <- data
  }
  
  df_miss_idx_fill <- df_miss_idx_fill |> as_tibble() 
  colnames(df_miss_idx_fill) <- colnames(df_miss)
  df_miss_idx_fill <- df_miss_idx_fill |>
    mutate(p4 = as.integer(p4),
           p5 = as.integer(p5))
  
  # merges both df and collapse those overlap completely
  map2_dfc(df_miss, df_miss_idx_fill, coalesce) 

}

# create 10 dataframes for each chain
for (data in 1:10) {
  df_miss_mice_i <- missing_fill_full(empty_df_mice |> filter(chain == data))
  assign(paste0("df_miss_mice_",data), df_miss_mice_i)
}

Wow, OK that’s a lot of code! Hopefully the comments above will be self-explainatory. Some of my code may be redundant. Now that we have 10 imputed dataset. The next recipe of the magic is fit each dataset to a model then pool them

Fit All!

# put all 10 dataset into a list (probably should've done this before for efficienty, oh well)
imp_list <- map(.x=1:10, ~ missing_fill_full(empty_df_mice |> filter(chain == .x)))

# fit one at a time and save
fits <- map(imp_list, ~ glm(outcome ~ p1 + p2 + p3 + p4 + p5,
                            data = .x, family = binomial))

Let’s take a peek at the first fit

fits[[1]] |> broom::tidy(conf.int = T)

Alright! Looks pretty ok.

Pool them

# assess the point estimate of coefficient
coef <- 1:6
fit_list <- 1:10
temp_coef <- c()
# temp_coef_var <- c()
model_coef <- c()
model_coef_var <- c()
for (i in coef) {
  for (j in fit_list) {
  temp_coef[j] <- fits[[j]]$coefficients[[i]]
  }
  model_coef[i] <- mean(temp_coef)
  model_coef_var[i] <- var(temp_coef)
}

# assess the pooled SE to construct 95% CI of the coefficients
se <- do.call(rbind, map(.x = fits, ~ tidy(.x) |> pull(std.error))) |> as_tibble()
within_var <- colMeans(se^2) 
between_var <- model_coef_var 
total_var <- within_var + between_var + between_var / length(fit_list)
final_se <- sqrt(total_var)

# degrees of freedom (Barnard–Rubin)
lambda <- (1 + 1/length(fit_list)) * between_var / total_var
df_com <- 1000 - 6
df_old <- (length(fit_list) - 1) / lambda^2
df_obs <- (df_com + 1) / (df_com + 3) * df_com * (1 - lambda)
degfree <- df_old * df_obs / (df_old + df_obs)

tval <- qt(0.975, degfree)
mice_coef <- tibble(term = fits[[1]] |> tidy() |> pull(term), estimate = model_coef, std.error = final_se, tval = tval,
       conf.low = estimate - tval * std.error, conf.high = estimate + tval * std.error) |>
  select(-tval) |>
  mutate(method = "our_mice")

To combine our 10 models, we use Rubin’s rules. For the estimate, we just average the 10 coefficients. For the standard error, we add two kinds of uncertainty: the usual uncertainty within each model, and how much the 10 coefficients disagree with each other because of the missing data. Since that disagreement is estimated from only 10 datasets, we use a t value with Barnard–Rubin degrees of freedom instead of z score of 1.96 for 95% CI.

library(mice)

imp  <- mice(df_miss, m = 10, maxit = 10, seed = 1, printFlag = FALSE)
fit  <- with(imp, glm(outcome ~ p1 + p2 + p3 + p4 + p5, family = binomial))
actual_mice_coef <- tidy(pool(fit), conf.int = TRUE) |>
  select(term, estimate, std.error, conf.low, conf.high) |>
  mutate(method = "actual_mice")

For completeness sake, let’s use the professional’s method and see if our method matches up.

Grab All the Coefficients From Different Models

## function to grab stats
grab_stats <- function(df, method) {
  model <- glm(outcome ~ p1 + p2 + p3 + p4 + p5, df, family = "binomial") 
  coef <- model |> tidy(conf.int = T) |> select(-statistic, -p.value) |> mutate(method = method)
  return(coef)
}

## full data
full_data_coef <- grab_stats(df, "full_data")

## complete case
cco_coef <- grab_stats(df_miss, "complete_case")

## complete case analysis
all_coef <- rbind(mice_coef, full_data_coef, cco_coef, actual_mice_coef)

## truth
truth <- tibble(term = c("(Intercept)", "p1", "p2", "p3", "p4", "p5"),
                truth = c(-0.5, 0.3, -0.2, 0.1, 0.5, 0.4))

## visualize all
all_coef |>
  ggplot(aes(x = estimate, y = term, color = method)) +
  geom_vline(xintercept = 0, linetype = "dotted", colour = "grey60") +
  geom_pointrange(aes(xmin = conf.low, xmax = conf.high),
                  position = position_dodge(width = 0.6)) +
  geom_point(data = truth, aes(x = truth, y = term),
             inherit.aes = FALSE, shape = 4, size = 4) +
  labs(x = "Estimate (log-odds)", y = NULL,
       caption = "× = true coefficient") +
  theme_bw()

Our from-scratch MICE (purple) lines up closely with the mice package (red) for every coefficient, and both stay near the full-data fit (teal). Their intervals are slightly wider than the full-data ones, reflecting the extra uncertainty from the missing data. Complete-case analysis (green) drops every row with any missing value, so its estimates wander and its intervals are much wider. The × marks show the true coefficients from our data-generating process.

For p4, both MICE versions land a bit above the full-data estimate. Since they agree, it’s likely a quirk of this dataset rather than a coding problem. p5 is MNAR, yet both versions look fine here. That may just be luck in a single run, technically it shouldn’t really work, right?

Hooray! We did it! Wow, that was a long journey! I started this since 10/9/26 on the procedure, draft codes etc. Met lots of hurdles over the weeks. I’d be honest with you, it was bumpy! My initial go was a wrong approach where I average out the 10 datasets imputed data (which was wrong), then realized that that is not the proper pooling. Thanks to LLM being consistent tutor and back and forth conversaiton to converge to the proper technique. Now, I know how it works under the hood!

Opportunities For Improvement

• next we should learn sampling from posterior instead of resampling
• Try predictive mean matching, mice’s default. It borrows real observed values instead of adding normal noise, which keeps imputations within realistic ranges.
• Make the code faster and tidier. Vectorise missing_fill(), keep everything in lists instead of global variables, and pass the random generator explicitly
• try out sensitivity analysis with mice

Lessons learnt

• learnt the machinery of mice!
• learnt about map2_dfc
• learnt about Barnard–Rubin degrees of freedom and how to pool using Rubin’s rule.
• learnt total variance formula

If you like this article:

• please feel free to send me a comment or visit my other blogs
• please feel free to follow me on BlueSky, twitter, GitHub or Mastodon
• if you would like collaborate please feel free to contact me

Le prix de l’intelligence artificielle baisse plus vite que celui de presque tout ce que l’industrie a connu. Cette baisse est réelle, mais elle ne mesure pas ce que vous payez. Un projet d’IA en entreprise n’achète pas des jetons : il achète des tâches réussies, et chaque tâche réussie porte le coût de celles qui ont échoué. Plus le jeton devient bon marché, plus la rentabilité d’un projet se joue ailleurs que sur la facture d’API : dans son taux de réussite, et dans le travail humain qu’il faut pour rattraper ses échecs.

Le jeton n’est pas un baril de pétrole

Les chiffres de baisse sont spectaculaires. En septembre 2026, l’institut de recherche Epoch AI estimait que le coût pour atteindre un niveau de performance donné avait baissé d’environ 47 % par trimestre depuis 2023, soit une division par 13 chaque année. C’est six fois plus rapide que la baisse du coût du calcul informatique lui-même.

Il y a pourtant un détail qui compte : Epoch AI ne mesure pas le prix du jeton. Elle mesure ce que coûte l’atteinte d’un résultat sur un banc d’essai, et prend soin de distinguer ces deux grandeurs. La nuance n’est pas académique, car le jeton n’a rien d’une matière première.

Il n’a pas de prix universel. Un baril de pétrole se compare à un autre baril. Le jeton d’un modèle ne vaut pas celui d’un autre : les tarifs varient d’un facteur dix ou plus selon le modèle, et les forfaits exprimés en « jetons inclus » sont à peu près illisibles pour un acheteur.

Sa quantité n’est pas prévisible. Personne ne sait, avant de lancer une tâche, combien de jetons elle consommera. Le même modèle peut produire une réponse parfaite du premier coup, ou s’égarer dans des développements inutiles et des relances. D’un modèle à l’autre, l’écart est considérable. En février 2026, le cabinet d’évaluation Artificial Analysis relevait qu’à tarif strictement identique, Claude Opus 4.6 avait consommé deux fois plus de jetons de sortie que son prédécesseur Opus 4.5 pour passer la même batterie de tests (58 millions contre 29 millions). Même prix unitaire, deux fois plus de jetons de sortie à payer. Dans le même relevé, GPT-5.2 consommait 130 millions de jetons de sortie, plus du double d’Opus 4.6, pour une facture totale pourtant comparable : 2 304 dollars contre 2 486. Le prix du jeton et le nombre de jetons varient ensemble, et dans des sens imprévisibles. Pris isolément, aucun des deux ne dit ce que coûte le travail.

Une image venue de l’atelier

Dans l’industrie, ce problème est résolu depuis longtemps. Aucun directeur d’usine ne juge une fonderie sur le prix du kilo de métal. Il la juge sur le coût de la bonne pièce : ce que coûte une pièce conforme, rebuts et retouches compris. Une fonderie qui vend son métal deux fois moins cher, mais rebute une pièce sur quatre et mobilise un régleur pour reprendre les autres, coûte plus cher que sa concurrente. L’indicateur qui compte est le rendement au premier passage : la part des pièces qui sortent conformes du premier coup, sans retouche.

L’IA générative relève exactement de cette logique. Le jeton est le kilo de métal. Ce que l’entreprise achète, c’est la bonne pièce.

Le cas : deux devis pour 1 000 notes de frais

Prenons un cas simulé, mais réaliste. Le directeur administratif et financier d’une entreprise de taille intermédiaire veut confier à un agent IA le contrôle des notes de frais : environ 1 000 par mois. L’agent enchaîne trois étapes : lire le justificatif, vérifier la conformité à la politique de remboursement, rédiger la décision.

Deux offres sont sur la table. Le modèle « économique » est affiché à 0,50 € le million de jetons ; le modèle « frontière » à 5 €, dix fois plus. Sur le papier, la décision paraît évidente.

Un pilote d’un mois corrige une première fois ce calcul. Le modèle économique, moins capable, boucle davantage, relit, se reprend et produit des réponses plus longues : il consomme en moyenne 60 000 jetons par note, contre 15 000 pour le modèle frontière. Le coût de calcul par note s’établit donc à 0,03 € contre 0,075 €. L’écart de prix, qui était de un à dix, n’est plus que de un à deux et demi.

Ce n’est pourtant que la moitié de l’histoire. La question qui décide de tout reste ouverte : combien de ces notes sont justes ?

Chaque réussite paie pour les échecs

Un agent IA est probabiliste. Il peut se tromper sans le signaler, et produire une décision fausse avec la même assurance qu’une décision juste.

Sur une chaîne de trois étapes, les erreurs s’accumulent.

Si le modèle économique réussit chaque étape dans 90 % des cas, les trois étapes réussissent ensemble dans environ 73 % des cas (0,90 × 0,90 × 0,90). Pour le modèle frontière, fiable à 97 % par étape, on obtient environ 91 %. Les industriels connaissent ce calcul sous le nom de rendement cumulé : la conformité d’une ligne est le produit de la conformité de chacun de ses postes. Ce n’est qu’un ordre de grandeur, valable si les erreurs sont indépendantes d’une étape à l’autre ; seul le pilote dit le chiffre réel. Mais il explique pourquoi une fiabilité qui semble élevée à chaque étape s’effondre sur la chaîne entière.

Les ingénieurs résument la conséquence économique d’une formule : le coût d’une tâche réussie est le coût d’une tentative divisé par la probabilité de succès. Traduit en langage de gestion : si une note sur quatre échoue, chaque note réussie paie sa propre tentative, plus une part de celles qui ont raté.

Reste à savoir comment on paie un échec. Il existe deux manières.

• Relancer la machine. Le coût est dérisoire : pour le modèle économique, 0,03 € divisé par 73 % donne environ 0,04 € par note réussie. Si l’histoire s’arrêtait là, le modèle économique gagnerait.
• Faire reprendre la note par un comptable. Une décision de remboursement erronée ne se relance pas à l’aveugle : quelqu’un doit reprendre le dossier, vérifier, trancher. Comptons un quart d’heure à 60 € de l’heure, soit 15 € par reprise.

Dans la plupart des processus métier, c’est la seconde voie qui s’impose. Et elle change tout.

Le verdict : le modèle le moins cher coûte trois fois plus

Le modèle dix fois moins cher au jeton revient trois fois plus cher par note acceptée. Les jetons pèsent moins de 1 % du coût pour le modèle économique, environ 5 % pour le modèle frontière. Tout le reste correspond à ce que la facture d’API ne montre jamais : le travail humain qui rattrape les échecs.

Ce tableau est même optimiste. Il suppose que chaque échec est repéré. En pratique, pour trouver les notes fausses, il faut un contrôle, et ce contrôle a lui aussi un coût. Il ne compte pas non plus les serveurs, l’orchestration ni la maintenance des règles. On ne paie pas un modèle : on paie un dispositif.

La règle de décision tient en une ligne :

Coût d’une tâche réussie = tout ce qu’on a dépensé pour l’obtenir (tentatives, relances, reprises) ÷ nombre de tâches réussies.

C’est le coût de la bonne pièce, appliqué à l’IA.

La dérive silencieuse : quand la baisse de prix masque une hausse de coût

Supposons que le DAF ait fait le bon choix et retenu le modèle frontière. Six mois plus tard, le fournisseur publie une nouvelle version et divise son prix par deux. La facture d’API passe de 75 € à 37,50 € par mois. Le tableau de bord affiche une économie de 50 %.

Mais un système probabiliste ne tombe pas en panne : il s’érode. La nouvelle version se comporte un peu différemment, les justificatifs changent de format, la politique de remboursement évolue. Sans qu’aucune alerte ne se déclenche, le taux de réussite passe de 91 % à 85 %. Cela fait 150 reprises au lieu de 87, soit 2 250 € de travail humain. Le coût par note acceptée passe d’environ 1,38 € à 2,29 € : +66 %, au moment même où le prix du jeton a baissé de moitié.

Ces chiffres sont une illustration, pas une loi : une dérive peut être plus lente, plus rapide, ou nulle. C’est précisément la raison pour laquelle on ne peut pas la supposer et qu’il faut la mesurer. Celui qui ne suit que la facture de jetons célébrera une économie. Seule la mesure régulière du taux de réussite révèle la hausse.

Ce qu’il faut exiger avant de signer

La baisse du prix du jeton est une bonne nouvelle pour l’industrie. Pour un décideur, elle déplace la question : plus le calcul devient bon marché, plus le coût d’une tâche se réduit au coût de ses échecs. Trois questions suffisent à remettre un projet sur ses pieds.

1. Quel est le coût d’une tâche réussie, échecs et reprises compris ? Non le prix du jeton, ni le prix d’une requête.
2. Quel taux de réussite a été mesuré sur nos propres dossiers ? Non celui d’un banc d’essai public, ni celui qu’annonce le fournisseur.
3. Qui mesure ce taux chaque mois, et en particulier après chaque mise à jour du modèle ?

Ne construisez plus vos dossiers d’investissement sur le coût théorique d’une requête. Évaluez le dispositif complet, et exigez la mesure du coût unitaire par tâche validée. Dans une usine, personne n’accepterait de piloter une ligne au prix du kilo de métal. Il n’y a aucune raison de le faire avec l’IA.

Pour aller plus loin, le livre blanc gratuit La facture fantôme de l’IA examine pourquoi les gains annoncés ne produisent pas toujours les économies attendues. Le Volume 1 de SAGA-IA permet d’éprouver la viabilité économique de votre propre projet, avec des cas et des exercices corrigés.

Le cas des notes de frais est une simulation construite pour cet article ; ses montants et ses taux sont des conventions explicites, non des observations.

Ask a large language model a question about your data, and you’ll probably get a confident answer. Ask it again, and you might get a different one. Sometimes, neither answer is correct.

That’s a problem if you’re building an AI agent other people need to use. How do you get dependable results from it, and what does it take to make it ready for production?

On 19 November 2026, join Pedro Silva, Principal Data Scientist at Jumping Rivers, for this free webinar.

Register for the webinar

From a Working Agent to Production

One way to make an agent more reliable is to keep the language model away from the arithmetic. Pedro will show how to build agents in R where the model decides which function to call and R does the calculations. The agent answers questions about your data with existing R functions, so the numbers come from R, not from the model.

Getting an agent to work in a demo is the easy part. Once other people start using it, new questions come up. How do you manage separate conversations for different users? How do you test an agent whose responses may vary between runs? How do you keep a record of its actions, and which safeguards can you rely on?

Pedro will work through each of these and show what it takes to turn a working prototype into a tool other people can depend on.

Who Should Attend?

The webinar is for data scientists, R developers and anyone looking at building and deploying AI agents.

It will be useful whether you’re trying out agents for the first time or have an application you want to open up to more users.

Continued Learning Benefits

If you join more of our webinars, you can also get discounts:

• Attend two webinars and receive 20% off tickets to the AI in Production 2027 conference
• Attend more than two webinars and receive 20% off any Jumping Rivers public training course

Webinar Details

• Topic: Agents in Production: From a Chatbot to a Tool You Can Rely On
• Speaker: Pedro Silva
• Date: Thursday, 19 November 2026
• Time: 1:15pm UK time
• Location: Online
• Cost: Free

Register for the webinar

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