Reposted from the original at https://blog.stephenturner.us/p/ctrlvee-extract-external-r-code-insert-inline-positron-rstudio-addin. 

Ever find yourself looking through a pkgdown page or a Quarto book, copying and pasting code chunks from your browser into your IDE? I do, and it’s a minor annoyance.1

My friend and colleague VP Nagraj published a new R package called ctrlvee that makes this a lot easier.

• CRAN: https://cran.r-project.org/package=ctrlvee

• GitHub: https://github.com/vpnagraj/ctrlvee

It does one thing. Put your cursor anywhere in an R script in Positron or RStudio, call the add-in, provide a URL, and a few milliseconds later you’ll have all the code from that page in your editor, separated by chunk boundaries (along with some metadata and a note to check the license!).

The package README provides a demonstration using the “Data Validation and QA” chapter of my Data Science Team Training book (dstt.stephenturner.us).

1. Install the package: install.packages("ctrlvee")

2. Run the add-in. In Positron you’ll open the command palette, search for Run RStudio Addin, then extract external R code and insert inline. You’ll get a modal asking you for a URL.

3. Paste one in. E.g., https://dstt.stephenturner.us/validation.html

4. The R code from the website appears in your editor

Here’s a demo.

Here’s what the extracted/inserted code looks like, from this source.

# -----------------------------------------------------------------
# Chunks fetched by ctrlvee from: https://dstt.stephenturner.us/validation.html
# Strategy: Rendered HTML page
# Date: 2026-05-16 05:14:44
# Chunks: 8
# NOTE: Check the source license before reusing this code.
# -----------------------------------------------------------------

flu <- data.frame(
    week = c(1, 2, 3, 4, 4),
    county = c("Fairfax", "Arlington", NA, "Loudoun", "Loudoun"),
    disease = c("Flu", "Flu", "Flu", "Flu", "Flu"),
    cases = c(23, 41, 18, -5, 12),
    rate = c(2.1, 3.8, 1.6, NA, 1.1)
)

flu

# ---- chunk boundary ----

if (any(flu$cases < 0, na.rm = TRUE)) {
    stop("Negative case counts detected. Inspect raw data before proceeding.")
}

# ---- chunk boundary ----

stopifnot(
    "Negative case counts" = all(flu$cases >= 0, na.rm = TRUE),
    "Missing county values" = !anyNA(flu$county),
    "Duplicate records" = !anyDuplicated(flu[, c("week", "county")])
)

# ---- chunk boundary ----

install.packages("pointblank")

# ---- chunk boundary ----

library(pointblank)

agent <- create_agent(tbl = flu, label = "Weekly flu surveillance") |>
    col_vals_gte(
        columns = cases,
        value = 0,
        label = "Case counts must be non-negative"
    ) |>
    col_vals_not_null(
        columns = c(week, county),
        label = "Week and county cannot be missing"
    ) |>
    rows_distinct(
        columns = c(week, county),
        label = "No duplicate week/county records"
    ) |>
    interrogate()

agent

# ---- chunk boundary ----

create_agent(tbl = flu, label = "Weekly flu surveillance — extended") |>
    col_is_numeric(
        columns = c(cases, rate),
        label = "Case count and rate must be numeric"
    ) |>
    col_vals_in_set(
        columns = disease,
        set = c("Flu", "COVID-19", "RSV"),
        label = "Disease must be from the approved list"
    ) |>
    col_vals_between(
        columns = week,
        left = 1,
        right = 52,
        label = "Week must be between 1 and 52"
    ) |>
    col_vals_gte(
        columns = rate,
        value = 0,
        na_pass = TRUE,
        label = "Rate must be non-negative (NAs allowed)"
    ) |>
    interrogate()

# ---- chunk boundary ----

if (!all_passed(agent)) {
    stop("Data validation failed. Review the agent report before proceeding.")
}

# ---- chunk boundary ----

library(readr)
library(pointblank)

flu <- read_csv("data/flu-2024.csv")

# Validate immediately after reading
agent <- create_agent(tbl = flu, label = "flu-2024 validation") |>
    col_vals_gte(columns = cases, value = 0, label = "No negative counts") |>
    col_vals_not_null(columns = c(week, county), label = "No missing keys") |>
    rows_distinct(columns = c(week, county), label = "No duplicate records") |>
    interrogate()

if (!all_passed(agent)) {
    stop("Validation failed — see agent report above.")
}

If you’ve read any of my past posts you know I like to program in several different languages, some of which I like more than others. Sometimes a problem calls for a particular language to be used, and with that comes adjusting one’s brain to thinking in that language and using the appropriate idioms to leverage that language’s features. But what if I don’t want to?

I don’t want to

The line between R and Python has been heavily blurred the last few years, particularly with {reticulate} enabling us to use Python within R code, RStudio rebranding as Posit and taking on a strong Python development effort, releasing Positron as a multi-language IDE, and Quarto being a multi-language rethink of Rmarkdown.

I occasionally need to use Python directly – an SDK wrapping an API exists and I don’t particularly want to spend a lot of time writing my own R version, especially before I know what I want to get out of the endpoints. At this point I tend to bump up against my muscle-memory from R and try to use functions I’m familiar with from R, but which don’t actually exist in Python. Now, that might sometimes be because the pattern I’m trying to encode simply has a different name in Python; instead of an sapply(x, f)

sapply(c(2, 3, 4, 5), \(x) x ^ 2)
## [1]  4  9 16 25

I should reach for map, in which case I am reminded that this produces a lazy iterator that doesn’t show me the results

map(lambda x: x ** 2, [2, 3, 4, 5])
## <map object at 0x10d7fbee0>

and so I need to wrap it into a list to get the values out

list(map(lambda x: x ** 2, [2, 3, 4, 5]))
## [4, 9, 16, 25]

Or, I could use a list comprehension which isn’t lazy

[v ** 2 for v in [2, 3, 4, 5]]
## [4, 9, 16, 25]

That’s the idiom that I should be reaching for. Sure.

Other times there’s a package I need to use and a slightly different way of approaching the problem. In R I love the table() function for getting histogram-like counts of the unique values of a vector

table(c("b", "a", "c", "a", "b", "a"))
## 
## a b c 
## 3 2 1

which in Python looks like

from collections import Counter

sorted(Counter(["b", "a", "c", "a", "b", "a"]).items())
## [('a', 3), ('b', 2), ('c', 1)]

Probably Pythonistas remember that idiom and the package to import and the .items() extractor and the fact that they maybe want to sort the result. But I kept coming back to a question I ask myself: what if I don’t want to? Why is there not a function that wraps this idiom? If there was, why not just call it “table”? Admittedly, it’s far from the catchiest, most memorable, or most useful name, but it’s immediately recognisable to an R user (ditto for “sapply”).

One approach I considered here was to just call R from Python. That can be done, but I doubt I or anyone else wants to deal with that every time we want to iterate over a list. There’s a package on the Python package index which seems to support this nicely: https://pypi.org/project/r-functions/ but it’s wrappers around individual R files, via RScript. I’m thinking more along the lines of ‘native Python with an R interface’.

Python is an object-oriented language, but it has functions, so why not make one

from collections import Counter

def table(x):
    return dict(sorted(Counter(x).items()))

table(["b", "a", "c", "a", "b", "a"])
## {'a': 3, 'b': 2, 'c': 1}
def sapply(x, func):
    return [func(v) for v in x]
  
sapply([2, 3, 4, 5], lambda x: x ** 2)
## [4, 9, 16, 25]

and have a nicer function interface to apply these idioms? I thought about this a bit longer, and realised there’s lots of functions I use in R that I wish I could use in Python. An idiom for finding the index of elements of a ‘vector’ (list in Python) which are true (TRUE in R, True in Python) is

[i for i, v in enumerate(x) if v]

but I just want to call which(x)

which(c(FALSE, FALSE, TRUE, FALSE , TRUE))
## [1] 3 5

so why not define this

def which(x):
    return [i for i, v in enumerate(x) if v]
  
which([False, False, True, False, True])
## [2, 4]

(remembering that Python is 0-indexed).

How far could one take this? Quite a long way!

I thought more about what differences would need to be accounted for, and one that immediately came to mind was that R is vectorised. If I was to recreate R’s character counting function nchar(s) as essentially len(s), I’d need to consider whether I wanted it to work on a single string or a ‘vector’ of strings

In R:

nchar(c("these", "all", "have", "different", "lengths"))
## [1] 5 3 4 9 7

But in Python, len() expects a single value, so it calculates the length of the list

len(["these", "all", "have", "different", "lengths"])
## 5

The ‘proper’ way to do it is to map over the list

[len(s) for s in ["these", "all", "have", "different", "lengths"]]
## [5, 3, 4, 9, 7]

but again, why do I need to use an idiom for this? What if I just made a decorator to change a regular function to a vectorised one by applying this list comprehension internally when it’s passed a list (or a tuple), and which otherwise just evaluates the function with the argument?

import functools

def make_vec(func):
    @functools.wraps(func)
    def wrapper(*args, **kwargs):
        if isinstance(args[0], (list, tuple)):
            return [func(xi, *args[1:], **kwargs) for xi in args[0]]
        return func(*args, **kwargs)
    return wrapper

@make_vec
def my_len(s):
    return len(s)

my_len(["these", "all", "have", "different", "lengths"])
## [5, 3, 4, 9, 7]

and I could name it… “nchar”!

The other use-case that came to mind was Elio venting (and referencing a post to which I also wrote a sort of response) that they needed to list the files in the current directory

Post by @eliocamp@mastodon.social

View on Mastodon

with the idiom

import os

[os.path.join(path, f) for f in os.listdir(path)]

The supplied suggestions included

from pathlib import Path

list(Path(path).iterdir())

(just rolls off the tongue, doesn’t it?) which returns a list of PosixPath() objects and is hardly easy to parse visually.

So, why not have a function?!?

import os

def list_files(path):
    return [os.path.join(path, f) for f in os.listdir(path)]

path = "path/to/files"

list_files(path)
## ['path/to/files/file1.txt', 'path/to/files/file2.txt', 'path/to/files/file3.csv']

I would have liked to call this list.files() but, since Python strictly uses the dot for method calling, it can’t be that.

This then raises the question of “should I support the arguments already in the R functions?” In this case, should it support a recursive argument? Yes, that adds complexity, but it’s surely do-able. At this point I reached for some AI assistance and had Claude help me to implement as many functions as we could think of, supporting as many common arguments as possible. This involved extending the decorator to support vectorising other arguments (which also need to be careful about dots).

On testing it out, it looked like we had something viable.

One last piece I wanted to support, though: the which() example above extracts the elements of a logical vector which are True, but in order to build that vector in the first place, I would naturally leverage R’s vectorisation as an array language. The two steps involved here are to first compute the comparison resulting in a logical vector, then to use which() to identify the indices of those which are true

which(c("c", "b", "a", "c", "a", "b") == "a")
## [1] 3 5

The vectorisation decorator above doesn’t help here, because it’s at the point of == that we want to vectorise

['c', 'b', 'a', 'c', 'a', 'b'] == 'a'
## False

This is False because the character 'a' is not equal to the given list.

The appropriate idiom is once again to use a list comprehension

which(x == 'a' for x in ['c', 'b', 'a', 'c', 'a', 'b'])
## [2, 4]

The solution I’m fond of is to create a new ‘Vec’ class which wraps binary operators with a list comprehension, again abstracting away this detail. This means implementing __eq__, __add__, __and__ and lots of other binary operations, but with that, and a wrapper to create such an object, the comparison operators can be vectorised

vals = vec(['c', 'b', 'a', 'c', 'a', 'b'])
which(vals == 'a')
## [2, 4]

Not pristine, but quite clean, if you ask me.

With all these pieces in place, adding implementations for common base R functions including most arguments and a way to vectorise lists, I wrapped everything up into a Python package (my first) to learn how to do it.

The workflow isn’t particularly painful, with my biggest complication being different versions of Python supporting different requirements in pyproject.toml, and so some GitHub Actions are failing because of that.

As part of building out the implementations I had Claude add tests for each of the functions with some expected values – if I do want to improve some of the idioms internally, I want to ensure I don’t change the values produced. That works for having any testing at all, but how can I be sure that I’m reproducing what I would get if I was working in R? One option was to just run all of the test functions by hand and confirm that the values look similar enough, accounting for list vs vector and 0 vs 1 indexing. Instead, Claude managed to write an adaptor for pytest which does the realignment of e.g. list_files to list.files (and similarly for arguments), realigns the indexing where needed, and runs all existing tests directly in R via rpy2 (skipping over some for which I don’t have tests yet). I’m disabling automated testing of this because I suspect it could get flaky dealing with both R and Python on GitHub Actions, but I can confirm that all the current tests pass.

I wanted to have a documentation website similar to what we have via {pkgdown} and came across quartodoc which is what the Python version of {pins} uses. Getting that to work required downgrading a specific Python dependency, but was otherwise painless.

I have a working package locally – how do I share it? This seemed like the perfect opportunity to learn what the release process looks like for Python. I have a handful of packages on CRAN and one on Bioconductor, and the process there is far from frictionless, with the side-effect that there’s some trust you can place on the interoperability of packages and minimal (automated) code checking. While Python is more ‘wild west’ in terms of what can be uploaded, it’s really nice to see that they do have an entirely separate test server where you can upload your package and see how it looks. I’m reminded of the quote

Everybody has a testing environment. Some people are lucky enough to have a totally separate environment to run production in.

Given that it’s not currently possible to run 100% of the CRAN checks locally (and even some that you can give a different result to what’s on their systems) this does make me a little jealous. I wonder whether the decrease in load from rejecting failing submissions would offset supporting a test submission server.

All went well pushing to the test server (via an authentication key) and I managed to build up the courage to push to the production instance… it’s live!

rfuns logo – R functions in Python… are fun

and the documentation site isn’t too bad, either (in my opinion).

This means that you can now run

uv add rfuns

(or the equivalent in whatever virtual environment management configuration you’re using, e.g. pip install rfuns) and start using some R functions directly in Python!

Depending on how you like to manage your imports, you can import everything

from rfuns import *

which([False, False, True, False, True])
## [2, 4]

or, if you prefer to namespace

import rfuns as r

r.which([False, False, True, False, True])
## [2, 4]

The list of functions currently imported, grouped into sections is:

nchar(["these", "all", "have", "different", "lengths"])
## [5, 3, 4, 9, 7]
grepl("ar", ["frog", "carpet", "basket", "dart"])
## [False, True, False, True]
sqrt([36, 81, 9])
## [6.0, 9.0, 3.0]

seq(5)
## [0, 1, 2, 3, 4]
seq(from_=0, to=10, by=2)
## [0, 2, 4, 6, 8, 10]

setdiff([5, 2, 4, 1], [2, 1])
## [5, 4]

set([5, 2, 4, 1]) -  set([2, 1])
## {4, 5}

