Building a sports betting model in R does not begin with machine learning or a complicated statistical formula. It begins with reliable data.

You need historical results, team or player statistics, bookmaker odds and a process for updating everything without manually downloading a new spreadsheet every day. Fortunately, R provides several packages and APIs that make it possible to build a reproducible sports betting data pipeline.

In this guide, you will learn how to obtain sports betting data in R, download current odds, organize historical information and prepare datasets for predictive modeling and backtesting.

What Data Do You Need for a Sports Betting Model?

A useful sports betting dataset normally combines two different types of information:

• Sports performance data: scores, schedules, team statistics, player statistics and play-by-play data.
• Betting market data: moneylines, point spreads, totals, bookmaker prices and historical closing odds.

The exact variables depend on the sport and market you want to predict. For example, an NFL point-spread model may use offensive EPA, defensive EPA, quarterback performance, home advantage, rest days and the bookmaker’s closing spread.

An NBA totals model could use pace, offensive rating, defensive rating, injuries, recent form and the market total.

Useful R Packages for Sports Data

The SportsDataverse ecosystem provides packages for several major sports:

• nflreadr and nflfastR for NFL data.
• hoopR for NBA and NCAA basketball.
• baseballr for MLB, college baseball and Statcast data.
• fastRhockey for NHL and hockey data.
• wehoop for WNBA and women’s college basketball.
• oddsapiR for current and historical sportsbook odds.

Install the core packages with:

install.packages(c(
  "tidyverse",
  "httr2",
  "jsonlite",
  "lubridate",
  "oddsapiR"
))

You do not necessarily need every sport-specific package. Install only the packages required for the leagues you intend to analyze.

Getting a Sports Odds API Key

One of the simplest ways to access bookmaker odds is The Odds API. It covers many sports, leagues, bookmakers and betting markets.

Create an account, obtain your API key and save it in your R environment. Avoid writing a private key directly inside a script that may later be shared online.

install.packages("usethis")
usethis::edit_r_environ()

Add the following line to the .Renviron file:

ODDS_API_KEY=YOUR_PRIVATE_API_KEY

Save the file and restart RStudio. You can then confirm that R can find the key:

Sys.getenv("ODDS_API_KEY")

Do not publish the result of this command or upload your key to GitHub.

Download Current Sports Betting Odds in R

The following example requests current NFL moneyline, spread and total prices from US bookmakers:

library(httr2)
library(jsonlite)
library(dplyr)
library(tidyr)
library(purrr)

api_key <- Sys.getenv("ODDS_API_KEY")

request_url <- paste0(
  "https://api.the-odds-api.com/v4/sports/",
  "americanfootball_nfl/odds"
)

response <- request(request_url) |>
  req_url_query(
    apiKey = api_key,
    regions = "us",
    markets = "h2h,spreads,totals",
    oddsFormat = "decimal",
    dateFormat = "iso"
  ) |>
  req_perform()

odds_raw <- resp_body_json(response, simplifyVector = FALSE)

The API response contains nested JSON because each event can include multiple bookmakers, markets and outcomes. A nested response is useful for storage, but it usually needs to be transformed before modeling.

Convert the API Response into Tidy Data

The following function converts the nested response into one row per event, bookmaker, market and outcome:

tidy_odds <- function(events) {

  map_dfr(events, function(event) {

    map_dfr(event$bookmakers, function(bookmaker) {

      map_dfr(bookmaker$markets, function(market) {

        map_dfr(market$outcomes, function(outcome) {

          tibble(
            event_id = event$id,
            sport = event$sport_title,
            commence_time = event$commence_time,
            home_team = event$home_team,
            away_team = event$away_team,
            bookmaker = bookmaker$title,
            market = market$key,
            outcome = outcome$name,
            odds = outcome$price,
            point = if (is.null(outcome$point)) NA_real_ else outcome$point,
            last_update = bookmaker$last_update
          )
        })
      })
    })
  })
}

odds_df <- tidy_odds(odds_raw)

glimpse(odds_df)

The resulting table can contain columns such as:

• home_team and away_team
• commence_time
• bookmaker
• market
• outcome
• odds
• point

Convert the timestamps into a proper date-time format before analyzing them:

library(lubridate)

odds_df <- odds_df |>
  mutate(
    commence_time = ymd_hms(commence_time),
    last_update = ymd_hms(last_update)
  )

Understanding Moneylines, Spreads and Totals

The API uses different market identifiers:

• h2h: head-to-head or moneyline betting.
• spreads: point-spread or handicap betting.
• totals: over/under markets.

You can filter the dataset to analyze a single market:

spread_odds <- odds_df |>
  filter(market == "spreads")

total_odds <- odds_df |>
  filter(market == "totals")

moneyline_odds <- odds_df |>
  filter(market == "h2h")

Convert Decimal Odds into Implied Probabilities

Decimal odds can be converted into raw implied probability using:

moneyline_odds <- moneyline_odds |>
  mutate(implied_probability = 1 / odds)

For example, decimal odds of 2.00 represent a raw implied probability of 50%. However, bookmaker probabilities normally add up to more than 100% because the prices include a margin, also known as vig or overround.

A simple way to remove this margin is to normalize the probabilities within each event and bookmaker:

fair_moneyline <- moneyline_odds |>
  group_by(event_id, bookmaker) |>
  mutate(
    raw_probability = 1 / odds,
    market_total = sum(raw_probability, na.rm = TRUE),
    fair_probability = raw_probability / market_total
  ) |>
  ungroup()

The resulting fair_probability column provides a basic no-vig market estimate that can be compared with probabilities generated by your model.

How to Collect Historical Betting Odds

A single snapshot is not enough for serious backtesting. You need to store odds repeatedly or use a provider that offers a historical odds endpoint.

Historical data should ideally include:

• The time when the odds were observed.
• The bookmaker.
• The opening price.
• Intermediate market prices.
• The closing price before the game started.
• The final score and betting result.

This distinction matters because a strategy tested against closing odds may produce very different results from one tested against prices available several hours before the game.

When saving a current snapshot, include the collection time:

odds_snapshot <- odds_df |>
  mutate(collected_at = Sys.time())

dir.create("data", showWarnings = FALSE)

file_name <- paste0(
  "data/odds_",
  format(Sys.time(), "%Y%m%d_%H%M%S"),
  ".csv"
)

readr::write_csv(odds_snapshot, file_name)

This creates a new timestamped file every time the script runs. For a larger project, a database such as SQLite or PostgreSQL is more efficient than storing hundreds of CSV files.

Combine Betting Odds with Sports Performance Data

Bookmaker odds become more useful when combined with historical results and predictive features. For NFL analysis, for example, you can use nflreadr to download play-by-play data:

install.packages("nflreadr")

library(nflreadr)
library(dplyr)

pbp <- load_pbp(2025)

team_features <- pbp |>
  filter(!is.na(posteam), !is.na(epa)) |>
  group_by(game_id, posteam) |>
  summarise(
    offensive_epa = mean(epa, na.rm = TRUE),
    success_rate = mean(success == 1, na.rm = TRUE),
    plays = n(),
    .groups = "drop"
  )

You can then aggregate these metrics before each game and join them to the odds table using team names, event dates or a custom event identifier.

For a complete introduction to NFL play-by-play data, EPA and win probability, see Football Analytics with R: NFL Data Science using nflfastR and nflverse.

Sports Data Sources for NFL, NBA, MLB and NHL

NFL Data

The nflreadr and nflfastR ecosystem provides schedules, rosters, player statistics and detailed play-by-play data. It is particularly useful for building features based on EPA, success rate, passing performance and win probability.

NBA Data

The hoopR package can be used to work with NBA and NCAA schedules, box scores and play-by-play information. Potential betting features include pace, offensive efficiency, defensive efficiency, shot profile and recent performance.

MLB Data

The baseballr package provides access to several baseball data sources. Useful variables may include starting pitcher performance, bullpen usage, park factors, batting metrics and Statcast information.

NHL Data

The fastRhockey ecosystem can help analysts access hockey schedules and play-by-play information. Common model features include expected goals, shot quality, goaltender performance, rest and special-teams efficiency.

Build a Simple Probability Model

After cleaning the data and creating features, you can begin with logistic regression. Suppose your dataset contains a binary variable called home_win and several pregame features:

model <- glm(
  home_win ~ home_rating_diff +
    rest_days_diff +
    recent_form_diff +
    market_probability,
  data = training_data,
  family = binomial()
)

test_data <- test_data |>
  mutate(
    predicted_probability = predict(
      model,
      newdata = test_data,
      type = "response"
    )
  )

This is only a baseline. It is usually better to begin with an interpretable model and a clean validation process before trying Random Forest, XGBoost or neural networks.

Identify Potential Value Bets

A potential value bet exists when your estimated probability is higher than the break-even probability implied by the available odds.

betting_candidates <- test_data |>
  mutate(
    break_even_probability = 1 / decimal_odds,
    expected_value = predicted_probability * decimal_odds - 1,
    model_edge = predicted_probability - break_even_probability
  ) |>
  filter(expected_value > 0)

A positive expected value in historical data does not guarantee future profit. Your probabilities must be calibrated, the backtest must avoid data leakage and the strategy must be tested on games that were not used to train the model.

Backtest the Model by Season

Randomly splitting individual games can accidentally allow future information to influence past predictions. A time-based split is generally more realistic.

training_data <- model_data |>
  filter(game_date < as.Date("2025-01-01"))

test_data <- model_data |>
  filter(game_date >= as.Date("2025-01-01"))

A useful backtest should report more than total profit. Consider tracking:

• Number of bets.
• Win rate.
• Return on investment.
• Maximum drawdown.
• Closing line value.
• Brier score.
• Log loss.
• Probability calibration.

If you want to learn how to use Elo ratings, Monte Carlo simulation and forecasting methods, explore Sports Prediction and Simulation with R: Monte Carlo, Elo Ratings, and Forecasting.

Using Bayesian Models for Sports Prediction

Bayesian models are especially useful in sports because team strength changes over time and the amount of available information varies between teams and players.

A Bayesian workflow can:

• Represent uncertainty with probability distributions.
• Update team estimates when new games are played.
• Use partial pooling to stabilize small samples.
• Estimate full predictive distributions instead of single values.
• Incorporate prior knowledge without treating it as certainty.

For a practical introduction to priors, posteriors, hierarchical models, prediction and model validation, see Bayesian Sports Analytics with R: Predictive Modeling for Betting & Performance.

Automate Daily Sports Data Updates

Once your script works, you can schedule it to run every day. A simple pipeline might perform the following steps:

1. Download the latest games and statistics.
2. Request current sportsbook odds.
3. Save a timestamped odds snapshot.
4. Update team and player features.
5. Generate probabilities for upcoming games.
6. Compare model probabilities with market prices.
7. Save a report containing potential opportunities.

On Windows, you can automate an R script with Task Scheduler. On Linux or a server, you can use a cron job. GitHub Actions can also run scheduled workflows, although private API keys should always be stored as encrypted secrets.

Common Sports Betting Backtesting Mistakes

Using Information That Was Not Available Before the Game

Every model feature must represent information available at the time the bet would have been placed. Season averages calculated using games played after the prediction date create data leakage.

Ignoring Changes in the Betting Line

Opening odds, morning odds and closing odds are not interchangeable. Record the exact timestamp and price that your strategy uses.

Testing Too Many Strategies

If you test hundreds of filters, one strategy may appear profitable by chance. Use an out-of-sample period that was not used to select the strategy.

Using Accuracy as the Only Metric

A model can predict many winners correctly and still lose money if it consistently selects overpriced favorites. Calibration and expected value are more relevant than accuracy alone.

Assuming a Small Positive Return Proves an Edge

Sports betting returns are noisy. A strategy needs enough independent bets and should be evaluated with uncertainty intervals, drawdowns and sensitivity tests.

From Raw Data to a Complete Betting System

A complete sports betting workflow can be summarized as:

1. Collect performance data and bookmaker odds.
2. Clean team names, dates and market identifiers.
3. Create features using only past information.
4. Train a probabilistic model.
5. Evaluate calibration on unseen games.
6. Compare predictions with no-vig market probabilities.
7. Backtest realistic prices and betting rules.
8. Monitor results and update the model over time.

For readers who want to connect probabilities with expected value, the Kelly criterion and bankroll management, Bayesian Sports Betting with R: Probability, Kelly Criterion and Betting Strategies provides a focused guide to data-driven betting decisions in R.

Frequently Asked Questions

Can I get sports betting data for free in R?

Yes. Several R packages provide free sports performance data, and some odds providers offer limited free API access. Historical betting odds and frequent API requests may require a paid plan.

What is the best R package for sports betting odds?

oddsapiR is a convenient option for accessing The Odds API from R. You can also call the API directly with packages such as httr2 and process its JSON response with R.

Can I obtain NFL, NBA, MLB and NHL data with R?

Yes. The R sports analytics ecosystem includes packages such as nflreadr, hoopR, baseballr and fastRhockey.

How many years of data do I need?

There is no universal minimum. More seasons provide a larger sample, but older data may describe a different competitive or betting environment. Time weighting and rolling training windows can help balance sample size and relevance.

Can a sports betting model guarantee profits?

No. Predictive models estimate probabilities under uncertainty. They can be evaluated and improved, but they cannot eliminate variance, bookmaker margins, model error or financial risk.

Conclusion

R provides the tools needed to build a complete sports betting data pipeline: data collection, cleaning, feature engineering, probability estimation, backtesting and automated updates.

The most important step is not choosing the most complicated algorithm. It is creating a reliable dataset that preserves the information and odds actually available before each event. Once that foundation is correct, you can compare logistic regression, Elo ratings, Bayesian models, machine learning and simulation methods in a realistic way.

Start with one sport and one betting market. Save every odds snapshot, build a simple baseline and evaluate it on a future season before adding more complexity.

This article is for educational and analytical purposes only. Sports betting involves financial risk. No model or strategy can guarantee a profit.

The post How to Get Sports Betting Data in R: Free APIs, Historical Odds and Daily Updates appeared first on R Programming Books.

It’s 2:47pm, the meeting has been running for 38 minutes, and someone has just said, “Can everyone see my screen?”

You’re trying to concentrate, your laptop fan sounds like it’s preparing for take-off, and somewhere outside the temperature has climbed to 34°C. You’ve already relocated once today in search of a patch of shade under a tree with a half-decent breeze, laptop balanced on your knees, one eye on the battery icon. Meanwhile, your calendar still has three more meetings in it.

Summer 2026 has a funny way of making us notice things like this. It’s shaping up to be the UK’s warmest summer on record, with several heatwaves already behind us. Rather than spending another afternoon staring at a spreadsheet, why not give R something more interesting to do?

Let’s talk about the weather

This year, the UK weather hasn’t so much been unpredictable as relentlessly, record-breakingly hot. Which, if you’re a data person, is still a great excuse to explore.

For example, imagine we have daily high temperatures for London, Manchester and Newcastle across this summer’s headline heatwaves. We could use R to quickly compare them, find the warmest days, and create a chart showing just how much hotter the south has been running than the north.

library(dplyr)

# Illustrative daily highs (°C) across three of this summer's heatwave events
summer_weather <- tribble(
 ~date, ~city, ~temperature,
 "2026-05-24", "London", 32.3,
 "2026-05-24", "Manchester", 25.8,
 "2026-05-24", "Newcastle", 21.9,
 "2026-05-25", "London", 34.8,
 "2026-05-25", "Manchester", 26.9,
 "2026-05-25", "Newcastle", 22.6,
 "2026-05-26", "London", 35.1,
 "2026-05-26", "Manchester", 27.4,
 "2026-05-26", "Newcastle", 23.1,
 "2026-07-28", "London", 29.6,
 "2026-07-28", "Manchester", 26.2,
 "2026-07-28", "Newcastle", 22.4,
 "2026-07-29", "London", 34.2,
 "2026-07-29", "Manchester", 28.9,
 "2026-07-29", "Newcastle", 23.8,
 "2026-08-13", "London", 37.0,
 "2026-08-13", "Manchester", 32.1,
 "2026-08-13", "Newcastle", 25.6,
 "2026-08-14", "London", 38.1,
 "2026-08-14", "Manchester", 32.8,
 "2026-08-14", "Newcastle", 26.3
) |>
 mutate(date = as.Date(date))

library(ggplot2)

ggplot(summer_weather, aes(x = date, y = temperature, colour = city)) +
 geom_line() +
 geom_point() +
 labs(
 title = "How warm has UK summer 2026 been?",
 subtitle = "Daily highs across three of this summer's headline heatwaves",
 x = NULL,
 y = "Temperature (°C)"
 ) +
 theme_minimal()

Suddenly, that spreadsheet of numbers becomes something you can actually explore. You can see, at a glance, just how much bigger the north-south gap gets once a heatwave really takes hold – London hit 38.1°C in mid-August while Newcastle stayed at 26.3°C on the same day.

And you don’t have to stop at temperature. You could look at rainfall (or the lack of it), reservoir levels, ice cream sales, train delays… or even investigate whether your team’s productivity mysteriously drops once the office hits 25°C.

A little summer upskilling

This is also where R becomes particularly useful.

Learning R isn’t just about knowing how to write code. It’s about becoming more comfortable taking a question, finding the right data, exploring it and turning the results into something that other people can understand.

Maybe you’ve been using R for a while but keep thinking, “There must be a better way to do this.”

There probably is.

Maybe you’re comfortable with the basics but want to get better at data visualisation, modelling, reproducible reporting or working with larger datasets.

That’s exactly where structured training can help.

At Jumping Rivers, we run practical training across R, Python, SQL, Quarto, Shiny and the wider Posit ecosystem. Our courses are designed around real-world work, so you can take what you learn straight back to your desk — preferably somewhere with a fan, or this year, maybe even air conditioning.

Make the summer count

The summer break can be a good time to step back from the usual routine and invest in skills that make the rest of the year easier.

You could spend a few hours learning a better way to manipulate data, finally get to grips with ggplot2, explore Quarto, or start building your first Shiny application.

And if you’re not sure what training would actually be useful for you or your team, that’s where we can help too.

We can look at the tools you’re currently using, the skills your team already has and where the gaps are, then recommend a training path that makes sense.

No complicated spreadsheets required.

Although, if you do have one – perhaps tracking this summer’s heatwaves – we’re always happy to help you turn it into something much more interesting.

Ready to make your next data project a little more enjoyable?

Explore our public training courses or get in touch with the Jumping Rivers team to talk about what would work best for you.

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

I continue to complete work on my next post on using the R wordcloud package.  As I normally do programming and wrting with my Lenovo desktop computer, I decided to experiment with installing R-base and RStudio on my Raspberry Pi Model 3B and tablet computer for those occasions when I desire to work while traveling.    The Raspberry Pi is running the latest version of Debian Trixie along with the Raspberry Pi desktop. The tablet is using the latest version of Android 16. My first observation relates to the availability of R for these platforms.

I. Installing R-base and RStudio on Raspberry Pi3/B,4,5

