CRAN, GitHub

TL;DR: {bdlnm} brings Bayesian Distributed Lag Non-Linear Models (B-DLNMs) to R using INLA, allowing to model complex DLNMs, quantify uncertainty, and produce rich visualizations.

install.packages("bdlnm")
library(bdlnm)

install.packages(
  "INLA",
  repos = c(
    getOption("repos"),
    INLA = "https://inla.r-inla-download.org/R/stable"
  ),
  dep = TRUE
)

# DLNMs and splines
library(dlnm)
library(splines)

# Data manipulation
library(dplyr)
library(reshape2)
library(stringr)
library(lubridate)

# Visualization
library(ggplot2)
library(gganimate)
library(ggnewscale)
library(patchwork)
library(scales)
library(plotly)

# Tables
library(gt)

# Execution time
library(tictoc)

col_mort <- "#2f2f2f"
col_temp <- "#8e44ad"

# Scaling parameters
a <- (max(london$mort_75plus) - min(london$mort_75plus)) /
  (max(london$tmean) - min(london$tmean))
b <- min(london$mort_75plus) - min(london$tmean) * a

p <- ggplot(london, aes(x = yday(date))) +
  geom_line(
    aes(y = a * tmean + b, color = "Mean Temperature"),
    linewidth = 0.4
  ) +
  geom_line(
    aes(y = mort_75plus, color = "Daily Mortality (+75 years)"),
    linewidth = 0.4
  ) +
  facet_wrap(~year, ncol = 3) +
  scale_y_continuous(
    name = "Daily Mortality (+75 years)",
    breaks = seq(0, 225, by = 50),
    sec.axis = sec_axis(
      name = "Mean Temperature (°C)",
      transform = ~ (. - b) / a,
      breaks = seq(-10, 30, by = 10)
    )
  ) +
  scale_x_continuous(
    breaks = yday(as.Date(paste0(
      "2000-",
      c("01", "03", "05", "07", "09", "11"),
      "-01"
    ))),
    labels = c("Jan", "Mar", "May", "Jul", "Sep", "Nov"),
    expand = c(0.01, 0)
  ) +
  scale_color_manual(
    values = c(
      "Daily Mortality (+75 years)" = col_mort,
      "Mean Temperature" = col_temp
    )
  ) +
  labs(x = NULL, color = NULL) +
  guides(color = "none") +
  theme_minimal() +
  theme(
    axis.title.y.left = element_text(
      color = col_mort,
      face = "bold",
      margin = margin(r = 8)
    ),
    axis.title.y.right = element_text(
      color = col_temp,
      face = "bold",
      margin = margin(l = 8)
    ),
    axis.text.y.left = element_text(color = col_mort),
    axis.text.y.right = element_text(color = col_temp)
  ) +
  transition_reveal(as.numeric(date))

animate(p, nframes = 300, fps = 10, end_pause = 100)

# Exposure-response and lag-response spline parameters
dlnm_var <- list(
  var_prc = c(10, 75, 90),
  var_fun = "ns",
  lag_fun = "ns",
  max_lag = 21,
  lagnk = 3
)

# Cross-basis parameters
argvar <- list(
  fun = dlnm_var$var_fun,
  knots = quantile(london$tmean, dlnm_var$var_prc / 100, na.rm = TRUE),
  Bound = range(london$tmean, na.rm = TRUE)
)

arglag <- list(
  fun = dlnm_var$lag_fun,
  knots = logknots(dlnm_var$max_lag, nk = dlnm_var$lagnk)
)

# Create crossbasis
cb <- crossbasis(london$tmean, lag = dlnm_var$max_lag, argvar, arglag)

seas <- ns(london$date, df = round(8 * length(london$date) / 365.25))

temp <- round(seq(min(london$tmean), max(london$tmean), by = 0.1), 1)

tictoc::tic()
mod <- bdlnm(
  mort_75plus ~ cb + factor(dow) + seas,
  data = london,
  family = "poisson",
  sample.arg = list(n = 1000, seed = 5243)
)
tictoc::toc()

8.33 sec elapsed

tictoc::tic()
mmt <- optimal_exposure(mod, exp_at = temp)
tictoc::toc()

7.3 sec elapsed

ggplot(as.data.frame(mmt$est), aes(x = mmt$est)) +
  geom_histogram(
    fill = "#A8C5DA",
    bins = length(unique(mmt$est)),
    alpha = 0.6,
    color = "white"
  ) +
  geom_density(
    aes(y = after_stat(density) * length(mmt$est) / length(unique(mmt$est))),
    color = "#2E5E7E",
    linewidth = 1.2,
    adjust = 2 # <-- key change: higher = smoother
  ) +
  geom_vline(
    xintercept = mmt$summary["0.5quant"],
    color = "#2E5E7E",
    linewidth = 1.1,
    linetype = "dashed"
  ) +
  scale_x_continuous(breaks = seq(min(mmt$est), max(mmt$est), by = 0.1)) +
  labs(x = "Temperature (°C)", y = "Frequency") +
  theme_minimal()

cen <- mmt$summary[["0.5quant"]]
tictoc::tic()
cpred <- bcrosspred(mod, exp_at = temp, cen = cen)
tictoc::toc()

6.83 sec elapsed

We can plot the full exposure-lag-response association using a 3-D surface:

matRRfit_median <- cpred$matRRfit.summary[,, "0.5quant"]
x <- rownames(matRRfit_median)
y <- colnames(matRRfit_median)
z <- t(matRRfit_median)

zmin <- min(z, na.rm = TRUE)
zmax <- max(z, na.rm = TRUE)
mid <- (1 - zmin) / (zmax - zmin)

plot_ly() |>
  add_surface(
    x = x,
    y = y,
    z = z,
    surfacecolor = z,
    cmin = zmin,
    cmax = zmax,
    colorscale = list(
      c(0, "#00696e"),
      c(mid * 0.5, "#80c8c8"),
      c(mid, "#f5f0e8"),
      c(mid + (1 - mid) * 0.5, "#c2714f"),
      c(1, "#6b1c1c")
    ),
    colorbar = list(title = "RR")
  ) |>
  add_surface(
    x = x,
    y = y,
    z = matrix(1, nrow = length(y), ncol = length(x)),
    colorscale = list(c(0, "black"), c(1, "black")),
    opacity = 0.4,
    showscale = FALSE
  ) |>
  layout(
    title = "Exposure-Lag-Response Surface",
    scene = list(
      xaxis = list(title = "Temperature (°C)"),
      yaxis = list(title = "Lag", tickvals = y, ticktext = gsub("lag", "", y)),
      zaxis = list(title = "RR"),
      camera = list(eye = list(x = 1.5, y = -1.8, z = 0.8))
    )
  )

The surface reveals two distinct risk regions. Hot temperatures produce a sharp, acute risk concentrated at the first lags, peaking at lag 0 and dissipating rapidly after the first lags. Cold temperatures produce a more modest and gradual increase in the first lags that does not fully disappear at longer lags. Intermediate temperatures near the MMT sit close to the RR = 1 reference plane across all lags.

The differential lag structure observed for heat- and cold-related mortality is consistent with known physiological mechanisms. Heat-related mortality tends to occur rapidly after exposure due to acute physiological stress, whereas cold-related mortality develops more gradually through delayed cardiovascular and respiratory effects, leading to increasing risk over longer lag periods.

matRRfit <- cbind(
  melt(cpred$matRRfit.summary[,, "0.5quant"], value.name = "RR"),
  RR_lci = melt(
    cpred$matRRfit.summary[,, "0.025quant"],
    value.name = "RR_lci"
  )$RR_lci,
  RR_uci = melt(
    cpred$matRRfit.summary[,, "0.975quant"],
    value.name = "RR_uci"
  )$RR_uci
) |>
  rename(temperature = Var1, lag = Var2) |>
  mutate(
    lag = as.numeric(gsub("lag", "", lag))
  )

temps <- cpred$exp_at

p <- ggplot() +
  # Lag-responses curves colored by temperature
  geom_line(
    data = matRRfit,
    aes(x = lag, y = RR, group = temperature, color = temperature),
    alpha = 0.35,
    linewidth = 0.35
  ) +
  scale_color_gradientn(
    colours = c(
      "#2166ac",
      "#4393c3",
      "#92c5de",
      "#d1e5f0",
      "#f7f7f7",
      "#fddbc7",
      "#f4a582",
      "#d6604d",
      "#b2182b"
    ),
    name = "Temperature"
  ) +
  # Start a new color scale for highlighted curves
  ggnewscale::new_scale_color() +
  # RR = 1 reference
  geom_hline(
    yintercept = 1,
    linetype = "dashed",
    color = "grey30",
    linewidth = 0.5
  ) +
  scale_x_continuous(breaks = cpred$lag_at) +
  scale_y_continuous(trans = "log10", breaks = pretty_breaks(6)) +
  labs(
    title = "Lag-response curves by temperature",
    x = "Lag (days)",
    y = "Relative Risk (RR)"
  ) +
  theme_minimal() +
  theme(legend.position = "top", panel.grid.minor.x = element_blank()) +
  transition_states(
    temperature,
    transition_length = 1,
    state_length = 0
  ) +
  shadow_mark(past = TRUE, future = FALSE, alpha = 0.6)

animate(p, nframes = 300, fps = 15, end_pause = 100)

allRRfit <- data.frame(
  temperature = as.numeric(rownames(cpred$allRRfit.summary)),
  lag = "overall",
  RR = cpred$allRRfit.summary[, "0.5quant"],
  RR_lci = cpred$allRRfit.summary[, "0.025quant"],
  RR_uci = cpred$allRRfit.summary[, "0.975quant"]
)

RRfit <- rbind(matRRfit, allRRfit)

# Split data
RRfit_lags <- RRfit |>
  filter(!lag %in% c("overall")) |>
  mutate(lag = as.numeric(lag))
RRfit_overall <- RRfit |>
  filter(lag %in% c("overall"))

temps <- cpred$exp_at
t_cold <- temps[which.min(abs(temps - quantile(temps, 0.01)))]
t_hot <- temps[which.min(abs(temps - quantile(temps, 0.99)))]

# Top plot: exposure-response curves for each lag and overall
p_main <- ggplot() +
  # Background: all lags, fading from vivid (small) to pale (large)
  geom_line(
    data = RRfit_lags,
    aes(x = temperature, y = RR, group = lag, color = lag),
    linewidth = 0.8
  ) +
  scale_color_gradientn(
    colours = c(
      "black",
      "#2b1d2f",
      "#4a2f5e",
      "#6a4c93",
      "#8b6bb8",
      "#b39cdb",
      "#d8c9f1",
      "#f3eef5"
    ),
    values = scales::rescale(c(0, 0.5, 1, 2, 3, 4, 5, 10, 20))
  ) +
  new_scale_color() +
  new_scale_fill() +
  # Credible intervals
  geom_ribbon(
    data = RRfit_overall,
    aes(
      x = temperature,
      ymin = RR_lci,
      ymax = RR_uci,
      fill = "1"
    ),
    alpha = 0.2
  ) +
  # Highlighted curves
  geom_line(
    data = RRfit_overall,
    aes(x = temperature, y = RR, color = "1"),
    linewidth = 1.2
  ) +
  geom_hline(
    yintercept = 1,
    linetype = "dashed"
  ) +
  scale_color_manual(values = "#a6761d", labels = "Overall (CrI95%)") +
  scale_fill_manual(values = "#a6761d", labels = "Overall (CrI95%)") +
  scale_y_continuous(
    transform = "log10",
    breaks = sort(c(0.8, pretty_breaks(5)(c(0.8, 4))))
  ) +
  labs(
    x = NULL,
    y = "Relative Risk (RR)",
    color = NULL,
    fill = NULL
  ) +
  theme_minimal() +
  theme(
    legend.position = "top",
    axis.text.x = element_blank(),
    plot.margin = margin(8, 8, 0, 8)
  )

# Bottom plot: histogram with percentile lines
p_hist <- ggplot(london, aes(x = tmean)) +
  geom_histogram(
    aes(y = after_stat(density), fill = after_stat(x)),
    binwidth = 0.5,
    color = "black",
    linewidth = 0.2
  ) +
  geom_vline(
    xintercept = t_cold,
    linetype = "dashed",
    color = "#053061",
    linewidth = 0.6
  ) +
  geom_vline(
    xintercept = t_hot,
    linetype = "dashed",
    color = "#67001f",
    linewidth = 0.6
  ) +
  geom_vline(
    xintercept = cen,
    linetype = "dashed",
    color = "grey20",
    linewidth = 0.6
  ) +
  annotate(
    "text",
    x = t_cold,
    y = Inf,
    label = "1st pct",
    vjust = 1.4,
    hjust = 1.1,
    size = 3.2,
    color = "#053061"
  ) +
  annotate(
    "text",
    x = t_hot,
    y = Inf,
    label = "99th pct",
    vjust = 1.4,
    hjust = -0.1,
    size = 3.2,
    color = "#67001f"
  ) +
  annotate(
    "text",
    x = cen,
    y = Inf,
    label = "MMT",
    vjust = 1.4,
    hjust = -0.1,
    size = 3.2,
    color = "grey20"
  ) +
  scale_x_continuous(limits = range(cpred$exp_at)) +
  scale_fill_gradientn(
    colours = c(
      "#053061",
      "#2166ac",
      "#4393c3",
      "#92c5de",
      "#d1e5f0",
      "#f7f7f7",
      "#fddbc7",
      "#f4a582",
      "#d6604d",
      "#b2182b",
      "#67001f"
    ),
    name = "Temperature"
  ) +
  labs(x = "Temperature (°C)", y = "Density") +
  theme_minimal() +
  theme(
    plot.margin = margin(20, 8, 8, 8),
    legend.position = "bottom"
  )

# Combine them:
p_main / p_hist + plot_layout(heights = c(3, 1))

tictoc::tic()
attr_forw <- attributable(
  mod,
  london,
  name_date = "date",
  name_exposure = "tmean",
  name_cases = "mort_75plus",
  cen = cen,
  dir = "forw"
)
tictoc::toc()

110.12 sec elapsed

col_af <- "black"

temp_colours <- c(
  "#053061",
  "#2166ac",
  "#4393c3",
  "#92c5de",
  "#d1e5f0",
  "#f7f7f7",
  "#fddbc7",
  "#f4a582",
  "#d6604d",
  "#b2182b",
  "#67001f"
)

af_med <- attr_forw$af.summary[, "0.5quant"]

# Pre-compute range once
af_min <- min(af_med, na.rm = TRUE) - 0.01
af_max <- max(af_med, na.rm = TRUE) + 0.01

df <- data.frame(
  date = london$date,
  x = yday(london$date),
  year = year(london$date),
  tmean = london$tmean,
  af = af_med
)

ggplot(df, aes(x = x)) +
  # Full-height temperature background per day
  geom_rect(
    aes(
      xmin = x - 0.5,
      xmax = x + 0.5,
      ymin = af_min,
      ymax = af_max,
      fill = tmean
    )
  ) +
  scale_fill_gradientn(
    colours = temp_colours,
    name = "Temperature (°C)"
  ) +
  # AF line on top
  geom_line(
    aes(y = af),
    color = col_af,
    linewidth = 0.7
  ) +
  scale_y_continuous(
    name = "Attributable Fraction (AF)",
    breaks = seq(0, 1, by = 0.1),
    limits = c(af_min, af_max),
    expand = c(0, 0)
  ) +
  scale_x_continuous(
    breaks = yday(as.Date(paste0(
      "2000-",
      c("01", "03", "05", "07", "09", "11"),
      "-01"
    ))),
    labels = c("Jan", "Mar", "May", "Jul", "Sep", "Nov"),
    expand = c(0, 0)
  ) +
  facet_wrap(~year, ncol = 3, axes = "all_x") +
  labs(x = NULL) +
  theme_minimal(base_size = 11) +
  theme(
    panel.spacing.x = unit(0, "pt"),
    strip.text = element_text(face = "bold", size = 10),
    legend.position = "top",
    legend.key.width = unit(2.5, "cm")
  )

rbind(
  "Attributable fraction" = attr_forw$aftotal.summary,
  "Attributable number" = attr_forw$antotal.summary
) |>
  as.data.frame() |>
  round(3) |>
  gt(rownames_to_stub = TRUE)

As we approach the decennial (10-year) anniversary since Jumping Rivers was founded in 2016, it’s a good time to reflect on what we have achieved in that time and share some lessons learned.

If you have read our blogs previously then you will be aware that Jumping Rivers is a consultancy and training provider in all things data science. But did you know that we offer over 50 different courses spanning R, Python, Git, SQL and more?

In this blog we will provide a glimpse into our internal process and share how we have streamlined the task of maintaining so many courses. Along the way we will share some good practices applicable to any big coding project, including packaging of source code and automated CI/CD.

The challenge

Let’s start by laying out the key challenges which face us.

1. Multilingual support

Our course catalogue consists of over 50 courses. The majority of these are either based on R or Python or both:

• 50% R
• 30% Python
• 5% R and Python
• 15% other (Git, SQL, Tableau, Posit and more)

At the very least, any solution that we come up with for standardising our courses must be compatible with both R and Python. Ideally it should also support some less taught languages including SQL and Git.

2. Maintenance

The world of R and Python is constantly changing. The languages themselves receive frequent updates, as do publicly available R packages on CRAN and Python packages on PyPI.

This has the consequence that code which worked one year ago (or even one day) may no longer be functional with the latest package versions. We will need some way to track this and ensure that the code examples covered in our courses remain relevant and error-free.

3. Demand

We deliver over 100 courses per year. For a relatively small team of data scientists, this can be a lot to juggle!

In an ideal world, the process of building the course materials, setting up the cloud environment for training, and managing all of the administration that goes along with this should be automated. That way, the trainer can focus on providing the highest quality experience for the attendees without having to worry about things going wrong on the day.

The solution

Our team is used to setting up data science workflows for clients, including automated reporting and migration of source code into packages. We have therefore applied these techniques in our internal processes, including training.

Automated reporting

You write a document which has to be updated on a regular basis; this might include a monthly presentation showing the latest company revenues. Does this scenario sound familiar?

We could regenerate the plots and data tables and manually copy and paste these into the report document. Even better, we can take advantage of free-to-use automated reporting frameworks including R Markdown and Quarto.

R Markdown and Quarto both work as follows:

• We provide a “YAML header” at the top of the report document with configuration and formatting options:

  ---
  title: "Introduction to Python"
  authors:
  - "Myles Mitchell"
  date: "2026-04-02"
  output: pdf
  ---
  

• The report body is formatted as Markdown and supports a mixture of plain text and code:

  ## Introduction
  At it's most basic, Python is essentially a calculator.
  We can run basic calculations as follows:
  ```{python}
  2 + 1
  ```
  We can also assign the output of a calculation to a
  variable so that it can be reused later:
  ```{python}
  x = 2 + 1
  print(x)
  ```
  

