I realized this week that I won’t be hiring a GIS focused research assistant the same way as I used to.

Claude Code is just too good at GIS in R now. I guess this happened around the release of Opus 4.7. I just hadn’t really pushed the coding agent until now.

This week I downloaded a global dataset of a coastal habitat of several gigabytes, expecting to spend my free time this week poking around getting a feel for how to do my intended analysis.

30 minutes later the analysis was done.

Claude solved what used to be significant and time consuming GIS issues with no help from me. All I gave it was the end goal and reference to a previous workflow for a similar project.

It smoothly addressed complex challenges like geometry self intersections and spatial intersections of raster to vector data. It even ran small-scale tests to identify processing bottlenecks and find the fastest solution for intersecting two large layers.

Attentive package authors already noticed something new in their r-universe page over the past few weeks: R-universe has quietly started building and checking packages for Windows ARM64.

• https://ropensci.r-universe.dev/pdftools#checktable
• https://cran.r-universe.dev/dplyr#checktable
• https://bioc-release.r-universe.dev/S4Vectors#checktable

R-universe has supported ARM64 on Linux and macOS for some time already, but not Windows. With this addition, packages are now built and checked on both x86_64 and ARM64 across all major operating systems, making the R ecosystem fully ARM ready.

Porting the R ecosystem to ARM Windows

The R-universe service is part of a broader community effort to bring the R ecosystem to Windows ARM64. This project builds on the work from Tomáš Kalibera, whose efforts laid much of the groundwork for modernizing R on Windows to support UCRT and ARM.

Among other things, the project provides daily builds of base R installers for Windows ARM64 (R-release, R-patched, and R-devel), built on GitHub Actions with the rtools45 toolchain. Because CRAN repositories do not have Windows ARM64 binaries yet, these installers are patched to install binary packages from R-universe instead, so install.packages() just works. There is nothing users need to configure.

What we had to fix along the way

Rolling this out on R-universe turned out to be a great way to find out which parts of the ecosystem were not quite ready for Windows on ARM. Some things that needed patching:

• The base R installers themselves, to pick up ARM64 binaries from the right repository path.
• External software such as JAGS, which needed some changes to compile rjags and runjags for arm.
• Rust integration: many R packages using Rust were hard-coding the x86_64-pc-windows-gnu target. The extendr project shipped a fix in extendr-api 0.8.2, and rustup now supports the aarch64-pc-windows-gnullvm target we need.
• Big compiled packages like arrow, plus a few dozen other packages with issues such as unconditionally passing -msse compiler flags (ARM chips do not do SSE) or autoconf scripts misdetecting the toolchain.

For broken CRAN and Bioconductor packages we maintain a set of install-time patches, and we are in the process of upstreaming these fixes to the package maintainers. With these patches in place, the entire R ecosystem now builds on Windows ARM64.

Testing in the real world

The fact that everything builds does not yet guarantee that everything behaves exactly as expected. We would love help from users to find any loose ends. Hence if you are running Windows on ARM hardware, or you are using a Windows VM on MacOS, we would love for you to take this for a spin.

For package developers there several ways to test packages with compiled code on Windows ARM64. As mentioned, if your package is already on CRAN or BioConductor, you can find binaries and check results on R-universe at https://cran.r-universe.dev/{pkg} or https://bioc.r-universe.dev/{pkg}.

In addition, GitHub Actions provides free windows-11-arm runners, and the standard workflows from r-lib/actions support them out of the box.

Helping the ecosystem evolve

This project is exactly one of the reasons for building the R-universe system. Besides providing an open space for developers to publish new software, it also serves as a versatile infrastructure that lets us experiment with new tooling and platforms to help move the ecosystem forward as a whole. Having an automated platform that continuously builds thousands of packages made it possible to find those missing pieces, fix them, and iterate quickly.

By identifying and fixing the main hurdles now, we want to help package authors address the problems before they hit end-users, and pave the way for a smooth introduction of Windows ARM64 support on mainstream repositories such as CRAN and Posit Package Manager when the time is ready.

FADE IN

Nine students in the M.S. in Psychology Research program at Brooklyn College, New York City (NYC), USA are taking the course Reproducible Psychological Research, where they are learning how to code for the first time. All of them have their cameras off. They are all learning R for the first time, and as an introduction to R, are being taught using the mtcars dataset.

Nine NYC graduate students who

• don’t own a car,
• don’t drive a car,
• and have no interest in cars.

Even with their screens off, with every additional mtcars column we talked about, they fell deeper and deeper asleep during the lectures.

If I wanted my students to learn anything in my class—not just R—I needed to make an immediate pivot.

Change needed to happen, and happen fast.

Relevance drives rigor

I didn’t take my students’ interests into consideration at the start of the course, and I was paying the price. How do I make sure that no matter what each individual student is interested in, they can be accommodated? Probably a big ask, but part of me felt like I was starting the semester in debt and needed to really pay up.

What do all my students have in common? They all live in NYC. Even if they’re not from NYC, they currently call it home.

Then came the idea:

What if I incorporate NYC Open Data into the course?

For anyone who doesn’t know, NYC does an incredible job of making public data available through the NYC Open Data Portal. Thousands of datasets, ranging from 311 calls to transportation data to slashings and stabbings (a favorite of one of my students), are available for free, for anyone, anywhere.

Immediately, there was a switch in my (virtual) classroom. Before, engagement was about as close to zero as you can get. Now, it skyrocketed. Students were asking more questions, which was great.

There were really only two options when it came to using NYC Open Data. The first, which I began using, was downloading static datasets as Excel or CSV files. While this helped with engagement, it also created a new set of problems: using outdated datasets (311, for instance, is updated daily), making sure students downloaded the correct files, and making sure they put them in the correct folders.

Luckily, the portal offers a fantastic solution to all of the problems associated with accessing static files: APIs!

However, this came with an even bigger problem:

In this class, I had to teach R, statistics, reproducible workflows, data analysis, data cleaning, data visualization, writing, interpretation, and really an innumerable number of additional skills. Did I really want to add APIs to this list?

Cue an idea.

nycOpenData: The beginning

At this point, I knew a few things:

1. Students need relevant data in order to stay engaged.
2. Students really enjoyed using NYC Open Data.
3. My students’ heads may explode if I teach them about APIs.

So I started thinking:

How do I incorporate NYC Open Data into the classroom moving forward without using APIs or static files?

What if I make a package?

Students are learning R, which means we’re also learning about packages—a crucial, yet simple, concept. What if there was a package dedicated to pulling NYC Open Data? That would fit right into the students’ learning path.

When I looked, there wasn’t one.

So I thought:

Why not make it myself?

Maybe a crazy idea. I had never developed an R package, never mind one that relied on APIs. What does the process look like? Do I have the skills? Will CRAN accept something that most likely only my nine students and I would ever use?

As a big believer in the mantra, “I’ll do it myself,” I already had the code to pull in the individual datasets in R, so… why not?

So, I went ahead and did it myself.

nycOpenData: The middle

Originally, I thought the easiest thing to do was build the package with one function: nyc_311(), which would simply pull in the 311 dataset. The perk was that it was one of the largest and most data-diverse datasets, so it could kill many birds with one stone.

But then I thought,

“No way will CRAN accept a package with just one function.”

So, I literally copied and pasted the same exact code 30 different times and adapted it for the 30 most popular datasets on the NYC Open Data Portal.

Through a lot of blood, sweat, and documentation, the first iteration of nycOpenData was officially created and accepted onto CRAN! I, a person who only installed packages, was now the creator of one? It was amazing to actually give back to the marvelous R community that provided me with so much for so many years.

This not only made teaching with NYC Open Data easier, but more fun too. The package did exactly what it needed to do: make working with NYC Open Data datasets incredibly easy in R.

With this version of nycOpenData, students were able to complete their final projects with one simple prompt:

Answer any question you have about NYC using open data.

Those projects eventually became the NYC Open Data Student Gallery Book, which was presented at NYC Open Data Week 2026. Each student’s project became its own chapter, and each student had the opportunity to present their work to the world. These same students who were once falling asleep while I basically sang them a lullaby, now were very much awake, displaying what they’ve learned and uncovered to an audience of professionals.

As an educator, my goal is to help my students be as successful as possible. Through nycOpenData, that goal became a reality.

All of the work I put into this package was worth it.

nycOpenData: Student contributions

The following semester, I had the opportunity to teach a follow-up class to the same cohort. Then came another crazy idea:

What if the same students who used nycOpenData to not only learn, but also conduct research, contributed to the same infrastructure they utilized?

Basically, what if I helped my students contribute to the nycOpenData package?

This would be similar to what a junior software developer might do. Using my code as a template, they could build their own functions for whatever dataset they were interested in. The package already had 30 functions—what’s another nine?

Through light introductions to GitHub, package development, and open-source software, each of my students successfully contributed their own function to the nycOpenData package. Suddenly, they weren’t just users of open-source software—they were contributors.

Good thing I didn’t do this last semester, as that would have guaranteed brain explosions from my students.

nycOpenData: rOpenSci

As time progressed, and as I assessed the student contributions, I became more familiar with package development. Wanting the best product possible, I decided to submit the package to rOpenSci through its public peer review process. I was a first-time package developer and wanted guidance from people who had far more experience than I did. I wanted anyone—whether my students or someone using the package halfway around the world—to have the best experience possible. After all, the goal was to make NYC Open Data more accessible for anyone. To go from helping just my students, to my students helping, to potentially helping the overall R community through rOpenSci—I would have never guessed this trajectory.

This ended up being one of the best experiences I’ve had as a programmer.

The review process didn’t simply determine whether the package was ready—it made the package better. Every suggestion strengthened the documentation, improved the user experience, and taught me something new about developing open-source software.

First of all, my editor and reviewers told me exactly what I needed to hear:

Having 39 functions that basically do the same thing is not only unmanageable, but bad code.

Which was 100% the truth.

And, if I’m being honest, the package itself came about while I was teaching a course on reproducibility. I was literally teaching my students that copy-and-paste was not something to do, while I was doing the exact same thing. Sure, I brushed it off because I needed to put something together quickly, but it was time to face the facts and fix it.

With that came good news and bad news. As for every action, there is an equal and opposite reaction:

• The good news: I was able to turn the 39 functions into just 3. Now, all 2,000+ datasets can be pulled from the portal. Easier to maintain, easier for users to understand, and overall a much better product.
• The bad news: All of my students’ work was now gone. (Sad face.)

Thanks to the thoughtful feedback and hard work of the editor, reviewers, and everyone involved in the review process, nycOpenData was successfully accepted into rOpenSci—an indescribable privilege. I honestly can’t believe the journey from where the package started to where it is today.

I am especially grateful to my editor, Ronny Hernandez Mora, and reviewers, Haolin Dong and Michael Pascale, for the time, care, and thoughtful feedback they invested throughout the review. As a first-time package developer, I couldn’t have asked for a better introduction to open-source software development. With the foundation finally built, another thought entered my mind:

Can we use this same infrastructure for other open data portals?

Turns out, you can!

The OpenData Ecosystem

With the infrastructure in place, another realization hit me: I never had to build it from scratch again.

The same foundation that powered nycOpenData could be adapted to other open data portals, making it possible to create an entire ecosystem of OpenData packages.

If students in NYC were motivated by using data regarding their own city, why wouldn’t the same hold true for students in other cities?

Since then, four additional packages have been built, turning one package into an ecosystem:

• nysOpenData — New York State Open Data
• mtaOpenData — Metropolitan Transportation Authority Open Data
• laOpenData — Los Angeles Open Data
• chiOpenData — Chicago Open Data
• ausOpenData — Austin Open Data

The ecosystem continues to grow. A former student and I recently collaborated on the ausOpenData package (which was just recently accepted to CRAN), and a current intern at the NYC Open Data Lab is building sfOpenData and SeattleOpenData. This not only gives others an opportunity to contribute, but also gives me a chance to continue involving students in the ecosystem after the original nycOpenData functions were “refactored.”

Hopefully, there are many more to come. The hard part is no longer building packages—it is deciding which open data portal comes next.

nycOpenData: The Final Product

What started as a single function has turned into six OpenData packages, one of which has been successfully peer-reviewed by rOpenSci, with the others built on the same foundation. Tens of billions of rows of public data, from across the country, are now easily accessible in R for anyone—from researchers and journalists to students and educators.

Looking back, I think this journey embodies the spirit of R.

One of the best things about R is that it is open source. Anyone with an idea can contribute, and those ideas have a way of growing into something much bigger than originally imagined. This entire ecosystem started because nine graduate students didn’t care about the mtcars dataset.

On top of that, the tools used to build this ecosystem were free and open source, and the packages themselves are freely available to anyone. In my family we say, “If it’s for free, it’s for me.”

What started in one virtual classroom ended up bringing together students, educators, reviewers, professionals, interns, and thousands of people who have since installed the packages.

Through this entire process, the biggest lesson learned was:

Relevance isn’t a reward for learning—it is what enables learning.

To my students

While this story has my name attached to it, it wouldn’t exist without students who were willing to try something new, contribute code, ask questions, and occasionally tolerate my crazy ideas. When we started the course, I asked you all to keep an open mind and try something new-and I never would ask you to do something that I wouldn’t.

Thank you for trusting me when I asked you to do something different.

FADE OUT

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

SAS is almost 50 years old. It predates the web, and it predates most of the open source tools data teams reach for today. While old, it is not obsolete. SAS remains central to clinical research, banking, insurance and other regulated industries, and those organisations do not keep it because it is new. They keep it because it is trusted, governed and built into critical processes.

So its age is beside the point. What matters is whether your organisation would benefit from open source for critical processes and whether it can be trusted. Some teams are already moving off it, and the useful question today is: Can regulated industries trust open source, and which of your workloads should adopt it?

Why SAS lasted

Software survives when organisations trust it, and that matters most where mistakes are expensive.

Clinical research is a clear example. The US FDA defines the formats it accepts for electronic study data, and those standards include SAS transport files alongside CDISC data models. The FDA does not require sponsors to run their analysis in SAS; the requirement applies to the format used to submit data. But submission standards shape working practices. Over time, organisations have built SAS-based processes around these requirements, including standard outputs, controlled datasets, audit trails and reproducible analysis. The FDA also asks sponsors to provide the software programs used for certain efficacy analyses, so reviewers can reproduce the results. In that setting, reproducibility sits at the heart of the process.

Banking and insurance follow the same pattern. Credit risk models and regulatory reports must be explained, tested and defended, and SAS became established because it supported reliable, governed workflows. That is why SAS lasted: it works in environments where trust matters more than novelty.

Why organisations keep SAS too long

Some organisations have strong reasons to keep SAS. Others keep it because replacing it looks harder than maintaining it.

Years of logic sit inside the code

A mature SAS system contains more than programs. It contains macros, formats, reporting rules and business decisions built up over decades. Some of those rules only exist in the code, the people who wrote them may have left, and the documentation may be incomplete. This creates grown logic: systems built gradually by different people, without one clear design.

Changing the system feels risky because nobody fully understands it. But leaving it unchanged increases the risk over time, as the system becomes more deeply embedded and fewer people retain the knowledge needed to maintain or replace it.

SAS is only a component

SAS is often connected to reporting pipelines, validation processes, internal tools and downstream systems, so removing one part can affect the rest. Banking teams have described cases where legacy systems are closely tied to daily operations, and replacing them requires more than rewriting code. It requires understanding how the whole process works.

Risks – Skills are getting harder to maintain

Organisations need people who understand the existing SAS system, and they also need people who can design and validate its replacement. These are often the same people. That creates a difficult position: the team needs specialist knowledge to maintain SAS and the same knowledge to move away from it. Waiting has been classed as the safest option by some, but waiting increases the knowledge gap and increases the risk.

