How long do wars last, on average? If a war such as that currently under way in Iran has lasted 74 days so far, how long do we expect it to last in total? For all sorts of reasons, inquiring minds are interested. Luckily there are some very well curated datasets out there, including the Correlates of War, that make it easy to answer these questions.

A caveat to all this applies that I am not a military historian, just an interested amateur. I’m very open to having mistakes of interpretation or method pointed out to me.

Distribution of wars’ durations

The Correlates of War data lets us see, for example, that this is the distribution (on a logarithmic scale) of durations of wars post-Napoleon:

You can see I’ve compared this to a log-normal distribution and found that it doesn’t have quite as fat tails as that. But that’s ok, I’m not too worried about the precise shape, because later on I’ll be using pretty straightforward empirical methods.

This data is only for inter-state wars, which are in contrast to intra-state (eg civil wars) and extra-state (eg with external non-state actors). As I’m interested in a reference population to compare the current USA-Israel-Iran war to, it’s the inter-state population I want.

The median length of a war is 139 days and the mean is 408 days.

The four day war in the dataset is the so-called “Football War” of 1969 between Honduras and El Salvador. The 3,734 day war was the much better-known “Vietnam War Phase II”, involving USA, Australia, Vietnam, Cambodia and others.

Here’s the code to import the data from the Correlates of War project and draw that first density plot:

library(tidyverse)
library(lubridate)
library(janitor)
library(glue)
library(ggrepel)
library(scales)

# https://correlatesofwar.org/data-sets/cow-war/


#----- import interstate war data----------------------

interstate <- read_csv("https://correlatesofwar.org/wp-content/uploads/Inter-StateWarData_v4.0.csv") |> 
  clean_names() |> 
  mutate(start_date = as.Date(sprintf("%04d-%02d-%02d", start_year1, start_month1, start_day1)),
         end_date = as.Date(sprintf("%04d-%02d-%02d", end_year1, end_month1, end_day1)))

interstate_wars <- interstate |> 
  group_by(war_num, war_name) |> 
  summarise(earliest_start= min(start_date),
            latest_end = max(end_date),
            bat_death = sum(bat_death)) |> 
  mutate(duration = as.numeric(latest_end - earliest_start),
         start_year = year(earliest_start)) |> 
  ungroup()

# what years covered? 1823 to 2003 at time of writing
range(interstate_wars$start_year)

#==========================plots=================
 
simple_caption <- "Source: Correlates of War, Inter-State War Data; analysis by freerangestats.info"

#-----------------distribution of duration------------
summary(interstate_wars$duration)

sim_norm <- data.frame(duration = 10 ^ (rnorm(1e6, 
                                        mean = log10(interstate_wars$duration), 
                                        sd = sd(log10(interstate_wars$duration)))))

interstate_wars |> 
  ggplot(aes(x = duration)) +
  geom_density() +
  geom_rug() +
  geom_density(data = sim_norm, colour = "orange") +
  annotate("text", x= 1, y = 0.18, label = "Simulated log-normal distribution", 
           colour = "orange", hjust = 0) +
  annotate("text", x= 300, y = 0.51, label = "Empirical distribution of war durations", 
           colour = "black", hjust = 0) +
  # carefully chosen labels for x axis:
  scale_x_log10(label = comma, breaks = c(range(interstate_wars$duration), 10, 100, 1000)) +
  labs(x = "Duration of wars (in days, logarithmic scale)",
       y = "Density",
       title = "Distribution of war durations, 1823 to 2003",
       subtitle = "More concentrated, less-fat tails than a log-normal distribution",
       caption = simple_caption) +
  # use coord to limit x axis so statistical calculations are all done on full data:
  coord_cartesian(xlim = c(1, 8000))

OK, so my main analytical task here is to work out the conditional expected duration of a war that has reached 74 days – the length so far of the USA-Israel-Iran war. Yes, I know there’s an incompletely observed ceasefire, but there’s also a blockade (or two), and that’s unambiguously an act of war under international law. So I’m counting the war as ongoing.

My chart to answer this question is this one:

What’s happening here is:

• the empirical cumulative distribution function of durations is the dark line – basically the cumulative frequency on the vertical axis, but expressed as a proportion.
• the grey line is a simple LOESS smoother of that cumulative frequency, useful for modelling values that aren’t exactly matched in the data.
• the red lines show the duration of the current war, and where it would fit in the distribution of 1823 to 2003 wars. It’s about 0.33 (defined in the code below as the variable current_cf), meaning that the current war is already longer than about 33% of wars.
• the horizontal blue line is half way in the vertical space between the horizontal red line and 1. Where it meets the smoothed line and drops a vertical blue line shows the expected median duration of a war that has gotten to this 0.33 point on the cumulative frequency.

So we see that of wars that get as long as 74 days, we expect the median total length to be 261 days. That’s a bit grim for those of us who think that even extending into June is going to be very bad indeed for the world economy, but it’s good to know. Of course, there’s plenty of wars that get to 74 days and then stop soon after, so there’s hope there too.

Here’s the code to do that bit of statistical inference and draw the chart:

#-------------------cumulative distribution--------------
interstate_cumulative <- interstate_wars |> 
  arrange(duration) |> 
  mutate(cumulative_freq = 1:n() / n()) 

# smoothed model of the cumulative distribution, including estimates of where
# the Iran war is on it:
model <- loess(cumulative_freq ~ log(duration), data = interstate_cumulative)
current_dur <- 74 # as at 13 May 2025 - war started 28 February 2026
current_cf <- predict(model, newdata = data.frame(duration = current_dur))

# inverse model to estimate duration given a cumulative frequency, useful for
# annotations on the chart:
inv_model <- loess(duration ~ x, 
                   data = data.frame(duration = interstate_cumulative$duration, 
                                     x = fitted(model)))

# of wars that last this long, what is the median cumulative frequency (i.e. half-way to 1):
conditional_median_freq <- (1 + current_cf) / 2
# of wars with that median cumulative frequency, convert it back into a duration,
conditional_median_dur <- predict(inv_model, data.frame(x = conditional_median_freq))