Notice that we have included chunks of Python code. By making use of chunk options we can configure code chunks to be executed when rendering the report. Any outputs from the code (plots, tables, summary statistics) can then be displayed.

By migrating the code logic into the report itself, we can update our report assets at the click of a button whenever the data changes.

We have taken inspiration from this approach with our course notes and presentation slides. This forces us to be rigorous with the code examples. Any runtime errors that are produced by faulty or outdated code would be visible in the course notes and by extension to the attendees of our courses.

Crucially for us, R Markdown and Quarto are both compatible with R and Python. They also support syntax highlighting for languages like Git and SQL, as well as a variety of output formats including HTML and PDF.

Internal R packages

So we have settled on a solution for building our course notes. But we have 50 different courses, and setting these up from scratch each time is going to get tedious!

A good practice in any coding project is to avoid duplication as much as possible. Instead of copying and pasting code, we should really be migrating code into functions which are self contained, reusable and easy to test. This will mean fewer places to debug when things inevitably go wrong.

Following a similar philosophy for our training infrastructure, we have migrated any reusable assets for our courses—including logos, template files and styling—into a collection of internal R packages.

When building a new course, the developer can now focus on the aspects that are unique to that course:

• Code examples
• Notes
• Exercises
• Presentation slides

Everything else is taken care of automatically:

• The appearance of the course notes and presentation slides.
• Build routines including converting the R Markdown / Quarto text files into HTML.

In addition to course templates, we also have internal packages for managing the administrative side of training, including:

• Calculating pricing quotes for clients.
• Generating post-course certificates.
• Spinning up a bespoke Posit Workbench environment for the course.
• Summarising attendee feedback.

And the list goes on!

GitLab CI/CD

With automated reporting and packaging of source code, we have created standardised routines that can be applied to any of our courses.

This does not change the fact that we have over 50 courses to maintain. We still need a way of testing our courses and tracking issues. This is where CI/CD (Continuous Integration / Continuous Development and Deployment) comes in.

CI/CD defines a framework for software development, including:

• Automated unit testing.
• Branching of source code and code review.
• Versioning and deployment of software.

If you maintain software then you have likely come across version control with Git. Cloud platforms like GitLab and GitHub provide tools for collaborative code development. Not only do they provide a cloud backup of your source code, they also provide the following features:

• CI/CD tools for automated testing, build and deployment.
• Branch rules for enforcing good practices like code review and unit testing.
• Versioning and tagging of source code.

Each of our courses is maintained via it’s own GitLab repository. The CI/CD pipelines for our courses are defined in a separate repository along with the internal R packages mentioned above.

When setting up a new course, the course repository will be automatically populated with the template CI/CD rules. All courses are therefore subject to the same stringent checks, including:

• Ensuring that the course notes build without errors.
• Enforcing code review of any course updates before these are merged into the main branch.
• Building and storing the artifacts (the rendered HTML notes and coding scripts) for the latest version of the course.

These checks are triggered by any updates to a course. We also schedule monthly CI/CD pipelines for all courses, with any issues immediately flagged to our trainers.

We have also taken advantage of GitLab’s folder-like structure for organising code repositories. Within the Jumping Rivers project on GitLab, we have a subproject called “training”. All of our course-related repositories are located “downstream” from this project. This means that any settings or environment variables defined at the “training” level are automatically applied to all of our courses.

In summary

The take-home lessons from this blog are applicable to any big coding project:

• Avoid duplication: migrate any reusable logic or assets into standalone packages.
• Utilise CI/CD workflows using GitLab, GitHub or similar.
• Focus on what matters by automating as much of the process as possible.

Our training infrastructure has taken 10 years to build and is still constantly evolving; we have not even covered the full process in this blog! For a deeper dive, check out this talk by Myles at SatRdays London 2024.

For more on automated reporting, check out:

• Quarto for the Python user.
• Parameterised presentations with Quarto.

For more on packaging of source code, check out:

• Writing a personal R package.
• Three-part series: Creating a Python package.
• Four-part series: R package quality.

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

I’ve posted about publication lag times previously. The “lag” refers to the time from submitting a paper and it appearing in a journal.

Publication lag times are still a frustration for researchers. Although preprints circumvent the delay in sharing science with others, publication is still king when it comes to evaluation. Contracts are short and publication delays can be long…

I recently saw a post comparing median publication lag times for genetics journals. This motivated me to update my code and rerun the analysis for cell biology journals to see what, if anything, has changed.

Methodology

I wrote an R package PubMedLagR which uses {rentrez} to retrieve the publication data from PubMed. Once we have this data, it is a matter of producing some plots which I have covered previously. To ensure that we are only looking at research papers, we use "journal article"[pt] as a search term, and also remove from the data anything else (Reviews, Commentaries etc.). Feel free to use it to look at other journals.

Caveats

Before we get started, there are some caveats. The analysis is only as good as the PubMed data. Not all journals submit their date information to PubMed, while for others it is incomplete (as we’ll see below). There are inaccuracies for sure. I found a paper that was supposedly submitted on 1970-01-01, more than 40 years before the journal started. Also, it’s well known that some journals “restart the clock” on a paper by rejecting it and allowing resubmission, then only counting the resubmitted version. So, comparison between journals is a little tricky, but it allows us to look at trends.

Let’s dive in

We can use the following code to grab the data we are interested in.

library(PubMedLagR)
jrnl_list <- c("J Cell Sci", "Mol Biol Cell", "J Cell Biol", "Nat Cell Biol",
               "EMBO J", "Biochem J", "Dev Cell", "FASEB J", "J Biol Chem",
               "Cells", "Front Cell Dev Biol", "Nature Communications",
               "Cell Reports", "Mol Cell", "Autophagy", "Cell Death Differ",
               "Cell Death Dis", "Cell Res", "Sci Adv", "Cell")
yrs <- 2006:2026
retrieve_journal_year_records(jrnl_list, yrs, batch_size = 250)
pprs <- pubmed_xmls_to_df()

The list of journals is somewhat arbitrary. I have included Nature Communications and Science Advances although they carry many other papers besides cell biology. I included Cell (even though there’s not much cell biology in there these days) and left out Nature and Science.

How many papers are in each of these journals and how has that changed over time?

There’s a huge increase in the number of papers published by Nature Communications, Science Advances and Cell Reports. Nature Communications is a behemoth, publishing ~12,500 papers in 2025.

There was a boom and bust in publications at Front Cell Dev Biol and Cells which could be due to reputational problems (like this). Other journals have declined. Most noticeably J Biol Chem, but others have taken a hit. The reasons behind these dynamics are discussed elsewhere.

Median publication lag time

We’ll use received-to-published as our measure of publication lag time. This is the time from submission to it appearing in the journal “in print”. It’s measured here in days.

For the last 5 years, the median lag time at Nature Cell Biology is over one year. Other outlets are close to this, e.g. Cell, Dev Cell; whereas others linger around 200 days, or lower in the case of J Cell Sci, FASEB J et al. The journal Autophagy is missing here because there are no data for it. Others, like Sci Adv, MBoC have only minimal data available. The shortest lag times were for Cells (which might not surprise some people).

The lion’s share of this lag time is the time from submission to acceptance (received-accepted), with a small contribution from the time taken to formally publish the article (accepted-published).

This is using the median time. Obviously, some papers whizz straight in, whereas others… don’t. Let’s have a look.

You can see many examples of papers that have lag times of 1, 2 or 3 years. I scaled all the plots to a maximum lag time of 3 years. There were several examples of papers with lag times of up to 7 years that looked to be genuine, but they distorted the view.

So what trends can we see?

The lag times are creeping up at some journals, but not at others. It’s possible to see some very short lag times (in amongst the longer ones) recently at Dev Cell, Mol Cell and EMBO J. I assume these are transfers in to the journal, leading to rapid publication. Although Cell also has many of those too, so perhaps there are other explanations.

To click though in more detail, I made some graphics for each journal which use ridgelines to look at the profile of lag times for each year.

And finally

Incidentally, nine of the top ten papers with the longest lag times were published in Nature Communications. The longest was this one (3263 days). Almost nine years! All papers have their battle stories and I’m sure this one has a tale to tell.

—

The post title comes from “Hold On Hope” by Guided by Voices.

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

This article is cross-posted on rOpenSci and R-Consortium blogs.

For more than two decades, the Bioconductor project has been a cornerstone of the R ecosystem, providing high-quality, peer-reviewed tools for bioinformatics and computational biology. Its curated repository model, rigorous review standards, and tightly coordinated release process have helped establish Bioconductor as one of the most trusted distribution channels in scientific computing.

However, the infrastructure that supports such a long-standing and large-scale project inevitably accumulates technical debt. Legacy build systems, bespoke tooling, and historically grown workflows add up to costly and unsustainable maintenance work. For this reason, Bioconductor is collaborating with R-universe to gradually modernize parts of its infrastructure, while accommodating the project’s scale, governance, and established processes. In turn, Bioconductor is helping R-universe expand and refine its features as we learn to serve the complex needs of the Bioconductor community.

This collaboration reflects a core principle of R-universe as an R Consortium Infrastructure Steering Committee (ISC) top-level project: supporting reviewed package repositories such as rOpenSci and Bioconductor, and providing modern, open, and reusable infrastructure that strengthens the broader R ecosystem.

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

For more than two decades, the Bioconductor project has been a cornerstone of the R ecosystem, providing high-quality, peer-reviewed tools for bioinformatics and computational biology. Its curated repository model, rigorous review standards, and tightly coordinated release process have helped establish Bioconductor as one of the most trusted distribution channels in scientific computing.

However, the infrastructure that supports such a long-standing and large-scale project inevitably accumulates technical debt. Legacy build systems, bespoke tooling, and historically grown workflows add up to costly and unsustainable maintenance work. For this reason, Bioconductor is collaborating with R-universe to gradually modernize parts of its infrastructure, while accommodating the project’s scale, governance, and established processes. In turn, Bioconductor is helping R-universe expand and refine its features as we learn to serve the complex needs of the Bioconductor community.

This collaboration reflects a core principle of R-universe as an R Consortium Infrastructure Steering Committee (ISC) top-level project: supporting reviewed package repositories such as rOpenSci and Bioconductor, and providing modern, open, and reusable infrastructure that strengthens the broader R ecosystem.

A shared mission: Tooling for managed repositories

R-universe was designed as a next-generation package distribution and build system for R. It provides:

• Continuous building and checking of R packages across platforms
• Binary packages for Windows, macOS, Linux, and WebAssembly
• Transparent and reproducible build environments managed via GitHub actions
• Dashboards and metadata APIs for monitoring ecosystem health and activity
• CRAN-like package repositories with discoverable metrics and documentation

From the outset, a key objective has been to support curated and reviewed communities — such as rOpenSci and Bioconductor — by offering modern infrastructure without requiring them to redesign their governance model or review processes.

For Bioconductor, this means incrementally introducing piece-wise functionality, with consideration for established release cycles and quality control mechanisms:

1. Setting up independent build and dashboard tooling, replicating processes from the current Bioconductor build systems on R-universe infrastructure
2. Mirroring Windows and macOS binaries produced on R-universe to Bioconductor
3. Exploring further integration of results and metadata produced by R-universe for Bioconductor health/activity monitoring and aiding the curation processes
4. Potential future steps toward deeper automation and harmonization

By taking small gradual steps towards adopting R-universe components, everyone gets the opportunity to experiment with new tooling and evaluate where adjustments may be needed in order to minimize disruption to existing practices.

An important milestone in this venture is that Bioconductor now uses R-universe to build the Windows and macOS binaries, which significantly reduces costs and the maintenance load on the Bioconductor team. Beyond binary distribution, we are currently exploring deeper integration of R-universe’s continuous check results into Bioconductor’s quality control and release processes.

Two Universes: Release and Development

Bioconductor maintains two distinct repositories:

• A release branch for stable packages
• A devel branch for ongoing development and the next release cycle

To mirror this structure, we currently operate two dedicated R-universe instances:

• Development branch: https://bioc.r-universe.dev
• Release branch: https://bioc-release.r-universe.dev

These universes integrate directly with Bioconductor’s existing Git infrastructure and provide continuous builds for packages in both branches.

Through the R-universe dashboard, package maintainers and users can:

• Inspect cross-platform check results
• Review extended BiocCheck diagnostics
• Monitor build logs and dependency graphs
• Explore rich package metadata and metrics
• Publish binary packages for Windows, macOS, and Linux

This provides a familiar yet modern interface for Bioconductor contributors, aligned with what users increasingly expect from contemporary R package infrastructure.

Information about each package is available on https://bioc.r-universe.dev/{pkgname}. For example, https://bioc.r-universe.dev/DESeq2 provides details on the DESeq2 package as shown below:

If this is your first time visiting R-universe, we recommend clicking the “Website Tour” button which will walk you through the most important information in 1 or 2 minutes.

Technical Documentation for Bioconductor Maintainers

The R-universe project maintains comprehensive technical documentation at https://docs.r-universe.dev. For Bioconductor specifically, we created a dedicated section summarizing the most relevant topics for developers to get started with R-universe: https://docs.r-universe.dev/bioconductor/

As the collaboration evolves and new components get introduced, the documentation will continue to be expanded. The goal is to provide Bioconductor maintainers with a clear reference point for understanding how R-universe fits into their development workflow, while maintaining compatibility with the established practices that have made Bioconductor a successful project within the R community.

Looking Ahead

Adopting new infrastructure inevitably involves adjustments. For Bioconductor developers, integrating with a new build and distribution system will likely require some changes to workflows, and time to become familiar with new or different package checks, build diagnostics, and binary distribution.

However, by gradually moving toward common infrastructure, the Bioconductor project will benefit from improvements that are being continuously developed and maintained for the broader R ecosystem. A system based on modern continuous integration (CI) will provide developers with improved tooling, and will give the core team more time to focus on community coordination and quality control, rather than on maintaining costly infrastructure. At the same time, the shared platform provided by R-universe can help to increase the visibility and accessibility of Bioconductor software to the greater R community.

We look forward to continuing this alliance and to working with the Bioconductor community to ensure that the next generation of infrastructure supports the project for many years to come.

I do not usually write posts that are calls to action. But sometimes, something important enough comes along that it would feel wrong to stay silent. This is one of those times.

How does the average marathoner pace their race? In this post, we’ll use R to have a look at a large dataset of marathon times to try to answer this question.

The ideal strategy would be to “even split” the race. This is where you run continually at the same pace from kilometre 0 to the finish. Let’s forget about “negative splitting”. This is where you speed up through the race, usually by running at a constant pace for the first half or three-quarters and then increasing the pace. Negative splits are for the pros not mere mortals! The difficulty with even-splitting the race is that it is very hard to know what pace you can maintain. The marathon gets hard for everyone after 30 km, so a slow down is almost inevitable. Certainly if you have started too fast you will fade. This situation is known as “positive splitting”.

Why is it so hard to know what pace you can maintain? Well, you can predict a pace based on existing races e.g. half marathon, and there are various ways to do this, but it is difficult to tell if you can hold that pace for the marathon. It’s such a brutal event that training up to run one takes time and it equally takes a while to recover, so experimentation is limited. Running a full marathon (at pace) in training, is not advised. So determining an ideal pace involves quite a bit of guesswork.

Let’s take a look at a big dataset of marathon times – we’ll use the New York City Marathon from 2025 – to see if we can understand how to pace a marathon. There’s an available dataset of chip times (meaning we don’t have to worry about dodgy GPS data) and the course has similar first and second half profiles, allowing us to use these times to understand negative/even/positive splitting. Let’s dive in.

You can skip to the code to play along or just see the analysis here.

First we can see using histograms of the difference between second half and first half of the marathon, that most runners positive split the marathon. There are very few runners who run a negative-split (blue bars, left of the dashed line). More runners even-split (yellow), but the majority run positive (red) split times.

For marathoners with finishing in times of below 3 h, the modal split is only +2 minutes. Over 21.1 km this is only a loss of 6 s per km. For marathoners with finishes of over three hours, this loss gets more severe. Those finishing outside of 5 h, ship 20 minutes or more in the second half.

At first glance this looks like better pace management by the faster runners, but these positive splits could be proportional to the paces being run. In other words, a slower runner should ship more time in the second half, because they’re running more slowly.

We can look at this data a different way and directly compare the first and second half times for each runner. Again this highlights just how few runners negative- or even-split the marathon. Most are positive splitting and are in the upper left half of the plot. We can also see that the data veers away from the ideal even-split (dashed line) with the slower paces. This veering looks linear (straight line).

We can fit a line to this data, and constrain it to go through (1,1) i.e. a 2 h marathoner even-splitting the race. To do this in R we can use lm(formula = I(y - 60) ~ I(x - 60) + 0, data = fitting) and this gives the coefficient for I(x – 60) as 1.24. This is essentially the fade co-efficient for the average runner in the 2025 edition of this race.

What does that mean? Well, for a runner achieving a 90 minute first half, their second half would most likely be: 60 + 1.239 * (90 – 60) = 97.17 minutes, so this would be a finish time of 3:07:10.

For anyone looking to run a 3 h New York Marathon, the average runner would therefore need to run 60 / 2.239 + 60 = 86.8 minutes for the first half to anticipate the fade. So 1:26:48 for the first half, and then 1:33:12 for the second half.

A more simple calculation is to take the mean of the ratio between the two half times for everyone in the dataset. This gives a fade coefficient of 1.13. The difference between these two fade co-efficients is due to the lack of constraint used in the fit. The ratio predicts a positive split being inevitable for the fastest runners, which is probably not true. Anyhow, this puts the first half time at 88 minutes for folks looking to run 3 h. These fade co-efficients are good predictors for a range of times, and I suspect would be similar at other marathon events with a similar profile. You can use them to calculate your ideal pace for a target finish time.

Finally, for the most accurate answer about sub-3 h pacing, we can look directly at runners finishing between 02:50:00 and 03:00:00 and see what they actually ran. The median first half time was 86.3 min (IQR = 84.4 – 87.87) and the second half was 89.62 (88.07 – 91.12). This gives a median finish time of 2:56:00. So running a 1:26:18 first half would give someone their best chance of finishing in under 3 h, allowing for the inevitable fade.

The takeaway message is: to finish within a goal time, do not assume even splits. That is, if you want to run 3 hours 30 min and bank on 90 minutes per half (4:59/km), you will most likely fail to hit the target. Build in a buffer of time to allow for the inevitable fade. A pace of 4:45/km is a better target pace (see below).

Good luck!

The code

This analysis was possible thanks to the uploader for making the chip time data available. Also, a shoutout to Nicola Rennie for sharing how to style social media handles in {ggplot2} graphics. This part of my code requires my {qBrand} library and should be skipped if you are running the code yourself (remove the caption = cap argument in the ggplot calls).