The vast majority of UK Universities haven’t taught SAS to undergraduate statistics students in over fifteen years!

Migration to open source is a full-stack problem that needs a full-stack solution

R and Python can perform the analysis, modelling and reporting that regulated organisations need. The main risk in a SAS to R migration lies in the process around the language; R itself is rarely the problem.

Regulated teams must prove that the new process works, and similar results are not enough. Outputs may need to match at row, value and format level, and teams must retain evidence showing how each comparison was completed. Undocumented macros create another risk: a macro may hold an important business rule that exists nowhere else, and a direct conversion can reproduce the code while losing the reasoning behind it. Migration means rebuilding the process itself, not simply translating code line by line.

A proven method for introducing open source: phased migration

Replacing the full SAS system at once creates unnecessary risk. A safer migration runs the old and new systems in parallel, so teams can compare outputs, test business rules, automate regression checks, document differences, maintain a rollback route, and validate each stage before release.

This takes longer than converting the code, but it produces a system the organisation can trust. A successful SAS to R migration is phased, governed and validated, with correctness and auditability built in from the start.

Is SAS worth keeping?

Sometimes SAS still has a clear role in regulated environments, supporting established processes, submission standards and governance models. But not every SAS workload is essential. Some remain because they are critical, and others remain because nobody has reviewed them.

Occasionally, the answer is both SAS and open source, not a choice between the two. Teams can use R or Python for new analysis, modelling and visualisation while keeping SAS for existing governance, reporting or submission processes, which lets organisations modernise without introducing unnecessary risk.

The first step is assessment. Before moving anything, identify which workloads are regulator-facing, which processes contain undocumented logic, which systems depend on SAS outputs, which workloads can move now, and what evidence will be required to validate the change. That gives the organisation a clear migration path.

How Jumping Rivers helps you achieve Open Source

We strategise for teams and systems from SAS to full implementation to R, treating migration as a controlled transformation rather than a code conversion exercise. That means understanding the existing system, identifying risk, validating outputs and moving workloads in stages.

If your organisation is deciding which SAS workloads to keep and which to move, get in touch. If you are still assessing whether R is the right choice, our Why Use R? course is designed for teams moving from SAS and SPSS. And if your organisation already uses R, Litmus provides audit-ready risk assessments for open source packages at package and function level.

Migration helps you adopt R. Validation helps you trust it.

Frequently asked questions

Is SAS still used in 2026?

Yes. SAS remains widely used in pharmaceutical research, banking, insurance and other regulated industries. Its value comes from governance, auditability and its place within established workflows.

Does the FDA require SAS?

No. The FDA accepts SAS transport files as a submission format, but this does not require sponsors to perform their analysis in SAS. These standards have, however, influenced the tools and processes used for clinical submissions.

Can R replace SAS in clinical trials?

Yes, in many workflows. The main challenge is validation: teams must show that R outputs are reproducible, correct and suitable for regulatory review. Many organisations use SAS and R together during the transition.

What are the main risks of migrating from SAS to R?

The main risks include undocumented business rules, incomplete knowledge of existing systems, differences between old and new outputs, and replacing too much at once. A phased migration with parallel testing and formal validation reduces these risks.

Can SAS and R be used together?

Yes. R can support new analysis, modelling and visualisation while SAS continues to support existing reporting, governance or submission processes. For many organisations, this is the safest route to modernisation.

References

• Techzine, 47 years of SAS: advanced age gives SAS an edge in the current AI landscape
• US FDA, Study Data Standards Resources
• US FDA, CDER Study Data Standards Research and Development
• US FDA, Study Data Technical Conformance Guide
• Banking.Vision, The challenge of legacy transformation from a software development perspective
• SDK.finance, How to Modernize Banking Legacy Systems via Cloud Migration
• pharmaverse blog, Faster Clinical Trial Reporting: Implementing CDISC SDTM and ADaM Standards with Open Source R Packages
• pharmaverse blog, Testing Containers and WebAssembly in Submissions to the FDA

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

Textbooks say LOOCV has the lowest bias but highest variance compared to 10 and 5-fold. Coded a K-Fold CV from scratch for learning to test that on simulated data — and at n=1000 it holds up. At n=100? Not so much.

The above image was generated via chatGPT. Uploaded all the text of this blog post and asked it to generate a cartoon. Very impressive! It used to be spelling error and gibberish of text in the past, but now cohesive words on image. Just wow.

Motivations

Crossvalidation is such a crucial step in Machine Learning (and traditional methods) that nowadays is incorporated in easy to use sklearn or tidymodels without us needing to build one from scratch. As with my other learning experience, the best way to learn the concept (other than learning the concept ) is to code it from the ground up and see how it works! In K-Fold CV, the training data is split into K chunks; the model is trained K times, each time holding out a different chunk. Performance is averaged across all K folds, giving a more stable estimate. A special case is Leave-One-Out CV (LOOCV), where each individual observation serves as its own validation set. It’s thorough but computationally expensive. I was told that, bias LOOCV < 10-fold < 5-fold; whereas variance LOOCV > 10-fold > 5-fold. Is that true? Also, what’s with the repeats, does that really reduce variance? Let’s check them out.

Objectives:

• Simulate data with a known data-generating process
• Implement K-Fold cross-validation from scratch
• Assessing RMSE
• Compare candidate models using CV RMSE
• Verify the best model on a held-out test set
• Opportunities For Improvement
• Lessons Learnt

Simulate Data

library(tidyverse)

set.seed(1)
n <- 1000
x <- rnorm(n)
w <- rnorm(n)
y <- 0.5*x^2 + -0.5*w + 0.3*w*x + rnorm(n)
df <- tibble(x,y,w)
idx <- sample(1:n, size=0.8*n)
train <- df[idx, ]
test <- df[-idx, ]

The above code simulates a dataset with 1000 observations, where the response variable y is generated based on a known data-generating process involving predictors x and w. The dataset is then split into a training set (80%) and a test set (20%). Let’s visualize.

df |>
  mutate(w_cut = cut_interval(w, n=5)) |>
  ggplot(aes(x=x, y=y, color=w_cut, group=w_cut)) +
  geom_point(alpha=0.5) +
  theme_bw() +
  geom_smooth(method = "gam", se=F)

Wow, very interesting visualization where the relationships are definitely not linear here. It’s some form of interaction between x and w. Let’s see if we can recover the underlying data-generating process using K-Fold Cross-Validation.

K-Fold Cross-Validation From Scratch

folds <- 5
segment_portion <- nrow(train)/folds
formula_list <- list(as.formula("y~x"),as.formula("y~I(x^2)"),as.formula("y~I(x^2)+w+w:x"),as.formula("y~I(x^3)+w+w:x"),
                     as.formula("y~w:x"),as.formula("y~w"),as.formula("y~x+w+x:w"),as.formula("y~I(x^2)+w:x"),
                     as.formula("y~I(x^2)+w"))

cv_log <- tibble()

for (formula in formula_list) {
print(formula)
predict_log <- y_log <- vector(mode="numeric",length=segment_portion*folds)
start <- 1
end <- segment_portion 

for (fold in 1:folds) {
    val_i <- train[start:end,]
    train_i <- train[-c(start:end),]
    model_i <- lm(formula,train_i)
    predict_i <- predict(model_i, val_i)
    predict_log[start:end] <- predict_i
    y_log[start:end] <- val_i$y
    start <- end + 1
    end <- start + segment_portion - 1
}

val_df <- tibble(predict=predict_log,y=y_log) |>
  mutate(formula=deparse(formula))
cv_log <- cv_log |>
  bind_rows(val_df)
}

## y ~ x
## y ~ I(x^2)
## y ~ I(x^2) + w + w:x
## y ~ I(x^3) + w + w:x
## y ~ w:x
## y ~ w
## y ~ x + w + x:w
## y ~ I(x^2) + w:x
## y ~ I(x^2) + w

Alright, what we’ve done above is a manual implementation of K-Fold Cross-Validation. We loop through each formula in our list, and for each formula, we split the training data into 5 folds. For each fold, we train the model on the other 4 folds and validate it on the current fold. We store the predictions and actual values for later evaluation.

We basically want to see which formula has the lowest RMSE across the folds. Let’s calculate that next. From the DGP formula, we know that the best model should be y~I(x^2)+w+w:x. Let’s see if we can recover that using K-Fold CV.

Assessing RMSE

cv_log |>
  group_by(formula) |>
  summarize(rmse = sqrt(mean((y-predict)^2))) |>
  arrange(rmse) |>
  mutate(rmse = format(rmse, digits = 8))

## # A tibble: 9 × 2
##   formula              rmse     
##   <chr>                <chr>    
## 1 y ~ I(x^2) + w + w:x 1.0397338
## 2 y ~ I(x^2) + w       1.1008693
## 3 y ~ I(x^2) + w:x     1.1607225
## 4 y ~ I(x^2)           1.2144857
## 5 y ~ x + w + x:w      1.2826723
## 6 y ~ I(x^3) + w + w:x 1.2912046
## 7 y ~ w                1.3578304
## 8 y ~ w:x              1.3732478
## 9 y ~ x                1.4451807

Here our loss function is RMSE since y is a continuous data and we’re trying to predict that. The formula with the lowest RMSE is indeed y~I(x^2)+w+w:x, which matches the underlying data-generating process. OK at least, right now we are able to recover the underlying DGP using 5-Fold Cross-Validation. But is there a difference between 5 fold, 10 fold, or even LOOCV? If there is a difference, how do we even assess that? In the past we were able to assess bias and variance based on a true ATE, but what on earth is a true RMSE !?!

To check whether the textbook claim (bias LOOCV < 10-fold < 5-fold; variance LOOCV > 10-fold > 5-fold) holds up, we ran a small simulation with help from Claude Sonnet 5. Since we control the data-generating process, we can compare the “correct formula” (assuming the correct formula has the lowest RMSE as above) with 500 different simulated dataset against a “true” RMSE estimated from a large test set (n=10000) — large enough, by the law of large numbers, to treat as ground truth. Averaging across simulations gives bias (how far off CV runs from the true error) and variance (how much CV’s estimate swings from sample to sample) for each method. Is this legit? If the textbook claim is correct, we should be able to observe bias LOOCV < 10-fold < 5-fold; variance LOOCV > 10-fold > 5-fold. Let’s see if we can observe that in the simulation below.

Compare Candidate Models

# set.seed(1)

# k-fold CV RMSE for a given formula and dataset (k = n gives LOOCV)
cv_rmse <- function(data, formula, k) {
  n <- nrow(data)
  folds <- sample(rep(1:k, length.out = n))
  preds <- numeric(n)
  for (i in 1:k) {
    train_i <- data[folds != i, ]
    val_i   <- data[folds == i, ]
    model_i <- lm(formula, train_i)
    preds[folds == i] <- predict(model_i, val_i)
  }
  sqrt(mean((data$y - preds)^2))
}

# "true" RMSE: fit on train, evaluate on a large fresh draw from the DGP
true_rmse <- function(train, formula, n_test = 10000) {
  x <- rnorm(n_test); w <- rnorm(n_test)
  y <- 0.5*x^2 - 0.5*w + 0.3*w*x + rnorm(n_test)
  test <- tibble(x, y, w)
  model <- lm(formula, train)
  sqrt(mean((test$y - predict(model, test))^2))
}

formula_true <- as.formula("y ~ I(x^2) + w + w:x")
n_sim   <- 500
n_train <- 100

results <- vector("list", n_sim)

for (s in 1:n_sim) {
  x <- rnorm(n_train); w <- rnorm(n_train)
  y <- 0.5*x^2 - 0.5*w + 0.3*w*x + rnorm(n_train)
  train_s <- tibble(x, y, w)

  results[[s]] <- tibble(
    sim      = s,
    true_err = true_rmse(train_s, formula_true),
    loocv    = cv_rmse(train_s, formula_true, k = n_train),
    cv5      = cv_rmse(train_s, formula_true, k = 5),
    cv10     = cv_rmse(train_s, formula_true, k = 10)
  )
}

sim_df <- bind_rows(results)

sim_long <- sim_df |>
  pivot_longer(cols = c(loocv, cv5, cv10), names_to = "method", values_to = "cv_estimate")

sim_long |>
  group_by(method) |>
  summarize(
    mean_cv_estimate = mean(cv_estimate),
    mean_true_error   = mean(true_err),
    bias     = mean(cv_estimate - true_err),
    variance = var(cv_estimate),
    .groups  = "drop"
  ) |>
  arrange(bias) |>
  mutate(variance = format(variance, digit = 8))

## # A tibble: 3 × 5
##   method mean_cv_estimate mean_true_error     bias variance     
##   <chr>             <dbl>           <dbl>    <dbl> <chr>        
## 1 loocv              1.00            1.00 0.000703 0.00053227867
## 2 cv10               1.00            1.00 0.000945 0.00053192050
## 3 cv5                1.00            1.00 0.00124  0.00053176441

Wow, looking at the results, we can see that the textbook claim holds up. LOOCV has the lowest bias, but the highest variance. 10-fold CV is in between, and 5-fold CV has the highest bias but lowest variance. But, noticed that we had to increase our digit to 8 to see the difference in variance, it’s really miniscule. Mainly because our n=1000 is already quite large, so the variance is already quite small. If we reduce n to 100, how would that look?

## # A tibble: 3 × 5
##   method mean_cv_estimate mean_true_error     bias variance    
##   <chr>             <dbl>           <dbl>    <dbl> <chr>       
## 1 loocv              1.02            1.02 -0.00400 0.0063905148
## 2 cv10               1.02            1.02 -0.00138 0.0065753049
## 3 cv5                1.03            1.02  0.00244 0.0068098389

!?!?!?!?! with n=100, the bias and variance order heuristics no longer hold up? Why is this? I don’t know. If you do, please let me know. I even increased the n_sim to 1000, but still same patern. Intersting how I had to push the n up to 850 in order to observe the textbook variance order again, though again it’s quite miniscule.

Verify On Test Set

test |>
  mutate(predict = predict(lm(as.formula("y ~ I(x^2) + w + w:x"), train), test)) |>
  mutate(residual = y-predict) |>
  mutate(res_square = residual^2) |>
  pull(res_square) |>
  mean() |>
  sqrt()

## [1] 1.018281

Alright! The RMSE on test set is quite similar to our average validation sets! Let’s visualize the predicted and actual y of the model on test set.

test |>
  mutate(predict = predict(lm(as.formula("y ~ I(x^2) + w + w:x"), train), test)) |>
  ggplot(aes(x=predict, y=y)) +
  geom_point(alpha=0.5) +
  theme_bw() +
  geom_smooth(method = "lm") +
  labs(title="Predicted vs Actual y on Test Set", x="Predicted y", y="Actual y")

summary(lm(y ~ predict, data=test |> mutate(predict = predict(lm(as.formula("y ~ I(x^2) + w + w:x"), train), test))))

## 
## Call:
## lm(formula = y ~ predict, data = mutate(test, predict = predict(lm(as.formula("y ~ I(x^2) + w + w:x"), 
##     train), test)))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.98185 -0.61410 -0.05514  0.62883  2.28837 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.08581    0.08390   1.023    0.308    
## predict      0.80885    0.07562  10.696   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.007 on 198 degrees of freedom
## Multiple R-squared:  0.3662,	Adjusted R-squared:  0.363 
## F-statistic: 114.4 on 1 and 198 DF,  p-value: < 2.2e-16

Slope here is 0.81, not 1 — even with the correct formula, coefficients are still estimated from a finite, noisy sample, so predict is really “true signal + estimation error.” Regressing y on a noisy version of itself pulls the slope below 1 — a known effect called attenuation (same idea as errors-in-variables bias). More training data shrinks that estimation error, so the slope should creep closer to 1. Intercept stays near 0 because the estimation error isn’t systematically biased in one direction — just noisy — so predictions still center correctly on average. If this is true, then if we increase our n to 10000, we should see slope is closer to 1 and intercept closer to 0. Let’s test this theory out.