# Draw chart of cumulative distribution:
interstate_cumulative |> 
  ggplot(aes(x = duration, y = cumulative_freq)) +
  geom_smooth(method = "loess", colour = "grey80") +
  geom_line() +
  # note that (seems a bit odd) need to manually do the scale transform to geom_segment here:
  geom_segment(x = log10(current_dur), xend = log10(current_dur), y = -Inf, yend = current_cf, colour = "red") +
  geom_segment(x = 0, xend = log10(current_dur), y = current_cf, yend = current_cf, colour = "red") +
  geom_segment(x = log10(conditional_median_dur), xend = log10(conditional_median_dur), y = -Inf, yend = conditional_median_freq, colour = "blue") +
  geom_segment(x = 0, xend = log10(conditional_median_dur), y = conditional_median_freq, yend = conditional_median_freq, colour = "blue") +
  
  annotate("text", x = current_dur * 0.95, y = 0.39, label = "Current Iran war", colour = "red", hjust = 1) +
  annotate("text", x = conditional_median_dur * 1.05, y = 0.62, colour = "blue", hjust = 0, vjust = 1, 
           label = glue("Median expectation conditional 
on at least {current_dur} days")) +
  scale_x_log10(label = comma, breaks = c(10, current_dur, 100, conditional_median_dur, 1000)) +
  labs(x = "Total duration of war (in days, logarithmic scale)",
       y = "Cumulative frequency of wars",
       title = "Expectations of duration of Iran war, based on modern inter-state wars' duration",
       subtitle = glue("Comparison to wars from 1823 to 2003. The median war that lasts {current_dur} days goes on to last {round(conditional_median_dur)} days."),
       caption = simple_caption)

We can use the same approach to calculate not just the median war duration (conditional on getting to 74 days) but other percentiles. For example, in the below we can construct an 80% prediction interval (between the 0.1 and 0.9 quantiles) of total duration of 94.9 and 1,752 days. To put this another way, from this 74 day point, only 10% of wars will have a total duration of 94.9 or less days (ie another 21 days).

All up, that’s a big range of course; the main thing it tells us is that wars last longer than many people would like, and there’s a big variation in wars’ duration.

# some prediction intervals, conditional on getting to 74 days:
probs <- c(0.05, 0.1, 0.5, 0.8, 0.9, 0.95)
more_freqs <- probs * (1 - current_cf) + current_cf
conditional_dur <- predict(inv_model, data.frame(x = more_freqs))
tibble(probability = probs, duration = conditional_dur)
# so 80% of wars that reach 74 days will have a total duration between 95 and 1,752 days

  probability duration
        <dbl>    <dbl>
1        0.05     82.3
2        0.1      94.9
3        0.5     261. 
4        0.8    1141. 
5        0.9    1752. 
6        0.95   2119. 

Duration and other factors

So I’d answered my main question but I was naturally curious about some other relationships too. Obviously one expects longer wars to have more deaths in battle; can we see this in the data? Yes we can:

I like this chart as presenting the scale of nearly two centuries of inter-state war in one easy visualisation.

We also see that if there’s a pattern in relationship between duration, deaths and when the war started (the starting year mapped to colour in the chart above) it’s not an obvious one. We’ll come back to that in the next chart, but first, here’s the code to create the scatter plot above.

#------------------Compare duration and number of deaths----------------
interstate_wars |> 
  ggplot(aes(x = duration, y = bat_death, label = war_name)) +
  geom_point(aes(colour = start_year), size = 3.5) +
  geom_text_repel(colour = "grey50", size = 2, seed = 123) +
  scale_y_log10(label = comma) +
  scale_x_log10(label = comma) +
  scale_colour_viridis_c() +
  labs(title = "Inter-state wars, 1823-2003",
       colour = "Starting year",
       x = "Duration in days",
       y = "Number of battle deaths",
       caption = simple_caption) +
  theme(legend.position = c(0.15, 0.8))

I was a bit worried about that “two centuries” thing. Are recent wars all much shorter, or perhaps much longer, than older wars? If so it would be a big limitation on my inference about likely war length. So I prepared one more plot to check out if there was an obvious relationship, more rigorously than just eye-balling colour on the previous plot. I was a bit surprised to see that actually there is no real growth or reduction in war duration over time:

I also quite like this chart as giving us an instant comparison of our current USA-Israel-Iran war with some of those in history. We can see that it is already longer than the Boxer Rebellion, but not quite as long as the Falkland Islands or the War for Kosovo (for all of these names I am using those provided by the Correlates of War project - I’m well aware that these are contested labels).

Here’s my final chunk of code drawing that last chart:

#------------Compare duration with when in history it happened---------------
interstate_wars |> 
  arrange(bat_death) |> 
  ggplot(aes(x = earliest_start, y = duration)) +
  geom_hline(yintercept = current_dur, colour = "red") +
  geom_point(aes(size = bat_death), shape = 1) +
  geom_text_repel(aes(label = war_name), colour = "steelblue", size = 3, seed = 123) +
  annotate("text", x= as.Date("1820-01-01"), y = current_dur + 8, hjust = 0,
           label = "Duration of 2026 US-Israel-Iran war so far", colour = "red") +
  scale_y_log10(label = comma) +
  scale_size_area(label = comma, max_size = 25) +
  labs(title = "Inter-state wars, 1823-2003",
       subtitle = glue("Compared to the USA-Israel-Iran war as at {Sys.Date()}"),
       x = "Start of war",
       y = "Duration of war (days)",
       size = "Number of batlle deaths:",
       caption = simple_caption)

That’s all folks. Stay safe out there.

A high can make a regression model look impressively accurate — but this number can be deceptive. If you want to understand why a high is not always a sign of a good model, read on!

In the post, Learning Data Science: Modelling Basics, we built a simple model to predict income from age. R printed a model summary containing something called R-squared, but we did not yet discuss what that value actually means.

At first sight, a high looks highly reassuring. In our example, the linear model achieved an close to 90%. That sounds impressive.

However, just as high classification accuracy can be misleading — as discussed in ZeroR: The Simplest Possible Classifier, or Why High Accuracy can be Misleading — a high can also create a false sense of confidence.

To understand why, it helps to examine the formula itself and then revisit the three models from the previous post: the mean model, the linear model, and the polynomial model.

The Meaning of

The coefficient of determination is defined as:

At first glance, the formula appears intimidating, but its basic idea is relatively simple.

The denominator

measures the total variation in the target variable. It quantifies how strongly the observed values differ from their mean.

The numerator

measures the remaining unexplained error after fitting the model.

Thus, measures the proportion of variation explained by the model.

An of:

• 0 means the model explains none of the variation,
• 1 means the model explains all variation perfectly.

This sounds straightforward enough. The difficulty is that perfectly explaining the observed data is not necessarily the same thing as building a useful predictive model.

The Mean Model

Let us begin with the simplest possible regression model.

Suppose we completely ignore age and simply predict the average income for every individual:

This is effectively the regression equivalent of ZeroR. The model does not learn any relationship at all.

In this case:

Therefore, the residual sum of squares becomes identical to the total sum of squares:

Substituting this into the formula gives:

The model explains none of the variation in the data.

This corresponds to the underfitting case discussed previously: the model is too simple to capture the underlying structure.

The Polynomial Model

Now consider the opposite extreme.

Instead of fitting a straight line, suppose we fit a polynomial of sufficiently high degree. In fact, if we have observations with distinct age values, a polynomial of degree up to can pass exactly through all observed data points:

In that case:

for all observations, implying:

and therefore:

The model achieves a perfect fit.

At first sight, this appears ideal. In practice, however, such a model often performs poorly on unseen data because it has adapted itself not only to the underlying relationship, but also to random fluctuations and noise within the training data.

This is the classical overfitting problem.

A perfect may therefore indicate not a particularly good model, but a model that has become too flexible.

The Linear Model

The linear model from the previous post lies between these two extremes.

It is simple enough to avoid memorizing every random fluctuation, yet flexible enough to capture a meaningful trend in the data.

This balance between simplicity and flexibility is one of the central themes in statistical learning.

The idea was summarized in the previous post with the following plot:

and by the famous observation attributed to George Box:

“All models are wrong, but some are useful.”

The objective in modelling is therefore not to maximize complexity or maximize , but to find a model that generalizes well beyond the observed sample.

Why Alone Is Insufficient

The key limitation of is that it evaluates fit on the observed data only.

It does not directly measure:

• predictive performance on unseen data,
• robustness,
• causal validity, or
• generalization ability.

As model complexity increases, almost always increases as well. A sufficiently flexible model can often achieve values very close to 1 even when its predictions on new data are poor.

For this reason, practical data science relies on additional evaluation methods such as:

• train-test splits,
• cross-validation,
• regularization,
• adjusted , and
• out-of-sample testing.

The goal is not to reproduce historical observations perfectly, but to construct models that remain useful when confronted with new data.

A high can therefore mean two very different things:

• the model has identified a genuine structure,
• or the model has merely adapted itself too closely to the training data.

Distinguishing between these possibilities is one of the central challenges of machine learning and statistical modelling.

Expected goals has become one of the most important concepts in modern football analytics. Instead of judging a team only by goals scored, xG helps us estimate the quality of the chances created. In this tutorial, we will build a practical expected goals model in R using football data, feature engineering, logistic regression, model evaluation, and visualization.

This is a hands-on guide for analysts who want to move beyond simple football statistics and start building reproducible soccer analytics workflows in R.

What Is Expected Goals?

Expected goals, usually written as xG, measures the probability that a shot becomes a goal. A shot from two meters in front of goal will usually have a high xG value, while a long-range shot from outside the box will usually have a low xG value.

An xG model can use variables such as:

• Shot distance
• Shot angle
• Body part used
• Game state
• Minute of the match
• Shot type
• Set-piece situation
• Home or away context

In this post, we will build a clean starter model using R. You can later extend it with richer event data, tracking data, or more advanced machine learning models.

Install and Load R Packages

# Core data science packages
install.packages(c(
  "tidyverse",
  "ggplot2",
  "dplyr",
  "readr",
  "janitor",
  "broom",
  "yardstick",
  "rsample",
  "pROC",
  "patchwork"
))

# Football data package
install.packages("worldfootballR")

library(tidyverse)
library(ggplot2)
library(dplyr)
library(readr)
library(janitor)
library(broom)
library(yardstick)
library(rsample)
library(pROC)
library(patchwork)
library(worldfootballR)

Create a Simple Shot Dataset

Different public football data sources structure shot data differently. To make this tutorial reproducible, we will first create a synthetic shot dataset that behaves like real football event data. Later, you can replace this with your own data from FBref, StatsBomb open data, Wyscout-style exports, or custom event feeds.

set.seed(123)

n_shots <- 5000

shots <- tibble(
  shot_id = 1:n_shots,
  player = sample(
    c("Player A", "Player B", "Player C", "Player D", "Player E"),
    n_shots,
    replace = TRUE
  ),
  team = sample(
    c("Team Red", "Team Blue", "Team Green", "Team White"),
    n_shots,
    replace = TRUE
  ),
  minute = sample(1:95, n_shots, replace = TRUE),
  x_location = runif(n_shots, min = 70, max = 120),
  y_location = runif(n_shots, min = 0, max = 80),
  body_part = sample(
    c("Right Foot", "Left Foot", "Header", "Other"),
    n_shots,
    replace = TRUE,
    prob = c(0.43, 0.32, 0.20, 0.05)
  ),
  situation = sample(
    c("Open Play", "Corner", "Free Kick", "Penalty", "Counter Attack"),
    n_shots,
    replace = TRUE,
    prob = c(0.68, 0.12, 0.08, 0.03, 0.09)
  ),
  home_away = sample(c("Home", "Away"), n_shots, replace = TRUE)
)

glimpse(shots)

Engineer Shot Distance and Angle

Distance and angle are two of the most important features in a basic xG model. We will assume the goal is centered at x = 120 and y = 40.

goal_x <- 120
goal_y <- 40

shots <- shots %>%
  mutate(
    distance_to_goal = sqrt(
      (goal_x - x_location)^2 + (goal_y - y_location)^2
    ),
    angle_to_goal = atan2(
      abs(goal_y - y_location),
      goal_x - x_location
    ),
    angle_degrees = angle_to_goal * 180 / pi
  )

shots %>%
  select(shot_id, x_location, y_location, distance_to_goal, angle_degrees) %>%
  head()

Create a Goal Outcome

For demonstration, we will simulate goals using realistic football logic. Shots closer to goal should be more likely to become goals. Penalties should have higher probability. Headers and long-range attempts should usually be harder.

shots <- shots %>%
  mutate(
    linear_probability =
      -2.8 -
      0.08 * distance_to_goal +
      0.025 * angle_degrees +
      if_else(body_part == "Header", -0.35, 0) +
      if_else(body_part == "Other", -0.60, 0) +
      if_else(situation == "Penalty", 3.00, 0) +
      if_else(situation == "Counter Attack", 0.35, 0) +
      if_else(situation == "Free Kick", -0.45, 0),
    
    goal_probability = plogis(linear_probability),
    goal = rbinom(n(), size = 1, prob = goal_probability)
  )

shots %>%
  summarise(
    total_shots = n(),
    total_goals = sum(goal),
    conversion_rate = mean(goal)
  )

Explore the Shot Data

shots %>%
  count(body_part, goal) %>%
  group_by(body_part) %>%
  mutate(rate = n / sum(n))
shots %>%
  group_by(situation) %>%
  summarise(
    shots = n(),
    goals = sum(goal),
    conversion_rate = mean(goal),
    avg_distance = mean(distance_to_goal),
    .groups = "drop"
  ) %>%
  arrange(desc(conversion_rate))

Visualize Shot Locations

ggplot(shots, aes(x = x_location, y = y_location, color = factor(goal))) +
  geom_point(alpha = 0.35) +
  coord_fixed() +
  labs(
    title = "Shot Map",
    x = "Pitch Length",
    y = "Pitch Width",
    color = "Goal"
  ) +
  theme_minimal()

Split Data into Training and Testing Sets

set.seed(123)

shot_split <- initial_split(shots, prop = 0.80, strata = goal)

train_data <- training(shot_split)
test_data  <- testing(shot_split)

nrow(train_data)
nrow(test_data)

Build a Logistic Regression xG Model

Expected goals is naturally suited to logistic regression because the outcome is binary: goal or no goal.

xg_model <- glm(
  goal ~ distance_to_goal +
    angle_degrees +
    body_part +
    situation +
    home_away +
    minute,
  data = train_data,
  family = binomial()
)

summary(xg_model)

Convert Model Output into xG Values

test_predictions <- test_data %>%
  mutate(
    xg = predict(xg_model, newdata = test_data, type = "response")
  )

test_predictions %>%
  select(player, team, goal, xg, distance_to_goal, angle_degrees) %>%
  head(10)

Evaluate the xG Model

A good xG model should not only predict goals, but also produce well-calibrated probabilities. If 100 shots each have an xG of 0.10, we would expect roughly 10 goals over a large enough sample.

test_predictions %>%
  summarise(
    actual_goals = sum(goal),
    expected_goals = sum(xg),
    avg_xg = mean(xg),
    actual_conversion = mean(goal)
  )

ROC AUC

roc_obj <- roc(
  response = test_predictions$goal,
  predictor = test_predictions$xg
)

auc(roc_obj)
plot(
  roc_obj,
  main = "ROC Curve for xG Model"
)

Brier Score

brier_score <- mean((test_predictions$xg - test_predictions$goal)^2)

brier_score

Create xG Buckets for Calibration

calibration_table <- test_predictions %>%
  mutate(
    xg_bucket = cut(
      xg,
      breaks = seq(0, 1, by = 0.05),
      include.lowest = TRUE
    )
  ) %>%
  group_by(xg_bucket) %>%
  summarise(
    shots = n(),
    avg_xg = mean(xg),
    actual_goal_rate = mean(goal),
    goals = sum(goal),
    .groups = "drop"
  ) %>%
  filter(shots >= 10)

calibration_table
ggplot(calibration_table, aes(x = avg_xg, y = actual_goal_rate)) +
  geom_point(size = 3) +
  geom_abline(intercept = 0, slope = 1, linetype = "dashed") +
  labs(
    title = "xG Model Calibration",
    x = "Average Predicted xG",
    y = "Actual Goal Rate"
  ) +
  theme_minimal()

Player-Level xG Analysis

Once every shot has an xG value, we can aggregate by player. This allows us to compare goals, expected goals, overperformance, and shot volume.

player_xg <- test_predictions %>%
  group_by(player) %>%
  summarise(
    shots = n(),
    goals = sum(goal),
    xg = sum(xg),
    goals_minus_xg = goals - xg,
    xg_per_shot = mean(xg),
    conversion_rate = mean(goal),
    .groups = "drop"
  ) %>%
  arrange(desc(xg))

player_xg
ggplot(player_xg, aes(x = reorder(player, xg), y = xg)) +
  geom_col() +
  coord_flip() +
  labs(
    title = "Expected Goals by Player",
    x = "Player",
    y = "Total xG"
  ) +
  theme_minimal()

Team-Level xG Analysis

team_xg <- test_predictions %>%
  group_by(team) %>%
  summarise(
    shots = n(),
    goals = sum(goal),
    xg = sum(xg),
    goals_minus_xg = goals - xg,
    avg_xg_per_shot = mean(xg),
    .groups = "drop"
  ) %>%
  arrange(desc(xg))

team_xg
ggplot(team_xg, aes(x = reorder(team, goals_minus_xg), y = goals_minus_xg)) +
  geom_col() +
  coord_flip() +
  labs(
    title = "Goals Minus xG by Team",
    x = "Team",
    y = "Goals - Expected Goals"
  ) +
  theme_minimal()

Shot Quality Distribution

ggplot(test_predictions, aes(x = xg)) +
  geom_histogram(bins = 40) +
  labs(
    title = "Distribution of Shot Quality",
    x = "Expected Goals",
    y = "Number of Shots"
  ) +
  theme_minimal()

Compare Goals and xG by Situation

situation_xg <- test_predictions %>%
  group_by(situation) %>%
  summarise(
    shots = n(),
    goals = sum(goal),
    xg = sum(xg),
    avg_xg = mean(xg),
    .groups = "drop"
  ) %>%
  arrange(desc(avg_xg))

situation_xg
situation_long <- situation_xg %>%
  select(situation, goals, xg) %>%
  pivot_longer(
    cols = c(goals, xg),
    names_to = "metric",
    values_to = "value"
  )

ggplot(situation_long, aes(x = reorder(situation, value), y = value, fill = metric)) +
  geom_col(position = "dodge") +
  coord_flip() +
  labs(
    title = "Goals vs Expected Goals by Situation",
    x = "Situation",
    y = "Value",
    fill = "Metric"
  ) +
  theme_minimal()

Build a More Advanced xG Model with Interactions

A simple model is useful, but football is full of interactions. For example, distance may affect headers differently than footed shots. We can include interaction terms in the model.

xg_model_interaction <- glm(
  goal ~ distance_to_goal * body_part +
    angle_degrees +
    situation +
    home_away +
    minute,
  data = train_data,
  family = binomial()
)

summary(xg_model_interaction)
test_predictions_interaction <- test_data %>%
  mutate(
    xg_interaction = predict(
      xg_model_interaction,
      newdata = test_data,
      type = "response"
    )
  )

mean((test_predictions_interaction$xg_interaction - test_predictions_interaction$goal)^2)

Compare Two xG Models

model_comparison <- tibble(
  model = c("Basic Logistic Regression", "Interaction Logistic Regression"),
  brier_score = c(
    mean((test_predictions$xg - test_predictions$goal)^2),
    mean((test_predictions_interaction$xg_interaction - test_predictions_interaction$goal)^2)
  ),
  total_predicted_goals = c(
    sum(test_predictions$xg),
    sum(test_predictions_interaction$xg_interaction)
  ),
  actual_goals = c(
    sum(test_predictions$goal),
    sum(test_predictions_interaction$goal)
  )
)

model_comparison

Create a Reusable xG Prediction Function

predict_xg <- function(model, new_shots) {
  new_shots %>%
    mutate(
      predicted_xg = predict(
        model,
        newdata = new_shots,
        type = "response"
      )
    )
}

new_predictions <- predict_xg(xg_model, test_data)

head(new_predictions)

Create a Custom Shot Example

custom_shot <- tibble(
  distance_to_goal = 12,
  angle_degrees = 28,
  body_part = "Right Foot",
  situation = "Open Play",
  home_away = "Home",
  minute = 62
)

predict(
  xg_model,
  newdata = custom_shot,
  type = "response"
)

Use worldfootballR for Real Football Workflows

For real projects, you can use packages such as worldfootballR to collect football data from public sources and build reproducible analysis pipelines. The exact available columns depend on the source and endpoint, so always inspect your data before modeling.

library(worldfootballR)
library(tidyverse)

# Example: get FBref match results
# Adjust country, gender, season_end_year, and tier depending on your project

premier_league_results <- fb_match_results(
  country = "ENG",
  gender = "M",
  season_end_year = 2025,
  tier = "1st"
)

glimpse(premier_league_results)
premier_league_results %>%
  clean_names() %>%
  head()

If you are building a full football analytics pipeline with FBref, Transfermarkt, and Understat-style workflows, a more structured project template can save a lot of time. I cover that type of end-to-end workflow in Mastering Football Data with worldfootballR, especially for readers who want reusable R scripts, clean folders, and practical football data examples.

Example: Clean Match Results Data

clean_results <- premier_league_results %>%
  clean_names()

clean_results %>%
  glimpse()
# Example structure will depend on the returned data
# Always check column names first

names(clean_results)

Build a Match-Level Team Summary

# This is an example pattern.
# You may need to adjust column names depending on your data source.

team_summary_example <- clean_results %>%
  summarise(
    matches = n()
  )

team_summary_example

Save Your xG Model

Once you have trained a model, save it so you can reuse it later in reports, dashboards, APIs, or automated pipelines.

saveRDS(xg_model, "xg_model_logistic_regression.rds")

loaded_xg_model <- readRDS("xg_model_logistic_regression.rds")

predict(
  loaded_xg_model,
  newdata = custom_shot,
  type = "response"
)

Create an xG Report Table

xg_report <- test_predictions %>%
  group_by(team, player) %>%
  summarise(
    shots = n(),
    goals = sum(goal),
    xg = round(sum(xg), 2),
    goals_minus_xg = round(goals - sum(xg), 2),
    xg_per_shot = round(mean(xg), 3),
    .groups = "drop"
  ) %>%
  arrange(desc(xg))

xg_report
write_csv(xg_report, "xg_player_report.csv")

Create an xG Shot Map

ggplot(test_predictions, aes(x = x_location, y = y_location)) +
  geom_point(aes(size = xg, alpha = xg)) +
  coord_fixed() +
  labs(
    title = "xG Shot Map",
    x = "Pitch Length",
    y = "Pitch Width",
    size = "xG",
    alpha = "xG"
  ) +
  theme_minimal()

Create a High-Value Chances Table

big_chances <- test_predictions %>%
  filter(xg >= 0.30) %>%
  arrange(desc(xg)) %>%
  select(
    player,
    team,
    minute,
    body_part,
    situation,
    distance_to_goal,
    angle_degrees,
    xg,
    goal
  )

big_chances %>%
  head(20)

Model Improvement Ideas

This starter xG model can be improved in many ways. A professional football analytics workflow may include:

• More accurate shot coordinates
• Goalkeeper position
• Defender pressure
• Pass type before the shot
• Through balls and cutbacks
• Shot speed
• First-time shots
• Game state
• Team strength
• Player finishing history

Train an XGBoost-Style Model Later

Logistic regression is interpretable and a good starting point. For higher predictive performance, you can later compare it with random forests, gradient boosting, or Bayesian models.

# Example packages for future model upgrades
# install.packages(c("xgboost", "ranger", "tidymodels"))

library(tidymodels)

# A future tidymodels workflow could look like this:

xg_recipe <- recipe(
  goal ~ distance_to_goal + angle_degrees + body_part + situation + home_away + minute,
  data = train_data
) %>%
  step_dummy(all_nominal_predictors()) %>%
  step_normalize(all_numeric_predictors())

xg_recipe

Build a Tidymodels Logistic Regression Workflow

logistic_spec <- logistic_reg() %>%
  set_engine("glm") %>%
  set_mode("classification")

xg_workflow <- workflow() %>%
  add_recipe(xg_recipe) %>%
  add_model(logistic_spec)

xg_fit <- fit(
  xg_workflow,
  data = train_data %>%
    mutate(goal = factor(goal, levels = c(0, 1)))
)

xg_fit
tidy(xg_fit)

Predict Probabilities with Tidymodels

tidy_predictions <- predict(
  xg_fit,
  new_data = test_data,
  type = "prob"
) %>%
  bind_cols(test_data %>% mutate(goal = factor(goal, levels = c(0, 1))))

head(tidy_predictions)
tidy_predictions %>%
  roc_auc(
    truth = goal,
    .pred_1
  )

Turn xG into Match Insights

The real value of expected goals is not just predicting whether one shot becomes a goal. The value comes from aggregation. Once every shot has a probability, you can create match-level and season-level insights.

match_shots <- test_predictions %>%
  mutate(
    match_id = sample(1:100, n(), replace = TRUE)
  )

match_xg <- match_shots %>%
  group_by(match_id, team) %>%
  summarise(
    shots = n(),
    goals = sum(goal),
    xg = sum(xg),
    .groups = "drop"
  )

match_xg %>%
  arrange(match_id, desc(xg)) %>%
  head(20)

Find Teams Creating Better Chances

team_chance_quality <- test_predictions %>%
  group_by(team) %>%
  summarise(
    shots = n(),
    total_xg = sum(xg),
    avg_xg_per_shot = mean(xg),
    big_chances = sum(xg >= 0.30),
    low_quality_shots = sum(xg <= 0.05),
    .groups = "drop"
  ) %>%
  arrange(desc(avg_xg_per_shot))

team_chance_quality

Final Thoughts

Building an expected goals model in R is one of the best ways to learn football analytics because it combines data cleaning, feature engineering, statistical modeling, visualization, and interpretation. A simple logistic regression model can already teach you a lot about shot quality, player performance, and team attacking style.

From here, the next steps are clear: use richer football data, improve your features, compare different models, evaluate calibration, and build repeatable workflows that can be updated every week during the season.

Expected goals is not the final answer to football analysis, but it is one of the best starting points for serious soccer data science in R.

The post How to Build an Expected Goals (xG) Model in R with worldfootballR appeared first on R Programming Books.

In the new version of unifiedml available on CRAN, you can benchmark different models using k-fold cross-validation (section 1 of this blog post), and there’s a unified interface for predicting model probabilities (section 2 of this blog post).

install.packages("unifiedml")

install.packages(c("e1071", "randomForest", "caret"))

install.packages("glmnet")

library(unifiedml)

1 – Benchmarking models

set.seed(123)

X <- iris[, 1:4]
y <- iris$Species

models <- list( # `Model` is exported from package 'unifiedml'
  glm  = Model$new(caret::train), # caret can be used (see https://topepo.github.io/caret/available-models.html)
  rf   = Model$new(randomForest::randomForest), # or a native pkg
  svm  = Model$new(e1071::svm) # or another pkg
)

params <- list(
  glm = list(method = "glmnet",
             tuneGrid = data.frame(alpha = 0, lambda = 0.01), # for caret model, all hyperparameters must be provided
             trControl = trainControl(method = "none")),
  rf  = list(ntree = 150), # Not necessarily needing to specify all hyperparameters
  svm = list(kernel = "radial",
             cost = 1,
             gamma = 0.1)
)

results <- unifiedml::benchmark(models, X, y, cv = 5, params = params)

[1/3] Fitting model: glm
Mean CV score for glm: 0.9533

[2/3] Fitting model: rf
Mean CV score for rf: 0.9600

[3/3] Fitting model: svm
Mean CV score for svm: 0.9733

print(results) # 5-fold cross-validation results

$glm
$glm$avg_score
[1] 0.9533333

$glm$scores
    fold1     fold2     fold3     fold4     fold5 
0.9333333 0.9666667 0.9333333 0.9333333 1.0000000 


$rf
$rf$avg_score
[1] 0.96

$rf$scores
    fold1     fold2     fold3     fold4     fold5 
0.9333333 1.0000000 0.9333333 0.9333333 1.0000000 


$svm
$svm$avg_score
[1] 0.9733333

$svm$scores
    fold1     fold2     fold3     fold4     fold5 
0.9666667 1.0000000 0.9666667 0.9333333 1.0000000 

# initialize empty vectors
model_vec <- c()
fold_vec  <- c()
score_vec <- c()

for (model in names(results)) {
  scores <- results[[model]]$scores

  model_vec <- c(model_vec, rep(model, length(scores)))
  fold_vec  <- c(fold_vec, names(scores))
  score_vec <- c(score_vec, as.numeric(scores))
}

df <- data.frame(
  model = model_vec,
  fold  = fold_vec,
  score = score_vec
)

library(ggplot2)

ggplot(df, aes(x = model, y = score, fill = model)) +
  geom_violin(trim = FALSE, alpha = 0.6) +
  geom_jitter(width = 0.08, size = 2) +
  theme_minimal() +
  labs(
    title = "Cross-validation score distribution",
    x = "Model",
    y = "Score"
  ) +
  theme(legend.position = "none")

2 - Unified interface for predicting probabilities

# Load required packages
library(unifiedml)
library(randomForest)
library(nnet)
library(e1071)

# Load iris dataset
data(iris)

# Setup reproducible data
set.seed(42)

# Create feature matrix (all 4 numeric features)
X <- as.matrix(iris[, 1:4])
colnames(X) <- c("Sepal.Length", "Sepal.Width", "Petal.Length", "Petal.Width")

# Target: Species (multi-class with 3 levels)
y_multiclass <- iris$Species

# Create binary classification target (Versicolor vs others)
y_binary <- factor(
  ifelse(iris$Species == "versicolor", "versicolor", "other"),
  levels = c("other", "versicolor")
)

# Split into train/test (75% train, 25% test)
set.seed(42)
train_idx <- sample(1:nrow(X), size = floor(0.75 * nrow(X)), replace = FALSE)
test_idx <- setdiff(1:nrow(X), train_idx)

X_train <- X[train_idx, ]
X_test <- X[test_idx, ]
y_train_multiclass <- y_multiclass[train_idx]
y_test_multiclass <- y_multiclass[test_idx]
y_train_binary <- y_binary[train_idx]
y_test_binary <- y_binary[test_idx]

cat("\n")
cat("============================================================================\n")
cat("IRIS DATASET - Summary\n")
cat("============================================================================\n")
cat(sprintf("Training samples: %d\n", nrow(X_train)))
cat(sprintf("Test samples: %d\n", nrow(X_test)))
cat(sprintf("Features: %d\n", ncol(X_train)))
cat(sprintf("Classes: %s\n", paste(levels(y_multiclass), collapse = ", ")))

# ============================================================================
# EXAMPLE 1: randomForest - Multi-class Classification on IRIS
# ============================================================================

cat("\n")
cat("============================================================================\n")
cat("EXAMPLE 1: randomForest - Multi-class Classification\n")
cat("============================================================================\n")

mod_rf <- Model$new(randomForest::randomForest)
mod_rf$fit(X_train, y_train_multiclass, ntree = 100)

cat("\nPredicting probabilities for first 5 test samples:\n")
probs_rf <- mod_rf$predict_proba(X_test[1:5, ])

cat("\nProbability matrix:\n")
print(round(probs_rf, 3))

cat("\nInterpretation:\n")
for(i in 1:5) {
  cat(sprintf("\nSample %d (Actual: %s):\n", i, as.character(y_test_multiclass[i])))
  cat(sprintf("  setosa:     %.1f%%\n", probs_rf[i, "setosa"] * 100))
  cat(sprintf("  versicolor: %.1f%%\n", probs_rf[i, "versicolor"] * 100))
  cat(sprintf("  virginica:  %.1f%%\n", probs_rf[i, "virginica"] * 100))
  cat(sprintf("  Predicted:  %s\n", colnames(probs_rf)[which.max(probs_rf[i, ])]))
}

# Get class predictions
pred_classes_rf <- mod_rf$predict(X_test[1:5, ], type = "class")
cat("\nPredicted classes (first 5):", as.character(pred_classes_rf), "\n")
cat("Actual classes (first 5):   ", as.character(y_test_multiclass[1:5]), "\n")

# Calculate accuracy on full test set
probs_all_rf <- mod_rf$predict_proba(X_test)
pred_all_rf <- colnames(probs_all_rf)[apply(probs_all_rf, 1, which.max)]
accuracy_rf <- mean(pred_all_rf == as.character(y_test_multiclass))
cat(sprintf("\nTest set accuracy: %.1f%%\n", accuracy_rf * 100))

# ============================================================================
# EXAMPLE 2: nnet - Multi-class Classification on IRIS
# ============================================================================

cat("\n")
cat("============================================================================\n")
cat("EXAMPLE 2: nnet - Multi-class Classification\n")
cat("============================================================================\n")

mod_nnet <- Model$new(nnet::nnet)
mod_nnet$fit(X_train, y_train_multiclass, size = 10, maxit = 200, trace = FALSE)

cat("\nPredicting probabilities for first 5 test samples:\n")
probs_nnet <- mod_nnet$predict_proba(X_test[1:5, ])

cat("\nProbability matrix (all 3 classes):\n")
print(round(probs_nnet, 3))

cat("\nDetailed predictions:\n")
for(i in 1:5) {
  cat(sprintf("\nSample %d (Actual: %s):\n", i, as.character(y_test_multiclass[i])))
  cat(sprintf("  setosa:     %.1f%%\n", probs_nnet[i, "setosa"] * 100))
  cat(sprintf("  versicolor: %.1f%%\n", probs_nnet[i, "versicolor"] * 100))
  cat(sprintf("  virginica:  %.1f%%\n", probs_nnet[i, "virginica"] * 100))
  cat(sprintf("  Predicted:  %s\n", colnames(probs_nnet)[which.max(probs_nnet[i, ])]))
}

# Get class predictions
pred_classes_nnet <- mod_nnet$predict(X_test[1:5, ], type = "class")
cat("\nPredicted classes (first 5):", as.character(pred_classes_nnet), "\n")
cat("Actual classes (first 5):   ", as.character(y_test_multiclass[1:5]), "\n")

# Calculate accuracy
probs_all_nnet <- mod_nnet$predict_proba(X_test)
pred_all_nnet <- colnames(probs_all_nnet)[apply(probs_all_nnet, 1, which.max)]
accuracy_nnet <- mean(pred_all_nnet == as.character(y_test_multiclass))
cat(sprintf("\nTest set accuracy: %.1f%%\n", accuracy_nnet * 100))

# ============================================================================
# EXAMPLE 3: SVM - Multi-class Classification on IRIS
# ============================================================================

cat("\n")
cat("============================================================================\n")
cat("EXAMPLE 3: SVM - Multi-class Classification\n")
cat("============================================================================\n")

mod_svm <- Model$new(e1071::svm)
mod_svm$fit(X_train, y_train_multiclass, probability = TRUE, kernel = "radial")

cat("\nPredicting probabilities for first 5 test samples:\n")
probs_svm <- mod_svm$predict_proba(X_test[1:5, ])

cat("\nProbability matrix:\n")
print(round(probs_svm, 4))

cat("\nDetailed predictions:\n")
for(i in 1:5) {
  cat(sprintf("\nSample %d (Actual: %s):\n", i, as.character(y_test_multiclass[i])))
  cat(sprintf("  setosa:     %.1f%%\n", probs_svm[i, "setosa"] * 100))
  cat(sprintf("  versicolor: %.1f%%\n", probs_svm[i, "versicolor"] * 100))
  cat(sprintf("  virginica:  %.1f%%\n", probs_svm[i, "virginica"] * 100))
  cat(sprintf("  Predicted:  %s\n", colnames(probs_svm)[which.max(probs_svm[i, ])]))
}

# Calculate accuracy
probs_all_svm <- mod_svm$predict_proba(X_test)
pred_all_svm <- colnames(probs_all_svm)[apply(probs_all_svm, 1, which.max)]
accuracy_svm <- mean(pred_all_svm == as.character(y_test_multiclass))
cat(sprintf("\nTest set accuracy: %.1f%%\n", accuracy_svm * 100))

============================================================================
IRIS DATASET - Summary
============================================================================
Training samples: 112
Test samples: 38
Features: 4
Classes: setosa, versicolor, virginica

============================================================================
EXAMPLE 1: randomForest - Multi-class Classification
============================================================================

Predicting probabilities for first 5 test samples:

Probability matrix:
  setosa versicolor virginica
1      1          0         0
2      1          0         0
3      1          0         0
4      1          0         0
5      1          0         0
attr(,"assign")
[1] 1 1 1
attr(,"contrasts")
attr(,"contrasts")$pred
[1] "contr.treatment"

attr(,"extraction_method")
[1] "fallback::1"
attr(,"model_class")
[1] "randomForest.formula"

Interpretation:

Sample 1 (Actual: setosa):
  setosa:     100.0%
  versicolor: 0.0%
  virginica:  0.0%
  Predicted:  setosa

Sample 2 (Actual: setosa):
  setosa:     100.0%
  versicolor: 0.0%
  virginica:  0.0%
  Predicted:  setosa

Sample 3 (Actual: setosa):
  setosa:     100.0%
  versicolor: 0.0%
  virginica:  0.0%
  Predicted:  setosa

Sample 4 (Actual: setosa):
  setosa:     100.0%
  versicolor: 0.0%
  virginica:  0.0%
  Predicted:  setosa

Sample 5 (Actual: setosa):
  setosa:     100.0%
  versicolor: 0.0%
  virginica:  0.0%
  Predicted:  setosa

Predicted classes (first 5): setosa setosa setosa setosa setosa 
Actual classes (first 5):    setosa setosa setosa setosa setosa 

Test set accuracy: 94.7%

============================================================================
EXAMPLE 2: nnet - Multi-class Classification
============================================================================

Predicting probabilities for first 5 test samples:

Probability matrix (all 3 classes):
  setosa versicolor virginica
1      1          0         0
2      1          0         0
3      1          0         0
4      1          0         0
5      1          0         0
attr(,"extraction_method")
[1] "fallback::5"
attr(,"model_class")
[1] "nnet.formula"

Detailed predictions:

Sample 1 (Actual: setosa):
  setosa:     100.0%
  versicolor: 0.0%
  virginica:  0.0%
  Predicted:  setosa

Sample 2 (Actual: setosa):
  setosa:     100.0%
  versicolor: 0.0%
  virginica:  0.0%
  Predicted:  setosa

Sample 3 (Actual: setosa):
  setosa:     100.0%
  versicolor: 0.0%
  virginica:  0.0%
  Predicted:  setosa

Sample 4 (Actual: setosa):
  setosa:     100.0%
  versicolor: 0.0%
  virginica:  0.0%
  Predicted:  setosa

Sample 5 (Actual: setosa):
  setosa:     100.0%
  versicolor: 0.0%
  virginica:  0.0%
  Predicted:  setosa

Predicted classes (first 5): setosa setosa setosa setosa setosa 
Actual classes (first 5):    setosa setosa setosa setosa setosa 

Test set accuracy: 97.4%

============================================================================
EXAMPLE 3: SVM - Multi-class Classification
============================================================================

Predicting probabilities for first 5 test samples:

Probability matrix:
  setosa versicolor virginica
1      1          0         0
2      1          0         0
3      1          0         0
4      1          0         0
5      1          0         0
attr(,"assign")
[1] 1 1 1
attr(,"contrasts")
attr(,"contrasts")$pred
[1] "contr.treatment"

attr(,"extraction_method")
[1] "fallback::1"
attr(,"model_class")
[1] "svm.formula"

Detailed predictions:

Sample 1 (Actual: setosa):
  setosa:     100.0%
  versicolor: 0.0%
  virginica:  0.0%
  Predicted:  setosa

Sample 2 (Actual: setosa):
  setosa:     100.0%
  versicolor: 0.0%
  virginica:  0.0%
  Predicted:  setosa

Sample 3 (Actual: setosa):
  setosa:     100.0%
  versicolor: 0.0%
  virginica:  0.0%
  Predicted:  setosa

Sample 4 (Actual: setosa):
  setosa:     100.0%
  versicolor: 0.0%
  virginica:  0.0%
  Predicted:  setosa

Sample 5 (Actual: setosa):
  setosa:     100.0%
  versicolor: 0.0%
  virginica:  0.0%
  Predicted:  setosa

Test set accuracy: 94.7%

Great strides in artificial intelligence development during the last five years produced agents that are now commonplace at work and home. It is humbling to note that virtually all frontier large language models today trace back to a preprint introducing the transformer neural network architecture¹ – a fifteen-page paper that profoundly rocked the world through waves of excitement and angst.

This paradigm shift in model design has also heavily influenced computer vision, leading to a surge in vision-language models (VLMs). Not only can such systems easily generalize across tasks such as segmentation, depth estimation and image generation or editing², they have also blown legacy models out of the water in object detection benchmarks, with little to no fine-tuning³.

However, it should not be lightly assumed that the transformer architecture is the only path forward to a more meaningful, cost-effective or even better-performing AI – not when we are still having trouble counting “r” in the word strawberry. Neuromorphic computation⁴, photonic neural networks⁵, JEPA⁶ and many other techniques have recently shown us different ways to design and implement intelligent systems that produce optimal solutions for a variety of problems.

Today I want to focus on a topic from a timeless book that inspired me to think differently and, particularly, to effectively apply a first principles approach to problem-solving. The topic is edge detection, and that book is Vision, by David Marr⁷. Just as On the Origin of Species⁸ and On Growth and Form⁹, this is yet another masterpiece that brought together different disciplines – in this case neurophysiology and computer vision – to revolutionise science.

In this blog post we will define and compare algorithms for image edge detection, and explore their remarkable similarity with neurophysiological readings.

Introduction

Modern computer vision is deeply rooted in Marr’s pioneering work. To understand any information-processing system, Marr argued, one must describe it at three interdependent levels of analysis:

• The computational level – what problem is being solved and why (e.g. edge detection)
• The algorithmic level – how it is solved, and what representations and procedures are used (e.g. the Laplacian transform)
• The implementational level – where it is physically realised (e.g. in vivo, in silico)

This layered thinking is what makes the book so enduring. Marr was not merely describing the visual system, he was arguing that to truly understand it you had to explain it at all three levels simultaneously. The book also features memorable passages on random dot stereograms¹⁰, binocular disparity and motion perception – overall, highly recommended for science enthusiasts.

Let us now introduce the key concepts underlying edge detection that leveraged this structured approach, to gain a better understanding of how it can be solved in practice.

Zero-crossing

From a computational perspective, edges are fundamentally spatial discontinuities in images. If for a brief moment we consider a simple greyscale image, edges are wherever dark-to-light and light-to-dark transitions occur, in whatever direction.

Because such transitions mathematically translate to local changes in pixel intensity, the most natural approach to identify edges is to compute image gradients, the two-dimensional equivalent of the derivative. The first derivative of image intensity evaluated across an edge produces a peak (for a dark-to-bright transition) or a trough (for a bright-to-dark transition), depending on the direction. However, the second derivative provides not only the means to identify both transition types, but also a beautifully simple detection mechanism: it crosses zero at the precise location of the edge. This is the essence of the zero-crossing.

First and second derivatives of a one-dimensional signal. Left: signal displaying low-to-high (dark-to-bright) transition. Middle: first derivative of the signal, capturing that transition as a peak. Right: second derivative of the signal, exhibiting the zero-crossing.

Marr and Hildreth formalised this insight by proposing the Laplacian of Gaussian (LoG) as the operator of choice¹¹. The Laplacian is the sum of second partial derivatives in both spatial dimensions:

Applied directly to a noisy image, the Laplacian amplifies every small intensity fluctuation. The Gaussian pre-filter solves this by smoothing the image at a chosen scale before differentiation. Because convolution is associative, the two steps can be combined into a single kernel – the LoG, also denoted :

This kernel, which resembles an inverted sombrero and is sometimes called the Mexican hat wavelet, produces a response that crosses zero exactly at an intensity edge. The width of the Gaussian determines the scale of detection: small preserves fine detail, large captures only coarse structure. Marr argued that the visual system operates simultaneously at multiple scales – an idea that would later resonate in scale-space theory and, much later, in the multi-scale feature maps of deep convolutional networks.

Marr went further and showed that the centre-surround organisation of retinal ganglion cell receptive fields – which he modelled as a Difference of Gaussians (DoG), the difference between a narrow excitatory Gaussian and a broader inhibitory one – is a close biological approximation of the LoG. Put differently, your retina is already computing zero-crossings before the signal ever reaches the visual cortex. The agreement between computational predictions and in vivo electrophysiological measurements, documented in Vision (p. 64), remains one of the most compelling examples of theory meeting experiment in all of neuroscience.

Zero-crossings are theoretically elegant, but the workhorse operators most computer-vision tools reach for – including, as we will see, the Canny detector itself – operate on first-derivative gradients. Let us look at those.

Image gradients

In practice, image gradients are computed using convolution filters – small kernels that slide across the image and produce a weighted local sum at each pixel, as illustrated in my post on convolutional neural networks. The two most widely used first-order gradient operators are:

Sobel: weights the central row and column more heavily, providing a modest degree of smoothing:

Prewitt: weights all neighbours equally:

In both cases the gradient magnitude at each pixel is , and the gradient direction is . Where the magnitude is large, a transition is occurring; where it is small, the neighbourhood is uniform.

The limitation of these operators is that they are sensitive to noise and produce thick, diffuse edges. Every pixel with a large gradient is flagged regardless of whether it truly lies on the edge or merely near it. This is precisely the problem that John Canny set out to solve.

Canny edge detection

Published in 1986, John Canny’s paper A Computational Approach to Edge Detection¹² remains one of the most cited works in computer vision. Canny framed edge detection as an explicit optimisation problem and derived a detector that simultaneously maximises three criteria: i) good detection (few missed edges, few “false alarms”), ii) good localisation (detected edges close to true edges), and iii) single response (one response per edge, not many). The resulting algorithm is a four-step pipeline outlined below:

Step 1 – Gaussian smoothing

As with the LoG, the first step is to suppress noise by convolving the image with a Gaussian kernel of width . The choice of directly governs the trade-off between noise suppression and fine detail preservation. A larger removes more noise but blurs genuine edges.

Step 2 – Gradient computation

The smoothed image is then differentiated, typically using Sobel kernels – to obtain the gradient magnitude and direction at every pixel.

Step 3 – Non-maximum suppression (NMS)

This step thins the edges. For each pixel, Canny checks whether its gradient magnitude is a local maximum along the gradient direction – that is, whether it is larger than its two neighbours in the direction . If it is not, it is suppressed to zero. The result is a set of thin, one-pixel-wide candidate edges.

Step 4 – Hysteresis thresholding

The final step uses two thresholds, and , to prune candidate pixels. A pixel whose gradient magnitude exceeds is accepted as an edge, and conversely a pixel below is rejected. A pixel between the two thresholds is accepted only if it is connected, directly or through other such pixels, to a strong edge pixel. This connectivity analysis – the defining feature of hysteresis – ensures that long, continuous edges are preserved even when their local gradient fluctuates, while isolated noise responses are discarded. For a more visual understanding of NMS and hysteresis I recommend reading the Canny edge detection documentation from OpenCV.

The elegance of Canny lies in how each step addresses a specific failure mode of earlier operators. In essence Gaussian smoothing tackles noise, NMS tackles thick edges and hysteresis tackles the false edge / broken edge trade-off that a single threshold cannot solve.

Let’s get started with Python

Time to practice! We will first build the separate components (Gaussian blur, Sobel gradients, LoG zero-crossings), then run the full Canny pipeline and explore how its parameters trade off recall against noise. We will use opencv and scikit-image alongside the usual suspects numpy and matplotlib. You can install all packages using the following shell command:

pip install opencv-python scikit-image matplotlib numpy

Image loading and preprocessing

We start by loading a greyscale image. For demonstration purposes I use a stock picture from scikit-image – feel free to use any other image of your choice.

import cv2import numpy as npimport matplotlib.pyplot as pltfrom skimage import datafrom scipy.ndimage import gaussian_laplace# Load a greyscale test image (uint8, values 0–255)image = data.camera() fig, ax = plt.subplots(figsize=(5, 5))ax.imshow(image, cmap='gray')ax.set_title('Original image')ax.axis('off')plt.tight_layout()plt.show()

The Laplacian of Gaussian and zero-crossings

Let us inspect the LoG response and its zero-crossings – the theoretical backbone we discussed earlier.

# LoG: positive sigma = apply Gaussian of that std, then Laplacianlog_response = gaussian_laplace(image.astype(float), sigma=2.0) # Zero-crossings: sign changes between neighbouring pixelsdef zero_crossings(log_img):    """Return a binary mask of zero-crossing locations."""    zc = np.zeros_like(log_img, dtype=bool)    # Check horizontal and vertical sign changes    for shift in [(0, 1), (1, 0)]:        shifted = np.roll(log_img, shift=shift, axis=(0, 1))        zc |= (np.sign(log_img) != np.sign(shifted))    return zc zc_mask = zero_crossings(log_response) fig, axes = plt.subplots(1, 2, figsize=(10, 4))axes[0].imshow(log_response, cmap='RdBu_r')axes[0].set_title('LoG response (σ=2.0)')axes[0].axis('off')axes[1].imshow(zc_mask, cmap='gray')axes[1].set_title('Zero-crossings')axes[1].axis('off')plt.tight_layout()plt.show()

Notice how the zero-crossing map already captures much of the scene’s edge structure, but it is sensitive to low-level noise and retains spurious responses in flat regions. This motivates the additional refinement steps of the Canny algorithm.

Gaussian smoothing and gradient computation

Before running the full Canny pipeline, it is instructive to inspect the intermediate steps. Here we apply a Gaussian blur and then compute Sobel gradients manually.

# Gaussian blur — sigma controlled by ksize (must be odd) and sigmaXblurred = cv2.GaussianBlur(image, ksize=(5, 5), sigmaX=1.4) # Sobel gradients in x and yGx = cv2.Sobel(blurred, cv2.CV_64F, dx=1, dy=0, ksize=3)Gy = cv2.Sobel(blurred, cv2.CV_64F, dx=0, dy=1, ksize=3) # Gradient magnitudemagnitude = np.sqrt(Gx**2 + Gy**2)magnitude = (magnitude / magnitude.max() * 255).astype(np.uint8) fig, axes = plt.subplots(1, 2, figsize=(10, 4))for ax, img, title in zip(axes, [blurred, magnitude], ['Gaussian blur (σ=1.4)', '|∇f| - Sobel magnitude']):    ax.imshow(img, cmap='gray')    ax.set_title(title)    ax.axis('off')plt.tight_layout()plt.show()

Canny edge detection with OpenCV

The OpenCV Canny() function accepts the image, the two hysteresis thresholds, and an optional aperture size for the Sobel operator. Crucially, the Gaussian smoothing step should be applied manually beforehand so you have full control over .

def run_canny(image, sigma, t_low, t_high, aperture=3):    """Apply Gaussian blur then Canny edge detection."""    # Kernel size: 2 * ceil(3*sigma) + 1 ensures the kernel covers ±3σ    ksize = 2 * int(np.ceil(3 * sigma)) + 1    blurred = cv2.GaussianBlur(image, (ksize, ksize), sigmaX=sigma)    edges = cv2.Canny(blurred, threshold1=t_low, threshold2=t_high, apertureSize=aperture)    return edgesedges = run_canny(image, sigma=1.4, t_low=50, t_high=150)fig, axes = plt.subplots(1, 2, figsize=(10, 4))axes[0].imshow(image, cmap='gray')axes[0].set_title('Original')axes[0].axis('off')axes[1].imshow(edges, cmap='gray')axes[1].set_title('Canny edges (σ=1.4, T_low=50, T_high=150)')axes[1].axis('off')plt.tight_layout()plt.show()

The effect of hysteresis thresholds

The ratio is at least as important as the absolute values. A common rule of thumb is to use a 2:1 or 3:1 ratio. Let us explore this now.

fig, axes = plt.subplots(2, 3, figsize=(14, 8))configs = [    (1.4, 20, 60, '↑ Recall, ↑ Noise\nσ=1.4, T=20/60'),    (1.4, 50, 150, 'Balanced\nσ=1.4, T=50/150'),    (1.4, 100, 200, '↓ Recall, ↓ Noise\nσ=1.4, T=100/200'),    (0.5, 50, 150, 'Fine scale\nσ=0.5, T=50/150'),    (2.0, 50, 150, 'Coarse scale\nσ=2.0, T=50/150'),    (4.0, 50, 150, 'Very coarse scale\nσ=4.0, T=50/150'),] for ax, (sigma, tl, th, title) in zip(axes.ravel(), configs):    result = run_canny(image, sigma=sigma, t_low=tl, t_high=th)    ax.imshow(result, cmap='gray')    ax.set_title(title, fontsize=9)    ax.axis('off') plt.tight_layout()plt.show()

The top row demonstrates the threshold effect – lower thresholds recover more edges but also more noise, higher thresholds yield cleaner output at the cost of broken contours. The bottom row shows the scale effect governed by – at small scales the detector responds to fine texture and noise, at large scales only the dominant structural boundaries survive.

Overlaying edges on the original image

A useful visualisation is overlaying detected edges on the original image. This facilitates the quality assessment of our workflow.

# Convert to RGB so we can draw edges in redoverlay = cv2.cvtColor(image, cv2.COLOR_GRAY2RGB)overlay[edges > 0] = [220, 40, 40] # red edges fig, ax = plt.subplots(figsize=(6, 6))ax.imshow(overlay)ax.set_title('Canny edges overlaid (σ=1.4, T=50/150)')ax.axis('off')plt.tight_layout()plt.show()

Similarly to the butterfly picture, here we find a solution that precisely identifies the sharpest edges from the image.

Conclusion

We have traced a line from the centre-surround receptive fields of retinal ganglion cells to the LoG operator and Marr’s zero-crossings, and from there to the Canny detector – one of the most popular algorithms in image processing. The key ideas are worth summarising:

• Edges are zero-crossings of the second derivative of image intensity, a principle Marr derived from first principles and validated against neurophysiology
• The LoG operator implements this computationally: a Gaussian pre-filter controls scale and suppresses noise, whereas the Laplacian finds sign changes
• Canny refines the idea with NMS for thin, well-localised edges, and hysteresis thresholding to preserve continuous contours without fragmenting them

Edge detection may seem like a solved problem in an era of end-to-end learned vision systems, but it remains the conceptual foundation of a surprisingly wide range of techniques. Some practical applications worth exploring on your own include:

• Hough transform for line and circle detection – it operates directly on edge maps
• Contour-based object detection – a classical pre-deep-learning approach that is still competitive in constrained domains
• Medical image segmentation – where edge-based pre-processing still complements learned models for thin-structure detection

That brings us to a close – thanks for reading, I hope this post was insightful and entertaining. Stay curious!

References

1. Vaswani, A., Shazeer, N., Parmar, N., Uszkoreit, J., Jones, L., Gomez, A. N., Kaiser, Ł., & Polosukhin, I. (2017). Attention Is All You Need. arXiv:1706.03762. ︎
2. Vision Banana team, Google DeepMind (2026). Image Generators are Generalist Vision Learners. arXiv:2604.20329. ︎
3. Robicheaux, P. et al. (2025). RF-DETR: Neural Architecture Search for Real-Time Detection Transformers. arXiv:2511.09554 (ICLR 2026). ︎
4. Kudithipudi, D. et al. (2025). Neuromorphic computing at scale. Nature 637, 801–812. ︎
5. Ashtiani, F., Idjadi, M. H., & Kim, K. (2026). Integrated photonic neural network with on-chip backpropagation training. Nature 651, 927–932. ︎
6. Assran, M. et al. (2023). Self-Supervised Learning from Images with a Joint-Embedding Predictive Architecture (I-JEPA). arXiv:2301.08243. ︎
7. Marr, D. (1982). Vision: A Computational Investigation into the Human Representation and Processing of Visual Information. W. H. Freeman; reissued by MIT Press (2010). ︎
8. Darwin, C. (1859). On the Origin of Species by Means of Natural Selection. John Murray, London. ︎
9. Thompson, D’A. W. (1917). On Growth and Form. Cambridge University Press. ︎
10. Julesz, B. (1971). Foundations of Cyclopean Perception. University of Chicago Press. ︎
11. Marr, D., & Hildreth, E. (1980). Theory of edge detection. Proceedings of the Royal Society of London B, 207(1167), 187–217. ︎
12. Canny, J. (1986). A Computational Approach to Edge Detection. IEEE Transactions on Pattern Analysis and Machine Intelligence, PAMI-8(6), 679–698. ︎

Read it in: Español.

We are excited to introduce the new team of mentors for the rOpenSci 2026 Champions Program! This year we have eleven individuals committed to open science, bringing together a rich diversity of backgrounds and perspectives. The team is made up of people joining the program for the first time, former Champions returning as mentors, and experienced mentors from previous cohorts returning to continue to strengthen this community.

This year’s mentors come from a variety of disciplines and countries, and are active voices in the R community in Latin America and beyond. With their guidance, the new group of Champions will not only develop their projects, but also grow as leaders in open science and research software development.

New mentors

Alber Hamersson Sánchez Ipia

Hi! I’m Alber and I’m going to be an rOpenSci mentor this year.

I was born in Colombia, in the department of Cauca, in one of the country’s most mountainous regions, called Tierradentro. I am a Cadastral and Geodetic Engineer from the Francisco José de Caldas District University in Colombia, where I earned a Master’s in Information and Communication Sciences; additionally, I completed another Master’s degree in Geoinformatics at the University of Münster in Germany, and later earned a PhD in Earth System Science at the National Institute for Space Research (INPE) in Brazil. Currently, I live and work in Brazil and serve as a research assistant at the same INPE.

Part of my daily work involves writing R code to process spatial data and ensure the reproducibility of the scientific experiment results, so I am familiar with R package development. Additionally, I am a co-author of the segmetric package, which is currently available on CRAN, and I maintain one of the Data Carpentry lessons, specifically the introduction to R for geospatial data.

I am interested in sharing the knowledge and experience I have accumulated so far with anyone who is going to write scientific or statistical software, particularly in Spanish. For this reason I am joining rOpenSci, where I hope to be part of and help build a community of developers.

Pablo Paccioretti

Hello! I am Pablo, Agricultural Engineer and PhD from the National University of Córdoba (Argentina), where I work as a researcher and teacher. Since my student years I have been interested in Statistics, which directed my work towards data analysis. In particular, I apply and develop methodologies and software tools to analyze georeferenced data from field trials and agricultural monitoring platforms.

I am interested in the development of open tools for data processing and analysis. I have developed scientific software, including R packages for georeferenced data analysis.

My participation in the Champions Program arises from an interest in strengthening the links between applied data analysis and programming, and promoting good practices in both areas. Through this program I hope to contribute to the community by sharing experiences and resources, while also learning from other professionals working in different contexts and disciplines.

Champions to mentors

Erick Navarro Delgado

Hello! My name is Erick Navarro, and I have a degree in biology from the Universidad Nacional Autónoma de México, and I am a PhD candidate in Bioinformatics at The University of British Columbia. I was born and raised in Mexico City, but currently live in Vancouver, Canada. My line of research is focused on developing computational tools to understand how genetic factors and environmental exposures/lived experiences act together or separately to shape our molecular landscape.

I am excited to participate in the rOpenSci Champions Program because I believe that open and accessible science is essential for conducting relevant research whose results benefit everyone in our society. In this program I hope to connect with new members of the open science community, share my programming skills, and drive software development in Latin America.

Guadalupe Pascal

Hello! My name is Guadalupe.

I am a researcher and project coordinator in optimization and data science for decision making in social systems, with transfeminist, open science and regional perspectives. I am also an associate professor of optimization and quantitative methods (UNLZ-UCA) and professor in data science and artificial intelligence courses (UGR). I have a Master’s in Decision Systems Engineering from URJC (Spain) and am an industrial engineer from UNLZ (Argentina), a PhD student in Information Technology and Engineering (URJC-UNLZ), and hold deplomas in Gender and Society (UNLZ), Cognitive Neuroscience (Neurotransmitting), and Education in the Age of Artificial Intelligence (UMET). I am part of the Matilda Latin American Open Chair and Women in Engineering, as a founding member and representative of the Gender Network of Engineering Faculties of Argentina.

I am also currently part of the rOpenSci community as a 2025-2026 cohort Champion, and I am very excited to be a mentor in this program for several reasons. On the one hand, I have a deep gratification from being engaged with the current cohort. From a simple point of view, the quality and rigor with which the program is implemented in all its instances have a direct impact on the quality and rigor of my own work. And from an holistic point of view, this serves as extremely valuable and compelling evidence of the synergy within communities of practice in developing skills and producing situated knowledge: the rOpenSci Champions Program is a concrete and real example of how communities share knowledge and, fundamentally, values, perspectives and embodied learning. On the other hand, I am looking forward to the challenge of being a mentor in this program because, although it is a role that I have played in other environments, I have never mentored the development of someone else’s R package. Finally, I would like to work in this role to share my experiences as both a mentor and mentee with the community. I believe that accompanying each other in a formative and transformative process is one of the most human dimensions of this ecosystem in which we work.

Andrea Gomez Vargas

I am Colombian by origin and Argentinean by choice, where today I live, develop my career, and actively participate in the community of R. I am a sociologist and work in the national statistics office of Argentina, in the area of social statistics, where I analyze information about the population to understand inequalities and living conditions.

The R community is my favorite space to share knowledge and build collectively. Currently, I am co-organizer of R in Buenos Aires and RSE Argentina, and I also participate in communities such as R-Ladies, LatinR and rOpenSci, contributing to the strengthening of networks at local, regional and global levels, promoting the learning and use of open tools in data science.

I was a Champion in the 2023-2024 cohort. where I developed {ARcenso} a package that facilitates access to historical census data for Argentina. I am motivated to continue in the program as a mentor to continue promoting open knowledge and to accompany other people in the development of projects with an impact on their communities.

Monika Avila Marquez

Hi, I am Monika, a postdoctoral researcher in statistics at the University of Geneva, where I work on causal inference and machine learning methods for panel data. I have a PhD in econometrics and my research focuses on the development of semi-parametric estimators that combine machine learning techniques with econometric foundations for estimating panel data models. I also work on mixed effects model selection and causal inference with interference.

I am co-organizer of the R-Ladies Geneva chapter, where I strive to build an inclusive community of practice for people using R in research.

I am participating as a mentor in this program because I want to give back for all that rOpenSci has given me. This community has accompanied me in my professional development – as a source of resources, as a learning space and as an example of what it means to do open science with rigor and generosity. Today I have the opportunity to offer that same support to others, and that excites me deeply.

Returning mentors

Luis D. Verde Arregoitia

Hi, I’m Luis D. Verde Arregoitia, a Mexican living in Xalapa, Mexico. Biologist and PhD in Biological Sciences, I am a mammal specialist with experience in R programming for data analysis, visualization and statistical modeling. I am also a certified instructor and author of several packages.

I was a mentor in two previous cohorts of the program where I have supported software developers in Latin America and I return with much enthusiasm to this new cohort.

Pao Corrales

I have a PhD in Atmospheric Sciences from the University of Buenos Aires (Argentina) and am currently working in Australia at the 21st Century Weather Centre as a research software engineer.

I actively participate in R-Ladies, R Forwards, The Carpentries, LatinR and rOpenSci, learning and sharing knowledge about R in the community. In 2023 I participated in the Champions Program as a Champion, submitting the agroclimate package to the rOpenSci peer review process. I learned a lot and connected with people from all over the world. Tt was an excellent experience!

I am passionate about teaching and helping other people grow in what they do, access new opportunities and develop professionally and as individuals. I am very excited to participate again this year as a mentor in the Latin America Champions Program.

Francisco Cardozo

My name is Francisco Cardozo. I am originally from Colombia and came to the United States to pursue my doctoral studies. I am currently working at the University of Miami as a postdoctoral researcher in the IMPAC research center, an institution dedicated to advancing our understanding of adolescent development. I have participated in the Champions Program on several occasions. Much of my professional work has focused on research design and the application of statistical methods, particularly through the use of the R software environment.

Milagros Mendoza

Hello, my name is Milagros. I am an ecologist and statistician driven by a desire to understand the complex systems that intertwine nature, society, and data. Throughout my career, I have worked with interdisciplinary data in the fields of climate, demography, and ecology, always striving to translate that data into knowledge that engages with reality and contributes to more informed decision-making. I am currently pursuing a postdoctoral fellowship at the Vale Institute of Technology in Brazil, where I am part of the research group on territories and natural resources.

I decided to serve as a mentor at rOpenSci because I am motivated to help more people develop confidence in using scientific tools, strengthen their critical thinking, and actively engage within the academic community. In this sense, I view mentoring as a learning space focused on dialogue and mutual growth.

Elio Campitelli

I am from Argentina but two years ago I moved to Australia because it is the only other country that starts with A and uses the same type of plug.

I am doing a postdoc at Monash University researching interactions between Antarctic sea ice and the atmosphere.

I have been a mentor to previous cohorts of the program. It was a great experience that I want to repeat once more.

What’s next

We are happy to have this diverse and talented team of mentors, who embody the values of collaboration and commitment to collective growth. Their support will be key to helping the new Champions move their ideas and projects forward and contribute to the development of a stronger and more diverse open science community.

The selection of Champions is now complete, and we’ll be announcing them soon.

Introduction

Differential Machine Learning (DML), as introduced in the recent arXiv paper (Differential Machine Learning for 0DTE Options with Stochastic Volatility and Jumps), extends supervised learning by incorporating not only function values but also their derivatives. In financial contexts, this often means sensitivities such as Greeks. However, when direct derivatives are unavailable, we can approximate market dynamics using volatility indicators.

In this project, we adapt DML to Bitcoin price forecasting. Instead of derivatives, we use RSI, MACD, and Bollinger Bands as proxies for volatility. These indicators capture momentum, trend strength, and price dispersion, providing a practical way to embed uncertainty into the learning process. To implement this, we design a twin-network architecture in Keras: one network learns price dynamics from time-based features, while the other learns volatility signals. Finally, we combine them via a stacking ensemble to achieve robust forecasts with confidence intervals.

Why Volatility Variables Instead of Derivatives?

• RSI (Relative Strength Index): Measures momentum and overbought/oversold conditions.
• MACD (Moving Average Convergence Divergence): Captures trend direction and strength.
• Bollinger Bands (upper/lower bands, %B): Quantifies price dispersion and volatility.

These indicators act as empirical substitutes for theoretical derivatives. While DML in its pure form requires sensitivities, in practice, these volatility proxies provide similar information about how prices respond to market forces.

Why Twin Networks?

The idea is to separate the learning tasks:

• The primary network models the continuous component of the price process.
• The auxiliary network models the volatility/jump component. Together, they mimic the decomposition found in stochastic models such as Bates or Heston, but implemented within a flexible neural framework.

Ensemble via Stacking

Once both networks are trained, their predictions are combined using a linear regression meta-model. This stacking ensemble learns the optimal weighting between the primary and auxiliary outputs. The result is a forecast that integrates both trend and volatility signals, significantly improving accuracy compared to either network alone.

Evaluation

• Metrics: RMSE and MAPE, computed with the yardstick package.
• Results:
  • Individual networks → RMSE ~76,000, MAPE ~99%.
  • Stacking ensemble → RMSE ~3,030, MAPE ~3.65.

This demonstrates the power of combining price and volatility signals in a unified framework.

Confidence Intervals

To quantify uncertainty, we compute residual-based confidence intervals around the point forecasts:

$${\widehat{y}}_t\pm 1.96\cdot {\sigma }_{\text{residuals}}$$

This approach uses the standard deviation of training residuals to generate 95% confidence bands. It provides interpretable uncertainty estimates without requiring explicit probabilistic modeling.

Visualization

The forecasts are visualized with ggplot2:

• Grey ribbon → confidence intervals.
• Red line → stacking ensemble forecast.
• Black line → actual BTC prices.

This design clearly communicates both the central forecast and the uncertainty range. The chart you will include at the end of the blog shows exactly this: a red forecast line, black actuals, and a grey confidence band, illustrating how the ensemble integrates volatility information into predictive intervals.

Keras3 in R: Flexible Deep Learning for Financial Forecasting

What is Keras3?

Keras3 is the modern R interface to the Keras deep learning library, built on top of TensorFlow. It allows R users to define, train, and evaluate neural networks with concise syntax while leveraging TensorFlow’s computational power. Unlike earlier versions, Keras3 is fully aligned with TensorFlow 2.x, ensuring long-term support and compatibility.

How We Used Keras3

In our workflow, Keras3 was the backbone for implementing the twin-network architecture:

Why ReLU?

• ReLU (Rectified Linear Unit) is the activation function used in hidden layers.
• Formula: $\text{ReLU}(x)=\max (0,x)$.
• Benefits:
  • Introduces non-linearity, enabling the network to learn complex relationships.
  • Efficient and helps avoid vanishing gradients.
  • Well-suited for financial data where signals can be sparse and directional.

Why Adam?

• Adam (Adaptive Moment Estimation) is the optimizer chosen.
• Combines momentum (using past gradients to accelerate learning) and adaptive learning rates (adjusting step sizes per parameter).
• Benefits:
  • Robust for noisy, non-stationary data like cryptocurrency prices.
  • Requires minimal tuning, making it ideal for plug-and-play workflows.
  • Widely adopted in both academic and applied machine learning.

Contribution to the R Ecosystem

Keras3 bridges the gap between R’s tidyverse/tidymodels ecosystem and modern deep learning:

• Integrates seamlessly with data preprocessing pipelines (recipes, timetk).
• Allows financial analysts and data scientists to stay within R while accessing TensorFlow’s deep learning capabilities.
• Encourages reproducibility: models can be defined, trained, and evaluated entirely in R, without switching to Python.
• Expands R’s role beyond traditional statistical modeling into state-of-the-art AI applications.

Why It Matters for DML

By using Keras3:

• We could separate learning tasks into a primary network (trend/seasonality) and an auxiliary network (volatility/momentum).
• Both networks were trained with ReLU activations and Adam optimization, ensuring stability and efficiency.
• Their outputs were combined in a stacking ensemble, yielding forecasts that integrate both price dynamics and volatility signals.

This demonstrates how Keras3 empowers R users to implement advanced architectures like twin networks, making Differential Machine Learning concepts practical in financial forecasting.

Conclusion

This case study demonstrates how Differential Machine Learning concepts can be adapted for financial forecasting in R:

• Volatility indicators serve as practical substitutes for derivatives.
• Twin-network architecture in Keras captures both trend and volatility.
• Stacking ensembles significantly improves predictive performance.
• Residual-based confidence intervals provide interpretable uncertainty estimates.

By combining academic ideas with reproducible R workflows, we can build robust forecasting pipelines that bridge theory and practice.

I tend to write a lot of functions that create specific graphics implemented with ggplot2. Although I try to pick graphic parameters (e.g. colors, text size, etc.) that are reasonable, I will typically define all relevant aesthetics as parameters to my function. As a result, my functions tend to have a lot of parameters. When I need to debug the function I need to have all those parameters set in the global environment which usually requires me highlighting each assignment and running it. This function automates this process. You can pass any function and it will attempt to set parameters to the given environment (the global environment by default). It will return a data frame with a column indicating if the variable was set and the value. This is useful to know what parameters don’t have a default value that need to be set yourself.