## ─ Session info ───────────────────────────────────────────────────────────────
##  setting  value
##  version  R version 4.5.3 (2026-03-11)
##  os       macOS Tahoe 26.3.1
##  system   aarch64, darwin20
##  ui       X11
##  language (EN)
##  collate  en_US.UTF-8
##  ctype    en_US.UTF-8
##  tz       Australia/Adelaide
##  date     2026-05-22
##  pandoc   3.6.3 @ /Applications/RStudio.app/Contents/Resources/app/quarto/bin/tools/aarch64/ (via rmarkdown)
##  quarto   1.7.31 @ /usr/local/bin/quarto
## 
## ─ Packages ───────────────────────────────────────────────────────────────────
##  package     * version date (UTC) lib source
##  blogdown      1.23    2026-01-18 [1] CRAN (R 4.5.2)
##  bookdown      0.46    2025-12-05 [1] CRAN (R 4.5.2)
##  bslib         0.10.0  2026-01-26 [1] CRAN (R 4.5.2)
##  cachem        1.1.0   2024-05-16 [1] CRAN (R 4.5.0)
##  cli           3.6.5   2025-04-23 [1] CRAN (R 4.5.0)
##  devtools      2.4.6   2025-10-03 [1] CRAN (R 4.5.0)
##  digest        0.6.39  2025-11-19 [1] CRAN (R 4.5.2)
##  ellipsis      0.3.2   2021-04-29 [1] CRAN (R 4.5.0)
##  evaluate      1.0.5   2025-08-27 [1] CRAN (R 4.5.0)
##  fastmap       1.2.0   2024-05-15 [1] CRAN (R 4.5.0)
##  fs            1.6.7   2026-03-06 [1] CRAN (R 4.5.2)
##  glue          1.8.1   2026-04-17 [1] CRAN (R 4.5.2)
##  htmltools     0.5.9   2025-12-04 [1] CRAN (R 4.5.2)
##  jquerylib     0.1.4   2021-04-26 [1] CRAN (R 4.5.0)
##  jsonlite      2.0.0   2025-03-27 [1] CRAN (R 4.5.0)
##  knitr         1.51    2025-12-20 [1] CRAN (R 4.5.2)
##  lattice       0.22-9  2026-02-09 [1] CRAN (R 4.5.3)
##  lifecycle     1.0.5   2026-01-08 [1] CRAN (R 4.5.2)
##  magrittr      2.0.4   2025-09-12 [1] CRAN (R 4.5.0)
##  Matrix        1.7-4   2025-08-28 [1] CRAN (R 4.5.3)
##  memoise       2.0.1   2021-11-26 [1] CRAN (R 4.5.0)
##  otel          0.2.0   2025-08-29 [1] CRAN (R 4.5.0)
##  pkgbuild      1.4.8   2025-05-26 [1] CRAN (R 4.5.0)
##  pkgload       1.5.0   2026-02-03 [1] CRAN (R 4.5.2)
##  png           0.1-9   2026-03-15 [1] CRAN (R 4.5.2)
##  purrr         1.2.2   2026-04-10 [1] CRAN (R 4.5.2)
##  R6            2.6.1   2025-02-15 [1] CRAN (R 4.5.0)
##  Rcpp          1.1.1   2026-01-10 [1] CRAN (R 4.5.2)
##  remotes       2.5.0   2024-03-17 [1] CRAN (R 4.5.0)
##  reticulate    1.45.0  2026-02-13 [1] CRAN (R 4.5.2)
##  rlang         1.1.7   2026-01-09 [1] CRAN (R 4.5.2)
##  rmarkdown     2.30    2025-09-28 [1] CRAN (R 4.5.0)
##  rstudioapi    0.18.0  2026-01-16 [1] CRAN (R 4.5.2)
##  sass          0.4.10  2025-04-11 [1] CRAN (R 4.5.0)
##  sessioninfo   1.2.3   2025-02-05 [1] CRAN (R 4.5.0)
##  usethis       3.2.1   2025-09-06 [1] CRAN (R 4.5.0)
##  vctrs         0.7.1   2026-01-23 [1] CRAN (R 4.5.2)
##  withr         3.0.2   2024-10-28 [1] CRAN (R 4.5.0)
##  xfun          0.56    2026-01-18 [1] CRAN (R 4.5.2)
##  yaml          2.3.12  2025-12-10 [1] CRAN (R 4.5.2)
## 
##  [1] /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library
## 
## ─ Python configuration ───────────────────────────────────────────────────────
##  python:         /Users/jono/.cache/uv/archive-v0/Q1veGTfRq3GBaNYBXjagV/bin/python
##  libpython:      /Users/jono/.local/share/uv/python/cpython-3.12.12-macos-aarch64-none/lib/libpython3.12.dylib
##  pythonhome:     /Users/jono/.cache/uv/archive-v0/Q1veGTfRq3GBaNYBXjagV:/Users/jono/.cache/uv/archive-v0/Q1veGTfRq3GBaNYBXjagV
##  virtualenv:     /Users/jono/.cache/uv/archive-v0/Q1veGTfRq3GBaNYBXjagV/bin/activate_this.py
##  version:        3.12.12 (main, Oct 28 2025, 11:52:25) [Clang 20.1.4 ]
##  numpy:          /Users/jono/.cache/uv/archive-v0/Q1veGTfRq3GBaNYBXjagV/lib/python3.12/site-packages/numpy
##  numpy_version:  2.4.6
##  
##  NOTE: Python version was forced by VIRTUAL_ENV
## 
## ──────────────────────────────────────────────────────────────────────────────

Over at Daring Fireball, John Gruber makes a passing observation about the Apple Sports app:

I’ve got some gripes about certain specific aspects of Apple Sports. Like, where does one even start to explain how much is wrong with their zero-sum visualization of team stats? Has anyone ever even seen a presentation like that before? Anyone?

That “Anyone” link lands over here. Hi everyone! The team stats image is quite confusing. It’s a summary of a game between the San Antonio Spurs and the Oklahoma City Thunder. I don’t know much about basketball, but I do know a bit about data visualization and in a pleasing coincidence my former student Josh Fink is the A-VP of Basketball Data Science for the Spurs. Here is the image that John objected to:

I just finished driving a very long way up the side of the country, so I’m kind of tired. But even allowing for that, boy, this way of representing things really is quite confusing. Not being an Apple Sports user I had to look at it for a bit to understand what was happening. But, now that it has given me a headache, I can kind of see why whoever designed this ended up in the undoubtedly bad place they did.

Before I get to why I have some sympathy for the designer, why did I find this representation of these numbers so disorienting? It’s not just just because I’ve been driving for nine hours. John is right to call the picture a “Zero Sum” representation. The design strongly suggests to the viewer that, within each row, we’re looking at each team’s share of a total. Each pair of black and blue lines seem to be vying for control of their whole row, with the longest line being the “winner” in each case.

This sort of representation would make perfect sense for a measure that really was zero sum. Take an example from a properly good sport, like rugby. There, like in basketball, to a first approximation a team either has the ball or it doesn’t.¹ But there’s no shot clock in rugby, and possession routinely gets turned over without the game stopping. So, knowing that Team A had 65% possession is not only informative, it also immediately entails that Team B had 35%. You could show that with a representation like one of the rows above.

Literally none of the measures in the Basketball data above are zero-sum in this way. Both teams could shoot 100% from the free throw line, or zero percent. But because the first three measures shown are percentages, this reinforces the zero-sum impression given by the lines. It certainly did that in my case. But then, starting with Assists, the remaining rows are just absolute numbers. When I started looking at the absolute numbers, I got confused a second time by the length of the lines. “Oh so it’s not a share, it’s the value” I thought—but no, they do correspond in terms of relative proportions to the teams share within each row. But they’re not really shares they’re just magnitudes. But they have to be shown in a fixed space and we want to make them relatively comparable somehow so … Argh.

It would be nice if there were One Weird Trick to fully fix this figure. But I’m not sure that there is. For example, at a minimum we could redraw these numbers to reflect the fact that they’re not zero-sum. Keep each measure as a row (i.e. on the y-axis) but have the lines, or columns, be side by side within each category instead of facing off. Like this:

This view at least lets you immediately see who “won” each measure. The viewer can just directly compare the length of the bars in each category. People are really good at doing that accurately. In that sense it’s much less confusing than the original. But there’s still a lot wrong with it. The core problem is that when we draw a graph like this, we’re usually putting the same kind of thing (e.g. countries, or religious groups, or sports teams) on the y-axis, and then seeing how different their scores are on some single measure (e.g. GDP, or number of adherents, or average points scored per game), which we put on the x-axis. Maybe we use color to break things out by some third measure as well.² In this case, I’ve just labeled the x-axis as generically as possible. “Value” covers the range of all the measures. The lowest value is 5, in Largest Lead. The highest is 88, in Free Throw %. But these numbers are not meaningfully comparable. The graph encourages us to compare across as well as within categories. But while within-category comparisons are meaningful, the between-category ones are not. There were way more Bench Points than Blocks in the game. But that is not a useful thing to know.

Knowing who won each measure isn’t nothing. It can be informative about how the game went, maybe especially when a team won the game but “lost” on a number of the measures. If you really wanted to lean in to that aspect, you could sort of justify the zero-sum view, and maybe look for a way to sort and order by “how much” a team “won” each category. But again, what’s the right denominator for those measures? For instance, do we care about a team’s share of all Defensive Rebounds in the game? Or do we care about the share of Defensive Rebounds a team won relative to every opportunity it had to make a Defensive Rebound? How meaningful is ordering our rows by those kinds of shares? Even worse, some measures (notably Fouls) are bad to “win”, so we’d have to do something about those.

Our fundamental problem is that we just have two cases (the teams) and fifteen different measures, or variables. Each variable, except for the three percentages, is in effect on its own scale. There’s no direct way to make comparisons across them. Sure, some of these measures are probably going to be associated with one another—e.g. Turnovers and Points Off Turnovers—but the numeric values aren’t directly comparable in general. If you know a lot about basketball you might have some informative rules of thumb about each one of these measures, or some of them in combination. But at that point the lines in this particular graph are not going to be doing any work for you; you’ll just end up looking directly at the numbers. If we had data on all these measures for every NBA game for a whole season then we could of course do much more with them, because then each measure would have a distribution across all games and across all teams.

As it is, the purpose of the “Stats” screen in Apple Sports is just to summarize information from a single game. The other thing I could think of to do with the numbers as kind of graph is something like this:

This is marginally more helpful than the one before just because, again, it gets rid of the unhelpful zero-sum look of the original. As I hope you can immediately see, it creates many other difficulties. It also doesn’t do away with the core problem. That problem is principally one of information design rather than data visualization. What I mean is that what we’re trying to organize is, in effect, fifteen pairs of related but fundamentally distinct numbers. If we had fifteen cases and two variables things would be simple. But with fifteen variables and two cases … well, this is not the kind of thing you can make a single effective and non-confusing graph out of. That’s why I kind of sympathize with the designer. In a constrained space they have to show thirty numbers (thirty two, including the score). Lots of information. A straight table seems like it would be boring. Surely there’s some way to thematically integrate the numbers in a visually appealing manner that brings out some of the relationships across the rows. That’s what graphs do; it seems like the right thing to reach for. But at its heart this information is not a graph. It just sort of looks like one, and that ends up confusing people.

A few days ago, I presented Conformalized TabPFN: Prediction Intervals for a Pretrained Transformer for Tabular Data in Python and R. Today, it’s about TabICL, another state-of-the-art tabular foundation model. TabICL requires no token, as you’ll notice in the following Python and R code.

1 – Python version

!pip install tabicl nnetsauce # scikit-learn matplotlib numpy

from sklearn.datasets import load_diabetes
from sklearn.model_selection import train_test_split
from sklearn.linear_model import RidgeCV
from sklearn.metrics import mean_squared_error
from tabicl import TabICLRegressor
import nnetsauce as ns
import numpy as np
import matplotlib.pyplot as plt
from time import time

# ── data ───────────────────────────────────────────────────
X, y = load_diabetes(return_X_y=True)
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, random_state=42
)

# ── base models ────────────────────────────────────────────
models = {
    "TabICL": TabICLRegressor(),
    "RidgeCV": RidgeCV(),
}

results = {}
for name, reg in models.items():
    start = time()
    conf = ns.PredictionInterval(reg, level=95)
    conf.fit(X_train, y_train)
    pi = conf.predict(X_test, return_pi=True)
    print(f"{name:10s}  time={time() - start:.1f}s")

    coverage = np.mean((pi.lower <= y_test) & (pi.upper >= y_test))
    width    = np.mean(pi.upper - pi.lower)
    rmse     = np.sqrt(mean_squared_error(y_test, pi.mean))

    results[name] = {"pi": pi, "coverage": coverage,
                     "width": width, "rmse": rmse}
    print(f"{name:10s}  RMSE={rmse:.1f}  "
          f"coverage={coverage:.3f}  avg_width={width:.1f}")

# ── plot side-by-side ──────────────────────────────────────
fig, axes = plt.subplots(1, 2, figsize=(12, 4), sharey=True)
colors = {"TabICL": "orange", "RidgeCV": "steelblue"}
max_idx = 50

for ax, (name, res) in zip(axes, results.items()):
    pi = res["pi"]
    x  = range(max_idx)
    ax.fill_between(x, pi.lower[:max_idx], pi.upper[:max_idx],
                     alpha=0.35, color=colors[name], label="95% PI")
    ax.plot(x, pi.mean[:max_idx], "k--", lw=1.5, label="predicted")
    ax.plot(x, y_test[:max_idx], "k.", ms=6, alpha=0.4, label="observed")
    ax.set_title(
        f"{name}  |  cov={res['coverage']:.3f}  width={res['width']:.1f}"
    )
    ax.legend(fontsize=8)

plt.suptitle("Conformalized TabICL vs RidgeCV — diabetes dataset")
plt.tight_layout()
plt.show()

Checkpoint 'tabicl-regressor-v2-20260212.ckpt' not cached.
 Downloading from Hugging Face Hub (jingang/TabICL).




tabicl-regressor-v2-20260212.ckpt:   0%|          | 0.00/114M [00:00<?, ?B/s]


TabICL      time=21.8s
TabICL      RMSE=54.4  coverage=0.955  avg_width=226.1
RidgeCV     time=0.0s
RidgeCV     RMSE=53.9  coverage=0.955  avg_width=211.5

2 - R version

 %load_ext rpy2.ipython # in a Colab notebook, use this

%R install.packages("reticulate")

%%R  # in Colab/Jupyter with rpy2; remove this line for pure R

library(reticulate)

# pip install tabicl nnetsauce scikit-learn matplotlib numpy

sklearn_ds  <- import("sklearn.datasets")
sklearn_ms  <- import("sklearn.model_selection")
sklearn_m   <- import("sklearn.metrics")
sklearn_lm  <- import("sklearn.linear_model")
tabicl      <- import("tabicl")
ns          <- import("nnetsauce")
np          <- import("numpy")
plt         <- import("matplotlib.pyplot")

# ── data ───────────────────────────────────────────────────
d       <- sklearn_ds$load_diabetes(return_X_y = TRUE)
X <- d[[1]]; y <- d[[2]]
sp      <- sklearn_ms$train_test_split(X, y,
             test_size = 0.2, random_state = 42L)
X_train <- sp[[1]]; X_test <- sp[[2]]
y_train <- sp[[3]]; y_test <- sp[[4]]

# ── helper: fit + evaluate ─────────────────────────────────
eval_model <- function(reg, name) {
  conf <- ns$PredictionInterval(reg, level = 95L)
  conf$fit(X_train, y_train)
  pi   <- conf$predict(X_test, return_pi = TRUE)

  cov  <- np$mean((pi$lower <= y_test) * (pi$upper >= y_test))
  wid  <- np$mean(pi$upper - pi$lower)
  rmse <- sqrt(sklearn_m$mean_squared_error(y_test, pi$mean))

  cat(sprintf("%-10s  RMSE=%.1f  coverage=%.3f  avg_width=%.1f\n",
              name, rmse, cov, wid))
  invisible(pi)
}

# ── run both models ────────────────────────────────────────
pi_tabicl  <- eval_model(tabicl$TabICLRegressor(),  "TabICL")
pi_ridge   <- eval_model(sklearn_lm$RidgeCV(),       "RidgeCV")

# ── plot ───────────────────────────────────────────────────
max_idx <- 50L
x_range <- np$array(0:(max_idx - 1))

plot_pi <- function(pi, title, col) {
  x_fill <- np$concatenate(list(x_range, x_range[max_idx:1]))
  y_fill <- np$concatenate(list(
    pi$upper[1:max_idx], pi$lower[max_idx:1]))
  plt$fill(x_fill, y_fill, alpha=0.35, fc=col, ec="None", label="95% PI")
  plt$plot(x_range, pi$mean[1:max_idx], "k--", lw=1.5, label="predicted")
  plt$plot(x_range, y_test[1:max_idx], "k.", ms=6L, alpha=0.4, label="observed")
  plt$title(title); plt$legend(fontsize=8L)
}