set.seed(1)
n <- 10000
x <- rnorm(n)
w <- rnorm(n)
y <- 0.5*x^2 + -0.5*w + 0.3*w*x + rnorm(n)
df <- tibble(x,y,w)
idx <- sample(1:n, size=0.8*n)
train <- df[idx, ]
test <- df[-idx, ]
formula <- as.formula("y ~ I(x^2) + w + w:x")
folds <- 10
segment_portion <- nrow(train)/folds
predict_log <- y_log <- vector(mode="numeric",length=segment_portion*folds)
start <- 1
end <- segment_portion 

for (fold in 1:folds) {
    val_i <- train[start:end,]
    train_i <- train[-c(start:end),]
    model_i <- lm(formula,train_i)
    predict_i <- predict(model_i, val_i)
    predict_log[start:end] <- predict_i
    y_log[start:end] <- val_i$y
    start <- end + 1
    end <- start + segment_portion - 1
}

val_df <- tibble(predict=predict_log,y=y_log) 

val_df |>
  summarize(rmse = sqrt(mean((y-predict)^2))) |>
  arrange(rmse) |>
  mutate(rmse = format(rmse, digits = 8))

## # A tibble: 1 × 1
##   rmse     
##   <chr>    
## 1 1.0126176

test |>
  mutate(predict = predict(lm(as.formula("y ~ I(x^2) + w + w:x"), train), test)) |>
  mutate(residual = y-predict) |>
  mutate(res_square = residual^2) |>
  pull(res_square) |>
  mean() |>
  sqrt()

## [1] 0.9884315

test |>
  mutate(predict = predict(lm(as.formula("y ~ I(x^2) + w + w:x"), train), test)) |>
  ggplot(aes(x=predict, y=y)) +
  geom_point(alpha=0.5) +
  theme_bw() +
  geom_smooth(method = "lm") +
  labs(title="Predicted vs Actual y on Test Set", x="Predicted y", y="Actual y")

summary(lm(y ~ predict, data=test |> mutate(predict = predict(lm(as.formula("y ~ I(x^2) + w + w:x"), train), test))))

## 
## Call:
## lm(formula = y ~ predict, data = mutate(test, predict = predict(lm(as.formula("y ~ I(x^2) + w + w:x"), 
##     train), test)))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -3.08711 -0.67088  0.01239  0.69015  2.86620 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.02778    0.02527   1.099    0.272    
## predict      0.96577    0.02369  40.762   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9884 on 1998 degrees of freedom
## Multiple R-squared:  0.454,	Adjusted R-squared:  0.4538 
## F-statistic:  1662 on 1 and 1998 DF,  p-value: < 2.2e-16

There you have it! The RMSE on test set is quite similar to our average validation sets! The slope is now 0.97, much closer to 1, and the intercept is 0.02, much closer to 0. This confirms our theory that with more training data, the estimation error decreases, leading to better predictions.

Opportunities For Improvement

• apply crossvalidation from scratch on hyperparameter tuning such as xgboost hyperparam, lasso/ridge/elasticnet/glmnet
• learn about nested resampling that uses an additional layer of resampling that separates the tuning activities from the process used to estimate the efficacy of the model. see here
• learn to code grouped and stratified cv from scratch
• learn a bit more on time-series CV
• need to try n_sim of 10000 and see if the bias and variance order holds up for n=100, since the simulation is quite slow, I didn’t run it to completion. Update: I did complete an n_sim of 10000 of n=100 and LOOCV still had the lowest variance.

Lessons learnt

• Learnt that we have to use I to indicate that we want to include polynomial terms in the formula, otherwise R will interpret them as interaction terms.
• Learnt deparse is useful to convert a formula object into a character string for logging purposes.

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

Graphical user interfaces for the R language are easy to use and getting more powerful all the time. Here is my updated comparison of jamovi, JASP, BlueSky Statistics (free & Pro), Rattle, RKWard, R-Instat, R AnalyticFlow, and R Commander.

With so many detailed reviews of Graphical User Interfaces (GUIs) for R available, which should you choose? It’s not too difficult to rate them based on the number of features they offer, so I’ll start there. Then, I’ll follow with a brief overview of each.

I’m basing the counts on the number of dialog boxes in each category of the following categories:

• Ease of Use
• General Usability
• Graphics
• Analytics
• Reproducibility

This data is trickier to collect than you might think. Some software has fewer menu choices, depending instead on more detailed dialog boxes. Studying every menu and dialog box is very time-consuming, but that is what I’ve tried to do to keep this comparison trustworthy.

Perhaps the biggest flaw in this methodology is that every feature adds only one point to each GUI’s total score. I encourage you to download the full dataset and consider which features are most important to you. If you decide to make your own graphs with a different weighting system, I’d love to hear from you in the comments below. Thomas Langkamp did an excellent comparison focused on teaching statistics. You can read it here.

Ease of Use

For ease of use, I’ve defined it primarily by how well each GUI meets its primary goal: avoiding code. They get one point for each of the following abilities, which include being able to install, start, and use the GUI to its maximum effect, including publication-quality output, without knowing anything about the R language itself. Figure 1 shows the result. R Commander is abbreviated Rcmdr, and R AnalyticFlow is abbreviated RAF. The commercial BlueSky Pro comes out on top by a slim margin, followed closely by JASP and RKWard. None of the GUIs achieved the highest possible score of 15, so there is room for improvement.

The ease of use features are:

• Installs without the use of R
• Starts without the use of R
• Remembers recent files
• Hides R code by default
• Uses its full capability without needing R code
• Data editor included
• Pub-quality tables w/out R code steps
• Simple menus that grow as needed
• Table of Contents to ease navigation
• Variable labels ease identification in the output
• Easy to move blocks of output
• Ease reading columns by freezing headers of long tables
• Accepts data pasted from the clipboard
• Easy to move header row of pasted data into the variable name field
• Includes AI interface to ease interpretation of output

General Usability

This category is dominated by data-wrangling capabilities, which data scientists and statisticians spend most of their time on. It also includes various types of data input and output. We see in Figure 2 that both BlueSky versions and R-Instat come out on top, not just due to their excellent selection of data-wrangling features but also for their use of the rio package for importing and exporting files. The rio package combines the import/export capabilities of many other packages, and it is easy to use. I expect the other GUIs will eventually adopt it, raising their scores by around 20 points.

The general usability features are:

• Operating systems (how many)
• Import data file types (how many)
• Import from databases (how many)
• Export data file types (how many)
• Languages displayable in UI (how many, besides English)
• Easy to repeat any step by groups (split-file)
• Multiple data files open at once
• Multiple output windows
• Multiple code windows
• Variable metadata view
• Variable types (how many)
• Variable search/filter in dialogs
• Variable sort by name
• Variable sort by type
• Variable move manually
• Model Builder (how many effect types)
• Magnify GUI for teaching
• R code editor
• Comment/uncomment blocks of code
• Package management (comes with R and all packages)
• Output: word processing features
• Output: R Markdown
• Output: LaTeX
• Data wrangling (how many)
• Transform across many variables at once (e.g., row mean)
• Transform down many variables at once (e.g., log, sqrt)
• Assign factor labels across many variables at once
• Project saves/loads data, dialogs, and notes in one file

Graphics

This category consists mainly of the number of plot types each software offers. However, the other items can be very important to completing your work. They should add more than one point to the graphics score, but I scored them one point since some will view them as very important while others might not need them at all. Be sure to see the full reviews or download the Excel file if those features are important to you. Figure 3 shows the total graphics score for each GUI. jamovi has a solid lead in this category. R-Instat’s ability to layer any geometric object, or “geom,” on top of another could substantially increase its count, but I don’t include that since it requires knowledge of R code.

When studying all the graphs in this article, it’s important to consider the difference between the relative and absolute performance. For example, relatively speaking, R Commander is not doing well here, but it does offer over 25 types of plots! That might be fine for your needs.

The graphics features are:

• Types of graphs (how many)
• Small multiples (i.e., facets)
• Large multiples (i.e., by group / split-file)
• Export graphics formats (how many)
• Interactive Plot Builder
• Interactive Plot Editor (e.g., add annotations, arrows, etc.)
• Able to combine any graphs into a single image

Continued…

The PROJ format for defining Coordinate Reference Systems (CRS) is becoming deprecated, as PROJ strings cannot accurately store modern datum transformations, coordinate epoch information, or complex axis order rules. However, that format is much easier and more practical for humans to type and read, and it is still often used as a legacy convenience. I, for one, use it often for centering a map projection on the mapped region, thus avoiding the distortion that occurs away from the projection center. For example, imagine I want to map New Zealand in a Lambert Azimuthal Equal-Area projection in R:

library(terra)library(geodata)nz <- gadm(country = "New Zealand", level = 0, path = tempdir())nz_prj <- project(nz, "+proj=laea")plot(nz_prj)

The map is heavily distorted, because it is far away from the center of the projection, which is 0,0 by default. With the PROJ format, we can easily set the center wherever we want, e.g.:

nz_centr <- crds(centroids(aggregate(nz)))nz_centr  # 172.6104 -41.82592nz_prj <- project(nz, "+proj=laea +lat_0=-41.8 +lon_0=172.6")# (this can be further automated to avoid the hard-coding,# but I keep it like this here for simplicity) plot(nz_prj)

Much better! We now have an equal-area map of New Zealand without visibly distorting the shape, and without having to look for specific equal-area CRSs normally used over there (and again every time we need to map another region of the world). This works as long as we use the WGS84 ellipsoid rather than a local datum, as those are currently not supported when using PROJ. So, I’m still using this regularly, and apprehensive about PROJ being deprecated… If anyone knows of a non-PROJ way of easily defining user-centered projections like this, please let me know!

After an extended beta period, JAGS 5.0.0 is now available. The release is more of a process than a single event, but we have started by uploading the source tarball and the Windows installer to SourceForge. The macOS installer will be along hopefully next week.

Thanks to Jeroen Ooms, Windows binary installers are available from the r-windows project on GitHub. This includes binaries for aarch64/Arm64.

Changes since the beta release