library(ggplot2)
library(ggtext)

sysfonts::font_add_google("Roboto", "roboto")
showtext::showtext_auto()

## data wrangling ----

# load csv file from url
url <- paste0("https://huggingface.co/datasets/donaldye8812/",
              "nyc-2025-marathon-splits/resolve/main/",
              "nyrr_marathon_2025_summary_56480_runners_WITH_SPLITS.csv")
df <- read.csv(url)

# the data frame is a long table
# we need to grab the time values where splitCode is "HALF" or "MAR"
df <- df[df$splitCode %in% c("HALF", "MAR"), c("RunnerID", "splitCode", "time")]
# reshape to wide format, values are in time
df <- reshape(df, idvar = "RunnerID", timevar = "splitCode", direction = "wide")
# calculate the split times in minutes
df$split_HALF <- as.numeric(
  as.difftime(df$time.HALF, format = "%H:%M:%S", units = "mins"))
df$split_MAR <- as.numeric(
  as.difftime(df$time.MAR, format = "%H:%M:%S", units = "mins"))
# calculate the second half time
df$split_SECOND_HALF <- df$split_MAR - df$split_HALF
# remove rows with NA values
df <- df[!is.na(df$split_SECOND_HALF), ]
# calculate the difference
df$Difference <- df$split_SECOND_HALF - df$split_HALF
# difference as a fraction of first half
df$Difference_Fraction <- df$Difference / df$split_HALF * 100
# classify into sub 3 hr, sub 4 hr, sub 5 hr, sub 6 hr, over 6 hr
df$Category <- cut(df$split_MAR,
                           breaks = c(0, 180, 210, 240, 300, Inf),
                           labels = c("Sub 3 h", "3:00-3:30", "3:30-4:00",
                                      "4:00-5:00", "Over 5 h"))

## plot styling ----

social <- qBrand::qSocial()
cap <-  paste0(
  "**Data:** New York City Marathon 2025 Results<br>**Graphic:** ",social
)

my_palette <- c("Sub 3 h" = "#cb2029",
                "3:00-3:30" = "#147f77",
                "3:30-4:00" = "#cf6d21",
                "4:00-5:00" = "#28a91b",
                "Over 5 h" = "#a31a6d")

## make the plots ----

ggplot(df, aes(x = Difference, fill = after_stat(x))) +
  # vertical line at x = 0
  geom_vline(xintercept = 0, linetype = "dashed", color = "black") +
  geom_histogram(breaks = seq(
    from = -59.5, to = 81.5, by = 1), color = "black") +
  scale_colour_gradient2(
    low = "#2b83ba",
    mid = "#ffffbf",
    high = "#d7191c",
    midpoint = 0,
    limits = c(-15,15),
    na.value = "#ffffffff",
    guide = "colourbar",
    aesthetics = "fill",
    oob = scales::squish
  ) +
  scale_x_continuous(breaks = seq(-45,90,15), limits = c(-40, 80)) +
  facet_wrap(~ Category, ncol = 1, scales = "free_y") +
  labs(caption = cap) +
  labs(title = "Most runners positive split the marathon",
       x = "Difference in minutes (Second Half - First Half)",
       y = "Number of Runners",
       caption = cap) +
  theme_classic() +
  # hide legend
  theme(legend.position = "none") +
  theme(
    plot.caption = element_textbox_simple(
      colour = "grey25",
      hjust = 0,
      halign = 0,
      margin = margin(b = 0, t = 5),
      size = rel(0.9)
    ),
    text = element_text(family = "roboto", size = 16),
    plot.title = element_text(size = rel(1.2),
                              face = "bold")
  )

ggsave("Output/Plots/nyc_marathon_2025_split_difference_histogram.png",
       width = 900, height = 1200, dpi = 72, units = "px", bg = "white")

ggplot() +
  geom_abline(slope = 1, linetype = "dashed", color = "black") +
  geom_point(data = df,
             aes(x = split_HALF, y = split_SECOND_HALF, colour = Category),
             shape = 16, size = 1.5, alpha = 0.1) +
  scale_x_continuous(breaks = seq(from = 0, to = 12 * 30, by = 30),
                     labels = seq(from = 0, to = 6, by = 0.5),
                     limits = c(1 * 60, 5 * 60)) +
  scale_y_continuous(breaks = seq(from = 0, to = 12 * 30, by = 30),
                     labels = seq(from = 0, to = 6, by = 0.5),
                     limits = c(1 * 60, 5 * 60)) +
  scale_colour_manual(values = my_palette) +
  labs(x = "First half time (h)",
       y = "Second half time (h)",
       caption = cap) +
  theme_bw() +
  theme(
    plot.caption = element_textbox_simple(
      colour = "grey25",
      hjust = 0,
      halign = 0,
      margin = margin(b = 0, t = 10),
      size = rel(0.9)
    ),
    text = element_text(family = "roboto", size = 16)
  ) +
  guides(colour = guide_legend(override.aes = list(alpha = 1)))

ggsave("Output/Plots/nyc_marathon_2025_split_difference_scatter.png",
       width = 1000, height = 800, dpi = 72, units = "px", bg = "white")

From this data we can also make some calculations to understand…

## fitting ----

# to fit, we'll constrain the line to go through (60,60), i.e. a
# 2 h marathoner who runs even splits
fitting <- data.frame(x = df$split_HALF,y = df$split_SECOND_HALF)
lm( I(y-60) ~ I(x-60) + 0, data = fitting)


# Call:
#   lm(formula = I(y - 60) ~ I(x - 60) + 0, data = fitting)
# 
# Coefficients:
#   I(x - 60)  
# 1.239  

# so for a 90 minute first half, second half would be:
# 60 + 1.239 * (90 - 60) = 97.17 minutes, a finish time of 3:07:10

# to run a 3 h New York Marathon, the average runner needs to run
# 60 / 2.239 + 60 = 86.8 minutes for the first half
# so 1:26:48 for the first half, and 1:33:12 for the second half

# a more simple approach is to calculate the mean of the ratios
mean_ratio <- mean(df$split_SECOND_HALF / df$split_HALF)
mean_ratio
# [1] 1.127581

# filter the df for finish times between 170 and 180 minutes
target <- df[df$split_MAR > 170 & df$split_MAR < 180,]
summary(target)


    RunnerID         time.HALF           time.MAR           split_HALF      split_MAR     split_SECOND_HALF   Difference     
 Min.   :48819892   Length:1289        Length:1289        Min.   :70.25   Min.   :170.0   Min.   : 82.70    Min.   :-7.6833  
 1st Qu.:48834548   Class :character   Class :character   1st Qu.:84.42   1st Qu.:173.6   1st Qu.: 88.07    1st Qu.: 0.7167  
 Median :48849752   Mode  :character   Mode  :character   Median :86.30   Median :176.0   Median : 89.62    Median : 3.0500  
 Mean   :48849498                                         Mean   :85.98   Mean   :175.7   Mean   : 89.73    Mean   : 3.7585  
 3rd Qu.:48864551                                         3rd Qu.:87.87   3rd Qu.:178.2   3rd Qu.: 91.12    3rd Qu.: 5.9000  
 Max.   :48878979                                         Max.   :92.87   Max.   :180.0   Max.   :106.02    Max.   :35.7667  
 Difference_Fraction      Category   
 Min.   :-8.5008     Sub 3 h  :1289  
 1st Qu.: 0.8051     3:00-3:30:   0  
 Median : 3.5390     3:30-4:00:   0  
 Mean   : 4.5178     4:00-5:00:   0  
 3rd Qu.: 7.0055     Over 5 h :   0  
 Max.   :50.9134   

—

The post title comes from “Marathon Man” by Ian Brown from his “My Way” album. He’s wearing a track suit on the cover but that’s not optimal wear for running a marathon.

A new version of unifiedml is out; available on CRAN. unifiedml is an effort to offer a unified interface to R’s machine learning models.

The main change in this version 0.2.1 is the removal of type (of prediction) from predict, and the use of... instead, which is more generic and flexible.

This post contains advanced examples of use of unifiedml for classification, with ranger and xgboost. More examples have been added to the package vignettes too.

install.packages("unifiedml")

install.packages(c("ranger"))

library("unifiedml")

Loading required package: doParallel

Loading required package: foreach

Loading required package: iterators

Loading required package: parallel

Loading required package: R6

1 – ranger example

library(ranger)



# 2 - 'ranger' classification ---------------------------

# -------------------------------
# S3 wrapper for ranger
# -------------------------------

# Fit function remains the same
my_ranger <- function(x, y, ...) {
  if (!is.data.frame(x)) x <- as.data.frame(x)
  y <- as.factor(y)
  colnames(x) <- paste0("X", seq_len(ncol(x)))
  df <- data.frame(y = y, x)
  fit <- ranger::ranger(y ~ ., data = df, probability = TRUE, ...)
  structure(list(fit = fit), class = "my_ranger")
}

# Predict only with newdata
predict.my_ranger <- function(object, newdata = NULL, newx = NULL, ...) {
  if (!is.null(newx)) newdata <- newx
  if (is.null(newdata)) stop("No data provided for prediction")
#  misc::debug_print(newx)
#  misc::debug_print(newdata)
  if (is.matrix(newdata)) newdata <- as.data.frame(newdata)
#  misc::debug_print(newdata)
  # Unconditionally rename to match training
  colnames(newdata) <- paste0("X", seq_len(ncol(newdata)))
#  misc::debug_print(newdata)
  preds <- predict(object$fit, data = newdata)$predictions
#  misc::debug_print(newdata)
  if (is.matrix(preds) && ncol(preds) == 2) {
    lvls <- colnames(preds)
    return(ifelse(preds[, 2] > 0.5, lvls[2], lvls[1]))
  }

  preds
}

# Print method
print.my_ranger <- function(x, ...) {
  cat("my_ranger model\n")
  print(x$fit)
}

# -------------------------------
# Example: Iris binary classification
# -------------------------------

set.seed(123)
iris_binary <- iris[iris$Species %in% c("setosa", "versicolor"), ]
X_binary <- iris_binary[, 1:4]
y_binary <- as.factor(as.character(iris_binary$Species))

# Train/test split
train_idx <- sample(seq_len(nrow(X_binary)), size = 0.7 * nrow(X_binary))
X_train <- X_binary[train_idx, ]
y_train <- y_binary[train_idx]
X_test <- X_binary[-train_idx, ]
y_test <- y_binary[-train_idx]

# Initialize and fit model
# Initialize model
mod <- Model$new(my_ranger)

# Fit on training data only
mod$fit(X_train, y_train, num.trees = 150L)

# Predict on test set
preds <- mod$predict(X_test)

# Evaluate
table(Predicted = preds, True =y_test)
mean(preds == y_test)  # Accuracy



# 5-fold cross-validation on training set
cv_scores <- cross_val_score(
  mod,
  X_train,
  y_train,
  num.trees = 150L,
  cv = 5L
)

cv_scores
mean(cv_scores)  # average CV accuracy


            True
Predicted    setosa versicolor
  setosa         15          0
  versicolor      0         15

1

  |======================================================================| 100%

1. 1
2. 1
3. 1
4. 1
5. 1

1

2 - xgboost example

library(xgboost)

my_xgboost <- function(x, y, ...) {
  
  # Convert to matrix safely
  if (!is.matrix(x)) {
    x <- as.matrix(x)
  }
  
  # Handle factors
  if (is.factor(y)) {
    y <- as.numeric(y) - 1
  }
  
  fit <- xgboost::xgboost(
    data = x,
    label = y,
    ...
  )
  
  structure(list(fit = fit), class = "my_xgboost")
}

predict.my_xgboost <- function(object, newdata, ...) {
  
  # Ensure matrix
  newdata <- as.matrix(newdata)
  
  preds <- predict(object$fit, newdata)
  
  # If binary classification → convert probs to class
  if (!is.null(object$fit$params$objective) &&
      grepl("binary", object$fit$params$objective)) {
    
    return(ifelse(preds > 0.5, 1, 0))
  }
  
  preds
}

predict.my_xgboost <- function(object, newdata = NULL, newx = NULL, ...) {
  
  # Accept both conventions
  if (!is.null(newx)) {
    newdata <- newx
  }
  
  newdata <- as.matrix(newdata)
  
  preds <- predict(object$fit, newdata)
  
  # Binary classification → class labels
  if (!is.null(object$fit$params$objective) &&
      grepl("binary", object$fit$params$objective)) {
    
    return(ifelse(preds > 0.5, 1, 0))
  }
  
  preds
}

print.my_xgboost <- function(x, ...) {
  cat("my_xgboost model\n")
  print(x$fit)
}


set.seed(123)  # for reproducibility

# Binary subset
iris_binary <- iris[iris$Species %in% c("setosa", "versicolor"), ]
X_binary <- as.matrix(iris_binary[, 1:4])
y_binary <- as.factor(as.character(iris_binary$Species))

# Split indices: 70% train, 30% test
train_idx <- sample(seq_len(nrow(X_binary)), size = 0.7 * nrow(X_binary))
X_train <- X_binary[train_idx, ]
y_train <- y_binary[train_idx]
X_test <- X_binary[-train_idx, ]
y_test <- y_binary[-train_idx]

# Initialize model
mod <- Model$new(my_xgboost)

# Fit on training data only
mod$fit(X_train, y_train, nrounds = 50, objective = "binary:logistic")

# Predict on test set
preds <- mod$predict(X_test)

# Evaluate
table(Predicted = preds, True =y_test)
mean(preds == y_test)  # Accuracy



# 5-fold cross-validation on training set
cv_scores <- cross_val_score(
  mod, 
  X_train, 
  y_train, 
  nrounds = 50, 
  objective = "binary:logistic", 
  cv = 5L
)

cv_scores
mean(cv_scores)  # average CV accuracy

Last summer, while relaxing on the beaches of Berck, a French town known for treating tuberculosis in kids by exposing them to the fresh maritime air (back in 19th century, they have antibiotics these days), I found myself daydreaming about building my own programming language.

Spoiler alert: I don’t know how to build programming languages, but I have developed extremely strong opinions over the years about the features a modern data science language should have. So could I use them fancy LLMs to build one?

Also, let’s get one question answered straight away: why create a new language instead of contributing to existing ones? I certainly do contribute, I maintain several R packages like {rix}, {rixpress}, and {chronicler}, and even have two Python packages (cronista and ryxpress), but I wanted a clean slate to build a system centered around a few non-negotiable principles and features I’ve implemented over the years in R:

• Reproducibility-First: A language where reproducibility isn’t a bolt-on afterthought managed by external tools, but the very foundation of the runtime.
• Aggressive Re-use: Instead of reinventing the wheel, this language would stand on the shoulders of giants. It’d use Nix for package management and environment isolation, and Apache Arrow as its high-performance backbone for data frames. R, Python, Julia and other languages would provide the algorithms and models.
• First-Class Pipelines: Scripts shouldn’t be a sequence of side-effects. In this language, pipelines would be mandatory and first-class citizens.
• Fail Early and Loudly: No silent type conversions or hidden NAs. If something is wrong, the language breaks immediately so you can fix it.
• Errors as Objects: Inspired by functional programming, errors are first-class values that can be inspected and handled gracefully.
• Two Pipes: I want two pipes, one for linear transformations, |>, and a maybe-pipe, ?|> for error recovery. Unlike the standard pipe, ?|> always forwards its value, including Errors, to the next function, allowing you to write handlers that inspect and potentially recover from them. Since Errors are just values, this composes naturally with the rest of the language.
• Polyglot by Design: Rather than re-implementing every statistical algorithm, this language would be designed to orchestrate and bridge R, Python, and Julia seamlessly.

Also, we’re in a post LLM world, and like them or not, they’re here to stay. They’re pretty useful to write boilerplate code and so any new language would be dead on arrival if it didn’t play nicely with LLMs. So such a new language would need to be written for LLMs primarily, because I don’t expect anyone to learn any new language. This is where the declarative nature of Nix is a huge advantage. Because environments are precisely described, it is much easier for LLMs to focus on generating code and not have to fight with environment setup. This is also the reason I took another radical decision: since Nix would be mandatory for setting up the environment, why bother building OS-specific binaries? I’d just build a Nix package for this language and let Nix handle the rest.

This architecture results in a DSL for orchestration, making it trivial to transfer data objects between different ecosystems without the usual FFI (Foreign Function Interface) friction.

With these ideas in mind, I started prompting Gemini to brainstorm and started by generating specification files. Very broad first, but as days went by, more and more focused. The way I went about it (and still go) is that I first brainstorm an idea with an LLM, then I ask it to generate a specification file, then I refine it, ask it to generate a new specification file, and so on. Once I’m happy with the spec, I ask an LLM to generate a minimal implementation of the spec. Usually writing the spec and a first implementation is a task shared between Claude and Gemini (through Antigravity). Then I open a pull request and ask GitHub Copilot to review it (usually with GPT-5.x). I repeat this process until I’m happy with the implementation. I always ask for documentation and unit tests (and golden tests when relevant, more on this later).

I started to really believe that I had something interesting, so I gave it a shot, and called it T. I had long joked that the natural successor to R should be called T (because R is the successor to S… and no, I’m not going to call it Q because that sounds like the word for ass in French).

Something else that made me confident I could succeed, besides my own hubris, was that I am pretty familiar with unit testing, test-driven development, trunk-based development and Nix. When you combine all these elements, it makes developing with LLMs quite safe.

So I just started prompting. And now I’m quite happy to announce that there is a beta version of T that you can use today!

By leveraging Nix as a build engine, T can treat complex data science workflows as buildable derivations. A typical T pipeline looks like this:

p = pipeline {
  -- 1. Python node: read data with pandas
  mtcars_pl = pyn(
    command = <{
import pandas as pd
pd.read_csv("data/mtcars.csv", sep="|")
    }>,
    include = ["data/mtcars.csv"],
    serializer = ^csv
  )

  -- 2. Python node: filter and serialize as CSV
  mtcars_pl_am = pyn(
    command = <{
mtcars_pl[mtcars_pl['am'] == 1]
    }>,
    deserializer = ^csv,
    serializer = ^csv
  )

  -- 3. R node: read CSV and take head using functions.R
  mtcars_head = rn(
    command = <{
my_head(mtcars_pl_am)
    }>,
    functions = ["src/functions.R"],
    deserializer = ^csv,
    serializer = ^csv
  )

  -- 4. R node: select column with dplyr
  mtcars_mpg = rn(
    command = <{
library(dplyr)
mtcars_head %>% select(mpg)
    }>,
    deserializer = ^csv,
    serializer = ^csv
  )

  -- Render Quarto report
  report = node(script = "src/report.qmd", runtime = Quarto)
}

-- Materialize the pipeline
populate_pipeline(p, build = true)
pipeline_copy() -- Copy the outputs from the Nix store to your working directory