#' Set function parameters to an environment.
#'
#' This function is designed to help debug functions. It will attempt to set all
#' the default parameter values to the specified environment (global environment
#' by default). This is useful for when you want to execute code within the 
#' function definition interactively but need the parameters set in the current 
#' environment.
#'
#' **Warning:** This function will modify the global environment and therefore 
#' violates CRAN policy
#' ["Packages should not modify the global environment (user’s workspace)"]
#' (https://cran.r-project.org/web/packages/policies.html#Source-packages).
#'
#' @param FUN the function to assign parameters to an environment.
#' @param envir the environment to assign the variables to. Defaults to the 
#'        global environment.
#' @param verbose whether to return the data frame invisibly or to print the results.
#' @return a data frame where row names correspond to the parameter name with 
#'        two columns: `set` which is logical indicating if the variable was set 
#'        and `value` with a character representation of the variable value.
set_function_params <- function(FUN, envir = globalenv(), verbose = interactive()) {
    params <- formals(FUN)
    params_set <- data.frame(row.names = names(params),
                             set = rep(FALSE, length(params)),
                             value = rep(NA_character_, length(params)))
    for(param in names(params)) {
        value <- params[[param]]
        if(!missing(value)) {
            if(is.character(value)) {
                assign(param, value, envir = envir)
                params_set[param,]$value <- value
            } else {
                assign(param, eval(value), envir = envir)
                params_set[param,]$value <- eval(value)
            }
            params_set[param,]$set <- TRUE
        }
    }
    if(verbose) {
        return(params_set)
    } else {
        invisible(params_set)
    }
}

Very recently I was trying to debug a function that creates profile plots for cluster analysis (clav::profile_plot(), documentation). This function has 23 parameters! Setting these all manually is pretty tedious.

# List objects in the current environment
ls()

[1] "set_function_params"

# Call the function
param_set_result <- set_function_params(clav::profile_plot)

# Check to see if the parameters are actually set
ls()

 [1] "bonferroni"          "center_alpha"        "center_band"        
 [4] "center_fill"         "cluster_label_hjust" "color_palette"      
 [7] "hjust"               "label_clusters"      "label_means"        
[10] "label_outcome_means" "label_profile_means" "param_set_result"   
[13] "point_size"          "se_factor"           "set_function_params"
[16] "standardize"         "text_size"           "title"              
[19] "ylab"               

We can examine the data frame which gives a summary of the parameters set (or not).

param_set_result

                      set               value
df                  FALSE                <NA>
clusters            FALSE                <NA>
df_dep              FALSE                <NA>
standardize          TRUE                TRUE
bonferroni           TRUE                TRUE
label_means          TRUE                TRUE
label_profile_means  TRUE                TRUE
label_outcome_means  TRUE                TRUE
center_band          TRUE                0.25
center_fill          TRUE             #f0f9e8
center_alpha         TRUE                 0.1
text_size            TRUE                   4
hjust                TRUE                 0.5
point_size           TRUE                   2
se_factor            TRUE                1.96
color_palette        TRUE                   2
cluster_labels      FALSE                <NA>
cluster_order       FALSE                <NA>
label_clusters       TRUE                TRUE
cluster_label_x     FALSE                <NA>
cluster_label_hjust  TRUE                   5
ylab                 TRUE Mean Standard Score
title                TRUE    Cluster Profiles

JAGS 5.0.0-beta is now available from SourceForge.

The beta release is for two groups of people:

• People who have written software depending on JAGS, in particular authors of R packages that depend on one of the four interfaces between R and JAGS – rjags, runjags, R2jags, and jagsUI. Currently some of these packages do not pass the CRAN tests with the new version of JAGS. Some time to fix these problems before the official release is helpful.
• Anyone who wants to try out the new version and find problems with it before the official release.

Please send feedback via the JAGS forums or file a bug report

The JAGS library

The following packages are available:

• Source tarball
• Windows binary installer (x86_64)
• macOS binary installer
  • There is a single macOS installer for both x86_64 and arm64.

The rjags package

In order to interface to JAGS 5.0.0 from R you will need rjags_5-1. This is not yet available from CRAN because some of the reverse dependencies do not yet work with version 5.0.0 of JAGS. The following packages are provided:

• Source tarball
• Windows binary (x86_64)
• macOS binaries
  • Separate binaries are provided for x86_64 and arm64 and for R version 4.5.3 and 4.6.0.

I’m getting more and more into data engineering these days and having used R for a long time, I’m seeing a lot of problems that look nail-shaped to my R-shaped hammer. The available tools to solve those problems exist for (presumably) very good reasons, so I wanted to take some time to dig into how to use them and compare their workflows to what I would otherwise naively do in R.

I should mention here that I’m currently open to data/code-related opportunities and am actively seeking a new role – if your organisation is looking for someone aligned with my skillset, please get in touch with me any way you can, e.g. contact at jcarroll.com.au.

I’m a firm believer in “you learn with your hands, not with your eyes” so I wanted to actually build something. I definitely could spin up Claude Code and have it produce the entire thing for me – and in a different project I might do that – but in this case I want to make the mistakes myself so I can learn where the complexity really lives and where my prior assumptions are misaligned. I did have Claude (the chat version, not the full coding agent) guide me through the steps to get this project running, and I did let it clean up my SQL; this project wasn’t about learning to better optimise my SQL, but understanding exactly what it produced will help me write a better version on my next iteration.

Thinking of a real-world project I could take for a spin, I decided to build some ingestion for my personal finances. I’ve used Quickbooks previously which connects up to my bank and helps categorise personal and business (as a freelance contractor) expenses. I decided I’ll build my own ‘slowbooks’ processing workflow based on some manual exports (I don’t think my bank has an API).

Both of the approaches I’ll compare here build on the idea of a Makefile which connects up commands to run based on dependencies, and only runs what is needed; if all the input dependencies of a step have not changed, there’s no need to re-run that step. From what I understand, you could largely get away with just writing some Makefiles (or the newer implementation just) but these two approaches help to better structure how that’s constructed.

This is a somewhat longer post than some of mine, so here’s some quick links to the sections

dbt

One tool that comes up frequently is ‘data build tool’ most commonly referred to as just dbt, though that full name doesn’t even show up on their website. Started in 2016, it’s released as a Python package (dbt-core) though if you do try to just install something called ‘dbt’ you get the cloud CLI tool which isn’t quite the same. Naming stuff is hard.

It’s a way to write code you can commit, which translates to SQL and performs data ingestion, processing, transformation, and storage in a structured way with relationships between various steps in the workflow. It adds macros on top of plain SQL to make the transformations easier, written in jinja, a template engine which enables writing something more like Python within SQL.

This episode of Data Science Lab from Posit walks through an example of using dbt, and while it’s a fantastic overview of what a project looks like, it can’t answer all of the ‘how would I do that?’ problems that will come up in a different project.

Like they did, I will use DuckDB for a database – I enjoyed reading through ‘DuckDB in Action’ with the DSLC.io book club and can definitely see the advantages over SQLite which I would previously have reached for in this case.

I installed dbt via uv – the official instructions use pip and I’ve been burned too many times with that tool; uv is much nicer. Nonetheless, I still encountered Python-related issues because it looks like dbt doesn’t yet support Python 3.14 and yet this isn’t mentioned in their instructions either. I got it working with this command, adding the dbt-duckdb extension I plan to use, as well as streamlit to make a dashboard later

uv init slowbooks --python 3.12
cd slowbooks
uv add dbt-duckdb duckdb streamlit

Adding a profiles.yml in the project root defining the database (DuckDB) I want to produce to store the tables

slowbooks:
  target: dev
  outputs:
    dev:
      type: duckdb
      path: slowbooks.duckdb
      schema: main

I can then initialise the project with

uv run dbt init . --skip-profile-setup

This creates the basic project structure, and there’s a lot going on.

I also needed to define a dependency in packages.yml so that I could use the macros

packages:
  - package: dbt-labs/dbt_utils
    version: [">=1.0.0"]

and ran

uv run dbt deps

I put my exported CSVs (several for my transaction/savings accounts and one for my credit card) in a new raw/ folder; my understanding is that the seeds/ folder is for static data, although that’s the folder used in the Posit tutorial above.

I also ran some pre-processing over my CSVs to categorise the merchants. My bank provides a ‘category’ and ‘subcategory’ for each item, but I wanted to be able to override some of those to more specific definitions so that I could group by them, e.g. ’total spent on books’ since I mainly buy those from just a couple of merchants. This produced a new CSV of patterns, resolved names, and classifications, since the ‘description’ of an item in my transactions might have, e.g.

Paypal *FruitShop 0401000000 Au

and I want to identify the ‘FruitShop’ part, so I can match against that pattern. This is a (fairly) static file (the source data will occasionally be extended), so that did go into seeds/.

{targets}

The whole time I’ve been learning about dbt I’ve had a voice in my head asking “can’t I just use {targets}?” Yes, it’s an R-specific tool, but it does a fantastic job at what it does. It’s not a new tool at all – this post from 2021 demonstrates the power of it, and Miles McBain has been singing the praises of it since at least as early as 2020 (along with the predecessor {drake}).

Rather than double up all of my inputs, I will just keep the {targets} implementation as a subdirectory of my dbt project and refer to the exact same source files. I will create a distinct database, though.

Installing {targets}, provided you already have a working R installation, is as straightforward as

install.packages("targets")

within an R session, be that in RStudio, Positron, Emacs, or a terminal.

As for the rest of the file structure, 100% of the R code here goes into a _targets.R file – much cleaner, albeit that’s a tradeoff in terms of separating different components.

Table of Contents

Comparing Workflows

For the actual processing I’m going to show both dbt and {targets} approaches in tabsets for switching back and forth.

For dbt a ‘model’ is a select statement producing a table, with the structure being models split out into three layers of increasingly production-ready data. From the dbt docs, these are defined as:

• Staging: Preparing atomic building blocks
• Intermediate: Purpose-built transformation steps
• Marts: Business-defined entities

and I’m trying to stick to that as best as I can.

Staging – Load Data

The first step was to ingest that into a ‘staging’ model. This is where the initial data loading happens. For this personal project I’ve exported the CSV files I need, and will do so again in the future, adding them to the same folder for de-duplication within the pipeline. In a more mature project these might be read from an API or a connection to a managed database, and both approaches can easily switch between different ’environments’ (dev, staging, prod, …) without adjusting much, certainly without having to rename all the dependency labels.

• 

  First I defined the sources in models/staging/sources.yml, leveraging DuckDB’s read_csv() to read all the CSV files in raw/

  version: 2
  
  sources:
    - name: all_raw
      schema: main
      meta:
        external_location: "read_csv('raw/*.csv', filename=true, union_by_name=true, header=true)"
      tables:
        - name: all_transactions
          description: "All raw CSV exports"
  

  Reading in these files occurs in a model models/staging/stg_bank.sql and models/staging/stg_cc.sql, the first of which is

  with source as (
      select * from {{ source('all_raw', 'all_transactions') }}
  ),
  
  filtered as (
      select * from source
      where filename not like '%visa_%'
  ),
  
  cleaned as (
      select
          -- Parse YYYYMMDD integer date format
          strptime(cast("Date" as varchar), '%Y%m%d')::date as date,
  
          -- Collapse whitespace runs in description
          regexp_replace(trim("Description"), '\s+', ' ', 'g') as description,
  
          -- Debit = spend (positive), Credit = refund/income (negative)
          coalesce("Debit", 0) - coalesce("Credit", 0) as amount_aud,
  
          "Category"    as raw_category,
          "SubCategory" as raw_subcategory,
          filename      as raw_source
      from filtered
      where "Date" is not null
  )
  
  select * from cleaned
  

  With a similar model for the credit card data. I’ve separated these based on matching to 'visa' in the filename, and I’ve named my files according to this.

  These are then combined in models/staging/stg_transactions.sql, referencing each of the dependencies with the ref() macro. A ‘surrogate key’ is created to uniquely identify rows, so that when I add more data, the duplicates will drop out. This does mean that any intentional duplicates: double records on the same date from the same merchant for the same amount – e.g. buying one ice-cream, dropping it, buying another – will also be dropped, but I’m considering that an edge-case and not worrying about it.

  with bank as (
      select * from {{ ref('stg_bank') }}
  ),
  
  cc as (
      select * from {{ ref('stg_cc') }}
  ),
  
  unioned as (
      select * from bank
      union all
      select * from cc
  ),
  
  with_surrogate_key as (
      select
          {{ dbt_utils.generate_surrogate_key(['date', 'description', 'amount_aud']) }} as transaction_id,
          date,
          description,
          amount_aud,
          raw_category,
          raw_subcategory,
          raw_source
      from unioned
      where "Description" not ilike '%Internet Withdrawal%' -- drop transfers between accounts
        -- this does include manual payments, but most of these are small
  ) 
  
  select * from with_surrogate_key
  

  I’ve also stripped out the ‘internet withdrawal’ records as these are mostly transfers between my own accounts. It also includes manual transfers to e.g. contractors or even some bills, but dealing with these didn’t seem worth the effort.

  One point worth noting here is that this processing is all in SQL; I definitely got the feeling after working with this tool that it was made for data folks who naturally reach for SQL when working with data. Personally, I prefer an abstraction on top of my SQL, so this felt limiting to me, but tastes will absolutely differ.

  The merchants seed file is automatically loaded with the name seed_merchants matching the file name.

• 

  The equivalent in {targets} uses tar_target() to identify dependencies and things to be output. I start by identifying the files I want to read in. A strict comparison would have been to do another grepv() with invert=TRUE but setdiff() works nicely here

  RAW_DIR <- "../raw"
  
  file_list <- list.files(RAW_DIR, full.names = TRUE)
  cc_list <- grepv("visa.*\\.csv$", file_list)
  bank_list <- setdiff(file_list, cc_list)
  

  (sidenote: ooh, yeah – I get to use that new grepv() added in R 4.5.0)

  Loading the data requires a function, but now we can leverage R and its abstractions, in this case {dplyr}, {stringr}, and {lubridate}

  stage_source <- function(files) {
    read_files(files) |>
      mutate(
        date = ymd(as.character(Date)),
        # Collapse whitespace runs — mirrors regexp_replace(trim(Description), '\s+', ' ', 'g')
        description = str_squish(Description),
        # Debit = spend (positive), Credit = refund/income (negative)
        amount_aud = coalesce(Debit, 0) - coalesce(Credit, 0),
        raw_category = Category,
        raw_subcategory = SubCategory,
        raw_source = filename
      ) |>
      filter(!is.na(Date)) |>
      select(
        date,
        description,
        amount_aud,
        raw_category,
        raw_subcategory,
        raw_source
      )
  }
  

  and combining the sources along with filtering out the transfers

  stg_transactions <- function(bank, cc) {
    bind_rows(bank, cc) |>
      # Drop inter-account transfers — mirrors WHERE description NOT ILIKE '%Internet Withdrawal%'
      filter(!str_detect(str_to_lower(description), "internet withdrawal")) |>
      surrogate_key(c("date", "description", "amount_aud")) |>
      select(
        transaction_id,
        date,
        description,
        amount_aud,
        raw_category,
        raw_subcategory,
        raw_source
      )
  }
  

  The surrogate_key function is something I did have to define, but Claude happily provided me with an equivalent to what’s in dbt

  surrogate_key <- function(df, cols) {
    df |>
      mutate(
        transaction_id = purrr::pmap_chr(pick(all_of(cols)), \(...) {
          vals <- list(...)
          parts <- purrr::map_chr(seq_along(vals), \(i) {
            v <- vals[[i]]
            if (is.na(v)) "^^NULL^^" else as.character(v)
          })
          digest::digest(
            paste(parts, collapse = "|"),
            algo = "md5",
            serialize = FALSE
          )
        })
      )
  }
  

  With those pieces, plus loading the merchants file, the full pipeline so far is

  list(
    tar_target(cc_files,   cc_list,   format = "file"),
    tar_target(bank_files, bank_list, format = "file"),
  
    # Staging
    tar_target(stg_bank, stage_source(bank_files)),
    tar_target(stg_cc, stage_source(cc_files)),
    tar_target(stg_txns, stg_transactions(stg_bank, stg_cc))
  )
  

Top of this section | Table of Contents

Intermediate - Joins and Enrichment

For this ‘simple’ example there won’t be a lot of difference between an ‘intermediate’ stage and a final ‘mart’ stage, but this is where the merging with the merchant categories occurs. The transactions from staging are loaded and joined according to the patterns I’ve defined in the seed file.

• 

  In this step the transaction items are categorised according to the merchant file, taking the entry with the most specific pattern match. In models/intermediate/int_transactions_categorised.sql

  with transactions as (
      select * from {{ ref('stg_transactions') }}
  ),
  
  merchants as (
      select * from {{ ref('seed_merchants') }}
  ),
  
  matched as (
      select
          t.*,
          m.keyword,
          m.merchant_name,
          m.merchant_category
      from transactions t
      left join merchants m
          on t.description ilike '%' || m.keyword || '%'
  
      -- Where multiple keywords match, take the longest (most specific)
      qualify row_number() over (
          partition by t.transaction_id
          order by length(m.keyword) desc
      ) = 1
  )
  
  select
      transaction_id,
      date,
      description,
      amount_aud,
      raw_category,
      raw_subcategory,
      raw_source,
      coalesce(merchant_name,     'Unknown')       as merchant_name,
      coalesce(merchant_category, 'Uncategorised') as merchant_category
  from matched
  

  From there, a monthly aggregation table is produced. The date_trunc() feature makes this fairly clean, and being able to sum() values is nice. In models/intermediate/int_monthly_balances.sql

  select
      date_trunc('month', date)::date as month,
      sum(amount_aud) as total_spend_aud,
      count(*) as transaction_count
  from {{ ref('int_transactions_categorised') }}
  group by 1
  

• 

  Rather than relying on the implementation of SQL, R can use the {fuzzyjoin} package to perform the matching

  categorise_transactions <- function(transactions, merchants) {
    # Mirrors: JOIN merchants ON description ILIKE '%' || keyword || '%'
    matched <- fuzzyjoin::fuzzy_left_join(
      transactions,
      merchants,
      by = c("description" = "keyword"),
      match_fun = \(x, y) {
        str_detect(str_to_lower(x), str_to_lower(y), negate = FALSE)
      }
    ) |>
      # Where multiple keywords match, prefer the longest — mirrors QUALIFY ROW_NUMBER() OVER (... ORDER BY length(keyword) DESC) = 1
      group_by(transaction_id) |>
      arrange(desc(str_length(keyword)), .by_group = TRUE) |>
      slice(1) |>
      ungroup() |>
      mutate(
        merchant_name = coalesce(merchant_name, "Unknown"),
        merchant_category = coalesce(merchant_category, "Uncategorised")
      ) |>
      select(
        transaction_id,
        date,
        description,
        amount_aud,
        merchant_name,
        merchant_category,
        raw_category,
        raw_subcategory,
        raw_source
      )
  
    matched
  }
  

  The month aggregations are a bread-and-butter problem for {dplyr}

  monthly_balances <- function(transactions_categorised) {
    transactions_categorised |>
      mutate(month = floor_date(date, "month")) |>
      group_by(month) |>
      summarise(
        total_spend_aud = sum(amount_aud),
        transaction_count = n(),
        .groups = "drop"
      )
    }
  

  and now the pipeline can include those steps

  list(
    tar_target(cc_files,   cc_list,   format = "file"),
    tar_target(bank_files, bank_list, format = "file"),
  
    tar_target(merchant_file, "../seeds/seed_merchants.csv", format = "file"),
    tar_target(merchants, readr::read_csv(merchant_file, show_col_types = FALSE)),
  
    # Staging
    tar_target(stg_bank, stage_source(bank_files)),
    tar_target(stg_cc, stage_source(cc_files)),
    tar_target(stg_txns, stg_transactions(stg_bank, stg_cc)),
  
    # Intermediate
    tar_target(int_categorised, categorise_transactions(stg_txns, merchants)),
    tar_target(int_monthly, monthly_balances(int_categorised))
  )
  

Top of this section | Table of Contents

Marts - Summaries and Outputs

I could create some definitive ‘data product’ results here, but for now this is very similar to the ‘intermediate’ stage with one additional grouping by merchant as well as month

• 

  This is basically just a select, but it does filter for uniqueness on the key. In models/marts/mart_transactions.sql

  {{
      config(
          materialized='incremental',
          unique_key='transaction_id'
      )
  }}
  
  select
      transaction_id,
      date,
      description,
      amount_aud,
      merchant_name,
      merchant_category,
      raw_category,
      raw_subcategory,
      raw_source
  from {{ ref('int_transactions_categorised') }}
  
  {% if is_incremental() %}
      where transaction_id not in (select transaction_id from {{ this }})
  {% endif %}
  

  and finally a month/category aggregation in models/marts/mart_category_trends.sql

  select
    date_trunc('month', date)::date as month,
    merchant_category,
    sum(amount_aud)                 as total_aud,
    count(*)                        as transaction_count
  from {{ ref('int_transactions_categorised') }}
  group by 1, 2
  

• 

  These are essentially the same as intermediate, but with an additional dimension for the monthly summary

  mart_transactions <- function(transactions_categorised) {
    # Equivalent to the incremental mart — deduplication by transaction_id
    transactions_categorised |>
      distinct(transaction_id, .keep_all = TRUE)
  }
  
  mart_monthly_summary <- function(mart_txns) {
    mart_txns |>
      mutate(month = floor_date(date, "month")) |>
      group_by(month, merchant_category) |>
      summarise(
        total_spend_aud = sum(amount_aud),
        transaction_count = n(),
        .groups = "drop"
      )
  }
  

Top of this section | Table of Contents

Tests / Validation

One ‘selling point’ I’ve seen for dbt is that it can additionally add validation tests within the workflow. That’s extremely useful to ensure that you’re not producing junk data inadvertently.

• 

  A schema.yml can be added to a model folder with details of tests to be run on the resulting object. This adds a test that the transaction_id is not null and is unique, and specifies the known values of merchant_category column

  version: 2
  
  models:
    - name: mart_transactions
      columns:
        - name: transaction_id
          tests:
            - not_null
            - unique
        - name: merchant_category
          tests:
          - accepted_values:
              name: merchant_category_is_valid
              values:
                - Accommodation
                - Business
                - Cash
                - Clothing
                - Dining & Drinks
                - Donations
                - Education
  [..truncated..]
  

  The tests/ folder can contain additional SQL tests to be run as part of the workflow. These just need to return a result that should be empty if all goes well, with some number of rows returned if they fail the validation - i.e. no news is good news. I will set up more of these as I figure out what else I want to define as my definition of ‘good quality’ but for now I’ll ensure that no records have the ‘Uncategorised’ category, which means I don’t have an entry for them in my seed_merchants.csv definition.

  I’ve set the option for this to 'warn' because while I do want to identify those missing categories, I don’t want it to stop the workflow entirely

  {{ config(severity='warn') }}
  
  select *
  from {{ ref('mart_transactions') }}
  where merchant_category = 'Uncategorised'
  

  In my case there are still some (34) uncategorised transactions (manual transfers), but the merchant_category_is_valid validation passes

  10 of 13 WARN 34 assert_all_transactions_categorised ........................... [WARN 34 in 0.01s]
  11 of 13 START test merchant_category_is_valid ................................. [RUN]
  11 of 13 PASS merchant_category_is_valid ....................................... [PASS in 0.01s]
  

  and otherwise (if I remove one of the ‘valid’ values)

  10 of 13 WARN 34 assert_all_transactions_categorised ........................... [WARN 34 in 0.02s]
  11 of 13 START test merchant_category_is_valid ................................. [RUN]
  11 of 13 FAIL 1 merchant_category_is_valid ..................................... [FAIL 1 in 0.02s]
  

  which in this case shows that one category didn’t match.

  The tests can also be run independently with

  uv run dbt test
  

• 

  {targets} doesn’t have a specific way to test results, but it does have a way to produce artefacts as part of the workflow in exactly the same way as we do for the data, so I can run arbitrary code including a full data validation.

  In my case I’ll use Appsilon’s
  {data.validator} which has an example for {targets} already, but one could just as easily use {pointblank} or {validate}.

  That leverages assertions from {assertr} and is again just another function

  run_tests <- function(mart_txns) {
    report <- data.validator::data_validation_report()
  
    data.validator::validate(mart_txns, name = "mart_transactions") |>
      data.validator::validate_cols(
        predicate = assertr::not_na,
        "transaction_id",
        description = "transaction_id is not null"
      ) |>
      data.validator::validate_cols(
        predicate = assertr::is_uniq,
        "transaction_id",
        description = "transaction_id is unique"
      ) |>
      data.validator::validate_cols(
        predicate = assertr::in_set(valid_categories),
        "merchant_category",
        description = "merchant_category is an accepted value"
      ) |>
      data.validator::add_results(report)
  
    data.validator::save_report(report, output_file = "validation_report.html")
    report
  }
  

  but it ties in nicely because errors here get reported correctly by {targets}.

  Once the workflow has run, the output file can be opened

  browseURL("validation_report.html")
  

  and if all went well it looks like this

  
  Successful {targets} + {data.validator} run

  otherwise a failure is reported

  
  Errored {targets} + {data.validator} run

  and clicking on ‘Show’ opens a table of the offending results.

Top of this section | Table of Contents

Analysis

What’s the point of organising this data if we’re not going to do something with it? This is where I start to really wonder if {targets} maybe has a bigger picture in mind when it connects up the data, because while dbt will do all of the processing in SQL, R will happily continue to do the analysis.

I think this is where a separation of concerns becomes necessary, and that depends on the scale of the data involved. While you or I working on a small project might be very happy to tie the analysis into the data preparation all in one place, Netflix probably wants to segregate the data processing and analysis steps into entirely different divisions, so tying a bow on the cleaned data and letting analysts pick it up from a database makes a lot more sense.

For my example, let’s say I’m interested in analysing which categories have out of the ordinary amounts of spend in a given month - have I spent more on groceries this month? To do that, I want to calculate the average spend in each category each month plus the variation and identify when the spend is more than a standard deviation away from the average.