• Better multi-core support on macOS (More on this in a future post).
• Fixed a delayed rejection bug in the MNormal sampler.
• Added a new elliptical slice sampler contributed by Marcel Jonker.
• The dic.samples function could throw a logic error (See ticket #66). I am unable to reproduce this but JAGS 5.0.0 includes a workaround that should avoid triggering this bug.
• Fixed a reproducibility bug in BaseRNGFactory.

I would not normally add a new sampler at this late stage, but the elliptical slice sampler was ready to go and performs surprisingly well for a non-adaptive sampler. If you experience any trouble with this new sampler then you can turn it of with

set.factory("bugs::EllipticalSlice", type="sampler", state=FALSE)

JAGS will then fall back on the adaptive random walk Metropolis-Hastings sampler MNormal.

The rjags package and CRAN migration

The current rjags package on CRAN (version 4-17) will not work with JAGS 5.0.0. You need version 5-x but for the time being this will only be available from SourceForge.

• Source tarball
• Windows binary (x86_64)
• macOS binary

The reason for not uploading the rjags package to CRAN is that it currently breaks a lot of reverse dependencies. I have been quite proactive in contacting the maintainers of these package and providing patches, but we still need to give them more time.

We are monitoring the reverse dependency situation and will be liaising with the CRAN maintainers.

Preview of upcoming posts

In order to help users and developers with the transition to JAGS version 5, a series of “What’s new in JAGS 5.0.0” blog posts are planned covering the following topics.

• Initial values
• Reproducibility
• Multicore support
• ABI compatibility
• New diagnostics

These will include guidance for package maintainers to make the transition to JAGS 5.0.0.

I’ve put together a short guide: How Do I Build a Model? A starting guide for new modellers. It’s free, online, and aimed squarely at the graduate student who has a research question and is staring at a blank script wondering where to even begin.

You can read it online at How Do I Build a Model?, or download a pdf.

I’ve never settled into one corner of ecological modelling. I’m not the Bayesian guy, or the foodweb model guy, or the machine learning guy, though I’ve dabbled in all three. That haphazard career has given me a broad view across a lot of different modelling traditions, and I keep noticing the same meta-principles showing up regardless of the discipline. This guide is my attempt to write those down.

It’s not a textbook on any particular method. Once you know what type of model you’re building, go read a book by a discipline expert (I’ve listed some favourites in the final chapter). This guide is for the step before that: turning a question into a model in the first place.

In their paper “A Synchronized Multi‑IMU Wearable System for Tracking of Joint‑Angles in Sports Motion Analysis” (arXiv:2607.26027v1), Samarasekera and colleagues set out to solve a very practical problem: how to reliably measure joint angles in dynamic sports movements using wearable IMUs. Their goal was to design a synchronized pipeline that could filter out noise, correct drift, and normalize ranges so that biomechanical signals could be trusted in real‑world conditions.

Our goal is different but conceptually parallel: to apply the same pipeline logic to financial time‑series data. Instead of knee flexion angles, we analyze the comparative dynamics of Gold and Bitcoin. The challenge is similar — financial returns are noisy, prone to drift, and difficult to compare directly. By adapting their pipeline, we aim to extract latent trends, mitigate drift, normalize ranges, and highlight relative divergence between assets.

The pipeline consists of four sequential stages:

1. Indirect Kalman Filtering (IKF) – extracting latent trends from noisy log returns.
   • State equation:

$$x_t=F\cdot x_{t-1}+w_t,w_t\sim N(0,Q)$$

• Observation equation:

$$y_t=H\cdot x_t+v_t,v_t\sim N(0,R)$$

• IKF separates signal from noise, with $Q$ and $R$ learned via maximum likelihood.

2. Drift Mitigation via High‑Pass Filtering – removing low‑frequency bias.
   • High‑pass operator:

$$x_t^{HP}=x_t-\frac{1}{N}\sum _{i=1}^N{x}_{t-i}$$

• This ensures long‑term drift does not dominate short‑term dynamics.

3. Range Normalization – mapping filtered values into a comparable scale.
   • Normalization function:

$$z_t=\frac{x_t^{HP}-\min (x^{HP})}{\max (x^{HP})-\min (x^{HP})}\cdot 100$$

• This allows Gold and Bitcoin to be compared on the same 0–100 scale.

4. Relative Notation (Directional Divergence) – highlighting when assets move in opposite directions.
   • Divergence indicator:

$$D_t=\text{sign}(r_t^{Gold})\ne \text{sign}(r_t^{BTC})$$

• If true, shaded regions show divergence; if false, alignment.

Why This Matters

• For the biomechanics researchers, the pipeline meant trustworthy joint‑angle tracking.
• For us, the same architecture means trustworthy comparative trend analysis between financial assets.
• In both cases, the pipeline’s strength lies in its modularity: IKF for latent signal extraction, high‑pass filtering for drift control, normalization for comparability, and relative notation for divergence detection.

✦ R Packages in Our Pipeline

To make this pipeline reproducible and transparent, we rely on a set of R packages. Each package plays a specific role in the architecture:

1. tidyverse
   • Purpose: Data manipulation and wrangling.
   • Functions: mutate, rename, inner_join, drop_na.
   • Why: Provides a clean grammar for transforming raw financial data into structured time‑series.
2. tidyquant
   • Purpose: Financial data acquisition and transformation.
   • Functions: tq_get (download asset prices), tq_transmute (compute log returns).
   • Why: Bridges tidyverse with financial APIs, enabling reproducible market data pipelines.
3. KFAS
   • Purpose: State‑space modeling and Kalman filtering.
   • Functions: SSModel (define state‑space model), fitSSM (MLE parameter estimation), KFS (Kalman smoothing).
   • Why: Implements the mathematical backbone of IKF, separating latent trends from noise.
4. zoo
   • Purpose: Rolling window operations.
   • Functions: rollmean.
   • Why: Used for high‑pass filtering by subtracting rolling averages, mitigating drift.
5. scales
   • Purpose: Normalization and rescaling.
   • Functions: rescale.
   • Why: Maps filtered trends into a 0–100 range, enabling comparability across assets.
6. ggbraid
   • Purpose: Divergence visualization.
   • Functions: geom_braid.
   • Why: Highlights periods of divergence between Gold and Bitcoin with shaded regions.
7. ggplot2 (via tidyverse)
   • Purpose: Visualization.
   • Functions: geom_line, labs, theme_classic.
   • Why: Provides the plotting framework for presenting normalized high‑pass Kalman trends.

library(tidyverse)
library(tidyquant)
library(KFAS)
library(ggbraid)

# --- Step 1: Log returns ---
gold <- tq_get("GC=F", from = "2024-01-01") %>%
  tq_transmute(select = close,
               mutate_fun = periodReturn,
               period = "daily",
               type = "log")

btc <- tq_get("BTC-USD", from = "2024-01-01") %>%
  tq_transmute(select = close,
               mutate_fun = periodReturn,
               period = "daily",
               type = "log")

df <- gold %>%
  rename(ret_gold = daily.returns) %>%
  inner_join(btc %>% rename(ret_btc = daily.returns), by="date") %>%
  drop_na()

# --- Step 2: IKF (Kalman filter) ---
model_gold <- SSModel(df$ret_gold ~ SSMtrend(1, Q = list(NA)), H = NA)
model_btc  <- SSModel(df$ret_btc  ~ SSMtrend(1, Q = list(NA)), H = NA)

fit_gold <- fitSSM(model_gold, inits = c(log(var(df$ret_gold)), log(var(df$ret_gold))))
fit_btc  <- fitSSM(model_btc,  inits = c(log(var(df$ret_btc)),  log(var(df$ret_btc))))

kf_gold <- KFS(fit_gold$model, smoothing = c("state"))
kf_btc  <- KFS(fit_btc$model,  smoothing = c("state"))

trend_gold <- as.numeric(kf_gold$alphahat)
trend_btc  <- as.numeric(kf_btc$alphahat)

# --- Step 3: High-pass filter (remove low-frequency drift) ---
hp_gold <- trend_gold - zoo::rollmean(trend_gold, k=30, fill=NA, align="right")
hp_btc  <- trend_btc  - zoo::rollmean(trend_btc,  k=30, fill=NA, align="right")

# --- Step 4: Range normalization ---
gold_norm <- scales::rescale(hp_gold, to=c(0,100))
btc_norm  <- scales::rescale(hp_btc,  to=c(0,100))

# --- Step 5: Relative notation (directional divergence) ---
df_norm <- df %>%
  mutate(
    gold_norm = gold_norm,
    btc_norm  = btc_norm,
    sign_gold = sign(ret_gold),
    sign_btc  = sign(ret_btc),
    divergence = sign_gold != sign_btc
  )

# --- Plot ---
ggplot(df_norm, aes(x=date)) +
  geom_line(aes(y=gold_norm, color="Gold"), linewidth=1.2) +
  geom_line(aes(y=btc_norm, color="Bitcoin"), linewidth=1.2) +
  geom_braid(aes(ymin=gold_norm, ymax=btc_norm, fill=divergence), alpha=0.3) +
  scale_color_manual(values=c("Gold"="#FFD700", "Bitcoin"="#1E3A8A"),
                     name = "") +
  scale_fill_manual(values=c("TRUE"="#FF4C4C", "FALSE"="#87CEFA"),
                    name = "Divergence") +
  scale_y_continuous(limits=c(0,100)) +
  labs(title="Gold vs Bitcoin: Kalman-Filtered High-Pass Trends",
       subtitle="Directional divergence highlighted with shaded regions",
       x="", y="Normalized High-Pass Kalman Trend") +
  theme_classic(base_family = "Roboto Slab") +
  theme(
    legend.position="bottom",
    legend.title=element_text(size=12, face="bold"),
    legend.text=element_text(size=10),
    plot.title=element_text(size=16, face="bold"),
    plot.subtitle=element_text(size=12),
    axis.title=element_text(size=12, face="bold"),
    axis.text=element_text(size=10),
    panel.grid.minor=element_blank()
  )

Conclusion

By adapting the IKF pipeline from the sports biomechanics paper, we demonstrate how state‑space filtering, drift mitigation, and normalization can reveal comparative dynamics among financial assets. This modular approach is reproducible and extensible to other domains where latent trends and divergence matter.

This is a post based on two other posts I wrote for LinkedIn. It has become quite long.

Run cell.
Inspect output.
Modify.
Run again.
Oops.
Restart kernel.
Run all.

for ...
    ...

We have all seen colourful signposts with great cities and their distances – how far each is from this current standpoint.

And fundamental question is, how correct these distances are? Or even, one should ask, how honest a particular signpost is, regarding the current standpoint.

So let’s put this utterly most fundamental question to the test And since it is vacation time, you can grab the R code and check it yourself

Image stumbling upon the post and wanted to check if this post is really relevant for this location. It can be a photo from other side of the world, hanging in the lobby of your holiday hotel And this sure ain’t honest signpost

The R code iterates through each city, finds matching points (where all circles meet) and searches for the nearest location. It uses OpenStreetMap (short: OSM) to extract all the geolocations and looks for possible intersections of all points

geocode_city <- function(city_name) {
 
  url <- modify_url(
    "https://nominatim.openstreetmap.org/search",
    query = list(
      q       = city_name,
      format  = "json",
      limit   = 1
    )
  )
  
  resp <- tryCatch(
    GET(url, user_agent("PointingSignFinder/1.0 (R script)")),
     error = function(e) {
       cat(" Error on endpoint")
       return(NULL)
     }
  )
  
  if (is.null(resp) || http_error(resp)) {
    cat(" Error on http\n")
    return(NULL)
  }
  
  result <- fromJSON(content(resp, as = "text", encoding = "UTF-8"))
  
  if (length(result) == 0) {
    cat("Error on result")
    return(NULL)
  }
  
  lat <- as.numeric(result$lat[1])
  lon <- as.numeric(result$lon[1])
  cat(sprintf(" found: %.4f°, %.4f°\n", lat, lon))
  
  Sys.sleep(1.1)
  
  list(lat = lat, lon = lon, display_name = result$display_name[1])
}


sign_location_finder <- function(cities,
                                 distances,
                                 tolerance = 50,
                                 coarse_res = 0.2,
                                 fine_res   = 0.02,
                                 nearby_radius = 200,
                                 nearby_min_pop = 100000) {
  
  stopifnot(length(cities) == length(distances))
  stopifnot(length(cities) >= 2)
  
  cat("Step 1: Geocoding cities\n")
  coords <- lapply(cities, geocode_city)
  
  failed <- which(sapply(coords, is.null))
  if (length(failed) > 0) {
    stop(sprintf("Location not founc: %s", paste(cities[failed], collapse = ", ")))
  }
  
  sign_data <- data.frame(
    city    = cities,
    lat     = sapply(coords, `[[`, "lat"),
    lon     = sapply(coords, `[[`, "lon"),
    dist_km = distances
  )


print(sign_data[, c("city", "lat", "lon", "dist_km")])
  
max_dist_deg <- max(sign_data$dist_km) / 111
lat_min <- min(sign_data$lat) - max_dist_deg - 5
lat_max <- max(sign_data$lat) + max_dist_deg + 5
lon_min <- min(sign_data$lon) - max_dist_deg * 2 - 5
lon_max <- max(sign_data$lon) + max_dist_deg * 2 + 5
  
 
lat_min <- max(lat_min, -85)
lat_max <- min(lat_max,  85)
lon_min <- max(lon_min, -180)
lon_max <- min(lon_max,  180)
  
cat(sprintf("\nStep 2: Search bounding box: lat [%.1f, %.1f], lon [%.1f, %.1f]\n",
              lat_min, lat_max, lon_min, lon_max))
  
  
  score_point <- function(plat, plon) {
    diffs <- sapply(seq_len(nrow(sign_data)), function(i) {
      d <- distHaversine(
        c(plon, plat),
        c(sign_data$lon[i], sign_data$lat[i])
      ) / 1000
      abs(d - sign_data$dist_km[i])
    })
    max(diffs)
  }
  
  
  cat(sprintf("\nStep 3: Search the grid (%.2f° resolution)...\n", coarse_res))
  grid <- expand.grid(
    lat = seq(lat_min, lat_max, by = coarse_res),
    lon = seq(lon_min, lon_max, by = coarse_res)
  )
  cat(sprintf("Checking %d grid points...\n", nrow(grid)))
  grid$score <- mapply(score_point, grid$lat, grid$lon)
  
  candidates <- grid[grid$score <= tolerance, ]
  cat(sprintf(" Found %d candidate cells within ±%d km tolerance.\n",
              nrow(candidates), tolerance))
  
  coarse_best <- grid[which.min(grid$score), ]
  
  fine_grid <- expand.grid(
    lat = seq(coarse_best$lat - 1, coarse_best$lat + 1, by = fine_res),
    lon = seq(coarse_best$lon - 1, coarse_best$lon + 1, by = fine_res)
  )
  fine_grid$score <- mapply(score_point, fine_grid$lat, fine_grid$lon)
  best <- fine_grid[which.min(fine_grid$score), ]
  
  cat(sprintf("\n>>> Estimated  locations:\n"))
  cat(sprintf("    Latitude  : %.4f°\n", best$lat))
  cat(sprintf("    Longitude : %.4f°\n", best$lon))
  cat(sprintf("    Max error : ±%.1f km\n", best$score))
  cat(sprintf("    Google Maps: https://www.google.com/maps?q=%.4f,%.4f\n", best$lat, best$lon))
  
  cat(sprintf("\nStep 5: Nearby cities (within %d km, pop > %s)...\n",
              nearby_radius, format(nearby_min_pop, big.mark = ",")))
  
  world_cities <- world.cities
  nearby <- world_cities |>
    filter(pop > nearby_min_pop) |>
    mutate(
      dist_to_sign = distHaversine(
        cbind(long, lat),
        c(best$lon, best$lat)
      ) / 1000
    ) |>
    filter(dist_to_sign <= nearby_radius) |>
    arrange(dist_to_sign) |>
    select(name, country.etc, lat, long, pop, dist_to_sign) |>
    head(10)
  
  if (nrow(nearby) > 0) {
    cat("\n  Cities found:\n")
    print(nearby, digits = 4)
    cat(sprintf("\n  Nearest: %s, %s (%.1f km away)\n",
                nearby$name[1], nearby$country.etc[1], nearby$dist_to_sign[1]))
  } else {
    cat("No major locations found in vicinity!")
  }
  

}

So now, that we understand the useless problem, we can put this to the test

# sample 2 - real with correct distances calculater from#  https://www.distancefromto.net/  "air distance"# result must be Ljubljana, Sloveniaresult <- sign_location_finder(  cities    = c("Koper", "Celje", "Maribor", "Kranj"),  distances = c(83, 61, 104, 24),  tolerance = 20) 

In this case I have entered the correct air distances and there should be a location present – which is capital of Slovenia – Ljubljana.

I have inserted the real values of air distances between Ljubljana all all four cities: Koper, Celje, Maribor, Kranj using air distance calculator.

And this where the circles should meet

So if you are enjoying your vacation and you stumble upon a sign, check it and let me know, if it corresponds to your location of not

As always, the complete code is available on GitHub in  Useless_R_function repository and the file Honest_geographical_location_signposts.R is here.

Check the repository for future updates!

Stay healthy, hydrated and happy R-coding!

The package garchf wraps the extensive GARCH modeling machinery in rugarch behind a forecast-style interface, so that GARCH-based forecasts can be produced, plotted, and evaluated with the same conventions used throughout the forecast ecosystem (in the spirit of functions like forecast::thetaf()). The workhorse function, xgarchf(), fits a GARCH model to a univariate time series – with a choice of variance model (sGARCH, eGARCH, gjrGARCH, apARCH, iGARCH, …), ARMA order for the conditional mean, GARCH order, and conditional distribution – and returns a standard forecast object containing point forecasts, prediction intervals, and simulated paths. Because the output integrates with forecast’s S3 methods, plot() and other generics work out of the box, and the simulated paths (sims) make it straightforward to compute richer probabilistic forecasting metrics (CRPS, pinball loss, Winkler scores, etc.) via cross-validation. Below, we illustrate xgarchf() on Google’s daily closing stock price returns, and then run a small cross-validation exercise comparing a few GARCH variants using the crossvalidation package.

library(garchf)

Registered S3 method overwritten by 'quantmod':
  method            from
  as.zoo.data.frame zoo 

1 – Graphs

y <- diff(log(fpp2::goog200))

fit1 <- xgarchf(y, h = 20, model = "eGARCH")
fit2 <- xgarchf(y, h = 20, model = "sGARCH")
fit3 <- xgarchf(y, h = 20, model = "gjrGARCH")
fit4 <- xgarchf(y, h = 20, model = "iGARCH")

print(summary(fitted(fit1)))

    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
0.001504 0.001504 0.001504 0.001504 0.001504 0.001504     
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
0.002783 0.002783 0.002783 0.002783 0.002783 0.002783     
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
0.001233 0.001233 0.001233 0.001233 0.001233 0.001233     
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
0.001485 0.001485 0.001485 0.001485 0.001485 0.001485 

    Ljung-Box test

data:  Residuals
Q* = 10.341, df = 10, p-value = 0.4111

Model df: 0.   Total lags used: 10

    Ljung-Box test

data:  Residuals
Q* = 10.341, df = 10, p-value = 0.4111

Model df: 0.   Total lags used: 10

    Ljung-Box test

data:  Residuals
Q* = 10.341, df = 10, p-value = 0.4111

Model df: 0.   Total lags used: 10

    Ljung-Box test

data:  Residuals
Q* = 10.341, df = 10, p-value = 0.4111

Model df: 0.   Total lags used: 10

2 - Cross-validation with the package crossvalidation

For crossvalidation, see https://github.com/Techtonique/crossvalidation.

spl_m5 <- function(predicted,
                   observed,
                   probs = c(0.005, 0.025, 0.165, 0.25,
                             0.5,
                             0.75, 0.835, 0.975, 0.995))
{
  if (is.null(predicted$sims)) {
    stop("predicted$sims is required")
  }

  # training series from forecast object
  train <- as.numeric(predicted$x)

  # M5 scaling denominator
  scale_denom <- mean(abs(diff(train)), na.rm = TRUE)

  if (scale_denom <= 0 || !is.finite(scale_denom)) {
    scale_denom <- 1
  }

  sims <- predicted$sims

  # quantile forecasts: rows = horizons, cols = probs
  qhat <- sapply(
    probs,
    function(p) apply(sims, 1, stats::quantile, probs = p, na.rm = TRUE)
  )

  observed <- as.numeric(observed)

  # pinball loss
  pinball_vec <- numeric(length(probs))

  for (j in seq_along(probs)) {

    u <- probs[j]
    q <- qhat[, j]

    pinball_vec[j] <- mean(
      ifelse(
        observed >= q,
        u * (observed - q),
        (1 - u) * (q - observed)
      ),
      na.rm = TRUE
    )
  }

  # scaled pinball loss for each quantile
  spl_vec <- pinball_vec / scale_denom

  # average scaled pinball loss
  mean_spl <- mean(spl_vec)

  out <- c(
    mean_SPL = mean_spl
  )

  names(spl_vec) <- paste0("SPL_", probs)

  c(out, spl_vec)
}

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

  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)
    )
  )

  spl <- spl_m5(predicted, observed)

  c(
    RMSE = rmse,
    MAE = mae,
    Coverage95 = coverage,
    Winkler95 = mean(winkler),
    CRPS = mean(scoringRules::crps_sample(observed, predicted$sims)),
    mean_SPL = as.numeric(spl["mean_SPL"])
  )
}