As you can see, each node has a command argument where you can write literal R or Python code. It is also possible to provide the path to a script instead. If packages need to be loaded for the code to work, you can just write the calls to load the required packages in the command argument as well.

While T is heavily inspired by the {targets} package in R, it takes the concept a step further by making pipelines first-class objects within the language itself. This means you can:

• Compose Pipelines: You can define small, modular pipelines and then merge them into larger ones using standard operators.
• Static Analysis: Because the DAG (Directed Acyclic Graph) is defined within the language, T can validate your entire workflow (checking for circular dependencies or missing data) before a single line of code even runs.
• Heterogeneous Execution: A single pipeline can effortlessly mix R, Python, and native T code. Data is passed between these nodes using built-in serializers like ^csv, ^arrow, or even specialized formats like ^pmml for traditional models and ^onnx for deep learning architectures. It is also possible to define your own serializers.
• Immutable State: Each node output is managed by Nix, meaning if you haven’t changed the code or the data for a specific node, T (via Nix) will simply pull the cached result from previous runs.

But don’t let the “orchestrator” label fool you; T is also a capable language in its own right. It features a selection of built-in packages inspired by the tidyverse for data manipulation. Thanks to its Arrow backend, it is surprisingly fast. I even maintain a CI benchmark running on NYC Taxi data to ensure performance remains competitive.

I made sure that T is pretty easy to use with LLMs by providing a file called summary.md in the root of the GitHub repository. This file is meant to be used by LLMs to quickly learn the language’s syntax and generate code accordingly. You could also provide the whole help documentation to the LLM (found in the repository under help/docs.json), but I found that a summary is usually enough. There is also another experimental feature I’m thinking about, called intent blocks. These blocks would essentially be first-class structured comments that would be used to anchor LLM’s behaviour and make it more deterministic. These blocks would be parsed by T and used to generate code accordingly. I have some ideas how these could look like, something like this:

intent {
  description: "Customer churn prediction",
  assumptions: ["Age > 18", "NA imputed with multiple imputation"],
  requires: ["dataset.csv"]
}

# 1. Initialize a new project
nix run github:b-rodrigues/tlang -- init --project my_t_project
cd my_t_project

[project]
name = "r_py_xgboost_t"
description = "A T data analysis project"

[dependencies]
# T packages this project depends on
# Format: package = { git = "repository-url", tag = "version" }
# Example:
# stats = { git = "https://github.com/t-lang/stats", tag = "v0.5.0" }

[r-dependencies]
packages = ["dplyr", "yardstick"]

[py-dependencies]
version = "python313"
packages = ["numpy", "pandas", "scikit-learn", "xgboost"]

[additional-tools]
packages = ["quarto"]

[t]
# Minimum T language version required
min_version = "0.51.2"

\(\)
packages = ["amsmath", "geometry", "hyperref", "biblatex"]

nix run github:b-rodrigues/tlang -- init --package my_package
cd my_package

[package]
name = "my_package"
version = "0.1.0"
description = "A brief description of what my_package does"
authors = ["brodriguesco"]
license = "EUPL-1.2"
homepage = ""
repository = ""

[dependencies]
# T packages this package depends on
# Format: package = { git = "repository-url", tag = "version" }

[t]
# Minimum T language version required
min_version = "0.5.0"

nix develop

Agentic AI tools like Claude Code can write and run code, fix its own errors, and produce a formatted report with figures. I wanted to know whether that translates into reliable ecological modelling, so we ran a test: three fisheries tasks, four AI models, ten independent runs each, scored against a rubric. The results are published in Fish and Fisheries.

We found agents can be genuinely useful, but only if you know how to use them well and only if you know enough about the analysis to catch what they miss.

A little bit less than two years ago, building on work by Jim Hester and Kevin Ushey, Davis Vaughan completed a very impactful JavaScript file for the R community: an R grammar for the Tree-sitter parsing generator. He even got a round of applause for it during a talk at the useR! 2024 conference! So, did he get cheered for… grammatical rules in a JavaScript file?

No, the audience was excited about the improved developer experience for R that this file unlocked. R tooling around Tree-sitter is how you get

• reformatting through Air and linting through Jarl;
• auto-completion or help on hover in the Positron IDE;
• better search for R on GitHub;
• and more!

In this post, we’ll explain what Tree-sitter is, and how tools built on Tree-sitter can benefit your R development workflow.

Code parsing: what is Tree-sitter?

Tree-sitter is a code parsing generator written in C, with bindings existing in several languages including Rust (and R!).

Let’s rewind a little bit. What does it mean to parse code?

Basically, given a string of code like

a <- mean(x, na.rm = TRUE)

How do you know that mean is a function name, na.rm an argument name, TRUE a logical? You have to parse that code into what’s called a parse tree. You do that in your head when reading R code.

R itself can obviously parse R code, thanks to its grammar. See for instance the commit that introduced R’s native pipe, which necessitated extending R’s syntax thus modifying its grammar.

You can use parse() and getParseData() to parse R code.

parse(
 text = "a <- mean(x, na.rm = TRUE)",
 keep.source = TRUE
) |>
 getParseData()
#> line1 col1 line2 col2 id parent token terminal text
#> 23 1 1 1 26 23 0 expr FALSE 
#> 1 1 1 1 1 1 3 SYMBOL TRUE a
#> 3 1 1 1 1 3 23 expr FALSE 
#> 2 1 3 1 4 2 23 LEFT_ASSIGN TRUE <-
#> 21 1 6 1 26 21 23 expr FALSE 
#> 4 1 6 1 9 4 6 SYMBOL_FUNCTION_CALL TRUE mean
#> 6 1 6 1 9 6 21 expr FALSE 
#> 5 1 10 1 10 5 21 '(' TRUE (
#> 7 1 11 1 11 7 9 SYMBOL TRUE x
#> 9 1 11 1 11 9 21 expr FALSE 
#> 8 1 12 1 12 8 21 ',' TRUE ,
#> 13 1 14 1 18 13 21 SYMBOL_SUB TRUE na.rm
#> 14 1 20 1 20 14 21 EQ_SUB TRUE =
#> 15 1 22 1 25 15 16 NUM_CONST TRUE TRUE
#> 16 1 22 1 25 16 21 expr FALSE 
#> 17 1 26 1 26 17 21 ')' TRUE )

Or you could transform that same data into XML using Gábor Csárdi’s {xmlparsedata}:

parse(
 text = "a <- mean(x, na.rm = TRUE)",
 keep.source = TRUE
) |>
 xmlparsedata::xml_parse_data(pretty = TRUE) |>
 xml2::read_xml() |>
 as.character() |>
 cat()
#> <?xml version="1.0" encoding="UTF-8" standalone="yes"?>
#> <exprlist>
#> <expr line1="1" col1="1" line2="1" col2="26" start="28" end="53">
#> <expr line1="1" col1="1" line2="1" col2="1" start="28" end="28">
#> <SYMBOL line1="1" col1="1" line2="1" col2="1" start="28" end="28">a</SYMBOL>
#> </expr>
#> <LEFT_ASSIGN line1="1" col1="3" line2="1" col2="4" start="30" end="31"><-</LEFT_ASSIGN>
#> <expr line1="1" col1="6" line2="1" col2="26" start="33" end="53">
#> <expr line1="1" col1="6" line2="1" col2="9" start="33" end="36">
#> <SYMBOL_FUNCTION_CALL line1="1" col1="6" line2="1" col2="9" start="33" end="36">mean</SYMBOL_FUNCTION_CALL>
#> </expr>
#> <OP-LEFT-PAREN line1="1" col1="10" line2="1" col2="10" start="37" end="37">(</OP-LEFT-PAREN>
#> <expr line1="1" col1="11" line2="1" col2="11" start="38" end="38">
#> <SYMBOL line1="1" col1="11" line2="1" col2="11" start="38" end="38">x</SYMBOL>
#> </expr>
#> <OP-COMMA line1="1" col1="12" line2="1" col2="12" start="39" end="39">,</OP-COMMA>
#> <SYMBOL_SUB line1="1" col1="14" line2="1" col2="18" start="41" end="45">na.rm</SYMBOL_SUB>
#> <EQ_SUB line1="1" col1="20" line2="1" col2="20" start="47" end="47">=</EQ_SUB>
#> <expr line1="1" col1="22" line2="1" col2="25" start="49" end="52">
#> <NUM_CONST line1="1" col1="22" line2="1" col2="25" start="49" end="52">TRUE</NUM_CONST>
#> </expr>
#> <OP-RIGHT-PAREN line1="1" col1="26" line2="1" col2="26" start="53" end="53">)</OP-RIGHT-PAREN>
#> </expr>
#> </expr>
#> </exprlist>

In both cases, you recognize words such as LEFT_ASSIGN or SYMBOL_FUNCTION_CALL. Parsing is an essential step before the code is actually executed, but parsed code can also be used for other purposes, such as analyzing code without brittle regular expressions (does it call a particular function?), navigating code (going from a function call to the definition of that function), or modifying code (replacing all occurrences of a function with another one).

Now, Tree-sitter performs this same code parsing but faster especially thanks to its support of incremental parsing – which is key to updating the syntax tree as you are typing in your editor for instance! Tree-sitter is agnostic in that it can parse any code as long as there is a grammar for it (think, Rosetta Stone plugins). It’s been used for many languages which means many tools have been built around it.

To have Tree-sitter “learn” a new language you need to give it a file containing the definition of the syntax of that language, what’s called a grammar. This is where the aforementioned JavaScript file by Davis Vaughan and collaborators comes into play! The treesitter-r repo, which provides a translation of the R grammar in the format expected by Tree-sitter, is the base of all tools presented in this post which use R code as their input.

Here’s how to use the {treesitter} R package for the same code as earlier. The {treesitter} R package allows us to use Tree-sitter from R. To parse R code with it, we need the language() function from {treesitter.r}¹.

library(treesitter)
#> 
#> Attaching package: 'treesitter'
#> The following object is masked from 'package:base':
#> 
#> range
language <- treesitter.r::language()
parser <- parser(language)
text <- "a <- mean(x, na.rm = TRUE)"
parser_parse(parser, text)
#> <tree_sitter_tree>
#> 
#> ── Text ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
#> a <- mean(x, na.rm = TRUE)
#> 
#> ── S-Expression ────────────────────────────────────────────────────────────────────────────────────────────────────────────────
#> (program [(0, 0), (0, 26)]
#> (binary_operator [(0, 0), (0, 26)]
#> lhs: (identifier [(0, 0), (0, 1)])
#> operator: "<-" [(0, 2), (0, 4)]
#> rhs: (call [(0, 5), (0, 26)]
#> function: (identifier [(0, 5), (0, 9)])
#> arguments: (arguments [(0, 9), (0, 26)]
#> open: "(" [(0, 9), (0, 10)]
#> argument: (argument [(0, 10), (0, 11)]
#> value: (identifier [(0, 10), (0, 11)])
#> )
#> (comma [(0, 11), (0, 12)])
#> argument: (argument [(0, 13), (0, 25)]
#> name: (identifier [(0, 13), (0, 18)])
#> "=" [(0, 19), (0, 20)]
#> value: (true [(0, 21), (0, 25)])
#> )
#> close: ")" [(0, 25), (0, 26)]
#> )
#> )
#> )
#> )

Tree-sitter is the workhorse of many tools, that are mentioned in the diagram below. All of them are dependent on tree-sitter and the R grammar provided to it. Some of them are command-line interfaces (CLIs), while others are R packages.

Browsing code interactively: Positron IDE, GitHub

The real reason why the audience applauded Davis Vaughan is that he explained how the R grammar for Tree-sitter had been deployed to GitHub so that we get almost as good experience browsing R code on GitHub as browsing, say, JS code. If we search for a function name in a repository for instance, its definition will be indicated in the search results. See Davis’ slides (also available in PDF), or refer to the video below showing how typing vetiver_model in the search bar from the R vetiver repo makes the function definition the first result, on which one can click to land into the definition.

Also very useful is the use of Tree-sitter by Ark, the R kernel used in the Positron IDE. Ark is how you get autocompletion and help on hover in Positron. The video below shows how you can extend the selection to further steps of a pipeline in Positron.

This use case of Tree-sitter is also featured in Davis’ slides. See also Lionel Henry’s and Davis Vaughan’s talk about Ark at posit conf 2024, especially the part about code assistance.

Other development environments such as Emacs also have support for Tree-sitter.

Searching/browsing code

You can parse and search R code using the {treesitter} R package and treesitter query syntax. The {treesitter} R package is a dependency of the {gander} package by Simon Couch, that is meant to be used for a better experience with LLMs when writing R code. Another use case of the {treesitter} R package is the {igraph.r2cdocs} extension to {roxygen2} for the {igraph} package, that parses all of igraph R code to then be able to identify, for each exported function, whether it (in)directly calls a function whose name ends with _impl, indicating a wrapper to a C igraph function whose docs can be then be linked from the manual of the R function.

The {pkgdepends} package calls Tree-sitter (C) to detect dependencies in files. Below we run it on the source of the saperlipopette R package.

pkgdepends::scan_deps(
 "../../../../../CHAMPIONS/saperlipopette",
 "../../../../../CHAMPIONS"
)
#> 
#> Dependencies:
#> + brio  @ R/blame.R, R/check-editor.R, R/clean-dir.R, R/committed-to-main.R, R/committed-to-wrong-branch.R, R/conflict…
#> + cli  @ inst/exo_bisect-Rprofile.en.R, inst/exo_bisect-Rprofile.es.R, inst/exo_bisect-Rprofile.fr.R, inst/exo_blame-…
#> + devtools  @ saperlipopette.Rproj
#> + fs  @ R/blame.R, R/check-editor.R, R/clean-dir.R, R/committed-to-main.R, R/committed-to-wrong-branch.R, R/conflict…
#> + gert  @ inst/exo_check_editor-Rprofile.en.R, inst/exo_check_editor-Rprofile.es.R, inst/exo_check_editor-Rprofile.fr.…
#> + knitr  @ README.Rmd
#> + parsedate  @ R/utils-git.R
#> + purrr  @ R/create-all.R, R/debug.R, R/log-deleted-file.R, R/log-deleted-line.R, R/revparse.R, R/roxygen2.R, R/worktre…
#> + rlang  @ R/create-all.R, R/roxygen2.R, R/utils-fs.R, R/utils-usethis.R, R/zzz.R
#> + rmarkdown  @ README.Rmd, vignettes/saperlipopette.qmd
#> + roxygen2  @ R/roxygen2.R, saperlipopette.Rproj
#> + saperlipopette @ README.Rmd, vignettes/saperlipopette.qmd
#> + tibble  @ R/roxygen2.R
#> + usethis  @ R/blame.R, R/check-editor.R, R/clean-dir.R, R/committed-to-main.R, R/committed-to-wrong-branch.R, R/conflict…
#> + vctrs  @ R/roxygen2.R
#> + withr  @ R/blame.R, R/check-editor.R, R/clean-dir.R, R/committed-to-main.R, R/committed-to-wrong-branch.R, R/conflict…
#> 
#> Test dependencies:
#> + fs  @ tests/testthat/test-blame.R, tests/testthat/test-check-editor.R, tests/testthat/test-clean-dir.R, tests/test…
#> + gert  @ tests/testthat/test-blame.R, tests/testthat/test-clean-dir.R, tests/testthat/test-committed-to-main.R, tests…
#> + rlang  @ tests/testthat/test-blame.R, tests/testthat/test-check-editor.R, tests/testthat/test-clean-dir.R, tests/test…
#> + saperlipopette @ tests/testthat.R
#> + testthat  @ tests/testthat.R
#> + withr  @ tests/testthat/test-blame.R, tests/testthat/test-check-editor.R, tests/testthat/test-clean-dir.R, tests/test…

ast-grep is a useful tool built on Tree-sitter for searching and re-writing code, with a clearer query syntax than Tree-sitter’s. Its name is reminiscent of grep, but with ast-grep we do not need to write brittle regular expressions . {astgrepr} by Etienne Bacher is an R wrapper to the Rust bindings of ast-grep, and is used in Etienne’s {flir} package for refactoring code.

The ast-grep command-line interface (CLI) itself is featured in a useful blog post by Emil Hvitfeldt where he explains how to document the usage of ast-grep for Claude.

Formatting and linting: Air, Jarl

Speaking of CLIs…

Air, by Davis Vaughan and Lionel Henry, is a CLI built on Tree-sitter, in Rust. It reformats code blazingly fast.

Jarl, by Etienne Bacher, is a CLI built on Air, therefore also on Tree-sitter, in Rust. It lints and fixes code, also blazingly fast. It can even detect unreachable code, unused functions and duplicated function definitions.

In both of these examples, the creation of CLIs wrapping Rust bindings was more efficient than the creation of R packages wrapping the {treesitter} R package, for several reasons:

• Rust CLIs can edit code very fast²;
• CLIs are integrated in extensions for popular IDEs (for instance Positron);
• a CLI is easier to install on CI than an R package that needs, well, an R installation.

More tools

A brief mention of some other interesting tools we’ve explored a bit less.

Configuring: {ts} for parsing JSON and TOML (not R!)

The {ts} package by Gábor Csárdi is the backbone of two R packages used for editing and manipulating:

• TOML {tstoml};
• JSON {tsjson}.

Compared to existing parsers in R for those formats, these two packages preserve comments.

Testing code: {muttest}

Mutation testing is a kind of testing where you, say, randomly swap + with - in your code (you mutate it) and you run your tests to see whether they catch the mutant. The {muttest} package by Jakub Sobolewski is an R package for mutation testing, that depends on the {treesitter} R package.

Diffing code: difftastic

The difftastic CLI by Wilfred Hughes is “a structural diff tool that understands syntax”. This means that difftastic doesn’t only compare line or “words” but actual syntax by looking at lines around the lines that changed (by default, 3). Even better, it understands R out of the box. See this blog post with examples of R code diffing.

Conclusion: more to come?

In this post, we’ve presented an overview of Tree-sitter based tooling for R or in R.

Note that this ecosystem of tools is very actively developed, so some tools might come and go. However, the idea that plugging the R grammar into a general parsing generator brings cool features to us R developers, will remain true. Maybe you will contribute to this ecosystem, either through an existing tool or by creating a new one?

IMPORTANT: The website https://www.techtonique.net is down until further notice.

https://www.techtonique.net contained an language-agnostic API for machine learning tasks (classification, regression, survival analysis, forecasting etc.).

As a result, do not buy the Gumroad tutorial then.

You can still use the packages https://github.com/Techtonique locally.

PS: It’s not an April’s fool joke.

Disclaimer: This blog contains opinions that are of the authors alone and do not necessarily reflect the strategy of their respective organizations.

install.packages("dataviewR")