R is available for a variety of UNIX, Windows, and MacOS systems. If you are running R on Windows, you are familiar with the 32- and 64-bit versions available for download and installation via an executable loader. While the R-base console has been available for the RPi platform, it has only been recently that the RStudio-server has been available for the ARM64 processor used in the RPi 3A/B, 4, and 5 models. I am currently using Debian Trixie 64-bit on my RPi 3/B. The R-base package is now available in Debian repositories so R can be installed via the RPi desktop menu rather than downloading binary builds or executable files.

For an RPi 3B or higher I would recommend the following:
-Make sure your microSD card is large enough. I am using a 32 GB card.
-Make sure your OS is up to date. Use the command line utility to run the following commands:
       sudo apt update (respond to prompts that follow)
       sudo apt full-upgrade (respond the prompts)
Depending on the model of RPi you are using, the memory card size and Debian version you are using this update could take quite some time.

Once the update is completed use the desktop menu to access the add/delete software option, search for the R package using r-base as the keyword and click on the appropriate icon to start the installation. When the installation is complete you should see the R icon in the desktop dropdown menu under the Programming or Science (or both) headings. The R-base console can now be run by clicking the menu icon, and R is now available for access by RStudio-server if it is installed.

Because the RPi uses an ARM processor, RStudio itself cannot be installed, but RStudio-server has been successfully ported to the platform. Additional information on R downloads can be found at the Posit web site RStudio IDE User Guide RStudio User Guide, and at the link RStudio Latest Builds. If you wish to install the RStudio-server from your RPi command line utility there are several steps, but the result is a working web-based interface with full RStudio-server build. Follow the steps listed below.
1. When installing new software run an update using:
     sudo apt update
2. The port of RStudio-server we are installing was designed for the Ubuntu OS so install dependencies needed for Debian using:
     sudo apt install gdebi-core libssl-dev libclang-dev
3. Get the build from the Posit Daily Builds library using 

wget https://dl.dailies.rstudio.com/server/jammy/arm64/rstudio-       server-2026.06.0-242-arm64.deb
4. Install the application using:
 sudo gdebi rstudio-server-2026.06.0-242-arm64.deb
5. When the installation is complete use the system service command to start RStudio-server with:
     sudo systemctl start rstudio-server (for the current  bootup)                                             and/or
     sudo systemctl enable rstudio-server (start at all bootups)
6. Open the Chromium or Firefox web browser from the desktop menu and access the RStudio-server by entering the URL:
      http:// IP address on your network>:8787
In my case I would enter http://192.168.4.115:8787
The screenshot shown below shows RStudio with the code from this article and the resulting output.

II. Using R on an Android Device

R and RStudio will not port directly to an Android based OS, but there are a few applications that will work with varying degrees of utility. I have a tablet that runs Android 16 and am using a free application, Rlytic,. Once installed from the Play Store users sign up with a username and password. When the program starts, a code entry console is displayed. Your code can be entered directly using the on-screen keyboard provided or can be loaded from your device file storage or cloud storage. The interface is easy to use. I have included a simple program example and some screenshots below.
Rlytic is free to use but is restricted to having only 2 programs active at a time. An unlimited version is available for purchase. I might also add that at the time of this writing Rlytic is running on R-base v.3 so users may run into some problems with more complex projects.

III. Using R on an iPhone

I currently use an iPhone 12 and was curious about any R applications that would work with it. I found an application called WebR which combines R-base 4.xxx with a text editor and browser interface. According to the program s author the application was designed for use by students in a classroom setting when learning statistics and/or R programming. It provides a highly mobile platform for Running R programs and quickly generates both text and graphics output. Once again, I will leave it to readers to engage the application s learning curve and will provide a simple example and screenshots below. The software is free and is available in the iPhone App Store.

IV. Sample Program: Raspberry Pi

The following code is a simple example of how R can be used to demonstrate the Central Limit Theorem in sampling from a population. The code uses the R-base rnorm function to generate randomly selected samples from a normally distributed population of values with a given population mean and standard deviation, finds the mean of each sample generated and graphs the sampling distribution. The code is shown below.

#population; sd=10; mean=65
#generate 25 samples of 25 observations
#calculate sample mean of each sample and plot distribution
###################################################
#code to generate samples and display all sample means
###################################################
Samples <- replicate(25, rnorm(25, mean=65, sd=10))
Samples #show the samples generated
##################################################
#code to calculate and display mean of each column of sample means
#################################################
SampleMeans <- colMeans(Samples)
SampleMeans #show the means of the samples generated
####################################################
#code to plot means of the sampling distribution
#####################################################
plot(density(SampleMeans),
main = “Density of Sample Means”,
xlab = “Sample Mean”)
The plot of the distribution of the sample means is shown below.

V. Sample Program: Rlytic
Here is the same code with the plot of the results for the Rlytic app on my Android 16 tablet. For brevity I have not included all the hashtag dialog from the RPi example. The screenshot and plot are shown.
#population; sd=10; mean=65
#generate 25 samples of 25 observations
#calculate sample mean of each sample and plot distribution
Samples <- replicate(25, rnorm(25, mean=65, sd=10))
#Samples
SampleMeans <- colMeans(Samples)
#SampleMeans
plot(density(SampleMeans),
main = “Density”,
xlab = Mean”)
The Rlytic Screen:

The Rlytic Plot: (Note Rlytic graphs are PDF format)

VI. Sample Program: WebR for iPhone
Here is a slightly modified version of the random sampling code entered into WebR on my iPhone 12.
x=rnorm(25, mean=65, sd=10)
plot(density(x))
Shown below is the resulting output. As in previous examples I did not print the output showing the randomly generated individual means.

I am still working on the next part of my tutorial on using wordcloud and related packages for the analysis of large, complex text files. Please look for my next post in the not-too-distant future.
D.M. Wiig
R Statistics and Programming
https://dmwiig.net

Deep learning for time series forecasting in R has historically faced a major architectural hurdle: Python overhead. Frameworks like modeltime.gluonts provide interface wrappers around AWS GluonTS, but they rely on a complex execution chain passing through reticulate, virtual environments, Python serialization, and MXNet/PyTorch backends.

To overcome the performance bottlenecks and dependency friction of cross-language bridging, we engineered a pure R-native DeepAR forecasting engine. Powered by the C++ LibTorch backend via R’s torch package, this architecture offers lightweight, in-memory execution without any Python or reticulate dependencies.

Architectural Comparison: Modeltime/GluonTS vs. Native R Torch

The architectural difference between traditional wrappers and our native C++ LibTorch binding lies in data marshalling and execution depth:

Deep Dive into the Code Architecture

Our R implementation mirrors the probabilistic depth of DeepAR while maintaining computational stability and clean visual interactivity.

1. Bounded Student-t Distribution Head

Financial time series, such as the SOXX ETF, exhibit heavy-tailed return distributions (“fat tails”) and sudden volatility shocks. Gaussian models often understate extreme risks or produce over-reactive prediction bands.

We implement a 3-head architecture off the LSTM hidden state:

• Location parameter (μ): Unconstrained linear output layer.
• Scale parameter (σ): Softplus activation layer with numerical stability offset.
• Degrees of freedom parameter (ν): Bounded dynamically between 4.0 and 30.0 using a scaled sigmoid:

Bounding ν ≥ 4.0 guarantees mathematically finite variance, preventing Monte Carlo variance explosion over multi-step autoregressive horizons.

2. Variance-Controlled Stochastic Monte Carlo Sampling

During the 10-day forecast horizon, we generate 100 autoregressive simulation paths. To align Monte Carlo variance with predicted σ, we scale Student-t samples by the theoretical variance factor:

$\text{scale\_factor} = \sqrt{\frac{\nu – 2}{\nu}}$

$y_{t} = \mu + \sigma \cdot \text{scale\_factor} \cdot t_{\nu}$

This ensures the trajectory bounds remain stable across multi-step autoregressive rollouts.

3. Granular Interactive Plotly Visualization

The frontend layer leverages ggplot2, ggtext, and plotly to deliver clean UI/UX interactivity:

• Embedded HTML Titles: Eliminates redundant legend boxes by color-coding series names directly inside the Markdown title using ggtext::element_markdown.
• Invisible Boundary Anchors: Invisible hover points (alpha = 0) are placed along the 95% confidence bounds (conf_hi and conf_lo). Users can inspect exact upper/lower boundary prices dynamically without cluttering the plot with extra lines.

Complete R Script

# ==============================================================================
# TITLE: Pure R-Native Torch DeepAR - Bounded Student-t Distribution Engine
# PATH: tool_nodes/forecasting/engine/evaluate_torch_deepar_student_t_bounded.R
# DEPLOYMENT TARGET: Native R Pipeline (Zero Python / Zero Reticulate Dependency)
# All code descriptions and labels are systematically maintained in English.
# ==============================================================================

if (!require("pacman")) install.packages("pacman")
pacman::p_load(tidyquant, tidyverse, timetk, torch, plotly, yardstick)

# 1. Fetch & Prepare Data from Yahoo Finance
df_dl <- tq_get("SOXX") %>%
  select(date, close) %>%
  filter(date >= last(date) - months(12)) %>%
  drop_na()

# Data Normalization Parameters
mean_close <- mean(df_dl$close)
sd_close   <- sd(df_dl$close)
df_dl      <- df_dl %>% mutate(close_scaled = (close - mean_close) / sd_close)

# Configuration Parameters
lookback_length   <- 20
prediction_length <- 10
num_paths         <- 100

train_data <- head(df_dl, nrow(df_dl) - prediction_length)
test_data  <- tail(df_dl, prediction_length)

# 2. Sequence Generator
create_sequences <- function(data_vector, lookback) {
  num_samples <- length(data_vector) - lookback
  x_mat <- matrix(0, nrow = num_samples, ncol = lookback)
  y_mat <- matrix(0, nrow = num_samples, ncol = 1)
  
  for (i in 1:num_samples) {
    x_mat[i, ]  <- data_vector[i:(i + lookback - 1)]
    y_mat[i, 1] <- data_vector[i + lookback]
  }
  
  list(
    x = torch_tensor(x_mat, dtype = torch_float())$unsqueeze(3),
    y = torch_tensor(y_mat, dtype = torch_float())
  )
}

seqs <- create_sequences(train_data$close_scaled, lookback_length)

# 3. Native Torch DeepAR Architecture with Bounded Student-t Head
deepar_student_net <- nn_module(
  "DeepARStudentNetBounded",
  initialize = function(input_size = 1, hidden_size = 32, num_layers = 2) {
    self$lstm     <- nn_lstm(input_size = input_size, hidden_size = hidden_size, 
                             num_layers = num_layers, batch_first = TRUE)
    self$fc_mu    <- nn_linear(hidden_size, 1)
    self$fc_sigma <- nn_linear(hidden_size, 1)
    self$fc_v     <- nn_linear(hidden_size, 1)
  },
  forward = function(x) {
    out <- self$lstm(x)
    last_hidden <- out[[1]][, dim(out[[1]])[2], ]
    
    mu    <- self$fc_mu(last_hidden)
    sigma <- nnf_softplus(self$fc_sigma(last_hidden)) + 1e-4
    
    # Bound degrees of freedom v between 4.0 and 30.0 to prevent explosive tails
    v     <- 4.0 + 26.0 * torch_sigmoid(self$fc_v(last_hidden))
    
    list(mu = mu, sigma = sigma, v = v)
  }
)

model     <- deepar_student_net()
optimizer <- optim_adam(model$parameters, lr = 0.003)

# Stable Student-t Negative Log-Likelihood Loss
student_t_nll_loss <- function(mu, sigma, v, y) {
  term1 <- torch_lgamma((v + 1) / 2)
  term2 <- torch_lgamma(v / 2)
  term3 <- 0.5 * torch_log(v * pi)
  term4 <- torch_log(sigma)
  
  residual <- (y - mu) / sigma
  term5 <- ((v + 1) / 2) * torch_log(1 + (residual$pow(2) / v))
  
  - (term1 - term2 - term3 - term4 - term5)
}

# 4. Training Loop
model$train()
for (epoch in 1:40) {
  optimizer$zero_grad()
  preds <- model(seqs$x)
  loss  <- student_t_nll_loss(preds$mu, preds$sigma, preds$v, seqs$y)$mean()
  loss$backward()
  
  # Gradient clipping for numerical stability
  nn_utils_clip_grad_norm_(model$parameters, max_norm = 1.0)
  optimizer$step()
}

# 5. Stochastic Monte Carlo Trajectory Sampling (Variance Variance-Controlled)
model$eval()
price_paths <- matrix(0, nrow = num_paths, ncol = prediction_length)
initial_input_seq <- tail(train_data$close_scaled, lookback_length)

with_no_grad({
  for (s in 1:num_paths) {
    curr_seq <- initial_input_seq
    
    for (t in 1:prediction_length) {
      curr_tensor <- torch_tensor(matrix(curr_seq, nrow = 1), dtype = torch_float())$unsqueeze(3)
      pred <- model(curr_tensor)
      
      mu    <- as.numeric(pred$mu)
      sigma <- as.numeric(pred$sigma)
      v_val <- as.numeric(pred$v)
      
      # Scaled Student-t sampling to strictly align variance with sigma
      scale_factor <- sqrt((v_val - 2) / v_val)
      sampled_scaled <- mu + sigma * scale_factor * rt(1, df = v_val)
      
      price_paths[s, t] <- sampled_scaled * sd_close + mean_close
      
      # Autoregressive slide
      curr_seq <- c(curr_seq[-1], sampled_scaled)
    }
  }
})

# 6. Extract Quantiles & Prepare Tidy Evaluation Data Frame
predicted_prices <- colMeans(price_paths)
lower_bound      <- apply(price_paths, 2, quantile, probs = 0.025)
upper_bound      <- apply(price_paths, 2, quantile, probs = 0.975)

df_eval <- tibble(
  date     = test_data$date,
  actual   = test_data$close,
  pred     = predicted_prices,
  conf_lo  = lower_bound,
  conf_hi  = upper_bound
)

# 7. Tidymodels / Yardstick Metric Engine
eval_metrics <- metric_set(mape, rmse, rsq)

metrics_summary <- df_eval %>%
  eval_metrics(truth = actual, estimate = pred) %>%
  select(.metric, .estimate) %>%
  rename(Metric = .metric, Value = .estimate)

print(metrics_summary)

mape_val <- metrics_summary %>% 
  filter(Metric == "mape") %>% 
  pull(Value)


# 8. Modern Interactive Plotly Visualization (Clean Lines & Clear Ribbon)

if (!require("pacman")) install.packages("pacman")
pacman::p_load(tidyquant, tidyverse, plotly, scales, glue, ggtext)

# 1. Prepare Dedicated Hover Text Layers
df_plot_actual <- df_eval %>% 
  select(date, actual) %>% 
  mutate(text_actual = glue::glue("<b>Actual Price:</b> ${round(actual, 2)}\n<b>Date:</b> {format(date, '%b %d, %Y')}"))

df_plot_pred <- df_eval %>% 
  select(date, pred) %>% 
  mutate(text_pred = glue::glue("<b>DeepAR Pred:</b> ${round(pred, 2)}\n<b>Date:</b> {format(date, '%b %d, %Y')}"))

df_plot_hi <- df_eval %>% 
  select(date, conf_hi) %>% 
  mutate(text_hi = glue::glue("<b>95% Upper Bound:</b> ${round(conf_hi, 2)}\n<b>Date:</b> {format(date, '%b %d, %Y')}"))

df_plot_lo <- df_eval %>% 
  select(date, conf_lo) %>% 
  mutate(text_lo = glue::glue("<b>95% Lower Bound:</b> ${round(conf_lo, 2)}\n<b>Date:</b> {format(date, '%b %d, %Y')}"))

# 2. Build GGPlot Spec with Invisible Boundary Anchors
p <- ggplot() +
  # Clean Background Ribbon
  geom_ribbon(
    data = df_eval,
    aes(x = date, ymin = conf_lo, ymax = conf_hi),
    fill  = "#808080",
    alpha = 0.20
  ) +
  # Invisible Upper Bound Hover Points (No Lines, Pure Hover)
  geom_point(
    data = df_plot_hi,
    aes(x = date, y = conf_hi, text = text_hi),
    color = "transparent",
    alpha = 0,
    size  = 3
  ) +
  # Invisible Lower Bound Hover Points (No Lines, Pure Hover)
  geom_point(
    data = df_plot_lo,
    aes(x = date, y = conf_lo, text = text_lo),
    color = "transparent",
    alpha = 0,
    size  = 3
  ) +
  # Actual Price: Solid Dark Line & Hover Points
  geom_line(
    data = df_plot_actual,
    aes(x = date, y = actual),
    color = "#2c3e50",
    linewidth = 1.2
  ) +
  geom_point(
    data = df_plot_actual,
    aes(x = date, y = actual, text = text_actual),
    color = "#2c3e50",
    size  = 2
  ) +
  # DeepAR Forecast: Dashed Red Line & Clean Hover Points
  geom_line(
    data = df_plot_pred,
    aes(x = date, y = pred),
    color = "#e74c3c",
    linetype = "dashed",
    linewidth = 1.2
  ) +
  geom_point(
    data = df_plot_pred,
    aes(x = date, y = pred, text = text_pred),
    color = "#e74c3c",
    size  = 2
  ) +
  # Formatting & Theme
  scale_y_continuous(labels = dollar_format(accuracy = 1)) +
  labs(
    x = "",
    y = "",
    title = paste0(
      "SOXX ETF <span style = 'color:#2c3e50'>Actual Prices</span> vs ",
      "<span style = 'color:#e74c3c'>Torch DeepAR Forecast</span><br>",
      "<span style='font-size:12px; color:#555555;'>10-Day Horizon | MAPE: ", round(mape_val, 2), "%</span>"
    )
  ) +
  theme_minimal() +
  theme(
    plot.title = element_markdown(
      hjust = 0.5, 
      face  = "bold"
    ),
    plot.background  = element_rect(fill = "#ffffff", color = NA),
    panel.background = element_rect(fill = "#ffffff", color = NA),
    panel.grid.minor = element_blank()
  )

# 3. Render Interactive Plotly Spec
font_family <- list(family = "Roboto Slab, Sans-Serif", size = 16)
label_font  <- list(font = list(family = "Roboto Slab, Sans-Serif", size = 13))

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

Conclusion

By implementing DeepAR directly in R via torch (LibTorch), we achieve a low-latency, zero-Python architecture that fits naturally into existing tidymodels workflows. The resulting pipeline delivers high-precision probabilistic predictions (achieving a MAPE of ~2.17% on a 10-day SOXX forecast horizon) with fast, in-memory performance suitable for production deployment.

Over the course of my life I’ve invested a lot of time and effort into the acquisition of different skills, abilities, competencies, specialized knowledge – for simplicity I will refer to all of it as “skills” in the following. Some of those skills, I would say, have even become an integral part of my identity, which isn’t unusual for a human, I guess. At my current age (sadly I’m not 17 anymore) I have accumulated a skill collection sizable enough that it needs to be managed systematically.

I need to decide: What do I want maintain? Maintain fully or partially, or merely keep-warm? Where do I want to improve? What’s the minimal maintenance dose required in each case? What new skills will I invest substantial efforts in? Which of my current skills do I choose to neglect (to balance everything out)?

Below are some directional thoughts on skill maintenance/improvement/neglect with some vague personal examples.