eval_metric80 <- 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.2

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

  spl <- spl_m5(predicted, observed)

  c(
    ME = me,
    RMSE = rmse,
    MAE = mae,
    Coverage80 = coverage,
    Winkler80 = mean(winkler),
    CRPS = mean(scoringRules::crps_sample(observed, predicted$sims)),
    mean_SPL = as.numeric(spl["mean_SPL"])
  )
}

res <- crossvalidation::crossval_ts(
  y = y,
  initial_window = 150,
  horizon = 10,
  fixed_window = FALSE,
  fcast_func = garchf::xgarchf,
  eval_metric = eval_metric95,
  fit_params = list(arma_order=c(0, 0), 
                    garch_order = c(1, 1),
                    model = "eGARCH",
                    level=95),
  show_progress = FALSE
)
print(summary(res))

res <- crossvalidation::crossval_ts(
  y = y,
  initial_window = 150,
  horizon = 10,
  fixed_window = FALSE,
  fcast_func = garchf::xgarchf,
  eval_metric = eval_metric95,
  fit_params = list(arma_order=c(0, 0), 
                    garch_order = c(1, 1),
                    model = "sGARCH",
                    level=95),
  show_progress = FALSE
)
print(summary(res))

res <- crossvalidation::crossval_ts(
  y = y,
  initial_window = 150,
  horizon = 10,
  fixed_window = FALSE,
  fcast_func = garchf::xgarchf,
  eval_metric = eval_metric95,
  fit_params = list(arma_order=c(0, 0), 
                    garch_order = c(1, 1),
                    model = "iGARCH",
                    level=95),
  show_progress = FALSE
)
print(summary(res))

      RMSE               MAE             Coverage95       Winkler95      
 Min.   :0.006459   Min.   :0.004896   Min.   :0.7000   Min.   :0.03441  
 1st Qu.:0.008001   1st Qu.:0.006266   1st Qu.:0.9000   1st Qu.:0.06153  
 Median :0.009135   Median :0.007061   Median :1.0000   Median :0.06466  
 Mean   :0.017189   Mean   :0.010613   Mean   :0.9475   Mean   :0.20090  
 3rd Qu.:0.022415   3rd Qu.:0.015066   3rd Qu.:1.0000   3rd Qu.:0.47706  
 Max.   :0.042741   Max.   :0.023054   Max.   :1.0000   Max.   :0.96625  
      CRPS             mean_SPL     
 Min.   :0.004154   Min.   :0.1170  
 1st Qu.:0.004922   1st Qu.:0.1362  
 Median :0.005311   Median :0.1481  
 Mean   :0.010039   Mean   :0.3380  
 3rd Qu.:0.017800   3rd Qu.:0.6763  
 Max.   :0.051416   Max.   :1.6607        

      RMSE               MAE             Coverage95       Winkler95      
 Min.   :0.006554   Min.   :0.005380   Min.   :0.7000   Min.   :0.03856  
 1st Qu.:0.008424   1st Qu.:0.006561   1st Qu.:0.9750   1st Qu.:0.06261  
 Median :0.009314   Median :0.007401   Median :1.0000   Median :0.07277  
 Mean   :0.017155   Mean   :0.010741   Mean   :0.9525   Mean   :0.17716  
 3rd Qu.:0.020081   3rd Qu.:0.013012   3rd Qu.:1.0000   3rd Qu.:0.34988  
 Max.   :0.042568   Max.   :0.022999   Max.   :1.0000   Max.   :0.49756  
      CRPS             mean_SPL     
 Min.   :0.004634   Min.   :0.1349  
 1st Qu.:0.005383   1st Qu.:0.1516  
 Median :0.005695   Median :0.1606  
 Mean   :0.008842   Mean   :0.2986  
 3rd Qu.:0.012903   3rd Qu.:0.4725  
 Max.   :0.019442   Max.   :0.7436        

    RMSE               MAE             Coverage95       Winkler95      
 Min.   :0.006556   Min.   :0.005514   Min.   :0.7000   Min.   :0.03580  
 1st Qu.:0.008306   1st Qu.:0.006239   1st Qu.:0.9750   1st Qu.:0.06307  
 Median :0.009119   Median :0.007007   Median :1.0000   Median :0.06478  
 Mean   :0.017029   Mean   :0.010484   Mean   :0.9525   Mean   :0.16725  
 3rd Qu.:0.018881   3rd Qu.:0.011711   3rd Qu.:1.0000   3rd Qu.:0.16760  
 Max.   :0.042713   Max.   :0.023064   Max.   :1.0000   Max.   :0.51546  
      CRPS             mean_SPL     
 Min.   :0.004851   Min.   :0.1363  
 1st Qu.:0.005317   1st Qu.:0.1510  
 Median :0.005700   Median :0.1591  
 Mean   :0.008720   Mean   :0.2928  
 3rd Qu.:0.009316   3rd Qu.:0.3126  
 Max.   :0.019649   Max.   :0.7586  

LLMs can produce a convincing media recommendation in a few seconds. The more useful question is whether the recommendation is correct: are the reach calculations right, are the assumptions visible, and does the plan fit the brief?

I work at Havas Media and I am the author of the Havas AI Media Quality Index (HAI-Q), a benchmark for media planning in the German market. I built it because generic knowledge tests are not enough for professional AI applications. A model can write a plausible paragraph about media strategy and still get a budget allocation, a target group definition, or a GRP calculation wrong.

This post has three parts. First, I will build a small media-planning evaluation in R using the vitals package. Then I will discuss why media-planning questions are unusually difficult to evaluate. Finally, I will focus on what the HAI-Q results tell us about general-purpose LLMs.

1. A small eval in R

At its simplest, an LLM eval has three components:

1. A dataset containing questions and targets.
2. A solver that answers the questions.
3. A scorer that evaluates the answers.

This is the basic structure used by vitals. The target column does not always need to contain one exact answer. For open-ended questions, it can contain grading guidance: what a good answer must include, which assumptions must be stated, and which errors should lose credit.

This is useful for testing the final answer of an LLM. It is less useful as a complete test of an AI agent. An agent consists of a harness and one or more models: it may retrieve data, call tools, run calculations, retry a failed step, maintain state, and decide when to ask for help. vitals mainly gives us a dataset, a solver, and a scorer. It does not automatically tell us whether the agent used the right tool, passed the right arguments, or recovered safely from an intermediate failure. It can be extended with custom solvers, but then the agent harness and trajectory checks still have to be built separately.

The questions

Here is a deliberately small evaluation set. The questions cover calculation, planning, and interpretation. The targets are written as instructions for a grader rather than as model answers.

library(dplyr)
library(tibble)

media_questions <- tribble(
  ~id, ~input, ~target, ~domain, ~task,

  "grp-calculation",
  "A campaign delivers 60% reach at an average frequency of 3. What is the campaign's GRP? State the assumption behind the calculation.",
  "A full-credit answer calculates 60 * 3 = 180 GRP. It explains that this is the standard simplified relationship GRP = reach (in percent) * average frequency. It must not describe GRP as the number of unique people reached.",
  "Reach and GRP", "Calculation",

  "budget-allocation",
  "A client has EUR 100,000 to reach adults aged 25-49 in Germany. There is no historical channel performance data. Propose a first-pass media allocation and explain how you would improve it after launch.",
  "There is no single correct allocation. A strong answer makes the uncertainty explicit, explains the role of each channel, avoids inventing CPMs or reach figures, and proposes a measurement and test-and-learn plan. It should distinguish a planning hypothesis from an observed result.",
  "Media strategy", "Recommendation",

  "reach-frequency-tradeoff",
  "Two plans have the same budget. Plan A has higher reach and lower frequency; Plan B has lower reach and higher frequency. How would you decide between them?",
  "A full-credit answer says that the choice depends on the communication objective, the selected channesls, the overall attention, audience size, purchase cycle, creative strength, and expected response. It should explain the reach-frequency trade-off and ask for missing information instead of declaring one plan universally better.",
  "Media strategy", "Reasoning"
)

media_questions

The questions are intentionally specific. “Create a good media plan” is difficult to score because almost any fluent answer can sound reasonable. A question about GRP, in contrast, gives us a concrete numerical check. The budget question tests something different: whether the model knows when it does not have enough information to justify a precise allocation.

In a real eval, I would expand this table with questions about channel roles, target-group definitions, reach curves, media costs, campaign objectives, and German market data. I would also add metadata such as difficulty, source, market, and the date on which the underlying data was valid.

The solver and scorer

The vitals introduction uses Task$new() to combine the dataset, a solver, and a scorer. The following example uses Claude as the solver and an LLM-based grader. An API key for the selected provider is required.

library(vitals)
library(ellmer)

media_task <- Task$new(
  dataset = media_questions,
  solver = generate(
    chat_openrouter(model = "claude-sonnet-4.6")
  ),
  scorer = model_graded_qa(partial_credit = TRUE),
  name = "A small media-planning eval"
)

# Run the solver and scorer.
media_task$eval()

# Bind the results into a tibble for analysis.
media_scores <- vitals_bind(media_task)
media_scores

media_scores |>
  count(score)

With partial_credit = TRUE, the scorer can return C for correct, P for partially correct, or I for incorrect. This is more useful than a single binary score for media planning. A response may use the right calculation but fail to explain its assumptions. Another answer may be directionally useful but invent a media cost. Those are different failures and should not be hidden in one number.

The result is not the end of the evaluation. I would inspect the individual answers and the grading decisions, especially on the first few runs. If the grader repeatedly gives full credit to a response that I consider unsafe, the problem may be the grading guidance rather than the model.

Comparing models

Once the task exists, comparing solvers is straightforward. For example, the same questions can be evaluated with an OpenAI model and then combined with the first run.

media_task_openai <- media_task$clone()

media_task_openai$eval(
  solver_chat = chat_openrouter(model = "gpt-5.4")
)

comparison <- vitals_bind(
  claude = media_task,
  openai = media_task_openai
) |>
  mutate(model = recode(
    task,
    claude = "Claude",
    openai = "OpenAI"
  ))

comparison |>
  count(model, score)

The exact model names will change over time. The important part is to keep the questions and scoring rules stable while changing one component at a time. Otherwise, a score difference is difficult to interpret.

This example is an LLM eval, not an agent eval. It compares the text returned at the end of the solver. For an agent, I would additionally log the tool calls, retrieved documents, intermediate calculations, retries, latency, cost, and final answer. The test should then assert properties of the whole trajectory, not only the final prose.

2. Why media-planning evals are hard

Media planning looks like a natural use case for LLMs. The work involves briefs, audience descriptions, channel recommendations, calculations, and written rationales. But the same combination makes the evaluation difficult.

There is rarely one correct plan

Many planning questions are underdetermined. A sensible allocation depends on the objective, budget, campaign length, audience size, creative assets, buying conditions, historical performance, and measurement setup. Two experienced planners can recommend different channel mixes and both be right.

This means that a target should often describe required properties rather than prescribe one number. For example, a full-credit answer might need to:

• state its assumptions;
• stay within the budget;
• connect channel choices to the objective;
• distinguish facts from hypotheses; and
• explain how the plan would be measured and revised.

The evaluator has to reward sound reasoning without turning one planner’s preference into a false ground truth.

Media facts are local and time-dependent

A media answer can be correct in one market and wrong in another. Reach, costs, inventory, audience definitions, platform availability, and panel methodology all vary by country and by date. A model that has learned general media vocabulary does not automatically have access to the latest German market data.

This is why a useful eval should record the market and the data vintage. It should also test whether the model provides a source or clearly labels an assumption when a current figure is unavailable. “The average CPM is X” is not a safe answer if the model cannot explain where X came from.

Definitions matter more than fluent prose

Media planning has many terms that are close enough to sound interchangeable but are not interchangeable in a calculation. Reach can mean gross or net reach. Frequency can be average or effective frequency. GRP is not the same as the number of unique people reached. Target-group percentages need a denominator.

These distinctions create easy failure modes. A response can contain the right formula but apply it to the wrong quantity. It can also produce a tidy table whose totals do not add up. Numerical checks and constraint checks therefore need to be part of the eval, not just a general assessment of writing quality.

LLM-as-a-judge is useful but not sufficient

The model_graded_qa() scorer in vitals makes open-ended evaluation scalable. It is also a model judging another model’s answer, which introduces its own errors. A grader may prefer confident language, overlook a subtle numerical mistake, or accept an answer because it resembles the supplied guidance.

For a serious benchmark, I would calibrate the grader against human-labelled examples and keep a sample for manual review. I would also separate different dimensions where possible:

• factual accuracy;
• numerical accuracy;
• constraint satisfaction;
• quality of reasoning; and
• usefulness of the recommendation.

One overall score is convenient, but it can hide which part of the planning process is failing.

vitals is not an agent harness

The distinction matters in practice. vitals is a useful way to ask an LLM a fixed set of questions and score the answers. It is not, by itself, a test harness for an AI agent. If an agent has access to a reach database, a calculator, or a planning API, the evaluation needs to check the interaction with those components as well.

For an end-to-end agent eval, I would test at least four layers: whether the harness selects the right tool, whether the tool call is valid, whether the intermediate result is interpreted correctly, and whether the final answer communicates uncertainty. A final response can look reasonable even when one of those earlier steps failed.

3. Results from the Havas HAI-Q

The Havas HAI-Q page describes a benchmark built for media planning in the German market. It contains 35 practical tasks across five areas:

1. Media strategy
2. Channel and budget allocation
3. Target-group evaluation
4. Reach and GRP
5. Media and effectiveness metrics

In this post I focus on the general-purpose LLM results:

The headline is uncomfortable: even GPT-5, the strongest general-purpose model in this comparison, is correct on fewer than half of the questions. It answers 16 of 35 questions correctly, leaving 19 answers that are not counted as correct. For practical media-planning work, that is not a safe error rate.

When I look at the individual HAI-Q questions, the errors are not evenly distributed. The models particularly struggle with numeric questions involving reach figures and audience sizes. They can produce a fluent explanation of a planning idea and still make a basic mistake in a denominator, a percentage, or the interpretation of a reach number. This is exactly the type of failure that is easy to miss if evaluation is based on plausibility rather than calculation.

There is also a clear improvement between the older and newer models. GPT-4o answers 4 of 35 questions correctly, compared with 16 of 35 for GPT-5. That is 12 additional correct answers and an improvement from 11.4% to 45.7%, or 34.3 percentage points. GPT-5 gets four times as many questions correct as GPT-4o on this benchmark. The improvement is substantial, even though the newer model is still not reliable enough for unsupervised media decisions.

Claude Sonnet 4.5 scores lower in this test, with 6 of 35 correct answers. This does not mean that the model is generally poor. It means that a general language benchmark and a domain-specific media-planning benchmark answer different questions. The HAI-Q is designed around the German market and a specific definition of correctness, including the numerical details that matter in planning.

The aggregate scores should still be read with care. They do not show the difficulty of every question, the variance across repeated runs, or the exact scoring decision for every answer. But the pattern is clear enough to be useful: current general-purpose LLMs are good at sounding like media planners, while still struggling to answer roughly half of these practical questions correctly.

Takeaways

An LLM eval for media planning should start with real questions from the workflow, not generic prompts about marketing. In R, a tibble with input and target is enough to get started, and vitals provides a practical way to connect questions, solvers, scorers, and result inspection.