• 

  Given that this needs the ‘final’ tables, it belongs in the models/marts folder. There is an analysis/ folder in the dbt project by default, but that’s for ad-hoc SQL queries that need to use ref() but don’t necessarily produce anything one wishes to persist.

  Calculating the standard deviations relies on DuckDB’s helpers, which I don’t even want to consider writing in bare SQL myself. In models/mart/mart_category_outliers.sql:

  with monthly as (
      select * from {{ ref('mart_category_trends') }}
  ),
  
  stats as (
      select
          merchant_category,
          avg(total_aud)                                    as mean_spend,
          stddev_samp(total_aud)                            as sd_spend,
          count(*)                                          as n_months
      from monthly
      where total_aud > 0
      group by 1
      having count(*) > 1  -- need >1 observation for stddev
  ),
  
  outliers as (
      select
          m.month,
          m.merchant_category,
          m.total_aud,
          s.mean_spend,
          m.total_aud - s.mean_spend                       as deviation,
          (m.total_aud - s.mean_spend) / s.sd_spend        as z_score
      from monthly m
      inner join stats s using (merchant_category)
      where m.total_aud > 0
        and abs((m.total_aud - s.mean_spend) / s.sd_spend) > 1
  )
  
  select * from outliers
  order by abs(z_score) desc
  

  This creates a new table in the database with the results.

• 

  The sd() function in R is no stranger to anyone who’s done stats, and it drops into this code cleanly because {dplyr} translates SQL, and supports DuckDB

  monthly_outliers <- function(mart_monthly) {
    spend <- mart_monthly |>
      filter(total_spend_aud > 0)
  
    stats <- spend |>
      group_by(merchant_category) |>
      summarise(
        mean_spend = mean(total_spend_aud),
        sd_spend = sd(total_spend_aud),
        n_months = n(),
        .groups = "drop"
      ) |>
      filter(n_months > 1) # need >1 observation for sd
  
    spend |>
      inner_join(stats, by = "merchant_category") |>
      mutate(
        z_score = (total_spend_aud - mean_spend) / sd_spend,
        deviation = total_spend_aud - mean_spend
      ) |>
      filter(abs(z_score) > 1) |>
      arrange(desc(abs(z_score))) |>
      select(
        month,
        merchant_category,
        total_spend_aud,
        mean_spend,
        deviation,
        z_score
      )
  }
  

  which one can examine by asking {dbplyr} to explain how it’s working

  dplyr::copy_to(
    DBI::dbConnect(duckdb::duckdb()), 
    data.frame(x = c(2, 3, 1, 5, 4)), 
    "example"
  ) |> 
    dplyr::summarise(sd_x = sd(x, na.rm = TRUE)) |>
    dplyr::show_query()
  
  ## <SQL>
  ## SELECT STDDEV(x) AS sd_x
  ## FROM example
  

  which shows that it uses the alias STDDEV(x).

Top of this section | Table of Contents

The Complete Workflow

That’s all the pieces I need to push data in the exported CSVs through the pipe and produce a database of monthly aggregated, categorised totals. Here’s how it looks in terms of the two tools.

• 

  The file structure is perhaps best seen in the docs website (see the next section), but essentially the files in models define the workflow

  models
  ├── intermediate
  │   ├── int_monthly_balances.sql
  │   └── int_transactions_categorised.sql
  ├── marts
  │   ├── mart_category_outliers.sql
  │   ├── mart_category_trends.sql
  │   ├── mart_transactions.sql
  │   └── schema.yml
  └── staging
      ├── sources.yml
      ├── stg_bank.sql
      ├── stg_cc.sql
      └── stg_transactions.sql
  

  Throughout the dbt steps I’ve detailed here, each .sql model produced a corresponding table in the DuckDB database. Running

  uv run dbt build
  

  runs through the DAG identifying what needs to be run before what, then runs all the steps in order. An example of this running well looks like

  01:51:09  Found 1 seed, 8 models, 4 data tests, 2 sources, 591 macros
  01:51:09  
  01:51:09  Concurrency: 1 threads (target='dev')
  01:51:09  
  01:51:09  1 of 13 START sql view model main.stg_bank ..................................... [RUN]
  01:51:09  1 of 13 OK created sql view model main.stg_bank ................................ [OK in 0.04s]
  01:51:09  2 of 13 START sql view model main.stg_cc ....................................... [RUN]
  01:51:09  2 of 13 OK created sql view model main.stg_cc .................................. [OK in 0.02s]
  01:51:09  3 of 13 START seed file main.seed_merchants .................................... [RUN]
  01:51:09  3 of 13 OK loaded seed file main.seed_merchants ................................ [INSERT 402 in 0.02s]
  01:51:09  4 of 13 START sql view model main.stg_transactions ............................. [RUN]
  01:51:09  4 of 13 OK created sql view model main.stg_transactions ........................ [OK in 0.04s]
  01:51:09  5 of 13 START sql view model main.int_transactions_categorised ................. [RUN]
  01:51:09  5 of 13 OK created sql view model main.int_transactions_categorised ............ [OK in 0.03s]
  01:51:09  6 of 13 START sql view model main.int_monthly_balances ......................... [RUN]
  01:51:09  6 of 13 OK created sql view model main.int_monthly_balances .................... [OK in 0.03s]
  01:51:09  7 of 13 START sql table model main.mart_category_trends ........................ [RUN]
  01:51:09  7 of 13 OK created sql table model main.mart_category_trends ................... [OK in 0.22s]
  01:51:09  8 of 13 START sql incremental model main.mart_transactions ..................... [RUN]
  01:51:09  8 of 13 OK created sql incremental model main.mart_transactions ................ [OK in 0.26s]
  01:51:09  9 of 13 START sql table model main.mart_category_outliers ...................... [RUN]
  01:51:09  9 of 13 OK created sql table model main.mart_category_outliers ................. [OK in 0.01s]
  01:51:09  10 of 13 START test assert_all_transactions_categorised ........................ [RUN]
  01:51:09  10 of 13 WARN 34 assert_all_transactions_categorised ........................... [WARN 34 in 0.01s]
  01:51:09  11 of 13 START test merchant_category_is_valid ................................. [RUN]
  01:51:09  11 of 13 PASS merchant_category_is_valid ....................................... [PASS in 0.01s]
  01:51:09  12 of 13 START test not_null_mart_transactions_transaction_id .................. [RUN]
  01:51:09  12 of 13 PASS not_null_mart_transactions_transaction_id ........................ [PASS in 0.01s]
  01:51:09  13 of 13 START test unique_mart_transactions_transaction_id .................... [RUN]
  01:51:09  13 of 13 PASS unique_mart_transactions_transaction_id .......................... [PASS in 0.01s]
  01:51:09  
  01:51:09  Finished running 1 incremental model, 1 seed, 2 table models, 4 data tests, 5 view models in 0 hours 0 minutes and 0.78 seconds (0.78s).
  01:51:09  
  01:51:09  Completed with 1 warning:
  01:51:09  
  01:51:09  Warning in test assert_all_transactions_categorised (tests/assert_all_transactions_categorised.sql)
  01:51:09  Got 34 results, configured to warn if != 0
  01:51:09  
  01:51:09    compiled code at target/compiled/slowbooks/tests/assert_all_transactions_categorised.sql
  01:51:09  
  01:51:09  Done. PASS=12 WARN=1 ERROR=0 SKIP=0 NO-OP=0 TOTAL=13
  

  The times on the left side are in UTC, and I can’t find a way to change that, which may be for the best - one can always convert to local after the fact if needed.

  Everything completed with an OK status except for the assertion which I’ve allowed to WARN because I haven’t categorised a handful of records.

  Looking at the resulting database, e.g. in a terminal, shows all the tables which have been created

  duckdb slowbooks.duckdb
  
  DuckDB v1.5.2 (Variegata)
  Enter ".help" for usage hints.
  slowbooks D show tables;
  ┌──────────────────────────────┐
  │             name             │
  │           varchar            │
  ├──────────────────────────────┤
  │ int_monthly_balances         │
  │ int_transactions_categorised │
  │ mart_category_outliers       │
  │ mart_category_trends         │
  │ mart_transactions            │
  │ seed_merchants               │
  │ stg_bank                     │
  │ stg_cc                       │
  │ stg_transactions             │
  └──────────────────────────────┘
  

• 

  The {targets} workflow is defined as a list within the _targets.R file

  list(
    # format="file" means targets re-runs downstream when file list or contents change
    tar_target(
      cc_files,
      cc_list,
      format = "file"
    ),
    tar_target(
      bank_files,
      bank_list,
      format = "file"
    ),
    tar_target(merchant_file, "../seeds/seed_merchants.csv", format = "file"),
  
    # Seeds
    tar_target(merchants, readr::read_csv(merchant_file, show_col_types = FALSE)),
  
    # Staging — passing file vectors directly so targets tracks dependencies correctly
    tar_target(stg_bank, stage_source(bank_files)),
    tar_target(stg_cc, stage_source(cc_files)),
    tar_target(stg_txns, stg_transactions(stg_bank, stg_cc)),
  
    # Intermediate
    tar_target(int_categorised, categorise_transactions(stg_txns, merchants)),
    tar_target(int_monthly, monthly_balances(int_categorised)),
  
    # Tests
    tar_target(validation, run_tests(int_categorised)),
    tar_target(is_violation, validation_violation(validation)),
  
    # Marts
    tar_target(mart_txns, mart_transactions(int_categorised)),
    tar_target(mart_monthly, mart_monthly_summary(mart_txns)),
  
    # Persist to DuckDB
    tar_target(
      db_mart_transactions,
      write_to_duckdb(mart_txns, "mart_transactions")
    ),
    tar_target(
      db_mart_monthly_summary,
      write_to_duckdb(mart_monthly, "mart_monthly_summary")
    ),
  
    # Analysis
    tar_target(outliers, monthly_outliers(mart_monthly))
  )
  

  and running tar_make() (or with the added validation, tar_make_catch()) from the working directory containing that file, the workflow is run

  + bank_files dispatched                      
  ✔ bank_files completed [281ms, 74.42 kB]
  + cc_files dispatched
  ✔ cc_files completed [0ms, 137.82 kB]
  + merchant_file dispatched
  ✔ merchant_file completed [1ms, 15.89 kB]
  + stg_bank dispatched
  ✔ stg_bank completed [88ms, 12.98 kB]
  + stg_cc dispatched
  ✔ stg_cc completed [84ms, 18.60 kB]
  + merchants dispatched
  ✔ merchants completed [137ms, 7.33 kB]
  + stg_txns dispatched
  ✔ stg_txns completed [101ms, 71.10 kB]
  + int_categorised dispatched
  ✔ int_categorised completed [871ms, 86.06 kB]                
  + mart_txns dispatched                                       
  ✔ mart_txns completed [0ms, 86.06 kB]                        
  + validation dispatched                                      
  ✔ validation completed [2.7s, 5.10 kB]                        
  + int_monthly dispatched                                      
  ✔ int_monthly completed [3ms, 650 B]                          
  + mart_monthly dispatched                                     
  ✔ mart_monthly completed [4ms, 2.78 kB]                       
  + db_mart_transactions dispatched                             
  ✔ db_mart_transactions completed [83ms, 70 B]                 
  + is_violation dispatched                                     
  ✔ is_violation completed [1ms, 48 B]                          
  + db_mart_monthly_summary dispatched                          
  ✔ db_mart_monthly_summary completed [79ms, 72 B]              
  + outliers dispatched                                         
  ✔ outliers completed [10ms, 2.30 kB]                          
  ✔ ended pipeline [4.6s, 16 completed, 0 skipped] 
  

  An important difference is that writing to the database only happened in two of the steps at the end, so the (distinct) database only contains those tables

  duckdb targets/slowbooks_r.duckdb
  
  DuckDB v1.5.2 (Variegata)
  Enter ".help" for usage hints.
  slowbooks_r D show tables;
  ┌──────────────────────┐
  │         name         │
  │       varchar        │
  ├──────────────────────┤
  │ mart_monthly_summary │
  │ mart_transactions    │
  └──────────────────────┘
  

  There’s nothing stopping me from also adding a write_to_duckdb() (a function which just opens the database, writes a table, then closes) call for any of the other steps, but I was satisfied that I am building the same thing in both cases.

Top of this section | Table of Contents

DAG / Visualisation / Docs

The similarity to a Makefile of these approaches depends on being able to determine what has ‘changed’ and what is the same, and this is where the two approaches differ. Both consider the workflow as a Directed Acyclic Graph (DAG) with steps taking dependencies on previous steps or data sources. This means I can visualise the workflow as a graph, but also makes for some important differences between how things actually run.

• 

  Every time I run dbt build the entire model is re-run. If the model is ‘incremental’ then it won’t need to do a full CREATE TABLE or re-categorise existing records every time, but all of the steps will be re-run. This also makes for more complexity if I do want to re-categorise existing records, in which case I need to add --full-refresh to the build step.

  With the model defined I can visualise it by generating and serving the documentation site with

  uv run dbt docs generate && uv run dbt docs serve 
  

  This builds a site locally hosted which includes all of the SQL code and a lineage graph showing how different pieces connect together

  
  dbt docs site (click to embiggen)

  Expanding this pane shows more of the DAG for the project, though not all of the connections

  
  dbt DAG for the whole slowbooks project (click to embiggen)

  This is really nice, and shows how the data flows from the raw data to the final summary.

• 

  This is where I think {targets} may have an advantage over dbt – since the workflow considers a hash of the data objects to determine what has changed, even if the code remains the same, it can identify which steps of the DAG are invalidated, and can skip over any steps which don’t need to be re-run.

  This is significant when the data you’re processing is no longer necessarily local to the machine running the pipeline. dbt performs the queries with SQL on the database (in this example the tables are written in DuckDB and materialised as views for downstream models), while the structure I’m using for {targets} here explicitly pulls in the data for R-native processing. I could make it more remote and use lazy tbl() operations via {dbplyr}, but it’s a trade-off one needs to consider.

  A full DAG for the project can be produced in an editor able to render HTML such as RStudio or Positron, with

  targets::tar_visnetwork()
  

  producing an interactive plot

  
  Up-to-date {targets} DAG visualisation (click to embiggen)

  If I change some of the code in the mart definition, and re-evaluate just that function, then re-running targets::tar_visnetwork() shows me which nodes are affected

  
  Invalidated {targets} DAG visualisation (click to embiggen)

  (note the different colour of some nodes). This is fantastic!

  What’s more, if I have a failing test during the validation, I can see what is downstream from that

  
  Failed validation in {targets} DAG (click to embiggen)

  That did require following the article in the {data.validator} docs to define

  validation_violation <- function(report) {
    violations_exist <- report$get_validations() %>%
      summarise(
        sum(num.violations, na.rm = TRUE) > 0
      ) %>%
      pull()
    if (isTRUE(violations_exist)) {
      rlang::abort(
        "Validation schema error",
        body = capture.output(report),
        class = "validation_violation"
      )
    }
    FALSE
  }
  

  and adding to the target

  tar_target(is_violation, validation_violation(validation)),
  

  and instead of running tar_make(), using a tryCatch()

  tar_make_catch <- function() {
      tryCatch(
          tar_make(),
          validation_violation = function(e) {
              print(e)
              tar_visnetwork()
          }
      )
  }
  

  Incredibly powerful stuff, right?

Top of this section | Table of Contents

Exploration

I only built the exploration dashboard as part of the dbt project because I’ve built plenty of shiny apps - I wanted to see what Claude could build based on this database data source. It built a streamlit app which shows the monthly spend broken down by category, and I had it add filters for the various categories, tables of transactions, and the monthly outliers.

The dashboard works great, albeit not perfectly. It looks something like this

There’s obvious issues with it - not least that the legend is incomplete, but for the sort of exploration I wanted to try out, it’s a great starting point.

It reads the summary tables directly from the database, so the analysis doesn’t need to happen within the app – a nice separation of business logic and visualisation.

Comparison

As a final sanity check, I’ll confirm that I get the same number of transactions in the monthly trend tables which are saved to both databases, albeit with different names

duckdb slowbooks.duckdb -c "select sum(transaction_count) from mart_category_trends;"
┌────────────────────────┐
│ sum(transaction_count) │
│         int128         │
├────────────────────────┤
│                   2028 │
└────────────────────────┘

duckdb targets/slowbooks_r.duckdb -c "select sum(transaction_count) from mart_monthly_summary;"
┌────────────────────────┐
│ sum(transaction_count) │
│         int128         │
├────────────────────────┤
│                   2028 │
└────────────────────────┘

As for what I like and don’t like about each approach:

• Language: I don’t reach for SQL as a primary language (it’s absolutely the second language everyone who codes with data should learn, in my opinion) so having to write everything in SQL myself doesn’t appeal so much to me. I’m very happy to be able to use {dplyr} more or less everywhere and have it write the SQL for me. On the other hand, I can see the value in moving to a language that’s closer to the data itself – the abstractions change over time ({dplyr} is notorious for this) so with fewer bells and whistles likely comes more stability. Provided the helper functions are available for things like basic statistics (e.g. in DuckDB) this doesn’t sound like too much of a downside. Doing a bit more research, it seems that dbt does support using Python for models, provided the adapter supports it (which dbt-duckdb does), so that’s a big win for those more familiar with Python, although I am under the impression that not everything works exactly the same for these models.

• Connection: I appreciate the massive leg-up that dbt offers in terms of handling connections to sources via extensions (e.g. dbt-duckdb). I’m sure if you’re not used to that then it looks like magic, but those familiar with working with databases via {dbplyr} and {DBI} it loses some of the wonder. Importantly, the dbt SQL code all runs within the database - downstream models rely on views so the data never really leaves the database. The R version could get closer to this, but I suspect the more common use-case is to actually pull down all of the data in which case it’s likely to fit within RAM.

• Version Control: For people not used to committing their work, again this seems like a huge step up, but R users get taught fairly early on to work with git and track their code, even if it’s just scripts. Someone used to just throwing SQL at a database from a terminal might be rightly amazed at the benefits opened up by tracking code this way, but for me it’s the default state.

• Layout: dbt has so many files for even a simple project that my VSCode file explorer runs off the screen. {targets} has everything in a single file. This could be organised more like dbt with a liberal use of source() calls at the top of _targets.R, say for each model and some utils.

• Interrogation: Perhaps there’s some more tooling I’m not aware of, but the {targets} visualisation of the DAG is a clear winner for me. Part of the tradeoff between ‘run everything locally’ and ‘run everything remotely’ is that I can inspect the intermediate data in the {targets} workflow with tar_read(id) and see what’s happening. I can read the generated table in the database, but for smallish data being able to just crack it open and have a look wins for me.

Other Solutions

While I’ve focused on this comparison between dbt and {targets}, these aren’t the only players in the game. I’m aware of Airflow, at least in the sense that it can ingest dbt pipelines and schedule them. For the Python folks there’s also prefect and dagster, the latter of which also has an R ingestion route in the form of dagsterpipes. A purely R solution is maestro which appears to target (pun intended) data coming from an API or database for which {targets} can’t identify the ‘up-to-date-ness’ (since that involves a hash of the file).

Table of Contents

Conclusion

I’ve vastly grown my understanding of both dbt and {targets} and have a much greater appreciation for what goes into using each of these to move and curate data. Plus, now I have a cool new toy I’ve built to explore my finances. I’m not sharing the code itself – partly so that I don’t risk committing my own finance data by accident, and partly because what I’ve done here isn’t anything you need to build on; if you’re interested in learning either or both of these tools, I recommend you do what I did and build a toy project.

I’m interested to hear what you think of this comparison – have I overlooked some significant difference or similarity? Some use-case where one of them would really shine over the other? Have I misrepresented something? I’m here to learn, so by all means please do let me know. And if you’re looking for someone with a history of programming and data who digs into projects this way, I’m on the market for opportunities.

As always, I can be found on Mastodon and the comment section below.

## ─ Session info ───────────────────────────────────────────────────────────────
##  setting  value
##  version  R version 4.5.3 (2026-03-11)
##  os       macOS Tahoe 26.3.1
##  system   aarch64, darwin20
##  ui       X11
##  language (EN)
##  collate  en_US.UTF-8
##  ctype    en_US.UTF-8
##  tz       Australia/Adelaide
##  date     2026-05-04
##  pandoc   3.6.3 @ /Applications/RStudio.app/Contents/Resources/app/quarto/bin/tools/aarch64/ (via rmarkdown)
##  quarto   1.7.31 @ /usr/local/bin/quarto
## 
## ─ Packages ───────────────────────────────────────────────────────────────────
##  package     * version date (UTC) lib source
##  blob          1.3.0   2026-01-14 [1] CRAN (R 4.5.2)
##  blogdown      1.23    2026-01-18 [1] CRAN (R 4.5.2)
##  bookdown      0.46    2025-12-05 [1] CRAN (R 4.5.2)
##  bslib         0.10.0  2026-01-26 [1] CRAN (R 4.5.2)
##  cachem        1.1.0   2024-05-16 [1] CRAN (R 4.5.0)
##  cli           3.6.5   2025-04-23 [1] CRAN (R 4.5.0)
##  DBI           1.3.0   2026-02-25 [1] CRAN (R 4.5.2)
##  dbplyr        2.5.2   2026-02-13 [1] CRAN (R 4.5.2)
##  devtools      2.4.6   2025-10-03 [1] CRAN (R 4.5.0)
##  digest        0.6.39  2025-11-19 [1] CRAN (R 4.5.2)
##  dplyr         1.2.1   2026-04-03 [1] CRAN (R 4.5.2)
##  duckdb        1.5.2   2026-04-13 [1] CRAN (R 4.5.2)
##  ellipsis      0.3.2   2021-04-29 [1] CRAN (R 4.5.0)
##  evaluate      1.0.5   2025-08-27 [1] CRAN (R 4.5.0)
##  fastmap       1.2.0   2024-05-15 [1] CRAN (R 4.5.0)
##  fs            1.6.7   2026-03-06 [1] CRAN (R 4.5.2)
##  generics      0.1.4   2025-05-09 [1] CRAN (R 4.5.0)
##  glue          1.8.1   2026-04-17 [1] CRAN (R 4.5.2)
##  htmltools     0.5.9   2025-12-04 [1] CRAN (R 4.5.2)
##  jquerylib     0.1.4   2021-04-26 [1] CRAN (R 4.5.0)
##  jsonlite      2.0.0   2025-03-27 [1] CRAN (R 4.5.0)
##  knitr         1.51    2025-12-20 [1] CRAN (R 4.5.2)
##  lifecycle     1.0.5   2026-01-08 [1] CRAN (R 4.5.2)
##  magrittr      2.0.4   2025-09-12 [1] CRAN (R 4.5.0)
##  memoise       2.0.1   2021-11-26 [1] CRAN (R 4.5.0)
##  otel          0.2.0   2025-08-29 [1] CRAN (R 4.5.0)
##  pillar        1.11.1  2025-09-17 [1] CRAN (R 4.5.0)
##  pkgbuild      1.4.8   2025-05-26 [1] CRAN (R 4.5.0)
##  pkgconfig     2.0.3   2019-09-22 [1] CRAN (R 4.5.0)
##  pkgload       1.5.0   2026-02-03 [1] CRAN (R 4.5.2)
##  purrr         1.2.2   2026-04-10 [1] CRAN (R 4.5.2)
##  R6            2.6.1   2025-02-15 [1] CRAN (R 4.5.0)
##  remotes       2.5.0   2024-03-17 [1] CRAN (R 4.5.0)
##  rlang         1.1.7   2026-01-09 [1] CRAN (R 4.5.2)
##  rmarkdown     2.30    2025-09-28 [1] CRAN (R 4.5.0)
##  rstudioapi    0.18.0  2026-01-16 [1] CRAN (R 4.5.2)
##  sass          0.4.10  2025-04-11 [1] CRAN (R 4.5.0)
##  sessioninfo   1.2.3   2025-02-05 [1] CRAN (R 4.5.0)
##  tibble        3.3.1   2026-01-11 [1] CRAN (R 4.5.2)
##  tidyselect    1.2.1   2024-03-11 [1] CRAN (R 4.5.0)
##  usethis       3.2.1   2025-09-06 [1] CRAN (R 4.5.0)
##  vctrs         0.7.1   2026-01-23 [1] CRAN (R 4.5.2)
##  withr         3.0.2   2024-10-28 [1] CRAN (R 4.5.0)
##  xfun          0.56    2026-01-18 [1] CRAN (R 4.5.2)
##  yaml          2.3.12  2025-12-10 [1] CRAN (R 4.5.2)
## 
##  [1] /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library
## 
## ──────────────────────────────────────────────────────────────────────────────

Have you ever looked at a freshly plotted scatter plot and immediately thought, “Ah, this is clearly a logarithmic curve with some heteroskedastic noise,” without running a single line of modeling code? How do you do that? You don’t perform gradient descent in your head. You use your intuition!

As an experienced data scientist, you have seen thousands of datasets in your career. When confronted with new data, your natural neural network (a.k.a. brain) simply draws on this vast library of past mathematical shapes and immediately recognizes the pattern. But what if an artificial neural network could do exactly the same thing? What if it could predict your data without actually being trained on it?

Welcome to the mind-bending world of In-Context Learning (ICL) for tabular data, brought to R via the incredible new TabPFN package (on CRAN).

The Transformer: From Text to Tables

To understand ICL, we have to talk about Large Language Models like ChatGPT (see also Building Your Own Mini-ChatGPT with R: From Markov Chains to Transformers!). When you give a chatbot an unfinished sentence, it doesn’t retrain its weights to guess the next word. It uses a Transformer architecture equipped with an attention mechanism (see also Attention! What lies at the Core of ChatGPT? (Also as a Video!)). It reads the words you provided, understands the dependencies between them (the grammar and context), and instantly extrapolates what comes next.