1. Skill maintenance generally requires far less time and/or effort than skill acquisition (or meaningful improvement). Therefore it makes sense to keep a skill sharp at a minimally acceptable level rather than letting it deteriorate. Of course, there are different levels of “maintain”, requiring different levels of effort. For instance, for now I merely keep-warm my music skills by playing my instruments only a few minutes per week.
2. Like the approach to acquiring a skill will differ between types of skills, such as motor skills, procedural skills, knowledge-based skills, etc., so will the approach to maintenance differ too. A few examples:
   • It seems that a memory-reliant skill will require active recall, e.g., spaced repetition (which I’ve been struggling with for almost 10 years now because it’s actually hard work to do the reviews). But, almost all non-trivial skills rely on memory to a significant degree.
   • Motor skill development and (later) maintenance appears to be based on regular practice of basically the same fundamental movements; for example, practicing scales on a musical instrument, or practicing the same kicks/punches over and over in a martial art.
   • Language skills rely primarily on regular exposure and use. Living in the US, this is something I need to actively/proactively seek out for my non-English languages. Reading books seems to work to some extent to maintain a language I already know, though it isn’t enough – at least it gives me an excuse to read fun but mediocre fiction (as in “I only read this to maintain my German/Russian!”) which I wouldn’t be reading otherwise.
3. A skill is made up of many component sub-skills that deteriorate at different rates. For example, fundamental hand and finger movements of playing an instrument are retained much better than individual music pieces or specific memorized chord progressions. Sometimes I may choose to maintain a specific sub-skill, and hope the best for the rest of the overall skill. For example, with the current AI wave, I chose to maintain only my competence in reviewing, rather than writing, code in certain my-non-main programming languages.
4. Despite point 1 above, many skills are easily restored if forgotten. There is not enough time in the day/week/month to dedicate to everything, and I need to make a conscious decision on what to neglect, eventually allowing it to atrophy. However, often skills come back quickly with a little use, even if they deteriorate quickly. A silly example: I lost the ability to type in QWERTY immediately after learning the Colemak keyboard layout; then, years later, easily relearned QWERTY in 1-2 weeks, but immediately lost the ability to type in Colemak; after using QWERTY exclusively for a year, I then easily relearned Colemak and have been maintaining both layouts in muscle memory ever since.

That’s the end of this blog post. So, obviously, a skill I should work on is writing.

Introduction

In previous posts I have described how to use the self-contained Python scripts feature in the uv Python package manager to create virtual environments to render Quarto documents using the Jupyter nbstata kernel and the python3 kernel. In this post I describe how to do the same for R scripts to render Quarto documents running R code using the knitr engine.

I recently discovered that there are now three uv-inspired package managers for R; ir, uvr, and rv (… maybe there are more?). I will concentrate on the first two because they allow defining self-contained R scripts. I find self-contained scripts a fast and lightweight way to define project dependencies, and I very rarely require a record of the exact package versions.

In the following examples I assume we are creating an R script, render.R, which contains one or more calls to quarto::quarto_render() for a lecture or tutorial. For the dependency R packages I include the packages the document itself needs, plus the quarto and knitr packages.

Example self-contained R script using ir

To define dependencies for ir, at the top of the script begin each comment line with #| then write a list under a packages key as follows – this is the list of packages I require for one of my practicals on missing data.

#| packages:
#|   - gtsummary
#|   - haven
#|   - tidyverse
#|   - VIM
#|   - quarto
#|   - knitr

# Rest of R code follows ...
# ... essentially one or sometimes multiple quarto::quarto_render() calls

This script can be run with

ir run render.R

Example self-contained R script using uvr

uvr follows the same dependency syntax as uv. Each line begins with a # comment, and the dependencies are defined as a TOML array of strings between # /// script and # ///. So the top of our render.R script looks as follows.

# /// script
# dependencies = [
#   "gtsummary",
#   "haven",
#   "tidyverse",
#   "VIM",
#   "quarto",
#   "knitr",
# ]
# ///

# Rest of R code follows ...
# ... essentially one or sometimes multiple quarto::quarto_render() calls

This script can be run with

uvr run render.R

Automation with just in a complex directory structure

For each course I teach I have the lecture or tutorial in a subdirectory. To run each script I could run the shell commands given above. To slightly improve efficiency I find that putting the following justfile at the top of the directory structure saves a bit of typing. The first recipe, render, uses my system R library, the others resolve packages via ir/uvr.

render dir=invocation_directory():
    cd "{{ dir }}" && Rscript render.R

ir dir=invocation_directory():
    cd "{{ dir }}" && ir run render.R

uvr dir=invocation_directory():
    cd "{{ dir }}" && uvr run render.R

I can simply type just ir or just uvr to render the lecture/tutorial given whichever directory I’m in.

Bonus 1 – Example self-contained Quarto document using ir

ir cleverly allows us to alternatively define the dependencies within the YAML header of a Quarto document, under an ir key. In this case we can remove the quarto package as we might assume we’d render this document by clicking the Render button in RStudio or using quarto render ... in the terminal.

---
title: My lecture/tutorial
ir:
  packages:
    - gtsummary
    - haven
    - tidyverse
    - VIM
    - knitr
---

Rest of Quarto document follows ...

Say this Quarto document is tutorial.qmd we would then render it with

ir render tutorial.qmd

More details are given in the ir Quarto docs.

Bonus 2 – Making the R script executable

With both ir and uvr (and indeed uv) we can optionally make the render.R script executable, say renaming to simply render, by adding the relevant shebang to the very top of the file.

For ir we add

#!/usr/bin/env -S ir run

and for uvr we add

#!/usr/bin/env -S uvr run

We then make the script executable

chmod +x render

and run it with

./render

Summary

I have shown how to make a self-contained, and optionally executable, R script to render Quarto documents using the knitr engine which automatically manages the required R packages. This functionality is provided by both the ir and uvr R package managers. This approach would also work for RMarkdown documents (of course one would need to swap the quarto package for the rmarkdown package in the list of dependencies).

Before the main content: I am creating an R Community on Google Groups. You can join the group using this form.

I’ve been working on a a modern dashboard framework for R using the beautiful Tabler Bootstrap theme. Furthermore, to render Tabler apps using a server I created Tabler Server alongside the process.

Installation

Old version, depends on Shiny:

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

New version, does not use Shiny:

# using the R-Universe
install.packages("tabler", repos = "https://pachadotdev.r-universe.dev")

# or using the remotes package
remotes::install_github("pachadotdev/tabler")

Quick Start

Single-script app

The following example uses the “combo” layout to recreate Shiny’s geyser example. The theme options can be adjusted from the code or the theme setting icon that can be hidden. See the example here.

I added a UI-only example to cover the different input elements and their options here.

Modular R package app

Create an R package with modular components:

library(tabler)

pkg_template("mydashboard")

See the package skeleton here. pkg_template() adds a DESCRIPTION and other components required for an R package to work.

For instance, Open Trade Statistics consists in a full dashboard that uses environment variables, SQL connections, caching, and D3 plots. Its code is here, and the result is here.

Loading/Progress bar

I added an example with a progress bar here. The progress bar hides the app while the new plots or other elements are computed.

Login page

This R package provides a login page that you can connect to a database or another system. The example here shows the dashboard after correctly typing the user “SpaceMariner” and password “IDDQD”. There is an example using RSQLite here.

I was thinking about adding a Google/Outlook/GitHub account login but I have no idea how to. If you know how and would like to contribute, please comment here.

Does it run Doom?

Yes. I tested the WebSocket compiling and running the Doom WASM version. See the code here.

Available Layouts

There are additional examples for each of the following layouts:

• Boxed (Default): Basic dashboard with top navbar and constrained width content area. This is the default layout.
• Combo: Combines vertical sidebar navigation with top header.
• Condensed: Compact layout with reduced padding/margins.
• Fluid: Full-width layout without container constraints.
• Fluid Vertical: Full-width layout with vertical sidebar.
• Horizontal: Layout with horizontal navigation menu.
• Navbar Dark: Layout with dark navbar theme.
• Navbar Overlap: Layout where content overlaps with navbar for a modern look.
• Navbar Sticky: Layout with sticky/fixed navbar that stays at the top when scrolling.
• RTL: Right-to-left layout for Hebrew/Arabic languages.
• Vertical: Vertical sidebar layout without top navbar.
• Vertical Right: Vertical sidebar positioned on the right side.
• Vertical Transparent: Vertical layout with transparent sidebar.

Note: tabler allows to pass layout = "navbar" and layout = "navbar-sticky-dark" which are wrappers for a light theme navbar layout and a dark theme sticky navbar layour, respectively.

Differences with Shiny

• Static plots (base, ggplot, tinyplot, etc.) render as SVG and can be downloaded with the right click button.
• URLs are of the form my.site/myapp?year=2000&country=gbr instead of my.site/myapp?year=2000&country=%22gbr%22

License

Apache License (>= 2)

Usually I don’t watch much youtube (a social phenomenon I have somewhat missed out on), but here’s an interesting video about another social phenomenon I have somewhat experienced myself .

With LLM coding, I feel like I can do anything and start working on too many things simultaneously / in multiple parallel sessions (pyfixest features and refactoring, a paper on our new demeaning algo, R bindings for the new algo, presentations on this and that, blog posts, a JOSS paper draft, etc etc).

I am paying for subscriptions with Codex, Claude, and OpenCode, and when I don’t spend all of my usage, I feel like I missed an opportunity, so sometimes I skip going to the gym to keep on working because I still have usage to spend or usage will reset in 30 minutes so I can keep on working prompting my agents to keep on working… This was particularly bad the weekend that Fable was temporarily available on the 20$ plan – I almost couldn’t leave my keyboard, the desire to keep on prompting was too strong. Needless to say that – because I prompted without much thinking – the results I got were not very good.