The difficult part is not running the model. It is defining what “correct” means when several plans can be reasonable, data changes by market and date, and a plausible answer can contain a subtle but consequential error. For agents, the evaluation also needs to cover the harness, tools, intermediate steps, and failure handling. A final answer score is not enough.

The HAI-Q results show that the latest general-purpose LLMs still answer more than half of these media-planning questions incorrectly. The largest weaknesses appear in numerical reasoning around reach and audience sizes. At the same time, the jump from GPT-4o to GPT-5 shows that model development is moving in the right direction.

For me, the practical conclusion is simple: use LLMs to accelerate media work, but test every numerical recommendation before relying on it. The next step is to combine model evals such as HAI-Q with end-to-end agent tests that measure tool use and the full planning workflow.

Resources

• Getting started with vitals
• vitals documentation
• HAI-Q: Havas AI Media Quality Index
• ellmer documentation

Like many people I have been keeping track of the slow-burn fuel crisis arising from the war in Iran, its possible expansion to the Red Sea, and exacerbation from the ongoing war between Russia and Ukraine. The latter seems to have reached a stage where serious damage is being done to Russia’s refining capacity and petro products’ exports.

I’ve now started a page on this website with a few charts on the fuel crisis that I’m going to be keeping up to date.

Stocks of crude are going down, but “how much is left”?

The motivation for today’s post was that I’d been seeing many news and opinion articles about how the US crude oil stocks—in total, and in the Strategic Petroleum Reserve in particular—have declined to modern record lows. But expert commentators are reluctant to show this decline in relation to the absolute amount left, as opposed to comparisons to seasonal or historic lows. This is understandable because no-one really knows (or perhaps tells) what the real “tank bottoms” are in a big complex system like the USA. A certain amount of oil, of uncertain volume but certainly in the tens (and possibly hundreds) of millions of barrels, is needed just to keep various bits of plumbing or storage from seizing up or collapsing.

So it’s hard to get a chart like this one, and I had to make it myself:

I’ve got no reputation as an oil expert to worry about, so I can safely draw those lines at a fairly arbitrary 24 (or 30) days of refinery throughput, and calculate the time it will take to get to those low levels “at the current rate”. It’s simplistic, but it’s an important reality check to see that there might be about 19 weeks left.

Fast as the stocks are being used up, the USA has a lot of oil available to it. But 19 weeks isn’t that many. Markets are betting that good things will happen in the weeks or months before the blue line reaches the red line. Because of that, oil is currently selling for around $90 per barrel instead of the $110+ that would be justified if people thought we had lost the traffic through the Strait of Hormuz indefinitely.

Not that I think 19 weeks is a realistic forecast. Long before US stocks reached that red line, there would be economic and political chaos. Efforts to keep retail prices artificially low would probably be overwhelmed and there would be material demand destruction. There would also be huge incentives for the US to give away anything to restore Persian Gulf and Red Sea transit routes. 19 weeks gives us a time frame within which something will have to change.

The value-add for me from this chart is showing the total oil available for the US in the event this becomes a real emergency for them; the comparison of that to a meaningful if arbitrary threshold in terms of days of refinery cover; and the calculation of weeks left at current draw-down rates to reach that high stress threshold.

Break out the SPR and also look at gasoline, distillates

The same data source has a lot of other detailed information, most of it far too detailed for me, but I’m interested enough in four key headline variables that are more detailed than “total stocks of crude oil”. In this chart I break down the total stocks into those in the Strategic Petroleum Reserve, and those in more normally accessed stocks. I also show the amounts held of two key refined products:

Code

Here’s the code that grabs the data from the US Energy Information Administration and draws that second plot:

library(tidyverse)
library(readxl)
library(scales)
library(glue)

fn <- "us-petrol-status-weekly.xls"

download.file(
  "https://ir.eia.gov/wpsr/psw01.xls",
  fn,
  mode = "wb"
)

us_stocks <- read_excel(fn, sheet = "Data 1", skip = 2) |>
  rename(
    crude = `Weekly U.S. Ending Stocks of Crude Oil  (Thousand Barrels)`,
    crude_spr = `Weekly U.S. Ending Stocks of Crude Oil in SPR  (Thousand Barrels)`,
    gasoline = `Weekly U.S. Ending Stocks of Total Gasoline  (Thousand Barrels)`,
    diesel = `Weekly U.S. Ending Stocks of Distillate Fuel Oil  (Thousand Barrels)`,
    date = Date
  )

#------------------plotting------------------

eia_caption <- glue("Source: Energy Information Administration (EIA). Accessed {format(Sys.Date(), '%d %B %Y')}")

# an annotation rectangle we can use in multiple charts
war_rect <- annotate(
  "rect",
  xmin = as.Date("2026-02-28"),
  xmax = Inf,
  ymin = -Inf,
  ymax = Inf,
  alpha = 0.5,
  fill = "grey80"
)

# Facet plot of total crude, SPR crude, diesel and gasoline stocks
us_stocks |>
  mutate(com_crude = crude - crude_spr) |>
  select(
    date,
    `Non-SPR crude oil` = com_crude,
    `Crude oil in SPR` = crude_spr,
    `Gasoline` = gasoline,
    `Distillate (mostly diesel)` = diesel
  ) |>
  gather(variable, value, -date) |>
  mutate(variable = fct_reorder(variable, value)) |>
  ggplot(aes(x = date, y = value / 1000)) +
  war_rect +
  facet_wrap(~variable, scales = "free_y") +
  geom_line(colour = "blue") +
  expand_limits(y = 0) +
  scale_y_continuous(label = comma) +
  labs(
    x = "",
    y = "Thousands of barrels",
    title = "US stocks of crude oil, gasoline, and distillate fuel oil (effectively diesel)",
    subtitle = "Showing both total crude oil stocks (crude) and those in the Strategic Petroleum Reserve (crude_spr)",
    caption = eia_caption
  )

As I side point, I’ve started using air to format my code. Nearly always it comes out looking better and easier to read than my manual formatting, and it’s certainly more consistent. Thanks to D. Vaughan and L. Henry and the tidyverse.org project.

And here’s the code to draw my highly polished main plot, the one of total crude oil stocks with the annotations. I’ve tried to do this in a way that it will keep working in the months ahead with minimal maintenance, we’ll see how that goes.

refinery_throughput <- 17.3 # as at 24 July, operating at 97% of US capacity. Daily refinery use.
plausible_stress <- 30 * refinery_throughput
plausible_high_stress <- 24 * refinery_throughput

crude_growth_summary <- us_stocks |>
  summarise(
    latest_crude = crude[date == max(date)],
    weeks = as.numeric(as.Date(max(date)) - as.Date("2026-04-03")) / 7,
    difference = (latest_crude -
      crude[date == as.Date("2026-04-03")]) /
      1000,
    ratio = latest_crude / crude[date == as.Date("2026-04-03")],
    # growth/decline per week in million barrels e.g. 10m barrels per week:
    difference_rate = difference / weeks,
    # growth/declien rate per week:
    growth_rate = 1 - exp(log(ratio) / weeks)
  ) |>
  mutate(
    weeks_at_this_rate = (plausible_stress - latest_crude / 1000) /
      difference_rate
  )

plot_us_stocks <- function(min_date = "2020-01-01", lab_x_diff = NULL) {
  if (!"Date" %in% class(min_date)) {
    min_date <- as.Date(min_date)
  }

  if (is.null(lab_x_diff)) {
    days_shown <- as.numeric(as.Date(max(us_stocks$date))) -
      as.numeric(min_date)
    lab_x_diff <- days_shown / 50
  }

  lv <- tail(us_stocks, 1)$crude / 1000

  p2 <- us_stocks |>
    filter(date >= min_date) |>
    ggplot(aes(x = date, y = crude / 1000)) +
    war_rect +
    geom_hline(yintercept = plausible_stress, colour = "darkred") +
    geom_hline(yintercept = plausible_high_stress, colour = "red") +
    annotate(
      "text",
      x = min_date + lab_x_diff,
      y = plausible_stress + 50,
      label = "Illustrative stress threshold:\n30 days of refinery cover",
      colour = "darkred",
      hjust = 0
    ) +
    annotate(
      "text",
      x = min_date + lab_x_diff,
      y = plausible_high_stress - 30,
      label = "Illustrative high stress threshold:\n24 days of refinery cover",
      colour = "red",
      hjust = 0,
      vjust = 1
    ) +
    annotate(
      "text",
      x = max(us_stocks$date) - 5e6,
      y = lv,
      label = glue(
        "{round(lv)} million barrels;\n{round(lv / refinery_throughput)} days of cover"
      ),
      size = 2.9,
      hjust = 1,
      vjust = 1,
      colour = "blue"
    ) +
    geom_line(colour = "blue") +
    expand_limits(y = 0) +
    scale_y_continuous(
      label = comma,
      sec.axis = sec_axis(
        ~ . / refinery_throughput,
        name = "Days of refinery throughput"
      )
    ) +
    labs(
      x = "",
      y = "Millions of barrels",
      title = "U.S. Total Crude Oil Stocks, including Strategic Petroleum Reserve",
      subtitle = glue(
        "Comparison of existing inventories with refinery throughput as at July 2026.
Decline since peak on 3 April 2026 is at {abs(round(crude_growth_summary$difference_rate, 1))} million barrels per week; {round(crude_growth_summary$weeks_at_this_rate)} weeks from stress threshold at this (hypothetical and linear) rate."
      ),
      caption = eia_caption
    ) +
    theme(axis.line.y = element_line(colour = "grey50"))

    svglite(
      glue("all-us-crude-from{min_date}.svg"),
      width = 10,
      height = 7
    )
    print(p2)
    dev.off()
}

plot_us_stocks("2000-01-01")
plot_us_stocks("2020-01-01")
plot_us_stocks("2025-01-01")

Because I did this plot in the form of a function, I have a couple of variants of the plot. One starting in 2020:

And the other zoomed right in to just January 2025 and onwards:

On the whole I prefer the plot I started with, showing the situation from 2000 and onwards. I’m trying to get some big picture perspective after all.

That’s all for today. The images on this page will stay the same as at the time of writing; but those on the fuel crisis monitoring page will be updated each week.

era v0.6.0 is now available on CRAN:

install.packages("era")
library(era)

This minor release adds functions for chronological comparison of year vectors:

x <- yr(c(200, 100, 300), "BCE")
yr_earlier_than(x, yr(150, "BCE"))
#> # BCE years <yr_lgl[3]>:
#> [1] FALSE FALSE  TRUE
#> # Era: Before Common Era (BCE): Gregorian years (365.2425 days), counted backwards from 1

yr_later_than(x, yr(150, "BCE"))
#> # BCE years <yr_lgl[3]>:
#> [1]  TRUE FALSE FALSE
#> # Era: Before Common Era (BCE): Gregorian years (365.2425 days), counted backwards from 1

And for calculating chronological differences between year vectors:

yr_difference(yr(300, "BCE"), yr(100, "BCE"))
#> # years <yr[1]>:
#> [1] 200
#> # Era: Gregorian years (365.2425 days)

These functions are aware of era directionality: “BCE” years count backwards, so 300 BCE is earlier than 100 BCE. They support work on tempo, a new package I am working on that represents chronological intervals and the temporal relations between them.

Links

• era on CRAN
• era package documentation
• Introductory vignette
• Source code (GitHub)

The goal of spuriouscorrelations is to keep alive the amazing examples from Tyler Vigen. Unfortunately, as of 2023-10-09, the website is down as my students noticed. Therefore, I decided to use the snapshot from the Internet Wayback Machine to save the datasets from 2023-06-07.

Installation

You can install the CRAN version of spuriouscorrelations with:

install.packages("spuriouscorrelations")

You can install the development version of spuriouscorrelations with:

remotes::install_github("pachadotdev/spuriouscorrelations")

Example

The package covers multiple examples for different spurious (and curious) correlations.

Here is a basic example which shows you how to plot a spurious correlation for the variables

• x: Number of people who drowned by falling into a pool
• y: Films Nicolas Cage appeared in

library(spuriouscorrelations)
library(tinyplot)

pool_drownings
   year   x y
1  1999 109 2
2  2000 102 2
3  2001 102 2
4  2002  98 3
5  2003  85 1
6  2004  95 1
7  2005  96 2
8  2006  98 3
9  2007 123 4
10 2008  94 1
11 2009 102 4
cor(pool_drownings$x, pool_drownings$y)
[1] 0.6660043
tinyplot(
  y ~ x,
  data = pool_drownings,
  main = sprintf("Correlation %s", round(cor(pool_drownings$x, pool_drownings$y), 3)),
  xlab = "Number of people who drowned by falling into a pool",
  ylab = "Films Nicolas Cage appeared in"
)

Converting the data to long format simplifies plotting both variables per year:

pool_drownings_2 <- reshape(
  pool_drownings, 
  varying = c("x", "y"),  # columns to collapse
  v.names = "value",      # name of the new value column
  timevar = "variable",   # name of the new ID column
  times = c("x", "y"),    # values to populate the ID column
  direction = "long"      # target format
)

pool_drownings_2
     year variable value id
1.x  1999        x   109  1
2.x  2000        x   102  2
3.x  2001        x   102  3
4.x  2002        x    98  4
5.x  2003        x    85  5
6.x  2004        x    95  6
7.x  2005        x    96  7
8.x  2006        x    98  8
9.x  2007        x   123  9
10.x 2008        x    94 10
11.x 2009        x   102 11
1.y  1999        y     2  1
2.y  2000        y     2  2
3.y  2001        y     2  3
4.y  2002        y     3  4
5.y  2003        y     1  5
6.y  2004        y     1  6
7.y  2005        y     2  7
8.y  2006        y     3  8
9.y  2007        y     4  9
10.y 2008        y     1 10
11.y 2009        y     4 11
tinyplot(
  value ~ year | variable, # "|" indicates the groping variable for the legend
  data = pool_drownings_2,
  main = sprintf("Correlation %s", round(cor(pool_drownings$x, pool_drownings$y), 3)),
  xlab = "Year",
  ylab = "Pooled observations",
  pch = 19 # solid dot shape
)

How about standarzing the variables to avoid the different scale visibility issue?

pool_drownings_3 <- pool_drownings

pool_drownings_3$x <- (pool_drownings_3$x - mean(pool_drownings_3$x)) / sd(pool_drownings_3$x)
pool_drownings_3$y <- (pool_drownings_3$y - mean(pool_drownings_3$y)) / sd(pool_drownings_3$y)

pool_drownings_3 <- reshape(
  pool_drownings_3, 
  varying = c("x", "y"),  # columns to collapse
  v.names = "value",      # name of the new value column
  timevar = "variable",   # name of the new ID column
  times = c("x", "y"),    # values to populate the ID column
  direction = "long"      # target format
)

pool_drownings_3
     year variable      value id
1.x  1999        x  0.8952867  1
2.x  2000        x  0.1696333  2
3.x  2001        x  0.1696333  3
4.x  2002        x -0.2450258  4
5.x  2003        x -1.5926679  5
6.x  2004        x -0.5560201  6
7.x  2005        x -0.4523554  7
8.x  2006        x -0.2450258  8
9.x  2007        x  2.3465935  9
10.x 2008        x -0.6596849 10
11.x 2009        x  0.1696333 11
1.y  1999        y -0.2470999  1
2.y  2000        y -0.2470999  2
3.y  2001        y -0.2470999  3
4.y  2002        y  0.6589330  4
5.y  2003        y -1.1531327  5
6.y  2004        y -1.1531327  6
7.y  2005        y -0.2470999  7
8.y  2006        y  0.6589330  8
9.y  2007        y  1.5649658  9
10.y 2008        y -1.1531327 10
11.y 2009        y  1.5649658 11
tinyplot(
  value ~ year | variable, # "|" indicates the groping variable for the legend
  data = pool_drownings_3,
  main = sprintf("Correlation %s", round(cor(pool_drownings$x, pool_drownings$y), 3)),
  xlab = "Year",
  ylab = "Pooled standardized observations",
  pch = 19, # solid dot shape
  type = "b" # add line to see the trend clearly
)