The genius of TabPFN is taking this exact architecture and applying it to spreadsheets. Instead of a sequence of words, the Transformer reads a sequence of data rows. It treats your features (X) and your target (Y) like the grammar of a language. By comparing all the rows and columns simultaneously in its “context window,” it figures out the dependencies in the table just like a language model figures out dependencies in text.

The model that arises is a foundation model for tabular data, or tabular foundation model for short.

This process is formally known as Few-Shot Learning. You aren’t giving the model an empty brain to train; you are “prompting” a pre-trained brain with a few dozen (or a few hundred) “shots” (rows) of your data to establish the pattern!

The Training Matrix: Learning the Shape of Maths

You might be wondering: If it isn’t training on my data, what exactly was it trained on?

This is where it gets incredibly cool. The researchers who built TabPFN didn’t train it on real-world datasets like housing prices or medical records. Instead, they wrote algorithms to generate millions of completely random, artificially created mathematical dependency structures.

They forced the network to practice on synthetic datasets containing every statistical quirk imaginable: linear trends, severe non-linearities, bizarre interaction effects, extreme missing data mechanisms, and sheer noise. Because it spent its entire training solving billions of abstract maths puzzles, the model learned the fundamental shape of causal mathematical dependencies. When it sees your real-world data, it’s just recognizing a pattern it has already solved synthetically a thousand times before.

Let’s see it in action

Let’s use the venerable iris dataset. Because iris is small and the mathematical boundaries are very clear, it’s the perfect candidate for few-shot learning. Notice how the code looks exactly like traditional machine learning, but under the hood, no training is actually happening!

# Load the package
library(tabpfn)

# 1. Prepare the Data
set.seed(42)
train_indices <- sample(seq_len(nrow(iris)), size = 0.7 * nrow(iris))

iris_train <- iris[train_indices, ]
iris_test  <- iris[-train_indices, ]

# 2. Fit the Model
cat("Generating embeddings...\n")
## Generating embeddings...
tab_fit <- tab_pfn(Species ~ ., data = iris_train)

# 3. Make Predictions
cat("Predicting...\n")
## Predicting...
predictions <- predict(tab_fit, new_data = iris_test)

# 4. Check the accuracy
accuracy <- sum(predictions$.pred_class == iris_test$Species) / nrow(iris_test)
cat("\nSuccess! Overall Accuracy:", round(accuracy * 100, 1), "%\n")
## 
## Success! Overall Accuracy: 97.8 %

When you run this, you will see an accuracy of 97.8%. The model looked at the few examples in iris_train, instantly recognized the multidimensional shapes separating the species using its synthetic intuition, and accurately classified the new test data without a single epoch of traditional backpropagation.

Conclusion

TabPFNis a paradigm shift. For small to medium tabular datasets, we no longer need to spend hours tuning hyperparameters for Random Forests or XGBoost. We can simply hand the data to an experienced, mathematically omniscient Transformer and let In-Context Learning do the heavy lifting.

Give it a try on your own data, and tell us about your experience in the comments below!

Recently I’ve been looking at hourly ridership data from the New York City Subway. Last time we learned that people go to work in the morning and come home in the evening, for example. (All together now: “Only in New York, baby!”) Today we’ll learn that bad weather makes people stay at home. Except, sometimes it doesn’t.

Regular readers will recall that the subway system carries a lot of passengers every day. The ridership data for the whole of 2025 represents just over 1.3 billion entries into the system via an OMNY tap or Metrocard. It’s available aggregated to hourly resolution by station complex. With that data in hand, we can calculate average hourly ridership for every day of the week. This gives us a profile of what, for example, a Monday or a Wednesday typically looks like, by hour. When calculating the average day-of-the week profile we exclude holidays and the like.

Meanwhile, the National Weather Service provides data on severe weather events that affected the New York City region in 2025. We could get more fine-grained if we wanted to, but for now we’ll just use the general list of events the NWS provides. Then we plot the Subway ridership profile for that specific date against the average profile for whatever day of the week the event happened on.

The gray lines are the baseline. The red ones are the bad weather day. The basic shape of the gray lines (and many of the red ones) is set by the rhythm of daily life. The sharp double-peak pattern is what someone I’ve shown too many of these graphs to has taken to calling “The Giant Cat-Ears of Employment”. The cat-ear shapes vary by work day (which might be the topic of another post), but are most sharply-contrasted with the weekends, which look more like little hillocks or drumlins.

We can see a few different cases in the panels. First are days when the weather event put no dent at all in people’s day. This is because of the incredible toughness resilience of New Yorkers, something they are surprisingly very modest about even though there was a weather event in the region that day, it just didn’t impinge on the city much, or at all. The light snow on February 11th or the heavy rain on March 6th are examples here. People just continued to go about their business.

Second are cases where there’s a lot of travel suppression but it’s not really—or not wholly—the weather that’s responsible. The winter storm on Friday December 26th is a case of this. That’s not a regular Friday. Many people are able to stay at home anyway, because it’s the day after Christmas.

Third are cases where the weather does seem to have suppressed travel. These are days like the snow on January 19th, or the shitty weather on Sunday February 16th. These events look like they made people stay at home. Some of these are more severe than others. The strongest example is the flash flooding on Thursday July 31st. That happened in the back half of the day and affected the evening commute directly.

Our fourth and final category is my favorite one. Sometimes snow makes no difference at all, especially if it’s on a workday. Sometimes it’s snowy on the weekend but you’re kind of sick of it, maybe because it’s late in the winter, so you’re either going about your business as usual or you’re just staying indoors. But there’s another kind of snow day.

The weekend of December 13th and 14th 2025 brought the city’s first measurable snow of the year, and in decent amounts, too—between four and eight inches of accumulation. Reports remarked on how long it had been in arriving. The result was that, over the weekend, ridership on the subway went up. Maybe on the Saturday it was to go out and buy the mandatory bread, milk, and eggs.¹ But maybe it was also just to be out in the snow. The next day, the people who didn’t have to go work slept in as usual. But that day, too, across the afternoon, more people than usual headed outside and took the subway somewhere. I’d like to think a bunch of them had a sled under their arm.

A note to myself on survival analysis — KM curves, log-rank tests & Cox models If I wrote it the way I understood it, maybe I’ll actually remember it

Motivations

We see survival analysis or more generally called time-to-event analysis almost all the time when we review journals articles on NEJM etc. Even though we understand the heuristic in interpreting some of the simpler result, I realized that I need to look at this a bit closer to full understand the works and math behind it. There was a recent project that made me feel that my understanding of this is not as competent as I had hoped after talking to one of my statistician colleagues, who also so happen to wrote this blog. Please take a look at Emily’s blog for a better, and a more accurate survival analysis tutorial. This blog is more for my learning so that I can refer back the fundamental when I need a refresher in the future. Also, if I were to write it the way I understood it, maybe that might increase the probably of me recollecting what I understood before. What are we waiting for? Let’s time-to-event this analysis!

Objectives:

= Time-to-event Analysis

• Survival function
• Let’s Calculate By Hand
• Simulation
  • Kaplan-Meier Estimator
  • Cox Proportional Hazard Model
• Acknowledgement
• Oppotunities For Improvement
• Lessons Learnt

Time-to-event Analysis

The name “survival analysis” is a bit misleading if you first encounter it outside of clinical research. The “survival” doesn’t necessarily mean staying alive — it means surviving without experiencing the event. But that even does not necessarily have to be mortality. It could be an unwanted event etc. Outside of clinical research, an event could be the time when a waymo arrives at your doorstep or someone flaked out. Hence time-to-event analysis appears to be more a better terminology, in my opinion.

This is different from the good ol regression is because time is the outcome, not only that it occurred or not (binary), but when! Now then if you’re like me, that’s just negative binomial regression, right? Not quite. Because, there is an additional special feature to time-to-event analysis called censoring.

Censoring can mean that the event did not occur, but it can also mean that we lost track of the patient, or the study ended before the event occurred.

I don’t know about you. But for me, censoring has a negative connotation. It sounds like we’re hiding something. But in survival analysis, censoring is actually a good thing. It means we have partial information about the time to event, even if we don’t know the exact time. So, think of censoring as the last time we noticed that the event DID NOT happen, and it’s usually coded as 0. In good ol regression, we usually will either do a complete case analysis (throw out missing data) or impute. But, imputing outcome is a tad odd, in most cases except here. The common censoring is right censoring, meaning we lose track of someone on the right side of the timeline.

Survival Function

The survival function, written S(t), answers one simple question: “What is the probability that a person has NOT yet experienced the event by time t?”. At t=0, everyone is event-free, so S(0) = 1 (100%). As time goes on, people experience the event, and S(t) decreases.

Let’s Calculate By Hand

Alright, the above looks quite self-explanatory. We have 5 patients, and we are tracking their time to event in months. Patient A experienced the event at 21 months, while patient B was censored at 30 months (we lost track of them). Patient C had the event at 51 months, patient D at 61 months, and patient E was censored at 80 months. Now let’s do some calculation.

Formula: $$ \hat{S}(t) = \prod_{i:, t_i \leq t} \left(1 – \frac{d_i}{n_i}\right) $$

• \(\hat{S}(t)\): The estimated survival function; the probability of surviving beyond time \(t\).
• \(\prod_{i:\, t_i \leq t}\): Product over all event times \(t_i\) that are less than or equal to \(t\).
• \(t_i\): The \(i\)-th observed event (death/failure) time.
• \(d_i\): The number of events (deaths/failures) that occurred at time \(t_i\).
• \(n_i\): The number of individuals at risk (still alive and under observation) just before time \(t_i\).
• \(\frac{d_i}{n_i}\): The estimated probability of the event occurring at time \(t_i\).
• \(1 - \frac{d_i}{n_i}\): The estimated probability of surviving through time \(t_i\).

or simplistically:

\(S(t) = S(t-1).(1-d/n)\)

Let’s calculate by hand:

That’s interesting! I don’t think I’ve calcualte these by hand before and that’s very helpful in just doing a simple example and observing the result. Alright, when we read articles, there typically is a group factor, how do they then use KM to generate 2 survival curves or survival estimate for each group? They just do the same thing but only for the subset of the data that belongs to that group. So, if we have a treatment and control group, we would calculate S(t) separately for each group, and then we can compare the two survival curves to see if there is a difference in survival between the groups. How? We can use the log-rank test to compare the survival curves, or we can use a Cox proportional hazards model to estimate the hazard ratio between the groups. Now, things are starting to look a tad more familiar. Let’s use R and some simple example and see if we can get to log-rank test on just simple KM etimator.

library(tidyverse)

simple_df <- tribble(
  ~time, ~status, ~treatment,
  5,1,1,
  2,1,0,
  6,0,1,
  1,1,0,
  2,0,1,
  2,0,0,
  7,1,1,
  3,1,0,
  7,1,1,
  2,1,0,
  1,1,0,
  6,1,1
) |>
  mutate(subject = row_number()) 

treatment_df <- simple_df |>
  filter(treatment == 1) |>
  arrange(time)

treatment <- tibble()
time <- treatment_df |> distinct(time) |> pull()
at_risk <- nrow(treatment_df)
S_t <- 1

for (i in time) {
  df_i <- treatment_df |>
    filter(time == i)
  status <- df_i |> pull(status) |> sum()
  n <- df_i |> pull(status) |> length()
  S_t <- S_t * (1 - status/at_risk)
  treatment <- treatment |>
    bind_rows(tibble(time=i,at_risk=at_risk,S_t=S_t,treatment=1))
  at_risk <- at_risk - n
}

(treatment)

## # A tibble: 4 × 4
##    time at_risk   S_t treatment
##   <dbl>   <int> <dbl>     <dbl>
## 1     2       6   1           1
## 2     5       5   0.8         1
## 3     6       4   0.6         1
## 4     7       2   0           1

no_treatment_df <- simple_df |>
  filter(treatment == 0) |>
  arrange(time)

no_treatment <- tibble()
time <- no_treatment_df |> distinct(time) |> pull()
at_risk <- nrow(no_treatment_df)
S_t <- 1

for (i in time) {
  df_i <- no_treatment_df |>
    filter(time == i)
  status <- df_i |> pull(status) |> sum()
  n <- df_i |> pull(status) |> length()
  S_t <- S_t * (1 - status/at_risk)
  no_treatment <- no_treatment |>
    bind_rows(tibble(time=i,at_risk=at_risk,S_t=S_t,treatment=0))
  at_risk <- at_risk - n
}

(no_treatment)

## # A tibble: 3 × 4
##    time at_risk   S_t treatment
##   <dbl>   <int> <dbl>     <dbl>
## 1     1       6 0.667         0
## 2     2       4 0.333         0
## 3     3       1 0             0

Visualize

rbind(treatment,no_treatment) |>
  bind_rows(tibble(
    time=c(0,0), status=c(0,0), treatment=c(1,0), subject=c(0,0), at_risk=c(6,6), S_t=c(1,1)  #add initial phase 
  )) |>
  ggplot(aes(x=time,y=S_t,color=as.factor(treatment))) +
  geom_step() +
  theme_bw()

Wow, since we created this simpel dataset, knowing treatment extended the time event, whereas no treatment didn’t, we can nicely see that when creating 2 KM plots stratified by the treatment group and plot it all they look very different. Now let’s quickly look at log rank test with survival package and then calculate by hand and see if we can reproduce the same p value.

## log-rank test
(log_rank_test <- survival::survdiff(Surv(time, status) ~ treatment, data = simple_df))

## Call:
## survival::survdiff(formula = Surv(time, status) ~ treatment, 
##     data = simple_df)
## 
##             N Observed Expected (O-E)^2/E (O-E)^2/V
## treatment=0 6        5     1.97      4.68      9.02
## treatment=1 6        4     7.03      1.31      9.02
## 
##  Chisq= 9  on 1 degrees of freedom, p= 0.003

Alright! It looks like a chi-square test and has a p val of 0.003. Let’s see if we can reproduce that. And it also looked like they use (O-E)^2/V as opposed to sum of (O-E)^2/E like the usual chi-square test to get the chisq statistic. Interesting.

V (variance) = n_0 * n_1 * d * (n - d) / (n^2 * (n - 1))

## log-rank test by hand
n0 <- 6
n1 <- 6

### time 1
simple_df |>
  filter(time == 1) |>
  group_by(treatment) |>
  summarize(n = n(),
            d = sum(status))

## # A tibble: 1 × 3
##   treatment     n     d
##       <dbl> <int> <dbl>
## 1         0     2     2

(t_1 <- tibble(n0=n0, n1=n1, d0=2, d1=0) |>
  mutate(n = n0 + n1) |>
  mutate(d = d0 + d1) |>
  mutate(V = n0 * n1 * d * (n - d) / (n^2 * (n - 1))) |>
  mutate(E0 = (n0 /n) * d) |>
  mutate(chi_i = (d0-E0)^2/V))

## # A tibble: 1 × 9
##      n0    n1    d0    d1     n     d     V    E0 chi_i
##   <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1     6     6     2     0    12     2 0.455     1   2.2

n0 <- n0 - 2
n1 <- n1 - 0

### time 2
simple_df |>
  filter(time == 2) |>
  group_by(treatment) |>
  summarize(n = n(),
            d = sum(status))

## # A tibble: 2 × 3
##   treatment     n     d
##       <dbl> <int> <dbl>
## 1         0     3     2
## 2         1     1     0

(t_2 <- tibble(n0=n0, n1=n1, d0=2, d1=0) |>
  mutate(n = n0 + n1) |>
  mutate(d = d0 + d1) |>
  mutate(V = n0 * n1 * d * (n - d) / (n^2 * (n - 1))) |>
  mutate(E0 = (n0 /n) * d) |>
  mutate(chi_i = (d0-E0)^2/V))

## # A tibble: 1 × 9
##      n0    n1    d0    d1     n     d     V    E0 chi_i
##   <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1     4     6     2     0    10     2 0.427   0.8  3.37

n0 <- n0 - 3
n1 <- n1 - 1

### time 3
simple_df |>
  filter(time == 3) |>
  group_by(treatment) |>
  summarize(n = n(),
            d = sum(status))

## # A tibble: 1 × 3
##   treatment     n     d
##       <dbl> <int> <dbl>
## 1         0     1     1

(t_3 <- tibble(n0=n0, n1=n1, d0=1, d1=0) |>
  mutate(n = n0 + n1) |>
  mutate(d = d0 + d1) |>
  mutate(V = n0 * n1 * d * (n - d) / (n^2 * (n - 1))) |>
  mutate(E0 = (n0 /n) * d) |>
  mutate(chi_i = (d0-E0)^2/V))

## # A tibble: 1 × 9
##      n0    n1    d0    d1     n     d     V    E0 chi_i
##   <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1     1     5     1     0     6     1 0.139 0.167     5

n0 <- n0 - 1
n1 <- n1 - 0

### time 5
simple_df |>
  filter(time == 5) |>
  group_by(treatment) |>
  summarize(n = n(),
            d = sum(status))

## # A tibble: 1 × 3
##   treatment     n     d
##       <dbl> <int> <dbl>
## 1         1     1     1

(t_5 <- tibble(n0=n0, n1=n1, d0=0, d1=1) |>
  mutate(n = n0 + n1) |>
  mutate(d = d0 + d1) |>
  mutate(V = n0 * n1 * d * (n - d) / (n^2 * (n - 1))) |>
  mutate(E0 = (n0 /n) * d) |>
  mutate(chi_i = (d0-E0)^2/V))

## # A tibble: 1 × 9
##      n0    n1    d0    d1     n     d     V    E0 chi_i
##   <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1     0     5     0     1     5     1     0     0   NaN

n1 <- n1 - 1

### time 6
simple_df |>
  filter(time == 6) |>
  group_by(treatment) |>
  summarize(n = n(),
            d = sum(status))

## # A tibble: 1 × 3
##   treatment     n     d
##       <dbl> <int> <dbl>
## 1         1     2     1

(t_6 <- tibble(n0=n0, n1=n1, d0=0, d1=1) |>
  mutate(n = n0 + n1) |>
  mutate(d = d0 + d1) |>
  mutate(V = n0 * n1 * d * (n - d) / (n^2 * (n - 1))) |>
  mutate(E0 = (n0 /n) * d) |>
  mutate(chi_i = (d0-E0)^2/V))

## # A tibble: 1 × 9
##      n0    n1    d0    d1     n     d     V    E0 chi_i
##   <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1     0     4     0     1     4     1     0     0   NaN

n1 <- n1 - 2

### time 7
simple_df |>
  filter(time == 7) |>
  group_by(treatment) |>
  summarize(n = n(),
            d = sum(status))

## # A tibble: 1 × 3
##   treatment     n     d
##       <dbl> <int> <dbl>
## 1         1     2     2

(t_7 <- tibble(n0=n0, n1=n1, d0=0, d1=2) |>
  mutate(n = n0 + n1) |>
  mutate(d = d0 + d1) |>
  mutate(V = n0 * n1 * d * (n - d) / (n^2 * (n - 1))) |>
  mutate(E0 = (n0 /n) * d) |>
  mutate(chi_i = (d0-E0)^2/V))

## # A tibble: 1 × 9
##      n0    n1    d0    d1     n     d     V    E0 chi_i
##   <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1     0     2     0     2     2     2     0     0   NaN

n1 <- n1 - 2

bind_rows(t_1, t_2, t_3, t_5, t_6, t_7) |>
  mutate(V = replace_na(V, 0),
         E0 = replace_na(E0, 0)) |>
  summarise(
    O0 = sum(d0),
    E0 = sum(E0),
    V  = sum(V)
  ) |>
  mutate(chi_sq = (O0 - E0)^2 / V)

## # A tibble: 1 × 4
##      O0    E0     V chi_sq
##   <dbl> <dbl> <dbl>  <dbl>
## 1     5  1.97  1.02   9.02

pchisq(q = 9.02, df = 1, lower.tail = F)

## [1] 0.002670414

We got it! If we round it up, it’s exactly 0.003 just like from survival.

Notice that we used E0, but you can use E1 and it would should return the same chi square statistic. Click below for details.

bind_rows(t_1, t_2, t_3, t_5, t_6, t_7) |>
  mutate(E1 = (n1 / n) * d) |>
  summarise(
    O1 = sum(d1),
    E1 = sum(E1),
    V = sum(V)
  ) |>
  mutate(chi_sq = (O1-E1)^2/V) 

## # A tibble: 1 × 4
##      O1    E1     V chi_sq
##   <dbl> <dbl> <dbl>  <dbl>
## 1     4  7.03  1.02   9.02

Take note that KM estimator can only estimate survival function and you can only compare the survival curves with log-rank test, but you can’t add more variables to adjust for confounding. Which means, we assume that there isn’t any confounding factors between treatment groups.

If any adjustment that’s needed, that’s where Cox proportional hazard model comes in. Now if we were to add age and run cox model, we would get a different hazard ratio and p value, but the log-rank test would still be the same because log-rank test is only comparing the survival curves without adjusting for any covariates. Let’s see that in action. Click below to expand, you’re going to see an interesting warning, complete separation.

simple_df <- tribble(
  ~time, ~status, ~treatment, ~age,
  5,1,1,30,
  2,1,0,80,
  6,0,1,35,
  1,1,0,85,
  2,0,1,32,
  2,0,0,30,
  7,1,1,25,
  3,1,0,90,
  7,1,1,98,
  2,1,0,98,
  1,1,0,89,
  6,1,1,20
) |>
  mutate(subject = row_number()) 

survival::coxph(Surv(time,status) ~ treatment+age, data = simple_df)

## Warning in coxph.fit(X, Y, istrat, offset, init, control, weights = weights, :
## Loglik converged before variable 1 ; coefficient may be infinite.

## Call:
## survival::coxph(formula = Surv(time, status) ~ treatment + age, 
##     data = simple_df)
## 
##                 coef  exp(coef)   se(coef)      z     p
## treatment -2.202e+01  2.729e-10  1.937e+04 -0.001 0.999
## age        2.275e-03  1.002e+00  1.192e-02  0.191 0.849
## 
## Likelihood ratio test=10.61  on 2 df, p=0.004959
## n= 12, number of events= 9

Let’s simulate the data so that we can estimate a more accurate hazard ratio with cox model and see how it compares to the true hazard ratio that we set in the simulation.

Simulation

library(survival)
library(survminer)

# simulate data of HR 0.55 (95%CI 0.442-0.674)
set.seed(1)
n <- 350
base_event <- 25
base_rate <- 1/base_event
treatment_event <- base_event + 20
treatment_rate <- 1/treatment_event
hr <- treatment_rate/base_rate
coef <- log(hr)
confounder <- rbinom(n,1,0.5)
treatment <- rbinom(n, 1, plogis(0.5*confounder))
true_time <- rexp(n, rate = base_rate*exp(coef*treatment+0.5*confounder))
cens_time <- runif(n, min = 0, max = treatment_event)         

df <- tibble(true_time, cens_time) |>
  mutate(time = pmin(true_time, cens_time),
         status = case_when(
           true_time <= cens_time ~ 1,
           TRUE ~ 0
         )) |>
  mutate(confounder = confounder) |>
  mutate(treatment = treatment |> as.factor())

head(df)

## # A tibble: 6 × 6
##   true_time cens_time  time status confounder treatment
##       <dbl>     <dbl> <dbl>  <dbl>      <int> <fct>    
## 1     13.0      42.3  13.0       1          0 0        
## 2      8.83      2.60  2.60      0          0 0        
## 3      4.49     13.7   4.49      1          1 0        
## 4     61.4      10.8  10.8       0          1 1        
## 5     47.7      17.3  17.3       0          0 1        
## 6     11.6      34.3  11.6       1          1 1

The above simulation, we set the true hazard ratio to be 0.55, which means that the treatment group has a 45% reduction in the hazard of the event compared to the control group. We then simulate the true time to event using an exponential distribution, and we also simulate a censoring time using a uniform distribution. The observed time is the minimum of the true time and the censoring time, and the status variable indicates whether the event was observed (1) or censored (0).

Simulating the above is helpful because we then know the true rate was derived from exponential function based on base rate multiplied by the hazard ratio, so we can then compare the estimated hazard ratio from the Cox model to the true hazard ratio we set in the simulation. The part that connects intuitively is how the exp(coef\*treatment+coef2\*confounder) is similar to the linear regression. If you noticed that we use base_rate*exp(...), it’s essentialy the same as exp(log(base_rate)+coef\*treatment+coef2\*confounder) which is the same as exp(intercept + coef\*treatment + coef2\*confounder), where the intercept is log(base_rate). So, in a way, the Cox model is modeling the log of the hazard function as a linear combination of the covariates, which is similar to how linear regression models the mean of the outcome as a linear combination of the covariates. The difference is that in Cox model, we are modeling the hazard function, which is the instantaneous rate of event occurrence at time t, whereas in linear regression, we are modeling the mean of the outcome variable.

Kaplan-Meier Estimator

(survdiff(Surv(time,status) ~ treatment, data = df))

## Call:
## survdiff(formula = Surv(time, status) ~ treatment, data = df)
## 
##               N Observed Expected (O-E)^2/E (O-E)^2/V
## treatment=0 149       74     62.3      2.20      3.61
## treatment=1 201       87     98.7      1.39      3.61
## 
##  Chisq= 3.6  on 1 degrees of freedom, p= 0.06

km_fit <- survfit(Surv(time, status) ~ treatment, data = df)

the log-rank test did show a significant difference between the two groups, which is expected since we set a true hazard ratio of 0.55 in the simulation. The Kaplan-Meier estimator will give us the estimated survival curves for each group, and we can visually compare them to see the difference in survival between the treatment and control groups.

Cox Proportional Hazard Model

cox_fit <- coxph(Surv(time, status) ~ treatment + confounder, data = df, x = T)
summary(cox_fit)