I usually start working on OSS after work around 21:00, and suddenly it’s 23:30 and I want to stop, but there is one last prompt to write and finally I stop and it is 00:30 but I feel no accomplishment because I haven’t been doing any deep work … but at least, I used all of my usage (#tokenmaxxing)!

Long story short, I think a more optimal strategy for me going forward is to reduce the number of subscriptions (fewer token limits to hit) and topics I work on in parallel and, to some degree, stop over-relying on an auto-mode-and-review workflow and switch back to a more “socratic” style of coding with LLMs.

When dockerizing an R Shiny application hosted via Shiny Server (built from shiny server docker image https://hub.docker.com/r/rocker/shiny), a common issue developers face is that environment variables set via ENV instructions in the Dockerfile (or passed at runtime via docker run -e) are completely missing inside the R Shiny app.

For example, if you set the following environment variable in your Dockerfile:

1

ENV SHOULD_IN_SHINY="From dockerfile"

When your app launches and runs Sys.getenv("SHOULD_IN_SHINY"), it returns an empty string "" instead of "From dockerfile". This is not a problem with Docker itself, but rather a consequence of how Shiny Server spawns R worker processes.

In this post, we will look into the root cause behind this behavior in Shiny Server and demonstrate the two recommended solutions to correctly expose environment variables to your R Shiny workers.

The Root Cause: How Shiny Server Spawns R Processes

The reason environment variables do not carry over to your R session lies in how Shiny Server executes R worker processes inside the container.

Shiny Server runs as a system service (typically as root or shiny). When launching an app instance, it re-executes R as the shiny unprivileged user using su. Specifically, the execution call combines two mutually exclusive flags:

1

su shiny --login --preserve-environment -c "... R ..."

Let’s break down what these flags request from su:

1. --login (-l): Starts a login shell. This intentionally resets the environment to a minimal whitelist (HOME, PATH, USER, TERM, etc.) and sources system startup profiles (/etc/profile, ~/.profile).
2. --preserve-environment (-p): Explicitly asks su to keep the current environment inherited from the caller (which includes Docker’s ENV variables).

Because Linux’s su utility cannot honor both conflicting behaviors, it chooses --login and prints a warning:

1

su: ignoring --preserve-environment, it's mutually exclusive with --login

Note: This warning is emitted by Shiny Server’s underlying process execution call, not by the R app itself. Attempting to suppress or patch it requires modifying and rebuilding Shiny Server C++ code, which is rarely practical.

Because the login shell wins, all custom environment variables passed to the container via ENV or docker run are wiped before R ever starts.

Solutions

Since patching Shiny Server is unnecessary, we can utilize the natural extension points provided by the login shell or R itself.

Solution 1: Use ~/.profile (Recommended for User-Level Shell Vars)

Since --login causes the shell to source ~/.profile for the shiny user, we can write our environment variables to /home/shiny/.profile during the Docker build phase.

In your Dockerfile:

1
2

# Append environment variable to the shiny user's profile
RUN echo 'export SHOULD_IN_SHINY="From dockerfile"' >> /home/shiny/.profile

When Shiny Server executes su shiny --login ..., the login shell will read /home/shiny/.profile and load SHOULD_IN_SHINY into the environment right before starting the R process.

Solution 2: Use ~/.Renviron (Recommended for R-Specific Configs)

Alternatively, R automatically inspects and loads ~/.Renviron on startup, right after the shell environment is initialized.

~/.Renviron is purpose-built for R:

• Scoped strictly to R processes.
• Does not use shell export keywords or $ variable expansions.
• Uses simple KEY=value key-value pairs.

In your Dockerfile, set up /home/shiny/.Renviron:

1
2
3

# Set up .Renviron for the shiny user
RUN echo 'SHOULD_IN_SHINY="From dockerfile"' >> /home/shiny/.Renviron \
    && chown shiny:shiny /home/shiny/.Renviron

In your Shiny app (app.R or server.R), you can access it reliably:

1
2

should_in_shiny <- Sys.getenv("SHOULD_IN_SHINY")
# Returns: "From dockerfile"

Summary

• The Problem: Shiny Server calls su shiny --login --preserve-environment, causing su to discard inherited environment variables (such as Docker ENV) in favor of a clean login shell.
• Solution 1 (~/.profile): Append export KEY="value" to /home/shiny/.profile in your Dockerfile.
• Solution 2 (~/.Renviron): Append KEY="value" to /home/shiny/.Renviron in your Dockerfile.

Both approaches integrate seamlessly with Docker builds and ensure your Shiny application receives all required configuration variables cleanly.

As for last year, I’ve got a poster accepted at COPA conference.

This year, the poster is about ‘Zero-Shot Probabilistic Stock Returns Forecasting with Pretrained RVFL Networks’.

Links to the paper and the code repository are below:

• Paper: https://www.researchgate.net/publication/412115839_Zero-Shot_Probabilistic_Stock_Returns_Forecasting_with_Pretrained_RVFL_Networks
• R Code: https://github.com/thierrymoudiki/2026_05_28_Pretrain_Ridge2_Stocks_Full_Pipeline

A few days ago I had a coffee in Paris with my long-time friend (and former Statistics & Computing editor) Gilles Celeux, and he mentioned me stopping solving and posting maths puzzles like those weekly published by Le Monde. They have indeed vanished with the retirement of the authors, but Gilles added that the arrival of LLMs would have made the exercise moot. I disagreed as (i) the fun of solving the puzzle on my own  has not gone away and (ii) the pedagogical appeal of the puzzle and its resolution remains. As the next Fiddler puzzle arrived in my mailbox, my resolution was put to the test (contrariwise to the previous entry, which did not require massive computations):

The Fiddler League consists of two teams. Over a season, they play each other 162 times. Each team has an equal chance of winning each game, and the results of games are independent. Over the season, on average, how many games would you expect the team with the better record to have won?

I started on the wrong foot with E[X|X≥81] when X is Bin(162,½), equal to 85.77 (either directly or with a Normal approximation), which differs from my second thought, E[max(X,162-X)]=86.07 (either directly or with a Normal approximation), which is larger because of the reflection produced by max. While the first computation was manageable, the second one seemed to involve simulation and I caved in prompting Claude, which provided the answer along with the connection

E[max(X,162−X)]=E[X∣X≥81](1+p₈₁​)−81p₈₁​

After some expansion, the League boasts 30 teams. Over a season, each team plays each other team five times. (Each team plays a total of 145 games.) Again, each team has an equal chance of winning each game, and the results of games are independent. Over the season, on average, how many games would you expect the team with the best record to have won?

The best record is max(X_i) with each of the 30 X_i‘s a sum of 29 Y_ij and the Y_ij=5-Y_ji distributed as Bin(5,½). The X_i‘s are thus Bin (145,½) but dependent. While I could not figure out a closed form answer for the expectation, a direct Monte Carlo resolution is obviously feasible, with Claude (rather than me) running it over 400 million repetitions, but a 30 dimensional Normal approximation exploiting the correlation of 1/29 between the components leads to roughly 85 as the expected value. (Again computed by a Claudicant simulation.)

While the conclusion that the Normal approximation is pretty accurate with so many terms in the Binomial variates is quite unsurprising, Claude saves me coding time without ruining the puzzle altogether. (And Gemini made me aware that the name of the café where Gilles and I regularly meet, L’Écir, is an Auvergne noun for a local, dangerous, mountain blizzard! Thus linking the place to the foundation of the café by Auvergne expatriates…)

Last month, Vicki Boykis recommended the book “Code that Fits in Your Head” by Mark Seemann. I was intrigued, especially by her takeaway that “Writing good software should be a slow and deliberate craft”.

The book reminded of The Pragmatic Programmer in the breadth of topics it covers. However, I much preferred its tone.

Sustainability

Throughout the book, the author makes it clear that he believes one should pay attention to the code one writes, to its internal quality. The following two sentences even got their little frame:

“The goal is not to write code fast. The goal is sustainable software.”

Arrange Act Assert

What useful phrase!

You know how in a test you can have something like:

test_that("Exception Bla handled", {
  withr::local_option("bla" = 1)
  test_thing <- 23

  x <- my_function(test_thing)

  expect_equal(x, 2)
})

The first two lines within the test_that() call prepare what’s needed. That’s the arrange phrase.

The line calling my_function() does the thing I want to test. That’s the act phrase.

The line calling expect_equal() checks that thing. That’s the assert phrase.

Arrange, act, assert! The book even mentions you can separate the three phases with empty lines.

Red Green Refactor

Another neat phrase! It describes how you could iteratively work on a part of your software.

• Red: you write a failing test for e.g. a feature.
• Green: you write the code to make the test pass.
• Refactor: you improve that code (but don’t make the test fail again!).

I appreciated how the author wrote that even the red phase isn’t easy, because you can get a passing test that you thought would fail. I think the point is made again elsewhere in the book: you can imagine a test will pass or fail under certain conditions, but better to actually prove it by running the test!

The Devil’s Advocate

A technique presented in the book is The Devil’s Advocate, in which you write wrong code on purpose, to see whether your tests will detect it. I suppose it’s a less random (but more difficult?) version of mutation testing, in which you evaluate how often your test suite detects mutants of your code. In R, mutation testing is provided by e.g. the mutator package.

Code reviews that do not block

There’s a short discussion of when to do code reviews (regularly) so as to not block your team mates. I enjoyed seeing this, as it’s important to not always be a rate-limiting factor. It shows the book’s advice is pragmatic.

Bundle smaller breaking changes

In the chapter “Augmenting code” that’s about working with existing code, the author – among other things – discusses whether to bundle or separate breaking changes into one or several releases. The decision criterion is how much work you create for “client developers” (maintainers of reverse dependencies).

It’s also one of the many places in the book where the author reminds us there’s no hard rule, that it’s the “art of software engineering”.

House or not House

The first chapter of the book presents and criticizes analogies of software engineering, like comparing developing software to building a house. That reminded of the talk Maintaining the house the tidyverse built by Hadley Wickham, that compared tidyverse maintenance to house maintenance.

The chapter’s thesis is that all analogies are useful and imperfect and you shouldn’t let them limit or cloud the way you view your work.

Justify exceptions

At some point the author writes that static code analysis brings false positives, but that if you disable specific rules, you should justify why. That’s how Jarl, a linter for R code, makes you suppress rules: you have to explain a reason.

Git is easy

I like these encouraging three sentences:

“Git isn’t the most user-friendly piece of technology on the planet, but you’re a programmer. You’ve managed to learn at least one programming language. Compared to that, learning the basics of Git is easy.”

Testing against databases

As an example of slow tests, the author mentions tests against databases. I want to use this as an excuse to plug the dittodb package to mock databases in tests.

Bisection without and with Git

In the chapter about Troubleshooting, the idea of bisection without Git is introduced. The author says it’s called bisection for lack of a better word. The idea being to go from a big piece of code with a bug to the smallest piece of code with the same bug, I suppose a better word is reprex.

“Being able to produce a minimal working example is a superpower in software engineering.”

Regarding Git bisection, here’s your friendly reminder that you can try it out using the saperlipopette R package.

Time boxing

The author extols the virtues of time boxing (e.g. the Pomodoro method). I’ve been using Entracte and would recommend giving it a go, if you have no similar setup yet!

When reading the book, between chunks of reading I’d crochet a few rounds.

Performance

I enjoyed the pragmatic discussion of performance: don’t forget to check the numbers and their significance.

Behavioural code analysis

I have never used behavioural code analysis, something that uses Git data, but its presentation reminded me of Git commands to get to know a project which is related (or the same?).

Code navigation

Yay to the mention of learning how to navigate code in your IDE. Such an important skill.

Conclusion

“Code that Fits In Your Head” is an ambitious (and rather long) book. I was glad to learn new phrases and to get the opportunity to get to know or reflect on so many topics.

Because the book was published in 2021, AI only makes an appearance as a brief footnote. The only moment I really thought AI would make a point less valid was a point about some kinds of refactoring taking ages: in the igraph R package, some changes have recently become feasible in an easier way thanks to using tools like Claude.

Some aspects of the book might be a tad annoying like the super simple diagrams or illustrations (think: a hammer drawing to illustrate the fact that to someone with a hammer, everything looks like a nail) but I suppose those play their role of breaking up pages.

Bonus: a reading list

If someone fairly new to software engineering were to ask me book recs, from looking at the stack near my desk I would recommend the following books, more specialized than “Code that Fits In Your Head”:

• The Art of Readable Code by Dustin Boswell and Trevor Foucher
• The Programmer’s Brain by Felienne Hermans (that I plan to re-read soon!)
• A Philosophy of Software Design by John Ousterhout
• A Git book, either Git in Practice by Mike McQuaid or Pro Git by Scott Chacon
• A book about team work. Maybe The Psychology of Software Teams by Cat Hicks, The Collective Edge by Colin M. Fisher, Dare to Lead by Brené Brown…
• A book about productivity/organization, but not a sanctimonious one about profound labor or the like. I enjoyed Make Time by Jake Knapp and John Zeratsky and Tranquility by Tuesday by Laura Vanderkam.
• A book critical of AI such as The AI Con by Emily M. Bender and Alex Hanna or Empire of AI by Karen Hao.

I also believe there are other ways to learn, like watching talks, but this post is about books.

What Is Posit Assistant?

Posit Assistant is an AI assistant built into your IDE, either Positron or RStudio. It’s built primarily for exploratory data analysis, and its pitch is session context and specialised, in-built skills.

If you follow Posit’s developments, you might have heard of their previous iterations of an integrated AI assistant. The naming has been a moving target:

• Positron Assistant came first, in early 2025, as a general coding assistant inside Positron
• Databot followed that August, aimed squarely at data exploration
• Posit Assistant replaces both, and is the first to be available in RStudio as well as Positron

In this blog, we will

• Take a look at a typical chat with Posit Assistant
• Explore what sets it apart from other integrated AI assistants
• Introduce the btw R package, for those who’d rather add R session context to their existing assistant instead

A Typical Chat

A great way to see it in action is to load some data and ask a broad question.

glimpse(happy)

happy |>
 summarise(across(everything(), ~sum(is.na(.)))) |>
 select(where(~. > 0))

drivers = select(happy,
 `GDP per capita`:`Perceptions of corruption`)

cor(happy$Score, drivers) |>
 round(2) |>
 t()
#> [,1]
#> GDP per capita 0.79
#> Social support 0.78
#> Healthy life expectancy 0.78
#> Freedom 0.57
#> Generosity 0.08
#> Corruption 0.39

Firstly, nothing in my request told it what the data represents, what columns existed or what the values looked like. Yet it came back with a description of the data, and a simple correlation analysis showing which variables drive happiness, and three suggestions of what you may want to explore next.

The other thing I like is that it answers by writing and then running code, which makes its decisions traceable. Language models are not designed to do maths (clue is in the “language” bit), but by writing code that does the maths, which we can scrutinise, copy and tweak, it gives us a much more trustworthy output.

How Does It Compare to Other Integrated Assistants?

Claude Code is probably the most popular alternative AI coding assistant, which is available inside Positron as an extension, and available anywhere you have a terminal through its CLI tool.

The two overlap on most of the basics. Both Posit Assistant and Claude Code:

• Can run on the same underlying model (though Posit Assistant also supports other providers)
• Can read and edit your files
• Generate and run code, from which they come to conclusions
• Can be given persistent project context files, through CLAUDE.md for Claude Code and AGENTS.md for Posit Assistant

So what sets them apart?

As we know now, Posit Assistant adds the session context. It knows the variables you have defined, the data frames you have loaded and the plot you are looking at, which is what made the exchange above possible.

Plots are one place where that makes a big difference. Posit Assistant can see the plots you previously made in the plot viewer, interpret them, and print plots back into the chat. That is worth a great deal in exploratory work, where much of what you learn comes from looking at a chart rather than a table. Claude Code can create a plot and save it to an image file, or ingest an image file, but the process is much clunkier.

It is not infallible, and Posit are candid about that: their research on how well LLMs interpret plots found models read straightforward charts well, but stumble when what is on screen contradicts what they expect to see.

It also brings the skills, which are bundles of instructions and code that the assistant pulls in by itself when it judges them relevant to what you asked. Thirteen ship built in, covering things like report authoring with Quarto and Jupyter Notebooks, and predictive modelling in R and Python, so it arrives at those tasks already knowing the conventions. You do not have to wait for it to spot the right one, though. Most skills are also exposed as slash commands, so typing /report, /quarto-authoring or /shiny-bslib runs that skill’s instructions there and then.

Claude Code is more autonomous, and will work through several steps before it checks in, whereas Posit Assistant stops for guidance more often. That is deliberate: it inherits Databot’s exploratory behaviour, so when it decides you are exploring data it runs a couple of tool calls, reports back and offers you some directions. The point is to keep you involved while the exploration is happening, rather than handing you a finished answer.

One thing to be aware of: Posit’s documentation lists “console history” as available context, which is a little misleading. Line by line detail of what has been run is not part of the context in either R or Python, so if you were hoping the assistant would know what code you’ve already run, or an error message that’s popped up, it will not.

Not ready to make the switch? Add context to your assistant with an MCP instead

If you would rather keep one assistant, it’s possible to give your AI assistant of choice some of the same session awareness with the btw R package, also developed by Posit. It works over the Model Context Protocol, and if MCP is new to you, Neal Richardson of Posit gave a good primer, “MCP, or Not MCP”, at our AI in Production conference.

The setup is straightforward, and the MCP setup guide walks you through it.

Once it is running, your assistant can inspect R objects and data frames in your session, check installed packages, and search for and describe CRAN packages.

There are trade-offs. It covers R only, so a Python session gets nothing from it. You also lose the skills, so your assistant comes to a Quarto document, Shiny app, or a modelling task without the guidance Posit Assistant would have loaded for you.

Closing Thoughts

Posit Assistant is easy to adopt, particularly if you already run the Claude Code extension in Positron or VS Code, because the workflow is familiar and the session context arrives without any configuration.

For exploratory data analysis it is an excellent tool. Being able to ask about the data you actually have loaded, and to read the code it wrote to answer you, changes the pace of that work more than I expected.

If you would rather keep to a single assistant, pair your go-to coding assistant with Posit’s btw package and you will get most of the way there.

This post grew out of “Improving Your Workflow with Positron and Claude”, a workshop we ran at AI in Production 2026. If sessions like that appeal, come and join us next year.

Find Out More

From us:

• A First Look at Positron and Posit Assistant – our free webinar on August 13th, 2026
• Positron vs RStudio – is it time to switch?
• “MCP, or Not MCP” – Neal Richardson (Posit) at AI in Production 2026

From Posit:

• A brief and biased history of Posit data science agents – Joe Cheng on why there were three of them
• Introducing AI in RStudio

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

It is important to know how to reshape a dataframe into a form that might be more suitable for the analyses we intend to perform. For me, it is also important to know how to do this with base R, without using any other packages. Yes, I know that other packages, such as dplyr and tidyr, are much better… Still, having the skills to perform most of our everyday tasks with base R may be a good idea. In particular, there are at least four routine tasks that we need to be able to perform:

1. subsetting
2. sorting
3. casting
4. melting

library(statforbiology)
students <- getAgroData("students")
head(students)
##   Id  Subject       Date Mark Year  HighSchool
## 1  1 AGRONOMY 10/06/2002   30 2001  HUMANITIES
## 2  2 AGRONOMY 08/07/2002   24 2001 AGRICULTURE
## 3  3 AGRONOMY 24/06/2002   30 2001 AGRICULTURE
## 4  4 AGRONOMY 24/06/2002   26 2001  HUMANITIES
## 5  5 AGRONOMY 23/01/2003   30 2001  HUMANITIES
## 6  6 AGRONOMY 09/09/2002   28 2001 AGRICULTURE

subData <- subset(students, Mark >= 28)
head(subData)
##    Id  Subject       Date Mark Year  HighSchool
## 1   1 AGRONOMY 10/06/2002   30 2001  HUMANITIES
## 3   3 AGRONOMY 24/06/2002   30 2001 AGRICULTURE
## 5   5 AGRONOMY 23/01/2003   30 2001  HUMANITIES
## 6   6 AGRONOMY 09/09/2002   28 2001 AGRICULTURE
## 11 11 AGRONOMY 09/09/2002   28 2001  SCIENTIFIC
## 17 17 AGRONOMY 10/06/2002   30 2001  HUMANITIES

subData <- subset(students, Mark <= 28 & Mark >= 26)
head(subData)
##    Id  Subject       Date Mark Year  HighSchool
## 4   4 AGRONOMY 24/06/2002   26 2001  HUMANITIES
## 6   6 AGRONOMY 09/09/2002   28 2001 AGRICULTURE
## 7   7 AGRONOMY 24/02/2003   26 2001  HUMANITIES
## 8   8 AGRONOMY 09/09/2002   26 2001  SCIENTIFIC
## 10 10 AGRONOMY 08/07/2002   27 2001  HUMANITIES
## 11 11 AGRONOMY 09/09/2002   28 2001  SCIENTIFIC

subData <- subset(students, Mark <= 28 & Mark >= 26 & 
                    Subject == "MATHS")
head(subData)
##      Id Subject       Date Mark Year  HighSchool
## 115 115   MATHS 15/07/2002   26 2001 AGRICULTURE
## 124 124   MATHS 16/09/2002   26 2001  SCIENTIFIC
## 138 138   MATHS 04/02/2002   27 2001  HUMANITIES
## 144 144   MATHS 10/02/2003   27 2001  HUMANITIES
## 145 145   MATHS 04/07/2003   27 2002  HUMANITIES
## 146 146   MATHS 28/02/2002   28 2001 AGRICULTURE

subData <- subset(students, Mark <= 28 & Mark >= 26 & 
                    (Subject == "MATHS" | 
                     Subject == "CHEMISTRY"))
head(subData)
##    Id   Subject       Date Mark Year   HighSchool
## 68 68 CHEMISTRY 04/06/2002   28 2001 OTHER SCHOOL
## 70 70 CHEMISTRY 04/06/2002   26 2001   ACCOUNTING
## 71 71 CHEMISTRY 04/06/2002   27 2001  AGRICULTURE
## 72 72 CHEMISTRY 23/01/2003   27 2001   SCIENTIFIC
## 75 75 CHEMISTRY 10/07/2002   27 2001  AGRICULTURE
## 81 81 CHEMISTRY 23/01/2003   28 2001  AGRICULTURE

subData <- subset(students, Mark <= 28 & Mark >= 26 & 
                    (Subject == "MATHS" | 
                     Subject == "CHEMISTRY"),
                  select = c(Subject, Mark, HighSchool))
head(subData)
##      Subject Mark   HighSchool
## 68 CHEMISTRY   28 OTHER SCHOOL
## 70 CHEMISTRY   26   ACCOUNTING
## 71 CHEMISTRY   27  AGRICULTURE
## 72 CHEMISTRY   27   SCIENTIFIC
## 75 CHEMISTRY   27  AGRICULTURE
## 81 CHEMISTRY   28  AGRICULTURE

subData <- subset(students, Mark <= 28 & Mark >= 26 & 
                    (Subject == "MATHS" | 
                     Subject == "CHEMISTRY"),
                  select = c(-Id, 
                             -Date,
                             -Year))
head(subData)
##      Subject Mark   HighSchool
## 68 CHEMISTRY   28 OTHER SCHOOL
## 70 CHEMISTRY   26   ACCOUNTING
## 71 CHEMISTRY   27  AGRICULTURE
## 72 CHEMISTRY   27   SCIENTIFIC
## 75 CHEMISTRY   27  AGRICULTURE
## 81 CHEMISTRY   28  AGRICULTURE

subData <- students[(students$Mark <= 29 & students$Mark >= 25),]
head(subData)
##    Id  Subject       Date Mark Year  HighSchool
## 4   4 AGRONOMY 24/06/2002   26 2001  HUMANITIES
## 6   6 AGRONOMY 09/09/2002   28 2001 AGRICULTURE
## 7   7 AGRONOMY 24/02/2003   26 2001  HUMANITIES
## 8   8 AGRONOMY 09/09/2002   26 2001  SCIENTIFIC
## 10 10 AGRONOMY 08/07/2002   27 2001  HUMANITIES
## 11 11 AGRONOMY 09/09/2002   28 2001  SCIENTIFIC

subData <- students
subData[(subData$Mark <= 29 & subData$Mark >= 25), "Mark"] <- NA
head(subData)
##   Id  Subject       Date Mark Year  HighSchool
## 1  1 AGRONOMY 10/06/2002   30 2001  HUMANITIES
## 2  2 AGRONOMY 08/07/2002   24 2001 AGRICULTURE
## 3  3 AGRONOMY 24/06/2002   30 2001 AGRICULTURE
## 4  4 AGRONOMY 24/06/2002   NA 2001  HUMANITIES
## 5  5 AGRONOMY 23/01/2003   30 2001  HUMANITIES
## 6  6 AGRONOMY 09/09/2002   NA 2001 AGRICULTURE

subData[is.na(subData$Mark), "Mark"] <- 0 
head(subData)
##   Id  Subject       Date Mark Year  HighSchool
## 1  1 AGRONOMY 10/06/2002   30 2001  HUMANITIES
## 2  2 AGRONOMY 08/07/2002   24 2001 AGRICULTURE
## 3  3 AGRONOMY 24/06/2002   30 2001 AGRICULTURE
## 4  4 AGRONOMY 24/06/2002    0 2001  HUMANITIES
## 5  5 AGRONOMY 23/01/2003   30 2001  HUMANITIES
## 6  6 AGRONOMY 09/09/2002    0 2001 AGRICULTURE

sortedData <- students[order(students$Mark), ]
head(sortedData)
##    Id   Subject       Date Mark Year   HighSchool
## 51 51   BIOLOGY 01/03/2002   18 2001   HUMANITIES
## 67 67 CHEMISTRY 20/02/2003   18 2002  AGRICULTURE
## 76 76 CHEMISTRY 24/02/2003   18 2002 OTHER SCHOOL
## 79 79 CHEMISTRY 18/06/2003   18 2002  AGRICULTURE
## 82 82 CHEMISTRY 18/07/2002   18 2001  AGRICULTURE
## 83 83 CHEMISTRY 23/01/2003   18 2001   SCIENTIFIC

sortedData <- students[order(-students$Mark), ]
head(sortedData)
##    Id  Subject       Date Mark Year  HighSchool
## 1   1 AGRONOMY 10/06/2002   30 2001  HUMANITIES
## 3   3 AGRONOMY 24/06/2002   30 2001 AGRICULTURE
## 5   5 AGRONOMY 23/01/2003   30 2001  HUMANITIES
## 17 17 AGRONOMY 10/06/2002   30 2001  HUMANITIES
## 18 18 AGRONOMY 10/06/2002   30 2001 AGRICULTURE
## 19 19 AGRONOMY 09/09/2002   30 2001 AGRICULTURE

sortedData <- students[order(-students$Mark, students$Subject), ]
head(sortedData)
##    Id  Subject       Date Mark Year  HighSchool
## 1   1 AGRONOMY 10/06/2002   30 2001  HUMANITIES
## 3   3 AGRONOMY 24/06/2002   30 2001 AGRICULTURE
## 5   5 AGRONOMY 23/01/2003   30 2001  HUMANITIES
## 17 17 AGRONOMY 10/06/2002   30 2001  HUMANITIES
## 18 18 AGRONOMY 10/06/2002   30 2001 AGRICULTURE
## 19 19 AGRONOMY 09/09/2002   30 2001 AGRICULTURE

sortedData <- students[order(-students$Mark, -xtfrm(students$Subject)), ]
head(sortedData)
##      Id Subject       Date Mark Year   HighSchool
## 116 116   MATHS 01/07/2002   30 2001 OTHER SCHOOL
## 117 117   MATHS 18/06/2002   30 2001   ACCOUNTING
## 118 118   MATHS 09/07/2002   30 2001  AGRICULTURE
## 121 121   MATHS 18/06/2002   30 2001   ACCOUNTING
## 123 123   MATHS 09/07/2002   30 2001   HUMANITIES
## 130 130   MATHS 07/02/2002   30 2001   SCIENTIFIC

rimsulfuron <- getAgroData("rimsulfuron")
head(rimsulfuron)
##                      Herbicide Plot Code Block Column WeedCover Yield
## 1             Rimsulfuron (40)    1    1     1      1      27.8 85.91
## 2             Rimsulfuron (45)    2    2     1      1      27.8 93.03
## 3             Rimsulfuron (50)    3    3     1      1      23.0 86.93
## 4             Rimsulfuron (60)    4    4     1      1      42.8 52.99
## 5    Rimsulfuron (50+30 split)    5    5     1      1      15.1 71.36
## 6 Rimsulfuron + thyfensulfuron    6    6     1      1      22.9 75.28

castData <- reshape(
  rimsulfuron[, c("Herbicide", "Block", "Yield")], 
  direction = "wide",
  idvar = "Herbicide",
  timevar = "Block"
)
castData
##                                     Herbicide Yield.1 Yield.2 Yield.3 Yield.4
## 1                            Rimsulfuron (40)   85.91   91.09  111.42   93.15
## 2                            Rimsulfuron (45)   93.03  105.00   89.19   79.04
## 3                            Rimsulfuron (50)   86.93  105.82  110.02   89.10
## 4                            Rimsulfuron (60)   52.99  102.86  100.62   97.04
## 5                   Rimsulfuron (50+30 split)   71.36   77.57  115.91   92.16
## 6                Rimsulfuron + thyfensulfuron   75.28   82.59   94.96   85.85
## 7                        Rimsulfuron + hoeing   73.22   86.06  118.01   98.32
## 8    Pendimethalin (pre) + rimsulfuron (post)   65.51   88.72   95.52   82.39
## 9  Pendimethalin (post) + rimsuulfuron (post)   94.82   87.72  102.05  101.94
## 10                        Rimsulfuron + Atred   94.11   89.86  104.34   99.63
## 11                             Thifensulfuron   78.47   42.32   62.52   24.34
## 12         Metolachlor + terbuthylazine (pre)   51.77   52.10   49.46   34.67
## 13                  Alachlor + terbuthylazine   12.06   49.58   41.34   16.37
## 14                                Hand-Weeded   77.58   92.08   86.59   99.63
## 15                                 Unweeded 1   10.88   31.77   23.92   20.85
## 16                                 Unweeded 2   27.58   51.55   25.13   38.61

WeedPop <- getAgroData("WeedPop")
WeedPop
##   Code POLLA CHEPO ECHCG AMARE XANST POLAV
## 1    A   0.1    33    11     0   0.1   0.1
## 2    B   0.1     3     3     0   0.1   0.0
## 3    C   7.0    19    19     4   7.0   1.0
## 4    D  18.0     3    28    19  12.0   6.0
## 5    E   5.0     7    28     3  10.0   1.0
## 6    F  11.0     9    33     7  10.0   6.0
## 7    G   8.0    13    33     6  15.0  15.0
## 8    H  18.0     5    33     4  19.0  12.0
## 9    I   6.0     6    38     3  10.0   6.0
mdati <- reshape(
  WeedPop,
  direction = "long",
  varying = colnames(WeedPop[2:7]),
  v.names = "WeedCover",
  times = colnames(WeedPop[2:7]),
  timevar = "Weed name",
  ids = WeedPop$Code,
  idvar = "Code",
)
head(mdati)
##         Code Weed name WeedCover
## A.POLLA    A     POLLA       0.1
## B.POLLA    B     POLLA       0.1
## C.POLLA    C     POLLA       7.0
## D.POLLA    D     POLLA      18.0
## E.POLLA    E     POLLA       5.0
## F.POLLA    F     POLLA      11.0

I gave the first of this year’s Ihaka Lectures at the University of Auckland this past July. It revisited and partially updated things I’ve been talking about for the past eighteen months or so on trust in data visualization (and by extension in scientific work) in the light of the many challenges we presently face.

AI in Production is about the reality of running AI and machine learning systems, not just building them. The talks come from people who have shipped something and lived with what came next: what worked, what surprised them, what they would do differently. The second edition returns to The Catalyst in Newcastle upon Tyne on 10 and 11 June 2027, and the call for speakers is open until 29 January.

What to expect

Hands-on workshops on Thursday 10 June, then a full day of talks and lightning talks on Friday 11 June. On the workshop day you take one morning and one afternoon session, with lunch included and a dinner and drinks reception in the Catalyst atrium afterwards. Topics across both days cover deployment, monitoring and observability, generative AI, retrieval-augmented generation, and responsible AI as it applies day to day. The talks we like best are the ones that cover what went wrong as well as what worked.

The shape of both days is on the conference site, and the FAQ answers most of the logistics questions.

Call for speakers

Proposals close on 29 January 2027. We want talks grounded in a system that exists and gets used by someone other than its author: what worked, what surprised you, what you would do differently. Deployment, observability, incident response, evaluation and guardrails, cost and latency, governance in daily practice, and the team lessons that come with keeping any of it running.

If you are unsure whether your talk fits, try this. Swap “AI” for any other kind of software system. If the talk still holds up, it is probably not for us. If it falls apart because the AI-specific parts are central to the story, that is what we are looking for.

First-time speakers are welcome. Our guide to submitting conference abstracts walks through the writing, and why submit to AI in Production addresses the usual worry about whether your work is ready. Send a title and up to 300 words through the submission form.

Key dates

• 8 January 2027: Super early bird deadline
• 29 January 2027: Call for speakers closes
• 5 March 2027: Early bird deadline
• 4 June 2027: General registration deadline
• 10 June 2027: Conference begins

Tickets

Book the conference day on its own, or a combined ticket that includes the workshop day.

Last year at AI in Production 2026

The 2027 line-up is still coming together, and there is room on the programme if you would like to speak. In the meantime, last year’s talks are a fair guide to the range: law, water, ecology, online gambling, central government, higher education, clinical research and enterprise software, as well as the software consultancies and tooling vendors you would expect.

Nathan Bilton of the law firm Weightmans covered who carries the liability when a chatbot gives a customer the wrong answer. Mac Misiura of Red Hat went through open source guardrails for securing large language model applications at scale. Grant Beasley of the online gambling operator tombola, with our own Myles Mitchell, showed how deep learning can identify players at risk of gambling harm. Seb Ringrose of Doubleword explained why low latency, low cost and high quality are a pick-two problem, and George Stagg of Posit closed the day on what building Posit Assistant taught him about agents. Our summary of the 2026 conference covers every talk, and you can watch the full 2026 playlist on YouTube.

Getting here

The Catalyst is at 3 Science Square, Newcastle Helix, Newcastle upon Tyne NE4 5TG, a 10 minute walk from Newcastle Central Station. Trains take about three hours from London and 90 minutes from Edinburgh. Newcastle International Airport is about 30 minutes from Central Station on the Metro. St James’ Park is a short walk away.

Make a trip of it

There is more to the city than the venue, and NewcastleGateshead has the full list. A few favourites:

• The Quayside, where the tilting Millennium Bridge and the Tyne Bridge sit side by side, with a Sunday market under the arches
• Grey Street and Grainger Town, Georgian streets that regularly top “best street in the UK” lists, plus the covered Grainger Market
• Newcastle Castle, the medieval keep that gave the city its name, with rooftop views over the Tyne
• BALTIC and The Glasshouse, contemporary art and music across the river in Gateshead
• Ouseburn Valley, the creative quarter: converted warehouses, street art, a city farm and some of the best independent pubs in the North East
• St James’ Park, home of Newcastle United, with stadium tours most days

Sponsorship

If your organisation would like to support the conference, email events@jumpingrivers.com.

Super early bird tickets are on sale until 8 January. Book your place.

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

Join our workshop on From Raw Data to Regulatory Results: Clinical Trial Programming in R using the pharmaverse,  which is a part of our workshops for Ukraine series! 

Here’s some more info: 

Title: From Raw Data to Regulatory Results: Clinical Trial Programming in R using the pharmaverse

Date: Thursday, September 17th, 18:00 – 20:00 CEST (Rome, Berlin, Paris timezone) 

Speaker: Edoardo Mancini, Edoardo is a Senior Data Scientist working at Roche in the UK. His work involves leading statistical programming activities in early and late-stage clinical trials across multiple therapeutic areas, including Ophthalmology and Hematology. Edoardo is also a passionate advocate for R and Open Source, pursuing these interests in his current role as maintainer for the open-source {admiral} R package, which is part of the wider pharmaverse and helps create ADaM datasets. 

Description: Have you ever wondered how raw patient data from a clinical trial transforms into the actual evidence used to approve new medicines? This two-hour, hands-on workshop aims to answer this question by giving you a glimpse into the clinical trial data pipeline. We’ll start with a crash course on industry data standards, tracing the journey from raw data to SDTM (Study Data Tabulation Model), into analysis-ready ADaM datasets, and finally into Tables, Listings, and Graphs (TLGs). From there, you will step into the shoes of a clinical statistical programmer: firstly, you will be guided through creating your own simple ADaM dataset(s) using the {admiral} package and test data. Then, you will bring your analyses to life using {ggplot2} and a tabulation package of your choice to build the kinds of polished visualizations and summary tables that drive real-world medical breakthroughs. No prior clinical trial experience is required, just bring your foundational R skills and your curiosity!

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!

Fit a cross-validated elastic net with cv.glmnet(x, y) on a design matrix that comfortably fits in memory, and the call can still fail. The reason is that the function needs several times the size of x while it runs, and nothing anywhere in your script says so. This post is about giving that requirement a home – a memory specification, written in R next to the code that knows it – first for a plain sequential call, then for the same call running in parallel.

TL;DR

Write down what the call needs, as one annotation on otherwise ordinary code;

fit <- cv.glmnet(x, y) |> resources(memory(4 * object.size(x)))

Better still, the memory specification need not stay at the call site. Once cv.glmnet() carries its own declaration, parallelization frameworks can make use of it, e.g.

fits <- lapply(xs, FUN = cv.glmnet, y = y) |> futurize()

Either way, that one declaration is meant to do two jobs;

1. checked before the work starts – parallel or not – it fails fast instead of minutes or hours later, and
2. handed to a parallel framework, it also decides how many tasks can run in parallel.

The goal of the Resource Project is to study, design, and implement these features in the Futureverse. Feedback and suggestions are welcome.

Problem: We have no way of declaring memory needs

We check arguments all the time. A stopifnot(is.matrix(x)) at the top of a function is second nature and in our muscle memory. It is useful because it helps functions fail fast when the wrong arguments are passed, and because we can give an informative error message. However, we have nothing for checking whether the machine is actually capable of running the function.

For example, you might not have enough memory available to perform a calculation. If you run out of memory, you might get:

> fit <- cv.glmnet(x, y)
Error: cannot allocate vector of size 1.7 Gb

Worse is when the operating system steps in first and kills the R process outright. There is no R error to catch, no traceback, and no chance to shut down gracefully:

$ Rscript fit-model.R
Killed
$

Neither one tells you how much was needed, and both arrive after the work is already underway.

After you have identified the problem to be a memory issue, it often becomes a trial-and-error game, and it can take a while to figure out how much memory you need. The finding might then end up as a source-code comment in the script, or as a rarely-read sentence in a package help page. Either way, since the memory needs are not in the code, R cannot check it, cannot size the parallelization from it, and cannot give informative error messages when you run out of memory.

Proposal: Check memory and fail fast

From rudimentary measurements and code inspections of glmnet, I found that cv.glmnet(x, y) needs roughly 3-4 times the memory of object.size(x) in addition to x itself. We can improve on this memory model¹, but for simplicity, let’s assume cv.glmnet(x, y) needs 4 * object.size(x) additional memory to be successful.

Let’s take a numeric design matrix of 250,000 rows and 1,000 columns, which clocks in around 1.9 GiB of memory, as an example:

x <- matrix(rnorm(250e3 * 1e3), nrow = 250e3, ncol = 1e3)
format(object.size(x), units = "GiB")
#> [1] "1.9 GiB"

Imagine a freeMemory() function for querying how much memory the current R process has available. With that, we could gatekeep our fit <- cv.glmnet(x, y) call as:

stopifnot(freeMemory() >= 4 * object.size(x))
#> Error: freeMemory() >= 4 * object.size(x) is not TRUE

That fails instantly and before attempting the model fit, if there is not enough memory available. We could also imagine a richer vocabulary that provides us with more informative error messages, e.g.

assert_resources(memory(4 * object.size(x)))
#> Error: UnmetResourceError: requires memory 7.5 GiB, available 3.1 GiB

In this project, we’re proposing an expr |> resources(...) syntax for declaring and asserting resource needs, including memory requirements, next to the code where it applies. In our example, it would look like:

fit <- cv.glmnet(x, y) |> resources(memory(4 * object.size(x)))
#> Error: UnmetResourceError: cv.glmnet(x, y) requires memory 7.5 GiB, available 3.1 GiB

This syntax preserves the original code and logic as-is, while allowing you to declare resource requirements that R can act on. In its most basic form, it effectively works like:

fit <- {
  assert_resources(memory(4 * object.size(x)))
  cv.glmnet(x, y)
}

If you wonder whether this is worth the extra code, consider who pays when it is missing. Even an experienced user might spend hours locating the reason for a “cannot allocate vector” error or an “OOM-killed” message. A researcher new to high-performance compute (HPC) clusters might lose a day to it, file a support ticket, or quietly conclude that R cannot handle their data.

The goal of this project is to not only reduce the amount of wasted compute resources, but also wasted human resources.

Batch, still sequential

Now suppose we are not fitting one model but ten. Consider a study of patients with some clinical endpoint (y) profiled in several ways, including expression, copy number, miRNA, methylation, somatic mutation, chromatin accessibility, proteomics, phosphoproteomics, metabolomics, and lipidomics. Assume we want to know which of these omics datasets predict the outcome the best. We have one y, and ten x design matrices in a list:

xs <- list(...)  # one design matrix per assay, ~8 GiB of them in total

Without memory protection, we would fit these models as:

fits <- lapply(xs, function(assay) cv.glmnet(assay, y))

The resolution of the technologies varies greatly, so the different x matrices vary greatly in size. Methylation is often of the highest resolution. Let’s assume its design matrix is 1.9 GiB in position xs[[4]], and the smallest is a tenth of that. Similarly to before, we could protect against memory overuse by using:

fits <- lapply(xs, function(assay) {
  cv.glmnet(assay, y) |> resources(memory(4 * object.size(assay)))
})
#> Error: UnmetResourceError: cv.glmnet(assay, y) requires memory 7.5 GiB, available 3.1 GiB

This tells you that one of the model fits would fail due to insufficient memory. Unfortunately, it does not fail instantly - it only fails when it tries to fit the too-large design matrix. Given that the largest matrix in this case happens to be in position four, you have already wasted efforts processing three cross-validation fits before failing. In the worst case, it could have processed nine out of the ten design matrices, before failing.

It would be better if the map-reduce call fails instantly, before attempting any model fits at all. We could make this happen if we move the declaration outside, so that the whole check happens prior to the model fits:

fits <- lapply(xs, function(assay) cv.glmnet(assay, y)) |>
  resources(function(assay) memory(4 * object.size(assay)))
#> Error: UnmetResourceError: cv.glmnet(assay, y) requires memory 7.5 GiB for xs[[4]], available 3.1 GiB

What is new is that the specification is a function rather than a fixed declaration. Somewhat simplified, this effectively checks the memory needs for all cross-validation fits first, before fitting them. Something like:

fits <- {
  assert_resources(lapply(xs, function(assay) memory(4 * object.size(assay))))
  lapply(xs, function(assay) cv.glmnet(assay, y))
}

This is the argument from the top of this post applied to a batch rather than to a call, and it is where failing fast saves the most. You avoid wasting all calls, which means less wasted compute resources, faster troubleshooting, and quicker fixes.

The declaration belongs on the function

In the above two examples, the caller had to declare the memory needs. That is a bit backwards. In order to do this, I had to investigate what cv.glmnet() needs, and nobody should have to repeat that exercise for every modeling function they call. It would be better if the author of the function could specify that. One approach would be to attach a resource-specification function² to the function itself;

resources(cv.glmnet) <- function(x, ...) memory(4 * object.size(x))

Technically, that sets attr(cv.glmnet, "resources") after having validated the function definition.

Next, we can have resources() look for such “resources” attributes and use them as the default, if found. So, if set, the memory-asserting call would become:

fit <- cv.glmnet(x, y) |> resources()
#> Error: UnmetResourceError: cv.glmnet(x, y) requires memory 7.5 GiB, available 3.1 GiB

Note how the burden is no longer on the caller, but on the function maintainer, to declare resources.

The map-reduce call works the same way, once we tell it which function is used:

fits <- lapply(xs, function(assay) cv.glmnet(assay, y)) |> resources(cv.glmnet)
#> Error: UnmetResourceError: cv.glmnet(assay, y) requires memory 7.5 GiB for xs[[4]], available 3.1 GiB

Handed a function that carries a declaration, resources() will use those declarations by default³.

Functions guarding themselves

With a resource function attached to the function, the cv.glmnet() function could easily guard the resources upfront by using⁴:

cv.glmnet <- function(x, y, ...) {
  assert_resources()
  # ... the real work
}

This would remove the need for the caller to use |> resources();

fit <- cv.glmnet(x, y)
#> Error: UnmetResourceError: cv.glmnet(x, y) requires memory 7.5 GiB, available 3.1 GiB

Attaching it to code you do not maintain

Even if a function does not carry a resource declaration, you can attach one yourself, e.g.

cv.glmnet <- glmnet::cv.glmnet
resources(cv.glmnet) <- function(x, ...) memory(4 * object.size(x))

That might be worth doing even for a single function in a single analysis script, preferably at the top of the script. It helps to gather such declarations in one place and avoids cluttering up the code, especially if the same function is used in multiple places.

Parallel: the same memory specification decides how many tasks run at once

So far, we’ve only discussed resource specification in sequential processing. Another goal of this project is to make use of them also in parallel processing. For instance, already today we can use the futurize package and its futurize() function to parallelize lapply(), purrr::map(), foreach(), and friends through that single pipe. In our case, we could use:

fits <- lapply(xs, function(assay) cv.glmnet(assay, y)) |>
  futurize()

Without memory protection, there is a great risk that you will run out of memory before you run out of CPU. Because of this, the goal is to support also:

fits <- lapply(xs, function(assay) cv.glmnet(assay, y)) |>
  resources(function(assay) memory(4 * object.size(assay))) |>
  futurize()

Here the resource specifications would not only protect against memory overuse, but also be used for scheduling the parallel tasks, e.g. by limiting the number of memory-hungry parallel tasks running concurrently on the same machine.

On a single machine

Consider using plan(multisession) where you have access to 32 GiB of memory and 16 CPU cores. That initiates 16 parallel workers. Each worker is a separate R process, and each task exports the matrix it is about to fit, so a task costs that matrix plus what fitting cv.glmnet() needs on top of that.

It is unlikely that you will be able to run 16 parallel model-fitting tasks with only 32 GiB. Even before turning to parallelization, you know that your design matrices in xs occupy 8 GiB of that memory, leaving you with at most 24 GiB for the parallel tasks. To fit the methylation assay, you need to export the 1.9 GiB matrix to the worker and then another 7.5 GiB to fit it, totaling 9.4 GiB. A task on one of the small assays costs a little under 4 GiB.

If you launched 16 parallel fits blindly, it’s quite likely that you would run out of memory and the operating system’s out-of-memory (OOM) killer may terminate your analysis. Given that none of this resource-specification framework exists today, the best you can do is to work the numbers yourself, and spin up only as many parallel workers as you can afford:

## 32 GiB RAM, ~8 GiB for the R session holding 'xs', 9.4 GiB for the largest task
plan(multisession, workers = parallelly::availableCores(max = (32 - 8) / 9.4))
#> 2 workers

Here we’re using availableCores() from the parallelly package, which respects common CPU allocations, while also limiting it manually via a custom memory-limit equation. Two concurrent parallel tasks is the best guess you have. However, with the above resource specifications, futurize() could probably do better and fit additional cross-validation models concurrently, especially the smaller ones.

On a compute cluster

On a high-performance compute (HPC) cluster, the job scheduler (e.g. Slurm and SGE) decides which execution nodes your job lands on, based on the amount of memory it requests. For example, Slurm declaration #SBATCH --mem=10G tells the scheduler that this job requires 10 GiB of RAM to run.

By declaring such memory needs within R;

plan(future.batchtools::batchtools_slurm)

fits <- lapply(xs, function(assay) cv.glmnet(assay, y)) |>
  resources(function(assay) memory(4 * object.size(assay))) |>
  futurize()

the future framework could work together with future.batchtools to translate each of the calculated resource needs into declarations understood by the job scheduler, which then can find appropriately sized slots on the cluster - all while maximizing the memory use but without ever running out of memory.

Ideally, everything is hidden away

Just as with resources(), if cv.glmnet() declares its own resource needs, futurize() can also take advantage of that. That would close the circle such that code existing already today, e.g.

fits <- lapply(xs, FUN = cv.glmnet, y = y) |> futurize()

would become resource aware, protect against overuse, and optimize scheduling overnight - all without code changes.

Outro

Phew, that was quite long, and yet, I only got to cover a tiny bit of what the Resource Project aims for. I discussed how we can manage memory from within R, but there are many other compute resources that limit us. For example, we also want to manage walltime, scratch space, GPU cores, and GPU memory.

If you have other thoughts or ideas, we’d love to hear from you. Please reach out on the Futureverse Discussions forum.

May the future be with you!

Henrik

I created R-Community for R users that are just starting with R or that have been using it for years.

I miss when I used Twitter and people had the “#rstats” hashtag that lead to a nice community. I closed my Twitter account when Mr. Elon Musk bought the platform, and similarly many users left to Mastodon, Bluesky, and other platforms. I thought Google Groups would be a good option as it does not require to configure an extra account as other technologies do (e.g., Slack, Discord, etc.). I also thought about moderation, as I received some trolling/bullying when I was a novice R user on Stackoverflow a few years ago.

This also follows from AI rise, which I find interesting but it replaces a lot of socialization in a world where smartphones already replaced parts of it in daily life. I tested a few AI engines things I have already worked on, and I observe that such engines often over-complicate answers.

If you are based in London, I am interested in reactivating the London R User Group as I think that PhD life sometimes puts you in an echo chamber that tends to over-specialization. I am willing to organize events to talk about R and socialize over pints or pizza.

You can join the group using this form.

Below is the group description.

Goal: to have a space to ask and answer R questions. No R preference dominates this space. Questions can be about base R, Tidyverse, data.table or anything in the wide R ecosystem.

Code of conduct: No bullying, hate, or rude acts of any kind. What is hateful to you, do not do to your fellow.

Community over Code: A strong, diverse community matters more than the software itself, because a healthy group can always improve bad code.

Earned Authority (Meritocracy): Influence and decision-making power are earned through active, public contributions.

Community of Peers: All individual contributors have equal weight in discussions, and hierarchical corporate structures do not apply here.

Openness and Transparency: Technical discussions happen in public view on this mailing list so everyone can participate conditional on the code of conduct.

Neutrality: Discussion about problem-solving, need for a new package or others are independent of any single vendor or corporate sponsor, protecting the long-term public good.

I have often found myself needing to plot a single graph with two y-axes having different scales. For example, this might be useful for representing temperature and rainfall data at a given location. Unfortunately, doing this with ggplot() is not straightforward.

While searching for possible solutions, I discovered that graphs with dual y-axes are generally not regarded as good tools for data visualisation. Hadley Wickham, the author of ggplot2 and several other important R packages, has given some good reasons for this (e.g., at this link). I do not intend to question these general arguments. Nonetheless, the idea that I cannot do with ggplot() something that I could easily do with Excel gets on my nerves. After all, such graphs can be quite useful in some specific circumstances, such as when displaying weather data.

I therefore started looking for a reasonable solution and eventually found my way. I would like to share it in this post.

Let’s consider the dataset DailyMeteoData.csv, which contains daily average temperature and rainfall data for two years at a location in my region (Umbria, central Italy). Let’s open the dataset and use dplyr to make a few useful transformations, such as:

1. converting the date character string into a date object;
2. adding variables for Year, Month, and Day of the Year (DOY).

library(dplyr)
fileName <- "https://www.casaonofri.it/_datasets/DailyMeteoData.csv"
dataMeteo <- read.csv(fileName) |>
  mutate(Date = as.Date(Date, format = "%d/%m/%Y"),
         Year = as.numeric(format(Date, format="%Y")),
         Month = as.numeric(format(Date, format="%m")),
         DOY = as.numeric(format(Date, format="%j")))
head(dataMeteo)

        Date Tavg Rain Year Month DOY
1 2011-01-01  5.9  0.0 2011     1   1
2 2011-01-02  6.2  0.6 2011     1   2
3 2011-01-03  4.5  0.0 2011     1   3
4 2011-01-04 -0.3  0.0 2011     1   4
5 2011-01-05  2.3  0.2 2011     1   5
6 2011-01-06  6.9  0.0 2011     1   6

Temperature is a continuous variable and can be easily represented using a line graph, whereas rainfall consists of discrete events and is therefore more usefully accumulated over periods such as ten-days or a month. We can use dplyr once more to create a new rainfall dataset containing the accumulated monthly rainfall, which is more suitable for our purposes. For the sake of simplicity, we assume that the year is divided into 12 months of equal length (approximately 30.4 days), and we calculate the DOY corresponding to the central day of each month, which will be used as the centre of the respective plot bar.

dataMeteo2 <- dataMeteo |>
  group_by(Year, Month) |>
  summarise(Rain = sum(Rain)) |>
  mutate(DOY = seq(15, 365, by = 365/12))
head(dataMeteo2)

# A tibble: 6 × 4
# Groups:   Year [1]
   Year Month  Rain   DOY
  <dbl> <dbl> <dbl> <dbl>
1  2011     1  40.4  15  
2  2011     2  32.8  45.4
3  2011     3 113.   75.8
4  2011     4  16.6 106. 
5  2011     5  27.4 137. 
6  2011     6  61.2 167. 

Now we can produce our first graph showing both rainfall and temperature. In the box below, I have played a little with the x-axis ticks and labels, and I have left the y-axis label empty for the moment.

library(ggplot2)
ggplot() +
  geom_bar(dataMeteo2, mapping = aes(x = DOY, y = Rain), fill = "grey",
           stat = "identity", width = 28) +
  geom_line(dataMeteo, mapping = aes(x = DOY, y = Tavg), col = "blue") +
  scale_x_continuous(breaks = c(365/12, 365/12*4, 365/12*8, 365) - 15, 
                     labels = c("Jan", "Apr", "Aug", "Dec"),
                     name = "") +
  scale_y_continuous(name = "") +
  facet_wrap(~Year) +
  theme_bw()

The previous graph does not work very well: the temperature line is hardly visible because its measurement scale is much smaller than that of rainfall. We therefore need to ‘scale’ the temperature variable so that it ranges approximately from 150 to 250. This will make the blue line clearly visible without interfering with the rainfall bars. The minimum and maximum temperature values are -4.2°C and 28.9°C, respectively; we want to map these original values to the new values 150 and 250, respectively, as shown in the figure below.

The previous figure tells us that we could transform the original temperature scale by using the equation of the straight line passing through the two points (-4.2, 150) and (28.9, 250). Thanks to what we remember from geometry courses, we can calculate the slope of such a straight line as:

while the intercept is:

Thus, the scaling equation is:

where is the new temperature scale, while is the original one. The reverse transformation is:

Now, we are ready to plot the graph. We transform the temperature to the new scale and plot it; furthermore, we include the second axis by using the sec.axis argument and the sec_axis() function, in which we specify the back-transformation function to the original scale.

dataMeteo <- dataMeteo %>% 
  mutate(newTavg = 162.69 + 3.02 * Tavg)

ggplot() +
  geom_bar(dataMeteo2, mapping = aes(x = DOY, y = Rain), fill = "grey",
           stat = "identity", width = 28) +
  geom_line(dataMeteo, mapping = aes(x = DOY, y = newTavg), col = "blue") +
  scale_x_continuous(breaks = c(365/12, 365/12*4, 365/12*8, 365) - 15, 
                     labels = c("Jan", "Apr", "Aug", "Dec"),
                     name = "") +
  scale_y_continuous(name = "Rain (mm)", 
                     sec.axis = sec_axis(~ (. - 162.69)/3.02, 
                                         name = "Daily Temperature (°C)",
                                         breaks = c(-10, 0, 10, 20, 30))) +
  facet_wrap(~Year) +
  theme_bw()

And we are done!

library(lubridate)
library(ggh4x)
dataMeteo1 <- dataMeteo |>
  mutate(
    Date = as.Date(Date, format = "%d/%m/%Y"),
    Year = lubridate::year(Date),
    Month = lubridate::month(Date))

dataMeteo2 <- dataMeteo1 |>
  group_by(Year, Month) |>
  summarise(
    Rain = sum(Rain),
    Temp = mean(Tavg, rm.na = TRUE)) |>
  mutate(Mth = factor(month.abb[Month], 
                      levels = c("Jan", "Feb", "Mar", "Apr",
                                 "May", "Jun", "Jul", "Aug",
                                 "Sep", "Oct", "Nov", "Dec")))
dataMeteo3 <- dataMeteo2 |>
  mutate(newTemp = 162.69 + 3.02 * Temp)

ggplot(dataMeteo3, mapping = aes(x = interaction(Mth, Year),
                                 group = 1)) +
  geom_col(aes(y = Rain), fill = "grey") +
  geom_point(aes(y = newTemp), size = 2, col = "blue") +
  geom_smooth(
    method = "gam", formula = y ~ s(x, bs = "cc"),
    aes(x = as.numeric(interaction(Mth, Year)), y = newTemp),
    col = "blue") +
  geom_point(aes(y = newTemp), size = 3, col = "blue", 
             fill = "white", shape = 21, stroke = 1) +
  scale_y_continuous(
    minor_breaks = scales::breaks_width(20),
    name = "Mean Monthly Total Rainfall (mm)",
    sec.axis = sec_axis(~ (. - 162.69) / 3.02,
      name = "Mean Monthly Daily Temperature (°C)",
      breaks = c(-10, 0, 10, 20, 30))) +
  guides(x = "axis_nested",
        y = guide_axis(minor.ticks=TRUE),
        y.sec = guide_axis(minor.ticks=TRUE)) +
  theme_bw(base_size = 18) +
  theme(
     panel.grid.minor = element_blank(),
     axis.title.x = element_blank(),
     axis.text.x = element_text(face = "bold", size = rel(0.5), angle = 90),
     axis.ticks = element_line(colour = "red"),
     ggh4x.axis.nestline.x = element_line(linewidth = 0.6),
     ggh4x.axis.nesttext.x = element_text(colour = "blue", 
                                          face = "bold", 
                                          size = rel(1.0))
  )

TLDR: This is the third post in a series examining relationships between the rule of law and other institutions as measured by the World Bank’s Worldwide Governance Indicators (WGI). Earlier posts demonstrated that the World Justice Project’s (WJP) measure can be proxied by the WGI Rule of Law (RoL) index, to take advantage of its wider coverage over a longer time-period.

This post builds on that analysis by asking how institutions should be conceptualized and analyzed: can the rule of law (and other institutions) sensibly be examined one at a time, or do they need to be treated as parts of an interdependent system? Using the six dimensions of the WGI as a proxy for institutional strength, results suggest dimensions are indeed bound together. Countries that score well on the rule of law also frequently score well on government effectiveness and corruption controls, year-to-year movements are positively linked across most dimensions and a common measure (or cause) accounts for the bulk of the variation between countries.

None of this will surprise anyone familiar with the literature and the analysis isn’t intended to prove institutions matter or that the WGI measure them. The point is narrower: to show why the rule of law can’t be examined in isolation and the risks of coming to the wrong conclusions when we try. It also sets up later posts in the series which will look at how the rule of law is connected with other national characteristics, such as economic growth.

Background

Institutions describe the constraints that structure political, economic and social interaction, such as how power is gained and exercised, how order is kept, and the way public resources are sourced, divided and used. Institutions shape incentives and opportunities: good institutions lead to more of the things we care about, like economic growth, peace and prosperity; while bad institutions lead to more of the thing we want less of, like poverty, inequality and conflict.

A sizable share of my work sits in economic development, so institutions are never far from mind. This series, though, started with a specific job: a client hired me to look at the links between the rule of law and economic growth, and I went looking for an accessible introduction to the topic I could share. What I found was this analysis from the Atlantic Council on why the rule of law is the key to prosperity.

Having noted how important institutions are and how difficult they are to define and measure, the authors move briskly to their conclusion that the rule of law is the single most influential factor behind long-term economic growth and societal wellbeing:

Source: Annie (Yu-Lin) Lee, Joseph Lemoine, 20 August 2025, Why the rule of law is the key to prosperity: Lessons from thirty years of data, Atlantic Council (link) [Accessed 5 August 2026]

That’s quite the claim and one of the inspirations for this post as from the looks of it, much of their evidence seems to come from the strength of a collection of pairwise correlations between the rule of law and their chosen proxies for well being (as measured by the Atlantic Council’s Freedom and Prosperity Index).

Now, I’m no big city lawyer, but it seems to me that pairwise correlations are a poor piece of evidence to provide for such a strong claim. Particularly when we’re wading into the murky waters of prosperity, institutions and the causal connection between the two. And the characteristics being analyzed are plausibly all part of the same interdependent system.

But, I’m not here to judge. Firstly, as the authors do mention that prosperity depends on the interplay of multiple institutional pillars, not just the rule of law. They’ve also tackled a complicated topic in an accessible way, which is an achievement in itself. And I’ve been there before too: as cross-country correlations are fun to explore and can provide a seemingly endless array of plausible policy interventions for making the world a better place. Also, the point of this post isn’t to criticize their analysis, but to help fill the gap I noticed when searching for resources on the topic.

Instead, this post attempts to fill a gap by presenting analysis demonstrating why the rule of law, and institutions more generally, are best conceptualized as interdependent pillars within a mutually reinforcing system. That interdependence is what makes the literature so conceptually interesting and so frustrating to analyze, since the little data available on the topic arrives bundled with multicollinearity, endogeneity, collider bias, omitted variables and measurement error.

The rule of law and institutions: a primer

To obnoxiously paraphrase research¹, the rule of law is thought to influence economic development and wider prosperity through a variety of avenues, such as enabling the enforcement of property rights, easing trade between unrelated parties and providing a peaceful means for resolving disputes. However, because the rule of law and its outcomes also depend on a wider set of institutions and a tangled coalition different interests, analyzing it is no simple task.²

On one level this is because institutions and their outcomes are interdependent. For instance, maybe the rule of law does drive prosperity, but perhaps prosperous countries are better able to invest in their legal systems too. Similarly, perhaps how the rule of law drives economic outcomes changes depending on the political system, geography and urbanization.

Defining and measuring amorphous concepts like the rule of law also happens to be hard, which means researchers lean heavily on surveys asking respondents for their perception of corruption, law and order and government effectiveness (etc). However, because a person’s opinion of a country’s performance in one area is likely to be influenced by their impression of it in others, these measures are likely to agree with one another for reasons that have little to do with the institutions themselves. A respondent who has watched a corruption scandal unfold is unlikely to rate the courts generously that same year.

Finally, it’s generally accepted that institutions are slow moving (or ‘sticky’), which means they don’t change much from one year to the next and will often exert their influence on wider outcomes indirectly. As a result, some of the most influential studies exploring the connection between institutions and prosperity analyse periods of a hundred years or more, such as Acemoglu, Johnson and Robinson who compared current income with a century old proxy for institutional quality (settler mortality).

Figure: Income vs. Settler mortality

Countries with higher GDP per capita now tended to be those with stronger colonial institutions (as proxied by settler mortality)

Source: Acemoglu, Daron, Simon Johnson, and James A. Robinson. 2001. “The Colonial Origins of Comparative Development: An Empirical Investigation.” American Economic Review 91 (5): 1369–1401. DOI: 10.1257/aer.91.5.1369

The Worldwide Governance Indicators

Although there are good reasons to question whether the WGI provides a reliable and holistic measure of institutions, (or even that it measure what it claims to) it is arguably one of the more reliable and well-tested attempts at measuring the cross-country quality of institutions. I also have a personal preference for the WGI as somebody who is hired to design, build and evaluate composite indicators, because it gets a core part of index design right: transparently sharing their methodology and data, and being open to making revisions to reflect feedback.

The WGI also has the practical advantage of being intuitive enough for outsiders to understand what it’s trying to measure. With governance reflecting the traditions and institutions by which authority in a country is exercised across six dimensions:

• Voice and Accountability: perceptions of the extent to which citizens can participate in selecting their government including electoral integrity, and of accountability mechanisms for citizens—reflected in the ability to access information, governmental oversight bodies, and a robust traditional/digital media landscape; and
• Political Stability: perceptions of the extent to which political power and governance are secure from destabilization, and of the likelihood that authority will be challenged or altered through violent, coercive, or unconstitutional means. The second aspect—government capacity to formulate and implement sound policies—includes:
• Government Effectiveness: perceptions of the quality of public services, the civil service, policy formulation and implementation, and the credibility of a government’s decisions; and
• Regulatory Quality: perceptions of the government’s ability to design and implement policies and regulations that promote private sector development. The third aspect—respect for institutions that govern economic and social interactions— comprises:
• Rule of Law: perceptions of the extent to which agents respect and follow the rules of society, including contract enforcement, property rights, the police, courts, and the likelihood of crime and violence; and
• Control of Corruption: perceptions of the extent to which public power is used for private gain, including both petty and grand corruption, as well as capture of the state by elites and private interests.

Source: World Bank, 2025, “The Worldwide Governance Indicators: Revised Methodology for Measuring Governance Using Perception Data December 2025”, link.

Note: For brevity, this post uses institutions, governance and the WGI interchangeably, in full knowledge that they aren’t the same thing. The WGI is the World Bank’s conceptualization of a particular set of institutions it considers useful, interesting and/or relevant. The WGI is not meant to be an all-encompassing measure suited to every use case, it’s just the most suitable measure for this analysis. The scores are also proxies for institutional quality rather than measures of it. Proxies come with the territory when trying to measure hard-to-define things like institutions and/or governance, as direct measurement is unavailable or unsuitable, particularly for cross-country comparisons.

Figure: Worldwide Governance Indicator Dimensions

The WGI attempts to measure governance across six dimensions

Source: my exceptional design skills

FYI: Daniel Kaufmann has also publicly responded to criticism of the WGI, which has resulted in several iterations of the methodology over time.

Project Setup and Data Cleaning

Once again, the WGI data can be downloaded here. Because the code analyzes all WGI indicators, fnc_read_wgi_sheets() load data from each sheet and combine it into a single dataframe.

Code: Plotting functions

#load the packages we'll probably need
library(tidyverse)
library(readxl)
library(janitor)
library(qgraph)
library(rstatix)   
library(broom)     
#Note: ppcor is also required, but called directly to avoid conflicts
#Install with install.packages("ppcor") if needed.


#define path
ref_wgi_path <- "./Data/wgidataset_with_sourcedata-2025.xlsx"

# sheet code -> output column name
ref_wgi_indices <- c(
  va = "wgi_voice",     # Voice and Accountability
  pv = "wgi_polstab",   # Political Stability & Absence of Violence
  ge = "wgi_goveff",    # Government Effectiveness
  rq = "wgi_regqual",   # Regulatory Quality
  rl = "wgi_rol",       # Rule of Law
  cc = "wgi_corrupt"    # Control of Corruption
)

# Read one WGI sheet and return iso3c, year, and the single renamed estimate.
# Renaming per sheet avoids the shared 'governance_estimate...' column colliding
# on join. Fails loudly if a sheet's estimate column is named differently.
fnc_read_wgi_sheet <- function(sheet_code, col_name, path) {
  read_excel(path, sheet = sheet_code) |>
    clean_names() |>
    rename(iso3c = economy_code,
           !!col_name := governance_estimate_approx_2_5_to_2_5) |>
    select(iso3c, year, all_of(col_name))
}

# read all six, then join on the common keys
dta_wgi_2025 <- ref_wgi_indices |>
  imap(\(col_name, sheet_code) fnc_read_wgi_sheet(sheet_code, col_name, ref_wgi_path)) |>
  reduce(full_join, by = c("iso3c", "year"))

#Add metadata (country + year-specific income bracket) taken once, from any sheet
tmp_dta_wgi_meta <- read_excel(ref_wgi_path, sheet = "rl") |>
  clean_names() |>
  rename(iso3c = economy_code, country = economy_name) |>
  select(iso3c, year, country, income_classification)

#add meta data
dta_wgi_2025 <- dta_wgi_2025 |>
  left_join(tmp_dta_wgi_meta, by = c("iso3c", "year")) |>
  relocate(iso3c, country, year, income_classification)

#drop temporary objects (ref_ objects are kept so the chunk can be re-run)
rm(tmp_dta_wgi_meta)

#define index column names and labels:
ref_wgi_cols   <- c("wgi_voice", "wgi_polstab", "wgi_goveff",
                    "wgi_regqual", "wgi_rol", "wgi_corrupt")

ref_wgi_labels <- c(wgi_voice   = "Voice & accountability",
                    wgi_polstab = "Political stability",
                    wgi_goveff  = "Govt effectiveness",
                    wgi_regqual = "Regulatory quality",
                    wgi_rol     = "Rule of law",
                    wgi_corrupt = "Control of corruption")

#shortened labels  
ref_wgi_short  <- c(wgi_voice   = "Voice & Acc.",
                    wgi_polstab = "Pol. stability",
                    wgi_goveff  = "Gov. effect.",
                    wgi_regqual = "Reg. quality",
                    wgi_rol     = "Rule of law",
                    wgi_corrupt = "Corruption")

# assumptions
ref_min_years <- 20    # minimum years of data before a country average is used
ref_pcor_cut  <- 0.05  # partial correlations below this aren't drawn in the network

The rule of law as part of an institutional portfolio

To demonstrate what a statistical minefield analyzing institutions is, we’ll start by applying the same technique that I criticized earlier: pairwise correlations. With the idea being to demonstrate that different measures of institutional strength tend to agree with one another. In addition to this being what you might expect if institutions were endogenously determined, interdependent and/or WGI dimensions measured something similar. It also illustrates why you can’t rely on associations alone to cleanly establish the causal connection between a particular institution and economic and social outcomes.

For comparing the shared strength of institutions we’ll use a country’s average score for each WGI dimension from 1996 to 2024.3 Although this reduces our sample size and ignores a variety of country-specific factors that might influence institutions, it provides a simple way to reduce year-to-year noise and focus on between-country effects. It’s also unlikely we’ll lose too much information, given institutions move slowly from year-to-year. But, averaging cuts both ways too: as stripping out year-to-year noise also strips out measurement error, making correlations mechanically higher within any single year.

Reflecting that WGI dimensions might just be unrelated series that are trending together, we’ll also examine pairwise correlations for annual changes in WGI scores across dimensions. Aside from checking whether the level associations are just shared trends, dimensions moving together year-to-year would be consistent with institutions being connected to one another. For instance, if the rule of law supports efforts to fight corruption, we might expect the quality of these two institutions to move in the same direction.

However, the use of averages for comparing levels and raw values for first-differences also warrants caution when making comparisons: As while the two describe the same countries, the different units of analysis and number of observations aren’t directly comparable.

Code: calculate country averages and first differences

#Country averages (levels).
#Countries are only kept where every dimension has at least ref_min_years of
#data
dta_wgi_country <- dta_wgi_2025 |>
  group_by(iso3c) |>
  filter(if_all(all_of(ref_wgi_cols), \(x) sum(!is.na(x)) >= ref_min_years)) |>
  summarise(across(all_of(ref_wgi_cols), \(x) mean(x, na.rm = TRUE)),
            nmb_years = n(),
            .groups = "drop")


#how many countries survive the coverage filter
chk_wgi_country_nmb <- nrow(dta_wgi_country)
chk_wgi_country_nmb

#first difference column names
ref_wgi_diff_cols <- paste0("d_", ref_wgi_cols)

#work out which rows sit next to a consecutive year
#(the WGI was published every other year before 2002, so we can't just lag blindly)
dta_wgi_gaps <- dta_wgi_2025 |>
  arrange(iso3c, year) |>
  group_by(iso3c) |>
  mutate(yr_gap = year - lag(year)) |>
  ungroup()

#then take the year-on-year change, but only where the gap is a single year
dta_wgi_diff <- dta_wgi_gaps |>
  group_by(iso3c) |>
  mutate(across(all_of(ref_wgi_cols),
                \(x) if_else(yr_gap == 1, x - lag(x), NA_real_),
                .names = "d_{.col}")) |>
  ungroup() |>
  select(iso3c, year, all_of(c(ref_wgi_cols, ref_wgi_diff_cols)))

#Complete-case frames used for BOTH the zero-order and partial correlations.
dta_wgi_lvl_cc <- dta_wgi_country |>
  select(all_of(ref_wgi_cols)) |>
  drop_na()

dta_wgi_dif_cc <- dta_wgi_diff |>
  select(iso3c, all_of(ref_wgi_diff_cols)) |>
  drop_na()

#sample sizes 
chk_sample_sizes <- tibble(
  series       = c("WGI Levels", "WGI First Differences"),
  nmb_obs      = c(nrow(dta_wgi_lvl_cc), nrow(dta_wgi_dif_cc)),
  nmb_country  = c(nrow(dta_wgi_lvl_cc), n_distinct(dta_wgi_dif_cc$iso3c)),
  unit         = c("country", "country-year")
)

Note: The WGI methodology has been updated in 2025 to improve comparability of governance scores over time.

Code: Correlation analysis

{r}
#create a function for calculating correlations 
fnc_cor_p <- function(dta, cols) {
  dta |>
    select(all_of(cols)) |>
    cor_test(vars = cols) |>
    rename(x = var1, y = var2) |>
    filter(match(x, cols) < match(y, cols)) |>
    mutate(across(c(x, y), \(v) sub("^d_", "", v))) |>
    select(x, y, cor, p)
}

#calculate association(s) on the shared complete-case frames
sum_wgi_cor_level <- fnc_cor_p(dta_wgi_lvl_cc, ref_wgi_cols)
sum_wgi_cor_diff  <- fnc_cor_p(dta_wgi_dif_cc, ref_wgi_diff_cols)

#Holm adjustment for testing 15 pairs at once. Note the p-values for the
#differences assume ~5,000 independent observations when they are really ~200
#countries observed repeatedly - so significance there is close to guaranteed
#and the effect sizes are what matter.
sum_wgi_cor_level <- sum_wgi_cor_level |>
  mutate(sig = p.adjust(p, "holm") < 0.05,
         r2  = cor^2)

sum_wgi_cor_diff <- sum_wgi_cor_diff |>
  mutate(sig = p.adjust(p, "holm") < 0.05,
         r2  = cor^2)

This code chunk just creates a basic bubble plot for visualizing the correlation matrix. This is perhaps more verbose than it needs to be, but it felt fair to give Claude the satisfaction of trying to apply my style guidelines.

Code: Plotting functions

#Bubble plot for a correlation matrix. value_col lets the same function handle
#zero-order (cor) and partial (pcor) tables.
fnc_plot_bubbles <- function(dta, value_col = "cor", wrap = 12) {
  
  ref_lbl <- setNames(str_wrap(ref_wgi_labels, wrap), ref_wgi_cols)
  
  dta <- dta |>
    rename(value = all_of(value_col)) |>
    mutate(x   = factor(x, ref_wgi_cols),
           y   = factor(y, ref_wgi_cols),
           lbl = ifelse(sig, sprintf("%.2f", value), sprintf("(%.2f)", value)))
  
  ggplot(dta, aes(x, y)) +
    geom_point(aes(size = abs(value), fill = value), shape = 21, colour = "white") +
    geom_text(aes(label = lbl, colour = abs(value) > 0.55), size = 3, fontface = "bold") +
    scale_fill_gradient2(low = "#922C40", mid = "#FFFFFF", high = "#16AF8E",
                         midpoint = 0, limits = c(-1, 1), name = "Correlation") +
    scale_colour_manual(values = c(`TRUE` = "white", `FALSE` = "#121212"), guide = "none") +
    scale_size_area(max_size = 16, limits = c(0, 1), guide = "none") +
    scale_x_discrete(labels = ref_lbl) +
    scale_y_discrete(labels = ref_lbl, limits = rev) +
    coord_fixed() +
    labs(x = NULL, y = NULL) +
    theme_minimal()
}

# Convert the long pairwise tibble back into a symmetric matrix for qgraph
fnc_pcor_matrix <- function(dta_pcor, col) {
  tmp_mat <- matrix(0, length(ref_wgi_cols), length(ref_wgi_cols),
                    dimnames = list(ref_wgi_cols, ref_wgi_cols))
  tmp_mat[cbind(dta_pcor$x, dta_pcor$y)] <- dta_pcor[[col]]
  tmp_mat[cbind(dta_pcor$y, dta_pcor$x)] <- dta_pcor[[col]]
  tmp_mat
}

#Network of direct links. Edges below `cut` in absolute size are NOT DRAWN, so a
#missing line means "smaller than the threshold", not "zero".
fnc_plot_network <- function(dta_pcor, cut = ref_pcor_cut, layout = "spring", seed = 123) {
  
  set.seed(seed)
  
  mat_pcor <- fnc_pcor_matrix(dta_pcor, "pcor")
  chk_sig  <- fnc_pcor_matrix(dta_pcor, "sig") == 1
  
  qgraph(mat_pcor,
         layout = layout, minimum = cut, maximum = 1, esize = 9, fade = TRUE,
         posCol = "#16AF8E", negCol = "#922C40",
         color = "white", border.color = "#E5E7EB", border.width = 2,
         labels = ref_wgi_short[ref_wgi_cols], label.cex = 0.8,
         label.color = "#121212", vsize = 11, shape = "circle",
         lty = ifelse(chk_sig, 1, 2),
         edge.labels = ifelse(chk_sig, sprintf("%.2f", mat_pcor),
                              sprintf("(%.2f)", mat_pcor)),
         edge.label.cex = 0.75, edge.label.bg = "white",
         edge.label.color = "#121212",  mar = rep(5, 4))
}

Strong institutions coincide with one another

In this section, pairwise correlations are used to test whether the strength of individual institutions tends to occur together.

It’ll come as no surprise that all of the WGI’s dimensions are strongly associated with each other, which suggests that on average a country scoring highly on one dimension probably scores highly on the others too. For the WGI’s Rule of Law measure, the pairwise associations are strong across the board, with the implied R squared statistic suggesting it accounts for somewhere between 65 and 90 percent of the variation in the other dimensions. Political Stability has the weakest associations, accounting for somewhere between 50 and 75 percent.

Figure: WGI pairwise correlations (levels)

The strength of institutions are strongly associated with one another across WGI dimensions

The network plot presents the same pairwise associations once the influence of the other dimensions is accounted for. Where partial correlations fall well below the zero-order correlations shown above, it points to dimensions carrying overlapping information. For the Rule of Law, the network diagram presents much weaker associations with other dimensions to the figure above, with moderate associations remaining for Voice and Accountability, Control of Corruption and Political Stability. Although it’s best not to take the implications of this analysis too far, one interpretation of these results that makes intuitive sense is that different institutions relate to one another differently. And sometimes this relationship might be indirect, such as government effectiveness indirectly influencing the rule of law via corruption controls.

Figure: WGI pairwise partial correlations (levels)

Pairwise correlations are lower once the influence of other WGI dimensions are accounted for

Code: Plotting pairwise correlations at levels

#produce and display the correlation matrix
plt_wgi_cor_level <- fnc_plot_bubbles(sum_wgi_cor_level)

plt_wgi_cor_level

#network plot for partial correlations 
plt_wgi_net_level <- fnc_plot_network(sum_wgi_pcor_level)

Institutions are connected to one another over time (albeit, loosely)

This section explores pairwise correlations between year-to-year movements of WGI dimension to explore evidence for institutions being interdependent.

Focusing on year-to-year movements in WGI dimensions rather than levels results in smaller correlation statistics across the board, which is to be expected given taking the first difference naturally inflates the influence of noise and measurement errors. Still, the associations point to a similar picture to associations at the levels: a positive and statistically significant association between all dimensions, with Political Stability being the weakest.

Having said this, statistical significance just suggests that shared movement between series isn’t zero, not that it’s particularly interesting. Added to this, given the first differences pool values so they are considered independent (despite being from the same countries), significance is close to guaranteed, which is yet another reason they should be interpreted with caution.

The implied explanatory power is also probably what matters more and it’s generally modest. For instance, the R squared statistic for Political stability suggests it can explain between 1 to 4 percent of the year-to-year variation in other dimensions. Whereas the Rule of Law, which holds the strongest pairwise explanatory power across dimensions (setting political stability aside), only explains somewhere between 9 to 17 percent of the variation, which while being nothing to sneeze at, still leaves a lot of unexplained movement.

Figure: WGI pairwise correlations (year-to-year changes)

Year-to-year movements in institutional quality are statistically associated with one another

Note: Associations between first differences are expected to be weaker for several reasons: institutions are slow moving; interdependence between institutions may operate with lags or in a non-linear fashion; and as differencing strips out persistent components of each series while retaining measurement errors intact, the signal to noise ratio is likely to be lower. The authors of the WGI also warn against taking too much stock of year-to-year changes in scores.

The second plot once again presents how year-to-year movements are associated with each other, but after accounting for influences outside the examined pair. Once again, the Rule of Law holds the strongest association with other dimensions, but its association with Control of Corruption is much lower, suggesting that much of the pairwise association between year-to-year changes above relate to other dimension outside the pair.

Figure: WGI pairwise partial correlations (year-to-year changes)

Pairwise correlations are reduced, but in most cases remain statistically significant, once the influence of other dimensions is accounted for

Code: plot first difference associations

#correlaton plot for correlations @ first difference
plt_wgi_cor_diff <- fnc_plot_bubbles(sum_wgi_cor_diff)

plt_wgi_cor_diff

#network plot for partial correlatons @ first difference
plt_wgi_net_diff <- fnc_plot_network(sum_wgi_pcor_diff,
                                     layout = plt_wgi_net_level$layout)

Six dimensions, one signal(?)

This section uses Principal Components Analysis (PCA) to explore whether all six dimensions of the WGI measure, or are caused by, the same thing.

A strong critique of the WGI is that dimensions more or less measure the same thing. From a conceptual standpoint this might have some weight, due to the perceived quality of one dimension influencing perceptions of another, the use of overlapping data sources and the fact that the definitions are not precisely defined. But, as somebody that works a lot with composite indices, I’d say this comes with the territory and I’m happy to leave it to the experts to argue among themselves.

However, from a statistical standpoint this might matter a lot, as it determines whether a measured link between the Rule of Law and outcomes like prosperity can be meaningfully interpreted as telling us anything about the Rule of Law.

This is explored in the code below by applying PCA to the WGI levels and first differences. PCA attempts to collapse collinear variables into a set of principal components (PC) that explain as much variability as possible. If the WGI’s dimensions are measuring or being driven by something common, we might expect a small number of PCs will explain a disproportional share of the variance.

Code: apply PCA to levels and first differences

#Fit a scaled PCA and pull out variance shares and PC1 loadings.
#Scaling so each index counts equally regardless of how spread out it is.
#Note the levels PCA runs on country averages (one row per country) while the
#differences PCA runs on pooled country-years - the shares are not measured on
#the same unit of analysis.
fnc_pca_tidy <- function(dta, cols, series) {
  
  dta_input <- dta |>
    select(all_of(cols)) |>
    drop_na()
  
  mod <- prcomp(dta_input, scale. = TRUE)
  
  #share of variance picked up by each component
  sum_var <- mod |>
    tidy(matrix = "eigenvalues") |>
    transmute(series, pc = paste0("PC", PC), var_pct = percent,
              nmb_obs = nrow(dta_input))
  
  #how strongly each index marks the first component
  sum_load <- mod |>
    tidy(matrix = "rotation") |>
    filter(PC == 1) |>
    transmute(series, index = sub("^d_", "", column), loading = abs(value))
  
  list(var = sum_var, load = sum_load)
}

ref_pca_series <- c("WGI Levels", "WGI First Differences")

tmp_pca <- list(fnc_pca_tidy(dta_wgi_lvl_cc, ref_wgi_cols,      ref_pca_series[1]),
                fnc_pca_tidy(dta_wgi_dif_cc, ref_wgi_diff_cols, ref_pca_series[2]))

sum_wgi_pca_var <- tmp_pca |>
  map("var") |>
  list_rbind() |>
  mutate(series = factor(series, ref_pca_series))

sum_wgi_pca_load <- tmp_pca |>
  map("load") |>
  list_rbind() |>
  mutate(series = factor(series, ref_pca_series))

#drop temporary objects
rm(tmp_pca)

The first plot presents the share of variance explained by each principal component. In the case of the the WGI’s levels, the PCA indicates that almost 90 percent of the measured variance in governance can be explained by a single principal component. Indicating that the majority of information presented by the six WGI dimensions could be efficiently described by a single measure. A result that supports the idea that either the WGI is measuring a similar thing and/or that a common factor is driving all six dimensions.

Both PCAs point in a similar direction, although the first differences are much less dramatic: with the first component picks up around 40 percent of the variance, against nearly 90 percent for the levels. Bear in mind these two numbers aren’t measured on the same unit of analysis, making them not directly comparable (i.e. approximately 200 country averages vs 5,000+ first differences). It’s therefore best not to read too much into comparisons, particularly given the authors of the WGI explicitly warn against analyzing year-to-year score movements.

Figure: Variance explained by principal component

The majority of variance can be explained by a single PC for WGI levels, while a larger number of PCs are required to provide a sufficient explanation of variability for first differences

Code: how concentrated is the common dimension?

{r}
plt_wgi_pca_scree <- sum_wgi_pca_var |>
  ggplot(aes(pc, var_pct)) +
  geom_col(fill = "#1E298D", width = 0.7) +
  geom_text(aes(label = scales::percent(var_pct, accuracy = 1)),
            vjust = -0.4, size = 3, colour = "grey40") +
  scale_y_continuous(labels = scales::percent, limits = c(0, 1)) +
  facet_wrap(~ series) +
  labs(x = NULL, y = "Variance explained (%)") +
  theme_minimal(base_size = 10) +
  theme(panel.grid.major.x = element_blank(),
        panel.grid.minor   = element_blank(),
        strip.text         = element_text(face = "bold", hjust = 0))

plt_wgi_pca_scree

The plot below presents PCA loadings, which measure how strongly each dimension contributes to PC1. Higher loadings mean a dimension is more closely tied to the component, and so shares more with the others. That the dimensions carry comparable loadings at both the levels and the first differences is again consistent with them measuring, or being influenced by, something common.

Political Stability sits lower than the rest, but I’d be reluctant to take my interpretation of that too far given how much variance remains unexplained. It’s also roughly what you might expect of a dimension intended to capture shocks rather than gradual shifts, or one behaving non-linearly (questions better answered in a separate post).

Figure: PCA Loadings

PCA loadings for both the levels and first differences are generally evenly spread

Code: PCA loadings

plt_wgi_pca_load <- sum_wgi_pca_load |>
  mutate(index = factor(index, ref_wgi_cols, ref_wgi_labels) |> fct_reorder(loading)) |>
  ggplot(aes(loading, index)) +
  geom_col(fill = "#1E298D", width = 0.65) +
  geom_text(aes(label = sprintf("%.2f", loading)),
            hjust = -0.25, size = 3, colour = "grey40") +
  scale_x_continuous(limits = c(0, 0.7), expand = expansion(c(0, 0.05))) +
  facet_wrap(~ series) +
  labs(x = "Absolute loading on PC1", y = NULL) +
  theme_minimal(base_size = 10) +
  theme(panel.grid.major.y = element_blank(),
        panel.grid.minor   = element_blank(),
        strip.text         = element_text(face = "bold", hjust = 0))

plt_wgi_pca_load

Summing up

Across three views of the same data, the six WGI dimensions behave like parts of one system rather than six separable things. Country averages correlate strongly across every pair; most associations remain even when accounting for the influence of other dimensions; and a single principal component accounts for a significant share of variation in governance scores between countries. Dimensions seem to move together too, albeit weakly, which is while not a slam dunk, does point to WGI scores having some interdependence.

But, it’s worth being clear about what this analysis doesn’t establish: Collinearity isn’t endogeneity and while a system of interdependent institutions might produce this pattern, so would WGI dimensions measuring the same thing.

But, the point of this post isn’t to support the legitimacy of the WGI or any other measure. And the distinction matters little to the point of this post. As whether dimensions move together because they co-determine each other, share a common driver, or share source data, the implications for anyone trying to understand or analyze institutions is the same: they have to be examined as a set as a correlation between one institution and an outcome may strong regardless of whether that pillar is doing the actual work or not.

How AI was used for this post: Claude was used heavily to refine the code and lightly leaned on to improve the accuracy and readability of the text. The former was mainly in an attempt to address Claude’s almost endless array of suggestions for adding more code, while the latter was mainly a result of having stared too long at my own writing.

No doubt errors remain as it’s quite the topic, which means I’ve made quite a number of edits to both the text and code. Feel free to contact me here.

Additional Note:

The motivation for writing this was the absence of descriptive analysis on the links between law and order and economic growth aimed at a general audience. Keeping the code in the post has probably cost it some readability, but my reasoning for sharing the analysis so openly was to make it easier for others to build on this post to fill the many gaps out there (including me in my future posts).

The inspiration for this series came from work I completed in 2025 for the Bingham Centre for the Rule of Law and the Law Society of England and Wales. A paper based on the work summarizing research on the topic is available here.

1. For a great summary of the research, see Dr Lopez-Gomez, L. (2026, The Rule of Law and the Institutional Roots of Economic Performance. The British Institute of International and Comparative Law, (link). Also see: Haggard, S., MacIntyre, A. and Tiede, L., 2008. The rule of law and economic development. Annu. Rev. Polit. Sci., 11(1), pp.205-234; and Besley, T., Bogart, D., Chapman, J. and Nuno, P., 2025. Justices of the peace: Legal foundations of the industrial revolution (No. 20214). Centre for Economic Policy Research. ︎
2. One of the better attempts I’ve seen to empirically test a causal connection between the two: Besley, T., Bogart, D., Chapman, J. and Nuno, P., 2025. Justices of the peace: Legal foundations of the industrial revolution (No. 20214). Centre for Economic Policy Research. ︎

The post Institutions Are a Package Deal: What a Correlation Can’t Tell You About the Rule of Law appeared first on Giles.

Read it in: Español.

Population censuses are a key source of information for understanding the composition of populations and how they change across regions within a country. They provide essential evidence for research, the design and evaluation of public policies, and informed decision-making.

Working with census data typically involves several preliminary steps: identifying what information is available, downloading the data, organizing it, and preparing it for analysis. When this information is scattered across multiple sources and provided in different formats, the process can become complex.

ARcenso is an R package developed as part of the rOpenSci Champions Program to facilitate access to Argentine census data and simplify its analysis.

In this article, the ARcenso development team demonstrates how to work with the historical census data available through the package in a simple, reproducible way.

The goal is to conduct a demographic analysis using official national data from the 1970 and 1980 censuses, combining visualisations and tables to explore the structure of the population—that is, how it is distributed by age and sex. This type of analysis helps characterise demographic changes over time, providing key information for research and the design of public policies.

Getting started

To run the code examples from this article on your own computer, make sure you have the following packages installed:

# If you don't have pak installed
install.packages("pak")

# Install ARcenso from GitHub
pak::pkg_install("soyandrea/arcenso")

# Install the required packages from CRAN
install.packages(c("dplyr", "tidyr", "ggplot2", "gt"))

Next, we load the packages needed to work with the census data and build indicators for the analysis:

library(arcenso) # access census data
library(dplyr) # data manipulation
library(tidyr) # data tidying and transformation
library(ggplot2) # data visualisation
library(gt) # table formatting

Accessing census data

Now that the required packages are installed and loaded, we can use ARcenso to access census data from the 1970 and 1980 censuses. The package allows users to retrieve information by census year, topic, and geographic area of interest, using the official geographic codes defined by the National Institute of Statistics and Censuses (INDEC) of Argentina .

The check_repository() function helps identify the datasets available in the package. The topic argument specifies the subject of interest (for example, population structure), while geo_code identifies the geographic area.

If you are unsure which values are available for these arguments, you can explore the package metadata. The geo_metadata object contains the available geographic areas and their corresponding codes, while census_metadata provides information about the topics and tables included in the package.

In this example, we use population structure data (topic = "estructura") for the entire country (geo_code = "00"). As a first step, we use check_repository() to see which datasets are available.

check_repository(topic = "estructura", geo_code = "00")

# A tibble: 4 × 3
id_cuadro anio titulo
<chr> <dbl> <chr>
1 1970_00_estructura_01 1970 Cuadro 1. Total del país. Población total, por gr…
2 1980_00_estructura_01 1980 Cuadro G3. Centros urbanos según tamaño y poblaci…
3 1980_00_estructura_02 1980 Cuadro G1. Total del país. Población total según …
4 1980_00_estructura_03 1980 Cuadro G2. Total del país. Población según sexo y…

The output of check_repository() lists the available years for the selected combination of topic and geographic area, together with the corresponding table identifiers (IDs) and their titles. Based on this information, we select the tables 1970_00_estructura_01 and 1980_00_estructura_03, which contain tabulations of the country’s total population by sex and age group and will be used throughout this analysis.

Tip: You can also explore the available data interactively using arcenso_app(). Once you find the table you need, simply copy its ID and use it in your analysis.

Data preparation

Once the tables of interest have been identified, we use get_census() to import them directly into the R session using their table identifiers.

As the structures of the two tables are not identical, we perform a series of transformations to harmonise the variables and create a consistent basis for comparison across censuses. In particular, we recode the age categories into five-year groups to ensure a common structure and add a column identifying the census year.

1970 Census

For the 1970 census, the data are already grouped into five-year age groups, so we only need to adjust the format of some labels and keep the variables relevant to the analysis.

# Country Total
poblacion_1970 <- get_census(id = "1970_00_estructura_01")

pob_1970 <- poblacion_1970 |>
 filter(sexo != "Total") |>
 mutate(
 censo = 1970,
 grupo_de_edad = case_when(
 grupo_de_edad == "0-4" ~ "00-04",
 grupo_de_edad == "5-9" ~ "05-09",
 TRUE ~ grupo_de_edad
 )
 ) |>
 rename(grupo_edad = grupo_de_edad) |>
 select(censo, sexo, grupo_edad, poblacion)

1980 Census

In contrast, the 1980 census table reports age as single years and includes additional information that is not relevant to this analysis. We therefore select the categories of interest and aggregate the data into five-year age groups to make them comparable with the 1970 census.

poblacion_1980 <- get_census(id = "1980_00_estructura_03")

pob_1980 <- poblacion_1980 |>
 filter(urbano_rural == "Total", sexo != "Total", edad != "Total") |>
 mutate(
 censo = 1980,
 edad_num = ifelse(edad == "85 y más", 85, as.numeric(edad)),
 grupo_edad = case_when(
 edad_num %in% c(0:4) ~ "00-04",
 edad_num %in% c(5:9) ~ "05-09",
 edad_num %in% c(10:14) ~ "10-14",
 edad_num %in% c(15:19) ~ "15-19",
 edad_num %in% c(20:24) ~ "20-24",
 edad_num %in% c(25:29) ~ "25-29",
 edad_num %in% c(30:34) ~ "30-34",
 edad_num %in% c(35:39) ~ "35-39",
 edad_num %in% c(40:44) ~ "40-44",
 edad_num %in% c(45:49) ~ "45-49",
 edad_num %in% c(50:54) ~ "50-54",
 edad_num %in% c(55:59) ~ "55-59",
 edad_num %in% c(60:64) ~ "60-64",
 edad_num %in% c(65:69) ~ "65-69",
 edad_num %in% c(70:74) ~ "70-74",
 edad_num %in% c(75:79) ~ "75-79",
 edad_num %in% c(80:84) ~ "80-84",
 TRUE ~ "85 y más"
 )
 ) |>
 select(censo, sexo, grupo_edad, poblacion)

Building the integrated database

Once both tables have been processed, we combine them into a single dataset and define the final structure of the variables, preparing the data for visualisation and the analysis of demographic indicators.

# We combine both geographic areas and prepare the variables for use
poblacion <-
 bind_rows(
 pob_1970,
 pob_1980
 ) |>
 mutate(
 poblacion = as.numeric(poblacion),
 sexo = factor(sexo, levels = c("Varones", "Mujeres")),
 grupo_edad = factor(
 grupo_edad,
 levels = c(
 "00-04",
 "05-09",
 "10-14",
 "15-19",
 "20-24",
 "25-29",
 "30-34",
 "35-39",
 "40-44",
 "45-49",
 "50-54",
 "55-59",
 "60-64",
 "65-69",
 "70-74",
 "75-79",
 "80-84",
 "85 y más"
 )
 )
 )

Population structure

Now that we have a dataset combining both censuses, with population totals organised by census year, sex, and age group, we can consistently compare the population composition across the two census years.

head(poblacion)


 # A tibble: 6 × 4
 censo sexo grupo_edad poblacion
 <dbl> <fct> <fct> <dbl>
 1 1970 Varones 00-04 1196950
 2 1970 Mujeres 00-04 1158350
 3 1970 Varones 05-09 1163050
 4 1970 Mujeres 05-09 1133950
 5 1970 Varones 10-14 1114300
 6 1970 Mujeres 10-14 1086850

Organising the population by sex and five-year age groups allows us to analyse its structure using a variety of visualisations. In this case, we use a population pyramid, which makes it easier to interpret both dimensions together.

Population pyramid

A population pyramid displays the distribution of the population by age and sex simultaneously, typically showing males on the left and females on the right. By representing the proportion of people in each five-year age group, it provides a clear visual summary of the population structure.

In this case, we calculated the relative distribution of the population within each census, allowing us to compare the population structure between 1970 and 1980 regardless of total population size. The shape of the pyramid also allows us to identify broad demographic patterns, such as a higher concentration of younger age groups or a relatively older population profile.

We then calculated the relative distribution of the population in each census and constructed a population pyramid to compare the population structure of the two censuses.

# Dataset for the Population Pyramid
piramide <- poblacion |>
 group_by(censo, sexo) |>
 mutate(
 poblacion_rel = if_else(
 sexo == "Varones",
 -poblacion / sum(poblacion),
 poblacion / sum(poblacion)
 )
 ) |>
 ungroup()

# Comparison Pyramid
piramide |>
 ggplot(aes(x = poblacion_rel, y = grupo_edad, fill = sexo)) +
 geom_col() +
 facet_wrap(~censo, ncol = 2) +
 scale_fill_manual(values = c("#00f59b", "#7014f2")) +
 scale_x_continuous(
 labels = function(x) paste0(abs(round(x * 100, 1)), "%"),
 limits = c(-0.15, 0.15),
 breaks = seq(-0.15, 0.15, by = 0.05)
 ) +
 labs(
 title = "Gráfico 1. Estructura de la población por sexo y grupo quinquenal de edad.",
 subtitle = "Argentina. Años 1970 y 1980",
 x = "Porcentaje",
 y = "Grupo quinquenal de edad",
 caption = "Fuente: INDEC, Censo Nacional de Población 1970 y 1980. Procesado con ARcenso.",
 fill = "Sexo"
 ) +
 theme_bw() +
 theme(
 legend.position = "bottom",
 strip.text = element_text(face = "bold", size = 12)
 )

Both censuses show a young population structure, with a high concentration of people in the younger age groups. However, by 1980, a slight shift towards the adult age groups is already apparent, suggesting the early stages of population aging.

Construction of demographic indicators

While the population pyramid provides an overall view of the population structure, demographic indicators offer summary measures that allow these patterns to be quantified and compared more precisely. In this section, we calculate two commonly used indicators to complement the visual analysis.

Aging index

The aging index compares the number of older adults (aged 65 years and over) with the number of children (aged 0–14 years). It provides a simple way to see whether the population is weighted more toward younger age groups or older age groups.

envejecimiento <- poblacion |>
 group_by(censo) |>
 summarise(
 poblacion_0a14 = sum(poblacion[
 grupo_edad %in% c("00-04", "05-09", "10-14")
 ]),
 poblacion_65ymas = sum(poblacion[
 grupo_edad %in% c("65-69", "70-74", "75-79", "80-84", "85 y más")
 ]),
 indice = round(poblacion_65ymas / poblacion_0a14 * 100, 0)
 )


gt(envejecimiento) |>
 tab_header(
 title = "Comparación del índice de envejecimiento",
 subtitle = "Argentina. Años 1970 y 1980"
 ) |>
 tab_spanner(
 label = "Población",
 columns = c(poblacion_0a14, poblacion_65ymas)
 ) |>
 fmt_number(
 columns = c(poblacion_0a14, poblacion_65ymas),
 decimals = 0,
 sep_mark = "."
 ) |>
 cols_label(
 poblacion_0a14 = "0 a 14 años",
 poblacion_65ymas = "65 años y más",
 indice = "Indice"
 ) |>
 tab_source_note(
 source_note = md(
 "**Fuente:** elaboración propia en base a datos de INDEC (Censos Nacionales de Población 1970 y 1980)."
 )
 )

The increase in the aging index reflects a change in the population age structure between 1970 and 1980, indicating a larger relative share of people aged 65 years and over compared with the younger population. This suggests a gradual aging of the population.

Female-to-Male ratio

This indicator shows how many women there are for every 100 men in a specific population group. Here, we calculate it for people aged 60 years and over, where differences between women and men tend to become more pronounced.

feminidad <- poblacion |>
 filter(
 grupo_edad %in% c("60-64", "65-69", "70-74", "75-79", "80-84", "85 y más")
 ) |>
 group_by(censo, grupo_edad, sexo) |>
 summarise(poblacion = sum(poblacion), .groups = "drop") |>
 pivot_wider(names_from = sexo, values_from = poblacion) |>
 mutate(
 indice_feminidad = round(Mujeres / Varones * 100, 0)
 ) |>
 select(-Varones, -Mujeres)


feminidad_plot <- feminidad |>
 pivot_wider(
 names_from = censo,
 values_from = indice_feminidad,
 names_prefix = "censo_"
 ) |>
 ggplot(aes(y = grupo_edad)) +
 geom_segment(
 aes(x = censo_1970, xend = censo_1980, yend = grupo_edad),
 color = "grey85",
 linewidth = 1
 ) +
 geom_point(
 aes(x = censo_1970),
 color = "#ff0f7b",
 size = 3
 ) +
 geom_point(
 aes(x = censo_1980),
 color = "#f89b29",
 size = 3
 ) +
 geom_text(
 aes(x = censo_1970, label = censo_1970),
 hjust = 1.4,
 size = 3
 ) +
 geom_text(
 aes(x = censo_1980, label = censo_1980),
 hjust = -0.4,
 size = 3
 ) +
 labs(
 x = "Mujeres por cada 100 varones",
 y = "Grupo de edad",
 title = "Cambio en el índice de feminidad de la población de 60 años y más",
 subtitle = "Argentina, Años 1970 y 1980",
 caption = "Fuente: INDEC, Censo Nacional de Población 1970 y 1980. Procesado con ARcenso."
 ) +
 theme_minimal()

feminidad_plot

The female-to-male ratio shows a higher proportion of women in the older age groups. This difference becomes more pronounced between 1970 and 1980 across all the age groups analysed, reflecting patterns of higher female life expectancy.

What does this analysis tell us?

In demographic analysis, the most time-consuming part is often not calculating indicators or creating visualisations, but the work that comes beforehand: identifying the available data for each census, understanding how it is organised, and building a consistent framework that enables comparisons across census years.

In this example, that process involved locating the relevant tables, understanding the structural differences between the 1970 and 1980 censuses, and harmonising the variables based on the available information to create a comparable dataset. This is precisely what ARcenso is designed to simplify. The check_repository() function helps identify the available tables, get_census() imports the data into R in a structured format, and arcenso_app() provides an interactive way to explore the repository.

Once the data have been organised, the analysis becomes more accessible, reproducible, and easier to extend to new questions.

ARcenso is still a work in progress. Future developments will incorporate additional census years and continue expanding the package’s analytical capabilities. The package was created to make Argentine census data easier to discover, access, and analyse, while building on the collaborative efforts of the open-source community.

If you would like to explore further,check out the vignettes in the package, which include additional examples to work with census data.