fig <- plt$figure(figsize=c(12, 4))
plt$subplot(1L, 2L, 1L); plot_pi(pi_tabicl, "Conformalized TabICL", "orange")
plt$subplot(1L, 2L, 2L); plot_pi(pi_ridge,  "Conformalized RidgeCV", "steelblue")
plt$suptitle("Conformalized TabICL vs RidgeCV — diabetes dataset")
plt$tight_layout()
plt$show()

    WARNING: The R package "reticulate" only fixed recently
    an issue that caused a segfault when used with rpy2:
    https://github.com/rstudio/reticulate/pull/1188
    Make sure that you use a version of that package that includes
    the fix.
    TabICL      RMSE=54.4  coverage=0.955  avg_width=226.1
RidgeCV     RMSE=53.9  coverage=0.955  avg_width=211.5

Probably a dataset that’s too easy for a Transformer. Conformalizing simple models helps them, in general, to obtain coverage rates close to the nominal level, as we see for RidgeCV here.

The manifold hypothesis, the idea that real-world high-dimensional data concentrates near a low-dimensional curved subspace, is foundational to modern machine learning. Many popular manifold learning methods such as UMAP, t-SNE, Isomap, and diffusion maps do achieve dimensionality reduction by embedding data into a flat Euclidean space , but they do not attempt to directly learn the underlying manifold. In contrast, the 2025 paper by Robinett et al., Atlas-based manifold representations for interpretable Riemannian machine learning, offers a proof of concept for directly tackling the manifold hypothesis based on fundamental ideas from differential geometry. It provides an algorithm for learning a low dimensional manifold from point data by constructing an atlas of charts. The paper is also notable for the design of an efficient data structure for working with the learned atlas and for the extensive supplementary materials that include a GitHub Repository containing several practical Python algorithms for doing calculations on manifolds, and an extraordinary amount of implementation detail.

Reading through Robinett et al., however, requires a fairly deep background in the theory of differential geometry. This post is an attempt to provide an on-ramp to Robinett et al. by discussing the relatively simple example of the two dimensional sphere, embedded in . It implements the Atlas-Learn data structures and algorithms in R, uses them to learn and then goes on to validate the Atlas-Learn algorithm for the sphere via three independent methods: 1) use numerical integration along the manifold to trace a great circle on the sphere, 3) recover the radius of curvature of the sphere from the atlas, and 4) verify the Gauss-Bonnet Theorem for the sphere.

The R code was mostly worked out by Claude Sonnet 4.3 in the context of participating in the Posit beta test for its AI Assistant. I found the integration of the AI engine into the RStudio IDE an effective means of communicating with Claude and managing the project workflow.

library(tidyverse)
library(cluster)    # pam() for k-medoids partitioning
library(RANN)       # nn2() for fast k-nearest-neighbor queries
library(plotly)     # interactive 3D visualization
library(purrr)      # map() / imap() for list operations

#| label: atlas-functions

# ===========================================================================
# PART 1: Quadratic feature helpers
# ---------------------------------------------------------------------------
# These functions implement the d=2 specialization of the general quadratic
# feature map. For general d, phi(xi) would have choose(d+1, 2) components.
# For d=2 it has exactly 3: [xi1^2, xi1*xi2, xi2^2].
# ===========================================================================

# Maps xi in R^2 to the three quadratic monomials used to model surface curvature.
# General d would give choose(d+1,2) monomials; for d=2 this is exactly 3.
quad_features <- function(xi) {
  c(xi[1]^2, xi[1] * xi[2], xi[2]^2)
}

# Jacobian of quad_features w.r.t. xi: a 3x2 matrix (d=2 specialization).
# Row i is the gradient of the i-th monomial:
#   d/dxi [xi1^2]     = [2*xi1,   0   ]
#   d/dxi [xi1*xi2]   = [xi2,     xi1 ]
#   d/dxi [xi2^2]     = [0,       2*xi2]
quad_jacobian <- function(xi) {
  matrix(
    c(2 * xi[1],   xi[2],         0,
              0,   xi[1],  2 * xi[2]),
    nrow = 3, ncol = 2, byrow = TRUE
  )
}

# ===========================================================================
# PART 2: Core chart operations
# ---------------------------------------------------------------------------
# Each chart stores:
#   mean  : R^3 — chart center (centroid of the cluster)
#   L     : R^{3x2} — orthonormal tangent basis (columns = v1, v2 from SVD)
#   M     : R^{3x1} — unit surface normal (v3 from SVD; scalar normal because D-d=1)
#   K     : R^{1x3} — quadratic curvature coefficients (1 row because D-d=1)
#   ell_A : R^{2x2} — ellipsoidal domain matrix (2x2 because d=2)
# ===========================================================================

# Evaluate the inverse chart map: xi in R^2  ->  point in R^3.
# Formula: x = mean + L*xi + M * (K * phi(xi))
#   - mean + L*xi : linear move along the tangent plane
#   - M*(K*phi)   : quadratic normal correction encoding surface curvature
chart_eval <- function(chart, xi) {
  q <- quad_features(xi)          # 3-vector: [xi1^2, xi1*xi2, xi2^2]
  as.vector(
    chart$mean +
    chart$L %*% xi +              # 3x2 * 2x1 = 3x1 tangent contribution
    chart$M * as.numeric(chart$K %*% q)   # 3x1 * scalar normal correction
  )
}

# Jacobian of the inverse chart map at xi: a 3x2 matrix.
# J(xi) = L + M * (K * dQ(xi))
#   dQ(xi) is the 3x2 Jacobian of the monomial map (quad_jacobian).
#   K * dQ  is 1x3 * 3x2 = 1x2; then M * (K*dQ) is 3x1 * 1x2 = 3x2.
# This is the key object for geodesic integration: it maps tangent vectors
# in R^2 chart coordinates to ambient R^3 velocity vectors.
chart_jacobian <- function(chart, xi) {
  dQ <- quad_jacobian(xi)                   # 3x2 monomial Jacobian
  chart$L + chart$M %*% (chart$K %*% dQ)   # 3x2 total Jacobian J(xi)
}

# Project an ambient R^3 point onto this chart's tangent-plane coordinates (R^2).
# The linear projection xi = L^T * (x - mean) is exact for points on the tangent
# plane; for points on the actual surface it is a first-order approximation.
chart_project <- function(chart, x) {
  as.vector(t(chart$L) %*% (x - chart$mean))
}

# Test whether tangent coords xi lie within the chart's ellipsoidal domain.
# The domain is {xi in R^2 : xi^T * A * xi <= 1}, a 2D Mahalanobis ball.
# This is O(d^2) = O(4) for d=2.
chart_in_domain <- function(chart, xi) {
  as.numeric(t(xi) %*% chart$ell_A %*% xi) <= 1
}

# ===========================================================================
# PART 3: atlas_learn() — the main learning function
# ---------------------------------------------------------------------------
# Input : X (N x 3 matrix of surface points), k (number of charts)
# Output: an S3 object of class "atlas" containing k atlas_chart objects
#
# The algorithm runs four steps per chart:
#   (a) k-medoids partitioning
#   (b) local PCA to find the tangent plane and normal
#   (c) quadratic regression for the curvature coefficients K
#   (d) ellipsoidal domain construction
#
# d=2, D=3 specializations that are hard-coded here:
#   - SVD takes nv=3 (the full 3x3 right-singular-vector matrix)
#   - L = V[,1:2] is 3x2; M = V[,3] is 3x1 (one normal vector, not a matrix)
#   - nu (normal offset) is a scalar per point, not a vector
#   - K is fitted as a 1x3 row vector (one row per normal direction)
#   - ell_A is 2x2 (domain is a 2D ellipse)
# ===========================================================================

atlas_learn <- function(X, k, ellipsoid_scale = 1.1) {
  # X              : N x 3 matrix of surface points
  # k              : number of charts
  # ellipsoid_scale: inflate domains by this factor so adjacent charts overlap

  message("Fitting k-medoids (k = ", k, ") ...")
  # PAM (Partitioning Around Medoids) is preferred over k-means for surfaces:
  # medoids are actual data points, making the partition robust to outliers.
  km <- cluster::pam(X, k, diss = FALSE)

  charts <- seq_len(k) |>
    purrr::map(\(j) {

      # --- (a) Extract the cluster and center it ----------------------------
      idx <- which(km$clustering == j)
      Xj  <- X[idx, , drop = FALSE]   # N_j x 3 matrix of cluster points
      m   <- colMeans(Xj)             # chart center (3-vector)
      Xc  <- sweep(Xj, 2, m)         # centered cluster: N_j x 3

      # --- (b) Local PCA: tangent plane and normal -------------------------
      # SVD of the centered cluster reveals local geometry:
      #   - first two right singular vectors (v1, v2) span the tangent plane
      #   - third right singular vector (v3) is the surface normal
      # d=2, D=3: we always take nv=3 because D=3 (fully determined).
      sv  <- svd(Xc, nu = 0, nv = 3)

      # L: 3x2 tangent basis (orthonormal columns). General form: R^{D x d}.
      L   <- sv$v[, 1:2, drop = FALSE]

      # M: 3x1 unit normal. d=2, D=3 specialization: D-d=1, so there is
      # exactly ONE normal direction and M is a column vector, not a matrix.
      M   <- sv$v[, 3, drop = FALSE]

      # Project each centered point into tangent / normal coordinates.
      # tau: N_j x 2 — tangent coordinates (2D because d=2)
      # nu:  N_j x 1 — normal offsets (scalar because D-d=1)
      tau <- Xc %*% L    # N_j x 2
      nu  <- Xc %*% M    # N_j x 1 (scalar per point; would be N_j x (D-d) in general)

      # --- (c) Quadratic curvature regression ------------------------------
      # Fit: nu ~ K * phi(tau)  where phi(tau) = [tau1^2, tau1*tau2, tau2^2]
      # K is 1x3 here (one output because D-d=1; would be (D-d)x3 in general).
      #
      # Design matrix Q: N_j x 3, each row is phi(tau_i)
      Q   <- t(apply(tau, 1, quad_features))   # N_j x 3

      # Ridge-regularized least squares: K = nu^T * Q * (Q^T*Q + esp*I)^{-1}
      # The ridge penalty eps=1e-10 prevents singularity when clusters are
      # nearly collinear in tangent coordinates.
      K   <- t(nu) %*% Q %*% solve(crossprod(Q) + 1e-10 * diag(3))  # 1 x 3

      # --- (d) Ellipsoidal domain ------------------------------------------
      # Domain: {xi in R^2 : xi^T A xi <= 1}, a Mahalanobis ball.
      # A_raw = Cov(tau)^{-1}: inverse covariance of the tangent coordinates.
      # This adapts the domain shape to the data spread in each direction.
      A_raw <- solve(cov(tau))

      # Scale A so the outermost point lands on the boundary, then inflate
      # by ellipsoid_scale (default 1.1) to ensure overlap with neighbors.
      qvals <- apply(tau, 1, \(xi) as.numeric(t(xi) %*% A_raw %*% xi))
      ell_A <- A_raw / (ellipsoid_scale * max(qvals))   # 2x2 positive-definite

      # Pack all chart data into an S3 object
      structure(
        list(mean = m, L = L, M = M, K = K, ell_A = ell_A, idx = idx,
             n_points = length(idx)),
        class = "atlas_chart"
      )
    })

  # Return the full atlas as an S3 object
  structure(
    list(
      charts        = charts,
      k             = k,
      clustering    = km$clustering,
      n_points      = nrow(X),
      ambient_dim   = ncol(X),   # D = 3
      intrinsic_dim = 2L         # d = 2 (hard-coded for surfaces in R^3)
    ),
    class = "atlas"
  )
}

# ===========================================================================
# PART 4: S3 print / summary methods
# ===========================================================================

print.atlas_chart <- function(x, ...) {
  cat(sprintf(
    "<atlas_chart>  %d points | mean [%s] | cond(ell_A) = %.1f\n",
    x$n_points,
    paste(round(x$mean, 3), collapse = ", "),
    kappa(x$ell_A)
  ))
  invisible(x)
}

# Returns a per-chart summary tibble
summary.atlas <- function(object, ...) {
  purrr::imap_dfr(object$charts, \(ch, i)
    tibble::tibble(
      chart      = i,
      n_points   = ch$n_points,
      mean_norm  = round(sqrt(sum(ch$mean^2)), 4),
      cond_ell_A = round(kappa(ch$ell_A), 1)
    )
  )
}

print.atlas <- function(x, ...) {
  cat(sprintf(
    "<atlas>  k = %d | ambient R^%d | intrinsic R^%d | %d points\n\n",
    x$k, x$ambient_dim, x$intrinsic_dim, x$n_points
  ))
  print(summary(x))
  invisible(x)
}

# ===========================================================================
# PART 5: Chart lookup
# ===========================================================================

# Return the index of the nearest chart whose domain contains x.
# Falls back to the globally nearest chart center if none qualify.
# Checking only the 6 nearest charts (by Euclidean center distance) is
# sufficient in practice and avoids an O(k) domain test at every step.
find_chart <- function(atlas, x) {
  dists      <- map_dbl(atlas$charts, \(ch) sum((x - ch$mean)^2))
  candidates <- order(dists)[seq_len(min(6L, atlas$k))]
  in_domain  <- purrr::keep(
    candidates,
    \(i) chart_in_domain(atlas$charts[[i]],
                         chart_project(atlas$charts[[i]], x))
  )
  if (length(in_domain) > 0L) in_domain[[1L]] else which.min(dists)
}

# ===========================================================================
# PART 6: Quasi-Euclidean retraction
# ---------------------------------------------------------------------------
# Advances one step from (chart_idx, xi) in the ambient direction tau_r3 (R^3).
#
# The pseudoinverse J^+ = (J^T J)^{-1} J^T is the key d=2, D=3 quantity:
#   - J is 3x2 (D x d)
#   - J^T J is 2x2 (d x d) — the pullback metric g; invertible for full-rank J
#   - J^+ is 2x3 (d x D) — pulls ambient vectors back to chart coordinates
#
# In general D > d, J^+ is the minimum-norm left inverse of J.
# ===========================================================================

retract_step <- function(atlas, chart_idx, xi, tau_r3) {
  chart  <- atlas$charts[[chart_idx]]
  J      <- chart_jacobian(chart, xi)        # 3x2 (D x d for d=2, D=3)
  Jp     <- solve(crossprod(J)) %*% t(J)    # 2x3 pseudoinverse: (J^T J)^{-1} J^T
  xi_new <- xi + as.vector(Jp %*% tau_r3)   # advance in chart coordinates

  if (chart_in_domain(chart, xi_new)) {
    # Still inside the same chart: evaluate and return
    return(list(
      chart_idx = chart_idx,
      xi        = xi_new,
      x         = chart_eval(chart, xi_new)
    ))
  }

  # Step crossed a chart boundary: find a new host chart and re-project
  x_cand <- chart_eval(chart, xi_new)       # approximate ambient position
  ci_new <- find_chart(atlas, x_cand)       # nearest chart that contains x_cand
  xi2    <- chart_project(atlas$charts[[ci_new]], x_cand)   # re-project
  list(
    chart_idx = ci_new,
    xi        = xi2,
    x         = chart_eval(atlas$charts[[ci_new]], xi2)
  )
}

# ===========================================================================
# PART 7: Geodesic path integration
# ---------------------------------------------------------------------------
# Traces a geodesic from x_start in direction direction_r3 (ambient R^3).
# Returns a tibble with columns x, y, z, step, chart_idx.
#
# At each step the ambient velocity vector is re-projected into the current
# chart's tangent plane and pushed back to ambient.  This keeps the direction
# of motion consistent across chart boundaries — without it, the fixed tau
# vector drifts after every chart transition and the path wanders off the
# intended geodesic.
# ===========================================================================