## Call:
## coxph(formula = Surv(time, status) ~ treatment + confounder, 
##     data = df, x = T)
## 
##   n= 350, number of events= 161 
## 
##               coef exp(coef) se(coef)      z Pr(>|z|)   
## treatment1 -0.3730    0.6887   0.1603 -2.327  0.01999 * 
## confounder  0.5024    1.6526   0.1605  3.130  0.00175 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##            exp(coef) exp(-coef) lower .95 upper .95
## treatment1    0.6887     1.4521     0.503    0.9429
## confounder    1.6526     0.6051     1.207    2.2635
## 
## Concordance= 0.582  (se = 0.024 )
## Likelihood ratio test= 13.34  on 2 df,   p=0.001
## Wald test            = 13.34  on 2 df,   p=0.001
## Score (logrank) test = 13.5  on 2 df,   p=0.001

# plot
ggsurvplot(
  fit = km_fit,
  data = df,
  # pval = TRUE,
  conf.int = TRUE,
  risk.table = TRUE,
  legend.labs = c("Control", "Treatment"),
  title    = "Kaplan-Meier Survival Curves"
)

Notice how the estimated hazard ratio from the Cox model is close to the true hazard ratio of 0.55 that we set in the simulation, our HR is 0.595 (95% CI 0.44-0.81). The Kaplan-Meier plot showed us the survival curves for each group, and we can visually see the difference in survival between the treatment and control groups.

Note: survdiff calculates log-rank test, survfit estimates the survival function, and coxph estimates the hazard ratio adjusting for covariates

There is an interesting article by Hernán that cautioned the use of unadjusted HR and the use unadjusted survival curve (which we did, because it’s based off KM estimator). He also mentioned that first, a single average HR across the entire follow-up can be misleading because the true effect may vary over time. Let’s see if we can apply that to our current plot. Let’s use adjustedCurves and see if it looks different.

library(adjustedCurves)

adjust_curve <- adjustedsurv(
  data = df, 
  ev_time = "time", 
  event = "status", 
  variable = "treatment", 
  method = "direct", 
  outcome_model = cox_fit, 
  conf_int = T)

plot(adjust_curve, conf_int = T)

Interesting. They do look different! Though the paper didn’t directly use this package. He also proposed a solution to more accurately estimate time-varying HR using pooled logistic regression and spline on time as feature. Let’s try this next time! So much to learn! On Emily Zabor’s blog, she did mention there is survival::cox.zph() function which allows us to check the assumption of proportional hazards.

Acknowledgement

Thanks to Emily Zabor’s tutorial and also personal advice on practical usage of survival analysis! Her blog contains so much more advanced topics and some functions and packages I’m planning to use in the future. It’s truly one of the more comprehensive and yet easy to understand tutorials I’ve seen on survival analysis. Thanks Emily!

Opportunities For Improvement

• learn competing risk analysis with fine-gray model
• learn to customize ggsurvplot
• use gtsummary::tbl_regression(exp = TRUE) to further beautify aHR
• test out Hernán’s proposed solution to calculate HR
• let’s test out other dataset such as BMT from SemiCompRisks, Melanoma from MASS,

Lessons learnt

• calculating by hand is helpful because I just realized we can’t just calculate at_risk with one row at a time because of the time occurs at the same time, it would be calculated at the same time, that’s why we used for loop for clarity
• survdiff calculates log-rank test, survfit estimates the survival function, and coxph estimates the hazard ratio adjusting for covariates
• censor is usually a good thing, but it could also mean lost to follow up.

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

Introduction

Random Vector Functional Link (RVFL) networks offer a simple yet powerful alternative to traditional neural networks for tabular data. Instead of learning hidden layers through backpropagation, RVFL generates them randomly (or not, if using a deterministic sequence of quasi-random numbers) and focuses all learning effort on a final, regularized linear model.

Formally, let

be the input data. RVFL networks (the ones described in this blog post) construct a set of nonlinear features by projecting (X) onto a random matrix

and applying an activation function (\(g(\cdot)\)):

These random nonlinear features are then concatenated with the original inputs to form an augmented design matrix:

The model prediction is obtained by fitting a linear model on this expanded space (hence, a nonlinear GLM):

Because (Z) can be high-dimensional and highly redundant, RVFL networks (the ones described in this blog post) rely on Elastic Net regularization (glmnet) to estimate the coefficients:

In this framework, randomness creates a rich pool of nonlinear transformations, while regularization selects and stabilizes the most useful ones. The result is a nonlinear model that combines the flexibility of neural networks with the efficiency and robustness of linear methods.

Of course, this blog post is not a proof of the title. It’s about R package rvflnet. But you can appreciate the high performance of RVFLs on regression, classification and survival analysis, an notably on the controversial Boston dataset (performs on par with Random Forest or Gradient Boosting).

0 – Install package

install.packages("survival", repos = "https://cran.r-project.org") # survival analysis

install.packages("remotes", repos = "https://cran.r-project.org")

devtools::install_github('thierrymoudiki/rvflnet') # Nonlinear glm (RVFL networks)

1 – Regression

set.seed(123)

library(glmnet)
data(Boston, package = "MASS")

# -------------------------
# Data
# -------------------------
X <- as.matrix(Boston[, -14])
y <- Boston$medv

n <- nrow(X)
idx <- sample(1:n, size = round(0.8 * n))

X_train <- X[idx, ]
y_train <- y[idx]

X_test <- X[-idx, ]
y_test <- y[-idx]

# -------------------------
# Grid
# -------------------------
grid <- expand.grid(
  n_hidden = c(175, 200, 225, 250),
  alpha = seq(0.1, 0.5, by=0.2),
  include_original = c(TRUE, FALSE),
  seed = 1,
  stringsAsFactors = FALSE
)

results <- vector("list", nrow(grid))

# -------------------------
# Loop
# -------------------------
for (i in seq_len(nrow(grid))) {

  params <- grid[i, ]

  #cat("\n========================================\n")
  #cat(sprintf("Run %d / %d\n", i, nrow(grid)))
  #print(params)

  # -------------------------
  # Fit model
  # -------------------------
  fit <- rvflnet::rvflnet(
    X_train, y_train,
    n_hidden = params$n_hidden,
    activation = "sigmoid",
    W_type = "gaussian",
    seed = params$seed,
    include_original = params$include_original, # direct link, skip connection or not
    alpha = params$alpha
  )

  # -------------------------
  # Evaluate full lambda path
  # -------------------------
  lambdas <- fit$fit$lambda

  preds <- predict(fit, newx = X_test, s = lambdas)

  rmse_path <- sqrt(colMeans((preds - y_test)^2))

  best_idx <- which.min(rmse_path)

  best_rmse <- rmse_path[best_idx]
  best_lambda <- lambdas[best_idx]

  # -------------------------
  # Sparsity
  # -------------------------
  coef_mat <- coef(fit, s = best_lambda)
  nonzero <- sum(coef_mat[-1, 1] != 0)

  # -------------------------
  # Verbose output
  # -------------------------
  #cat(sprintf("Best RMSE: %.4f\n", best_rmse))
  #cat(sprintf("Best lambda: %.6f\n", best_lambda))
  #cat(sprintf("Non-zero coeffs: %d\n", nonzero))

  # -------------------------
  # Store
  # -------------------------
  results[[i]] <- data.frame(
    n_hidden = params$n_hidden,
    alpha = params$alpha,
    include_original = params$include_original,
    seed = params$seed,
    rmse = best_rmse,
    lambda = best_lambda,
    nonzero = nonzero
  )
}

# -------------------------
# Aggregate
# -------------------------
results_df <- do.call(rbind, results)
results_df <- results_df[order(results_df$rmse), ]
print(head(results_df))

Loading required package: Matrix

Loaded glmnet 4.1-10



               n_hidden alpha include_original seed     rmse     lambda nonzero
s= 0.027561759      200   0.1             TRUE    1 2.881935 0.02756176     190
s= 0.017620327      200   0.3             TRUE    1 2.884739 0.01762033     167
s= 0.012734248      200   0.5             TRUE    1 2.889339 0.01273425     158
s= 0.036435024      175   0.1             TRUE    1 2.920012 0.03643502     165
s= 0.016833926      175   0.5             TRUE    1 2.938472 0.01683393     136
s= 0.023293035      175   0.3             TRUE    1 2.941267 0.02329304     144

An RMSE of 2.88 is on par with Random Forest or Gradient Boosting, with a significantly faster computation time.

2 - Classification

2 - 1 Binary Classification

set.seed(123)

data(iris)

# Binary classification: setosa vs others
y <- ifelse(iris$Species == "setosa", 1, 0)
X <- as.matrix(iris[, 1:4])

# Train/test split
n <- nrow(X)
idx <- sample(1:n, size = round(0.8 * n))

X_train <- X[idx, ]
y_train <- y[idx]

X_test <- X[-idx, ]
y_test <- y[-idx]

# -------------------------
# Fit model
# -------------------------
cv_model <- rvflnet::cv.rvflnet(
  X_train, y_train,
  n_hidden = 50,
  activation = "relu",
  W_type = "gaussian",
  family = "binomial",
  nfolds = 5
)

# -------------------------
# Predictions (probabilities)
# -------------------------
(probs <- predict(cv_model, X_test, type = "response"))

# Convert to class
y_pred <- ifelse(probs > 0.5, 1, 0)

all.equal(as.numeric(y_pred), as.numeric(predict(cv_model, X_test, type="class")))

# -------------------------
# Diagnostics
# -------------------------

# Accuracy
acc <- mean(drop(y_pred) == y_test)
cat("Accuracy:", acc, "\n")

# Confusion matrix
table(Predicted = y_pred, Actual = y_test)

TRUE

Accuracy: 1 



         Actual
Predicted  0  1
        0 20  0
        1  0 10

2 - 2 Multiclass Classification

set.seed(123)

data(iris)

y <- as.numeric(iris$Species)
X <- as.matrix(iris[, 1:4])

# Train/test split
n <- nrow(X)
idx <- sample(1:n, size = round(0.8 * n))

X_train <- X[idx, ]
y_train <- y[idx]

X_test <- X[-idx, ]
y_test <- y[-idx]

# -------------------------
# Fit model
# -------------------------
cv_model <- rvflnet::rvflnet(
  X_train, y_train,
  n_hidden = 50,
  activation = "relu",
  W_type = "gaussian",
  family = "multinomial",
  nlambda = 25,
  nfolds = 5
)

# -------------------------
# Diagnostics
# -------------------------

# Accuracy
acc <- colMeans(predict(cv_model, X_test, type="class") == y_test)
cat("Accuracies:", acc, "\n") # consider other metrics

Accuracies: 0.1666667 0.7666667 0.9333333 0.9666667 0.9666667 0.9666667 0.9666667 0.9666667 0.9666667 0.9666667 0.9666667 0.9666667 0.9666667 0.9666667 0.9666667 0.9666667 0.9666667 0.9666667 0.9666667 0.9666667 0.9666667 0.9666667 0.9666667 

3 - Nonlinear Cox survival analysis

3 - 1 Example 1

library(survival)
library(rvflnet)

data(ovarian)

X <- as.matrix(ovarian[, c("age", "resid.ds", "rx", "ecog.ps")])
y <- Surv(ovarian$futime, ovarian$fustat)

set.seed(123)
n <- nrow(X)
train_idx <- sample(1:n, size = round(0.8 * n))

X_train <- X[train_idx, ]
X_test  <- X[-train_idx, ]
y_train <- y[train_idx]
y_test  <- y[-train_idx]

# -------------------------
# Fit model
# -------------------------
cv_fit <- rvflnet::cv.rvflnet(
  X_train, y_train,
  family = "cox",
  nfolds = 5,
  type.measure = "C"
)

plot(cv_fit)

# Out-of-sample C-index
print(glmnet::Cindex(pred = predict(cv_fit, X_test), y = y_test))


Warning message in data(ovarian):
“data set ‘ovarian’ not found”


[1] 0.8571429

3 - 2 Example 2

library(glmnet)
library(survival)

data(pbc)
pbc2       <- pbc[!is.na(pbc$trt), ]
pbc2$event <- as.integer(pbc$status[!is.na(pbc$trt)] == 2)
pbc2$sex_n <- as.integer(pbc2$sex == "f")

feat_cols <- c("trt","age","sex_n","ascites","hepato","spiders","edema",
               "bili","chol","albumin","copper","alk.phos","ast",
               "trig","platelet","protime","stage")

df <- pbc2[, c("time", "event", feat_cols)]
for (col in feat_cols)
  if (any(is.na(df[[col]])))
    df[[col]][is.na(df[[col]])] <- median(df[[col]], na.rm = TRUE)

set.seed(42)
idx_train <- sample(nrow(df), floor(0.75 * nrow(df)))
train <- df[idx_train, ]; test <- df[-idx_train, ]
X_tr  <- as.matrix(train[, feat_cols])
X_te  <- as.matrix(test[,  feat_cols])
y_tr   <- Surv(train$time, train$event)

fit <- rvflnet::rvflnet(
  X_tr, y_tr,
  family = "cox",
  alpha=0.1, lambda=0.1 # not recommended
)

y_te   <- Surv(test$time, test$event)
ci <- glmnet::Cindex(predict(fit, X_te), y_te)

cat("\n=== Test-set C-index ===\n")
print(ci)


=== Test-set C-index ===
[1] 0.8218117

fit <- rvflnet::rvflnet(
  X_tr, y_tr,
  family = "cox",
  alpha=0.1, nlambda=50
)

y_te   <- Surv(test$time, test$event)

(cis <- apply(predict(fit, X_te), 2, function(x) glmnet::Cindex(x, y_te)))

#cat("\n=== Test-set C-index ===\n")
plot(log(fit$fit$lambda), cis, type = 'l')
abline(h=0.8, lty=2, col="red")

s0

0.5

s1

0.762812872467223

s2

0.802145411203814

s3

0.811084624553039

s4

0.811680572109654

s5

0.814064362336114

s6

0.815852205005959

s7

0.817640047675805

s8

0.820023837902265

s9

0.81942789034565

s10

0.817640047675805

s11

0.81823599523242

s12

0.81823599523242

s13

0.815852205005959

s14

0.814660309892729

s15

0.813468414779499

s16

0.813468414779499

s17

0.815852205005959

s18

0.814660309892729

s19

0.82061978545888

s20

0.81942789034565

s21

0.82181168057211

s22

0.82061978545888

s23

0.817044100119189

s24

0.817640047675805

s25

0.81823599523242

s26

0.814660309892729

s27

0.810488676996424

s28

0.803933253873659

s29

0.802145411203814

s30

0.799761620977354

s31

0.793206197854589

s32

0.789034564958284

s33

0.777711561382598

s34

0.771156138259833

s35

0.766984505363528

s36

0.756853396901073

s37

0.748510131108462

s38

0.743146603098927

s39

0.735399284862932

s40

0.728843861740167

s41

0.721692491060787

s42

0.718116805721096

s43

0.717520858164482

s44

0.716924910607867

s45

0.716924910607867

s46

0.715733015494636

s47

0.716328963051251

s48

0.715137067938021

s49

0.713945172824791

Frank Harrell’s Regression Modeling Strategies online seminar will take place May 14, 15, 18, and 19.

This workshop covers principled strategies for building, validating, and interpreting multivariable regression models for a wide range of outcomes, with emphasis on predictive accuracy, avoiding overfitting, and interpreting estimated effects. It explores spline methods, data reduction, benefits of Bayesian modeling, robust semiparametric ordinal, longitudinal, and survival models, and rigorous resampling-based validation, illustrated with applied case studies and R examples. More details here.

Along with the 1-day Introduction to R, Regression, and the rms Package, these virtual seminars are offered through Instats, in association with the A.S.A.

Join our workshop on Reactive Shiny Apps and Deployment with Google Cloud Run: Intermediate R Shiny Workshop,  which is a part of our workshops for Ukraine series! 

Here’s some more info: 

Title: Reactive Shiny Apps and Deployment with Google Cloud Run: Intermediate R Shiny Workshop 

Date: Thursday, May 21st, 18:00 – 20:00 CET (Rome, Berlin, Paris timezone) 

Speaker:Alfredo Hernández Sánchez is a Marie Skłodowska Curie ERA Postdoctoral Fellow at Vilnius University, where he leads the FIRSA project on financial regulation and innovation in Europe. His work combines applied research, data analysis, and reproducible computational methods, with a strong interest in turning research outputs into accessible digital tools such as dashboards and interactive web applications. He works extensively with R, Quarto, and Shiny in academic and policy-oriented settings.

www.alfredohs.com 

Description:  This workshop is designed for people who already know the basics of Shiny and want to build apps that are more robust, more reactive, and easier to maintain. We will look at practical reactive patterns, app structure, and some common choices that make Shiny dashboards easier to develop as they become more complex. In the second part of the workshop, I will show how a Shiny app can move from local development to a public deployment on Google Cloud Run, using a real dashboard project as an example. The session will give participants a practical introduction to a cloud-based workflow for publishing and maintaining Shiny applications in highly customizable environment. Basic familiarity with Shiny is assumed, and some previous experience building simple apps will help participants get the most out of the session.

Minimal registration fee: 20 euro (or 20 USD or 800 UAH)

Please note that the registration confirmation is sent 1 day before the workshop to all registered participants rather than immediately after registration

How can I register?

• Go to https://bit.ly/3wvwMA6 or https://bit.ly/4aD5LMC  or  https://bit.ly/3PFxtNA and donate at least 20 euro. Feel free to donate more if you can, all proceeds go directly to support Ukraine.

• Save your donation receipt (after the donation is processed, there is an option to enter your email address on the website to which the donation receipt is sent)

• Fill in the registration form, attaching a screenshot of a donation receipt (please attach the screenshot of the donation receipt that was emailed to you rather than the page you see after donation).

If you are not personally interested in attending, you can also contribute by sponsoring a participation of a student, who will then be able to participate for free. If you choose to sponsor a student, all proceeds will also go directly to organisations working in Ukraine. You can either sponsor a particular student or you can leave it up to us so that we can allocate the sponsored place to students who have signed up for the waiting list.

How can I sponsor a student?

• Go to https://bit.ly/3wvwMA6 or https://bit.ly/4aD5LMC  or https://bit.ly/3PFxtNA and donate at least 20 euro (or 17 GBP or 20 USD or 800 UAH). Feel free to donate more if you can, all proceeds go to support Ukraine!

• Save your donation receipt (after the donation is processed, there is an option to enter your email address on the website to which the donation receipt is sent)

• Fill in the sponsorship form, attaching the screenshot of the donation receipt (please attach the screenshot of the donation receipt that was emailed to you rather than the page you see after the donation). You can indicate whether you want to sponsor a particular student or we can allocate this spot ourselves to the students from the waiting list. You can also indicate whether you prefer us to prioritize students from developing countries when assigning place(s) that you sponsored.

If you are a university student and cannot afford the registration fee, you can also sign up for the waiting list here. (Note that you are not guaranteed to participate by signing up for the waiting list).

You can also find more information about this workshop series,  a schedule of our future workshops as well as a list of our past workshops which you can get the recordings & materials here.

Looking forward to seeing you during the workshop!

Dear rOpenSci friends, it’s time for our monthly news roundup! You can read this post on our blog. Now let’s dive into the activity at and around rOpenSci!

Tomáš Kalibera (1978–2026)

The rOpenSci team is deeply saddened at the loss of Tomáš Kalibera, a member of R-Core and contributor to the R community for almost 10 years. Tomáš passed away on April 1.  Our thoughts are with Tomáš’s friends and family at this time. 

Read Jan Vitek’s remembrance of Tomáš.

rOpenSci HQ

New editors Alec Robitaille and Lucy D’Agostino McGowan

We’re excited to welcome Alec Robitaille and Lucy D’Agostino McGowan as new editors. Alec joins our general review team, and Lucy our statistical software review team. Read more in the post introducing them!

Champions Program update

We’re excited to share that we’ve finished selecting the new cohort for the rOpenSci Champions Program! This was not an easy process. The quality of the proposals was exceptionally high, which made the selection both challenging and inspiring. We’re grateful to everyone who applied and shared their ideas with us. Please join us in welcoming our new Champions: Bastián Olea Herrera, Durga Valentina Linares Herrera, José Daniel Conejeros, Denisse Fierro Arcos, Evelia Lorena Coss Navarrete, Gladys Choque Ulloa, Linda Jazmín Cabrera Orellana, Patricia Andrea Loto, Marina Cecilia Cock, María Florencia Tames, and Estefanía Torrejón.

Over the coming months, they will contribute to the R Community through developing new packages, reviewing packages, and submitting packages for peer review. We’re looking forward to working with this amazing group and supporting their projects!

Collaborating between Bioconductor and R-universe on Development of Common Infrastructure

Bioconductor is collaborating with R-universe to gradually modernize parts of its infrastructure, while accommodating the project’s scale, governance, and established processes. In turn, Bioconductor is helping R-universe expand and refine its features as we learn to serve the complex needs of the Bioconductor community. Read more in the blog post.

rOpenSci Staff presentations

Yanina Bellini Saibene at R/Medicine 2026

Yani will deliver her keynote talk “Software Sustainability and Community Management” on Thursday May 7th, 11:15AM–12:15PM ET.

Jeroen Ooms at “Where Do R Packages Live?”

Jeroen will take part in an online discussion panel on Wednesday 20 May at 5:00 PM – 6:00 PM (AEST). The panel is organized by the Statistical Computing and Visualisation section of the Statistical Society of Australia (details).

Updates to the goodpractice R package

We have long recommended the goodpractice package, which identifies issues with R packages, and advises how to fix them. Thanks to a huge amount of work by valued community member Athanasia Mo Mowinckel, goodpractice has been extended and improved to include entirely new suites of checks, and improved ability to control which checks are run. A blog post describing the updates will be published soon, but in the meantime, we encourage you to install the current development version and try it out yourself.

Analyse your targets pipeline

Following our blog post about tree-sitter, Tyler Morgan Wall was inspired to create a static analysis tool for targets pipelines (as a VS Code extension).

Coworking

Read all about coworking!

• Tuesday May 5th 2026, 9:00 Australia Western (01:00 UTC) “Code Review with rOpenSci” with Steffi LaZerte and cohost Liz Hare.
  • Explore resources for Code Review
  • Sign up to volunteer to do software peer-review at rOpenSci
  • Meet cohost, Liz Hare, and discuss resources for Code Review with rOpenSci.
• Tuesday June 2nd 2026, 14:00 Europe Central (12:00 UTC) [theme to be determined], with Steffi LaZerte and cohost to be determined.
  • Explore resources related to the theme
  • Meet the cohost, and other attendees, and discuss the theme or other topics.

And remember, you can always cowork independently on work related to R, work on packages that tend to be neglected, or work on what ever you need to get done!

Editors’ Office Hours

We are exploring hosting a new event, Editors’ Office Hours, where you can drop in to ask questions about rOpenSci Software Peer Review on or near the third Tuesday of each month, alternating among timezones to accommodate different parts of the world.

Upcoming office hours:

• Tuesday May 19, 16:00-17:00 Europe Central (14:00-15:00 UTC) (event)

useR! 2026 Diversity Scholarship Program

useR! 2026 is offering diversity scholarships to support participation from people in underrepresented or historically marginalized groups within the R community. The program includes both registration fee waivers and full needs-based scholarships, which cover conference registration as well as travel and lodging (via reimbursement). Applications are open to eligible participants worldwide until May 10, and will be reviewed based on need, eligibility, and potential impact by a committee from Forwards, RLadies+ Global, and rOpenSci.

Found all the details and important links on the conference website: https://user2026.r-project.org/additional/diversity_scholarship.html

Software

New packages

The following package recently became a part of our software suite:

• reviser, developed by Marc Burri together with Philipp Wegmueller: Analyzes revisions in real-time time series vintages. The package converts between wide revision triangles and tidy long vintages, extracts selected releases, computes revision series, visualizes vintage paths, and summarizes revision properties such as bias, dispersion, autocorrelation, and news-noise diagnostics. It also identifies efficient releases and estimates state-space models for revision nowcasting. Methods are based on Howrey (1978) https://doi.org/10.2307/1924972, Jacobs and Van Norden (2011) https://doi.org/10.1016/j.jeconom.2010.04.010, and Kishor and Koenig (2012) https://doi.org/10.1198/jbes.2010.08169. It has been reviewed.

Discover more packages, read more about Software Peer Review.

New versions

The following ten packages have had an update since the last newsletter: osmextract (v0.6.0), Athlytics (v1.0.5), emodnet.wfs (v2.1.2), fellingdater (v1.2.1), readODS (v2.3.5), git2rdata (v0.5.2), weatherOz (v3.0.0), promoutils (v0.5.0), allcontributors (v0.2.3), and reviser (v0.1.1).

The writexl package has a new maintainer, Bill Denney. NLMR is now maintained by Jakub Nowosad.

Software Peer Review

There are eighteen recently closed and active submissions and 4 submissions on hold. Issues are at different stages:

• One at ‘6/approved’:

  • reviser, Tools for Studying Revision Properties in Real-Time Time Series Vintages. Submitted by Marc Burri. (Stats).

• Two at ‘5/awaiting-reviewer(s)-response’:

  • saperlipopette, Create Example Git Messes. Submitted by Maëlle Salmon.

  • pkgmatch, Find R Packages Matching Either Descriptions or Other R Packages. Submitted by mark padgham.

• Four at ‘4/review(s)-in-awaiting-changes’:

  • pvEBayes, Empirical Bayes Methods for Pharmacovigilance. Submitted by Yihao Tan. (Stats).

  • logolink, An Interface for Running NetLogo Simulations. Submitted by Daniel Vartanian.

  • ActiGlobe, Wearable Recording Processor for Time Shift Adjustment and Data Analysis. Submitted by C. William Yao.

  • galamm, Generalized Additive Latent and Mixed Models. Submitted by Øystein Sørensen. (Stats).

• Five at ‘3/reviewer(s)-assigned’:

  • nycOpenData, Convenient Access to NYC Open Data API Endpoints. Submitted by Christian Martinez.

  • RAMEN, RAMEN: Regional Association of Methylome variability with the Exposome and geNome. Submitted by Erick Navarro-Delgado.

  • ernest, A Toolkit for Nested Sampling. Submitted by Kyle Dewsnap. (Stats).

  • rcrisp, Automate the Delineation of Urban River Spaces. Submitted by Claudiu Forgaci. (Stats).

  • priorsense, Prior Diagnostics and Sensitivity Analysis. Submitted by Noa Kallioinen. (Stats).