SEX == "F" & AGE > 65 & TRT01P == "Xanomeline High Dose"

library(dataviewR)
library(pharmaverseadam)

dataviewer(adsl, adlb)

@online{pujari2026,
  author = {Pujari, Siddhesh and Kumar N, Madhan and S, Gomathi and
    Haight, Mackenzie},
  title = {Meet {dataviewR:} {The} {View()} {You} {Always} {Wanted}},
  date = {2026-03-31},
  url = {https://pharmaverse.github.io/blog/posts/2026-03-29-meet-dataviewr-the/meet-dataviewr-the-view-you-always-wanted.html},
  langid = {en}
}

The Emergency departments leading the transformation of Alzheimer’s and dementia care (ED-LEAD) study, which I have written about in the past, is approaching the end of its third year. This multifactorial design evaluates three independent, yet potentially synergistic, interventions aimed at improving care for persons living with dementia (PLWD) and their caregivers.

To estimate intervention effects, we are using what I’ve called the HEx-factor model, a Bayesian hierarchical exchangeable factorial model. The original plan was to conduct all analyses using Stan. However, we’ve run into a bit of a snafu. I’ve been working through the problem, and thought I’d share here.

The challenge turns out to be a computational one. Because the ED-LEAD analyses must be conducted on National Institute on Aging (NIA) Data LINKAGE servers, we are working in a somewhat restricted software environment, at least with respect to Bayesian data analysis. In particular, we have not been able to install or run Stan, which was our analytic engine of choice. This forced us to consider alternatives, and we turned to JAGS, which is available in the Linkage environment and certainly is well-suited for Bayesian hierarchical modeling.

At first glance, this might seem like a straightforward substitution. Both Stan and JAGS allow us to specify the same likelihood and priors. However, I quickly noticed that the models were not performing as well in JAGS as they had in Stan. It turns out that the samplers used in JAGS are more sensitive to posterior dependence than the Hamiltonian Monte Carlo (HMC) methods implemented in Stan.

I set out to understand and fix the problem, and found that a simple reparameterization—re-coding the binary treatment indicators—made a substantial difference. With this change, the JAGS sampler was able to explore the posterior distribution much more efficiently, yielding results comparable to those obtained with Stan.

To understand why this happens, I ran a series of simple simulations comparing the original and reparameterized versions of a basic two-way factorial model. That is what I present here.

library(simstudy)
library(data.table)
library(ggplot2)
library(rjags)
library(coda)
library(posterior)
library(broom)
library(gt)

RNGkind("Mersenne-Twister", "Inversion", "Rejection")
set.seed(824)

s_gen <- function(n = 2000,
                    alpha = -0.8,
                    beta_a = 0.5,
                    beta_b = 0.9,
                    beta_ab = -0.3) {
  
  def <- 
    defData(varname = "A", formula = 0.5, dist = "binary") |>
    defData(varname = "B", formula = 0.5, dist = "binary") |>
    defData(varname = "AB", formula = "A*B", dist = "nonrandom") |>
    defData(varname = "A_c", formula = "A - 0.5", dist = "nonrandom") |>
    defData(varname = "B_c", formula = "B - 0.5", dist = "nonrandom") |>
    defData(varname = "AB_c", formula = "A_c * B_c", dist = "nonrandom") |>
    defData(
      varname = "Y", 
      formula = "..alpha + ..beta_a * A + ..beta_b * B + ..beta_ab * AB",
      dist = "binary", link = "logit"
    )
    
  genData(n, def)
  
}

dd <- s_gen()

fit_01 <- glm(Y ~ A * B, data = dd, family = binomial)
fit_c  <- glm(Y ~ A_c * B_c, data = dd, family = binomial)

tidy(fit_01)
## # A tibble: 4 × 5
##   term        estimate std.error statistic  p.value
##   <chr>          <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept)   -1.03     0.0995    -10.4  3.90e-25
## 2 A              0.835    0.135       6.18 6.34e-10
## 3 B              1.14     0.134       8.46 2.73e-17
## 4 A:B           -0.670    0.186      -3.61 3.11e- 4
tidy(fit_c)
## # A tibble: 4 × 5
##   term        estimate std.error statistic  p.value
##   <chr>          <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept)   -0.212    0.0464     -4.57 4.91e- 6
## 2 A_c            0.501    0.0929      5.39 7.05e- 8
## 3 B_c            0.802    0.0929      8.63 6.07e-18
## 4 A_c:B_c       -0.670    0.186      -3.61 3.11e- 4

model_01 <- "
model {
  for (i in 1:N) {
    Y[i] ~ dbern(p[i])
    logit(p[i]) <- alpha + beta_a * A[i] + beta_b * B[i] + beta_ab * AB[i]
  }
  
  alpha   ~ dnorm(0, 0.25)
  beta_a  ~ dnorm(0, 0.25)
  beta_b  ~ dnorm(0, 0.25)
  beta_ab ~ dnorm(0, 25)
}
"

model_c <- "
model {
  for (i in 1:N) {
    Y[i] ~ dbern(p[i])
    logit(p[i]) <- alpha + gamma_a * A_c[i] + gamma_b * B_c[i] + gamma_ab * AB_c[i]
  }
  
  alpha   ~ dnorm(0, 0.25)
  gamma_a  ~ dnorm(0, 0.25)
  gamma_b  ~ dnorm(0, 0.25)
  gamma_ab ~ dnorm(0, 25)
}
"

fit_jags <- function(dat, model_string, centered = FALSE,
                     n_chains = 3, burn = 2000, n_iter = 5000) {
  
  # jdat <- as.list(dat[, .(Y, A, B, AB, A_c, B_c, AB_c)])
  if (centered) {
    jdat <- as.list(dat[, .(Y, A_c, B_c, AB_c)])
    vars <- c("alpha", "gamma_a", "gamma_b", "gamma_ab")
  } else {
    jdat <- as.list(dat[, .(Y, A, B, AB)])
    vars <- c("alpha", "beta_a", "beta_b", "beta_ab")
  }
  jdat$N <- nrow(dat)
  
  mod <- jags.model(
    textConnection(model_string),
    data = jdat,
    n.chains = n_chains,
    quiet = TRUE
  )
  
  update(mod, burn, progress.bar = "none")
  
  samp <- coda.samples(
    mod,
    variable.names = vars,
    n.iter = n_iter,
    progress.bar = "none"
  )
  
  samp
}

samp_01 <- fit_jags(dd, model_01, centered = FALSE)
samp_c  <- fit_jags(dd, model_c, centered = TRUE)

diag_tbl <- function(samp, model_name) {
  post <- as_draws_df(samp)
  summ <- summarise_draws(post)
  out <- as.data.table(summ)
  out[, model := model_name]
  out[]
}

diag_01 <- diag_tbl(samp_01, "0/1-coded")
diag_c  <- diag_tbl(samp_c, "centered")

get_lor_summary <- function(samp, model_name) {
  dt <- as.data.table(as_draws_df(samp))
  
  if (model_name == "0/1-coded") {
    dt[, lOR_A := beta_a]
    dt[, lOR_B := beta_b]
    dt[, lOR_AB := beta_a + beta_b + beta_ab]
  } else {
    dt[, lOR_A := gamma_a - 0.5 * gamma_ab]
    dt[, lOR_B := gamma_b - 0.5 * gamma_ab]
    dt[, lOR_AB := gamma_a + gamma_b]
  }
  
  dt[, .(
    mean_A = mean(lOR_A),
    mean_B = mean(lOR_B),
    mean_AB = mean(lOR_AB),
    sd_A = sd(lOR_A),
    sd_B = sd(lOR_B),
    sd_AB = sd(lOR_AB)
  )]
}

lor_01 <- get_lor_summary(samp_01, "0/1-coded")
lor_c  <- get_lor_summary(samp_c,  "centered")

one_run <- function(
  n = 2000,
  truth = c(alpha = -0.8, beta_a = 0.5, beta_b = 0.9, beta_ab = -0.3),
  n_chains = 3,
  burn = 1000,
  n_iter = 3000
) {
  
  dd <- s_gen(
    n = n,
    alpha = truth["alpha"],
    beta_a = truth["beta_a"],
    beta_b = truth["beta_b"],
    beta_ab = truth["beta_ab"]
  )
  
  samp_01 <- fit_jags(
    dd, model_01, centered = FALSE, 
    n_chains = n_chains, burn = burn, n_iter = n_iter
  )
  
  samp_c <- fit_jags(
    dd, model_c, centered = TRUE, 
    n_chains = n_chains, burn = burn, n_iter = n_iter)
  
  get_metrics <- function(samp, model_name) {
    post <- as_draws_df(samp)
    summ <- as.data.table(summarise_draws(post))
    summ[, model := model_name]
    summ[]
  }
  
  out <- rbindlist(list(
    get_metrics(samp_01, "0/1-coded"),
    get_metrics(samp_c,  "centered")
  ))
  
  out[]
}

nsim <- 500

sim_res <- rbindlist(mclapply(seq_len(nsim), function(i) {
  out <- one_run()
  out[, sim := i]
  out[]
}, mc.cores = 5))

I’m currently on a quest to better know and understand treesitter-based tooling for R. To make it short, treesitter is a tool for parsing code, for instance recognizing what is a function, an argument, a logical in a string of code. With tools built upon treesitter you can search, reformat, lint and fix, etc. your code. Exciting stuff, running locally and deterministically on your machine.

Speaking of “etc.”, Etienne Bacher helpfully suggested I also look at treesitter-based tooling for other languages to see what’s still missing in our ecosystem. This is how I stumbled upon difftastic by Wilfred Hughes, “a structural diff tool that understands syntax”. This means that difftastic doesn’t only compare line or “words” but actual syntax by looking at lines around the lines that changed (by default, 3), Even better, it understands R out of the box¹.

Many thanks to Etienne Bacher not only for making me discover difftastic but also for useful feedback on this post!

Installing difftastic

To install difftastic I downloaded a binary file for my system from the releases of the GitHub repository, as documented in the manual.

difftastic on two files

You can run difftastic on two files, a bit like you would use the waldo R package on two objects.

Let’s compare:

a <- gsub("bad", "good", x)

to

a <- stringr::str_replace(x, "bad", "good")

respectedly saved in old.R and new.R. The CLI is called difft not difftastic. I use the “inline” display rather than the two columns default in order to save horizontal space.

difft old.R new.R --display inline

We’d get to this nice looking diff:

The parentheses and "bad" and "good" arguments are ignored.

We can also get the JSON version of this diff, which is an unstable feature which usage requires setting an environment variable:

export DFT_UNSTABLE=yes
difft old.R new.R --display json

This gets us

{"aligned_lines":[[0,0],[1,1]],"chunks":[[{"lhs":{"line_number":0,"changes":[{"start":5,"end":9,"content":"gsub","highlight":"normal"},{"start":23,"end":24,"content":",","highlight":"normal"},{"start":25,"end":26,"content":"x","highlight":"normal"}]},"rhs":{"line_number":0,"changes":[{"start":5,"end":12,"content":"stringr","highlight":"normal"},{"start":12,"end":14,"content":"::","highlight":"keyword"},{"start":14,"end":25,"content":"str_replace","highlight":"normal"},{"start":26,"end":27,"content":"x","highlight":"normal"},{"start":27,"end":28,"content":",","highlight":"normal"}]}}]],"language":"R","path":"content/post/2026-03-26-difftastic/new.R","status":"changed"}

Now, none of this isn’t very useful because I would never compare files in this way… I use version control!

difftastic with Git

We can set difftastic as the external diff tool for Git globally or for the current project.

For instance with the gert R package, to set it locally:

gert::git_config_set("diff.external", "difft")

If I want to use the inline display I’d set:

gert::git_config_set("diff.external", "difft --display inline")

Then git diff will by default use difftastic. Most interestingly for me, git show --ext-diff will use difftastic. I never use git diff directly but I do look at more or less recent commits a lot.

Say I am interested in the commit that removed roxygen2’s dependency on stringi, I’ll run:

git show 7a1dd39866699a2b0a034bb15244c07698a1e2e7 --ext-diff

and get:

This isn’t spectacular because this is a small diff, but I enjoy the highlighting of the parentheses of the removed nested call, and of the logical.

Cool features of difftastic

Building on two examples of the difftastic homepage…

Ignoring formatting changes

Since formatters can so helpfully apply your formatting preferences, reviewing formatting changes in a patch that’s about something else entirely is useless and annoying. Imagine having a function definition that fits on a single line, then adding one argument to it.

Going from

f <- function(myarg1 = foo, myarg2 = bar) {}

to

f <- function(
  myarg1 = foo,
  myarg2 = bar,
  myarg3 = baz
) {}

Because the definition is now longer than 80 characters, your formatter might switch the definition to be on multiple lines. But the actually interesting change is the addition of one argument.

Native Git diff² would show:

Git with difftastic would show:

The matching of delimiters is why I found the difftastic’s display of the roxygen2 commit more pleasing.

Matching delimiters in wrappers

The Git diff can look a bit ugly when you simply move code from one function to the other.

Say we go from

f <- function() {
  1 + 1
}

to

f <- function() {
  g()
}

g <- function() {
  1 + 1
}

Git diff would show:

Whereas Git with difftastic would show:

Will I use difftastic?

I really like the concept behind difftastic and the few Git commits I looked at with it rendered nicely. Now, what’s missing for me to use difftastic a lot is its integration with the tools where I actually use Git:

• Positron including the GitLens extension;
• GitHub Pull Request Files tab.

In any case, I’ll continue learning about tools based on treesitter, some of which like Air and Jarl I can already use directly from my IDE.

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

rOpenSci HQ

rOpenSci Dev Guide 1.0.0: Trilingual and Improved

rOpenSci Software Peer Review’s guidance is gathered in an online book that keeps improving! It is now available in English, Spanish and Portuguese. Read more in the release announcement

Champions Program Update

We are still going through the Champions selection process, and we’re excited to share that the new group of mentors has already been selected and is now actively reviewing Champions applications.

This cohort brings together a wonderful mix of returning Champions stepping into mentorship roles, mentors continuing their contributions, and new members joining the program. The 2026 mentors are Andrea Gómez Vargas, Pablo Paccioretti, Alber Hamersson Sánchez Ipia, Erick Isaac Navarro Delgado, Francisco Cardozo, Luis Verde Arregoitia, Monika Ávila Márquez, Guadalupe Pascal, Pao Corrales, and Elio Campitelli. Together, they represent a diverse and vibrant community across Colombia, Mexico, Argentina, Brazil, and Bolivia, with some currently based in Switzerland, Canada, the United States, and Australia. We’re very happy to see this growing, interconnected network supporting the next cohort of Champions.

R-Universe update

You can now download artifacts and log files from R-Universe without being logged in with a GitHub account, for example https://ropensci.r-universe.dev/opencv#checktable.

Software review and usage of AI tools

Authors submitting new software for peer review are now required to explain potential usage of generative AI tools in their package development. All submission templates now include a mandatory check-box:

- [ ] Generative AI tools were used to produce some of the material in this submission.
If so, please describe usage, and include links to any relevant aspects of your repository.

This is the start of our updates to accommodate generative AI tools in package development, as described in our recent blog post. The next phase will involve updates to our Dev Guide, explaining requirements and recommendations for authors, reviewers, and editors. All updates are intended to permit generative AI tools to be used in any useful way, while minimising the burden on those who volunteer their own time to keep our software peer review service running.

Software review bot updates

The ropensci-review-bot now delivers an initial report to all new software pre-submissions and submissions, identifying the five most similar packages from both all rOpenSci packages, and all CRAN packages. The matches are generated by our ropensci-review-tools/pkgmatch package (itself reviewed in this review issue). Matching is based on an “term frequency-inverse document frequency” algorithm, using inverse document frequencies from all rOpenSci and CRAN packages. Similar package reports can also be manually triggered (by editors only) with @ropensci-review-bot similar packages, like in this example for the pkgmatch package itself.

Coworking

Read all about coworking!

• Tuesday April 7th 2026, 9:00 Americas Pacific (16:00 UTC) “Getting to know the CSID Network” with Steffi LaZerte and cohosts Irene Ramos and Adamu Saleh Saidu.
  • Learn more about the CSID Network
  • Meet cohosts, Irene Ramos and Adamu Saleh Saidu, and learn more about the CSID Network and how you might get involved.
• Tuesday May 5th 2026, 9:00 Australia Western (01:00 UTC) “Code Review with rOpenSci” with Steffi LaZerte and cohost Liz Hare.
  • Explore resources for Code Review
  • Sign up to volunteer to do software peer-review at rOpenSci
  • Meet cohost, Liz Hare, and discuss resources for Code Review with rOpenSci.

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

Software

New packages

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

• suwo, developed by Marcelo Araya-Salas together with Jorge Elizondo-Calvo and Alejandro Rico-Guevara: Streamline searching/downloading of nature media files (e.g. audios, photos) from online repositories. The package offers functions for obtaining media metadata from online repositories, downloading associated media files and updating data sets with new records. It has been reviewed by Eric R. Scott and Hugo Gruson.

Discover more packages, read more about Software Peer Review.

New versions

The following eleven packages have had an update since the last newsletter: cffr (v1.3.0), pkgmatch (v0.5.2), tarchetypes (0.14.1), rgbif (v3.8.5), saperlipopette (v0.1.1), gutenbergr (v0.5.0), trud (v0.2.1), naijR (v0.7.0), sasquatch (v0.1.3), lingtypology (v1.1.25), and rerddap (v1.2.3).

Post on dfms release: Releasing dfms 1.0: Fast and Feature-Rich Estimation of Dynamic Factor Models in R.

Software Peer Review

There are fifteen recently closed and active submissions and 5 submissions on hold. Issues are at different stages:

• One at ‘6/approved’:

  • suwo, Access Nature Media Repositories Through R. Submitted by Marcelo Araya-Salas.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

• Three at ‘1/editor-checks’:

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

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

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

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

On the blog

Software Review

• Software Review in the Era of AI: What We Are Testing at rOpenSci by Mark Padgham, Noam Ross, Maëlle Salmon, Yanina Bellini Saibene, Mauro Lepore, Emily Riederer, Jouni Helske, and Francisco Rodriguez-Sanchez. rOpenSci is testing preliminary policies on the use of generative AI tools, with proposed updates to documentation and procedures for authors submitting software for review, for editors, and for reviewers.

• rOpenSci Dev Guide 1.0.0: Trilingual and Improved by Maëlle Salmon, Mark Padgham, and Noam Ross. Updates in version 1.0.0 of the online book ‘rOpenSci Packages: Development, Maintenance, and Peer Review’. Other languages: rOpenSci Dev Guide 1.0.0: Trilingüe y mejorada (es), Guia de desenvolvimento da rOpenSci 1.0.0: trilíngue e aprimorado (pt).

• Breaking Release of the patentsview R Package by Russ Allen and Chris Baker. Breaking Release of the patentsview R Package.

Calls for contributions

Calls for maintainers

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

• NLMR, R package to simulate neutral landscape models. Issue for volunteering.

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

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

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

Calls for contributions

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

Package development corner

Some useful tips for R package developers.

A new R core member!

The R Foundation announced that Heather Turner has joined the R Core Team!

How to browse the R mailing lists

The official mailing lists of the R project like R-package-devel are full of important and useful information. How to browse them, given that the default website is not easy to search? You can use the mail-archive website (thanks to Hugo Gruson for the reminder!) or a new project by James Balamuta: the R Mailing Lists Archive!

“Claude Code: Setting up ast-grep with R support”

Thanks to Mauro Lepore for sharing this blog post by Emil Hvitfeldt: “Claude Code: Setting up ast-grep with R support”. ast-grep is a tool for querying code by syntax rather than brittle regular expressions. The blog post describes how to add R support to this tool, and how to take advantage of it when using Claude.

On muffling messages from packages

A follow-up on our post “Please Shut Up! Verbosity Control in Packages”.

• With the {cli} R package you can change the default handler for messages. See the docs. It seems mostly used to muffle messages, e.g. in flir.
• Here’s how the usethis R package muffles gert message selectively.

Last words

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

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

With the conflict in Iran causing worldwide disruption to energy markets, I have both a work and personal interest in energy supply in Pacific islands, which led me to this blog post. Here I look at just two aspects of energy: electricity generation, and household cooking. Nothing fancy here, just accessing some data and drawing a couple of plots.

Electricity generation

Here is the source of electricity for Pacific island countries, plus Australia and New Zealand, collated by Our World In Data from Energy Institute data that ultimately comes from government estimates:

There’s a pretty obvious story here: most of the Pacific is very dependent on “oil” (in the form of diesel) for generation of most of its electricity. There are some small steps towards renewables happening in recent years, but the vulnerability to a price or availability shock for diesel is pretty obvious.

Here’s the code for producing that, using the valuable owidapi R package to access the Our World in Data API.

#---------------Set up-----------------
library(owidapi)
library(tidyverse)
library(countrycode)
library(WDI)
library(jsonlite)
library(janitor)
library(httr2)


pic_codes <- 
  c(
    "ASM", "COK", "FSM", "FJI", "PYF", "GUM", "KIR", "MHL", "NRU", "NCL",
    "NIU", "MNP", "PLW", "PNG", "PCN", "WSM", "SLB", "TKL", "TON", "TUV",
    "VUT", "WLF", "AUS", "NZL"
  )
stopifnot(length(pic_codes) == 24)

# visual check we've got the right country codes for the Pacific:
countrycode::countrycode(pic_codes, origin = "iso3c", destination = "country.name.en")

#=======================electricity source===================

palette <- c(
  coal = "brown",
  gas = "magenta",
  oil = "red",
  kerosene = "red",
  electricity = "purple",
  solar = "yellow",
  wind = "steelblue",
  hydro = "darkblue",
  bioenergy = "lightgreen",
  charcoal = "grey",
  biomass = "darkgreen",
  'other renewables' = "darkgreen"
)

#-------------------electricity mix-----------------
elec_mix <- owid_get(
  chart_id = "share-elec-by-source",
  entities  = pic_codes
)

elec_data <- elec_mix |> 
   rename(country = entity_name) |> 
   select(-entity_id) |> 
   gather(variable, value, -country, -year) |>
   filter(value != 0) |> 
   filter(year > 2001) |> 
   mutate(variable = gsub("_share_of_electricity__pct", "", variable, fixed = TRUE),
          variable =gsub("_", " ", variable),
          variable =gsub(" excluding bioenergy", "", variable),
          variable = fct_drop(variable)) |> 
   mutate(variable = fct_relevel(variable, c("bioenergy", "hydro", "other renewables"), after = Inf)) |> 
   mutate(country = fct_relevel(country, c("Australia", "New Zealand"), after = Inf)) |> 
  group_by(country) |> 
  mutate(prop_pc= sum(value[variable %in% c("oil", "gas") & year == max(year)]) 
         / sum(value[year == max(year)])) |> 
  ungroup() |> 
  mutate(country = fct_reorder(country, prop_pc))

# Draw chart
elec_data |> 
  ggplot(aes(x = year, y = value, fill = variable)) +
 facet_wrap(~country, ncol = 5) +
  geom_col() +
  scale_fill_manual(values = palette) +
  scale_y_continuous(label = percent_format(scale = 1)) +
   labs(y = "Percentage of electricity",
        fill = "Source:",
        title = "Share of electricity by source",
        subtitle = "Countries shown in increasing order of vulnerability of electricity to a petrochemicals price or availability crisis.",
        x = "",
        caption = "Source: Ember (2026); Energy Institute - Statistical Review of World Energy (2025). Data processed by Our World In Data.") +
   theme(axis.text.x = element_text(angle = 45, hjust = 1))

Cooking fuel

OK, so electricity generation could be threatened by a lack of diesel. What about household cooking? This next chart draws on the definitive World Health Organization Household Energy Database, which models (based on what household survey data that is available) what households are using to cook:

Again, we see a lot of reliance on petrochemical products, particularly kerosene and liquid natural gas. The latter has been promoted as a relatively clean and healthy fuel to cook with compared to burning biomass (e.g. wood, coconuts, etc).

The larger Melanesian countries, with high rural populations, are those with the greatest use still of biomass for cooking. Most Pacific island countries do most of their cooking with oil or gas derived energy (remembering from the first chart, that ‘electricity’ often means diesel, ultimately).

Here’s the code to produce that chart. I used an LLM (I forget which) for the code to access the API itself, but I tested it and tweaked it to match my style, and the chart of course is all my own code.

#--------------------cooking-------------------
# The definitive source is the WHO  WHO Household Energy Database 
# which draws on various household surveys
# See https://www.who.int/data/gho/data/themes/air-pollution/cooking-fuel-and-technology-database-by-fuel-category

# next half dozen lines of code were supplied by Co-pilot and minimally
# tweaked by me for my style
indicator_code <- "PHE_HHAIR_PROP_POP_CATEGORY_FUELS"  # % by fuel type [3](https://millenniumindicators.un.org/wiki/spaces/SDGeHandbook/pages/35291272/Indicator+7.1.2)
url <- paste0("https://ghoapi.azureedge.net/api/", indicator_code)

resp <- request(url) |> 
  req_headers(`Accept` = "application/json")  |> 
  req_perform()

cooking_data <- fromJSON(resp_body_string(resp), flatten = TRUE)$value |>
  as_tibble() |> 
  clean_names()


pic_cooking_data <-cooking_data |> 
  filter(spatial_dim %in% pic_codes) |> 
  filter(dim1 == "RESIDENCEAREATYPE_TOTL") |> 
  mutate(fuel_type = tolower(gsub("HOUSEHOLDCOOKINGFUEL_FUEL_", "", dim2))) |> 
  mutate(year = as.numeric(time_dimension_value)) |> 
   select(value = numeric_value, 
          iso3_code = spatial_dim,
          value = numeric_value,
          year,
          fuel_type) |> 
  mutate(country = countrycode(iso3_code, origin = "iso3c", destination = "country.name.en"),
         country = gsub("Federated States", "Fed St", country)) |>
  group_by(country) |> 
  mutate(prop_gke = sum(value[fuel_type %in% c("gas", "kerosene", "electricity") & year == max(year)]) 
         / sum(value[year == max(year)])) |> 
  ungroup() |> 
  mutate(country = fct_reorder(country, prop_gke))

# Draw chart
pic_cooking_data |> 
  ggplot(aes(y = value, x = year, fill = fuel_type)) +
  facet_wrap(~country, ncol = 5) +
  # the numbers don't add up to 100 always, due to being modelled estimates
  #, not fully MECE, not counting dual fuels, etc. Good practice advice
  # is to not force them to add to 100%
  geom_area() +
  scale_fill_manual(values = palette) +
  scale_y_continuous(label = percent_format(scale = 1)) +
  labs(title = "Household primary fuel used for cooking",
       subtitle = "Estimates are modelled by WHO, and not adding up to 100% is a known limitation.
Countries shown in increasing order of vulnerability of cooking to a petrochemicals price or availability crisis.",
       x = "",
       fill = "Fuel type:",
       y ="Proportion of households",
       caption = "Source: WHO Household Energy Database")

That’s all, just a quick one today.

Dipping my toes into financial statements — income, balance sheet & cash flow. Still don’t fully get it, but slowly piecing together the story these numbers tell. Warren Buffett makes it look easy Baby steps!

Motivations

I’ve always wanted to learn financial statement, what it means, what it tells us, what Warren Buffett sees in them. Following the book Warren Buffett and the Interpretation of Financial Statements: The Search for the Company with a Durable Competitive Advantage and Datacamp: Analyzing Financial Statement in Python, I’ve made some notes for myself and also create the metrics functions, so that I can use and view them easily in the future. I’ll be honest, I still don’t fully understand it, but at least I can refer back to this as I look at these statements more frequently.

Disclaimer:

This is purely for educational purposes. This is not a financial advice, nor am I a financial advisor. This is a note to myself. If you find any mistakes or error, please let me know. Thanks! Also, there are a lot of information on each section, I won’t be covering all of them, just mostly the metrics from the book and also points I found interesting.

Objectives:

• The Skeleton of Financial Statements
• Let’s Take An Example
• Income Statement
  • Gross Profit Margin
  • Depreciation
  • Interest Payment to Operating Income
  • Income Before Tax
  • Income After Tax
  • Net Earnings
  • Per Share Earning
  • Operating Margin
  • Metrics
• Balance Sheet
  • Current Assets
  • Net Receivables To Gross Sale Ratio
  • The Current Ratio
  • Property, Plant, and Equipment
  • Short Term Debt
  • Long Term Debt
  • Retained Earnings
  • Treasury Stock
  • Return On Shareholders’ Equity
  • Metrics
• Cash Flow
  • Operating Income
  • Capital Expenditure
  • Stock Buyback
  • Metrics
• Combine All Metrics
• Let’s Look At Another Examples
• Oppotunities For Improvement
• Lessons Learnt

The Skeleton of Financial Statements

A financial statement is a formal record that shows a company’s financial activities and position, typically consisting of three core components: the income statement (which shows revenues earned and expenses incurred to calculate profit or loss over a period), the balance sheet (which presents what the company owns as assets, what it owes as liabilities, and the difference between them as equity at a specific point in time), and the cash flow statement (which tracks the actual movement of cash in and out of the business through operating activities, investing activities, and financing activities). The income statement reveals profitability, and the cash flow statement shows liquidity and how money actually moves through the business.

If we were to think of a kid’s lemonaid shop, the income statement would show how much money the shop made from selling lemon tea and how much it spent on ingredients, paying Johnny hourly to sell (salary) to calculate the profit. The balance sheet would list the shop’s assets (like cash in the register, inventory of lemons, sugar, and any equipment) and liabilities (like loans or unpaid bills – money your parent you borrowed from to buy all of the above) to show the net worth of the business at a given moment. The cash flow statement would track the actual cash coming in from customers and going out for expenses, giving insight into whether the shop has enough liquidity to cover its day-to-day operations.

It sounds simple, in the big picture, but these are just the basic skeleton of financial statements. There are many nuances and details that we need to understand to really grasp the story that these statements are telling us.Each section has its own items and some of these items are good at forming different metrics to tell the story of how the lemonaid business is doing. Below is just a snapshot of Apple’s financial statement.

Income Statement

Balance Sheet

Cash Flow

Let’s Take An Example

Let’s go to Alpha Vantage and create a free api key and then pull Apple’s 10 year financial statement and go through as an exercise.

https://www.sec.gov/Archives/edgar/data/320193/000032019325000079/aapl-20250927.htm#i719388195b384d85a4e238ad88eba90a_181

library(httr)
library(jsonlite)
library(tidyverse)

api_key <- "your_api_key_here" # Use .Renviron to be safer like below combine_all code

## Create a function to pull data
get_data <- function(fx,ticker) {
  raw <- GET(paste0(
    "https://www.alphavantage.co/query?function=",fx,
    "&symbol=",ticker,"&apikey=", api_key
  )) %>%
    content(as = "text", encoding = "UTF-8") %>%
    fromJSON()
  
  df <- raw$annualReports |> 
    as_tibble() |> 
    mutate(across(-c(fiscalDateEnding, reportedCurrency), as.numeric)) |>
    mutate(fiscalDateEnding = as.Date(fiscalDateEnding)) |>
    arrange(fiscalDateEnding)

  return(df)
}

## financial statement
income <- get_data("INCOME_STATEMENT","AAPL")
balance <- get_data("BALANCE_SHEET","AAPL")
cashflow <- get_data("CASH_FLOW","AAPL")

Let’s visualize the income statement

Income Statement

Alright, here we go! Income statement tells the story of how the goods and services are doing in the market, how much it costs to produce and sell them, and how much profit is left after all expenses are accounted for. It is a dynamic statement that shows the flow of money over a period of time, typically a quarter or a year. It is like a movie that tells the story of the company’s operations and profitability. Below are some ratios and heuristics of which companies have a durable competitive advantage, according to Warren Buffett.

Gross Profit Margin

\(Gross Profit Margin = Gross Profit / Total Revenue\)

In the book, we are looking for consistent, as a general rule, above 40% to consider a company to have a durable competitive advantage. Gross profit here is basically Total Revenue - Cost of Revenue, does not take into account of R&D, admin etc. Gross profit margin is a measure of how much profit a company makes from its core operations, before accounting for other expenses. A higher gross profit margin indicates that the company has a strong competitive position in the market and is able to generate more profit from its sales.

So on Apple, you can see that in the sales section, there is product and service. I assume product is the hardware and service is like the cloud storage etc. The cost of revenue section has the same sections, that depicts how much cost to make these products/services. Again, remember this is all just about the goods, does not take into account of the R&D, office that manages these, administrative etc, I think.

Depreciation

Apparently this is a non-cash expense that reflects the reduction in value of a company’s assets over time. It is an accounting method used to allocate the cost of tangible assets (like machinery, equipment, buildings) and intangible assets (like patents, copyrights) over their useful lives. Depreciation allows companies to spread out the expense of an asset over several years, rather than recognizing the entire cost in the year it was purchased.

A quick ratio in the book is

\(Depreciation / Gross Profit\)

which should be low ~5-7% which indicates that the company is not heavily reliant on physical assets that may lose value over time. Warren stated that EBITDA (Earnings Before Depreciation, Taxes, and Amortization) is something Wall Street loves, but in Warren’s eyes, depreciation is a real expanse. Why is it that he said EBITDA is something Wall Street love? Because it shows a higher profit by excluding non-cash expenses like depreciation, which can make a company look more profitable than it actually is… interesting.

Interest Payment to Operating Income

\(Interest Payment / Operating Income\)

should be less than 15%. This ratio indicates that the company is not overly burdened by debt and has a healthy balance between its operating income and interest expenses. A lower ratio suggests that the company is generating sufficient operating income to cover its interest payments, which is a positive sign of financial stability.

Operating Income = Total Revenue - Cost of Revenue - Operating Expenses. And Operating Expenses = R&D + Selling, General and Administrative Expenses. Wow, so many terminologies and I still don’t fully understand them, but at least I can refer back to this when I look at these statements.

Income Before Tax

This is a value that Warren uses when he calculates the return that he is getting when he a whole business.

Income After Tax

This is a value is the truth test of a business. What is reported with SEC should reflect the pre-tax income reported on the income statement. Apparently some companies like to report higher values than the the truth? was what the book said

Net Earnings

\(Net Earning / Total Revenue\)

The heuristic for this is we should look for more than 20% which indicates that the company is able to generate a significant amount of profit from its total revenue, hence a long term competitive advantage.

Net Earnings = Operating Income - Interest Expense - Taxes. This is the bottom line of the income statement, which represents the company’s total profit after all expenses have been deducted from total revenue. It is a key indicator of a company’s profitability and financial performance. A higher net earnings figure indicates that the company is generating more profit from its operations, which can be a sign of a strong competitive position in the market.Interest Expense - Taxes`.

Per Share Earning

This is something Warren wants to see consistent increment over time.

This is calculated by: \(EPS = Net Income/Outstanding Share\)

Does that mean if net income is negative, we’ll have negative share?

Operating Margin

\(Operating Margin = Operating Income / Total Revenue\)

This is a measure of a company’s profitability that indicates how much profit it generates from its operations relative to its total revenue. A higher operating margin suggests that the company is more efficient at converting revenue into profit, which can be a sign of a strong competitive position in the market. Warren looks for more than 10% which indicates that the company has a durable competitive advantage. ?????

Metrics

income_metric <- function(df) {
  df_i <- df |>
    mutate(grossProfitMargin = grossProfit / totalRevenue,
           depreciationToGrossProfit = depreciationAndAmortization / grossProfit,
           interestExpenseToOperatingIncome = interestExpense / operatingIncome,
           netEarningMargin = netIncome / totalRevenue,
           operatingMargin = operatingIncome / totalRevenue)
  
  hline <- tribble(
    ~param, ~hline_value, ~color2,
    "grossProfitMargin", 0.4, "red", 
    "depreciationToGrossProfit", 0.07, "blue",
    "interestExpenseToOperatingIncome", 0.15, "blue",
    "netEarningMargin", 0.2, "red",
    "operatingMargin", 0.1, "red"
  ) |>
    mutate(
      param = factor(param),
      ymin = ifelse(color2 == "red", -Inf, hline_value),
      ymax = ifelse(color2 == "red", hline_value, Inf)
    )
  
  columns <- hline$param
  
  plot <- df_i |>
    select(fiscalDateEnding, grossProfitMargin:operatingMargin) |>
    pivot_longer(cols = c(grossProfitMargin:operatingMargin), 
                 names_to = "param", values_to = "values") |>
    mutate(param = factor(param, levels = columns)) |> 
    ggplot(aes(x = fiscalDateEnding, y = values)) +
    geom_rect(
      data = hline,
      aes(ymin = ymin, ymax = ymax, fill = color2),
      xmin = -Inf, xmax = Inf,
      alpha = 0.15,
      inherit.aes = FALSE
    ) +
    geom_line() +
    facet_wrap(. ~ param, scale = "free_y") +
    scale_fill_identity() +
    theme_bw()
  
  return(plot)
}

income_metric(income) 

Net margin

What is EBTDA? It stands for Earnings Before Depreciation, Taxes, and Amortization. It is a measure of a company’s operating performance that excludes non-operating expenses such as depreciation, taxes, and amortization. It is often used as a proxy for cash flow from operations, as it focuses on the core profitability of the business before accounting for non-cash expenses and tax obligations.

Balance Sheet

\(Asset = Liability + Shareholder Equity\)

Reminded me of Full Metal Alchemist famous phrase “tōka kōkan”, equivalent exchange.

Current Assets

Current asset is any asset that can be reasonably expected to be converted into cash within one year. This includes cash and cash equivalents, accounts receivable, inventory, and other short-term assets. Current assets are important because they provide insight into a company’s liquidity and ability to meet its short-term obligations.

Net Receivables To Gross Sale Ratio

\(Net Receivables / Gross Sale\)

If a company consistently has a lower ratio (? 5% or less), it may indicate that the company is efficient in collecting payments from customers and has a lower risk of bad debts. A higher ratio may suggest that the company is having difficulty collecting payments, which could lead to cash flow issues and potential losses from uncollected receivables.

The Current Ratio

\(Current Ratio = Total Assets / Total Liabilities\)

The current ratio is a liquidity ratio that measures a company’s ability to pay off its short-term liabilities with its short-term assets. A current ratio of 1 or higher is generally considered good, indicating that the company has enough assets to cover its liabilities. A current ratio below 1 may indicate that the company may have difficulty meeting its short-term obligations. This makes sense.

Property, Plant, and Equipment

This is the value of a company’s physical assets, such as land, buildings, machinery, and equipment. It is important to consider the value of PPE when evaluating a company’s financial health and potential for growth. A company with a significant amount of PPE may have a competitive advantage in its industry, as it may be able to produce goods or services more efficiently than its competitors. However, it is also important to consider the age and condition of the PPE, as well as any potential liabilities associated with it. In the title of this chapter, it says For Warren Not Having Them Is A Good Thing

Short Term Debt

From Warren’s perspective, when it comes to investing in financial institutions, he’s always shied away from companies who are bigger borrowers in short-term than long-term.

\(Short Term Debt/ Long Term Debt\)

Not really sure what the heuristic threshold is, but let’s use less than 1 as a good indicator?

Long Term Debt

As a general rule, a company with durable competitve edge, will have no or little long term debt to maintain their business operation. Let’s ensure this is not an uptrend on visualization

Retained Earnings

Retained earnings is the portion of a company’s net income that is retained and not distributed as dividends to shareholders. It represents the accumulated profits that a company has reinvested in its business over time. Retained earnings can be used for various purposes, such as funding research and development, expanding operations, paying off debt, or acquiring other companies. It is an important metric for investors to consider when evaluating a company’s financial health and growth potential, as it indicates how much profit the company has generated and how it has been utilized to support its long-term success.

\(Growth Rate = (Ending Retained Earnings/Beginning Retained Earnings)^{1/years}-1\)

Treasury Stock

Treasury stock refers to shares that a company has repurchased from its shareholders. These shares are held in the company’s treasury and are not considered outstanding shares. Treasury stock can be used for various purposes, such as to increase shareholder value, to have shares available for employee compensation plans, or to prevent hostile takeovers. When a company repurchases its own shares, it reduces the number of outstanding shares in the market, which can increase the value of the remaining shares and potentially boost earnings per share (EPS). However, it is important for investors to consider the reasons behind a company’s decision to buy back its own stock and how it may impact the company’s financial health and long-term growth prospects.

Return On Shareholders’ Equity

\(Return On Shareholders' Equity = Net Income / Shareholders' Equity\)

According to the book, high ROSE means come play How about stop and smell the rose? The book did not tell us exactly what the threshold is, but the companies of choice has about ~30-35%

Metrics

balance_metric <- function(df, income) {
  df_b <- df |>
    mutate(currentRatio = totalCurrentAssets / totalCurrentLiabilities,
           netReceivablesToGrossSale = currentNetReceivables / income$grossProfit,
           shortToLongTermDebt = shortTermDebt / longTermDebt,
           growthRateRetainedEarnings = (retainedEarnings / lag(retainedEarnings))^(1/n()) - 1,
           returnOnShareholdersEquity = income$netIncome / totalShareholderEquity)
  
  hline <- tribble(
    ~param, ~hline_value, ~color2,
    "currentRatio", 1, "red", 
    "netReceivablesToGrossSale", 0.05, "blue",
    "shortToLongTermDebt", 1, "blue",
    "growthRateRetainedEarnings", 0.05, "red",
    "returnOnShareholdersEquity", 0.3, "red"
  ) |>
    mutate(
      param = factor(param),
      ymin = ifelse(color2 == "red", -Inf, hline_value),
      ymax = ifelse(color2 == "red", hline_value, Inf)
    )
  
  columns <- hline$param
  
  plot <- df_b |>
    select(fiscalDateEnding, currentRatio:returnOnShareholdersEquity) |>
    pivot_longer(cols = c(currentRatio:returnOnShareholdersEquity), 
                 names_to = "param", values_to = "values") |>
    mutate(param = factor(param, levels = columns)) |> 
    ggplot(aes(x = fiscalDateEnding, y = values)) +
    geom_rect(
      data = hline,
      aes(ymin = ymin, ymax = ymax, fill = color2),
      xmin = -Inf, xmax = Inf,
      alpha = 0.15,
      inherit.aes = FALSE
    ) +
    geom_line() +
    facet_wrap(. ~ param, scale = "free_y") +
    scale_fill_identity() +
    theme_bw()
  
  return(plot)
}

balance_metric(balance, income)

Cash Flow

Operating Income

Cash flow from operating income starts with net income and then add back depreciation and amortization. This is because depreciation and amortization are non-cash expenses that reduce net income but do not actually involve a cash outflow. By adding them back, we can get a better picture of the actual cash generated by the company’s operations.

Capital Expenditure

Also known as CapEx, refers to the funds that a company uses to acquire, upgrade, and maintain physical assets such as property, buildings, technology, equipment, or machinery. CapEx is an important metric for investors to consider when evaluating a company’s financial health and growth potential, as it indicates how much the company is investing in its long-term success.

If we were to look at Apple’s cash flow statement, CapEx is payments for acquisition of proptery, plant, and equipment.

\(Capital Expenditure/Net Earning\)

And the heuristic is ~50% or less.

Stock Buyback

Stock buyback, also known as share repurchase, refers to a company’s practice of buying back its own shares from the open market. This can be done for various reasons, such as to increase shareholder value, to have shares available for employee compensation plans, or to prevent hostile takeovers. When a company repurchases its own shares, it reduces the number of outstanding shares in the market, which can increase the value of the remaining shares and potentially boost earnings per share (EPS).

On the Apple cash flow statement, stock buyback is listed as repurchase of common stocks. Others might use issuance of (retirement) stocks.

Metrics

cashflow_metric <- function(df) {
  df_c <- df |>
    mutate(operatingCashFlow = netIncome + depreciationDepletionAndAmortization,
           capitalExpenditureToNetEarning = capitalExpenditures / netIncome,
           stockBuybackToNetEarning = abs(proceedsFromRepurchaseOfEquity / netIncome))
  
  hline <- tribble(
    ~param, ~hline_value, ~color2,
    "operatingCashFlow", 0, "red", 
    "capitalExpenditureToNetEarning", 0.5, "blue",
    "stockBuybackToNetEarning", 0.5, "red"
  ) |>
    mutate(
      param = factor(param),
      ymin = ifelse(color2 == "red", -Inf, hline_value),
      ymax = ifelse(color2 == "red", hline_value, Inf)
    )
  
  columns <- hline$param
  
  plot <- df_c |>
    select(fiscalDateEnding, operatingCashFlow, capitalExpenditureToNetEarning, stockBuybackToNetEarning) |>
    pivot_longer(cols = c(operatingCashFlow, capitalExpenditureToNetEarning, stockBuybackToNetEarning), 
                 names_to = "param", values_to = "values") |>
    mutate(param = factor(param, levels = columns)) |> 
    ggplot(aes(x = fiscalDateEnding, y = values)) +
    geom_rect(
      data = hline,
      aes(ymin = ymin, ymax = ymax, fill = color2),
      xmin = -Inf, xmax = Inf,
      alpha = 0.15,
      inherit.aes = FALSE
    ) +
    geom_line() +
    facet_wrap(. ~ param, scale = "free_y") +
    scale_fill_identity() +
    theme_bw()
  
  return(plot)
}

cashflow_metric(cashflow)

Use key valuation metrics (P/E, EV/EBITDA, P/B, P/S, etc) to determine how cheap or expensive a stock is.

Combine All Metrics

library(httr)
library(jsonlite)
library(tidyverse)
library(ggpubr)

api_key <- Sys.getenv("avkey")

## Create a function to pull data
get_data <- function(fx,ticker,share=F) {
  raw <- GET(paste0(
    "https://www.alphavantage.co/query?function=",fx,
    "&symbol=",ticker,"&apikey=", api_key
  )) %>%
    content(as = "text", encoding = "UTF-8") %>%
    fromJSON()
  
  if (share==T) {
    df <- raw$annualEarnings |>
      select(fiscalDateEnding, reportedEPS) |>
      mutate(fiscalDateEnding = ymd(fiscalDateEnding))
  } else {
  df <- raw$annualReports |> 
    as_tibble() |>
    mutate(fiscalDateEnding = as.Date(fiscalDateEnding)) |>   
    mutate(across(where(is.character), as.numeric)) |>        
    arrange(fiscalDateEnding)
  }

  return(df)
}

income_metric <- function(df) {
  df_i <- df |>
    mutate(grossProfitMargin = grossProfit / totalRevenue,
           depreciationToGrossProfit = depreciationAndAmortization / grossProfit,
           interestExpenseToOperatingIncome = interestExpense / operatingIncome,
           netEarningMargin = netIncome / totalRevenue,
           operatingMargin = operatingIncome / totalRevenue)
  
  hline <- tribble(
    ~param, ~hline_value, ~color2,
    "grossProfitMargin", 0.4, "red", 
    "depreciationToGrossProfit", 0.07, "blue",
    "interestExpenseToOperatingIncome", 0.15, "blue",
    "netEarningMargin", 0.2, "red",
    "operatingMargin", 0.1, "red"
  ) |>
    mutate(
      param = factor(param),
      ymin = ifelse(color2 == "red", -Inf, hline_value),
      ymax = ifelse(color2 == "red", hline_value, Inf)
    )
  
  columns <- hline$param
  
  plot <- df_i |>
    select(fiscalDateEnding, grossProfitMargin:operatingMargin) |>
    pivot_longer(cols = c(grossProfitMargin:operatingMargin), 
                 names_to = "param", values_to = "values") |>
    mutate(param = factor(param, levels = columns)) |> 
    ggplot(aes(x = fiscalDateEnding, y = values)) +
    geom_rect(
      data = hline,
      aes(ymin = ymin, ymax = ymax, fill = color2),
      xmin = -Inf, xmax = Inf,
      alpha = 0.15,
      inherit.aes = FALSE
    ) +
    geom_line() +
    facet_wrap(. ~ param, scale = "free_y") +
    scale_fill_identity() +
    theme_bw() +
    ggtitle("Income Statement Metrics")
  
  return(plot)
}

share_metric <- function(df) {
    plot <- df |>
      mutate(reportedEPS = as.numeric(reportedEPS)) |>
      ggplot(aes(x = fiscalDateEnding, y = reportedEPS)) +
      geom_line() +
      geom_smooth() +
      theme_bw() +
      ggtitle("Share Metrics")
    return(plot)
}

balance_metric <- function(df, income) {
  df_b <- df |>
    mutate(currentRatio = totalCurrentAssets / totalCurrentLiabilities,
           netReceivablesToGrossSale = currentNetReceivables / income$grossProfit,
           shortToLongTermDebt = shortTermDebt / longTermDebt,
           growthRateRetainedEarnings = (retainedEarnings / lag(retainedEarnings))^(1/n()) - 1,
           returnOnShareholdersEquity = income$netIncome / totalShareholderEquity)
  
  hline <- tribble(
    ~param, ~hline_value, ~color2,
    "currentRatio", 1, "red", 
    "netReceivablesToGrossSale", 0.3, "blue",
    "shortToLongTermDebt", 1, "blue",
    "growthRateRetainedEarnings", 0.05, "red",
    "returnOnShareholdersEquity", 0.3, "red"
  ) |>
    mutate(
      param = factor(param),
      ymin = ifelse(color2 == "red", -Inf, hline_value),
      ymax = ifelse(color2 == "red", hline_value, Inf)
    )
  
  columns <- hline$param
  
  plot <- df_b |>
    select(fiscalDateEnding, currentRatio:returnOnShareholdersEquity) |>
    pivot_longer(cols = c(currentRatio:returnOnShareholdersEquity), 
                 names_to = "param", values_to = "values") |>
    mutate(param = factor(param, levels = columns)) |>
    ggplot(aes(x = fiscalDateEnding, y = values)) +
    geom_rect(
      data = hline,
      aes(ymin = ymin, ymax = ymax, fill = color2),
      xmin = -Inf, xmax = Inf,
      alpha = 0.15,
      inherit.aes = FALSE
    ) +
    geom_line() +
    facet_wrap(. ~ param, scale = "free_y") +
    scale_fill_identity() +
    theme_bw() +
    ggtitle("Balance Sheet Metrics")
  
  return(plot)
}

cashflow_metric <- function(df, income) {
  df_c <- df |>
    mutate(operatingCashFlow = netIncome + depreciationDepletionAndAmortization,
           capitalExpenditureToNetEarning = capitalExpenditures / netIncome,
           stockBuybackToNetEarning = case_when(
             proceedsFromRepurchaseOfEquity < 0 ~ -proceedsFromRepurchaseOfEquity / income$netIncome,
             income$netIncome < 0 & proceedsFromRepurchaseOfEquity > 0 ~ NA_real_,
             income$netIncome < 0 & proceedsFromRepurchaseOfEquity < 0 ~-proceedsFromRepurchaseOfEquity / income$netIncome))
  
  hline <- tribble(
    ~param, ~hline_value, ~color2,
    "operatingCashFlow", 0, "red", 
    "capitalExpenditureToNetEarning", 0.5, "blue",
    "stockBuybackToNetEarning", 0.5, "red"
  ) |>
    mutate(
      param = factor(param),
      ymin = ifelse(color2 == "red", -Inf, hline_value),
      ymax = ifelse(color2 == "red", hline_value, Inf)
    )
  
  columns <- hline$param
  
  plot <- df_c |>
    select(fiscalDateEnding, operatingCashFlow, capitalExpenditureToNetEarning, stockBuybackToNetEarning) |>
    pivot_longer(cols = c(operatingCashFlow, capitalExpenditureToNetEarning, stockBuybackToNetEarning), 
                 names_to = "param", values_to = "values") |>
    mutate(param = factor(param, levels = columns)) |> 
    ggplot(aes(x = fiscalDateEnding, y = values)) +
    geom_rect(
      data = hline,
      aes(ymin = ymin, ymax = ymax, fill = color2),
      xmin = -Inf, xmax = Inf,
      alpha = 0.15,
      inherit.aes = FALSE
    ) +
    geom_line() +
    facet_wrap(. ~ param, scale = "free_y") +
    scale_fill_identity() +
    theme_bw() +
    ggtitle("Cash Flow Metrics")
  
  return(plot)
}

show_all <- function(name) {
  income <- get_data("INCOME_STATEMENT",name)
  share <- get_data("EARNINGS",name,share=T)
  balance <- get_data("BALANCE_SHEET",name)
  cashflow <- get_data("CASH_FLOW",name)
  plot1 <- income_metric(income) 
  plot2 <- share_metric(share)
  plot3 <- balance_metric(balance, income)
  plot4 <- cashflow_metric(cashflow, income)
  stacked_plot <- ggarrange(plot1,plot3,nrow=2)
  squished_plot <- ggarrange(plot4,plot2,ncol=2, widths = c(2,1))
  combineplot <- ggarrange(plotlist = list(stacked_plot,squished_plot), nrow = 2, heights = c(2,1))
  return(combineplot)
}

Let’s Look At Another Example

As a heuristics, we’ve made our dataviz where red is like lava, you want to stay above it. Blue is like the sky, you want to stay below it. The sweet zone is in between. When there is a loess curve, we are trying to see if the share is consistently increasing.

Apple

show_all("AAPL")

Looking the income metrics, Apple seems to be doing pretty good. Gross profit margin, net earning margin, operating income margin are all above thresholds. Depreciation to gross profit ratio is appropriately below the threshold, which means there is no recent purchase of property, plants, machines, etc where these are the ones that will add to depreciation. There is also a low to no interest expense, which is good! No debt? Next we move on to balance sheet metrics, current ratio is not great, liquidity is not great, net receivable is ?OK, I guess it make sense, if your products are popular and provide some sort of finance option, you might have some receivables. Short term debt is much lower than long term debt, which is good. The ROSE is smelling pretty good there too! In terms of cash flow metrics, it’s looking really good, it has high cash flowing in, low CapEx, and consistently buying its own share from 2020 onwards.

Microsoft

show_all("MSFT")

As for Microsoft, now that we’re a bit more familiar with the sections, let’s just string our read instead of separating them. Good consistent profit for all profit, operating and net income. Interestingly, depreciation is high along with CapEx, did they buy property, plants or machine? Ah, they built data center? Interest expense is low, which is good. Good liquidity factor with current ratio above threshold, net receivables is OK same as Apple (a competitor), short/long debt is good as well. ROSE is close to the threshold (worth watching). Cash is flowing, not much stock buyback, but overall EPS is consistently increasing. With its investment in data centers for Azure, next few years we should be seeing a persistent depreciation to gross profit ratio, along with CapEx, right? I read somewhere where you can’t place all depreciation in a single year.

NVIDIA

show_all("NVDA")

Now let’s take a look at NVIDIA. Wow, gross profit, operating income, and net earning margin all really good, better than Apple and Microsoft! Well, make sense with all these data centers, they need GPUs from NVIDIA. Look at that depreciation and CapEx, barely any! That’s great, less to maintain and the current plants they have are adequate to supply the demand. Low interest expense too, not much debt interest. Liquidity is great with high current ratio! Net receivables is good too, these big companies are paying NVIDIA back! ROSE is also very fragrant! Cash is flowing through the roof. Not much stock buyback. Also exponential EPS! You know, Warren did say that be mindful of companies with R&D cost.

Intel

show_all("INTC")

Wow, first glance with the profits, we’re seeing lava red. We then see high rise in depreciation and CapEx. Did they buy more property, plants, or machines? Intel announces plans for an investment of more than $28 billion for two new chip factories in Licking County, Ohio. And maybe Germany and Arizona too? Interesting. Notice that the interest expense towards the end of 2025 was in 2 digit but negative? It’s because the the operating income was negative and a high interest expense. I wonder if I should make this absolute as opposed to keeping the negative. Anyway, this means that not only was Intel in the red but also paying quite a bit of interest from the debt, I assume. The current ratio looks pretty good, same goes with the net receivables. Short term debt is also much lower compared to long term. ROSE not so good. As for the cash flow, when we add depreciation and demortization back to net income, it brought the red back up. Notice how the EPS dropped significantly lately, and notice that the buyback of share is missing? We code it where if it’s a positive value of buyback, which means it’s issuing stock instead of repurchasing, hence NA. Wow, what do you think? Will they be able to turn this around? Do data centers use Intel chips or AMD’s threadripper?

AMD

show_all("AMD")

Interesting income statement metric results. Good gross profit margin, but the net earning is barely. Interesting depreciation trend, what happened in 2024, where CapEx didn’t budge? Great liquidity, pretty high net receivables. Wow, very high short-to-long term debt in ?2021-2022. ROSE is essentially 0. With cash flowing good towards 2025, and also high buybacks 2022 ish. Finally, uptrending EPS! This is an odd one. Pasting my observation onto Claude and wow, these findings are due to Xilinx acquisition. That makes sense! Didn’t build new plants, old property, plants, machines have already been amortized, hence low depreciation. The acquisition is funded by debt, hence high short-to-long term debt ratio. Very interesting, indeed! So, it could be true that data centers prefer AMD over Intel given the good profit?

Google

show_all("GOOG")

Very good profit of all 3 metrics. After 2020, depreciation downtrended which is good and so did CapEx. Not much debt. Downtrending current ratio but still good. Very good net receivables! ROSE is coming up. Great cash flow. There is buyback of equity. Rising EPS. All great signs for Google! Being an all-service (no hardware?) company, this is quite good and healthy!

3M

show_all("MMM")

Now on to 3M, not a tech company, but let’s see if our workflow of looking at financial statement will help us tell a story. Profit is not as great as other tech companies we looked at, though gross profit is above threshold, but the net earnings is in the lava. Notice both net earnings and operating margins dipped quite low in 2024? And then a pretty high interest expense to operating income ratio in 2025? I wonder if because of losing money, they borrowed long term, hence short-to-long term debt is not too high? Current ratio is acceptable, but why is net receivables so high for 3M? Clients were not able to pay 3M? ROSE is pretty good in 2024 and 2025 even when profits weren’t. But why? Cash flow with high variance last 3 years. There is stock buyback in 2025. And even though EPS downtrended past 2 years but still remained high. This is a very interesting one as well. 3M should be a company where there’s competitive advantage because they make lots of daily usables which doesn’t need a whole lot of R&D. But why the anomaly in 2024 requiring debt? What do you think?

Disney

show_all("DIS")

Alright, what about Disney? Profit doesn’t seem as good as I expected. Both gross profit and net earning margins were below threshold, except for operating income margin. Interesting, in 2021, there is increase in depreciation, CapEx, and interest expense with low short-to-long debt ratio. It almost looked like they borrowed some money to purchase something new? Did they rebuild a park or something? In terms of liquidity, it’s on the lava zone in 2025 and seemed like it downtrended for the past 3 years as well. Net receivables is OK? Though I would think it should be lower? ROSE is red. Cash flow looks good, this is interesting because whatever that caused the depreciation when added back now they have cash. Also interesting to note that from 2020 to 2024 they were issuing stock rather than buying back. Their EPS appear to be quite volatile for the pat 3 years.

Boeing

show_all("BA")

Last but not least, let’s look at Boeing. Three profit margin metrics in the lava zone. Spike of interest expense in 2021 along with depreciation. Something physical was bought with borrowed money that’s paid long term. Current ratio is good. Essentially no net receivables. Consistent floating ROSE on the red sea. Operating cash flow in the red for a few years then positive in 2025. CapEx spiked in 2025, ?what was spent in capital where it doesn’t really depreciate? Buyback stock is interesting, missing several past few years, essentially issuing stocks when we have NA data. EPS downtrend to the negatives.

Opportunities For Improvement

• We need to recheck growth rate again, not sure if the code is correct
• perhaps build this into an R package, and develop further
• Learn about valuation metrics such as P/E, EV/EBITDA, P/B, P/S
• Not really sure if short-to-long term debt ratio of 1 is a good threshold, seems too high.
• include actual numbers on geom_label with ggrepel and reduce font size

Lessons learnt

• Depreciation and CapEx seem to have strong correlation, which makes sense.
• learnt Alpha Vantage API, quite straight forward
• learnt to look at financial statement and its metrics

If you like this article:

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

Following the post on exact Shapley values for time series explainability, this post illustrates an example of how to use sensitivity analysis to explain time-series forecasts, based on the ahead::dynrmf model and external regressors. What is sensitivity analysis in this context? It’s about evaluating the impact of changes in the external regressors on the time-series forecast.

The post uses the ahead::dynrmf_sensi function to compute the sensitivities, and the ahead::plot_dynrmf_sensitivity function to plot the results.

First, install the package:

devtools::install_github("Techtonique/ahead")

Then, run the following code:

# devtools::install_github("Techtonique/ahead")
# install.packages(c("fpp2", "e1071", "patchwork"))

library(ahead)
library(fpp2)
library(patchwork)
library(e1071)

#' # Example 1: US Consumption vs Income
sensitivity_results_auto <- ahead::dynrmf_sensi(
y = fpp2::uschange[, "Consumption"],
xreg = fpp2::uschange[, "Income"],
h = 10
)

plot1 <- ahead::plot_dynrmf_sensitivity(sensitivity_results_auto, 
                           title = "Sensitivity of Consumption to Income (Ridge)",
                           y_label = "Effect (ΔConsumption / ΔIncome)")

#' # Example 1: US Consumption vs Income
sensitivity_results_auto_svm <- ahead::dynrmf_sensi(
  y = fpp2::uschange[, "Consumption"],
  xreg = fpp2::uschange[, "Income"],
  h = 10, 
  fit_func = e1071::svm # additional parameter passed to ahead::dynrmf
)

plot2 <- ahead::plot_dynrmf_sensitivity(sensitivity_results_auto_svm, 
                                        title = "Sensitivity of Consumption to Income (SVM)",
                                        y_label = "Effect (ΔConsumption / ΔIncome)")

 
# Example 2: TV Advertising vs Insurance Quotes
sensitivity_results_tv <- ahead::dynrmf_sensi(
 y = fpp2::insurance[, "Quotes"],
   xreg = fpp2::insurance[, "TV.advert"],
   h = 8
 )

plot3 <- ahead::plot_dynrmf_sensitivity(sensitivity_results_tv,
                           title = "Sensitivity of Insurance Quotes to TV Advertising (Ridge)",
                           y_label = "Effect (ΔQuotes / ΔTV.advert)")

sensitivity_results_tv_svm <- ahead::dynrmf_sensi(
  y = fpp2::insurance[, "Quotes"],
  xreg = fpp2::insurance[, "TV.advert"],
  h = 8, 
  fit_func = e1071::svm # additional parameter passed to ahead::dynrmf
)

plot4 <- ahead::plot_dynrmf_sensitivity(sensitivity_results_tv_svm,
                                        title = "Sensitivity of Insurance Quotes to TV Advertising (SVM)",
                                        y_label = "Effect (ΔQuotes / ΔTV.advert)")

(plot1+plot2)

(plot3+plot4)

LLM-enhanced momentum investing combines traditional momentum signals with real-time news interpretation by large language models (LLMs). The idea is straightforward: stocks with strong past returns are candidates for momentum portfolios, but their inclusion and weight are refined by LLM-generated sentiment scores derived from firm-specific news. This hybrid approach improves risk-adjusted returns (Sharpe, Sortino) and is particularly effective in concentrated, high-conviction portfolios.

Key Parameters

1. Lookback Window (k)

• Definition: The number of past days of news considered for sentiment analysis.
• Role: Determines how much recent information the LLM uses to judge momentum continuation.
• Example: If k = 5, the model analyzes the last 5 business days of headlines and summaries for each stock.

2. Forecast Horizon (l)

• Definition: The period over which momentum continuation is predicted.
• Role: Sets the “target” for the LLM’s forecast — how far into the future the model should judge momentum persistence.
• Example: If l = 5, the LLM predicts whether momentum will continue for the next 5 trading days.
• Connection to Rebalancing: The forecast horizon typically aligns with the rebalancing cycle. For weekly rebalancing, the horizon is 5 days; for monthly rebalancing, it’s ~21 days.

3. Portfolio Size (m)

• Definition: The number of stocks selected after LLM scoring.
• Role: Controls how concentrated or diversified the portfolio is.
• Example: From the top 20 YTD performers, you might select the top 10 after sentiment scoring.

4. Rebalancing Frequency (T)

• Definition: How often the portfolio is updated with new signals.
• Role: Sets the rhythm of portfolio refresh — weekly, monthly, or quarterly.
• Example: Weekly rebalancing means recalculating momentum and sentiment scores every 5 trading days.

Concept

The strategy begins with a classic momentum screen: select the top 20 S&P 500 companies by year-to-date (YTD) performance. Instead of stopping there, the approach integrates large language model (LLM) sentiment analysis of firm-specific news. By analyzing the last 5 business days of headlines and summaries, the LLM produces a score indicating whether momentum is likely to continue.

These scores are then used to re-weight the portfolio, tilting allocations toward companies with stronger news sentiment. Finally, the portfolio is narrowed to the top 10 conviction stocks.

Selected Parameters

• Lookback Window: 5 days of firm-specific news.
• Rebalancing Frequency: Weekly updates of the portfolio.
• Forecast Horizon: 5 trading days (aligned with the rebalancing cycle).

This setup ensures that the LLM is asked to judge whether momentum will persist until the next rebalance, making the signals both short-term and actionable.

What the Code Does Step by Step

1. Fetching the Data
   • The R script first pulls all S&P 500 tickers.
   • It calculates YTD returns for each stock.
   • The top 20 stocks by performance are selected as momentum candidates.
2. News Sentiment Analysis with LLM
   • For each of these 20 stocks, the code queries Bing News for recent headlines with Azure AI Services.
   • The last 5 business days of news are collected.
   • These headlines and summaries are sent to a Microsoft Foundry-hosted LLM.
   • The LLM outputs a score (0–1) indicating whether sentiment supports momentum continuation or signals reversal.
3. Portfolio Tilting
   • LLM scores are normalized to [-1, +1].
   • Baseline equal weights are tilted according to these scores.
   • The top 10 stocks by adjusted weight form the final portfolio

4. Visualization

• A styled table is created using the gt package.
• Adjusted weights are color-coded (red–green gradient).
• The final portfolio is saved as an image (top10.png).

Strategic Insight

This approach mirrors the methodology in the Swiss Finance Institute paper:

• Momentum ranking provides the baseline.
• LLM sentiment scoring refines stock selection and weighting.
• Portfolio tilting integrates qualitative news signals into quantitative allocation.

library(httr)
library(jsonlite)
library(tidyquant)
library(tidyverse)
library(lubridate)
library(gt)
library(gtExtras)
library(scales)
library(showtext)
library(webshot2)

# 1. Environment & Auth Setup
sysfonts::font_add_google("Roboto Slab", "roboto_slab")
showtext_auto()

# Azure & Bing Credentials
bing_key           <- "<your-bing-key>"
bing_endpoint      <- "<your-bing-endpoint>" 
azure_llm_key      <- "<your-llm-key>"
azure_llm_endpoint <- "<your-llm-endpoint>"

# R Part (first): S&P 500 Screening with Monthly Returns
sp500_tickers <- 
  tq_index("SP500") %>% 
  select(symbol, company)

# Calculate YTD change
momentum_df <- 
  sp500_tickers %>%
  tq_get(get = "stock.prices", from = floor_date(today(), "year")) %>%
  group_by(symbol) %>%
  arrange(date) %>% # Ensure chronological order for first/last functions
  summarize(
    total_return = (last(adjusted) / first(adjusted)) - 1, 
    .groups = "drop"
  ) %>%
  inner_join(sp500_tickers, by = "symbol") %>%
  slice_max(total_return, n = 20) %>%
  select(symbol, company)

# R Part (second): News Search and LLM Analysis
analyze_momentum_continuation <- function(ticker, company_name) {
  
  # Construct the Bing News Search query: ticker + ' stock'
  query_str <- paste0(ticker, " stock")
  news_url <- paste0(bing_endpoint, "v7.0/news/search")
  
  # Call Bing News Search API
  news_res <- GET(news_url, add_headers(`Ocp-Apim-Subscription-Key` = bing_key), 
                  query = list(q = query_str, count = 5, freshness = "Day"))
  
  # Short pause to throttle calls and avoid rate limit errors
  Sys.sleep(1)
  
  news_text <- ""
  if (status_code(news_res) == 200) {
    content <- fromJSON(content(news_res, "text", encoding = "UTF-8"))
    if (length(content$value) > 0) {
      news_text <- paste(content$value$name, content$value$description, collapse = " | ")
    }
  }
  
  # Construct the LLM payload for Azure AI Foundry
  prompt_payload <- list(
    messages = list(
      list(role = "system", content = "You are an LLM Enhanced Momentum Investing Agent."),
      list(role = "user", content = paste0(
        "Headlines + Summaries for ", company_name, " (", ticker, "): ", news_text,
        "\nPerform sentiment analysis based on the last 5 days of news (lookback=5, horizon=5). ",
        "Infer whether sentiment supports momentum continuation or signals reversal. ",
        "Return a JSON object with two fields: 'subsector' (string) and 'llm_score' (string: probability 0-1)."))
    ),
    temperature = 0.1
  )
  
  llm_res <- POST(url = azure_llm_endpoint, 
                  add_headers(`api-key` = azure_llm_key, `Content-Type` = "application/json"),
                  body = prompt_payload, encode = "json")
  
  if (status_code(llm_res) == 200) {
    llm_out <- fromJSON(content(llm_res, "text", encoding = "UTF-8"))
    # Parse JSON response without using regex
    llm_json_data <- fromJSON(llm_out$choices$message$content)
    return(as.data.frame(llm_json_data))
  } else {
    return(data.frame(subsector = "N/A", llm_score = "0.5"))
  }
}

# Execute Analysis: Merge top 20 tickers with LLM scores
news_scores_df <- momentum_df %>%
  mutate(analysis = map2(symbol, company, analyze_momentum_continuation)) %>%
  unnest(analysis) %>%
  mutate(llm_score = as.numeric(llm_score))

# R Part (final): Portfolio Tilting and Visualization
# Normalize scores to [-1, +1] and tilt weights
tilted_portfolio <- news_scores_df %>%
  mutate(
    norm_score = rescale(llm_score, to = c(-1, 1), from = c(0, 1)),
    base_weight = 1 / n(),
    adj_weight = base_weight * (1 + norm_score)
  ) %>%
  mutate(adj_weight = adj_weight / sum(adj_weight)) %>%
  slice_max(adj_weight, n = 10)

# Create gt visualization using original column names
final_table <- tilted_portfolio %>%
  select(company, subsector, adj_weight) %>%
  gt() %>%
  tab_header(title = "Top 10 Tilted S&P 500 Momentum Stocks") %>%
  # Use cols_label for human-readable labels without renaming underlying columns
  cols_label(
    company = "Company", 
    subsector = "Subsector", 
    adj_weight = "Adjusted Weight "
  ) %>%
  # Apply color intensity with scales::col_numeric
  data_color(
    columns = adj_weight, 
    colors = col_numeric(palette = c("red", "green"), domain = NULL)
  ) %>%
  fmt_percent(
    columns = contains("adj_weight"), 
    decimals = 2,
    locale = "en" 
  ) %>% 
  cols_align(align = "center") %>%
  opt_table_font(font = google_font("Roboto Slab"))

# Save the visualization as top10.png using webshot
gtsave(final_table, "top10.png")

Final Observation

Looking at the resulting portfolio, one striking feature is the dominance of energy and petroleum companies. Firms such as ConocoPhillips, EOG Resources, ExxonMobil, Occidental Petroleum, Marathon Petroleum, Valero, Phillips 66, Chevron, and Baker Hughes all appear prominently.

This heavy tilt toward energy is not random—it reflects how geopolitical tensions (Iran and U.S.–Israel war dynamics) have amplified the importance of oil and gas in global markets. News sentiment around these companies has been strongly supportive of continued momentum, pushing them into the top allocation slots.

Over the course of my career as a Data Scientist, I’ve worked on projects ranging from simple code reviews, to large application builds. For the most part, I have used R to do this.

If you’re getting into coding or data science, one question you’re probably asking yourself is “Which language should I learn?”

This blog aims to show you why R might be a good decision.

R was built for data (not just programming)

Unlike general purpose languages (such as Python), R was designed specifically for statistics and data analysis.

That means:

• Built in statistical tools
• Powerful visualisation capabilities
• Research level methods available immediately

With packages like the tidyverse, you can clean, analyse, and visualise data with surprisingly little code.

High demand in analytics, research, and healthcare

R is especially popular in many sectors such as:

• Healthcare & biostats
• Academic research
• Government departments
• Finance & risk modeling
• Pharmaceutical companies

Here are some examples of R in production use:

• The {bbplot} R package. Yes, the BBC use R to create graphics for their website!
• Health and wellbeing profiling app for the NHS
• During the Covid-19 pandemic, the financial times had a stats tracker in which the graphs were built with R.

Knowing some R will give you a competitive edge if you’re looking at working within these sectors.

Open source with the backing of Posit

R is open source. This means that:

• It’s free, and always will be!
• Anyone can view the source code the makes up R, there are.
• Each R package (a folder containing code) has to live on GitHub.com, for everyone to see.
• It has a large community of contributors. There are great forums to get help such as Stack Overflow, Posit Community and the R weekly newsletter and tonnes more.
• There are thousands more available functionalities compared to paid softwares such as SPSS, SAS or Excel.

Posit, who maintain the free to use RStudio and Positron IDEs (integrated development environment), have many full time staff working solely on maintaining and creating new functionality within R. This means we get:

• Defined accountability
• Predictable release cycles
• Bugs can be solved quicker

Incredible data visualisation possibilities

Being able to communicate your findings with stakeholders is very important in data science, and one of R’s biggest strengths is visualisation and reporting.

With the {ggplot2} package, you can create publication ready charts with very little code. The R Graph Gallery has some amazing examples of what is possible with {ggplot2}.

With the {quarto} and {shiny} packages, you are able to build reproducible reports and interactive dashboards. All this without needing to know any HTML, CSS or JavaScript.

Beginner friendly learning curve

This is very much my own opinion. Compared to other languages, I think R is fairly intuitive and feels rewarding much earlier on in the journey. It also has (in my opinion), the most beginner friendly programme to code in, called RStudio.

Most people attend only two days worth of training with Jumping Rivers, and say they feel ready to start tackling their own data problems.

So… is R worth learning in 2026?

I think so. If you want pure software engineering or large-scale production systems, you may need Python. But for becoming a strong data thinker, and giving you an edge in your analysis, R is one of the best starting points.

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