atlas_geodesic <- function(atlas, x_start, direction_r3,
                           n_steps = 100L, step_size = 0.02) {
  ci  <- find_chart(atlas, x_start)
  ch  <- atlas$charts[[ci]]
  xi  <- chart_project(ch, x_start)

  # Project the initial direction onto the starting chart's tangent plane,
  # then normalize to unit speed.  This ensures tau_r3 is always in T_x M.
  J      <- chart_jacobian(ch, xi)
  Jp     <- solve(crossprod(J)) %*% t(J)          # 2x3 pseudoinverse
  tau_r3 <- as.vector(J %*% (Jp %*% direction_r3))  # project onto tangent plane
  tau_r3 <- tau_r3 / sqrt(sum(tau_r3^2)) * step_size  # unit-speed scaling

  steps     <- vector("list", n_steps + 1L)
  chart_ids <- integer(n_steps + 1L)
  steps[[1L]]     <- chart_eval(ch, xi)
  chart_ids[[1L]] <- ci

  for (i in seq_len(n_steps)) {
    res    <- retract_step(atlas, ci, xi, tau_r3)
    ci     <- res$chart_idx
    xi     <- res$xi
    steps[[i + 1L]]     <- res$x
    chart_ids[[i + 1L]] <- ci

    # Re-project tau_r3 into the current chart's tangent plane and normalize.
    # This is the identity-transport approximation: it holds the direction
    # constant in the current chart's frame rather than applying the
    # Christoffel correction that true parallel transport would require.
    J_new  <- chart_jacobian(atlas$charts[[ci]], xi)
    Jp_new <- solve(crossprod(J_new)) %*% t(J_new)
    tau_r3 <- as.vector(J_new %*% (Jp_new %*% tau_r3))
    tau_r3 <- tau_r3 / sqrt(sum(tau_r3^2)) * step_size
    J      <- J_new
  }

  do.call(rbind, steps) |>
    as_tibble(.name_repair = "minimal") |>
    set_names(c("x", "y", "z")) |>
    mutate(step = row_number(), chart_idx = chart_ids)
}

# ===========================================================================
# PART 8: Sphere-specific helpers
# ---------------------------------------------------------------------------
# These are ground-truth functions for S^2 used only for validation;
# they play no role in the Atlas-Learn algorithm itself.
# ===========================================================================

# Great-circle geodesic distance between two unit vectors (in radians).
# d(p, q) = arccos(p . q), clamped to [-1, 1] to avoid NaN from floating point.
sphere_dist <- function(p, q) {
  acos(pmax(pmin(sum(p * q), 1.0), -1.0))
}

# Unit tangent vector at p pointing toward q along the great circle.
# Computed by projecting q onto the tangent plane at p and normalizing.
geodesic_direction_sphere <- function(p, q) {
  v   <- q - sum(p * q) * p   # remove radial component
  nrm <- sqrt(sum(v^2))
  if (nrm < 1e-12) return(NULL)
  v / nrm
}

# Run the Atlas geodesic from p toward q for the exact number of steps
# needed to cover the true great-circle distance, then measure endpoint error.
atlas_endpoint_error <- function(atlas, p, q, step_size = 0.02) {
  d   <- sphere_dist(p, q)
  dir <- geodesic_direction_sphere(p, q)

  if (d < 1e-6 || is.null(dir)) {
    return(tibble(true_dist = d, endpoint_error = 0, n_steps = 0L))
  }

  n_steps  <- max(1L, round(d / step_size))
  path     <- atlas_geodesic(atlas, p, dir,
                              n_steps = n_steps, step_size = step_size)

  endpoint <- as.numeric(path[nrow(path), c("x", "y", "z")])
  # Re-project endpoint onto S^2 before measuring angle error.
  # This separates geodesic-direction error from on-manifold drift.
  endpoint <- endpoint / sqrt(sum(endpoint^2))

  tibble(
    true_dist      = d,
    endpoint_error = sphere_dist(endpoint, q),
    n_steps        = as.integer(n_steps)
  )
}

#| code-summary: "Show the code"

#| label: generate-sphere

set.seed(4471)
N     <- 2000L
theta <- runif(N, 0, 2 * pi)
phi   <- acos(runif(N, -1, 1))   # uniform area measure on S^2

sphere_pts <- tibble(
  x = sin(phi) * cos(theta),
  y = sin(phi) * sin(theta),
  z = cos(phi)
)

X_sphere <- as.matrix(sphere_pts)

#| label: learn-atlas

atlas_cache <- "atlas_s2_cache.rds"

if (file.exists(atlas_cache)) {
  message("Loading cached atlas from ", atlas_cache)
  atlas_s2 <- readRDS(atlas_cache)
} else {
  message("Cache not found — fitting atlas (k = 25) ...")
  atlas_s2 <- atlas_learn(X_sphere, k = 25L)
  saveRDS(atlas_s2, atlas_cache)
  message("Atlas saved to ", atlas_cache)
}

# Attach chart labels to the point-cloud tibble for visualization
sphere_pts <- sphere_pts |>
  mutate(chart = factor(atlas_s2$clustering))

atlas_s2
├── k              integer          number of charts (25)
├── n_points       integer          total data points (2000)
├── ambient_dim    integer          dimension of embedding space (3 = D)
├── intrinsic_dim  integer          dimension of the manifold (2 = d)
├── clustering     integer[2000]    chart label for every data point
└── charts         list[25]         one atlas_chart per cluster
    └── charts[[j]]
        ├── mean     numeric[3]     centroid of cluster j in R³
        ├── L        numeric[3×2]   orthonormal tangent basis (columns = v₁, v₂)
        ├── M        numeric[3×1]   unit surface normal = v₃  [For d=2: single vector]
        ├── K        numeric[1×3]   quadratic curvature coefficients  [for d=2: 1 row]
        ├── ell_A    numeric[2×2]   ellipsoidal domain matrix  [d=2 specialisation: 2×2]
        ├── idx      integer[n_j]   row indices into the N×3 data matrix
        └── n_points integer        cluster size n_j

purrr::imap_dfr(atlas_s2$charts, function(ch, j) {
  tibble(
    Chart  = j,
    `n pts` = ch$n_points,
    `centre (x,y,z)` = paste(round(ch$mean, 2), collapse = ", "),
    # For the unit sphere, M should align with the radial direction (M·r̂ ≈ 1)
    `M · r̂`  = round(abs(sum(ch$M * ch$mean / sqrt(sum(ch$mean^2)))), 3),
    # Quadratic curvature coefficients
    k1 = round(ch$K[1], 3),
    k2 = round(ch$K[2], 3),
    k3 = round(ch$K[3], 3),
    # Semi-axes of the 2D ellipsoidal domain (1/sqrt of eigenvalues of ell_A)
    `ell axis a` = round(1 / sqrt(eigen(ch$ell_A, only.values=TRUE)$values[2]), 3),
    `ell axis b` = round(1 / sqrt(eigen(ch$ell_A, only.values=TRUE)$values[1]), 3)
  )
}) |> print(n = 25)

# A tibble: 25 × 9
   Chart `n pts` `centre (x,y,z)`    `M · r̂`     k1     k2     k3 `ell axis a`
   <int>   <int> <chr>                 <dbl>  <dbl>  <dbl>  <dbl>        <dbl>
 1     1      82 0.91, 0.3, -0.16      0.999 -0.216 -0.078 -0.068        0.498
 2     2      70 0.35, -0.9, -0.06     1     -0.136  0     -0.163        0.486
 3     3      84 -0.64, 0.71, -0.13    1      0.182 -0.031  0.101        0.44 
 4     4      73 -0.04, 0.95, -0.12    1      0.162 -0.024  0.112        0.463
 5     5      81 -0.75, 0.4, 0.46      1      0.115 -0.056  0.184        0.431
 6     6      73 0.28, 0.6, -0.7       0.999 -0.203  0.016 -0.114        0.518
 7     7      78 -0.93, 0.07, -0.22    1      0.169 -0.087  0.055        0.46 
 8     8      86 -0.08, -0.3, 0.91     1     -0.237  0.051  0.099        0.53 
 9     9      73 -0.18, -0.87, -0.38   1      0.111 -0.036  0.127        0.421
10    10      85 0.87, 0.07, 0.4       1     -0.148 -0.002 -0.13         0.471
11    11      75 -0.36, 0.21, 0.87     1      0.2   -0.034  0.054        0.556
12    12      67 -0.79, -0.52, -0.21   1      0.178 -0.001  0.109        0.431
13    13      84 0.04, 0.05, -0.96     1      0.182  0.012  0.081        0.446
14    14      66 -0.39, -0.83, 0.31    1     -0.189 -0.006 -0.066        0.511
15    15      84 0.28, -0.55, -0.74    1      0.132 -0.066  0.105        0.431
16    16      87 -0.52, 0.43, -0.68    1     -0.251  0.009  0.026        0.511
17    17     100 0.83, -0.44, -0.2     1      0.144  0.082  0.133        0.526
18    18      68 0.07, -0.8, 0.55      1     -0.24   0.072 -0.088        0.474
19    19      91 0.59, 0.71, 0.24      0.999  0.229 -0.002 -0.003        0.537
20    20      93 -0.81, -0.26, 0.45    1      0.17   0.084  0.115        0.455
21    21      95 0.41, 0.27, 0.83      1      0.155 -0.087  0.056        0.532
22    22      66 -0.46, -0.32, -0.78   1     -0.179 -0.025 -0.072        0.525
23    23      70 -0.12, 0.78, 0.54     1      0.183  0.045  0.1          0.495
24    24      93 0.67, 0.07, -0.69     1      0.192 -0.001  0.116        0.453
25    25      76 0.62, -0.45, 0.58     1     -0.215 -0.039  0.028        0.512
# ℹ 1 more variable: `ell axis b` <dbl>

#| label: atlas-chart-plot

n_charts <- length(unique(sphere_pts$chart))
chart_colors <- colorRampPalette(RColorBrewer::brewer.pal(9, "Set1"))(n_charts)

plot_ly(
  data   = sphere_pts,
  x = ~x, y = ~y, z = ~z,
  type   = "scatter3d",
  mode   = "markers",
  color  = ~chart,
  colors = chart_colors,
  marker = list(size = 2.5, opacity = 0.6),
  text   = ~paste0("Chart ", chart),
  hoverinfo = "text"
) |>
  layout(
    title  = "S² coloured by Atlas-Learn chart membership (k = 25)",
    scene  = list(
      xaxis = list(title = "x"),
      yaxis = list(title = "y"),
      zaxis = list(title = "z"),
      aspectmode = "cube"
    ),
    showlegend = FALSE
  )

Here I am again, writing about crochet and programming! I’ve continued creating my tons of shitty stitches cute creatures.

New Git analogy! The crochet lifeline

First Git/crochet analogy.

I read a great crochet book about creating your own amigurumi patterns: Créer ses propres modèles d’amigurumis au crochet by Clotilde Massot and Lise Grandjonc. One of its authors (Clotilde Massot) is a software developer who, among other things, published an octocat pattern that I bought and used… Anyway in the book they explain that when you create a pattern, you will probably have to undo your work several times. Undoing a round is easy if you have a stitch marker in the first stitch of the current/last round: you undo until you hit that stitch marker. But what about undoing several rounds? In that case, you’ll be better off if instead of using a stitch marker at the beginning of the current round, you use contrasting yarn stuck under the first stitch of each round. That yarn creates what the authors of the book call your lifeline!

Now, if that’s not a great analogy for commits and the ability to reset…

Communities of practice

If you’ve had the opportunity to attend or watch the wonderful useR! 2025 keynote talk of my rOpenSci colleague Yanina Bellini Saibene, you’ve heard of communities of practice. In her talk, Yani quoted Etienne Wenger who defined communities of practice as “groups of people who share a passion for something that they know how to do, and who interact regularly in order to learn how to do it better“. Yani mentioned her swimming team and English conversation club. Well, I found a community of practice for crochet: a stitch club at a local café! Participants meet up to crochet side by side, talking a lot about crochet: comparing yarns, exchanging tips, etc. The first time I went, we even started with a round of introductions where we said what each of us would work on during that meeting, which reminded me of rOpenSci coworking sessions where participants do exactly that.

Usefulness of seeing others’ work

I’ve worked on some patterns by Yan Schenkel a.k.a. Pica Pau. One cool aspect of the patterns is that they include a link to a gallery where anyone can upload pics of their take on each creature. So you can look at them, maybe seeing a crucial (to you) angle that’s absent from the pattern, noticing whether some “flaw” of your own project is present in others’ projects, comparing variations and picking what you prefer before you start, etc. For instance I stared at many pics of Alberto Seagull before making the legs for mine.

The usefulness of seeing others’ crocheted animals reminds me of how reading and reviewing open-source, or our colleagues’, code helps us (and LLMs, I suppose) learn how to do, or not to do, some things, how it help us refine our taste! Ironically, I do publish most of my code, but I haven’t had interest in doing that for my crocheting yet.

Conclusion

In summary, I keep finding excuses to talk about crochet. My pink octocat approves!

Note for myself and others: how to remove the “Leaflet | ” prefix in map attribution using the R package {leaflet}.

According to some sources we could write leaflet(options = leafletOptions(attributionPrefix = "")) but it doesn’t work in my case.

So instead we can execute some javascript:

library(leaflet)

leaflet() |>
  addTiles(urlTemplate = "", attribution = "Only my data") |> 
  htmlwidgets::onRender("function(el, x) {
    // Remove Leaflet attribution prefix
    this.attributionControl.setPrefix('');
    }")

Introduction

rOpenSci’s R-universe system is an open source platform allowing users to create their own CRAN-like universe of R packages.

It is absolutely fantastic. It is particularly useful in one area I research, Mendelian randomization (at the interface of Epidemiology and Genetic Epidemiology), because a lot of the packages are GitHub/GitLab-only.

Therefore, I setup and maintain https://mrcieu.r-universe.dev/ to include both packages from our MRCIEU GitHub organisation (from the MRC Integrative Epidemiology Unit at the University of Bristol, UK), and as many of the GitHub-only packages for Mendelian randomization I could find.

It is difficult to overstate how useful this is. For the first time, not only do researchers have a list of the Mendelian randomization packages in one place but they can install binaries – without having to go through the hassle – especially on (Ubuntu) Linux – of remotes::install_github(). Researchers can also see how often packages are updated and R-universe checks for changes in packages approximately every hour, keeping it always up to date.

This post gives five tips I have developed to help manage my R-universe.

Tip 1: Referring to a package from a pull request instead of from a branch on a fork

In the Mendelian randomization field many of these GitHub-only packages are not well written or abandoned once the PhD student/researcher leaves. Often when I add a package to our R-universe I find that their build fails, or they have R CMD check errors and warnings, or after several months their build fails because they are not maintained. I sometimes look into the failed builds and check problems. If it’s clear just a few fixes are required to rectify the situation I often open a pull request. Often that pull request is not responded to.

Previously, for such cases I would switch the source of the package entry in packages.json to be from the relevant branch on my fork. However, I have always felt a bit uneasy about this. I wondered if GitHub had a way to refer to the pull request branch without having to switch the repository. It turns out that it does. The format of pull request branch names is refs/pull/{number}/head where {number} is the number assigned once the PR is opened. Therefore, when I open a PR on a package I now add the "branch" field to the package entry in packages.json as follows.

  {
    "package": "GWASBrewer",
    "url": "https://github.com/jean997/GWASBrewer",
    "branch": "refs/pull/18/head"
  },

I switch back to the default branch if the PR is merged.

Tip 2: Justfile recipe for adding a package to packages.json

I regularly find that I need to add or remove a package. Manually editing the packages.json file is not hard, but I have found the following Justfile (Just is like Make, but specifically designed for running commands and has a much friendlier syntax) recipes helpful for doing this quickly.

These recipes require uv and just to be installed and on your PATH (uv automatically installs the required version of Python and creates/destroys/manages any required virtual environments). To use them, copy them into a text file named justfile at the top level of your R-universe registry repository and follow the instructions.

This recipe adds a package to your packages.json in alphabetical order. It has one required argument and 3 optional arguments.