• Three at ‘2/seeking-reviewer(s)’:

  • EpiStrainDynamics, Infer temporal trends of multiple pathogens. Submitted by Saras Windecker. (Stats).

  • lakefetch, Calculate Fetch and Wave Exposure for Lake Sampling Points. Submitted by jeremylfarrell.

  • fcmconfr, Fuzzy Cognitive Map Analysis in R. Submitted by benroston. (Stats).

• Three at ‘1/editor-checks’:

  • ciecl, International Classification of Diseases ICD-10/ICD-11 for Chile. Submitted by Rodolfo Tasso.

  • RAQSAPI, A Simple Interface to the US EPA Air Quality System Data Mart API. Submitted by mccroweyclinton-EPA.

  • coevolve, Fit Bayesian Generalized Dynamic Phylogenetic Models using Stan. Submitted by Scott Claessens. (Stats).

Find out more about Software Peer Review and how to get involved.

On the blog

Software Review

• reviser: Analyzing Real-Time Data Revisions in R by Marc Burri. A short introduction to reviser for analyzing real-time data vintages and revisions in R.

• Expanding the Editorial Team: Alec Robitaille and Lucy D’Agostino McGowan Join as Editors by Alec Robitaille, Lucy D’Agostino McGowan, and Yanina Bellini Saibene. Introducing two new editors for rOpenSci statistical software peer review.

• Collaborating between Bioconductor and R-universe on Development of Common Infrastructure by The rOpenSci Team and The Bioconductor Team. R-consortium ISC top-level project R-universe is working with Bioconductor to help gradually modernize parts of its infrastructure, while accommodating the project’s scale, governance, and established processes.

Tech Notes

• A Better R Programming Experience Thanks to Tree-sitter by Maëlle Salmon. Modern tooling for parsing, searching, formatting, editing R code, just like for other programming languages.

Calls for contributions

Calls for maintainers

If you’re interested in maintaining any of the R packages below, you might enjoy reading our blog post What Does It Mean to Maintain a Package?.

• landscapetools, R package for some of the less-glamorous tasks involved in landscape analysis. Issue for volunteering.

• hddtools, Tools to discover hydrological data, accessing catalogues and databases from various data providers. Issue for volunteering.

• qualtRics, download Qualtrics survey data. Issue for volunteering.

Calls for contributions

Refer to our help wanted page – before opening a PR, we recommend asking in the issue whether help is still needed.

Package development corner

Some useful information for R package developers.

Useless code, duplicated function? Refactoring with Jarl

The Jarl CLI by Étienne Bacher received several useful new features for package developers:

• unused_function
• duplicated_function_definition.

They are a nice complement to unreachable_code.

Read more in the release announcement.

Git commands to get to know a project

Ally Piechowski wrote an insightful post entitled “The Git Commands I Run Before Reading Any Code”, which suggests Git commands which are useful for understanding the code base. For instance a command to determine which files recently changed the most!

To complement this post, Garrick Aden-Buie wrote “a little git-recon bash script that runs them in series, complete with some ascii bar plots”.

Another R mailing list archive

Last month we listed several ways to browse the R mailing lists. Newsletter reader Florian Kohrt kindly wrote to us to mention the plain text backup maintained by Michael Chirico.

Will R run out of random seeds? Useful seed explainer

Andrew Heiss published an useful and interesting deep dive into random seeds.

Enforcing the coalesce operator

Are you enjoying the coalesce operator %||% introduced in R 4.4.0?

L %||% R newly in base is an expressive idiom for the phrases if(!is.null(L)) L else R or if(is.null(L)) R else L.

Consider enforcing it via Jarl or lintr.

As a reminder, the operator can be used in older versions of R through the backports R package or by importing it from rlang.

Positive AI reading

If you’re feeling some AI dread, you might enjoy:

• The blog post “Antidote” by Vicky Boykis;
• The talk “Is my degree worthless” by Davis Vaughan.

Last words

Thanks for reading! If you want to get involved with rOpenSci, check out our Contributing Guide that can help direct you to the right place, whether you want to make code contributions, non-code contributions, or contribute in other ways like sharing use cases. You can also support our work through donations.

If you haven’t subscribed to our newsletter yet, you can do so via a form. Until it’s time for our next newsletter, you can keep in touch with us via our website and Mastodon account.

dplyr verbs are descriptive: let’s make them more verbose!

Yet another pipe for R.

Repost for better image handling on r-bloggers.

Motivation

In SAS, every DATA step prints a log:

NOTE: There were 120000 observations read from WORK.SALES.
NOTE: 7153 observations were deleted.
NOTE: The data set WORK.RESULT has 112847 observations and 11 variables.

R’s dplyr pipelines are silent. logrittr fills that gap with %>=%, a drop-in pipe that logs row counts, column counts, added/dropped columns, and timing at every step, with no function masking.

With Fira Code ligatures, %>=% renders as a single wide arrow visually similar to %>% with an underline added, like a subtitle or, say, to read between the lines of a pipeline (what happened).

Multiples contexts

Things happens:

NOTE: There were 120000 observations read from WORK.SALES.
NOTE: 120000 observations were deleted.
NOTE: The data set WORK.RESULT has 0 observations and 11 variables.

“It’s here we’ve lost all rows in script execution”.

Pro

Reading this a long time after execution of a script helps you see:

• what happened at each stage of data processing without having to rerun the code, for example in a production environment where the input data is constantly changing
• monitor key processes
• Make sure you can explain what happened (an audit, for example)

In professional contexts it’s often needed.

Educational

This will also be clearer thanks to a console log for those with little experience with the tidyverse: people who are taking their first steps in programming by following a tutorial or teaching themselves.

Installation

install.packages('logrittr', repos = 'https://guillaumepressiat.r-universe.dev')

# or from github
# devtools::install_github("GuillaumePressiat/logrittr")

See github or r-universe.

Usage

library(logrittr)
library(dplyr)

iris %>=%
  as_tibble() %>=%
  filter(Sepal.Length < 5)  %>=%
  mutate(rn = row_number()) %>=%
  semi_join(
    iris %>% as_tibble() %>=%
      filter(Species == "setosa"),
    by = "Species"
  )  %>=%
  group_by(Species) %>=%
  summarise(n = n_distinct(rn))

── iris  [rows:       150  cols:    5] ─────────────────────────────────────────────────────
ℹ as_tibble()                            rows:       150 +0        cols:    5 +0    [   0.0 ms]
ℹ filter(Sepal.Length < 5)               rows:        22 -128      cols:    5 +0    [   3.0 ms]
ℹ mutate(rn = row_number())              rows:        22 +0        cols:    6 +1    [   1.0 ms]
  added: rn
ℹ > filter(Species == "setosa")          rows:        50 -100      cols:    5 +0    [   1.0 ms]
ℹ semi_join(iris %>% as_tibble() %>=%    rows:        20 -2        cols:    6 +0    [   5.0 ms]
  filter(Species == "setosa"), by =
  "Species")
ℹ group_by(Species)                      rows:        20 +0        cols:    6 +0    [   3.0 ms]
ℹ summarise(n = n_distinct(rn))          rows:         1 -19       cols:    2 -4    [   2.0 ms]
  dropped: Sepal.Length, Sepal.Width, Petal.Length, Petal.Width, rn
  added: n

Screenshot

library(dplyr)
library(logrittr)

logrittr_options(lang = "en", big_mark = ",", wrap_width = NULL, max_cols = 3)

nycflights13::flights %>=% 
  as_tibble() %>=%
  group_by(year, month, day) %>=% 
  count() %>=% 
  tidyr::pivot_wider(values_from = "n", names_from = "day") %>=% 
  glimpse()

Related package: tidylog

tidylog is a really neat package that gives me motivation for this one. tidylog works by masking dplyr functions, which doesn’t seem ideal to me.

Anyway this also was a moment for me to test a new programmer tool that is used a lot for programming at this time.

logrittr uses a custom pipe operator and never touches the dplyr namespace. Its console output is colorful and informative thanks to the cli package.

Working with lumberjack

If you already know the lumberjack package, compatibility is available with logrittr (timings are approximates).

Calling logrittr_logger$new():

library(lumberjack)
library(dplyr)

l <- logrittr_logger$new(verbose = TRUE)
logfile <- tempfile(fileext=".-r.log.csv")

iris %L>%
   start_log(log = l, label = "iris step") %L>%
   as_tibble() %L>%
   filter(Sepal.Length < 5) %L>%
   mutate(rn = row_number()) %L>%
   group_by(Species) %L>%
   summarise(n = n_distinct(rn)) %L>%
   dump_log(file=logfile, stop = FALSE)
   

mtcars %>% 
  start_log(log = l, label = "mtcars step") %L>%
   count() %L>%
   dump_log(file=logfile, stop = TRUE)


logdata <- read.csv(logfile)

Will write logrittr log content of multiple data steps in the same csv file.

Limitations

• Like tidylog, logrittr only works with dplyr pipelines on R data.frames (in memory) and is not able to do so with dbplyr pipelines from databases (remote/lazy table).

• Join cardinalities nicely done in tidylog are difficult to have from the pipe as join is already done, at this time we only show N row and N col evolution (before / after).

• Yes it’s another pipe, not ideal. We can dream of a with_logging(TRUE) context that will activate behaviour of logrittr pipe in |> or in %>%.

Take another pipe for a spin

logrittr prioritizes the user experience with a structured and colorful display in the console.

For now, this package is just a proof of concept that gave me a chance to experiment a bit with the cli package and few other things. But I think there’s a need for that in R, in a specific area where SAS outputs are so informative.

• https://guillaumepressiat.r-universe.dev/logrittr
• https://github.com/guillaumepressiat/logrittr

Introduction

Optimisation Models

grouper provides two distinct integer linear programming optimisation models.

library(grouper)
library(ompr)
library(ompr.roi)
library(ROI.plugin.glpk)

Preference-Based Assignment

The Preference-Based Assignment (PBA) model allows educators to assign student groups to topics to maximise overall student preferences for those topics. The topics can be viewed as project titles. The model allows for repetitions of each project title. This formulation also allows each project team to comprise multiple sub-groups. This is useful in cases where the project requires teams with different functionality to work together, e.g. where one team works on a front-end while the other develops a back-end model.

To execute the optimisation routine, an instructor prepares:

1. 1. 1. A group composition table listing the member students within each self-formed group
      2. A preference matrix containing the preference that each self-formed group has for each topic.
      3. A YAML file defining the remaining parameters of the model.

Examples

pba_gc_ex002
#>   id grouping
#> 1  1        1
#> 2  2        1
#> 3  3        2
#> 4  4        2
#> 5  5        3
#> 6  6        3
#> 7  7        4
#> 8  8        4

Each student is in a self-formed group of size 2, indicated via the grouping column. Suppose that, for this set of students, the instructor wishes to assign students into two topics, with each topic having two sub-groups. This requires the preference matrix to have 4 columns – one for each topic-subgroup combination. Remember that the ordering of topics/subtopics in the preference matrix should be:

Topic1-Subtopic1, Topic2-Subtopic1, Topic1-Subtopic2, Topic2-Subtopic2

Thus there should be 4 rows in the preference matrix – one for each self-formed group.

pba_prefmat_ex002
#>      col1 col2 col3 col4
#> [1,]    4    3    2    1
#> [2,]    3    4    2    1
#> [3,]    1    2    4    3
#> [4,]    1    2    3    4

The YAML file for this model contains the following parameters:

n_topics: 2
B: 2
R: 1
nmin: 2
nmax: 2
rmin: 1
rmax: 1

B corresponds to the number of sub-topics per topic, while r_min and r_max denote the minimum and maximum number of repetitions of each topic. n_min and n_max denote the minimum and maximum number of members in each sub-topic group.

It is possible to assign each self-formed group to its optimal choice of topic-subtopic combination. In our solution, we should see that group 1 is assigned to subtopic 1 of topic 1, group 2 is assigned to sub-topic 1 of topic 2, and so on.

df_ex002_list <- extract_student_info(pba_gc_ex002, "preference", 
                                     self_formed_groups = 2, 
                                     pref_mat = pba_prefmat_ex002)
yaml_ex002_list <- extract_params_yaml(system.file("extdata", 
                                         "pba_params_ex002.yml",  
                                          package = "grouper"),
                                      "preference")
m2 <- prepare_model(df_ex002_list, yaml_ex002_list, "preference")
result2 <- solve_model(m2, with_ROI(solver="glpk"))

assign_groups(result2, assignment = "preference", 
              dframe=pba_gc_ex002, yaml_ex002_list, 
              group_names="grouping")
#>   topic2 subtopic rep group size
#> 1      1        1   1     1    2
#> 2      2        1   1     2    2
#> 3      1        2   1     3    2
#> 4      2        2   1     4    2

Diversity-Based Assignment

The Diversity-Based Assignment (DBA) model enables educators to assign students to groups and topics with the dual, but weighted, aims of maximising diversity (based on student attributes) within groups and balancing specific skill levels across different groups.

To execute the DBA optimisation routine, the instructor prepares:

1. 1. 1. A group composition table containing:
         1. the member students within each self-formed group,
         2. the demographics that will be used to compute pairwise dissimilarity between students, and
         3. a numeric measure of each student’s skill.
      2. A YAML file defining the remaining parameters of the model.

Examples

dba_gc_ex001
#>   id major skill groups
#> 1  1     A     1      1
#> 2  2     A     1      2
#> 3  3     B     3      3
#> 4  4     B     3      4

It is intuitive that an assignment into two groups of size two, based on the diversity of majors alone, should assign students 1 and 2 into the first group and the remaining two students into another group.

The corresponding YAML dba_gc_ex001.yml file for this exercise consists of the following lines:

n_topics:  2
R:  1
nmin: 2
nmax: 2
rmin: 1
rmax: 1

To run the assignment, we can use the following commands. We can use either the gurobi solver, or the glpk solver for this example. Both are equally fast.

# Indicate appropriate columns using integer ids.
df_ex001_list <- extract_student_info(dba_gc_ex001, "diversity",
                                      demographic_cols = 2, 
                                      skills = 3, 
                                      self_formed_groups = 4)

yaml_ex001_list <- extract_params_yaml(system.file("extdata", 
                                         "dba_params_ex001.yml",  
                                         package = "grouper"),
                                       "diversity")
m1 <- prepare_model(df_ex001_list, yaml_ex001_list,
                    assignment="diversity",w1=0.5, w2=0.5)

result3 <- solve_model(m1, with_ROI(solver="glpk"))
assign_groups(result3, assignment = "diversity", 
              dframe=dba_gc_ex001, 
              group_names="groups")
#>   topic rep group id major skill
#> 1     1   1     2  2     A     1
#> 2     1   1     3  3     B     3
#> 3     2   1     1  1     A     1
#> 4     2   1     4  4     B     3

We can see that students 2 and 3 have been assigned to topic 1, repetition 1. Students 1 and 4 have been assigned to topic 2, repetition 1. w₁ and w₂ both have weights 0.5, which means the skills and demographic inputs are given equal weight in the optimisation.

At present, the routines use the daisy function from the cluster package to compute a pairwise dissimilarity matrix between students. However, it is also possible to supply your own custom dissimilarity matrix. Consider the following dataset of 4 students:

dba_gc_ex003
#>   year   major self_groups id
#> 1    1    math           1  1
#> 2    2 history           2  2
#> 3    3    dsds           3  3
#> 4    4    elts           4  4

Now consider a situation where we wish to consider years 1 and 2 different from years 3 and 4, and math and dsds (STEM majors) to be different from elts and history (non-STEM majors). For each difference, we assign a score of 1. This means that students 1 and 2 would have a dissimilarity score of 1 due to their difference in majors. Students 1 and 3 would also have a score of 1, but due to their difference in years. Students 1 and 4 would have score of 2, due to their differences in majors and in years. The overall dissimilarity matrix would be:

d_mat <- matrix(c(0, 1, 1, 2,
                  1, 0, 2, 1,
                  1, 2, 0, 1,
                  2, 1, 1, 0), nrow=4, byrow = TRUE)

To run the optimisation for this model, we can execute the following code:

df_ex003_list <- extract_student_info(dba_gc_ex003, "diversity",
                                       skills = NULL,
                                       self_formed_groups = 3,
                                       d_mat=d_mat)
yaml_ex003_list <- extract_params_yaml(system.file("extdata",   
                                         "dba_params_ex003.yml",
                                         package = "grouper"), 
                                       "diversity")
m3 <- prepare_model(df_ex003_list, yaml_ex003_list, w1=1.0, w2=0.0)
result <- solve_model(m3, with_ROI(solver="glpk")

assign_groups(result, "diversity", dba_gc_ex003,
              group_names="self_groups")
#>   topic rep group year   major id
#> 1     1   1     1    1    math  1
#> 2     1   1     4    4    elts  4
#> 3     2   1     2    2 history  2
#> 4     2   1     3    3    dsds  3

As you can see, the members of the two groups have maximal difference between them – they differ in terms of their year, and in terms of their major. Notice that we specified

skills = NULL

and

w2 = 0.0

This ensures that no skills columns were taken into account in this optimisation.

Gurobi Optimiser

Shiny Applications

library(shiny)
runApp(appDir=system.file("shiny", "dbaWebApp", package="grouper"))

# Analogous code for PBA app:
# runApp(appDir=system.file("shiny", "pbaWebApp", package="grouper"))

Here is a screen shot of the diversity-based shiny application.



The system folders with the shiny apps also contain example csv files for use with the apps.

More Details

• • Define the number of repetitions for each topic.
  • Define the max. and min. number of group members for each topic.

• • CRAN page
  • Github repository

Introduction to describe()

The describe() function from R’s psych package (Revelle, 2023) provides a comprehensive statistical summary of your dataset. Unlike R’s base summary() function, it includes additional metrics that are particularly useful for data exploration and assumption checking.

library(psych)
describe(your_data)

Breaking Down the Output Columns

Here’s what each column in the output represents:

Key Statistics and Their Interpretation

Central Tendency

• Mean vs. Median: Differences indicate skewness
• Trimmed Mean: Removes influence of outliers (default drops top/bottom 10%)

Variability

• SD vs. MAD: Use MAD when outliers are present
• Range: Simple but outlier-sensitive

Distribution Shape

• Skewness:
  • >0: Right-tailed
  • <0: Left-tailed
  • 0: Symmetric
• Kurtosis (Excess):
  • >0: Heavy-tailed (more outliers than normal)
  • <0: Light-tailed

Practical Examples

Example 1: MPG from mtcars

describe(mtcars$mpg)

Output Interpretation:

   vars  n   mean    sd median trimmed   mad min  max range skew kurtosis   se
1     1 32 20.09 6.03   19.2   19.70 5.41 10.4 33.9  23.5 0.61    -0.37 1.07

• Right-skewed (mean > median, positive skew)
• Light-tailed (negative kurtosis)
• SD (6.03) > MAD (5.41): Suggests some outlier influence

When to Use Which Statistic

Extending the Functionality

Adding IQR

The default describe() doesn’t show IQR, but you can add it:

library(dplyr)
describe(mtcars) %>% 
  mutate(IQR = apply(mtcars, 2, IQR, na.rm = TRUE))

Comparing Groups

Use describeBy() for grouped statistics:

describeBy(mtcars$mpg, group = mtcars$cyl)

Conclusion

R’s describe() function provides a powerful starting point for exploratory data analysis. By understanding each statistic it provides, you can:

• Detect data quality issues
• Choose appropriate analysis methods
• Understand your variables’ distributions
• Make informed decisions about data transformations

For formal reporting, consider supplementing these metrics with visualization and statistical tests.

Pro Tip: Always visualize your data alongside these statistics – numbers tell part of the story, but plots reveal the full picture!

Happy coding!

—
Reference:
Revelle, W. (2023). psych: Procedures for Psychological, Psychometric, and Personality Research. Northwestern University.

What Are Coefficients?

In statistics and data analysis, coefficients are numerical measures that quantify relationships between variables or characteristics of data distributions. They serve as fundamental indicators in statistical modeling and data interpretation.

1. Regression Coefficient

Definition

The regression coefficient measures the relationship between an independent variable (X) and a dependent variable (Y).

Formula

For linear model Y = aX + b:

• a: Regression coefficient (change in Y per unit change in X)
• b: Intercept

R Implementation

# Linear regression example
model <- lm(mpg ~ wt, data = mtcars)
summary(model)

# Extract coefficients
coef(model)

Interpretation

A coefficient of -5.34 for vehicle weight (wt) means each additional ton reduces mileage by 5.34 mpg on average.

2. Coefficient of Determination (R²)

Definition

R-squared represents the proportion of variance in the dependent variable explained by the model (0-1 scale).

R Code

# Get R-squared value
summary(model)$r.squared

Guidelines

• R² = 0.75 → Model explains 75% of data variation
• Higher values indicate better model fit

3. Coefficient of Variation (CV)

Definition

CV is a standardized measure of dispersion expressed as percentage of the mean.

Formula

CV% = (Standard Deviation / Mean) × 100%

R Function

# Calculate CV
cv <- function(x) {
  (sd(x, na.rm = TRUE)/mean(x, na.rm = TRUE)) * 100
}

# Example usage
cv(mtcars$mpg)

Interpretation Benchmarks

• CV < 15%: Low variability
• 15-30%: Moderate variability
• >30%: High variability

4. Correlation Coefficient

Definition

Measures the strength and direction of linear relationship between two variables (-1 to 1).

R Implementation

# Calculate correlation
cor(mtcars$mpg, mtcars$wt)

# Correlation matrix
cor(mtcars[, c("mpg", "wt", "hp")])

Interpretation

• 1: Perfect positive correlation
• -1: Perfect negative correlation
• 0: No linear correlation

Other Common Coefficients

Implementation in R

Comprehensive Analysis

library(psych)

# Descriptive statistics (includes multiple coefficients)
describe(mtcars)

# Full regression output
summary(lm(mpg ~ ., data = mtcars))

Custom Coefficient Calculations

# Multi-coefficient function
data_analysis <- function(x) {
  list(
    mean = mean(x),
    sd = sd(x),
    cv = cv(x),
    skewness = moments::skewness(x),
    kurtosis = moments::kurtosis(x)
  )
}

lapply(mtcars[, 1:4], data_analysis)

Visualization

library(ggplot2)
ggplot(mtcars, aes(wt, mpg)) + 
  geom_point() + 
  geom_smooth(method = "lm") +
  labs(title = "MPG vs Weight with Regression Line",
       x = "Weight (tons)",
       y = "Miles per Gallon")

Key Takeaways

1. Select coefficients based on analytical goals:
   • Variable relationships → Regression/Correlation coefficients
   • Model evaluation → R-squared
   • Variability comparison → CV
2. R advantages:
   • Built-in functions for all major coefficients
   • Seamless integration of statistical and visual analysis
3. Best practices:
   • Understand assumptions behind each coefficient
   • Combine statistical results with domain knowledge
   • Clearly distinguish between different coefficients
4. Advanced applications:

   # Robust regression (for outlier-resistant coefficients)
   library(MASS)
   rlm(mpg ~ wt, data = mtcars)
   
   # Standardized coefficients
   library(lm.beta)
   lm.beta(model)

By mastering these statistical coefficients and their R implementations, you’ll be equipped to conduct more rigorous data analysis and communicate results effectively. Remember that coefficients are tools – their proper interpretation always depends on context and research questions.

Happy coding!

Summer School in Science Mapping (SSSM) 2025 – I International Edition

Title: Producing Systematic Literature Reviews with Bibliometrix R and Biblioshiny
Date & Location: June 9-13, 2025 – Naples, ITA

We are pleased to announce the upcoming Summer School in Science Mapping (SSSM) 2025 – I International Edition, an intensive training program focused on conducting Systematic Literature Reviews using the Bibliometrix R package and its shiny-app Biblioshiny.

Organized by the academic spin-off K-Synth in collaboration with the Department of Economics and Statistics at the University of Naples Federico II, the school will be held in Naples, Italy, from June 9 to June 13, 2025.

Aim and Scope

The SSSM 2025 is an intensive training program tailored for early-career researchers and academics seeking to enhance their expertise in bibliometric methods and scientific mapping. By integrating theoretical foundations with practical sessions, the school equips participants with robust skills in citation analysis, co-citation techniques, science mapping, and reproducible workflows for scholarly evaluation. Designed as both a learning and networking opportunity, SSSM 2025 fosters methodological development and international collaboration in a dynamic, research-oriented environment.

The school’s content covers:

– Overview of bibliometric theory and methods
– Query design and data retrieval from major scientific databases
– Descriptive, relational, and structural bibliometric analyses in R
– Practical training in Bibliometrix R package and Biblioshiny app
– Applications to real-world research review cases

Lecturers and Guest Speakers

The school will be led by Professors Massimo Aria and Corrado Cuccurullo, the developers of Bibliometrix and Biblioshiny.

Additionally, the 2025 edition will feature the following keynotes by distinguished international scholars in the field of scientometrics:

– Nicolas Robinson-Garcia (University of Granada), Scientific Director of the Computational Social Sciences and Humanities Unit (U-CHASS)
– Manuel Jesús Cobo Martín (University of Cádiz), lead developer of the SciMAT software
– Nicola De Bellis (University of Modena and Reggio Emilia), Coordinator of the Bibliometric Office and author of influential studies in the evaluation of scientific research

Target Audience and Prerequisites

This Summer School is designed for PhD students, postdoctoral researchers, and academics affiliated with universities or research institutions. Participants are expected to have a basic knowledge of R programming and be familiar with RStudio.

Registration and Fees

Registration is open on the official Bibliometrix website (check the Summer School section):

https://www.bibliometrix.org/sssm/

For any inquiries, feel free to contact the organizing committee at: sssm@bibliometrix.com