Join our workshop on Introduction to All Things About bpvars, the R package for Forecasting with Bayesian Panel Vector Autoregressions,  which is a part of our workshops for Ukraine series! 

Here’s some more info: 

Title: All Things About bpvars, the R package for Forecasting with Bayesian Panel Vector Autoregressions

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

Speaker: Tomasz Woźniak is a Bayesian econometrician developing new econometric methods for applied macroeconomic research. He has been a specialised R user for eighteen years, and has recently joined The R Journal as an Associate Editor. He is the author of several R packages, available at https://bsvars.org/, that combine blazingly fast algorithms written in C++ with the convenience of data analysis in R. He works as a senior lecturer at the University of Melbourne, where he has an extensive research, teaching, and engagement portfolio. Tomasz supports Ukraine. Slava Ukrainie!

Description: This session is a unique opportunity to learn to use the bpvars package for forecasting with Bayesian panel vector autoregressions in a wide range of applications spanning:

* data preparation,

* specification of panel vector autoregressions,

* estimation of the models,

* estimation with missing observations,

* forecasting labour market outcomes globally,

* plotting and reporting predictions of different qualities, including point and density forecasts, marginal, conditional, or restricted forecasts,

* performing recursive expanding window forecasting,

* generating forecasting performance reports.

A sequence of hands-on exercises that fully prepares attendees to work with the package supports all of this. You are welcome to browse existing resources, including:

* bpvars website: https://bsvars.org/bpvars/

* bpvars CRAN profile: https://cran.r-project.org/package=bpvars

* package vignette: https://doi.org/10.48550/arXiv.2606.14143

* other materials: https://bsvars.org/bpvars/#resources

See you soon!

Preparation:

This session assumes the attendees know the basics of time series analysis. Install and load the bpvars package, then check the documentation by running the code below. Run the example you will find there in R. If it works, you’re ready to go. Thanks!