# add a package entry to packages.json in alphabetical order
[arg("branch", short="b")]
[arg("pkgname", short="p")]
[arg("subdir", short="s")]
add url pkgname="" branch="" subdir="":
    #!/usr/bin/env -S uv run --python 3.14 python3
    import json, re, sys
    url = "{{ url }}"
    if re.fullmatch(r'[^/]+/[^/]+', url):
        url = f"https://github.com/{url}"
    pkgname = "{{ pkgname }}" or url.rstrip("/").split("/")[-1]
    branch = "{{ branch }}"
    subdir = "{{ subdir }}"
    with open("packages.json") as f:
        packages = json.load(f)
    if any(p["package"] == pkgname for p in packages):
        print(f"Error: '{pkgname}' already exists in packages.json", file=sys.stderr)
        sys.exit(1)
    entry = {"package": pkgname, "url": url}
    if branch:
        entry["branch"] = branch
    if subdir:
        entry["subdir"] = subdir
    packages.append(entry)
    packages.sort(key=lambda p: p["package"].lower())
    with open("packages.json", "w") as f:
        json.dump(packages, f, indent=2)
        f.write("\n")
    print(f"Added {pkgname}")

Where url is say https://github.com/MRCIEU/TwoSampleMR, except that for GitHub packages you can specify this as MRCIEU/TwoSampleMR.

To add a GitHub package whose name matches its repository name, simply run

just add username/reponame

You can inspect the recipe’s arguments and options with

just --usage add

Usage: just add [OPTIONS] url

Arguments:
  url

Options:
  -p pkgname [default: ""]
  -b branch [default: ""]
  -s subdir [default: ""]

The 3 optional arguments allow you to specify the package name (-p pkgname), branch (-b branchname), or subdirectory (-s subdirectory) the package is in. For example, to add a GitHub package whose package name does not match its repository name run

just add username/reponame -p pkgname

Tip 3: Justfile recipe for removing a package from packages.json

This recipe removes a package from your packages.json.

# remove a package entry from packages.json
remove pkgname:
    #!/usr/bin/env -S uv run --python 3.14 python3
    import json, sys
    pkgname = "{{ pkgname }}"
    with open("packages.json") as f:
        packages = json.load(f)
    filtered = [p for p in packages if p["package"] != pkgname]
    if len(filtered) == len(packages):
        print(f"Error: '{pkgname}' not found in packages.json", file=sys.stderr)
        sys.exit(1)
    with open("packages.json", "w") as f:
        json.dump(filtered, f, indent=2)
        f.write("\n")
    print(f"Removed {pkgname}")

Run it with

just remove pkgname

Tip 4: Justfile recipe for checking packages.json is valid

When manually editing packages.json it is very easy to forget a comma or to miss a closing bracket or quotation mark. This recipe checks your JSON is valid.

# check packages.json
check:
    uv run --python 3.14 -m json.tool packages.json > /dev/null && echo "JSON check passed"

Run it with

just check

Tip 5: Conveniently view a package’s dependencies

Knowing a package’s full strong dependency list is useful — for example, when a breaking change somewhere in the chain causes unexpected build failures. While there are several ways to determine this in R, R-universe shows you the full list immediately.

Navigate to the R-universe page for the package you are interested in, say https://mrcieu.r-universe.dev/TwoSampleMR and click the dependencies pill.

It expands showing the full dependency list.

Summary

In summary, I have shown five tips I find useful to manage a large R-universe.

Knowing a model’s prediction is useful. Knowing how confident that prediction is, even more so. Conformal prediction provides exactly that: statistically valid prediction intervals with guaranteed coverage (under certain conditions), regardless of the underlying model or data distribution.

In this post, we pair two powerful tools: TabPFN, a pretrained transformer for tabular data, and nnetsauce’s PredictionInterval (which implements Split Conformal Prediction), which wraps any scikit-learn-compatible regressor into a conformal predictor. We demonstrate the full pipeline on the diabetes dataset, first in Python, then in R via reticulate. Both versions produce identical results: a coverage rate of 96.7% at a nominal 95% level.

1 – Python version

!pip install tabpfn tabpfn_client

!pip install nnetsauce

import tabpfn_client

API_TOKEN = "" # <- Paste your TabPFN token here (from https://priorlabs.ai/tabpfn)


tabpfn_client.set_access_token(API_TOKEN)

from sklearn.datasets import load_diabetes
from sklearn.model_selection import train_test_split
from tabpfn_client import TabPFNRegressor
from sklearn.metrics import mean_squared_error
import numpy as np

reg = TabPFNRegressor()

X, y = load_diabetes(return_X_y=True)
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

reg.fit(X_train, y_train)
preds = reg.predict(X_test)
rmse = np.sqrt(mean_squared_error(y_test, preds))
print(-rmse)

00:00 Fitting... |

WARNING:tabpfn_client.client:The provided train set hashes match previously uploaded train sets.


00:00 Fitting... Done!
00:00 Predicting... -

WARNING:tabpfn_client.client:The provided test set hash matches a previously uploaded test set.


00:01 Predicting... Done!
-51.559912022529886

import nnetsauce as ns

reg_conformal = ns.PredictionInterval(reg, level=95)
reg_conformal.fit(X_train, y_train)
preds = reg_conformal.predict(X_test, return_pi=True)

00:00 Fitting... |

WARNING:tabpfn_client.client:The provided train set hashes match previously uploaded train sets.


00:00 Fitting... Done!
00:00 Predicting... -

WARNING:tabpfn_client.client:The provided test set hash matches a previously uploaded test set.


00:01 Predicting... Done!
00:00 Predicting... -

WARNING:tabpfn_client.client:The provided test set hash matches a previously uploaded test set.


00:01 Predicting... Done!
00:00 Predicting... -

WARNING:tabpfn_client.client:The provided test set hash matches a previously uploaded test set.


00:01 Predicting... Done!

print(f"coverage_rate: {np.mean((preds.lower<=y_test)*(preds.upper>=y_test))}")

coverage_rate: 0.9662921348314607

import warnings
import matplotlib.pyplot as plt


warnings.filterwarnings('ignore')

split_color = 'green'
split_color2 = 'orange'
local_color = 'gray'

def plot_func(x,
              y,
              y_u=None,
              y_l=None,
              pred=None,
              shade_color="",
              method_name="",
              title=""):

    fig = plt.figure()

    plt.plot(x, y, 'k.', alpha=.3, markersize=10,
             fillstyle='full', label=u'Test set observations')

    if (y_u is not None) and (y_l is not None):
        plt.fill(np.concatenate([x, x[::-1]]),
                 np.concatenate([y_u, y_l[::-1]]),
                 alpha=.3, fc=shade_color, ec='None',
                 label = method_name + ' Prediction interval')

    if pred is not None:
        plt.plot(x, pred, 'k--', lw=2, alpha=0.9,
                 label=u'Predicted value')

    #plt.ylim([-2.5, 7])
    plt.xlabel('$X$')
    plt.ylabel('$Y$')
    plt.legend(loc='upper right')
    plt.title(title)

    plt.show()


max_idx = 50
plot_func(x = range(max_idx),
          y = y_test[0:max_idx],
          y_u = preds.upper[0:max_idx],
          y_l = preds.lower[0:max_idx],
          pred = preds.mean[0:max_idx],
          shade_color=split_color2,
          title = f"conformalized TabPFN ({max_idx} first points in test set)")

2 – R version

For this R version, I used R in the same notebook as Python, in Google Colab.

%load_ext rpy2.ipython

%R install.packages("reticulate")

%%R

# Conformalized TabPFN in R via reticulate

library(reticulate)

# ── 0. Python environment ──────────────────────────────────────────────────────
# Use your preferred Python env. Uncomment one (automatic on Google Colab):
# use_python("/usr/bin/python3")
# use_virtualenv("r-tabpfn")
# use_condaenv("r-tabpfn")

# Install required packages into the active Python env (run once)
# py_install(c("tabpfn", "tabpfn_client", "nnetsauce", "scikit-learn",
#              "matplotlib", "numpy"), pip = TRUE)

# ── 1. Imports ─────────────────────────────────────────────────────────────────
sklearn_datasets  <- import("sklearn.datasets")
sklearn_model_sel <- import("sklearn.model_selection")
sklearn_metrics   <- import("sklearn.metrics")
tabpfn_client     <- import("tabpfn_client")
ns                <- import("nnetsauce")
np                <- import("numpy")
plt               <- import("matplotlib.pyplot")
warnings          <- import("warnings")

# ── 2. TabPFN API token ────────────────────────────────────────────────────────
API_TOKEN <- ""   # <-- paste your TabPFN token here (from https://priorlabs.ai/tabpfn)
tabpfn_client$set_access_token(API_TOKEN)

TabPFNRegressor <- tabpfn_client$TabPFNRegressor

# ── 3. Data ────────────────────────────────────────────────────────────────────
diabetes   <- sklearn_datasets$load_diabetes(return_X_y = TRUE)
X          <- diabetes[[1]]
y          <- diabetes[[2]]

split      <- sklearn_model_sel$train_test_split(X, y, test_size = 0.2, random_state = 42L)
X_train    <- split[[1]]
X_test     <- split[[2]]
y_train    <- split[[3]]
y_test     <- split[[4]]

# ── 4. Fit TabPFN regressor ────────────────────────────────────────────────────
reg   <- TabPFNRegressor()
reg$fit(X_train, y_train)
preds_plain <- reg$predict(X_test)

rmse <- sqrt(sklearn_metrics$mean_squared_error(y_test, preds_plain))
cat(sprintf("TabPFN RMSE: %.4f\n", rmse))

# ── 5. Conformal prediction with nnetsauce ─────────────────────────────────────
reg_conformal <- ns$PredictionInterval(reg, level = 95L)
reg_conformal$fit(X_train, y_train)
preds <- reg_conformal$predict(X_test, return_pi = TRUE)

coverage <- np$mean((preds$lower <= y_test) * (preds$upper >= y_test))
cat(sprintf("Coverage rate: %.4f\n", coverage))

# ── 6. Plot (first 50 test points) ────────────────────────────────────────────
warnings$filterwarnings("ignore")

max_idx    <- 50L
x_range    <- np$array(0:(max_idx - 1))   # numeric index
y_obs      <- y_test[1:max_idx]
y_upper    <- preds$upper[1:max_idx]
y_lower    <- preds$lower[1:max_idx]
y_pred     <- preds$mean[1:max_idx]

# Build the filled polygon (matplotlib-style concatenation)
x_fill <- np$concatenate(list(x_range, x_range[max_idx:1]))
y_fill <- np$concatenate(list(y_upper, y_lower[max_idx:1]))

fig <- plt$figure()
plt$plot(x_range, y_obs,  "k.", alpha = 0.3, markersize = 10L,
         label = "Test set observations")
plt$fill(x_fill, y_fill, alpha = 0.3, fc = "orange", ec = "None",
         label = "Conformal Prediction interval")
plt$plot(x_range, y_pred, "k--", lw = 2L, alpha = 0.9,
         label = "Predicted value")
plt$xlabel("Index")
plt$ylabel("Y")
plt$legend(loc = "upper right")
plt$title(sprintf("Conformalized TabPFN (first %d points in test set)", max_idx))
plt$tight_layout()
plt$show()
# To save instead: plt$savefig("conformalized_tabpfn.png", dpi = 150L)

00:02 Fitting... Done!
00:02 Predicting... Done!
TabPFN RMSE: 51.5599
00:01 Fitting... Done!
00:02 Predicting... Done!
00:00 Predicting... -

WARNING:tabpfn_client.client:The provided test set hash matches a previously uploaded test set.


00:01 Predicting... Done!
00:02 Predicting... Done!
Coverage rate: 0.9663

Exploring the CovR/S two-component system in Group A Strep — from genome annotation with Bakta & BaktFold, to AlphaFold confidence metrics, and a first attempt at protein docking with Haddock3. Learning as we go!

Motivations

Since the last ampC adventure, I’m really curious about the mechanism of some of these bacterial virulence. Remember how chromosomal ampC organisms use ampG, ampD, and then ampR to repress class C beta lactamase gene? It’s such an orchestrated endeavor. What about streptococcus pyogenes and its virulence? How can it be a colonizer on one end and then virulence on the other that caused a number of devastating infection? Let’s learn a bit of the mechanism, and of course why not use this opportunity too to learn some other bioinformatic tools along the way? And see if we can use existing knowledge to make it more fun and educational! I’m looking forward to this! Join me in exploring the mechanism of the CovR/S two-component system in streptococcus pyogenes, aka Group A strep!

Objectives:

• What is CovR/S Two-component System?
• Let’s Look A Where Does It Show in NCBI
• What If We Have WGS? How to Annotate?
• What Are Hypothetical?
  • What Is An Acceptable AlphaFold Confidence?
  • A New Tool Called BaktFold
• Do All Streptococcus Pyogenes Have CovR/S Two-component System?
• Do Other Streptococcus species Have CovR/S Two-component System
• What Does CovR Look Like?
  • What would a Phosphorylated CovR Look like?
• Opportunities For Improvement
• Lessons Learnt

What is CovR/S Two-component System?

The CovR/S (Control of Virulence) system functions as a sophisticated environmental sensor that integrates multiple host-derived signals to orchestrate the transition from colonization to invasive disease. The evidence reveals three primary environmental triggers that modulate this master regulatory system. Simplistically, CovS (Sensor) senses first -> gets activated -> trigger CovR (Regulator) -> downstream repression. Breakage of such system will de-repress the virulence factors.

1. Magnesium Levels: The Baseline Sensor High extracellular magnesium concentrations (typical of healthy tissue) activate CovS kinase activity, leading to increased CovR phosphorylation and repression of virulence genes. This creates a colonization-friendly state where GAS maintains low virulence factor expression suitable for asymptomatic carriage.

2. LL-37 Antimicrobial Peptide: The Invasion Signal LL-37 cathelicidin peptide — released by neutrophils and epithelial cells during inflammation — directly binds to the extracellular domain of CovS and inhibits its kinase activity. This creates a paradoxical host-pathogen interaction where the host’s antimicrobial defense actually triggers bacterial virulence. LL-37 binding to CovS reduces CovR phosphorylation, leading to derepression of multiple virulence factors including pyrogenic exotoxin A, DNase Sda1, streptolysin O, and hyaluronic acid capsule. Critically, LL-37 signaling converts GAS from a colonizing to an invasive phenotype, with marked increases in resistance to opsonophagocytic killing by human leukocytes.

3. Acidic Stress: The Tissue Environment Sensor Acidic conditions (pH < 7.0, typical of infected or inflamed tissue) enhance CovR/S-dependent gene repression through activation of the covR/S promoter itself. This creates a negative feedback loop where tissue acidosis increases CovR/S expression, which then more strongly represses virulence factors.

Below is an image referenced directly from source that depicts the mechanism of CovR/S system in streptococcus pyogenes.

Let’s Look A Where Does It Show in NCBI

Let’s go to here. I picked out streptococcus pyogenes reference genome annotation and search for cov and this popped up.

There you go! CovS and CovR. Sometimes these system can also be known as CsrR/CsrS (Capsule Synthesis Regulator). There may have been 2 different research groups discovered these identical gene and called it differently?

It’s also so interesting that these 2 genes are so close to each other.

What If We Have WGS? How to Annotate?

Alright, let’s pick a random streptococcus pyogenes and see if we can use bakta to help us annotate. Let’s look at this one.

Install Bakta, see here

#make sure you use the environment name you created
conda activate bakta_env 

bakta \
  --db /path/to/bakta_db \
  --output rabdom_bakta_output \
  --prefix random \
  --threads 8 \
  --skip-crispr \
  --force \
  random.fna

After it is done, when you look in the folder, you will see something like this

When we look at the annotated gff3 file, we can see that there are 2 features annotated as two-component system response regulator and two-component system sensor histidine kinase. These are likely to be CovR and CovS.

library(ape)
library(tidyverse)

readLines("random.gff3") |> str_detect("FASTA") |> which() #found fasta, apparently bakta has ###FASTA inserted and ape cannot handle

## [1] 2025

tmp <- tempfile()
readLines("random.gff3")[1:2024] |> writeLines(tmp)
gff <- read.gff(tmp, GFF3 = T)
gff |>
  filter(str_detect(attributes, "[Cc][Oo][Vv]"))