“`

install.packages(“bpvars”)

library(bpvars)

?bpvars

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

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

How can I register?

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

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

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

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

How can I sponsor a student?

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

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

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

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

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

Looking forward to seeing you during the workshop!

I’m a big fan of Github Copilots ghost text and next edit suggestions. As a data analyst they give you much more control over the pace and direction of analysis code than a full AI agent. They also let you practice the core skill an agent needs from you: writing a clear specification.

I’ll walk through it in R, but the same ideas apply in any language. I’m assuming you’re on VScode with the Copilot extension installed. This should work with the free or paid plans (though I haven’t checked free myself).

library(dplyr)
library(ggplot2)
dat <- data.frame(
    time = rep(1:10, 3),
    species = rep(c("A", "B", "C"), each = 10),
    value = c(rnorm(10, mean = 5), rnorm(10, mean = 10), rnorm(10, mean = 15))
)

datA <- dat |> filter(species == "A")
m1 <- lm(value ~ time, data = datA)
coef(m1)["time"]

ggplot(datA, aes(x = time, y = value)) +
    geom_point() +
    geom_smooth(method = "lm", se = FALSE) +
    labs(
        title = "Species A Time Series",
        x = "Time",
        y = "Value"
    )

fit_fun <- function(species_name) {
    datA <- dat |> filter(species == species_name)
    m1 <- lm(value ~ time, data = datA)
    coef(m1)["time"]

    ggplot(datA, aes(x = time, y = value)) +
        geom_point() +
        geom_smooth(method = "lm", se = FALSE) +
        labs(
            title = paste("Species", species_name, "Time Series"),
            x = "Time",
            y = "Value"
        )
}

fit_fun("A")

# Simulate a new dataset of abundance at x-y coordinates

# plot a 2D map

# fit a model with interaction between x and y

Four hundred twelve new packages were submitted to CRAN in June. Here are my Top 40 picks in nineteen categories: Bioarchaeology, Biology, Climate Studies, Computational Methods, Ecology, Epidemiology, Finance, Functional Data Analysis, Machine Learning, Medical Statistics, Networks, Pharmacokinetics, Probability, Programming, Psychometrics, Risk Analysis, Statistics, Time Series, and Utilities.

This post describes the R version of Python’s GPopt (https://docs.techtonique.net/GPopt/GPopt.html), a package for ‘Bayesian’ optimization of black-box functions and machine learning hyperparameter tuning, using Gaussian Process Regression and other conformalized surrogates. The package is available on GitHub, and through the R universe.

Keep in mind that this package is for Machine Learning hyperparameter tuning: the global minimum won’t always be found, but this isn’t an issue, since it means you aren’t overfitting the training set.

It’s ported the same way as nnetsauce for R was: with uv to create an isolated Python virtual environment containing the Python GPopt package, and reticulate to call into it from R. Every function in this R package is a thin wrapper that returns the underlying Python object; the general rule is: object accesses with .’s in Python are replaced by $’s in R.

See this post for the technique: Finally figured out a way to port python packages to R using uv and reticulate.

Install

1. Create a Python virtual environment with uv

# pip install uv # if necessary
uv venv venv
source venv/bin/activate      # on Windows: venv\Scripts\activate
uv pip install pip GPopt

Keep track of where venv/ lives – you’ll pass its path as venv_path to every function in this package.

2. Install the R package

install.packages("remotes")
remotes::install_github("Techtonique/GPopt_r") 

reticulate will be installed automatically as a dependency.

Examples

Minimizing the Branin function

This is a standard test function for optimization algorithms. GPOpt is more suitable for expensive black-box functions, but this is a good example to illustrate the usage of the package.

library(GPopt)

branin <- function(x) {
  x1 <- x[1]; x2 <- x[2]
  term1 <- (x2 - (5.1 * x1^2) / (4 * pi^2) + (5 * x1) / pi - 6)^2
  term2 <- 10 * (1 - 1 / (8 * pi)) * cos(x1)
  term1 + term2 + 10
}

opt <- GPOpt(
  lower_bound = c(-5, 0),
  upper_bound = c(10, 15),
  objective_func = branin,
  n_init = 10,
  n_iter = 40,
  venv_path = "./venv"
)

opt$optimize(verbose = 1L)
print(opt$x_min)  # best parameters
print(opt$y_min)  # best objective value

Tuning a scikit-learn model’s hyperparameters

library(GPopt)

sklearn <- get_sklearn(venv_path = "./venv")
RandomForestClassifier <- sklearn$ensemble$RandomForestClassifier

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

mlopt <- MLOptimizer(scoring = "accuracy", cv = 5, venv_path = "./venv")

param_config <- list(
  n_estimators = list(bounds = c(10, 300), dtype = "int"),
  max_depth    = list(bounds = c(1, 20),   dtype = "int")
)

mlopt$optimize(
  X_train = X, y_train = y,
  estimator_class = RandomForestClassifier(),
  param_config = param_config,
  verbose = 1L
)

print(mlopt$get_best_parameters())
print(mlopt$get_best_score())

Bayesian optimization with early stopping

library(GPopt)

opt <- BOstopping(
  f = branin,
  bounds = rbind(c(-5, 10), c(0, 15)),
  venv_path = "./venv"
)
result <- opt$optimize(n_iter = 100L)

Using a custom (conformalized) surrogate model

library(GPopt)

sklearn <- get_sklearn(venv_path = "./venv")
ns <- get_nnetsauce(venv_path = "./venv")

opt <- GPOpt(
  lower_bound = c(-5, 0),
  upper_bound = c(10, 15),
  objective_func = branin,
  acquisition="ucb",
  method="splitconformal",
  surrogate_obj = ns$PredictionInterval(sklearn$ensemble$RandomForestRegressor()),
  venv_path = "./venv"
)
opt$optimize(verbose = 1L)

Biological phenomena are often studied by examining how a numerical variable, usually called the response (e.g., the weight or height of an organism), is affected by another variable, usually called the predictor (e.g., time or fertiliser application rate). These ‘predictor-response’ relationships are commonly described using models, which express the response as a mathematical function of the predictor. In general, a model can be written as:

where is the response, is the predictor, and is a collection of parameters, usually denoted by symbols, such as the letters of the Roman or Greek alphabet. The element is the function, which determines the shape of the relationship when plotted on an x–y graph.

Because mathematical modelling plays such a central role in biology and many other scientific disciplines, biologists need to be familiar with the most important mathematical functions. More importantly, they need to be able to “read” these functions and use their parameters to describe, interpret, and quantify biological processes. With this aim in mind, I have compiled a collection of the mathematical functions most commonly encountered in biology, explaining the meaning of their parameters, with particular emphasis on their biological interpretation rather than their mathematical properties.

Curve shapes

Functions are often classified according to the shape they exhibit when plotted on an x–y graph. This approach is adopted, for example, in Ratkowsky (1990), and the classification presented below is largely based on that work.

• Polynomials
  1. Straight line function
  2. Quadratic polynomial function
• Concave/Convex curves (no inflection)
  1. Exponential function
  2. Asymptotic function / Negative exponential function
  3. Power function
  4. Logarithmic function
  5. Rectangular hyperbola
• Sigmoidal curves
  1. Logistic function
  2. Gompertz function
  3. Modified Gompertz function
  4. Log-logistic function
  5. Weibull (type-1) function
  6. Weibull (type-2) function
• Curves with maxima/minima
  1. Peaked sigmoidal function
  2. Bragg function
  3. Lorentz function
  4. Beta function

I have chosen to include a relatively large number of functions, so this post is necessarily rather long. To make navigation easier, you can first inspect Figure 1 to identify the type of curve of interest and then use the links above to jump directly to the corresponding section.

Throughout this post, denotes the predictor, denotes the response, and the letters of the Roman alphabet (, , , …) denote the model parameters, which determine the shape of the function.

D(expr = expression(a + b*X + c*X^2), name = "X")

b + c * (2 * X)

Concave/Convex curves (no inflection points)

Exponential function

The most common parameterisation of an exponential function is:

Depending on the sign of , Eq. 3 describes a concave-up monotonically increasing shape (exponential growth; ) or a concave-up monotonically decreasing shape (exponential decay; ).

The parameter represents the value of when , whereas the meaning of can be understood by considering the first derivative:

D(expression(a * exp(k * X)), "X")

a * (exp(k * X) * k)

In algebraic terms, it is:

and thus:

Therefore, represents the relative rate of change (increase/decrease), throughout the domain, which is often known in growth analysis as the Relative Growth Rate (RGR).

Other equivalent parameterisations are:

and:

Equations 3-5 are equivalent, as can be shown by setting and :

Another, slightly different, parameterisation is commonly used in bioassay studies, mainly to describe exponential decay processes:

where corresponds to in the previous models and .

For all the exponential decay equations presented above, as . A lower asymptote, , can be incorporated into Eq. 6 to account for situations in which the response does not approach zero as the predictor tends to infinity:

Exponential functions are often used to describe the growth of populations under non-limiting environmental conditions or the degradation of xenobiotics in the environment (first-order degradation kinetics). In both cases, represents time and is therefore restricted to non-negative values, while (or, equivalently, ) must be different from zero. Check the post at this link for an example of fitting.

Asymptotic function

This function appears in several different parameterisations and is also known as the monomolecular growth function, the Mitscherlich law, or the von Bertalanffy law. Owing to its biological interpretation, the most widely used parameterisation is:

It describes a monotonically increasing, concave-down curve approaching a horizontal asymptote as tends to infinity. The parameter represents the maximum attainable value of (the plateau), while is the value of when (the initial value). The parameter is proportional to the rate at which the response approaches the plateau. Indeed, the first derivative is

D(expression(a - (a - b) * exp (- m * X)), "X")

(a - b) * (exp(-m * X) * m)

Considering Eq.8, we can write:

and then:

The previous expression shows that the relative rate of change in the response (the RGR in growth analysis) is not constant, as it is for the exponential function. Instead, it depends on the attained value of . In particular, the RGR is greatest at the beginning of the process, when is smallest, and gradually approaches zero as approaches the plateau .

Another closely related parameterisation is often encountered, in which is replaced by , by , and by . Some simple algebraic manipulations are also introduced to make the model consistent with the parameterisation commonly adopted for biological assay models:

In some applications, is reparameterised through its logarithm to facilitate model fitting:

For all these equations, setting (or equivalently ) yields a curve passing through the origin. In this case, the most common parameterisation is:

which is known as the negative exponential function.

The asymptotic function is widely used in biology, not only for growth analysis but also, for example, in weed competition studies. The negative exponential function has also been used to model the absorbed Photosynthetically Active Radiation () as a function of leaf area index (). In this context, represents the incident PAR (), while is the light extinction coefficient. Check the post at this link for an example of fitting.

Power function

The power function is also known as Freundlich function or allometric function. The most common parameterisation is:

Its shape is highly flexible and is determined by the value of the parameter (see Figure 1). When , the response increases with and the curve is concave down. When , decreases as increases and the curve is concave up. Finally, when , increases with and the curve is concave up. These three cases are illustrated in Figure 1 using different colours. The function is defined only for and has no horizontal asymptote as .

The biological interpretation of the parameters and is not straightforward. Both influence the slope of the curve, as can be seen by examining the first derivative:

D(expression(a * X^b), "X")

a * (X^(b - 1) * b)

The power function is mathematically equivalent to an exponential function of , since:

The power function (named as the Freundlich equation) is widely used in agricultural chemistry, for example to model the sorption of xenobiotics in soil. It is also used to describe the relationship between the number of plant species and the sampled area (the species–area relationship). Check the post at this link for an example of fitting.

Logarithmic function

This function is linear in :

Because of the logarithmic transformation, the function is only defined for . The parameter determines the shape of the curve. When , the response increases with and the curve is concave down. Conversely, when , decreases as increases and the curve is concave up, as shown in Figure 1.

The interpretation of the parameters is fairly straightforward. The parameter is the value of the response when . Changing while keeping constant simply shifts the curve vertically without altering its shape, so that the difference between any two curves remains constant for every value of (Fig. 2).

The parameter controls the slope of the curve and, more specifically, it is equal to the slope at , as can be seen from the expression for the first derivative:

D(expression(a + b*log(X)), "X")

b * (1/X)

Changing while keeping constant produces curves that all intersect at (Fig. 2).

In biology, logarithmic functions are used, for example, to describe species–area relationships in ecology, enzyme kinetics in biochemistry, and sensory perception in neurobiology. Check the post at this link for an example of fitting.

Rectangular hyperbola

The rectangular hyperbola is commonly known as the Michaelis-Menten function and is usually parameterised as

The curve is monotonically increasing and concave down, approaching the horizontal asymptote (plateau) . It passes through the origin of the axes ( when ). The parameter is the value of that produces a response equal to . Indeed,

which is readily solved to obtain .

The first derivative is:

D(expression( (a*X) / (b + X) ), "X")

a/(b + X) - (a * X)/(b + X)^2

From this expression, we can see that the initial slope (at ) is .

Equation 14 is not defined for and has no biological meaning for . In most biological applications, however, both and are positive.

An equivalent parameterisation is:

In this form, the parameters and are often replaced by and , respectively. Although the response is not defined for , it tends to zero as .

Another common parameterisation includes the initial slope as an explicit parameter because of its biological relevance. This is obtained by dividing both the numerator and denominator by and noting that , so that :

The Michaelis-Menten function is widely used in pesticide chemistry, enzyme kinetics (Eq. 14), biological assay models (Eq. 15), and weed competition studies, where it is used to describe crop yield loss as a function of weed density (Eq. 16). An example of fitting is shown in another post, at this link.

Sigmoidal functions

Sigmoidal functions are S-shaped, with two horizontal asymptotes (a lower and an upper asymptote) and an inflection point. They can be parameterised in countless ways, which may be confusing when selecting the most appropriate model for a particular biological process. In this post, I will use the parameterisations proposed in Ritz et al (2019), which is very consistent and puts all the curves on the same ground. Furthermore, the fitting properties are good, which contribute to good convergence and reliable estimation.

For the sake of simplicity, it is useful to recognise that most sigmoidal models belong to one of three basic families, which differ only in the position of the inflection point relative to the two asymptotes:

1. Logistic: the ordinate of the inflection point lies halfway between the lower and upper asymptotes; consequently, the S-shape is symmetric (the concave-up and the concave-down portions ‘mirror’ each other).
2. Gompertz: the ordinate of the inflection point lies closer to the lower asymptote, and, thus, the concave-up portions is ‘shorter’ than the concave-down portion.
3. Modified Gompertz: the ordinate of the inflection point lies closer to the upper asymptote, and, thus, the concave-up portion is ‘longer’ than than the concave-down portion.

Each of these three equations can describe either increasing or decreasing responses, depending on the sign of one of the model parameters (see below). Their different shapes are compared side by side in Fig. 3.

Each of these three families (logistic, Gompertz and modified Gompertz) can be transformed by replacing with , giving raise to three corresponding families:

1. Log-logistic
2. Type-I Weibull
3. Type-II Weibull

Like their linear-scale counterparts, these three families can describe either increasing or decreasing responses, depending on the sign of the same model parameter (see below). In the end, we have twelve different sigmoidal functions to choose from. The choice among them depends primarily on the expected symmetry of the response and on whether the predictor is more naturally interpreted on a linear or logarithmic scale.

Thanks to their versatility, sigmoidal functions have been used for an uncountable number of biological applications. Just to mention a few examples in agriculture, I would like to cite plant growth, dose-response curves in biological assays and the time-course of seed germination. For some examples, take a look at my post at this link.

Logistic function

The logistic curve is derived from the cumulative logistic distribution function. It is symmetric about its inflection point and can be parameterised as:

where is the upper asymptote, is the lower asymptote, is the value of at the inflection point, and is proportional to the slope at the inflection point. Because the curve is symmetric, also represents the value of that produces a response halfway between and (commonly referred to as the ED50 in biological assays). The parameter may be either positive or negative and, consequently, the response may either increase or decrease as increases.

This equation is known as the four-parameter logistic model. If appropriate, constraints can be imposed on the parameter values. For example, may be fixed at 0, yielding the three-parameter logistic model. If, in addition, is fixed at 1, the model reduces to the two-parameter logistic model.

Gompertz function

The Gompertz curve can be parameterised in many different ways. I prefer a parameterisation that closely resembles that of the logistic function:

Unlike the logistic function, the Gompertz curve is not symmetric about its inflection point. It exhibits a longer lag phase at the beginning, followed by a progressively steeper increase before gradually approaching the upper asymptote. The parameters have essentially the same interpretation as those of the logistic function, except that , the abscissa of the inflection point, does not correspond to the value of producing a response halfway between and .

As with the logistic function, four-, three-, and two-parameter Gompertz models can be obtained by constraining one or both asymptotes.

Modified Gompertz function

We have seen that the logistic curve is symmetric about its inflection point, whereas the Gompertz curve exhibits a longer lag phase at the beginning, followed by a progressively steeper increase. A different asymmetric pattern can be obtained by modifying the Gompertz function as follows:

The resulting curve increases rapidly at the beginning but gradually slows down as it approaches the upper asymptote. As with the logistic and Gompertz functions, one or both asymptotes can be constrained ( and/or ), giving rise to four-, three-, and two-parameter versions of the modified Gompertz model.

Log-logistic function

The log-logistic curve is symmetric with respect to . A log-normal curve has a very similar shape, although it is used much less frequently. In biological assays (and also in germination studies), the log-logistic function is commonly parameterised as

The parameters have the same interpretation as in the logistic function. In particular, represents the value of that produces a response halfway between and (the ED50). It is easy to show that the above equation is equivalent to

Like the logistic function, the log-logistic model can be fitted in four-, three-, or two-parameter versions by constraining one or both asymptotes. It is widely used to describe dose-response relationships in biological assays, seed germination, and crop growth.

Weibull function (type 1)

The Type I Weibull function is the logarithmic counterpart of the Gompertz function, being defined on rather than on . It is parameterised as

The parameters have essentially the same interpretation as those of the other sigmoidal functions presented above. In particular, and are the lower and upper asymptotes, respectively, while is the value of corresponding to the inflection point. Unlike the log-logistic function, however, does not correspond to the ED50.

Weibull function (type 2)

The Type II Weibull function is closely related to the Type I Weibull function, but it describes a different type of asymmetry, analogous to that of the modified Gompertz function:

The parameters have the same interpretation as those of the Type I Weibull function. In particular, and are the lower and upper asymptotes, respectively, while is the value of at the inflection point. As with the Type I Weibull function, does not correspond to the ED50. One or both asymptotes may be constrained, giving rise to four-, three-, and two-parameter versions of the model.

Another ‘flexible’ sigmoid

I would also like to mention another sigmoidal function that has been widely used in biology because of its supposed flexibility: the Richards function:

In this model, the degree of asymmetry is controlled by the parameter . When , the function reduces to the logistic model, whereas, as decreases below 1, the curve progressively resembles a Gompertz function. Figure 4 illustrates the different shapes obtained for different values of .

Despite this apparent flexibility, the Richards function has poor statistical properties for parameter estimation. In particular, the additional parameter often causes strong correlations among the estimates, making the fitting process unstable and the parameter estimates difficult to interpret. Ratkowsky (1990) described it as having “more undesirable nonlinear regression behaviour than almost any nonlinear regression model in common use.” For these reasons, the Richards function will not be considered further in this blog.

      As a kid, my baseball hero was Mickey Mantle. In his 18-year Major League career, he amassed tremendous statistics, and surely he was one of the game’s all-time greats. When Mantle retired, he was third on the all-time career home run list with 536 home runs, trailing only Babe Ruth (714) and Willie Mays (587). Unfortunately, he suffered numerous injuries, leaving fans to wonder what his numbers would look like if he had enjoyed a healthier career.

      For decades many people considered (and still consider) Babe Ruth to be the greatest all-around baseball player, because he was a star pitcher and a star hitter. The counter-argument to the Mantle “what-if” is: what if Ruth had not spent his first six years as a pitcher and instead had been a full-time hitter? (Similarly, what if Ted Williams had not spent some of his prime years in the military.)

      Nowadays Shohei Ohtani is certainly a once-in-a-generation two-way star. However, he spent his early prime years in Japan’s professional baseball, and he may not have enough time left in US baseball to amass these historic lifetime stats. Aaron Judge is another exceptional modern hitter; hopefully, he stays healthy so we can watch his career numbers climb.

      To compare Mantle and Ruth’s lifetime stats objectively, I turned to data. The Sean Lahman baseball dataset contains Major League player stats back to 1871, and is available in the R library lahman. This was a good opportunity to practice data manipulation using dplyr.

      Ruth played 102 more games than Mantle over a 22-year career (2,503 versus 2,401). His lifetime statistics eclipse Mantle’s in every offensive category except stolen bases — though Ruth did steal home 10 times and hit nearly twice as many triples, suggesting he was faster than most realize. His lifetime batting average was an impressive .342, sitting just behind Tris Speaker and Ted Williams.

      Lineup protection played a massive role for both men. Ruth generally batted directly before Lou Gehrig, and Mantle batted directly before Yogi Berra. While both Ruth and Mantle drew plenty of walks, presumably pitchers rarely chose to walk them intentionally just to face Gehrig or Berra.

      I was also curious about their defensive metrics. For a fair comparison, I filtered the data to isolate only their outfield appearances. There is no value in comparing an outfield throw to an assist Ruth made on a comeback grounder while pitching. Similarly, I excluded Mantle’s infield appearances; he primarily played first base in his final two seasons and filled in briefly at other infield spots early in his career.

      Ruth played 222 more outfield games than Mantle (2,241 versus 2,019). Their total outfield putouts were nearly identical, likely because Mantle played centerfield and covered more ground. However, Ruth recorded nearly twice as many outfield assists (204 versus 117), which aligns with the arm strength expected of a former pitcher. While Ruth’s 204 assists don’t match stars such as Roberto Clemente’s 266, it remains a highly respectable number.

      Ultimately, the data shows that even if you completely ignore his pitching career, Babe Ruth built an outstanding, standalone career as both a hitter and a fielder.

              
          MANTLE  RUTH
games       2401  2503
at_bats     8102  8398
runs        1677  2174
hits        2415  2873
doubles      344   506
triples       72   136
home_runs    536   714
BA         0.298 0.342
rbi         1509  2217
sb           153   123
bb          1733  2062
so          1710  1330
games_of    2019  2241
putouts     4438  4444
assists      117   204

      Here is my R code:

library(Lahman)
library(tidyverse)

data(People)
which(People$nameLast == "Mantle")   # 11773,  playerID = mantlmi01
which(People$nameLast == "Ruth")   # 16674 , playerID = ruthba01 
data(Batting)    
data(Fielding)   # for fielding, want games in outfield POS == 'OF'

df_fielding <- Fielding %>% 
  filter(playerID %in% c('mantlmi01', 'ruthba01'), POS == 'OF') %>% 
  group_by(playerID) %>% 
  summarize(games_of = sum(G, na.rm = TRUE), putouts = sum(PO, na.rm = TRUE),
    assists = sum(A, na.rm = TRUE)
  )

df <- Batting %>% 
  filter(playerID %in% c('mantlmi01', 'ruthba01')) %>% 
   group_by(playerID) %>% 
   summarize(games = sum(G, na.rm = TRUE), at_bats = sum(AB, na.rm = TRUE),
                runs = sum(R, na.rm = TRUE), hits = sum(H, na.rm = TRUE),
                doubles = sum(X2B, na.rm = TRUE), 
                triples = sum(X3B, na.rm = TRUE), home_runs = sum(HR, na.rm = TRUE),
                BA = round(hits / at_bats,3), 
                rbi = sum(RBI, na.rm = TRUE),  sb = sum(SB, na.rm = TRUE),
                bb = sum(BB, na.rm = TRUE),  so = sum(SO, na.rm = TRUE),
                )  %>%
  merge(df_fielding, by = "playerID", all.x = TRUE) 

df_transposed <- as.data.frame(t(df))
colnames(df_transposed) <- c("MANTLE", "RUTH")
df_transposed <- df_transposed[-1, ]
df_transposed

End

Remember the “AI rat penis” incident at Frontiers in Cell & Developmental Biology?

Well, that was two and a half years ago. That’s enough water under the bridge to now check in and see how this incident affected the journal. If we look at how many papers published, we should be able to see whether the incident affected submissions.

I will say upfront that this analysis did not reveal what I thought it would. It’s well known that Frontiers journals had huge growth some years ago and I also knew that submissions had decreased in the last few years. So, I expected that this drop off was due to authors not wanting to publish in the journal after the AI figure debacle or perhaps a combination of that and the study by Mark Hanson and colleagues showing suspiciously fast turnaround times at journals, including titles from this publisher.

The plots can be quickly generated using {PubMedLagR} which is available here.

Jump to the code, or just look at the plots.

If we pull all the articles in PubMed for this journal and see how many were published per month, we get a plot like this. The dotted line shows the date of the rat penis incident.

So there was a huge boom in publications in the COVID-19 pandemic era and then a decrease. This decrease predates the rat penis incident. If anything, the number of papers published began to stabilise and perhaps even recover following the incident!

It’s a bit hard to see if the composition of articles has changed over time from the plot above. So let’s replot these data and look at the fraction of all articles that are papers versus other types. I did some spot-checking and there were several papers that are clearly Commentary, Editorial and Review types that are not properly tagged on PubMed. So take this with a pinch of salt. The articles are split 50:50 between papers and other types and this has not really changed. It could be argued that commissioned content might be unaffected whereas as directly submitted papers might drop off if authors were concerned about a journal. But there’s no evidence of that here.

This is looking at published papers. Obviously this metric lags behind the author behaviour we’re interested in. Did the rat penis incident affect submissions to the journal? We can look at the received dates of published papers:

Pretty much the same story: month-on-moth submission are steadily rising after the incident.

Obviously we only have submission data for papers that were eventually published in the journal. We don’t know what the rejection rate was over this period. It would be possible for a beleaguered journal to reduce the rejection rate to stabilise the numbers of papers it publishes. But if we assume the rejection rate is constant, then submissions follow similar dynamics.

Conclusion: the AI debacle didn’t negatively affect submissions at this journal.

Maybe potential authors were unaware of this issue (although there was plenty of publicity), or maybe they view the incident as a one-off, or perhaps they’re not concerned about science integrity issues.

It’s no secret that the biggest drivers of journal choice for authors are 1) the journal impact factor and 2) the turnaround time/hassle to publish the work. The journals with the highest impact factor and the fastest turnaround win big. A scandal like this doesn’t seem to have affected journal choice. Of course, this is just a simple analysis and there are many factors at play here.

Isn’t it ironic though that authors want to publish in a journals with high impact factors? They want their paper to be viewed as potentially highly citable because of where it is published; but if that means that their paper appeared in a journal that has a scientific integrity issue, the cry is “judge the paper on its own merits, not where it’s published”…

The code

library(PubMedLagR)
library(ggplot2)

retrieve_journal_year_records("Front Cell Dev Biol", 2015:2026, batch_size = 200, papers_only = FALSE)
all <- pubmed_xmls_to_df(clean = FALSE)
# remove duplicate rows
all <- all[!duplicated(all),]
# clear Data/ and then
retrieve_journal_year_records("Front Cell Dev Biol", 2015:2026, batch_size = 200)
pprs <- pubmed_xmls_to_df()
# classify the papers in all
all$paper <- ifelse(all$pmid %in% unique(pprs$pmid),"paper", "other")
# convert publication date to the first of the month
all$year_month <- as.Date(paste0(substr(all$pubdate,1,7),"-01"))

# frontiers colours for fun
x <- c("242 130 91", "24 157 196")
fcolours <- sapply(strsplit(x, " "), function(x)
  rgb(x[1], x[2], x[3], maxColorValue=255))

# ggplot of number of papers per month-year
p1 <- ggplot(all, aes(x = year_month, fill = paper)) +
  geom_bar() +
  scale_fill_manual(values = fcolours) +
  geom_vline(xintercept = as.POSIXct(as.Date("2024-02-14")), linetype=2) +
  labs(x = "", y = "Papers per month") +
  theme_classic()
p1

# stacked plot version
p2 <- ggplot(all, aes(x = year_month, fill = paper)) +
  geom_bar(position = "fill") +
  scale_fill_manual(values = fcolours) +
  geom_vline(xintercept = as.POSIXct(as.Date("2024-02-14")), linetype=2) +
  labs(x = "", y = "Composition") +
  theme_classic()
p2

# what about submissions?
all$sub_year_month <- as.Date(paste0(substr(all$recdate,1,7),"-01"))
p3 <- ggplot(all, aes(x = sub_year_month, fill = paper)) +
  geom_bar() +
  scale_fill_manual(values = fcolours) +
  geom_vline(xintercept = as.POSIXct(as.Date("2024-02-14")), linetype=2) +
  labs(x = "", y = "Papers per month") +
  theme_classic()
p3

I added some styling to the plots to post them here, but the code will give the same three plots as in the post.

—

The post title comes from “New Frontier” by Donald Fagen from his 1982 “The Nightfly” LP.