##        seqid    source type  start    end score strand phase
## 331 contig_1 Pyrodigal  CDS 303144 303830    NA      +     0
## 332 contig_1 Pyrodigal  CDS 303836 305338    NA      +     0
##                                                                                                                                                                                                                             attributes
## 331           ID=EEBJGP_00327;Name=two-component system response regulator CovR;locus_tag=EEBJGP_00327;product=two-component system response regulator CovR;Dbxref=BlastRules:WP_002991052,SO:0001217,UniRef:UniRef50_Q49XM7;gene=covR
## 332 ID=EEBJGP_00328;Name=two-component system sensor histidine kinase CovS;locus_tag=EEBJGP_00328;product=two-component system sensor histidine kinase CovS;Dbxref=BlastRules:WP_002991036,SO:0001217,UniRef:UniRef50_D3KVE8;gene=covS

There you go! We found them after annotation. Wait a minute… what are those hypotheticals on our folder? Why are they there? Are they important? Let’s find out.

What Are Hypothetical?

Hypotheticals are those proteins that we don’t know what they do. They are annotated as “hypothetical protein” because they are predicted to be proteins based on the DNA sequence, but we have no experimental evidence of their function. They are often annotated as “hypothetical” because they have no known homologs in other organisms, or because they have no known domains or motifs that can be used to predict their function.

When we take a peek at the random.hypotheticals.tsv, it looks like this:

Let’s check how many hypotheticals we have here for this genome

tmp <- tempfile()
readLines("random.hypotheticals.tsv")[3:194] |> writeLines(tmp)
hypo <- read_tsv(tmp)
nrow(hypo)

## [1] 191

OK we have 191 of hypotheticals. Let’s see if a new tool on the block will be able to add some annotation to these hypotheticals and see if we can find anything interesting. We could also use filter and see if we can see those hypotheticals

gff |>
  filter(str_detect(attributes,"hypothetical protein")) |>
  nrow()

## [1] 192

Hmm.. they don’t tally. But let’s move on.

A New Tool Called BaktFold

BaktFold is a new tool that uses AlphaFold to predict the structure of proteins and then uses that structure to predict the function of the protein. It is a very powerful tool that can be used to annotate hypothetical proteins. Check out their paper

baktfold run \
  -i random.json \
  -o random_baktfold_output \
  -d /path/to/baktfold_db \
  -t 8 \
  -f 

readLines("baktfold.gff3") |> str_detect("FASTA") |> which() #found fasta, apparently bakta has ###FASTA inserted and ape cannot handle

## [1] 2025

tmp <- tempfile()
readLines("baktfold.gff3")[1:2024] |> writeLines(tmp)
gff_baktfold <- read.gff(tmp, GFF3 = T)
gff_baktfold |>
  filter(str_detect(attributes,"hypothetical protein")) |>
  nrow()

## [1] 106

Wow, this is really cool! We can see that we have much less hypotheticals! About 86 less! Let’s take a look what were previous hypotheticals and what they are annotated now and how?

gff_hypo <- gff |>
  filter(str_detect(attributes,"hypothetical protein")) |>
  pull(start)

new_df <- tibble(start = gff_hypo, temp = NA)

gff_baktfold_hypo <- gff_baktfold |>
  right_join(new_df) |>
  mutate(temp = case_when(
    !str_detect(attributes,"hypothetical protein") ~ attributes,
    TRUE ~ temp
  ))

## Joining with `by = join_by(start)`

Take a look at the image above, if you see temp column with NA, that would be hypotheticals from bakta. If you see it filled, it means baktfold was able to identify it across one or more databases. It’s interesting how baktfold work, it conducts sequential protein structure-based searches against four complementary structure databases (SwissProt, Alphafold Cluster Database, PDB, CATH). Protein sequences are transformed into Foldseek 3Di tokens via the ProstT5 protein language model and subsequently searched against structure databases via Foldseek. Pretty cool! Also, interestingly, when comparing a few of the baktfold predicted functional proteins with NCBI’s annotation, we sometimes do see some baktfold-annotated functions whereas NCBI labeled them as uncharacterized gene. This is not an exhaustive or thorough comparison by any means, but interesting to note.

Speaking of AlphaFold, we’ve always wanted to know a bit more about AlphaFold confidence. When we look at the predicted structure of a protein, how do we know if we can trust it? What is an acceptable AlphaFold confidence? Let’s learn a bit more.

What Is An Acceptable AlphaFold Confidence?

**These are edited responses from back-and-forth Claude Sonnet 4.6 query, **

AlphaFold reports per-residue confidence as per-residue confidence as pLDDT (predicted local distance difference test), scored 0–100:

It tells you how well-placed AlphaFold thinks each amino acid is relative to nearby residues. Crucially, it is not a measure of experimental validation — it is the model’s self-assessed confidence.

For docking, pLDDT is essentially a proxy for how much you can trust the binding pocket geometry. Active site residues need pLDDT ≥ 90 ideally, with ≥ 70 as a minimum. Second-shell residues within ~8 Å should also clear 70, and any low-confidence loops capping the binding site entrance are a red flag even if the catalytic residues themselves look fine. Meaning, might be a good idea to visualize with B factor in ChimeraX to ensure the binding sites are acceptable.

For MD, it seems to be a bit more forgiving. Regions scoring 70–90 will generally equilibrate fine — the force field redistributes strain and lets uncertain side chains settle. Regions in the 50–70 band need longer equilibration (100+ ns) and staged restraint release to avoid unphysical collapse (interesting area to explore). Below 50, MD often can’t rescue the geometry — these regions should either be truncated if non-essential, cross-checked against other databases, or explored with enhanced sampling methods. OK what about PAE on their website?

PAE (predicted aligned error) is the second major confidence metric AlphaFold produces, and it tells you something fundamentally different from pLDDT. Where pLDDT is a per-residue score asking “how confident am I in this residue’s local geometry,” PAE is a pairwise score asking “how confident am I in the relative position and orientation of residue A with respect to residue B.” It’s an N×N matrix where every cell (i,j) contains the expected position error in Å for residue j when residue i is used as the alignment reference.

Why it matters? pLDDT can look great across an entire protein — every residue scores above 80 — but if the PAE between two domains is high, that confident-looking structure is misleading. The two domains are individually well-folded, but AlphaFold is telling you it has no idea how they pack against each other. For docking, check the PAE within the domain containing your binding site — you want a dark block there, confirming the domain’s internal geometry is reliable as a unit. For MD, high inter-domain PAE is a heads-up that you may need enhanced sampling to explore the conformational space between domains rather than assuming the AlphaFold pose is the dominant one.

Two metrics before trusting AlphaFold. pLDDT is local — want ≥ 70 at active site, ≥ 90 for catalytic residues. PAE is pairwise — dark green means confident relative positioning between any two residues. Single-chain: check diagonal at binding site. Multi-chain: off-diagonal blocks tell you if the predicted interface is real. Check both before docking or MD.

Do All Streptococcus Pyogenes Have CovR/S Two-component System?

Interestingly, after downloaded 2840 Streptococcus pyogenes refseq annotation feature (gff3), I found 99.96% (2839/2840) of these contain CovR/S genes. Why not 100%? Turns out to be this GCF_005472355.1 that doesn’t have CovR/S listed in the annotation. I used Bakta to annotate it, and still no luck. Then used Baktfold to further annotate, couldn’t find it either. Used tblastn to look for CovR/S protein, the best return was 42% identity. Wow, does this isolate really have absent gene for those 2?

Note: Expand below to see the utility of using –dehydrate and then rehydrate for downloading and annotating large number of genomes.

datasets download genome taxon "Streptococcus pyogenes" \
  --annotated \
  --assembly-source refseq \
  --include gff3 \
  --dehydrated \
  --filename strep_pyogenes_refseq_gff3.zip

datasets rehydrate --directory strep_pyogenes_refseq_gff3/ --max-workers 10

Do Other Streptococcus species Have CovR/S Two-component System?

We’ve seen what other streptococcus can do clinically, I wonder if they have similar system when compared to streptococcus pyogenes. Let’ download all reference gene of streptococcus genus and see if we can find it in their annotations.

datasets download genome taxon "Streptococcus" \
  --reference \
  --include gff3 \
  --assembly-source RefSeq \
  --dehydrated \
  --filename strep_dehydrated.zip

unzip strep_dehydrated.zip

datasets rehydrate --directory strep_dehydrated

df_cov <- read_csv("covr_covs_annotation.csv")

df_cov_single <- df_cov |>
  filter(covR_present == "yes" & covS_present == "no") |> 
  pull(organism) 

df_cov |> head(10)

## # A tibble: 10 × 8
##    accession  organism covR_present covS_present covR_gene_names covS_gene_names
##    <chr>      <chr>    <chr>        <chr>        <chr>           <chr>          
##  1 GCF_00002… Strepto… no           no           <NA>            <NA>           
##  2 GCF_00016… Strepto… no           no           <NA>            <NA>           
##  3 GCF_00018… Strepto… no           no           <NA>            <NA>           
##  4 GCF_00018… Strepto… no           no           <NA>            <NA>           
##  5 GCF_00018… Strepto… yes          no           covR            <NA>           
##  6 GCF_00018… Strepto… yes          no           covR            <NA>           
##  7 GCF_00022… Strepto… no           no           <NA>            <NA>           
##  8 GCF_00025… Strepto… no           no           <NA>            <NA>           
##  9 GCF_00037… Strepto… no           no           <NA>            <NA>           
## 10 GCF_00037… Strepto… no           no           <NA>            <NA>           
## # ℹ 2 more variables: covR_annotation <chr>, covS_annotation <chr>

Wow, amonng the streptococcus species (n=153) streptococcus pyogenes appears to be the ONLY species that contains both CovR/S two-component system?

Interestingly, there were 17 that have covR but not covS.

These are Streptococcus parauberis NCFD 2020, Streptococcus ictaluri 707-05, Streptococcus didelphis DSM 15616, Streptococcus castoreus DSM 17536, Streptococcus iniae, Streptococcus phocae, Streptococcus bovimastitidis, Streptococcus catagoni, Streptococcus equi subsp. zooepidemicus, Streptococcus dysgalactiae, Streptococcus halichoeri, Streptococcus penaeicida, Streptococcus hongkongensis, Streptococcus porcinus, Streptococcus uberis, Streptococcus canis, Streptococcus pseudoporcinus. How curious!

What Does CovR Look Like?

Let’s take a look at Alpha Fold Database. Looking at CovR active site, mutation of D53A showed no dimerization of CoVR, which means that is a phosphylation site. Dimerization means that another phosphorylated CovR molecule is binding to another phosphorylated CovR, which is the active form of CovR. Which then binds to the DNA and represses the virulence genes.

The above shows both the PAE, plDDT and the predicted structure. The PAE is pretty good, the plDDT is also pretty good, with most of the residues above 70. The predicted structure also looks pretty reasonable, with the active site D53 highlighted in red.

What would a Phosphorylated CovR Look like?

Apparently for a CovR with D53E mutation (so it mimics the phosphorylated state), the structure is quite different from the wild type. The D53E mutation causes a conformational change that allows CovR to dimerize and bind to DNA, even in the absence of phosphorylation. Let’s visualize a D53E CovR predicted by AF3.

We can see that the highlighted portion is glutamine which is E. I’m not sure if I can tell the difference between wild-type and D53E. Let’s plug it into R and see if we can see differences in RMSD

library(bio3d)

covr_wt <- read.pdb("AF-D3KVK6-F1-model_v6.pdb")
covr_d53e <- read.pdb("covr_d53e_chainA.pdb")
covr_wt_idx <- atom.select(covr_wt, "calpha")
covr_d53e_idx <- atom.select(covr_d53e, "calpha")
fit <- fit.xyz(fixed = covr_wt$xyz, 
               mobile = covr_d53e$xyz, 
               fixed.inds  = covr_wt_idx$xyz,
                mobile.inds = covr_d53e_idx$xyz)

n_res <- length(covr_wt_idx$atom)
n_res_seq <- seq(1,n_res*3,3)

rmsd_vec <- vector(mode = "numeric", length=length(covr_wt_idx$atom))

for (i in 1:n_res) {
  idx <- c(n_res_seq[i]:(n_res_seq[i]+2))
  rmsd_vec[i] <- rmsd(covr_wt$xyz[covr_wt_idx$xyz][idx], fit[covr_d53e_idx$xyz][idx])
}

df <- tibble(
  residue = covr_wt$atom$resno[covr_wt_idx$atom],
  aa = covr_wt$seqres,
  rms = rmsd_vec
)

df |>
  ggplot(aes(x=residue,y=rms)) +
  geom_line() +
  geom_label(aes(label=aa),size=2) +
  theme_bw()

Alright, what this tells us is that certain residues actually didn’t vary much, but if we start seeing higher rmsd such as residue 180s, those are really high rmsd. Later on after hoddock, we’ll see if we can correlate these high RMSD residues to residues in protein-protein complex that interacted.

Next we’ll take a look at using Haddock3 for the first time and we see what we get. Our assessment is to evaluate the haddock score and body surface area (BSA) to see if a dimerization of 2 of the same covr (D53 vs D53E) would be more favorable in comparison. Since we’ve never done this before, we had Claude Code set up for us. This is mainly for my notes purposes so in the future when I review back I can reference and modify accordingly. The steps are:

1. Prepare the input files: Duplicate the CovR pdb but change the chain to A and B.
2. Set up the haddock input files: Create a text file of active and passive residue, then use haddodck3-restratins to genereate tbl file. The text file will look something like this

we can create the actpass_A.txt file with the following content:

create actpass_A.txt

87 88 89 90 91 98 99 100 101 106 107 108 109 110 111 112 113 114 115 116 117 118 120 121
81 82 83 84 85 86 93 94 102 103 104 105 119 122 123 124 125 126

With the above of alpha4-beta5-alpha5 the first row are active residues - solve-exposed residues; the 2nd row is passive residues - surface-exposed neighbors surrounding the active patch. Then we use haddock-restraint

Generate tbl file

conda run --name haddock3 haddock3-restraints active_passive_to_ambig \
      actpass_A.txt actpass_A.txt \
      --segid-one A --segid-two B \
      > ambig_covr_dimer.tbl

This basically tells haddock3 that residue 87 in chain A should interact with residue 87, 88, etc. in chain B with a distance of 2.0 Å and a weight of 1.0. which looks something like this

assign (resi 87 and segid A)
  (
         (resi 87 and segid B)
          or
         (resi 88 and segid B)
          or
         ...
  ) 2.0 2.0 0.0

3. Create a covr_dimer_workflow.toml file to specify the parameters for the docking run.

Create covr_dimer_workflow.toml

  run_dir = "run_covr_dimer"
  mode = "local"
  ncores = 4
  postprocess = true

  molecules = [
      "covr_d53e_chainA.pdb",
      "covr_d53e_chainB.pdb",
  ]

  [topoaa]
  autotoppar = false
  delenph = true

  [rigidbody]
  ambig_fname = "ambig_covr_dimer.tbl"
  sampling = 200
  sym_on = true
  nc2sym = 1
  c2sym_sta1_1 = 1
  c2sym_end1_1 = 228
  c2sym_seg1_1 = "A"
  c2sym_sta2_1 = 1
  c2sym_end2_1 = 228
  c2sym_seg2_1 = "B"

  [seletop]
  select = 100
  
  [flexref]
  ambig_fname = "ambig_covr_dimer.tbl"
  sym_on = true
  nc2sym = 1
  c2sym_sta1_1 = 1
  c2sym_end1_1 = 228
  c2sym_seg1_1 = "A"
  c2sym_sta2_1 = 1
  c2sym_end2_1 = 228
  c2sym_seg2_1 = "B"

  [emref]
  ambig_fname = "ambig_covr_dimer.tbl"
  sym_on = true
  nc2sym = 1
  c2sym_sta1_1 = 1
  c2sym_end1_1 = 228
  c2sym_seg1_1 = "A"
  c2sym_sta2_1 = 1
  c2sym_end2_1 = 228
  c2sym_seg2_1 = "B"

  [clustfcc]
  plot_matrix = true

  [seletopclusts]
  top_models = 4

  [emscoring]

The above you would have to change some of the settings and param including your tbl, pdb files. And change your ncore accordingly

3. Run Haddock: Use the command line to run Haddock with your input files. For example:

conda install - c bioconda -n haddock3 haddock3
conda activate haddock3
haddock3 -i it1.pdb -r it1.tbl -o haddock_output

What does Haddock3 do? Haddock3 is a flexible docking software that allows you to model the interaction between two or more biomolecules (like proteins) based on experimental data or predicted interactions. It uses a combination of rigid-body docking, semi-flexible refinement, and scoring to predict the most likely binding modes between the molecules. The workflow of Haddock3 goes like this rigid pose -> flexible refinement -> scoring -> clustering.

4. Look at the analysis results

result_d53 <- read_tsv("capri_ss_d53.tsv") |> arrange(caprieval_rank) |> mutate(prot="d53") |> slice_head(n=5)
result_d53e <- read_tsv("capri_ss_d53e.tsv") |> arrange(caprieval_rank) |> mutate(prot="d53e") |> slice_head(n=5)
compare_result <- rbind(result_d53,result_d53e) |>
  select(score,bsa,total,prot) |>
  pivot_longer(cols = c(score:total), names_to = "param", values_to = "values") 

compare_result |>
  ggplot(aes(x=prot,y=values,fill=prot)) +
  geom_boxplot(alpha=0.2,width=0.3) +
  geom_violin(alpha=0.5) +
  facet_wrap(.~param, scale="free_y") +
  theme_bw()

Among the top 5 ranked docking models, the phosphomimetic D53E mutant demonstrated consistently superior HADDOCK3 scores (median −138.1 vs −117.4 kcal/mol) and substantially larger buried surface area (median 2,567 vs 2,161 Ų, +18.8%), supporting enhanced dimerization propensity. Energy decomposition revealed that D53E gains its advantage primarily through van der Waals interactions (median −86.3 vs −32.7 kcal/mol, ~2.6×), despite weaker electrostatic contributions (median −335.7 vs −514.1 kcal/mol). Not really sure what this means. But let’s visualize our rank 1 pdb!

Note to self: BSA is calculated by taking the sum of the solvent-accessible surface areas (SASA) of the two individual proteins and subtracting the SASA of the complex. A larger BSA indicates a more extensive interface between the two proteins, which often correlates with stronger binding affinity. In this case, the D53E mutant’s larger BSA suggests it forms a more stable dimer compared to the wild-type D53, consistent with its role as a phosphomimetic that promotes dimerization and activation of CovR.

Let’s see which residues of these 2 chains interact!

haddock_d53e_d53e <- read.pdb("emscoring_1.pdb")

chainA <- atom.select(haddock_d53e_d53e, chain = "A")
chainB <- atom.select(haddock_d53e_d53e, chain = "B")
dist_matrix <- dist.xyz(chainA$xyz, chainB$xyz)
  
get_interface <- function(pdb, cutoff = 5.0) {
  
  chainA <- atom.select(pdb, chain = "A")
  chainB <- atom.select(pdb, chain = "B")
  
  # Get coordinates as matrices
  coordA <- matrix(pdb$xyz[chainA$xyz], ncol = 3, byrow = TRUE)
  coordB <- matrix(pdb$xyz[chainB$xyz], ncol = 3, byrow = TRUE)
  
  # Find contacting atom indices
  contact_i <- c()
  contact_j <- c()
  
  for (i in 1:nrow(coordA)) {
    dists <- sqrt(rowSums(sweep(coordB, 2, coordA[i,])^2))
    hits <- which(dists < cutoff)
    if (length(hits) > 0) {
      contact_i <- c(contact_i, rep(i, length(hits)))
      contact_j <- c(contact_j, hits)
    }
  }
  
  # Map back to residues
  resA <- pdb$atom[chainA$atom[unique(contact_i)], ] |>
    as_tibble() |>
    distinct(resno, resid) |>
    mutate(chain = "A")
  
  resB <- pdb$atom[chainB$atom[unique(contact_j)], ] |>
    as_tibble() |>
    distinct(resno, resid) |>
    mutate(chain = "B")
  
  bind_rows(resA, resB)
}

cont <- get_interface(haddock_d53e_d53e) |>
  filter(chain == "A") |> 
  select(residue=resno, aa=resid) |>
  mutate(contact = 1)

df |>
  left_join(cont, by = c("residue","aa")) |>
  mutate(contact = case_when(
    is.na(contact) ~ 0,
    TRUE ~ contact
  )) |>
  ggplot(aes(x=residue,y=rms)) +
  geom_line() +
  ggrepel::geom_label_repel(aes(label=aa,fill=as.factor(contact)),alpha=0.5,size=2) +
  theme_bw() +
  theme(legend.position = "none")

I couldn’t easily find a function in bio3d to look for closer contact between the 2 chains, hence got Claude to produce a code to filter out anything less 5 Angstrom and map it back to our RMSD df. Wow, interesting! Not all high rmsd are close contact residues, and not all close contact residues are high in rmsd.

Final Thoughts

Wow, all of the above took a long time! But it was quite interesting to learn several things here. We learnt quite a few things here. It’s our first time running a protein-protein docking! Even though I don’t deeply understand the method and the results, it’s a good start! Let’s keep at it and learn some more next time! If you notice something wrong here, please feel free to let me know!

Opportunities For Improvement

• Need to explore Prots5 and Foldseek in the future, the concept is quite interesting. Turning 3D coordinates into one dimention, very latent spacey.
• Need to understand haddock3 a bit more, the method and its output
• Need to understand how alpha helices and beta sheets occur, the math behind it and see if we can reproduce that from scratch
• Need to figure out how to reproduce a phosphorylated CovR structure instead of using a phosphomimetic
• Need to learn pymol

Lessons learnt

• learnt CovR/S two-component system in streptococcus pyogenes
• learnt what hypotheticals are, and what other methods we can use to further identify these hypotheticals
• learnt strep pyogenes is the only strep species (reference gene only) that has CovR/S, some other strep species have CovR but no CovS.
• learnt Baktfold
• learnt the bare basics of haddock3
• learnt RMSD, BSA, haddock score, angstrom unit (just eucleadian distance of xyz)

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

A previous post introduced the crossvalidation package for R. This time, the focus is on probabilistic forecasting — evaluating not just how accurate point forecasts are, but how well-calibrated prediction intervals are, using empirical coverage rates and Winkler scores – and crossvalidation.

install.packages("remotes")

install.packages("forecast")

remotes::install_github("Techtonique/crossvalidation")

library(crossvalidation)

Example 1

require(forecast)
data("AirPassengers")



eval_metric <- function(predicted, observed)
{
  error <- observed - predicted$mean

  me <- mean(error)
  rmse <- sqrt(mean(error^2))
  mae <- mean(abs(error))

  # ----- 80% interval -----

  lower80 <- predicted$lower[, 1]
  upper80 <- predicted$upper[, 1]

  coverage80 <- mean(
    observed >= lower80 & observed <= upper80
  )

  alpha80 <- 0.20

  winkler80 <- ifelse(
    observed < lower80,
    (upper80 - lower80) + (2 / alpha80) * (lower80 - observed),
    ifelse(
      observed > upper80,
      (upper80 - lower80) + (2 / alpha80) * (observed - upper80),
      (upper80 - lower80)
    )
  )

  # ----- 95% interval -----

  lower95 <- predicted$lower[, 2]
  upper95 <- predicted$upper[, 2]

  coverage95 <- mean(
    observed >= lower95 & observed <= upper95
  )

  alpha95 <- 0.05

  winkler95 <- ifelse(
    observed < lower95,
    (upper95 - lower95) + (2 / alpha95) * (lower95 - observed),
    ifelse(
      observed > upper95,
      (upper95 - lower95) + (2 / alpha95) * (observed - upper95),
      (upper95 - lower95)
    )
  )

  c(
    ME = me,
    RMSE = rmse,
    MAE = mae,
    Coverage80 = coverage80,
    Winkler80 = mean(winkler80),
    Coverage95 = coverage95,
    Winkler95 = mean(winkler95)
  )
}

(res <- crossval_ts(y=AirPassengers, initial_window = 10,
horizon = 3, fcast_func = forecast::thetaf, eval_metric = eval_metric))
print(colMeans(res))


Loading required package: forecast



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

         ME        RMSE         MAE  Coverage80   Winkler80  Coverage95 
  2.6570822  51.4271704  46.5118747   0.6590909 218.4527816   0.8459596 
  Winkler95 
312.1383104 

Example 2

eval_metric <- function(predicted, observed)
{
  error <- observed - predicted$mean

  me <- mean(error)
  rmse <- sqrt(mean(error^2))
  mae <- mean(abs(error))

  # Only one interval returned
  lower <- predicted$lower
  upper <- predicted$upper

  coverage <- mean(
    observed >= lower & observed <= upper
  )

  alpha <- 0.05

  winkler <- ifelse(
    observed < lower,
    (upper - lower) + (2 / alpha) * (lower - observed),
    ifelse(
      observed > upper,
      (upper - lower) + (2 / alpha) * (observed - upper),
      (upper - lower)
    )
  )

  c(
    ME = me,
    RMSE = rmse,
    MAE = mae,
    Coverage95 = coverage,
    Winkler95 = mean(winkler)
  )
}

fcast_func <- function(y, h, ...)
{
  forecast::thetaf(
    y,
    h = h,
    level = 95
  )
}

res <- crossval_ts(
  y = AirPassengers,
  initial_window = 10,
  horizon = 3,
  fcast_func = fcast_func,
  eval_metric = eval_metric
)

print(colMeans(res))

  |======================================================================| 100%
         ME        RMSE         MAE  Coverage95   Winkler95 
  2.6570822  51.4271704  46.5118747   0.8459596 312.1383104 

boxplot(res[, "Coverage95"])

Your tests pass. Coverage is high. Everything looks fine — until someone finds a bug in production that your didn’t catch – all because of a poor assertion.

Code coverage tells you which lines ran. It says nothing about whether those lines are actually tested. You can delete every assertion in your test suite, run covr, and still see 100%. Coverage is a measure of execution, not correctness. That gap is exactly what {muttest} was built to close — and 0.2.0 makes it much more capable than the previous version.

See the full changelog here.

What Is Mutation Testing?

Mutation testing asks a harder question than coverage: if this code were subtly wrong, would your tests notice?

It works by making small, deliberate changes to your source code — swapping > for >=, flipping TRUE to FALSE, replacing && with || — and then running your test suite against each modified version. Each modified version is called a mutant. If your tests fail, the mutant is killed: your tests noticed the change. If your tests pass, the mutant survived: your tests are blind to that kind of bug.

The result is a mutation score:

$\text{Mutation Score} = \frac{\text{Killed Mutants}}{\text{Total Mutants}} \times 100\%$

• 0% — Your tests pass no matter what the code does. Assertions are missing or trivial.
• 100% — Every mutation triggers a test failure. Your tests are tight.

Unlike coverage, this score reflects assertion quality, not just execution. A test suite full of expect_true(is.numeric(x)) checks will hit 100% coverage while missing every meaningful failure. Mutation testing exposes that.

Why You Should Care

Here is the canonical example. The function is_adult has a boundary condition:

# R/is_adult.R
is_adult <- function(age) {
  age >= 18
}

And these tests give 100% coverage:

# tests/testthat/test-is_adult.R
test_that("is_adult returns TRUE for adults", {
  expect_true(is_adult(25))
})

test_that("is_adult returns FALSE for minors", {
  expect_false(is_adult(10))
})

Both tests pass. Both would still pass if >= were accidentally replaced with >. The boundary value 18 is never tested, so neither mutant is killed:

#' R/is_adult.R — mutant 1: ">=" → ">"
is_adult <- function(age) {
  age > 18
}

Imagine this bug makes it to production. A 18 year old user tries to sign up, and the system rejects them. The bug is real, but your tests never saw it coming.

Running muttest exposes this immediately:

library(muttest)

plan <- muttest_plan(
  mutators = comparison_operators()
)
muttest(plan)

The progress table shows one survivor. The fix is a single test:

test_that("is_adult returns TRUE at the boundary age", {
  expect_true(is_adult(18))  # kills the >= → > mutant
})

This surviving mutant is not a problem to fix — it’s a specification you forgot to write.

The LLM Test Problem

Many developers now use LLMs to generate tests. Who likes to write tests themselves anyway?

LLMs are fast and produce syntactically correct code, but they may produce obvious cases, miss boundaries or just test properties of the code. The is_adult test suite above is what a language model might produce: structurally fine, semantically incomplete.

Mutation testing gives you an objective signal for how strong tests actually are, whether you wrote them yourself or they were generated by an LLM. A low mutation score doesn’t mean the LLM did a bad job — it means you now know exactly where to strengthen the assertions. LLM-generated tests need external validation just as much as human-written tests do.

muttest provides tools to help with this validation.

What’s New in 0.2.0

Expanded Mutator Library

The biggest addition in this release is a full roster of new mutators, organized into individual mutators and ready-made preset collections.

New individual mutators:

• boolean_literal("TRUE", "FALSE") — flips boolean constants: TRUE → FALSE
• na_literal("NA", "NULL") — swaps NA variants and NULL: NA → NULL
• negate_condition() — wraps if conditions in !(...): if (x > 0) → if (!(x > 0))
• remove_condition_negation() — strips leading ! from conditions: if (!done) → if (done)
• numeric_increment() / numeric_decrement() — shifts numeric constants by one: 5 → 6, 5 → 4
• index_increment() / index_decrement() — shifts subscript indices: x[i] → x[i + 1L]
• string_empty() — replaces non-empty strings with "": "hello" → ""
• string_fill() — replaces empty strings with "mutant": "" → "mutant"
• call_name("any", "all") — swaps function names: any(x) → all(x)
• remove_negation() — removes ! anywhere: !is.na(x) → is.na(x)
• replace_return_value() — replaces explicit return values with NULL: return(x) → return(NULL)
• delete_statement() — removes assignments and standalone calls one at a time, catching untested side effects and dead assignments

New preset collections — pass a single call and get the full set of relevant mutators:

• boolean_literals() — TRUE ↔ FALSE, T ↔ F
• na_literals() — NA ↔ NULL, NA ↔ NA_real_, NA ↔ NA_integer_, NA ↔ NA_character_
• numeric_literals() — combines numeric_increment() and numeric_decrement()
• index_mutations() — combines index_increment() and index_decrement()
• string_literals() — combines string_empty() and string_fill()
• condition_mutations() — combines negate_condition() and remove_condition_negation()

The three operator presets from 0.1.0 are still there — arithmetic_operators(), comparison_operators(), logical_operators() — and now they have company.

A practical starting configuration covers most of what you’d want to catch in business logic:

plan <- muttest_plan(
  source_files = "R/my_file.R",
  mutators = c(
    arithmetic_operators(),
    comparison_operators(),
    logical_operators(),
    condition_mutations(),
    numeric_literals(),
    list(remove_negation())
  )
)

Layer in boolean_literals(), na_literals(), string_literals(), or index_mutations() based on what your code actually does.

Mutators Are Now Parametrized

Individual mutators accept configuration arguments. operator("+", "-") and boolean_literal("TRUE", "FALSE") let you define exactly which token to replace and with what — so you can express the mutations that matter for your domain without writing a custom mutator from scratch. The Mutator base class is also now exported for cases where you want to go further and build an entirely custom mutator.

Survived Mutants Are Now Reported

The ProgressMutationReporter previously showed you only killed and total mutant counts. In 0.2.0, it now reports survived mutants — the ones your tests missed.

This is the signal that matters. Survivors are not noise; each one represents a real gap in your test suite. Seeing them surfaced directly in the progress output makes the feedback loop tighter: run muttest, read the survivors, add a test, repeat.

i Mutation Testing
  |   K |   S |   E |   T |   % | Mutator  | File
v |   1 |   0 |   0 |   1 | 100 | > → <    | shipping.R
x |   1 |   1 |   0 |   2 |  50 | > → >=   | shipping.R
-- Survived Mutants -----------------------------------------------
shipping.R  > → >=
2-   if (weight_kg > 5) 15.00 else 5.00
2+   if (weight_kg >= 5) 15.00 else 5.00
-- Results --------------------------------------------------------
[ KILLED 1 | SURVIVED 1 | ERRORS 0 | TOTAL 2 | SCORE 50.0% ]

Timeouts and Improved Error Handling

Mutation testing works by running your test suite once per mutant. Some mutations produce code that hangs — an infinite loop, a blocking call, a computation that never completes. In 0.1.0 that would stall your entire run.

In 0.2.0, muttest() supports per-mutant timeouts. Set a timeout and any mutant whose test run exceeds it is marked as errored. The rest of the run continues unaffected.

Error handling in general has been improved. When test execution fails unexpectedly, errors are now captured and reported cleanly rather than surfacing as unhandled conditions that stop the whole run. This makes mutation testing more robust in real projects where test environments are not always perfectly controlled.

Parallel Execution

The 0.1.0 release ran mutants sequentially. In large files with many mutants, that adds up. muttest() now supports parallel execution with {mirai} under the hood: mutants can be run concurrently across multiple workers, cutting run time on larger repositories.

Getting Started

Install from CRAN:

install.packages("muttest")

Pick one file with meaningful logic — branching, comparisons, arithmetic. Define a plan:

library(muttest)

plan <- muttest_plan(
  source_files = "R/your_file.R",
  mutators = comparison_operators()
)

muttest(plan)

Read the output. Find the survivors. Add the tests they imply. Repeat.

Start with one file and one mutator preset. Aim for a meaningful score improvement each iteration rather than chasing 100% immediately. A score of 80%+ on critical business logic is a strong starting target.

Try it on a file where you suspect the tests are weak. The survivors will tell you exactly what to add.

I’d Love to Hear From You

{muttest} is still fresh and its features and interface might change. The new mutator library covers a wide range of patterns, but there are certainly mutations specific to your domain that aren’t covered yet. If you run into a case where the right mutation is missing, an existing mutator behaves unexpectedly, or something in the output is hard to interpret, please open an issue on GitHub.

Feature requests are equally welcome. If there’s a kind of code change you’d want to test for and there’s no good way to express it yet, please drop an issue in the repository.

I came across a post recently by a machine learning engineer who made the bold claim that logistic regression is the worst name for an algorithm ever, or something along those lines¹. Many statisticians of the more old-school type seemed to disagree. This led me to think a bit more deeply about the subject. I’ve already written several posts on bad terminology in statistics (see confidence level, line of best fit, r squared) so I might have been expected to agree with the machine learning view, but in this case I agree with the statisticians, and I would like to explain why.

What data scientists think regression is

In data science classes, students are taught that there are two kinds of predictive modelling. In both cases, the aim is to predict a response $Y$ given a vector of features $X$. If $Y$ is real-valued (numeric in R terminology) then it’s a regression problem. If $Y$ is categorical then it’s a classification problem. I’m not sure where this terminology originated, but it’s certainly been propogated very widely by Hastie and Tibshirani’s classic The Elements of Statistical Learning.

In logistic regression, your data consists of some feature values $X$ and a response $Y \in \lbrace 0, 1 \rbrace$. In this case, the response is definitely categorical, so someone trained in data science would indeed call this a classification problem. But if you look more closely at the output produced by logistic regression, its predicted values are numbers, namely the probability of each data point being in the class labelled $1$. You need to do something to these numbers (for example, use a cutoff) in order to get a predicted class.

For example, in R:

set.seed(100)
N <- 100
a <- -1
b <- 1
x <- 2 * rnorm(N)

# simulated binary data
y <- rbinom(N, 1, 1/(1 + exp(-a -b * x)))

# plot observed values in grey
plot(x, y, pch=19, xlab="x", ylab="y",
     col=rgb(0, 0, 0, 0.3), las=1)

# fit logistic regression
model <- glm(y ~ x, family="binomial")

# plot predicted values in red
points(x, 
       predict(model, data.frame(x=x), 
                  type="response"),
       col=rgb(1, 0, 0, 0.3),
       pch=19)

In fact, it’s quite hard to think of a machine learning algorithm which directly predicts class membership rather than some sort of measure of how strongly a data point is a member of a class. Even Naive Bayes is making some sort of attempt to predict the probability of class membership. The simplest algorithm which directly predicts the class instead of the probability of class membership is the 1-nearest neighbour algorithm. (But if you used a larger number of neighbours, say 20, you would get some sort of estimate of how confident you were in your prediction.)

What statisticians think regression is

The term regression comes from Galton’s idea of regression to the mean (which I have written about here). Originally this was the observation that tall parents tend to have children who are shorter than them, and vice versa. The heights of children seem to regress towards the mean of the whole population.

More generally, the values of the response $Y$ corresponding to some fixed value of the features $x_0$ will follow some probability distribution. The mean of this distribution is $E[Y \vert x_0]$. The observed values of $Y$ will cluster around this mean. If you repeatedly draw values of $Y$, a large value will tend to be followed by a smaller value, and vice-versa. Thus, $E[Y \vert X]$ will tend to be smaller than $Y$ if $Y$ is unusually large, and larger than $Y$ if $Y$ is unusually small². You can see this if you use linear regression to predict $Y$ given $X$, as in the following example.

set.seed(100)
N <- 500

x <- rnorm(N)
y <- 0.4 * x + 0.8 * rnorm(N)
plot(x, y)
abline(coef(lm(y~x)), col="red")

(Note how the slope of the regression line is shallower than the “slope” which the eye perceives in the cloud of data points, which is the principal axis.)

But some algorithms don’t give you any regression effect. For example, an overfitted decision tree (a.k.a 1-NN regressor) will not show any regression to the mean, as in the following example. Note that the blue line does not under- or over-predict for the extreme values of $x$.

x <- c(1:9)
y <- c(-10, seq(-1,1, length=7), 10)
pred_nn <- function(xx) y[which.min(abs(xx - x))[1]]

plot(x, y)
abline(coef(lm(y~x)), col="red")
xx <- seq(1, 9, length=1000)
lines(xx, sapply(xx, pred_nn), type="s", lty=2, col="blue")

In this case, you have an algorithm which is predicting a numerical value, so data scientists would call it a regression, but it’s not actually exhibiting any regression. How annoying!

What regression actually is

Although it’s too late to rewrite the textbooks, maybe it could be argued that regression and classification should have been defined in the following way. If a predictive model directly predicts a response $Y$ given features $X$, then it should be called a classification model (even if $Y$ is numeric, as in the previous example). But if the model predicts $E[Y \vert X]$, then it should be called a regression model.

What about logistic regression? In this case, the model is predicting $P(Y=1 \vert X)$ which is just $E[Y \vert X]$. So the statisticians were right in the first place! Logistic regression is a regression model. It only becomes a classification model if you apply a second model to it. Usually this takes the form of a decision tree which predicts $Y=1$ if $E[Y \vert X] > p_0$ for some choice of $p_0$ and $Y=0$ otherwise. This decision tree is a classification model. But logistic regression itself isn’t.

1: I’m a little wary of calling myself a data scientist these days, partly because I think the profession has been devalued by various attempts to cash in on its popularity (leading to a glut of people with high confidence and low experience) and partly because I think data science is becoming a bit of a toxic brand with all the real-world harm being done by AI, data centres, mass surveillance, etc.

2: Anecdote time: at one of my old jobs we had to entertain a vendor who was basically selling a Kaggle-style workflow as a software-as-a-service product. The sales rep built a model on some of our data and presented it. In their write-up they included the observation that “interestingly, we noticed that the model tends to underpredict for large values of $x$ and overpredict for small values of $x$”. Well, that’s not very surprising because that’s what every predictive model does!

Digging through our memory box, we came across a conversation from which we tried to piece together when it all began with rOpenSci.

On July 13, 2011, an email was sent with the idea of a shared blog, a clever domain name, and a way to connect R package developers who cared about open science. The name “rOpenSci” appear in that email. A few months before that, the first commits had already been pushed to what would become taxize and treeBASE, two packages that quietly planted the seed of something much bigger.

That was 15 years ago. This year, we celebrate.

Quinceañera time

Fifteen years¹ is a milestone worth marking properly, that is why we want to celebrate with our community. We have a full year of activities planned, and we want you along for all of it.

Expect several diverse events, retrospectives and a few surprises we’re still stitching together. We’ll be reflecting on what we’ve built, highlighting the work of contributors old and new, and dreaming out loud about the next 15 years.

Stay tuned to this blog and our newsletter for announcements as the year unfolds.

First up: co-working session and casual virtual community celebration

We’re kicking things off with one co-working session on Tuesday, June 2 and two virtual celebrations on Wednesday, June 10 and Wednesday, June 17 in different timezones, so as many people as possible can join.

Each 90-minute celebration is built around rotating small-group conversations. You’ll meet community members you may not know yet, and together you’ll dig into questions to reflect on rOpenSci’s past and future.

Old friends and new faces alike, we would love to share this celebrations with you. No registration needed.

Thank you

To everyone who has contributed a package, reviewed code, written a blog post, helped in the forum, showed up to a community call, or simply used our tools in your research, thank you!

Here’s to 15 years and to whatever comes next.

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

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

Distribution of wars’ durations

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

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

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

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

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

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

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

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


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

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

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

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

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

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

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

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

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

My chart to answer this question is this one:

What’s happening here is:

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

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

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

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

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

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

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

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

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

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

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

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

Duration and other factors

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

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

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

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

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

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

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

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

That’s all folks. Stay safe out there.

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

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

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

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

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

The Meaning of

The coefficient of determination is defined as:

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

The denominator

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

The numerator

measures the remaining unexplained error after fitting the model.

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

An of:

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

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

The Mean Model

Let us begin with the simplest possible regression model.

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

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

In this case:

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

Substituting this into the formula gives:

The model explains none of the variation in the data.

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

The Polynomial Model

Now consider the opposite extreme.

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

In that case:

for all observations, implying:

and therefore:

The model achieves a perfect fit.

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

This is the classical overfitting problem.

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

The Linear Model

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

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

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

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

and by the famous observation attributed to George Box:

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

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

Why Alone Is Insufficient

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

It does not directly measure:

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

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

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

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

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

A high can therefore mean two very different things:

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

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

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

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

What Is Expected Goals?

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

An xG model can use variables such as:

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

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

Install and Load R Packages

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

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

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

Create a Simple Shot Dataset

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

set.seed(123)

n_shots <- 5000

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

glimpse(shots)

Engineer Shot Distance and Angle

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

goal_x <- 120
goal_y <- 40

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

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

Create a Goal Outcome

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

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

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

Explore the Shot Data

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

Visualize Shot Locations

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

Split Data into Training and Testing Sets

set.seed(123)

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

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

nrow(train_data)
nrow(test_data)

Build a Logistic Regression xG Model

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

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

summary(xg_model)

Convert Model Output into xG Values

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

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

Evaluate the xG Model

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

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

ROC AUC

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

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

Brier Score

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

brier_score

Create xG Buckets for Calibration

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

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

Player-Level xG Analysis

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

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

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

Team-Level xG Analysis

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

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

Shot Quality Distribution

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

Compare Goals and xG by Situation

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

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

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

Build a More Advanced xG Model with Interactions

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

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

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

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

Compare Two xG Models

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

model_comparison

Create a Reusable xG Prediction Function

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

new_predictions <- predict_xg(xg_model, test_data)

head(new_predictions)

Create a Custom Shot Example

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

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

Use worldfootballR for Real Football Workflows

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

library(worldfootballR)
library(tidyverse)

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

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

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

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

Example: Clean Match Results Data

clean_results <- premier_league_results %>%
  clean_names()

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

names(clean_results)

Build a Match-Level Team Summary

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

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

team_summary_example

Save Your xG Model

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

saveRDS(xg_model, "xg_model_logistic_regression.rds")

loaded_xg_model <- readRDS("xg_model_logistic_regression.rds")

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

Create an xG Report Table

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

xg_report
write_csv(xg_report, "xg_player_report.csv")

Create an xG Shot Map

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

Create a High-Value Chances Table

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

big_chances %>%
  head(20)

Model Improvement Ideas

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

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

Train an XGBoost-Style Model Later

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

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

library(tidymodels)

# A future tidymodels workflow could look like this:

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

xg_recipe

Build a Tidymodels Logistic Regression Workflow

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

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

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

xg_fit
tidy(xg_fit)

Predict Probabilities with Tidymodels

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

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

Turn xG into Match Insights

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

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

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

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

Find Teams Creating Better Chances

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

team_chance_quality

Final Thoughts

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

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

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

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

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

install.packages("unifiedml")

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

install.packages("glmnet")

library(unifiedml)

1 – Benchmarking models

set.seed(123)

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

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

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

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

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

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

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

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

$glm
$glm$avg_score
[1] 0.9533333

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


$rf
$rf$avg_score
[1] 0.96

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


$svm
$svm$avg_score
[1] 0.9733333

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

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

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

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

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

library(ggplot2)

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

2 - Unified interface for predicting probabilities

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

# Load iris dataset
data(iris)

# Setup reproducible data
set.seed(42)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Predicting probabilities for first 5 test samples:

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

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

Interpretation:

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

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

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

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

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

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

Test set accuracy: 94.7%

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

Predicting probabilities for first 5 test samples:

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

Detailed predictions:

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

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

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

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

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

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

Test set accuracy: 97.4%

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

Predicting probabilities for first 5 test samples:

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

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

Detailed predictions:

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

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

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

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

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

Test set accuracy: 94.7%

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

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

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

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

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

Introduction

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

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

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

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

Zero-crossing

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

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

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

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

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

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

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

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

Image gradients

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

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

Prewitt: weights all neighbours equally:

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

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

Canny edge detection

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

Step 1 – Gaussian smoothing

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

Step 2 – Gradient computation

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

Step 3 – Non-maximum suppression (NMS)

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

Step 4 – Hysteresis thresholding

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

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

Let’s get started with Python

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

pip install opencv-python scikit-image matplotlib numpy

Image loading and preprocessing

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

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

The Laplacian of Gaussian and zero-crossings

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

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

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

Gaussian smoothing and gradient computation

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

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

Canny edge detection with OpenCV

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

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

The effect of hysteresis thresholds

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

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

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

Overlaying edges on the original image

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

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

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

Conclusion

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

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

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

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

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

References

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

Read it in: Español.

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

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

New mentors

Alber Hamersson Sánchez Ipia

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

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

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

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

Pablo Paccioretti

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

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

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

Champions to mentors

Erick Navarro Delgado

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

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

Guadalupe Pascal

Hello! My name is Guadalupe.

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

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

Andrea Gomez Vargas

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

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

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

Monika Avila Marquez

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

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

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

Returning mentors

Luis D. Verde Arregoitia

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

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

Pao Corrales

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

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

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

Francisco Cardozo

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

Milagros Mendoza

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

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

Elio Campitelli

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

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

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

What’s next

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

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

Introduction

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

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

Why Volatility Variables Instead of Derivatives?

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

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

Why Twin Networks?

The idea is to separate the learning tasks:

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

Ensemble via Stacking

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

Evaluation

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

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

Confidence Intervals

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

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

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

Visualization

The forecasts are visualized with ggplot2:

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

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

Keras3 in R: Flexible Deep Learning for Financial Forecasting

What is Keras3?

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

How We Used Keras3

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

Why ReLU?

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

Why Adam?

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

Contribution to the R Ecosystem

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

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

Why It Matters for DML

By using Keras3:

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

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

Conclusion

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

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

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

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

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

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

# List objects in the current environment
ls()

[1] "set_function_params"

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

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

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

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

param_set_result

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

JAGS 5.0.0-beta is now available from SourceForge.

The beta release is for two groups of people:

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

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

The JAGS library

The following packages are available:

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

The rjags package

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

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