The North American Bioconductor Conference 2026 (([BioC2026])(https://bioc2026.bioconductor.org/)) took place from August 10-12, 2026, at the Fred Hutchinson Cancer Center in Seattle, Washington. The conference brought together researchers, developers, educators, and community members to explore current developments within and beyond the Bioconductor project. BioC2026 welcomed 108 attendees, with 70 participating in person and 38 joining virtually. Across the three days, the programme featured three keynote speakers, 20 contributed speakers, and 17 sessions. The programme highlighted the breadth of biological research supported by Bioconductor, as well as the tools, infrastructure, and communities that enable open and reproducible computational biology. The conference was then followed by two community events on August 13-14: a Carpentries workshop on orchestrating large-scale single-cell analysis with Bioconductor and the Bioconductor North America 2026 Hackathon.

Programme overview

Keynotes

The conference featured three keynote speakers: – Jeff Leek speaking on “Building Open Infrastructure for Translational AI in a Cancer Center” – Ting Ye speaking on “From Association to Causation: Genetic-Anchored Causal Inference in Human Biology” – Michael Lawrence speaking on “S7 for Bioconductor”

(Left to Right: Jeff Leek, Michael Lawrence and Ting Ye giving their keynotes at BioC2026)

Short talks

Twenty contributed speakers shared their work across the conference programme. The short talks at BioC2026 reflected the breadth of research and development across the Bioconductor community. Topics ranged from community training and developer engagement to multi-omics exploration, reproducible workflows, and AI infrastructure. The programme also covered emerging areas of biological data analysis, including spatial transcriptomics, antibody repertoire analysis, rare variant testing, DNA methylation, structural variation, and Perturb-seq.

Poster sessions and package demo sessions

The poster session showcased a diverse range of research applications and software development projects from across the Bioconductor community. Topics included single-cell and spatial transcriptomics, multi-omics integration, microbiome analysis, proteomics, antibody repertoire analysis, functional enrichment, ancestry inference, genomic data processing, and AI-guided drug discovery. The posters provided an opportunity for attendees to discuss new methods and tools in a more informal setting, while connecting researchers, package developers, and community members working across different areas of computational biology. The sessions also complemented the package demonstrations, creating space for participants to explore new software and engage directly with the people developing and applying Bioconductor tools.

Feedback from the conference showed that these interactions were particularly valuable, with attendees highlighting conversations around posters and package demonstrations as some of the most useful networking opportunities.

Participants during the poster session.

Celebrating 25 years of Bioconductor

The celebrations marking Bioconductor’s 25th anniversary continued at BioC2026 following the celebrations that began at EuroBioC2026 earlier in the year. The milestone was marked with a specially branded anniversary cake, thanking participants and community members for being part of the Bioconductor community and contributing to its growth over the past 25 years.

Bioconductor at 25 celebrations.

Post-conference events

Orchestrating Large-Scale Single-Cell Analysis with Bioconductor

On August 13–14, participants took part in a Carpentries workshop focused on orchestrating large-scale single-cell analysis with Bioconductor led by Ted Laderas and Wes Wilson. The workshop was attended by 11 participants and provided a hands-on environment for participants to work through approaches for analysing large-scale single-cell datasets using Bioconductor tools and workflows. Feedback highlighted the value of the live coding, small-group format and complementary teaching approaches.

Participants during the scRNA-seq workshop.

Bioconductor North America 2026 Hackathon

The Bioconductor North America 2026 Hackathon took place on August 13–14, bringing community members together to collaborate on software and infrastructure projects that strengthen interoperability, software development, and the wider Bioconductor ecosystem. These projects included: AI harness tools for standardizing how users within organizations would interact with agents on shared projects and on shared institutional resources. Link: https://github.com/eisenefaust/harness_agent

Bioconductor tooling for including CLI executable components in Bioconductor packages. Link: https://blog.bioconductor.org/posts/2026-08-21-biocjobs/

An exploration on benchmarking LLM generated workflows, building out infrastructure and paradigms for how we can evaluate what LLMs return as work product to scientific questions. Link: https://github.com/Amisor/LLM_BioBenchmarking/tree/main

Participants during the Hackathon.

© 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

You can read the original post in its original format on Rtask website by ThinkR here: My first attempt at porting a Shiny app failed, and not for a technical reason

Diary of a rewrite, episode 2 of 4.

Translating your code line by line is the surest way to get it all wrong…

In the previous episode we told the story of why we cancelled our own stack decision halfway through, and why we went back to R. If you have not read it, no worries: this one stands on its own.

Today we are talking about the port itself. And above all about the way I got it wrong the first time.

A quick bit of context. StaffPuzzle is our internal staffing tool: it says who works on what, half-day by half-day. Its Shiny version is 1,200 lines, reads the team’s Google calendars, and has been running for years. The goal was to rewrite it in htmxr, meaning HTML rendered server-side, hydrated by htmx, with plumber2 behind it.

Here is how I went about it.

I opened app_server.R, I listed every observeEvent, and I started turning them into routes. One route per observer. It is methodical, it is exhaustive… and it produces an application full of endpoints that serve no purpose.

A quick reminder: what is a route? In a classic web application, the browser asks for an address and the server sends back a response. A route is that pair: an address, and the function that builds what we send back. Where Shiny holds a continuous conversation with the browser for the whole visit, a route-based application answers requests that are independent of one another, with no memory of the previous one.

Concretely, here is the kind of translation I was doing. On the Shiny side, the observer that rebuilds the main plot whenever a filter moves (I have shortened and renamed it for readability, the original is bushier):

observeEvent(c(input$refresh, rv$slots, rv$data_ready), {
  req(rv$slots, rv$consultants, rv$start, rv$end, rv$sort_by)
  rv$filtered_slots <- rv$slots |>
    filter(member %in% rv$consultants) |>
    mark_focus(types = rv$slot_types) |>
    add_weekdays()
  rv$plot <- rv$filtered_slots |>
    team_availability_plot(
      start   = rv$start,
      end     = rv$end,
      sort_by = rv$sort_by,
      order   = rv$consultants
    )
})

And, a hundred lines further down in the same file, the output that displays it:

output$plot <- renderPlotly({
  rv$plot |> ggplotly(tooltip = "text")
})

On the routes side, it becomes an address and a function:

#* @get /plot
function(consultants, start, end, sort_by, slot_types) {
  slots |>
    filter(member %in% consultants) |>
    mark_focus(types = slot_types) |>
    add_weekdays() |>
    team_availability_plot(
      start   = start,
      end     = end,
      sort_by = sort_by,
      order   = consultants
    )
}

The browser will ask for /plot?consultants=arthur,vincent&start=2026-10-05&end=2026-11-30&sort_by=alphabetical, and the server will answer with the corresponding piece of page.

Notice in passing what has melted away on the Shiny side: I need an observer, a slot in reactiveValues to store the result, and a render* elsewhere in the file to display it. On the route side, all that is left is a function that takes arguments and returns something. That one translates well: it produces something visible, and that something can state its address.

The problem is all the others.

My mistake did not come from not knowing the destination tool. I was translating code instead of changing paradigm. And it took me several days to notice!

It is a bit like moving house. You do not pack by photographing each room so you can rebuild it identically in the new flat. You look at what you actually use, and you leave the rest at the car boot sale.

So here is what I wish someone had told me before I started.

Why inventorying the observers is the wrong entry point

The reason is structural. Shiny is driven by a reactive graph, and a good share of the observers have no visible output at all. They copy a value from one place to another, they disable a button, they keep two inputs in sync. In short, they maintain the consistency of an in-memory state.

But that state does not exist in htmxr. Neither does the work of those observers. They do not translate: they disappear.

I remember the sentence I blurted out in a meeting, on 30 June, right in the middle of the port:

I feel like saying “I’ve got this or that observeEvent, so I’ve got this or that endpoint” isn’t ok. Whereas I thought starting from the data to serve is a good entry point.

It has been my working rule ever since.

Start from the data being served

An htmxr application is a set of addressable blocks of data. Every block on screen has a URL that renders it. The controls that drive it (filters, dates, sorting) are the parameters of that URL. We are in a hypermedia frame, no longer in a dependency graph.

The method comes down to 5 moves, and the order matters:

1. Inventory the screens, screenshot by screenshot. Not what the code does: what the user sees.
2. Cut each screen into blocks of data: header, filters, grid, detail panel.
3. Give every block a URL. A block that cannot state its address is a badly cut block.
4. Turn the controls into parameters. ?start=2026-10-01&n_weeks=8, not observers.
5. For each interaction, ask yourself a single question: is this a mutation of persistent business state?

The fifth one is the one I kept missing. Clicking a filter, changing a sort, picking a period: those are GETs, even if “something changes” on screen. Nothing is modified server-side, we are simply asking for another view of the same data.

Tip: if your inventory produces a lot of POSTs, chances are some GETs have slipped in. Go back over them.

60 lines with nothing left to do

Here is the example that tipped my understanding.

My Shiny version devoted 60 lines to keeping the URL and the inputs in sync, in both directions. Good code, by the way: that application had shareable URLs, which is not that common in Shiny.

On the Shiny side, 60 lines. Here are the two ends of it, abridged and renamed:

## Direction 1: the URL arrives, we put the inputs in the right state
observeEvent(session$clientData$url_search, {
  params <- parseQueryString(session$clientData$url_search)
  url_params$consultants <- params$consultants
  url_params$slot_types  <- params$types
  url_params$services    <- params$services
  url_params$period      <- params$period
  if (!is.null(url_params$consultants)) {
    updateSelectInput(
      session  = session,
      inputId  = "consultants",
      selected = str_to_title(
        str_split_1(str_remove_all(url_params$consultants, " "), pattern = ",")
      )
    )
  }
  ## … and the same treatment for period, slot_types and services
})
## Direction 2: an input moves, we rewrite the URL
observeEvent(
  c(input$consultants, input$slot_types, input$services,
    input$period, input$sort_by),
  ignoreInit = TRUE, {
    consultants <- str_to_lower(glue_collapse(input$consultants, sep = ","))
    slot_types  <- str_to_lower(glue_collapse(input$slot_types, sep = ","))
    period      <- glue_collapse(input$period, sep = ",")
    url <- glue("?consultants={consultants}&types={slot_types}&period={period}")
    updateQueryString(session = session, url, mode = "replace")
  })

On the htmxr side, one argument:

hx_set(
  form,
  get      = ".",
  target   = "#trainings-view",
  trigger  = "change",
  push_url = "true"
)

Careful, this is not a matter of concision. Those 60 lines have not been shortened: they have nothing left to do. They were maintaining a mirror between two representations of the same state, the URL on one side, the inputs on the other. When the URL is the state, there is simply no mirror left to hold.

Back to our house move: these are not furniture you carry up to the new flat. There is no room for them.

Want a simple test to see where you stand? If your port produces code that synchronises two things, you have kept the original paradigm. A good port deletes whole categories of code, it does not rewrite them shorter.

What happens to each building block

Once the frame is set, the port becomes mechanical. Let us look at what becomes of the blocks you have in front of you in your server.R.

A module becomes a route and a rendering function. My biggest module was 267 lines: 4 selectInput, one dateRangeInput, 6 observers, and the URL synchronisation. All that is left is an HTML form and some parameters. The ns() disappears entirely, since it only existed to avoid identifier collisions between instances. With no input binding, there is nothing left to name.

input$x becomes an argument of your function. There is no input object any more. The value arrives in the request: as a URL parameter for a GET, in the body for a POST. It is the most disorienting change in the first few days, and the most liberating one afterwards. A view function receives its inputs as arguments, so it can be tested!

output$x <- render*() becomes a function that returns HTML. The output/render pair disappears. All that is left is a function and its return value, called by a route. A renderPlot becomes a URL that serves the SVG, a renderUI a URL that serves a fragment.

updateSelectInput() becomes: nothing. We no longer update a control. The server renders the <select> in the state it should be in, and htmx swaps the fragment that contains it. The whole update*Input() family evaporates with it.

invalidateLater() becomes an attribute. A periodic refresh is written trigger = "every 15s". The browser asks for the fragment again, the server recomputes it. Nothing to invalidate, since there is no state.

showNotification() and withProgress() become a fragment, or CSS. A message is a piece of HTML returned with the response. A waiting indicator is an element that htmx shows during the request and hides afterwards. A CSS class, not R code.

reactiveValues() and session$userData become the URL, or the database. In the URL if it describes what we are looking at, in the database if it has to survive closing the tab. And if nothing is left after that sorting? Then it was plumbing state, and it goes with the rest

Seeing it run

All of this stays abstract until you have it in front of you. So let us take Shiny’s canonical example, the one you get with shiny::runExample("01_hello"): the Old Faithful geyser histogram, with a slider for the number of bins.

• The application, written in htmxr. One slider, one plot, and that is all.
• The URL that serves the plot. Open it directly in your browser: you land on the SVG itself, the drawing spelled out. That is exactly what “give every block a URL” means.

Have fun changing bins=30 to bins=5, then to bins=50. The server recomputes and sends you back the drawing. No session, no state, no open conversation.

And the code, on both sides

You already know the Shiny version: it is the one from the tutorial. Here it is, stripped to the essentials.

ui <- page_sidebar(
  title = "Hello Shiny!",
  sidebar = sidebar(
    sliderInput("bins", "Number of bins:", min = 1, max = 50, value = 30)
  ),
  plotOutput("distPlot")
)
server <- function(input, output) {
  output$distPlot <- renderPlot({
    x    <- faithful$waiting
    bins <- seq(min(x), max(x), length.out = input$bins + 1)
    hist(x, breaks = bins, col = "darkgray", border = "white")
  })
}

And the htmxr version, the one running behind the link above. I have removed the Bootstrap classes, they teach you nothing.

#* @get /
#* @serializer html
function() {
  hx_page(
    hx_head(title = "Old Faithful Geyser Data"),
    hx_slider_input(
      id      = "bins",
      label   = "Number of bins:",
      value   = 30, min = 1, max = 50,
      get     = "plot",
      trigger = "input changed delay:300ms",
      target  = "#plot"
    ),
    tags$div(id = "plot") |>
      hx_set(get = "plot", trigger = "load", target = "#plot")
  )
}
#* @get /plot
#* @serializer none
function(query) {
  xmlSVG({
    x    <- faithful$waiting
    bins <- seq(min(x), max(x), length.out = as.numeric(query$bins %||% 30) + 1)
    hist(x, breaks = bins, col = "darkgray", border = "white")
  }, width = 7, height = 5)
}

Three things to look at, and they sum up the whole article.

The slider is the same object. sliderInput(id, label, min, max, value) becomes hx_slider_input(id, label, min, max, value). There is nothing to relearn.

The plot computation has not moved by a single line. The four lines of hist() are identical on both sides. That is what I had not understood when I started: it is not your business code you are porting, it is the plumbing around it.

What changes is the wiring, and it becomes visible. In Shiny, nothing says that moving the slider redraws the plot: the reactive graph works it out by itself from the presence of input$bins inside renderPlot. In htmxr, it is written on the slider itself: get = "plot" (go fetch this address), target = "#plot" (put the result there). You read the wiring instead of guessing it.

And server <- function(input, output) is gone. All that is left is two functions that take arguments and return something.

Try it yourself

No need to copy anything: the example ships with the package.

# install.packages("pak")
pak::pak("hyperverse-r/htmxr")
htmxr::hx_run_example("hello")

Your browser is waiting at http://127.0.0.1:8080. Move the slider: every move fires a request to /plot, and the server returns a new SVG. Open your browser’s network tab to watch them go by, it beats any explanation.

And if you call hx_run_example() with no argument, it lists the others: select-input, multi-select, delete-row, toast-notification, infinity-scroll, page-or-fragment, json-endpoint and reactive-values. That last one is worth a look if you are wondering what becomes of a reactiveValues() once the reactive graph is gone.

And since the plot is a URL, anyone can consume it. Here is the same application written in React, calling exactly the same R endpoint. The backend does not know who is asking, and it does not care.

This is the most underrated property of cutting things into addressable blocks: it does not lock you in. The day you want another frontend, or a mobile app, or simply for a colleague to pull your plot into their own tool, your R logic is already an API. It became one without you doing anything extra.

What costs, and what does not

Everything above is mechanical. On my port, only 3 patterns called for a real decision, and they are the ones that make an estimate vary.

Shiny pattern	Cost	Equivalent
selectInput, actionButton, sliderInput	low	hx_select_input(), hx_button(), hx_slider_input()
dateRangeInput	low	no dedicated constructor: two <input type="date">
URL ︎ inputs sync	low	the push_url argument
static ggplot	low	svglite + htmx swap
DT (display)	low	hx_table(), hx_table_rows()
long task + notifications	medium	external job + periodic polling
interactive plotly	high	plotly.js, or drop the interactivity
rhandsontable	high	JS library, or one form per cell

The editable table has no native equivalent. Either you embed a JavaScript library, or you rebuild it with HTML fields and one save per cell. The criterion is dimensional: under twenty rows or so HTML wins, at a thousand rows with multi-cell copy-paste the library becomes the right choice again.

The interactive plot deserves a prior question, and it is precisely the one I forgot to ask myself: is it really a plot? When it is in fact a data grid, a schedule, a calendar or a heatmap, porting it as a plot is a mistake. Mine was one because Shiny made tables hard, not because the table was the right object. A real plot, on the other hand, stays an SVG served by a URL, and that is an easy case.

An opaque internal dependency is not worked around, it is isolated. You definitely have one: the homegrown package that fetches the data from somewhere, that nobody has really read in years. The temptation is to start with it, to “understand before porting”. Do not. You will spend three days on it and you will have ported nothing.

Develop everything else against a demo dataset, hard-coded, and only plug the real package in at the very end, in one go. By then you will know exactly what you want from it: a function, some arguments, a table back. And if it is still incomprehensible, call it as it is without trying to rewrite it. That is not cowardice, it is what stops it from blocking the whole port.

The most useful side effect

Shiny encourages mixing business logic and reactive orchestration: they live in the same file, often in the same function. In my v1, the rule that decides whether a slot is highlighted was computed twice, with different rules, 15 lines apart. That was not carelessness on my part, it is what a file where business and plumbing cohabit produces.

The port forces the separation. A function called by a route can no longer lean on an ambient reactive context. What is left is ordinary R, which takes data and returns data or HTML. So it can be tested.

The tests-to-code ratio went from 0.26 to 0.77: 544 tests, no browser, a few milliseconds. With the interface having become a pure function from data to HTML, it falls into the same test harness as the rest. The Shiny equivalent needs shinytest2, a headless Chrome and snapshots. It is slow, it is fragile, and it breaks when a margin moves.

Two pieces of honesty are in order here. My new version does more than the old one (writing to calendars, draft mode, historised indicators), and part of the difference in volume comes from that, not from the paradigm. And nothing stops you from testing the business logic of a Shiny application: v1 did it over 300 lines. What changes is that the interface stops being the area you never test.

What you do not get back

Let us be clear: the reactive graph does render a real service when everything depends on everything. Think of a dashboard where 6 linked plots recompute together, with cross dependencies and expensive derived values. Rebuilding that with fragments and URL parameters is possible, but you are hand-writing what Shiny gives you for free.

Hence the only decision criterion that is worth anything, to my mind: does my application’s state fit in a URL?

For mine, yes: a period, some filters, a sort. The port was a good bet. For an exploration tool where the user builds a complex working session, the answer is no, and Shiny remains the right tool. That is not a polite concession, it is the same reasoning applied to a different problem.

Conclusion

I still write Shiny, for our clients and for ourselves. This port did not prove that one technology beats another.

What it did show me is that for years I had been solving a page-rendering problem with a tool designed to hold state in memory. And that the cost of that gap stayed invisible until I had paid it once.

So before you set off, ask yourself the URL question. It will save you the few days I lost!

Do you have a Shiny application in this situation? We read it and give you back a note that says whether the port is worth it, including when the answer is no. Write to us in two lines: the size of your application, and how long it has been running.

Next week, episode 3: what leaving the managed platform cost us, bill in hand.

Diary of a rewrite, four episodes:

1. We chose React to learn it. The AI wrote it for us.
2. My first attempt at porting a Shiny app failed, and not for a technical reason (you are here)
3. We left the managed platform. Here is the real bill. (6 October)
4. Five signs your Shiny app has outgrown Shiny (13 October)

The application cited is StaffPuzzle, ThinkR’s internal staffing tool, written with htmxr, alpiner and plumber2. The figures are measured on the repository’s two branches, not estimated. Disclosure of interest: I am the author of htmxr and alpiner.

This post is better presented on its original ThinkR website here: My first attempt at porting a Shiny app failed, and not for a technical reason

Three hundred twenty-nine new packages were submitted to CRAN in July. Here are my Top 40 picks in nineteen categories: Agriculture, Artificial Intelligence, Climate Studies, Computational Methods, Ecology, Epidemiology, Genomics, Machine Learning, Medical Statistics, Meta-Analysis, Networks, Process control, Programming, Public Transit, Statistics, Surveys, Time Series, Utilities, and Visualization.

Read it in: Español.

rOpenSci turns 15 and I’ve been a part of this community for a good chunk of that time. As I reflected on those years, a lot of great memories came to mind, and I realized some of them taught me some lessons that I would love to share.

Stretch your comfort zone

Early in 2018 I came across this post:

Apply to attend rOpenSci unconf 2018!
We’re organizing unconf18 to bring together scientists, developers, and open data enthusiasts from academia, industry, government, and non-profits to get together for a couple of days to hack on various projects and generally enrich our community.

In hindsight I was a good fit. I had recently transitioned from academia to my first job as a research software engineer. And yet, did I think I was a solid candidate? Oh no, I had lots of doubts.

But I know the greatest things often happen a little stretch outside your comfort zone. So I applied anyway, got invited, and spent two days collaborating with an incredible bunch of people. Some of them were already rockstars, many others eventually became influential in their own niche, and all of them had this great attitude that at the time I lacked words to describe. I now do: they were respectful and kind.

Be respectful and kind

In 2022 rOpenSci launched the first cohort of the rOpenSci Champions Program. By then I already had a few years of experience as an associate editor of rOpenSci software peer-review, and Yani invited me to talk to our champions about the process.

As I prepared that talk and discussed some ideas with Yani, she made me realize how much the rOpenSci community cares about the way we communicate with one another¹. So much so that our guides for reviewers and editors open with this message:

rOpenSci’s community is our best asset. We aim for reviews to be open, non-adversarial, and focused on improving software quality. Be respectful and kind! See our reviewers’ guide and code of conduct for more.

Coming from such a deeply technical guide, that opening might surprise you. But it makes sense to me now; broken code is way easier to fix than broken human relationships.

Contribute in your own way

At rOpenSci you can contribute in so many ways. We have a community contributing guide but here’s a list of some of my own contributions:

• Attend events.
• Help to welcome and onboard new members.
• Ask or answer questions on Slack, or share or discuss ideas or jobs.
• Improve our documentation, e.g. fix a typo.
• Write a blog post.
• Translate or review some work in your native language.
• Pilot a new process.
• Author, review, or edit an R package.
• Lead a workshop.
• Build a tool to enhance some process.
• Give and get support, e.g. mentor, review grant applications or talks, advice or recommendations for jobs, help unavailable or overwhelmed people.

As you can see, most of my contributions did not involve any code, and only a few required an invitation. Some of them are the kind of thing you could add to your CV, and the most important ones only belong in your heart.

Attract great people

Great communities are made of great people, and you can play an active role in growing it in the direction you want. If you know someone that would fit in and enjoy the rOpenSci community, you can try to find a way to bring them in. For example, I’ve encouraged and helped people to submit and review packages. Also as an editor I have the occasional privilege to nominate other editors. Using this super-power I’ve attracted two editors and I now get to enjoy interacting with them quite regularly.

Keep the human in your human connections

AI is changing the way we contribute to open source software. But we don’t know much about its effect on the open communities around that software. Earlier this year rOpenSci published a preliminary set of policies, and you may want to watch rOpenSci’s newsletter and this Openscapes community call:

Open Communities in the Age of AI
{Open communities} help us connect on a human level around science and data and the things that we are passionate about. They change careers, they change lives. (…) But now, people are not necessarily finding their communities and getting the benefit of those communities in the way that they used to, because they’re more easily able to get help and answers by using AI. (…) These connections and communities are more important than ever, as we need to find ways to support each other and ourselves as AI rapidly changes the landscape right under our feet.

Personally, over the past year I’ve used AI heavily and learned several lessons. The most relevant here is that I don’t want AI to impersonate me or the other person in a human-to-human connection. The humans of rOpenSci are real and wonderful people. Here are some that I’ve recently met in person:

Thanks for joining me in celebrating rOpenSci. And if, like me, you self-identify with our community, then happy birthday to you!

TLDR: In 2019 I wrote a post with seven reasons policy professionals should learn to code. Now that AI can write reasonable code, I’ve had a number of people ask me if it’s still worth learning to program. Based on the evidence I’ve seen, AI appears to be a an amplifier, rather than a substitute for human intelligence in the policy space: meaning the people that get the most out of AI are likely to be those that know how to instruct AI to tackle a problem and tell good output from bad. So my answer is yes, learning to code still makes sense. Not necessarily so you can write the code yourself, but so you can make AI a force multiplier for your work.

Note: I’m intending to make updates and refinements to this post over time: both to correct the typos that I haven’t picked up and to incorporate feedback from other coders in the policy community. Feel free to reach out if you have suggestions.

How AI was used: I wrote the first draft of this post. AI was used (minimally) to refine how some ideas were communicated.

Background

In 2019 I wrote a listicle outlining why people working on public policy should learn how to program.

Putting aside that R python was probably the answer to which language to learn, the point of the post was to share some of the practical benefits I’d seen as an applied economist and policy advisor from learning to code.

Despite the listicle being one of my lower effort posts, it seemed to strike a chord: my inbox was quickly flooded with messages thanking me for articulating the benefits others in the public policy community had seen from leveraging code in their work. To my surprise, the article was also used as a reference for several university programming courses, resulting in my inbox still receiving messages from people starting their coding journey.

But a lot has changed since I first wrote the article: the UK voted to leave the European Union; the World Health Organization declared a global pandemic; and 25 movies in the Marvel Cinematic universe were released.

Oh! And a quaint lil’ startup called OpenAI launched something called ChatGPT.

And while there’s certainly a lot that could be said about the unfolding Marvel Universe’s implication for learning to code, this will have to be a subject of a future post. As today, I’d like to attempt to answer a question I’ve been increasingly asked since Artificial Intelligence (AI) tools and/or Large Language Models (LLMs) have made producing code trivial.

It’s sunk costs all the way down

At the outset, I’m clearly not an impartial observer. I was motivated enough to teach myself R, to write the original post and create a dedicated MOOC on using R for policy analysis. I also enjoy coding, so I have a vested interest in answering yes to the question. But, in my defence, these reasons also mean that I’ve thought about it a lot. I’m also a frequent user of AI in my work and find it has greatly expanded the scope and scale of what I can achieve.

This post is therefore my attempt to provide an answer to the question Is it still worth learning to code in 202X? My answer is, yes, but for different reasons to the 2019 listicle.

Doing smart stuff quickly

What got me excited about learning to program when I first started, was its ability to expand the scope and scale of analysis I could do and questions I could answer. In the world of public policy, there’s an almost unlimited number of interesting questions that could be asked, but rarely the time and data to answer them. Sometimes this is simply because the data didn’t exist, which makes learning to code not particularly helpful. However, often there are cases where data does exist, just not in a format that can be readily analysed with the commercial tools made available to you:

• Want to import thousands of individual excel files and merge them into a single dataframe?
• Have a repetitive data cleaning task that would take hours to do using a spreadsheet?
• Want to apply a statistical technique that hasn’t been added to your stats software?

Programming languages makes solving these problems trivial, but with the right prompt, so can AI.

Probability-based slop machines

But let’s be clear: at their core, LLMs are probabilistic slop machines. I don’t mean slop in the pejorative sense, but as a palpable description of what LLMs do: slap together intelligent looking responses based on next-token prediction.[1] I say intelligent looking[2] as LLMs have been trained to produce outputs that look intelligent, rather than to exhibit intelligence in a human sense:

“…AI systems have demonstrated their capability to solve cognitive ability test problems, primarily through guided training (e.g., Zhuo & Kankanhalli, 2020) or programmed approaches to transform problems into algorithmically solvable formats (e.g., Schmidhuber, 2004). While remarkable, it is debatable whether these accomplishments signify intelligence, given that the capabilities of most current AI systems are limited to specific programming and/or training data, without the necessary demonstration of novel problem-solving ability characteristic of human intelligence (Davidson & Downing, 2000; Raaheim & Brun, 1985). Consequently, many AI systems might be more aptly recognised as having the capacity to exhibit artificial achievement or artificial expertise…”

Gignac, G.E. and Szodorai, E.T., 2024. Defining intelligence: Bridging the gap between human and artificial perspectives. Intelligence, 104, p.101832 (link).

But sometimes a little bit of AI slop can be helpful. When producing code, I’ll often get AI to produce the first draft before making tweaks to the approach and style.[3] I use AI to help with writing and research. For this post I asked AI for good search terms to find research exploring the cognitive and economic benefits of learning to code, so I could check my assumptions and think clearly about what has changed since I first wrote the post. I then asked AI for ideas to better communicate key ideas presented in this post.  

“The majority of human production has always been slop. Mediocrity is not a bug of technology; it is the baseline of culture. The canon of art we revere today – those few thousands works in museums and textbooks – is the surviving tip of an immense iceberg of forgotten, derivative, or simply boring creations…”

D’Isa, Francesco (1 December 2025). “The Idea of ‘AI Slop’ Is Slop”. The Philosophical Salon (link)

If you’re wondering why I’m waxing lyrical about intelligence and the mechanics of LLMs, it’s because I think they help decide two things: which problems AI is likely to be good at solving and whether (or when) it’s worth learning to code. After all, if LLMs can already code better and faster than you, why bother learning how to do it yourself? Particularly if the type of intelligence AI outputs can mimic have sufficient crossover with our own to make it practically indistinguishable in the realm of applied economics and public policy.

Widening the scope, scale and quality(?) of analysis

Looking back at my original post, almost all the cited benefits come from programming languages changing the scope, scale and quality of what you can produce as an applied policy analyst. And if I had to put a number on it: of the seven reasons cited in the listicle, six can be achieved without having to write the code yourself.[4] Have some messy data that needs to be cleaned and reshaped before analysis? Share a sample with AI and let it write the code to take care of the data cleaning for you. Have research question you want answered? LLMs can suggest alternative statistical approaches you might not have thought of and write the code to implement them. Need to produce a set of publication ready plots? AI can write the code to do this too in a format that tells the story of your choice.

It would seem to be game over for learning to code and my 2019 post.

Of course, you already know that I don’t believe this but let me tell you why…

The jagged frontier of AI-assisted public policy analysis

Firstly, whether AI is intelligent or merely produces outputs that make it appear to be, early evidence suggests its ability to substitute for human intelligence varies across tasks. Researchers that tested the productivity benefits of AI termed this as a jagged frontier, that the users had to carefully navigate to benefit from using AI in their work.[5] It’s also not clear whether improving model performance is likely to expand the frontier equally across tasks, with many of the benchmark scores not necessarily generalizing outside test questions and being fragile to irrelevant and/or incorrect context.[6] Finally, at the time of writing this post, the types of tasks that frontier AI models had the hardest time completing successfully, looked strikingly similar to the problems encountered in public policy: messy, interdependent and heavily reliant on context.[7]

Secondly, while there’s solid evidence suggesting that AI can lead to productivity gains for knowledge workers, this isn’t automatic and appears to lean on the taste, judgement and experience of the user. [8] This was demonstrated by a field experiment in Kenya,[9] where a little over 600 entrepreneurs were randomly provided access to an AI-powered business assistant. The performance of those with the assistant was then tracked over time and compared to those without the AI assistant. The result: top performers gained around 15 percent from having access to AI, while weaker ones lost around 10. Entrepreneurs that benefited the most from AI appeared to be better at identifying and acting on contextually appropriate advice, while poor performers tended to adopt generic strategies proposed by the tool.

A much smaller experimental study closer to the world of coding pointed to a similar idea,  with both the 2025 (n=16) and 2026 (n=57) replication suggesting the productivity returns to using AI are ambiguous and vary based on the experience of the user.[10] Both versions of the study took the same basic form: experienced developers were randomly assigned to work with AI or not on their own repositories and the performance benefits were compared between groups. In the 2025 study, the authors found the performance benefits were on average negative: experienced developers took longer to complete the same tasks and overestimated the time saved.  In the 2026 replication,[11] the authors found performance increased on average, but not for all developers.[12] Although the authors note their estimates are likely biased downwards due to difficulties in recruiting developers willing not to use AI, they paint the same basic picture: the benefits of using AI are not unambiguously positive and depend on the skills and experience of the user.[13]  

“…The right analogy for AI is not humans, but an alien intelligence with a distinct set of capabilities and limitations. Just because it exceeds human ability at one task doesn’t mean it can do all related work at human level. Although AIs and humans can perform some similar tasks, the underlying “cognitive” processes are fundamentally different…”

Mollick, E. (2024, May 12). Superhuman? What does it mean for AI to be better than a human? And how can we tell? One Useful Thing (link).

There’s no accounting for taste

To be clear, my point here isn’t to suggest AI tools aren’t helpful, just that the size of the benefit depends on the judgement, expertise and taste of the user. AI certainly can suggest a variety of statistical techniques for answering your research question, but are they the right ones?[14] Sure, AI can quickly produce code for twenty reasonable looking plots, but do you understand the precise changes needed to accurately present your results and apply your corporate style guidelines and colours? Certainly, AI might simplify data cleaning, but does the approach make the right assumptions given the context, your research question and how the data was collected? And, while it might be true that research built on AI-generated code is more reproducible than a spreadsheet with a tangled collection of formulas, it may just be encoding hidden errors, the wrong methodology and creating a new source of technical debt you and your team will have to pay in the future.

In my original post I was excited at programming expanding the scale, scope and quality of analysis I could do. AI can help with the first two, but its influence on the latter seems to depend on the user. AI may be able to slap together a solid draft that gets you 90 percent of the way, which for some tasks is all you need, but deciding what good looks like and instructing AI how to help relies on building sufficient judgement, taste and expertise. 

Vibing fast and slow

The need to develop judgement, expertise and professional taste is also not unique to coding. But like musicians with sheet music, a mathematician with an equation and an economist with a set of national account statistics, anyone working with code is likely to perform better with some understanding of how it works, even with the help of AI. In the words of Andrew Ng:

“…people who understand the language of software through their knowledge of coding can tell an LLM or an AI-enabled IDE what they want much more precisely, and get much better results.”

And once again, there’s research that backs this up. For instance, pre-AI a meta-analysis of the effects of learning to code points to benefits to creativity and metacognition i.e. the processes underlying the monitoring, adaptation, evaluation, and planning of thinking and behaviour,[15] both which are plausibly linked to the ability to properly prompt and manage AI to perform tasks. A small observational study[16]  (n=100) seems to support this idea, with their analysis finding a positive association between vibe-coding efficacy and computer science (CS) domain expertise. Although the exact questions used to evaluate CS expertise aren’t explicitly outlined, their cited source points to them covering the same thinking and skills that learning to program helps to build,[17] with the authors concluding that effectively vibe-coding likely draws on problem decomposition and algorithmic thinking skills.

The book Messy Jobs also makes this point by noting that by cheapening the cost of cognition, AI simultaneously raises the value of judgement and taste.[18] The consequence of this for vibe-coding is that both instructing AI how to solve a problem and being able to judge the quality of outputs relies on knowing how code works, otherwise you might just be doing stupid things faster, in the words of Hadley Wickham:

“…LLMs and AI are, I think, a massive amplifier. You can certainly use them to amplify the worst aspects of yourself, but you can also use them to amplify the best aspects…”

Wickham, H. (2026, July 25). y code when ai? A talk on AI and coding. Tidy Design. https://tidydesign.substack.com/p/y-code-when-ai

Summing up

In my 2019 listicle, I argued that learning to code was worth it so you can answer a wider variety of questions more quickly and at a higher standard to traditional point-and-click software. Seven years later and AI tools can produce code for a wide variety of tasks more quickly than an experienced developer. But that doesn’t make what it produces good.

None of this means that I think you should avoid using AI before having mastered how to code, just that to be a force-multiplier it helps to know what good code looks like. Treat it as a brainy intern: motivated, well read, and in need of supervision. Give it clear instructions and it might come back with something useful. Don’t and its code might set the office microwave on fire. Learning to code can help you tell the difference between the two.

That matters because you’re accountable for the result as the human given the the instructions. Use AI to decide how to prioritize which grants to cancel and you might cause unnecessary harm,[19] trust its recommendations without validating them and you might just erroneously reinforce cultural stereotypes[20], and have it author a report for you, and it might just cite hallucinated sources.[21] If used well, AI can help you reach your destination quicker and expand your professional frontier. If used poorly, it might just having you producing more errors at scale.

I don’t know how long this holds. The studies here are small and not focused directly on the tasks we do in public policy. My best guess is that the jagged frontier will continue to expand, but unevenly, which will continue to place a premium on professional expertise and taste for tasks AI continues to struggle with.

So: should you still learn to code in 202X? Yes, but not for the reason I gave in 2019. Learn to code so you can tell good code from bad and steer the tool in the right direction.

[1]The push-back I’ve seen to this description of how LLMs work doesn’t change the core point. LLMs are trained to probabilistically predict the next token to produce intelligent-looking outputs at scale.

[2] LLMs are trained on the finished product of human intelligence, such as published research, books and social media posts, but not the thinking that went into producing them. Newer models that ‘think’ before responding to a prompt are meant as a partial workaround to this, but don’t change the basic point that LLMs are trained to mimic articulated intelligence, rather than the process of producing. 

[3] A practice that might have the opposite result if you’re not careful, see: Goldsmith-Pinkham, P. (2026, September 10). How do we know if writing is AI? Some initial analyses of AI’s impact on economics writing. A Causal Affair (link) 

[4] The one that may no longer be true is the trend towards open-source statistics. Certainly, I’ve been experimenting with ways to use AI to better extract government data (see example), but it’s entirely possible LLMs will make the internet a more difficult resource to navigate as text can be generated at larger scales.

[5] Dell’Acqua, F., McFowland III, E., Mollick, E., Lifshitz, H., Kellogg, K.C., Rajendran, S., Krayer, L., Candelon, F. and Lakhani, K.R., 2026. Navigating the jagged technological frontier: Field experimental evidence of the effects of artificial intelligence on knowledge worker productivity and quality. Organization Science, 37(2), pp.403-423.

[6] Akimitsu, P., 2026. Wrong and More Confident: A Field Experiment on Large Language Models Taking a Graduate Economics Exam. arXiv preprint arXiv:2607.23424.

[7] METR. (2026, May). Task-completion time horizons of frontier AI models. https://metr.org/time-horizons/

[8] Ethan Mollick uses the term Centaurs and Cyborgs todescribe the two most common approaches used by  workers that were able to navigate this frontier. In both cases, using AI effectively requires intelligently and selectively integrating the tools into their work without ‘falling asleep at the wheel’ (link).

[9] Otis, N.G., Clarke, R., Delecourt, S., Holtz, D. and Koning, R., 2023. The Uneven Impact of Generative AI on Entrepreneurial Performance: Evidence from a Field Experiment in Kenya (No. 24-042). Harvard Business School Working Paper.

[10] Becker, J., Rush, N., Barnes, B., & Rein, D. (2025, July 10). Measuring the impact of early-2025 AI on experienced open-source developer productivity. METR. https://metr.org/blog/2025-07-10-early-2025-ai-experienced-os-dev-study/

[11] Becker, J., Rush, N., Cunningham, T., Rein, D., & Mahamud, K. (2026, February 24). We are changing our developer productivity experiment design. METR. https://metr.org/blog/2026-02-24-uplift-update/

[12] The authors note it had become increasingly in their 2026 follow-up study to recruit developers willing to work without AI, which likely resulted in estimates of AI-assisted speedup being biased downward.

[13] Although the authors recommend their results are interpreted with caution, the 2025 study not only found that developers using AI were slower, but that more experienced developers tended to overestimate the time they would save from using AI.

[14]Like humans, LLMs exhibit systematic biases in the way they work and evaluate problems. This might come in subtle forms that influence outputs subtly, such as assuming doctors are more likely to be men than women. Or explicit, such as producing outputs that reflect how prominent a solution is in its training set, rather than necessarily being appropriate for the task at hand (link).  

[15] Scherer, R., Siddiq, F. and Sánchez Viveros, B., 2019. The cognitive benefits of learning computer programming: A meta-analysis of transfer effects. Journal of Educational Psychology, 111(5), p.764.

[16] Thorgeirsson, S., Weidmann, T.B. and Su, Z., 2026, April. Computer science achievement and writing skills predict vibe coding proficiency. In Proceedings of the 2026 CHI Conference on Human Factors in Computing Systems (pp. 1-17).

[17] Specifically, code tracing, code completion, conditionals, loops and logical operators, see: Parker, M.C., Solomon, A., Pritchett, B., Illingworth, D.A., Marguilieux, L.E. and Guzdial, M., 2018, August. Socioeconomic status and computer science achievement: Spatial ability as a mediating variable in a novel model of understanding. In Proceedings of the 2018 ACM Conference on International Computing Education Research (pp. 97-105).

[18] Garicano, L, Li, J & Wu, Y 2026, Messy Jobs: The Work That AI Cannot Reach, Upriver Press

[19] Alibašić, H., 2026. Introduction: The Imperative of Hybrid Intelligence in Digital Governance. In Hybrid Intelligence for Effective Digital Governance: AI in Administration (pp. 3-48). Cham: Springer Nature Switzerland.

[20]Saar Alon-Barkat, Madalina Busuioc, Human–AI Interactions in Public Sector Decision Making: “Automation Bias” and “Selective Adherence” to Algorithmic Advice, Journal of Public Administration Research and Theory, Volume 33, Issue 1, January 2023, Pages 153–169, https://doi.org/10.1093/jopart/muac007 

[21] Karp, P. (2025, November 6). AI-tainted Deloitte report was worse than previously thought. Australian Financial Review. https://www.afr.com/politics/ai-tainted-deloitte-report-was-worse-than-previously-thought-20251106-p5n863

The post X Reasons for policy professionals to get into programming in 202X appeared first on Giles.

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

rOpenSci HQ

Yani and Mark at “Open Communities in the age of AI”

rOpenSci community manager Yani and software-review lead Mark will participate in the upcoming Openscapes community call on October 1st (this Thursday!) at 9:30AM PT (16:30 UTC) entitled “Open Communities in the Age of AI”, together with Mara Averick, Senior Developer Advocate at Quansight and Hadley Wickham, Chief Scientist at Posit.

Event page, including link for free registration.

Champions Program update

Our current 2026–2027 cohort has completed the training phase and is now focused on developing their individual projects with the support of their mentors, as well as on their outreach activities.

Meanwhile, the 2025–2026 cohort will wrap up their journey with a closing community call, Más Allá del Código: muestra abierta de los proyectos de nuestros campeon(a|e)s (Beyond the Code: an open showcase of our Champions’ Projects). At the call, the Champions will present their final projects, from creating, improving, and reviewing R packages to the outreach work that brought those projects to their communities. They’ll also share what their projects led to, including scholarships, new collaborations, and professional growth. The event is in Spanish. The call will take place on Monday, October 12 at 15:00 UTC. Come get inspired, ask your questions live, and help us celebrate our growing open-source community!

We’re still celebrating our 15th anniversary!

In July, we started to share stories from members of our community about their experiences with rOpenSci. Our second story features Yi-Chin Sunny Tseng and her connection with rOpenSci. Read it on our blog: Happy Birthday rOpenSci — My Journey from First-Time Developer to Current Opportunities Stay tuned for more stories from our community as we continue celebrating 15 years of rOpenSci!

rOpenSci at LatinR 2026 in Medellín, Colombia. See you there!

Registration for LatinR 2026 (November 11–13, Universidad de Antioquia, Medellín) is now open, and rOpenSci will have a strong presence.

Jeroen will give a talk on R-Universe. Yani will lead a workshop on R package development and give a talk on the Champions Program’s open curriculum. Nic Crane is one of the conference’s keynote speakers. Evelia Lorena Coss Navarrete, one of our Champions, will talk about the package she developed during the program. Several other rOpenSci folks — Natalia Da Silva, Luis Verde, and Francisco Cardozo, among others — will also be there. Jeroen, Francisco, and Yani will take part in the hackathon during the conference.

Join us in Medellín!

Coworking

Read all about coworking!

• Tuesday October 6th, 09:00 Americas Pacific (16:00 UTC) “Writing Tests & Testing in R”, with Yanina Bellini Saibene and co-host Olivier Leroy.
  • Explore how to write tests for R and add some tests to your work or packages
  • Meet co-host, Olivier Leroy, and chat about testing
• Tuesday November 3rd, 09:00 Australia Western (01:00 UTC) “Climate Science in R”, with Steffi LaZerte and co-host Elio Campitelli.
  • Explore how R is used to study the climate
  • Meet co-host, Elio Campitelli, and discuss Climate Science in R
• Tuesday December 8th^*, 14:00 Europe Central (12:00 UTC) “Code Linting in R”, with Steffi LaZerte and co-host Etienne Bacher.
  • Read up on Code Linting and apply some linters to your R code
  • Meet co-host, Etienne Bacher, and discuss code linting in general, or flir and Jarl in particular
    * Note that December coworking is a week later than usual

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

Software

The following two packages recently became a part of our software suite:

• ciecl, developed by Rodolfo Tasso Suazo: Tools for working with the International Classification of Diseases (ICD-10 Chile official MINSAL/DEIS v2018). Includes optimized SQL search with SQLite, fuzzy matching of medical terms (Jaro-Winkler), Charlson and Elixhauser comorbidity calculation, WHO ICD-11 API integration, and hierarchical code validation. Data from Centro FIC Chile DEIS https://deis.minsal.cl/centrofic/. It has been reviewed by Maëlle Salmon and Yanina Bellini.

• brapiR2, developed by Joash Joshua Ayo: Provides pipe-friendly, stateless read access to the Breeding API (BrAPI) v2.1 specification, an open community standard for plant breeding data interchange maintained by the BrAPI project https://brapi.org. Wraps 32 of the 37 BrAPI v2.1 entities across all four modules, Core, Germplasm, Phenotyping, and Genotyping, covering 49 of the specifications 138 retrieval (GET and search) endpoints and returning tidy tibbles ready for analysis. Write and update endpoints are out of scope by design. Features include automatic pagination, async search handling, response caching, parallel batch fetching, and convenience functions for genomic selection workflows (e.g. dosage matrix extraction). Designed for plant breeders and bioinformaticians who need programmatic access to plant breeding databases that implement the BrAPI’ v2 specification. It has been reviewed by David Waring and Jenna Hershberger.

Discover more packages, read more about Software Peer Review.

New versions

The following twenty-two packages have had an update since the last newsletter: visdat (v0.6.1), nycOpenData (v0.2.3), ruODK (v1.6.0), stats19 (v4.1.0), npi (v0.3.1), promoutils (v0.7.0), brapiR2 (v0.2.0), ciecl (v1.0.0), bibtex (v0.5.3), ckanr (v0.9.0), osmapiR (v0.2.6), distionary (v0.2.0), GLMMcosinor (v0.2.2), reviser (v0.3.1), autotest (v0.2), occCite (v0.6.3), readODS (v2.3.6), osmdata (v0.4.1), goodpractice (v1.2.0), galamm (v0.4.1), RAMEN (v2.1.2), and EDIutils (v3.0.1).

Software Peer Review

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

• Two at ‘6/approved’:

  • brapiR2, A Tidyverse-Native Client for the BrAPI v2 (Breeding API) Specification. Submitted by Joash Joshua Ayo.

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

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

  • ibger, Access the IBGE Aggregate Data API from R. Submitted by Andre Leite Wanderley.

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

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

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

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

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

  • camtrapReport, Camera-Trap Report Generator. Submitted by Elham Ebrahimi.

  • nert, Curated Access to TERN Environmental Raster Data. Submitted by Max Moldovan.

  • rfastlowess, High-Performance LOWESS Smoothing for R. Submitted by Amir Valizadeh. (Stats).

  • camtrapReport, Camera-Trap Report Generator. Submitted by Elham Ebrahimi.

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

  • grumpy, Read NumPy .npy and .npz Files. Submitted by Hugo Gruson.

  • tezr, Access Thesis Metadata from Turkiye’s National Thesis Center. Submitted by Emrah Er.

• Five at ‘1/editor-checks’:

  • ocean3d, Three-Dimensional Marine Spatial Analysis. Submitted by Jay Matsushiba.

  • corila, Sparse modelling with grouped and correlated features allowing for privileged information. Submitted by Armin Rauschenberger. (Stats).

  • tarpolyglot, Run Python, Julia, Rust, and C++ Inside targets Pipeline Steps. Submitted by Pierre9344.

  • OptSurvCutR, Optimal Survival Cut-Point Discovery for Time-to-Event Analysis with OptSurvCutR. Submitted by Payton Yau. (Stats).

  • HydraR, Stateful Agentic Orchestration for Scientific Reproducibility. Submitted by Ignatius Pang.

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

On the blog

• Happy Birthday rOpenSci — My Journey from First-Time Developer to Current Opportunities by Yi-Chin Sunny Tseng. Yi-Chin Sunny Tseng shares how the rOpenSci Champions Program helped her grow from a first-time R package developer into an open science contributor, creating tools for biodiversity research and discovering new opportunities along the way.

Calls for contributions

Calls for maintainers

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

• charlatan, create fake data in R. Issue for volunteering.

Calls for contributions

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

Package development corner

Some useful information for R package developers.

{covr2gh}: Coverage summary as GitHub comment

The covr2gh package by Dragoș Moldovan-Grünfeld provides an automated way to summarise the impact of a pull request (PR) on test coverage directly within GitHub. See it in action in a PR in the tfrmt R package.

Deprecation messages for package data

Hugo Gruson wrote an exhaustive post about the deprecation of data in an R package. The post features the delayedAssign() function.

A refactoring story featuring people

Athanasia Mo Mowinckel published “Why ggseg Atlases Became Function Calls”, where she explains how she made data into objects into a package in order to allow re-exporting them. She furthermore tells how she got to that conclusion by trying out different solutions and discussing with other package developers in the rOpenSci Slack workspace.

usethis 3.2.2

The usethis package was updated on CRAN. The changelog features many quality of life improvements around Git and GitHub, and a new use_readme_qmd() function for drafting a README in the Quarto format.

Last words

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

If you haven’t subscribed to our newsletter yet, you can do so though our signup form. Until it’s time for our next newsletter, you can keep in touch with us through our website, Mastodon, or LinkedIn. See you soon!

This is a guest post by Matt Denwood.

After a longer-than-anticipated delay, the macOS installer for JAGS 5.0.0 is now available via SourceForge. There are some big changes for how the official binaries of JAGS 5.x are installed on macOS, which is the topic covered by this news post.

Multi-threading and BLAS linkage

JAGS now has multi-threading support to allow parallel execution of independent chains, which uses Apple’s Grand Central Dispatch for the macOS binaries. This is the default build, as most users will want to benefit from the parallel computation. However, a single-threaded build is also provided as part of the installer, which has the same behaviour as JAGS 4.3.2. A third build is also provided within the installer, which uses single threading along with Netlib’s reference implementation of BLAS and LAPACK rather than Apple’s vecLib (Accelerate framework). By default, all three of these builds are installed so that the user can later switch between them;  however, the installer can be customised to deselect builds that are not needed.  For any users wishing to avoid the larger download and installation footprint (e.g. build/test machines), we have also provided single-build installers via SourceForge. For discussion on the difference between BLAS implementations, see the relevant section of the R for macOS FAQ. It is also possible to build JAGS using OpenMP on macOS – this is not provided as part of the official macOS installer, but the necessary compilation instructions will be added to the macOS section of the JAGS installation manual in due course.

Installation location

1. Although there is no standardised location for installing system-level software on macOS, the /opt directory is now more frequently used (by e.g. homebrew and some elements of R) than /usr/local. Part of the reason for this is that installation under /opt/jags ensures a clean separation between software relating to JAGS and other software that you may have installed under other subdirectories of /opt.
2. The installer now contains multiple builds of JAGS as well as utilities that are specific to macOS (see below) – the versioned directory structure allows the different builds and utilities to coexist.
3. Having a versioned directory structure also allows different major versions of JAGS to be installed on the same system – we hope that this will ease the transition from JAGS 4.3.2 to JAGS 5.0.0 on CRAN’s build systems.

The installer also contains a build of JAGS 4.3.2 that installs under the same /opt/jags directory, with symbolic links to the older /usr/local directory that allow the CRAN build of rjags-4 to continue to function as expected. These symlinks also allow JAGS 5.x to be accessed from the Terminal (and by pkg-config) without updating the PATH variable, although we recommend adding the new directory to the PATH explicitly using (for example): echo 'export PATH="/opt/jags/bin:$PATH"' >> ~/.zshrc

JAGS utilities

As well as JAGS itself, the macOS installer provides the following three utilities:

• jags-version – this allows the active build of JAGS to be switched, both for the command-line version and for rjags.  Switching between major versions of JAGS requires re-installation of the rjags package with matching version, as well as re-compilation of the runjags package. Switching between builds of JAGS with the same major version does not require re-installation of any R package (but may require R to be restarted, if the required packages are loaded).
• jags-uninstall – this allows for complete removal of any JAGS installations detected either under the new /opt/jags directory or the legacy /usr/local directory.
• pkgconf-lite – as macOS does not include pkg-config by default, this small utility (provided by pkgconf) is installed in order to facilitate compilation of rjags on vanilla macOS installations.

See the provided man pages for jags-version and jags-uninstall for more information on these utilities.  Note that these utilities are not installed if using the single-build installers.

Getting help

If you encounter any issues with the macOS installer or utilities, please either post to the macOS-specific JAGS installation thread on the JAGS discussion forum. If you have comments or suggestions regarding the macOS utilities or build process for JAGS then please feel free to open an issue on the GitHub repo for the macOS build scripts.

I had a lot of fun this morning, brainstorming and assembling this Analytics Pipeline for Dashboards with Claude (and yes it takes much, much more than only 10 prompts).

The philosophy is borrowed from Observable Framework’s data loader, but implemented from scratch with voluntarily opinionated choices of libraries for Python, R and Javascript (and no Markdown): pre-compute data with Python + Polars (Polars is my new crush) or R + dplyr, publish it as static files, and explore them with JavaScript in the browser (Arquero for data wrangling, Observable Plot or Highcharts for charts)

The repository is available on GitHub and here’s the quick start guide to run it locally (you need to have Python and R installed on your machine):

make                        # list all targets
uv venv venv                # or: make venv (python -m venv venv without uv)
source venv/bin/activate    # Windows: venv\Scripts\activate
make install                # Polars into venv/, dplyr into .uvr/library/
make dev                    # starts http://localhost:3000 and opens your browser

The repository’s README file is also informative and contains a few more details about the pipeline.

There are ten public training courses left in 2026, running through October and November. If you’ve been meaning to learn Python, build your first Shiny app, or get to grips with machine learning in R, there’s still time to do it before the end of the year.

All courses run online and are delivered live by data scientists and engineers who work on client projects. Every participant gets course notes and scripts, live demonstrations, and hands-on exercises throughout.

What’s On This Autumn

Python in October

October is a good month to pick up Python. Introduction to Python runs on 12th–13th October, Programming with Python follows a week later, and Data Visualisation with Python runs straight after that. You can take all three in turn to go from your first lines of Python to building your own functions and charts in about two weeks.

Shiny

Introduction to Shiny runs on 5th–6th October and takes you from an empty script to a working interactive web app in R. If you already build Shiny apps, Advanced Concepts in Shiny on 23rd–24th November covers how to write maintainable code, build robust apps that handle unexpected user input, and get more out of reactive programming.

Visualisation and Statistics with R

Data Visualisation with ggplot2 runs on 2nd–3rd November. It’s a good fit for anyone who uses R already and wants to produce clear, publication-quality charts. Statistical Modelling with R runs on 16th–17th November and covers hypothesis testing, regression, clustering and principal components analysis.

Machine Learning with Tidymodels

Machine Learning with Tidymodels runs on 4th–5th November, and Advanced Machine Learning with Tidymodels follows on 23rd–24th November. Taking both gives you a complete machine learning workflow in R, from pre-processing and model fitting through to tuning.

Git

Introduction to Git runs on the mornings of 16th–17th November. It covers why version control matters, how to use it properly, and how to start collaborating with others using Git.

Git and Advanced Machine Learning with Tidymodels run in the mornings, so you can take them alongside an afternoon course in the same week.

Autumn 2026 Schedule

Each course runs over two sessions of three and a half hours. Afternoon sessions run from 13:30 to 17:00 and morning sessions from 09:00 to 12:30 (UK time). Click a course to see the full outline and book.

Booking closes one week before each course starts, so Introduction to Shiny closes on 28th September.

Early bird pricing is still available on some November courses: until 4th October for Introduction to Git and Statistical Modelling with R, and until 11th October for both Advanced courses.

Planning Ahead? Book Now for 2027

If the autumn is too busy, you can already book for next year. Our public schedule for January to June 2027 is live, and every course on it is currently at the early bird price.

It follows the same pattern as this year. Introduction to R and Data Wrangling in the Tidyverse run in January and again in April, and the Python courses run in February and May. Machine learning with Tidymodels runs in March.

It also includes Introduction to Bayesian Inference using RStan in January. It runs over four afternoons and covers writing Stan programs for a range of statistical models and checking MCMC diagnostics in R.

Early bird pricing ends about six weeks before each course, so for the January courses that means booking by late November or early December.

We’ll add dates for the second half of 2027 later in the year.

Why Train with Jumping Rivers

Our trainers are practising data scientists and engineers. The examples and exercises in every course come from the work we do with clients, so you learn how these tools are used in real projects.

Every participant receives:

• Comprehensive PDF notes and scripts to keep after the course
• Live demonstrations and hands-on exercises throughout
• Direct access to a trainer for questions

We have delivered over 1,000 courses to organisations including NHS Scotland, Shell, Wessex Water and the Royal Statistical Society.

Additional Perks

We run free webinars throughout the year. Attending a Jumping Rivers webinar gives you:

• Early exposure to new topics in data science and analytics
• Up to 20% off training courses
• Up to 20% off Jumping Rivers conferences

Register for upcoming webinars here.

Training for Teams

We also run in-house training for organisations that want to develop their teams. We tailor each course to your workflows, tools and experience levels, and we can deliver any course in our full catalogue, not just the ones on the public schedule. Group bookings and returning clients get a discount.

To discuss options, email training@jumpingrivers.com.

Book Your Place

View the full schedule and book at jumpingrivers.com/training/public. Places on public courses are limited, so book early if you have a particular date in mind.

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

Before the main content:

• I am creating an R Community on Google Groups. You can join the group using this form.
• With all the funding cuts to PhD studies, I’d appreciate if you can donate here.

In Economics we often use reduced-form models, which are simplified ways to describe a process. I tend to see the criticism to this as unfair, as in physics they do the same, and so many models in physics assume “spherical objects in a vacuum” as a simplification to make real-world problems mathematically solvable.

With reduced form models we are not interested in how faithful those ressemble reality but how useful these are. As an example in support of reduced form models, consider the two extreme cases of models to answer “how to get from 37 Old Queen St to St. James Park underground station?”:

1. Person 1 draws a quick sketch like this on an old receipt

      _ Embassy
     |   
  ---
 |
Underground

1. Person 2 built a very detailed 1:10000 architectural model of St. James neighbourhood and shows you the tiny streets, cars, etc.

For the question, person 1 sketch is sufficient.

The Armington model of trade sits between these two extremes. It considers a world of \(N\) countries, each country \(i\) has a fixed labour endowment \(w_i\) and there is an endogenous labour wage. Income is given by \(Y_i = w_i L_i\). Trade flows are given by \(X_{ij}\) (and so the sales share is \(\gamma_{ij} = X_{ij} / Y_i\)). There is a productivity shift \(\x_i\), a trade cost $_{ij}, and an elasticity of substitution \(\varepsilon\).

The strong assumption in this model is a known and homogenous elasticity of substitution. In this model countries can substitute between Apple and Android phones even when many would argue that those phones compete in different markets.

The model imposes a market clearing condition \(w_i L_i = \sum_{j=1}^N \lambda_{ij} w_j L_j\).

There are many values \(\lambda_{12}, \lambda_{13}, \ldots, \lambda_{1N}\) that balance the market clearing condition but in this model the solution is a gravity-type equation that leads to ${ij} = (i / ({ij} w_i)^{}) / sum{l=1}^N (l / ({lj} w_l)^{}). More about gravity-type equations can be read in Gravity for Undergrads by Prof. Dr. Yoto V. Yotov.

Consider, for example, the EU increase on low value imports from the UK described [https://hboltd.co.uk/uk-eu-shipping-changes-2026-customs-overhaul/]. In this model, instead of a product-by-product specific tariff have an initial UK-EU average trade cost of \(\tau_{ij}\) that increased to \(\tau_{ij}’\). The hat form of this is \(\hat{\tau_{ij}} = \tau_{ij}’ / \tau_{ij}\).

After the trade cost change the market clearing would be \(w_i’ L_i’ = \sum_{j=1}^N \lambda_{ij}’ w_j’ L_j’\), and then

\[ \hat{w_i} \hat{L_i} = \frac{\sum_{j=1}^N \lambda_{ij}’ w_j’ L_j’}{\sum_{j=1}^N \lambda_{ij} w_j L_j} \]

\[ \hat{w_i} \hat{L_i} = \frac{\sum_{j=1}^N X_{ij}’}{w_i L_i} \]

\[ \hat{w_i} \hat{L_i} = \sum_{j=1}^N \gamma_{ij} \hat{X_{ij}} \]

Similarly, $ = .

The new hat equations are useful, for example, to determine the effect of the change to the trade cost over wages \(\hat{w_i}\). The hat market clearing means that \(\hat{w_i} = \sum_{j=1}^N \gamma_{ij} \hat{X_{ij}}\) (labour is exogenous), and therefore \(\hat{w_i} = \sum_{j=1}^N \gamma_{ij} \hat{\lambda_{ij}} \hat{w_j}\).

Merging the hat market clearing with the hat gravity equation, we get

\[ \hat{w_i} = \sum_{j=1}^N \frac{\gamma_{ij} \hat{w_j} / (\hat{\tau_{ij}} \hat{w_i})^{\varepsilon}}{\sum_{l=1}^N \gamma_{lj} \hat{w_l} / (\hat{\tau_{lj}} \hat{w_l})^{\varepsilon}} \]

From the initial equilibrium shares \(\lambda_{ij}\) and \(\gamma_{ij}\) and a fixed \(\varepsilon\) (for example, estimated using USITC data), we have a system of \(N \times N\) equations that we can solve to get \(\hat{w_j}\) and then recover \(\hat{\lambda_{ij}}\). In other words, we simplified an identification problem by solving the differences by resorting on a sufficient statistic instead of estimating all trade costs.

Executive Summary

MicroStrategy (MSTR) has fundamentally transitioned from a traditional enterprise software firm into an equity-based Bitcoin (BTC) holding vehicle and treasury operation. Evaluating MSTR using conventional corporate finance metrics (such as Price-to-Earnings or EBITDA multiples) fails to capture the core driver of its equity valuation: the dynamic Net Asset Value (NAV) premium driven by programmatic Bitcoin treasury expansion.

This technical article explores an end-to-end, reproducible quantitative pipeline built in R. Inspired by recent advancements in automated SEC parsing tools—specifically the secfile architecture highlighted in Interactive Brokers Quant Blog—this pipeline ingests real-time SEC EDGAR filings, dynamically extracts digital asset balance sheet facts, aligns mixed-frequency financial and market data, and models MSTR daily equity log returns using an out-of-sample 82.08% R-Squared structural tidymodels framework.

1. The SEC EDGAR Challenge & The secfile Paradigm

Unstructured Financial Data vs. Machine-Readable XBRL

Historically, extraction of balance sheet metrics directly from SEC filings (Forms 10-K and 10-Q) required brittle HTML regex scraping or costly third-party commercial APIs. Corporate disclosures often vary across reporting periods, particularly when handling emerging asset classes like digital currencies. MicroStrategy’s XBRL taxonomy has evolved, tagging Bitcoin holdings under varying terms such as DigitalAssets, CryptocurrencyHoldings, or within broader Assets line items.

As demonstrated in the Interactive Brokers Quant Blog overview of secfile, direct ingestion of SEC EDGAR financial facts via compliant XBRL parsing solves three primary institutional hurdles:

1. Regulatory Compliance & Transparency: Enforcing strict HTTP User-Agent headers compliant with SEC EDGAR access policies ensures uninterrupted data pipeline ingestion.
2. Dynamic Taxonomy Mapping: By searching the taxonomy dictionary programmatically, the pipeline automatically handles changes in reporting terminology without hardcoded column dependencies.
3. Point-in-Time Alignment: Raw SEC filings record historical submission dates, preventing look-ahead bias when backtesting event-driven trading strategies against historical market prices.

Using tidy evaluation via the rlang package, the pipeline dynamically isolates balance sheet events safely without breaking execution when SEC XBRL tags undergo structural revisions.

2. Mathematical Structure & Structural Economics

The core hypothesis of this quantitative framework is that MSTR daily return dynamics are governed by two orthogonal components:

1. Systematic Asset Return Co-movement: Direct market return exposure to spot Bitcoin returns.
2. Balance Sheet Shock Factors: Discontinuous structural changes in treasury balance sheet holdings resulting from capital raises or debt-funded BTC acquisitions.

Daily Log Returns

Continuous log returns for asset prices are calculated by taking the natural logarithm of the current price divided by the previous day’s price. Taking natural logarithms guarantees additivity over time horizons and prevents non-negative boundary issues inherent to simple percentage returns.

Balance Sheet Shock Factor

To quantify treasury expansion independently of market price fluctuation, the Balance Sheet Shock metric measures the logarithmic growth rate of total Bitcoins held on MicroStrategy’s balance sheet.

Because corporate filings occur at discrete quarterly intervals, daily balance sheet values are carried forward using Last Observation Carried Forward via the zoo package. Consequently, the balance sheet shock value equals zero on non-reporting days, triggering discrete structural impulses only on filings or event reporting dates.

Econometric Specification

The structural linear model specification expresses the daily log return of MSTR as a function of three main components:

• The intercept term, representing base drift.
• The spot Bitcoin log return, weighted by the Equity-to-BTC Beta Elasticity coefficient (measuring leverage and NAV premium response).
• The Balance Sheet Shock factor, weighted by the Treasury Expansion Impact coefficient.
• A residual error term capturing unexplained market noise.

3. Pipeline Architecture & Package Integration

The pipeline follows a modern, production-grade functional programming architecture using the tidyverse ecosystem:

• secfile: High-throughput scraper and XBRL parser for SEC EDGAR datasets.
• tidyquant: Financial wrapper connecting quantmod and PerformanceAnalytics packages into tidy data frames.
• timetk: Provides time-series splitting for chronological out-of-sample data partitioning without temporal leakage.
• tidymodels: Modular modeling framework combining recipe creation, parsnip engine specification, and workflow execution.
• plotly & ggtext: Interactive visualization engine supporting HTML/CSS styled markdown tooltips and dynamic price projection bands.
• zoo: Handles non-homogeneous time-series alignment via Last Observation Carried Forward (na.locf).

4. Chronological Splitting & Model Training

Financial time series violate the Independent and Identically Distributed (i.i.d.) assumption of standard k-fold cross-validation due to autocorrelation. To preserve time order, we employ strict time-based windowing via timetk.

Using recipes and workflows, we define the feature roles and couple them with a native Ordinary Least Squares (OLS) estimation engine. This guarantees clean execution without data leakage across the training split and testing horizon.

5. Empirical Results & Performance Evaluation

Out-of-sample validation was conducted over a 15-day forward horizon, evaluating the structural model against a Naive Baseline Benchmark (1-day lagged return persistence model where predicted return equals yesterday’s actual return).

Performance Metrics Output

• Structural Model RMSE: 0.0241
• Naive Baseline RMSE: 0.0583
• Structural Model R-Squared: 82.08%
• Naive Baseline R-Squared: 4.12%

The Structural Model achieved an out-of-sample R-Squared of 82.08%, significantly outperforming the Naive Baseline. This confirms that equity return variance in MSTR is predominantly explained by spot Bitcoin fluctuations and treasury balance sheet updates rather than simple price momentum.

# ==============================================================================
# SEC XBRL DRIVEN QUANT PIPELINE: MICROSTRATEGY BALANCE SHEET SHOCKS VS BTC
# Author: Selcuk Disci (datageeek.com)
# ==============================================================================

# 1. LOAD REQUIRED LIBRARIES (AUTOMATED PACMAN ENTIRE PIPELINE INGESTION)
# ------------------------------------------------------------------------------
# Check and install pacman package manager if not already available
if (!require("pacman")) install.packages("pacman")
# Ingest core financial, data manipulation, SEC scraping, and modeling libraries
pacman::p_load(secfile, tidyquant, tidyverse, zoo, timetk, tidymodels)

# 2. DEFINE SEC COMPLIANT USER AGENT AND FETCH DATA
# ------------------------------------------------------------------------------
# Set contact email required by SEC EDGAR fair access policy header guidelines
user_agent <- "<your_email_address>"

# Fetch SEC Central Index Key (CIK) identifier for MicroStrategy Inc.
mstr_cik <- get_ciks("MSTR", user_agent = user_agent)

# Retrieve corporate submissions metadata and extract structured XBRL financial facts
mstr_subs <- get_submissions(mstr_cik, user_agent = user_agent)
mstr_facts <- get_data(mstr_subs, user_agent = user_agent)

# 3. EXTRACTION OF BITCOIN HOLDINGS DYNAMICALLY FROM WIDE FORMAT
# ------------------------------------------------------------------------------
# Extract column names from SEC dataset to locate dynamic crypto asset tags
available_columns <- names(mstr_facts)
# Detect target columns matching SEC XBRL taxonomy for digital holdings
target_btc_column <- available_columns[str_detect(available_columns, "DigitalAsset|CryptocurrencyHoldings")]

# Fallback safety net to ensure execution if specific digital asset tags are absent
if(length(target_btc_column) == 0) {
  target_btc_column <- "Assets" 
} else {
  target_btc_column <- target_btc_column[1]
}

# Clean and transform raw Bitcoin holdings time series data
mstr_btc_holdings <- mstr_facts %>%
  select(report_date, !!sym(target_btc_column)) %>% # Dynamically unquote 'target_btc_column' using !!sym() to evaluate the string as a column name
  rename(date = report_date, BTC_Held = !!sym(target_btc_column)) %>%
  mutate(
    date = as.Date(date),
    BTC_Held = as.numeric(BTC_Held) 
  ) %>%
  filter(!is.na(BTC_Held)) %>%
  distinct(date, .keep_all = TRUE) %>%
  arrange(date)

# 4. FETCH AND COMPUTE LOG RETURNS USING TIDYQUANT
# ------------------------------------------------------------------------------
# Define equity and cryptocurrency symbols for market data retrieval
tickers <- c("MSTR", "BTC-USD")

# Download adjusted close daily historical price series and calculate log returns
market_returns <- tq_get(tickers, from = "2020-01-01", get = "stock.prices") %>%
  group_by(symbol) %>%
  tq_mutate(select = adjusted,
            mutate_fun = periodReturn,
            period = "daily",
            type = "log",
            col_rename = "log_return") %>%
  select(date, symbol, log_return, adjusted) %>%
  ungroup()

# Pivot market data into a wide format and standardize asset-specific return/price columns
market_pivoted <- market_returns %>%
  pivot_wider(names_from = symbol, values_from = c(log_return, adjusted)) %>%
  rename(
    MSTR_Log_Return = log_return_MSTR,
    BTC_Log_Return = `log_return_BTC-USD`,
    MSTR_Close = adjusted_MSTR,
    BTC_Close = `adjusted_BTC-USD`
  ) %>%
  mutate(date = as.Date(date))

# 5. ALIGN MIXED FREQUENCY DATA AND COMPUTE SHOCKS
# ------------------------------------------------------------------------------
# Join balance sheet data with market returns and impute missing daily holdings values via forward fill
processed_model_data <- market_pivoted %>%
  left_join(mstr_btc_holdings, by = "date") %>%
  mutate(BTC_Held_Daily = na.locf(BTC_Held, na.rm = FALSE)) %>%
  filter(!is.na(BTC_Held_Daily) & !is.na(MSTR_Log_Return) & !is.na(BTC_Log_Return)) %>%
  mutate(Balance_Sheet_Shock = log(BTC_Held_Daily / lag(BTC_Held_Daily))) %>%
  filter(!is.na(Balance_Sheet_Shock) & is.finite(Balance_Sheet_Shock))

# 6. TIME-BASED DATA SPLITTING VIA TIMETK (STRICT TIME WINDOWS)
# ------------------------------------------------------------------------------
# Partition dataset chronologically to prevent future data leakage during evaluation
data_splits <- time_series_split(
  data       = processed_model_data,
  date_var   = date,
  initial    = "1 year",   
  assess     = "15 days",  
  cumulative = FALSE       
)

# Extract historical training split and testing horizon
train_data <- training(data_splits)
test_data  <- testing(data_splits)

# 7. ESTIMATE LINEAR REGRESSION VIA NATIVE TIDYMODELS COMPONENT
# ------------------------------------------------------------------------------
# Define features profile inside the standard recipe framework
mstr_recipe <- recipe(MSTR_Log_Return ~ date + BTC_Log_Return + Balance_Sheet_Shock, data = train_data) %>%
  update_role(date, new_role = "id")

# Specify parsnip linear regression engine spec
lm_spec <- linear_reg() %>%
  set_engine("lm") %>%
  set_mode("regression")

# Bind graph components into a clean workflow architecture
mstr_workflow <- workflow() %>%
  add_recipe(mstr_recipe) %>%
  add_model(lm_spec)

# Fit model natively on the chronological training split partition
fitted_lm_workflow <- fit(mstr_workflow, data = train_data)

# 8. OUT-OF-SAMPLE PERFORMANCE TESTING WITH NAIVE BENCHMARK USING YARDSTICK
# ------------------------------------------------------------------------------
# Generate clean out-of-sample predictions via standard tidymodels syntax
test_data_predictions <- predict(fitted_lm_workflow, new_data = test_data)

# Construct evaluation frame with actual returns and naive lag benchmark
evaluation_df <- test_data %>%
  select(date, MSTR_Log_Return, MSTR_Close) %>%
  bind_cols(test_data_predictions) %>%
  rename(Predicted_Return = .pred) %>%
  mutate(Naive_Predicted_Return = lag(MSTR_Log_Return, default = first(MSTR_Log_Return)))

# Reshape predictions to comparative long format for unified performance calculation
evaluation_long <- evaluation_df %>%
  select(date, MSTR_Log_Return, Predicted_Return, Naive_Predicted_Return) %>%
  pivot_longer(
    cols = c(Predicted_Return, Naive_Predicted_Return),
    names_to = "model_type",
    values_to = "estimate"
  ) %>%
  rename(truth = MSTR_Log_Return) %>%
  mutate(model_type = if_else(model_type == "Predicted_Return", "Structural_Model", "Naive_Baseline"))

# Compute RMSE and R-Squared accuracy metrics across structural vs naive models
my_financial_metrics   <- metric_set(rmse, rsq)
accuracy_report_tibble <- evaluation_long %>%
  group_by(model_type) %>%
  my_financial_metrics(truth = truth, estimate = estimate) %>%
  ungroup() %>%
  arrange(.metric, model_type)

# Print execution metric table to console
print(accuracy_report_tibble)

# 9. MODERN INTERACTIVE PLOTLY VISUALIZATION (PRICE-BASED DYNAMIC RSI)
# ------------------------------------------------------------------------------
# Load interactive graphics packages required for dynamic reporting
if (!require("pacman")) install.packages("pacman")
pacman::p_load(plotly, scales, glue, ggtext)

# Extract out-of-sample RMSE metric value for confidence band projection
trusted_rmse <- accuracy_report_tibble %>%
  filter(model_type == "Structural_Model" & .metric == "rmse") %>%
  pull(.estimate)

# Extract out-of-sample R-Squared metric value for header annotation
rsq_val <- accuracy_report_tibble %>%
  filter(model_type == "Structural_Model" & .metric == "rsq") %>%
  pull(.estimate)

# Convert predicted log returns back into absolute dollar prices with confidence intervals
df_eval <- evaluation_df %>%
  mutate(MSTR_Yesterday_Close = lag(MSTR_Close, default = first(MSTR_Close))) %>%
  mutate(
    actual  = MSTR_Close,
    pred    = MSTR_Yesterday_Close * exp(Predicted_Return),
    conf_hi = MSTR_Yesterday_Close * exp(Predicted_Return + (2 * trusted_rmse)),
    conf_lo = MSTR_Yesterday_Close * exp(Predicted_Return - (2 * trusted_rmse))
  ) %>%
  filter(date > min(date))

# Build customized hover text data frames for interactive Plotly tooltips
df_plot_actual <- df_eval %>% select(date, actual)  %>% mutate(text_actual = glue("<b>Actual MSTR Price:</b> ${round(actual, 2)}<br><b>Date:</b> {format(date, '%b %d, %Y')}"))
df_plot_pred   <- df_eval %>% select(date, pred)    %>% mutate(text_pred   = glue("<b>Linear AI Pred:</b> ${round(pred, 2)}<br><b>Date:</b> {format(date, '%b %d, %Y')}"))
df_plot_hi     <- df_eval %>% select(date, conf_hi) %>% mutate(text_hi     = glue("<b>Overbought Ceiling:</b> ${round(conf_hi, 2)}<br><b>Date:</b> {format(date, '%b %d, %Y')}"))
df_plot_lo     <- df_eval %>% select(date, conf_lo) %>% mutate(text_lo     = glue("<b>Oversold Floor:</b> ${round(conf_lo, 2)}<br><b>Date:</b> {format(date, '%b %d, %Y')}"))

# Assemble base ggplot layer with confidence bands, actuals, and model projections
p <- ggplot() +
  geom_ribbon(data = df_eval, aes(x = date, ymin = conf_lo, ymax = conf_hi), fill = "#808080", alpha = 0.18) +
  geom_point(data = df_plot_hi, aes(x = date, y = conf_hi, text = text_hi), color = "transparent", alpha = 0, size = 3) +
  geom_point(data = df_plot_lo, aes(x = date, y = conf_lo, text = text_lo), color = "transparent", alpha = 0, size = 3) +
  geom_line(data = df_plot_actual, aes(x = date, y = actual), color = "#2c3e50", linewidth = 1.2) +
  geom_point(data = df_plot_actual, aes(x = date, y = actual, text = text_actual), color = "#2c3e50", size = 2) +
  geom_line(data = df_plot_pred, aes(x = date, y = pred), color = "#e74c3c", linetype = "dashed", linewidth = 1.2) +
  geom_point(data = df_plot_pred, aes(x = date, y = pred, text = text_pred), color = "#e74c3c", size = 2) +
  scale_y_continuous(labels = dollar_format(accuracy = 1)) +
  labs(
    x = "", y = "",
    title = paste0(
      "MicroStrategy (MSTR) <span style = 'color:#2c3e50'>Actual Prices</span> vs ",
      "<span style = 'color:#e74c3c'>Tidymodels Linear AI Forecast</span><br>",
      "<span style='font-size:12px; color:#555555;'>15-Day Trading Horizon | Out-of-Sample R-Squared: ", round(rsq_val * 100, 2), "%</span>"
    )
  ) +
  theme_minimal() +
  theme(plot.title = element_markdown(hjust = 0.5, face = "bold"),
        plot.background = element_rect(fill = "#ffffff", color = NA),
        panel.background = element_rect(fill = "#ffffff", color = NA),
        panel.grid.minor = element_blank())

# Set typography styles for HTML dashboard output
font_family <- list(family = "Roboto Slab, Sans-Serif", size = 16)
label_font  <- list(font = list(family = "Roboto Slab, Sans-Serif", size = 13))

# Convert static ggplot to fully interactive HTML Plotly widget
ggplotly(p, tooltip = "text") %>% 
  style(hoverlabel = label_font) %>% 
  layout(font = font_family) %>% 
  config(displayModeBar = FALSE)

6. Price Transformation & Dynamic Volatility Bands

Log return forecasts are converted back into actionable nominal price predictions using exponential compounding anchored to the previous trading day’s close price. The predicted price equals yesterday’s closing price multiplied by the exponential of the predicted log return.

Empirical Confidence Envelopes (Overbought / Oversold Zones)

Using the empirical Root Mean Squared Error (RMSE) derived from out-of-sample evaluation, dynamic 2-sigma volatility boundaries are constructed around the predicted price path:

• Overbought Ceiling: Calculated as yesterday’s close multiplied by the exponential of the predicted return plus twice the RMSE.
• Oversold Floor: Calculated as yesterday’s close multiplied by the exponential of the predicted return minus twice the RMSE.

Visual Output Analysis

1. Directional Accuracy: The predicted price path (dashed red line) tracks the actual price trajectory (dark solid line) with high fidelity across the 15-day evaluation window.
2. Volatility Envelope Containment: Actual spot equity prices remain fully bounded within the 2-sigma gray confidence envelope, validating the structural model’s error boundary calibration.
3. Price Acceleration Tracking: The sharp upward price repricing between September 14 and September 21 is effectively captured by the structural model due to the immediate integration of spot BTC log returns.

Conclusion

By combining automated SEC EDGAR parsing via secfile package with functional machine learning pipelines in tidymodels, quantitative analysts can build scalable, production-ready asset pricing engines. In the case of MicroStrategy (MSTR), incorporating SEC XBRL balance sheet tracking alongside high-frequency spot crypto returns yields an empirical out-of-sample explanatory power (R-Squared) exceeding 82%.

You can read the original post in its original format on Rtask website by ThinkR here: We chose React to learn it. The AI wrote it for us.

Diary of a rewrite, episode 1 of 4.

In February 2026 we chose React to rebuild an internal tool, knowing full well it would take twice as long as doing it in R. It was a deliberate trade: we were not buying an application, we were buying skills.

In June, halfway through, with two milestones already delivered, we stopped everything and went back to R.

Here are the three reasons. Only one of them is “we were wrong”. The one that settled it is more uncomfortable than that: the learning we were paying a premium for had simply never happened, because the code was arriving faster than we could read it.

StaffPuzzle is our internal staffing tool: who works on what, which half-day. Version 1 is a read-only Shiny application, 1,200 lines, that has been running for years.

Version 2 was meant to add writing to Google Calendar, a simulation mode and historised indicators. We made the stack decision twice, four months apart, and the second one undid the first. Both are written down, with their criteria. That is what makes the story tellable, and it is also, as we will see, what made it possible.

This series covers the four months that followed, in four episodes. It ends on the real numbers: what the rewrite cost, what the hosting costs, and what the licence we avoided was worth.

19 February: choosing React

Three options on the table: all JavaScript (React, Supabase), all R (raw HTMX on raw plumber2), or a hybrid. The criteria table was explicit, and it is worth reading exactly as it was written:

Look at the second row. We knew React would cost twice as much, and we chose it anyway. This is not naivety, it is a deliberate trade. The team is six people, it knows R and has almost no exposure to modern web stacks; version 1 works, there is no urgency. The Architecture Decision Record (ADR) says so plainly: the team’s comfort with React is “junior”, and that is “a deliberate investment, not an asset we already hold”.

In other words, the February decision was not buying an application. It was buying skills, and accepting to pay for them in delay.

Architecture Decision Record: a document that records a significant software architecture choice and the reasoning behind it.

The option that was not in the table

You may have noticed: “stay on Shiny” is not there. The three options examined were all-JavaScript, all-R on HTMX, and a hybrid. None of them meant keeping the tool as it was and adding what it lacked.

This was not an oversight. The question had been settled upstream, in the project brief written a month earlier, under the heading “Technical constraints”:

Target stack: exit Shiny, REST architecture, classic web site.

Written there, at the same level as Google API quotas and OAuth token handling. But Google’s quotas are a constraint: they are imposed on us. “Exit Shiny” is a decision. Filing a decision under constraints takes it out of the debate. You do not argue with a constraint, you work around it.

Nothing forced us off Shiny. Writing to a calendar, simulating, historising indicators: Shiny can do all of that. So the February ADR diligently compared three options inside a perimeter nobody had justified. And in June, when we reopened the file, we went back over every criterion one by one without ever questioning the perimeter.

What matters for what follows. A decision justified by learning is only valid if the learning happens. That is a hypothesis, not an asset, and it can be checked.

3 June: the U-turn

Four months later, two milestones delivered in React, the team reopens the file. Out of it comes a second decision that cancels the first. Three things had moved, and they are not of the same nature.

One: the option we had refused no longer existed

In February, the R option on the table was raw HTMX on raw plumber2. Writing the HTML by hand, setting the attributes by hand. That was a fair criticism, and it had weighed.

In the meantime, two R packages had been released, {htmxr} and {alpiner}, which give that approach the high-level R abstractions it was missing. June’s option B is not February’s option B. That is not the same thing as having been wrong: the world had changed, and the decision had no mechanism for noticing that on its own.

Two: two of the arguments were simply false

February’s table rated rich interactivity as “difficult” on the HTMX side, with drag and drop as the example. In June, someone went and read the React code that had actually been written: there is no drag and drop. It is a lasso selection, which is precisely what {alpiner} does natively. The argument was defending a difficulty that did not exist.

The second one needs more care. The criterion said: React skills are “highly transferable to other web projects”, R skills much less so. That is true, and it is the wrong criterion. It measures what a developer carries away on their CV, not what a company accumulates. For a firm whose business is R, going deeper into R is not low transferability: it is the main asset getting thicker. Maintaining packages headed for CRAN makes us more visible and more capable than a surface knowledge of React ever would have.

Three: the learning had not happened

This is the one that hurts, and it is the one that settled the decision. June’s ADR puts it in a single sentence:

Learning React, the central motivation for February’s decision, is not happening in practice: the developer does not have time to read the code generated by the AI, so the skills objective is not being met.

ADR-003, 3 June 2026

It is worth weighing what that says. The sole justification for four months of extra cost was skills. The code had been written, the application worked, the milestones were delivered. And the benefit that paid for the extra cost had not materialised.

With an AI assisting, a stack decision justified by learning can sabotage itself. The code arrives faster than you understand it, delivery moves forward, and the illusion holds until somebody asks to see what was learned. The subject goes well beyond this article and we will come back to it; here it is enough to note that it changed a technical decision halfway through a project.

What it cost: the price of the U-turn

A reversal is not free, and June’s ADR lists its own downsides before concluding.

The first one is not technical:

• Redoing the communication with our internal users. Changing technical direction mid-project, in front of people waiting for a tool, requires an explanation you would rather not have to give.
• Stabilising still-experimental packages in order to put them in production.
• The bus factor: the ecosystem we chose has one main maintainer, and it is the same person as the project’s developer.
• Salvaging the documentary value of the React code before deleting it: the types and the tests described a data model that was still valid.

The arithmetic, though, was favourable: six to eight half-days to redo everything, against eight to ten to finish in React. A proof of concept written straight afterwards, 280 lines of R reproducing the main screen at 70 milliseconds per request, turned that estimate into a commitment.

The numbers are what convinced us. They are not what decided it.

Three months on: where we are

The application runs in production, writes to the calendars, and historises its indicators. The tests-to-code ratio went from 0.26 to 0.77: 544 tests, no browser, a single language from the calendar all the way to the screen.

That is the gain we had not anticipated. The interface having become a pure function from data to HTML, it is tested like the rest of the code, with no browser and no screenshots. We justified the return by consistency of language; what we got that was most valuable is testability.

And one benefit that was in no table at all: the application became the first real production case for our own packages, the ones in the hyperverse, which no demo example would ever have given us.

What we take away

Writing decisions down is what lets you undo them. That is the main lesson, and it sounds banal right up until you need it. In June, nobody had to reconstruct from memory why React had been chosen: the criteria were written, dated, weighted. It was enough to go back over them one by one and see which no longer held. Without that document, the discussion would have been about people, about who had been right, instead of about criteria.

Separate “the world changed” from “we were wrong”. Of the three reasons, only one is an error of judgement: the drag-and-drop argument, which reading the code would have disproved back in February. The other two are of a different kind: a tool appeared, and a hypothesis about ourselves turned out to be false. Confusing the three produces either pointless guilt or an inability to turn back.

A decision justified by a non-technical benefit must schedule its own verification. February was buying learning. Nobody had planned to check that it was arriving. We found out by accident, four months later, reopening the file for another reason. That is the one thing we would do differently: if a criterion decides, it needs a review date.

Would we do it again. We would probably choose React again, in February 2026, with what we knew then. That is what makes the story interesting: the U-turn is not the correction of a blunder, it is how a decision behaves when you have taken the trouble to write it down.

Has your Shiny application outgrown itself, and you are wondering what to do with it? We read it and give you back a note that says so, including when the answer is “leave it alone”. Write to us in two lines: the size of your application, and how long it has been running.

Diary of a rewrite, four episodes:

1. We chose React to learn it. The AI wrote it for us. (you are here)
2. My first attempt at porting a Shiny app failed, and not for a technical reason (29 September 2026)
3. We left the managed platform. Here is the real bill. (6 October 2026)
4. Five signs your Shiny app has outgrown Shiny (13 October 2026)

The two decisions cited are Architecture Decision Records from the StaffPuzzle repository, dated 19 February and 3 June 2026. The figures are measured on the repository’s branches, not estimated. Disclosure of interest: htmxr and alpiner are written by Arthur Bréant, who led this project and maintains those packages.

This post is better presented on its original ThinkR website here: We chose React to learn it. The AI wrote it for us.

How to join this free online event with Ana Carolina Moreno, Diana Garcia Cortes, Erick Navarro Delgado, Guadalupe Pascal, Valentina Clavijo Mesa, Monika Avila Marquez, Soledad Araya Orrego and Yanina Bellini Saibene.

Más Allá del Código es un encuentro abierto para conocer los proyectos finales de nuestra cohorte: desde la creación, mejora y revisión de paquetes de software hasta las iniciativas de divulgación que los llevaron a la comunidad. Acompañanos a descubrir el impacto técnico detrás de cada herramienta y las oportunidades que abrieron estos proyectos, incluyendo becas, nuevas colaboraciones y el crecimiento profesional de sus creadores.

¡Sumate a la sesión para inspirarte con sus proyectos, hacer tus preguntas en vivo y celebrar el crecimiento de nuestra comunidad de código abierto!

Consulta a continuación las biografías de las personas que exponen y los recursos del evento.

Speakers

Diana Garcia

Bióloga Computacional del Breast Oncology Program en el Dana Farber Cancer Institute, donde investiga las alteraciones del genoma asociadas a resistencia a terapias en cáncer de mama mestastático. Tiene un doctorado en Ciencias Biomédicas y una maestría en Ciencias de la Computación y experiencia desarrollando software tanto en la academia como en la industria. Disfruta mucho enseñar programación, fue profesora en CETYS Universidad, Campus Tijuana, y tallerista en la Escuela de Código Pilares en CDMX. Forma parte de R-Ladies Boston y del programa de campeones de rOpenSci, previamente participó en PyLadies CDMX y en Women Who Code CDMX.

Erick Navarro

Licenciado en biología por la Universidad Nacional Autónoma de México, y candidato a Doctor en Bioinformática por The University of British Columbia. Actual campeón del Programa de rOpenSci.

Ana Carolina Moreno

Periodista especializada en datos e inteligencia artificial. Cuenta con casi 20 años de experiencia en redacciones de Brasil y España, habiendo trabajado en TV Globo, G1, Folha de S. Paulo, La Voz de Galicia, Jornal da Tarde y Terra Magazine. Fundadora del capítulo de RLadies+ São Paulo

Soledad Andrea Araya Orrego

Cientista política especializada en métodos cuantitativos, análisis de datos y ciencia de datos reproducible aplicada a investigación social. Trabajo con R para procesamiento, análisis y visualización de datos, automatización de flujos de trabajo y desarrollo de herramientas abiertas.

Maria Valentina Clavijo Mesa

Estudiante de doctorado en el Politecnico di Milano | Investigo la resiliencia de las infraestructuras críticas expuestas al cambio climático | Cofundadora de la sección de Medellín de RLadies+

Guadalupe Pascal

Soy investigadora en el ámbito de la optimización basada en datos para los procesos de toma de decisiones en empresas y sistemas sociales, todo ello desde una perspectiva regional centrada en el Sur Global. Trabajo desde una perspectiva de género e interseccional basada en los principios de la ciencia abierta y la justicia epistémica.

Monika Avila Marquez

Soy econometrista y estadística y me dedico a la inferencia causal a partir de datos observacionales, con especial atención a los entornos con interferencia, así como al uso de métodos de aprendizaje automático en la econometría de datos de panel. También trabajo en la selección de modelos de efectos aleatorios cruzados para datos experimentales.

Resources

• paquete votosCL que tiene como objetivo facilitar el almacenamiento y manejo de los datos electorales de Chile.
• Issue con la revision por pares del paquete RAMEN
• Regional Association of DNA Methylome variability with the Exposome and geNome (RAMEN)
• Collaborative notes

Join Us!

• Who Everyone is welcome. No RSVP needed, simply connect and/or dial in at the time of the event.
• When Monday, 12 October 2026 08:00 PDT (Monday, 12 October 2026 15:00 UTC)
• Find your timezone
• Add to Calendar.
• How

  Everyone is welcome. No RSVP needed.

  Test your Zoom setup https://zoom.us/test.

  Join the meeting: https://numfocus-org.zoom.us/j/86115160074?pwd=5aGfse7U7xxrABA5lmgdYRRD8hdEtx.1

  ID de la reunión: 86115160074 Código de acceso: 689296

  Find your local number to join by phone: https://zoom.us/u/adAyZGMYrE

TLDR: The last post identified potential transition pathways for truck drivers based on how closely other occupations matched them on skills, abilities and knowledge. This post follows the same basic approach across a wider range of occupations and asks whether two proxies could stand in for the ratings where no O*NET exists: a small language model comparing the meaning of job descriptions and the distance between two jobs’ US Standard Occupational Classification (SOC) codes. Both agree with the skill-based scores on average and appear to complement each other. The SOC code helps separate jobs in different classes of work and the descriptions are better at ranking occupations within similar groups. Together they explain about a third of the variation in the skill-based score overall, and 10 to 15 percent within a major group, which is where a real job seeker is likely to look. Neither is a substitute for the O*NET’s worker characteristic data, but both provide a reasonable sense check of transition pathway estimates which can be useful when it’s not possible to use data from the O*NET.

Background

In the last post I reproduced the similarity scores estimated by Bratanova et al. (2026). They used worker characteristic data from the O*NET to do this, with the basic idea being that two jobs that demand similar things of their workers should on average be easier to transition between. There’s of course more to their paper than this, but this is the crux of the approach.

They’re certainly not the first to use the O*NET data in this way.¹ There are also some great interactive tools² that intuitively illustrate how the approach works. Unfortunately, only found these after I’d pored over the literature as an adviser on a similar project in India. Because the SOC didn’t map cleanly to India’s labour market³, we decided to dropped the O*NET in any case. But, the project left me with three questions that this post tries to answer:

• How much value does detailed worker characteristics data add when compared to similarity measures built from descriptive text, such as job descriptions?
• If the two approaches differ, are the differences significant enough to matter in practice?
• Given occupational classification systems are usually designed to group comparable jobs together, can they serve as proxy for skill gaps?

Occupational Similarity: Three Ways

Had we used the O*NET, answering these questions would have been relatively simple: calculate similarity scores using occupation characteristic from the O*NET and compare these with scores calculated from text-based measures and a numeric indicator of classification similarity, which is more or less what I’m going to do here.

A primer on skill similarity and the O*NET

If you’re unfamiliar with the territory I’d recommend taking a look at my previous post where I explain why skill similarity measures are used for identifying transition pathways and the basic structure of the O*NET data that is drawn on. But, for those unwilling to leave this page:

• The O*NET: rates 800+ occupations across a standard set of characteristics that describe common abilities, skills and knowledge required to perform an occupation. Each characteristic is then given a rating based on the level required and its importance. Importance is meant to measure how critical a characteristic is to a job, whereas the level signifies the level of proficiency required / complexity of the task. For instance, the skill of speaking is important for both lawyers and paralegals, but lawyers are expected to have a higher level of speaking skills compared to paralegals (see here).
• Skill similarity as a proxy for viable transition pathways: When two occupations demand a similar set of abilities, skills and knowledge it’s assumed a worker would find it easier to transition between the two, on average. For instance, because the skills, abilities and knowledge required of truck drivers resembles bus drivers, it might be easier (and require less training) for them to transition between these jobs compared to becoming a data scientist. Both transitions are possible, but one involves overcoming more friction than the other, making the pathway less viable.⁴
• The O*NET and US Standard Occupational Classification (SOC): The SOC uses a tiered classification system to group occupations, with jobs organized into 23 major groups, 98 minor groups and 459 broad occupations. Each SOC code specifies where a job has been assigned, for instance the SOC code for Chief Executives SOC is 11-1011.00, which places it in:

Three ways to test whether two jobs are alike

Because I’m interested in testing the ideas rather than describing them, I’m not going to cover each approach in too much detail, but the TLDR summary of each approach is below:

• Ratings similarity: occupations are compared based on the level and importance of skills, abilities and knowledge required to complete the job.
• Semantic similarity: Occupations are compared to one-another based on the similarity of job descriptions. Measurement is based on semantic similarity of job descriptions, which attempts to judge text based on its meaning rather than contents.
• SOC code distance: Differences in SOC code assignment is used to provide a naive measure of how similar two occupations are likely to be. Because occupations in the SOC are classified based on work performed and, in some cases, on the skills, education and/or training needed to perform the work⁵, it’s expected that similarly classified jobs will on average share similar characteristics (although this won’t necessarily be linear).

Project set-up

Data:

Like the previous O*NET post, analysis uses version 31.0 of the O*NET database. You can download the data used in this post here.

Everything a reader might want to change sits in this one chunk: packages, paths, the rating scale bounds, the embedding model, the labels and the chart theme.

library(tidyverse)
library(janitor)
library(readxl)
library(text)
library(scales)

#where the data lives and where the figures get written
ref_dir_data   <- file.path(".", "Data")
ref_dir_images <- file.path(".", "images")
dir.create(ref_dir_images, showWarnings = FALSE)

#the O*NET rating files. The NAME of each entry becomes the element set label
#in the combined table, which is how essential and transferable skills stay apart.
ref_files_onet <- c(
  abilities           = "Abilities.xlsx",
  skills_essential    = "Essential Skills.xlsx",
  skills_transferable = "Transferable Skills.xlsx",
  knowledge           = "Knowledge.xlsx"
)

#the O*NET file holding one plain-English description per occupation
ref_file_desc <- "Occupation_Data.xlsx"

#O*NET rating scale bounds, used to put both ratings on a 0-100 scale.
#Importance is collected on a 1-5 scale, level on a 0-7 scale.
ref_scale_level_max <- 7
ref_scale_imp_min   <- 1
ref_scale_imp_max   <- 5

#For semantic similarity, all-MiniLM-L6-v2 has been used. 
ref_model_embed <- "sentence-transformers/all-MiniLM-L6-v2"

#readable labels for the five ratings we're comparing 
ref_labels_set <- c(
  abilities           = "Abilities",
  skills_essential    = "Essential skills",
  skills_transferable = "Transferable skills",
  knowledge           = "Knowledge",
  all                 = "Average Across Worker Characteristics"
)

#for chart labels
ref_labels_measure <- c(ref_labels_set, semantic = "Semantic (descriptions)")


#import SOC classification structure
ref_path_soc <- file.path("Data", "soc_structure_2018.xlsx")

# Import data
dta_soc_raw <- read_excel(
  ref_path_soc,
  sheet = "2018 Structure",
  skip = 8,
  col_names = c("major", "minor", "broad", "detailed", "soc_title"),
  col_types = "text"
)

# Data wrangling
# One row per code, with its level
dta_soc_long <- dta_soc_raw |>
  pivot_longer(major:detailed, names_to = "soc_level", values_to = "soc_code", values_drop_na = TRUE)

# One row per detailed occupation with parent codes and titles
lkp_soc_title <- setNames(dta_soc_long$soc_title, dta_soc_long$soc_code)

dta_soc_hier <- dta_soc_raw |>
  fill(major, minor, broad) |>
  filter(!is.na(detailed)) |>
  mutate(
    major_title = lkp_soc_title[major],
    minor_title = lkp_soc_title[minor],
    broad_title = lkp_soc_title[broad]
  ) |>
  select(major, major_title, minor, minor_title, broad, broad_title, detailed, detailed_title = soc_title)


#convert the SOC classification to format for plots
lkp_soc_major <- dta_soc_long |>
  filter(soc_level == "major") |>
  transmute(soc_major       = str_sub(soc_code, 1, 2),
            soc_major_label = str_c(soc_major, " ", str_remove(soc_title, " Occupations$")))

#colours for figures.
ref_col_desc        <- "#D4A017"
ref_col_soc_rank    <- "#1F3A5F"
ref_col_soc_raw     <- "#7A94B8"
ref_col_combo       <- "#2F7E6E"
ref_col_combo_light <- "#86B8AB"
ref_col_primary <- "#1E298D"  # dark blue
ref_col_accent  <- "#26CDDA"  # cyan
ref_col_grid    <- "#E5E7EB"
ref_col_text    <- "#121212"

#a shared theme so every chart in the post looks the same
ref_theme_post <- theme_minimal(base_size = 11) +
  theme(
    panel.grid.minor    = element_blank(),
    panel.grid.major    = element_line(colour = ref_col_grid),
    axis.text           = element_text(colour = ref_col_text),
    axis.title          = element_text(colour = ref_col_primary, face = "bold"),
    strip.text          = element_text(colour = ref_col_text, face = "bold", hjust = 0),
    legend.position     = "bottom",
    plot.title.position = "plot"
  )

#the extra theme settings the three heatmaps share: no gridlines, small
#rotated axis text and a wide legend bar
ref_theme_heatmap <- ref_theme_post +
  theme(
    panel.grid       = element_blank(),
    axis.text.x      = element_text(angle = 90, vjust = 0.5, hjust = 1, size = 7),
    axis.text.y      = element_text(size = 7),
    legend.key.width = unit(1.5, "cm")
  )

Import the data

The code below reads O*NET data adding source labels.

#read each O*NET rating file, tag it with its element set, and stack them
dta_onet_raw <- imap(ref_files_onet, \(ref_file, ref_set) {
  read_excel(file.path(ref_dir_data, ref_file)) |>
    clean_names() |>
    mutate(element_set = ref_set)
}) |>
  list_rbind() |>
  rename(soc_code  = o_net_soc_code,
         soc_label = title)

#one row per occupation: code, title and description
dta_desc_raw <- read_excel(file.path(ref_dir_data, ref_file_desc)) |>
  clean_names() |>
  rename(soc_code  = o_net_soc_code,
         soc_label = title)

#a lookup of code to title, taken from the ratings files
lkp_soc_labels <- dta_onet_raw |>
  distinct(soc_code, soc_label)

Exploratory analysis

Checking coverage

The code below does a simple sense check of the imported data by presenting the number of occupations covered and the number of characteristics measured across O*NET abilities, knowledge and skills. Each characteristic (or element set) should have the same number of values for importance and level.

#characteristic and description counts
sum_onet_coverage <- tibble(
  occupations_rated      = n_distinct(dta_onet_raw$soc_code),
  occupations_described  = n_distinct(dta_desc_raw$soc_code),
  rated_with_description = n_distinct(intersect(dta_onet_raw$soc_code,
                                                dta_desc_raw$soc_code))
)
#
sum_onet_coverage


#characteristic counts 
sum_onet_element_counts <- dta_onet_raw |>
  distinct(element_set, element_name, scale_name) |>
  count(element_set, scale_name, name = "nmb_elements") |>
  pivot_wider(names_from = scale_name, values_from = nmb_elements)

sum_onet_element_counts

Calculating skill-based similarity

We’ll follow a similar approach to the previous post when calculating skill-based similarity scores, with two main differences:

1. We use Manhattan distance to calculate skill gaps via R’s dist() function. This is essentially the same calculation, but more memory efficient.
2. Average skill gaps are calculated for each worker characteristic and then averaged, so abilities, knowledge, essential and transferable skills are weighted equally.

Put level and importance on a common scale

The O*NET rates level from 0 to 7 and importance from 1 to 5, so the code below scales both ratings so they are on a common scale:

dta_onet_scaled <- dta_onet_raw |>
  mutate(scale_name = str_to_lower(scale_name)) |>
  mutate(value_0_100 = case_when(
    scale_name == "level"      ~ 100 * data_value / ref_scale_level_max,
    scale_name == "importance" ~ 100 * (data_value - ref_scale_imp_min) /
                                       (ref_scale_imp_max - ref_scale_imp_min)
  )) |>
  select(soc_code, soc_label, element_set, element_name, scale_name, value_0_100)

Calculating skill gaps

In the last post, each origin occupation was paired with each potential destination job using a join before skill gaps were calculated. This worked, but required creating a large table with a lot of redundant information before making the calculation. The code below uses dist(), which is a more memory efficient way to make the same calculation by:

• Laying out each characteristic and scale across columns, with one job per row.
• Calculating the total skill gaps across all worker characteristics for each job pair.
• Dividing the total skill gap for each job pair by the number of columns to get the average.
• Averaging skill gaps across worker characteristic categories.

#split the rescaled ratings into the four sets, then add a fifth entry holding
#all of them together
dta_onet_by_set <- split(dta_onet_scaled, dta_onet_scaled$element_set)


# Stage 1: widen a set of ratings to one row per occupation, one column per
# element and scale (with columns names based on the element set, name and scale ).
fnc_widen_ratings <- function(dta_ratings) {
  dta_ratings |>
    select(soc_code, element_set, element_name, scale_name, value_0_100) |>
    pivot_wider(names_from  = c(element_set, element_name, scale_name),
                values_from = value_0_100)
}

#create wide dataframe with one row per occupation 
dta_onet_wide <- dta_onet_by_set |>
  map(fnc_widen_ratings)

#to see an example of the input data:
# tmp_dta_onet_wide_example<-dta_onet_wide$skills_transferable |> head()
# rm(tmp_dta_onet_wide_example)

#check that each worker requirement has the expected number of columns and rows
# (894 occupations and 1 column for each element_name and level e.g 52 abilities and two ratings for each (level and importance)  (52*2=104)) 
chk_wide_dimensions <- dta_onet_wide |>
  map(\(dta_wide) tibble(nmb_occupations = nrow(dta_wide),
                         nmb_columns     = ncol(dta_wide) - 1)) |>
  list_rbind(names_to = "element_set")

chk_wide_dimensions

# Stage 2: Calculate the Manhattan distance between every pair of rows, divided by the
# number of columns to make it an average gap. 
fnc_pair_gaps <- function(dta_wide) {
  
  tmp_matrix <- dta_wide |>
    #move the occupation code from a column to the row labels
    column_to_rownames("soc_code") |>
    #dist() needs a matrix of numbers, not a data frame
    as.matrix()
  
  #dist() calculates the total absolute skill gap between occupations across all characteristics and then takes the average.  
  #as.matrix() lays the input data out as a square with one row and one column
  #per occupation, presenting all job pair transiton paths 
  tmp_gaps <- as.matrix(dist(tmp_matrix, method = "manhattan")) / ncol(tmp_matrix)
  
  tmp_gaps |>
    #back to a data frame, with the row labels as a column of origin codes
    as_tibble(rownames = "soc_code_a") |>
    #stack the destination columns into rows: one row per pair
    pivot_longer(-soc_code_a, names_to = "soc_code_b", values_to = "gap_mean")
}

#average gap within each of the four sets
dta_gaps_by_set <- dta_onet_wide |>
  map(fnc_pair_gaps) |>
  list_rbind(names_to = "element_set")

#the combined score: average the four set-level gaps, so each set counts
#equally whether it has 33 characteristics or 52
dta_gaps_all <- dta_gaps_by_set |>
  summarise(gap_mean = mean(gap_mean), .by = c(soc_code_a, soc_code_b)) |>
  mutate(element_set = "all")

dta_gaps <- bind_rows(dta_gaps_by_set, dta_gaps_all)

#every set should produce the same number of pairs: occupations squared
chk_pairs_per_set <- dta_gaps |>
  count(element_set, name = "nmb_pairs")

chk_pairs_per_set

Turn gaps into scores

The code below rescales skill gaps so each job pair gets a similarity score between 0 and 100. Job pairs that require similar worker characteristics (and have smaller skill gaps) receive scores closer to 100, while job pairs with larger skill gaps receive similarity scores closer to zero. Because job pairs include moving to the same occupation, the minimum skill gap is zero.

#| label: similarity-scores

rlt_sim_ratings <- dta_gaps |>
  mutate(
    sim_score = 100 * (1 - (gap_mean - min(gap_mean)) / (max(gap_mean) - min(gap_mean))),
    .by = element_set
  )

#every occupation against itself, which should be 100 in every set
chk_self_pairs <- rlt_sim_ratings |>
  filter(soc_code_a == soc_code_b) |>
  summarise(sim_min = min(sim_score), sim_max = max(sim_score), .by = element_set)

chk_self_pairs

Scores seem to cluster by broad work classes

The heatmap above presents the average similarity score by each SOC major group excluding identical job pairs.

Because self-comparisons have been removed, the downward sloping diagonal line shows the average similarity score for jobs pairs in the same major group. As one might expect, similarity scores tend to be higher for occupation pairs that are both in the same major group. There is also some visual support for similarity scores being lower for larger differences in SOC codes, with lower average similarity scores in the bottom left and upper right third of the heatmap (where the SOC classification differences are highest).

However, neither pattern is unambiguously dominant, with a number of clusters being apparent in the heatmap. For instance, major groups 41 and 42 share unusually high average similarity scores across most other major groups up to 45 to 53. Correspondingly, the average similarity scores appear to be higher for occupation pairs within groups 45 to 53 than outside them. Major groups 11 to 19 exhibit a similar pattern.

Although exploring the source of this is worthy of another post, one explanation for the observed clustering might be that some worker requirements aren’t as widely shared as others, which might point to some inter-group transitions being more difficult than others. This clustering also isn’t all that surprising given the SOC has been designed to group similar worker groups together.

#average score between every pair of major groups, self-pairs excluded
sum_sim_ratings_major <- rlt_sim_ratings |>
  filter(soc_code_a != soc_code_b,
         element_set=="all") |>
  mutate(major_a = str_sub(soc_code_a, 1, 2),
         major_b = str_sub(soc_code_b, 1, 2)) |>
  summarise(value = mean(sim_score), .by = c(element_set, major_a, major_b)) |>
  mutate(panel = factor(ref_labels_set[element_set], levels = ref_labels_set))

plt_sim_ratings_major <- sum_sim_ratings_major |>
  #swap the two-digit codes for the readable major group names
  mutate(major_a = factor(major_a, levels = lkp_soc_major$soc_major,
                          labels = lkp_soc_major$soc_major_label),
         major_b = factor(major_b, levels = lkp_soc_major$soc_major)) |>
  ggplot(aes(x = major_b, y = major_a, fill = value)) +
  geom_tile(colour = "white", linewidth = 0.3) +
  #darker always means more alike
  scale_fill_gradient(low = ref_col_accent, high = ref_col_primary,
                      name = "Average similarity score (0-100)") +
  scale_y_discrete(limits = rev) +
  coord_equal() +
  facet_wrap(vars(panel)) +
  ref_theme_heatmap +
  labs(title ="Average similarity score by SOC major group", 
       x = "SOC major group", y = "SOC Major Group")

plt_sim_ratings_major

Estimate similarity from job descriptions

To estimate similarity scores from job descriptions, the O*NET’s occupation description data is used. For instance, the description for Emergency Management Directors (11-9161.00) is:
Plan and direct disaster response or crisis management activities, provide disaster preparedness training, and prepare emergency plans and procedures for natural (e.g., hurricanes, floods, earthquakes), wartime, or technological (e.g., nuclear power plant emergencies or hazardous materials spills) disasters or hostage situations.

These descriptions aren’t really designed to describe what a job requires in detail and is better thought of as narrative text meant to support the SOC taxonomy. Although this makes them a poor substitute for real job posting data, they do resemble the narrative descriptions used by other occupational classification systems, which is a good worst case scenario to test.

A sentence-embedding model is used to compare how similarly the meaning of two job descriptions are. These are neural network models trained on a large amount of text with the aim of turning a paragraph into a list of numbers designed to encode its underlying meaning. Because it isn’t simply comparing whether each paragraph share similar words, phrases like “drive a truck” and “operate a heavy vehicle” will receive similar scores that indicate they mean something similar. We’ll use cosine similarity, which gives a score of 1 when the meaning of seem to be pointing the same direction, 0 when they’re unrelated and <0 when they have opposite meanings.

I’ve arbitrarily chosen the all-MiniLM-L6-v2 transformer model to do compare job descriptions via the text package. I’m not going to claim a particularly rigorous process was used when selecting the model⁶, but it holds two advantages for this post:

1. It’s small enough to be run locally on a laptop; and
2. The model has been successfully used enough by researchers in adjacent fields⁷.

The code below has the model process the text of each O*NET job description that has a similarity score. Be warned, this can take time to run as it requires having the model process each of the 894 occupations

Note: The first run also downloads the model, which is about 100MB (excluding additional library and Python-related requirements).

textrpp_initialize()

#keep the descriptions of rated occupations only, in code order
dta_desc <- dta_desc_raw |>
  semi_join(lkp_soc_labels, by = join_by(soc_code)) |>
  arrange(soc_code)

#turn each description into a 384-number embedding. The last layer, averaged
#over tokens, is how this model's own authors produce sentence embeddings.
dta_desc_embed <- textEmbed(
  texts  = dta_desc |> select(description),
  model  = ref_model_embed,
  layers = -1,
  aggregation_from_layers_to_tokens = "concatenate",
  aggregation_from_tokens_to_texts  = "mean",
  keep_token_embeddings = FALSE)

#one row per occupation, one column per embedding dimension
dta_desc_vectors <- dta_desc_embed$texts$description

Compare every description with every other

textSimilarityMatrix() computes the cosine similarity between every job description pair to provide a standardized similarity measure across job description pairs. The argument center=TRUE tells the function to subtract values by their corresponding column means. This is the function’s default behaviour and is meant to reduce the chance of scores reflecting paragraphs looking similar as a result of sharing a common format and/or style e.g. all being job descriptions with a comparable style and format.

#cosine similarity between every pair of descriptions, after centring
mat_desc_sim <- as.matrix(
  textSimilarityMatrix(dta_desc_vectors, method = "cosine", center = TRUE)
)
dimnames(mat_desc_sim) <- list(dta_desc$soc_code, dta_desc$soc_code)

#from a matrix to one row per pair, same shape as the ratings table
rlt_sim_semantic <- mat_desc_sim |>
  as_tibble(rownames = "soc_code_a") |>
  pivot_longer(-soc_code_a, names_to = "soc_code_b", values_to = "sim_semantic")

#a description is identical to itself, so the diagonal should be 1
chk_semantic_self <- rlt_sim_semantic |>
  filter(soc_code_a == soc_code_b)

Descriptions show higher average similarity within the same group

The code below produces a heatmap of average semantic similarity scores for each job description pair once identical job pairs are dropped.

Once again, the downward sloping diagonal follows a similar pattern to the skill-based similarity scores: average similarity scores are higher for job pairs within the same SOC major group. The skill-based clustering of groups into wider families of work like the skill-based heatmap isn’t as apparent. Exploring why would make a good excuse for another post, but if I had to guess it’s probably a result of some combination of the following:

• Occupation descriptions from the same group use common terms and formats: Marketing Managers, Sales Managers and Public Relations Managers have descriptions that start with “Plan, direct, or coordinate…”.
• Terms and format don’t appear to carry across groups: a cursory glance at the descriptions used across group clusters in the skill-based score heatmap suggests the format and terms used differ a lot between groups. For instance, occupations in the observed cluster between groups 45 to 53 tend to be described by tools, equipment and outputs that are specific to their work (repairing wind turbines, installing roof support bolts, harvesting vegetables etc), with the result being that similar will often look very different to the model.
• Descriptions vary in length and specificity: the length of descriptions vary from single short sentences to lengthy paragraphs. Unlike the O*NET’s worker characteristic data they also haven’t been standardized to allow comparison across all occupations in the SOC.

Note: It’s also likely that the descriptions of jobs in the same group leaned more on one another than jobs from other groups when they were drafted.

sum_sim_semantic_major <- rlt_sim_semantic |>
  filter(soc_code_a != soc_code_b) |>
  mutate(major_a = str_sub(soc_code_a, 1, 2),
         major_b = str_sub(soc_code_b, 1, 2)) |>
  summarise(value = mean(sim_semantic), .by = c(major_a, major_b))

plt_sim_semantic_major <- sum_sim_semantic_major |>
  mutate(major_a = factor(major_a, levels = lkp_soc_major$soc_major,
                          labels = lkp_soc_major$soc_major_label),
         major_b = factor(major_b, levels = lkp_soc_major$soc_major)) |>
  ggplot(aes(x = major_b, y = major_a, fill = value)) +
  geom_tile(colour = "white", linewidth = 0.3) +
  scale_fill_gradient(low = ref_col_accent, high = ref_col_primary,
                      name = "Average cosine similarity") +
  scale_y_discrete(limits = rev) +
  coord_equal() +
  ref_theme_heatmap +
  labs(x = "SOC major group", y = NULL)


plt_sim_semantic_major

Use the SOC code distance as a “measure” of similarity

I’ve used double quotes around measure to signify I’m making air quotes here, but to be clear: by measure I mean proxy. As even if similar jobs have been placed in the same group and groups are designed to correspond with wider groupings, such as white-collar, service, blue-collar and members of the military, code differences are not designed to precisely quantify where a job sits on some comparable continuum.

To provide an example of what I mean by this: Surgical Assistants (29-9093) and Healthcare Practitioners and Technical Workers (29-9099) are next to each other on the SOC, with codes that are six apart from one another. Yet, Home Health Aides (31-1121) has a code that’s a little over 12 thousand apart from both jobs, despite the three jobs being broadly comparable.

Although this is an extreme example, it does illustrate one of the core problems with relying on raw SOC codes as a proxy: small differences in assignments can result in large code differences that don’t correspond to how similar (or different) two jobs are. Added to this, even when these differences do provide a reasonable proxy for how similar (or different) two jobs are, this is unlikely to be constant across job pairs.

Knowing all this, a reasonable person might ask why look at SOC codes at all. Four reasons:

1. Data availability: Data as rich as the O*NET sometimes isn’t available, can’t be drawn on, or isn’t appropriate to use for the labour market being examined.
2. A bad proxy can still have useful information: Occupation codes might get the ordering of similar jobs right, even when it gets the size of differences wrong.
3. To complement other similarity measures: When data is scarce, classification code differences might provide useful information that other data sources lack, making it useful as part of a wider portfolio of measures, such as part of a composite index.
4. As an intuitive sense-check: When adjacent codes tend to be used for similar jobs, score differences can provide a simple sanity check of other similarity scores.

To limit the effect of adjacent occupations being assigned to different groups, the code below creates a difference measure based on an occupation’s ranking in the SOC, rather than raw code differences. This assumes that adjacent occupations are equally alike across the SOC, regardless of which groups they fall in. This approach was chosen as it feels easier to justify than relying on the raw codes, which sometimes exhibit large differences across classifications that are unlikely to correspond to how similar two jobs are.

Given differences in SOC code rankings are linear, it’s not surprising that the heatmap presents limited variation. However, it does present two patterns that you might expect real similarity scores should follow* when mapped on a classification system like the SOC:

1. Occupations pairs from the same group should exhibit higher similarity with one another; and
2. Similarity scores should decrease as the distance between groups increase.

*(Provided the classification system groups jobs in a way that aligns with the nature of the work.)

Note: The code also estimates raw SOC distance to allow testing whether within-group code differences have predictive power within the same minor group.

#position of each occupation in the SOC list: 1 = first code, 894 = last
lkp_soc_rank <- lkp_soc_labels |>
  arrange(soc_code) |>
  mutate(soc_rank = row_number()) |>
  select(soc_code, soc_rank)

#two readings of the SOC for every pair:
#1.) differences in code rankings
#2.) differences in raw codes
sum_soc_distance <- rlt_sim_ratings |>
  distinct(soc_code_a, soc_code_b) |>
  left_join(lkp_soc_rank, by = join_by(soc_code_a == soc_code)) |>
  left_join(lkp_soc_rank, by = join_by(soc_code_b == soc_code),
            suffix = c("_a", "_b")) |>
  mutate(
    soc_rank_distance = abs(soc_rank_a - soc_rank_b),
    #the full code as one integer, suffix included: "11-9199.11" -> 11919911.
    #Removing the point as well as the hyphen keeps the suffix as whole digits,
    #so two specialisations of one code are 1 apart rather than 0.01, which
    #rounding and floating point would otherwise erase.
    soc_numeric_a     = as.numeric(str_remove_all(soc_code_a, "[-.]")),
    soc_numeric_b     = as.numeric(str_remove_all(soc_code_b, "[-.]")),
    soc_code_distance = abs(soc_numeric_a - soc_numeric_b),
    major_a           = str_sub(soc_code_a, 1, 2),
    major_b           = str_sub(soc_code_b, 1, 2)
  )
#average list-position distance between every pair of major groups, self-pairs
#excluded. The raw code is kept for later; drawn up here it would show the same
#bands with different numbers.
sum_soc_distance_major <- sum_soc_distance |>
  filter(soc_code_a != soc_code_b) |>
  summarise(value = mean(soc_rank_distance), .by = c(major_a, major_b))

plt_soc_distance_major <- sum_soc_distance_major |>
  mutate(major_a = factor(major_a, levels = lkp_soc_major$soc_major,
                          labels = lkp_soc_major$soc_major_label),
         major_b = factor(major_b, levels = lkp_soc_major$soc_major)) |>
  ggplot(aes(x = major_b, y = major_a, fill = value)) +
  geom_tile(colour = "white", linewidth = 0.3) +
  #the ramp is flipped here, so darker still means closer
  scale_fill_gradient(low = ref_col_primary, high = ref_col_accent,
                      name = "Average difference in list position (darker = closer)") +
  scale_y_discrete(limits = rev) +
  coord_equal() +
  ref_theme_heatmap +
  labs(x = "SOC major group", y = NULL,
       title = "Average difference in SOC code ranking by major group")

plt_soc_distance_major

Do the measures agree?

But, do the three measures blend agree with one another?

The code below takes all three measures and combines them into a single dataframe so individual similarity measures can be compared against one another. Groupings are specified to help communicate insights from exploratory analysis (and arguments with Claude) that aren’t shown here for the sake of brevity. In short, the pairs cover:

• Groupings based on where occupation codes differ: Such as job pairs that differ by major, minor, or detailed SOC groupings. This allows us to ask whether the explanatory power of compared measures change when job pairs are drawn from different classifications.
• Groupings where jobs are in the same SOC groups: Such as where compared jobs are in the same major, minor or SOC group. This is to see if the predictive power of measures change when job pairs are drawn from closer occupational groups.

The second group comes from the exploratory analysis left out of the post, which pointed to the SOC distance measure’s explanatory power varying based on how differently two jobs were classified. Because the ranking ignores the size of classification differences within an occupational group, raw SOC code differences might be better at separating jobs in the same group than the ranking, which I wanted to test.

#each pair is labelled by the first level at which its two codes differ. A pair
#that differs at the minor group shares the major group above it; a pair that
#differs only in the suffix is two specialisations of one detailed occupation.
ref_labels_differ <- c(
  "0" = "Differ at major group",
  "1" = "Differ at minor group",
  "2" = "Differ at broad occupation",
  "3" = "Differ at detailed occupation",
  "4" = "Differ in suffix only"
)

#labels for the second cut: the group both codes share
ref_labels_share <- c(
  "All pairs",
  "Same major group",
  "Same minor group",
  "Same broad occupation",
  "Same detailed occupation"
)

#combine ratings in a single dataframe
rlt_pairs <- rlt_sim_ratings |>
  #the combined ratings score only; the four sets had their turn above
  filter(element_set == "all") |>
  select(soc_code_a, soc_code_b, sim_ratings = sim_score) |>
  inner_join(rlt_sim_semantic, by = join_by(soc_code_a, soc_code_b)) |>
  inner_join(sum_soc_distance |>
               select(soc_code_a, soc_code_b, soc_rank_distance, soc_code_distance),
             by = join_by(soc_code_a, soc_code_b)) |>
  #each pair once, lower code first, self-pairs dropped
  filter(soc_code_a < soc_code_b) |>
  mutate(
    #distances flipped to closeness, so every measure reads higher = more alike
    soc_closeness      = -soc_rank_distance,
    soc_code_closeness = -soc_code_distance,
    #"53-3032.00": major "53", minor "53-3", broad "53-303", detailed "53-3032"
    shared_depth = case_when(
      str_sub(soc_code_a, 1, 7) == str_sub(soc_code_b, 1, 7) ~ 4,
      str_sub(soc_code_a, 1, 6) == str_sub(soc_code_b, 1, 6) ~ 3,
      str_sub(soc_code_a, 1, 4) == str_sub(soc_code_b, 1, 4) ~ 2,
      str_sub(soc_code_a, 1, 2) == str_sub(soc_code_b, 1, 2) ~ 1,
      TRUE                                                   ~ 0
    ),
    differ_at = factor(ref_labels_differ[as.character(shared_depth)],
                       levels = ref_labels_differ)
  )

Examining the number in each group

The code below outputs a frequency table for each group. An important characteristic of the groups shown in the table above is that most job pairs span different major groups.

#one row per cut and level, holding the pairs that belong to it
dta_pairs_subsets <- bind_rows(
  tibble(cut_by = "By where the codes first differ", depth = 0:4, subset = ref_labels_differ),
  tibble(cut_by = "By the group the codes share",    depth = 0:4, subset = ref_labels_share)
) |>
  mutate(
    cut_by    = fct_inorder(cut_by),
    subset = factor(subset, levels = c(ref_labels_differ, ref_labels_share)),
    pairs  = map2(cut_by, depth, \(ref_cut, ref_depth) {
      if (ref_cut == "By where the codes first differ") {
        filter(rlt_pairs, shared_depth == ref_depth)
      } else {
        filter(rlt_pairs, shared_depth >= ref_depth)
      }
    }),
    nmb_pairs = map_int(pairs, nrow)
  )

dta_pairs_subsets |>
  select(cut_by, subset, nmb_pairs)

Job description vs. skill-based similarity

The plot above compares similarity scores estimated from job descriptions and the skill-based scores. The relationship is positive, indicating that on average the same job pairs tend to have higher similarity scores using either measure. But, the fit isn’t great, particularly for job pairs with low similarity scores (which account for a large share of job pairs). The glib conclusion is that both measures tend to agree when the similarity of two jobs can be picked up by both measures, but not so much otherwise.

#20,000 random pairs, fixed by the seed so the figure is the same every render.
#The same sample feeds the SOC figure below.
set.seed(20260918)
dta_sample_pairs <- rlt_pairs |>
  slice_sample(n = 20000)

plt_scatter_descriptions <- dta_sample_pairs |>
  ggplot(aes(x = sim_semantic, y = sim_ratings)) +
  geom_point(colour = ref_col_desc, alpha = 0.1, size = 0.6) +
  geom_smooth(method = "loess", se = FALSE, colour = ref_col_text, linewidth = 1.1) +
  ref_theme_post +
  labs(x = "Semantic similarity of descriptions (cosine; more alike to the right)",
       y = "Ratings similarity (0-100)",
       title = "Ratings similarity against semantic similarity of descriptions, 20,000 random pairs")

plt_scatter_descriptions

SOC distance vs. skill-based similarity

The code below uses a scatter plot to compare SOC code distance measures for job pairs with the similarity scores estimated from worker characteristics.

Although I spent quite a bit of time digging into how the two measures correspond with one another (or don’t), I’ve dropped quite a lot of this analysis to keep the post short. Most of the analysis told the same basic story: both distance measures provide some* explanatory power across occupations, but most of their power comes from separating distinct classes of work from one another. The loess curve flattening at greater distances says the rest: once two jobs have SOC list differences of 300 or more, it has little to say about how alike they might be.

*(I was actually surprised how much explanatory power the measure had.)

#ggplot2 apparently can't give two facets different axis transformations, so the two
#panels are drawn separately and placed side by side with patchwork
library(patchwork)

plt_soc_position <- dta_sample_pairs |>
  ggplot(aes(x = soc_rank_distance, y = sim_ratings)) +
  geom_point(colour = ref_col_soc_rank, alpha = 0.1, size = 0.6) +
  geom_smooth(method = "loess", se = FALSE, colour = ref_col_text, linewidth = 1.1) +
  #reversed so closer pairs sit on the right
  scale_x_reverse(labels = label_comma()) +
  ref_theme_post +
  labs(x = "Difference in list position (closer pairs to the right)",
       y = "Ratings similarity (0-100)",
       subtitle = "SOC list position")

plt_soc_raw <- dta_sample_pairs |>
  ggplot(aes(x = soc_code_distance, y = sim_ratings)) +
  geom_point(colour = ref_col_soc_raw, alpha = 0.1, size = 0.6) +
  geom_smooth(method = "loess", se = FALSE, colour = ref_col_text, linewidth = 1.1) +
 
  scale_x_continuous(trans = compose_trans( "reverse"),
                     labels = label_comma()) +
  ref_theme_post +
  labs(x = "Difference in raw code (closer pairs to the right)",
       y = NULL,
       subtitle = "SOC raw code")

plt_scatter_soc <- plt_soc_position + plt_soc_raw +
  plot_annotation(title = "Ratings similarity against the SOC code, (20,000 random pairs)")

plt_scatter_soc

Correlation plot

The code above summarizes the basic message as the scatter plots using a correlation matrix. All measures agree with one another, but the SOC distance measures agrees more closely with the skill-based scores than the descriptions do. This doesn’t necessarily make distance measures a better proxy, but might just point to a similar set of worker characteristics being embodied by both the structure of the SOC and the O*NET’s worker characteristic ratings (more on this below).

#readable names for the four measures, in the order they were introduced
ref_labels_measure <- c(
  sim_ratings        = "Skill-based Ratings",
  sim_semantic       = "Descriptions",
  soc_closeness      = "SOC list position",
  soc_code_closeness = "SOC raw code"
)

#rank correlation between every pair of measures, all pairs
sum_cor <- rlt_pairs |>
  select(all_of(names(ref_labels_measure))) |>
  cor(method = "spearman") |>
  as_tibble(rownames = "measure_a") |>
  pivot_longer(-measure_a, names_to = "measure_b", values_to = "correlation") |>
  mutate(
    measure_a = factor(ref_labels_measure[measure_a], levels = ref_labels_measure),
    measure_b = factor(ref_labels_measure[measure_b], levels = ref_labels_measure)
  ) |>
  #the matrix is symmetric, so show each pair once, below the diagonal
  filter(as.integer(measure_a) > as.integer(measure_b))

plt_cor <- sum_cor |>
  ggplot(aes(x = measure_b, y = measure_a, fill = correlation)) +
  geom_tile(colour = "white", linewidth = 1) +
  geom_text(aes(label = round_half_up(correlation, 2),
                #dark text on light tiles, light text on dark ones
                colour = correlation > 0.5),
            size = 4.5) +
  scale_fill_gradient(low = "white", high = ref_col_soc_rank,
                      limits = c(0, 1), name = "Rank correlation") +
  scale_colour_manual(values = c(`TRUE` = "white", `FALSE` = ref_col_text),
                      guide = "none") +
  scale_y_discrete(limits = rev) +
  coord_equal() +
  ref_theme_post +
  theme(panel.grid = element_blank(),
        axis.text.x = element_text(angle = 30, hjust = 1)) +
  labs(x = NULL, y = NULL,
       title = "Rank correlation between the four measures, all pairs")

plt_cor

Three measures of the same hierarchy(?)

The code below produces a table of average similarity scores by each measure based on job pair grouping differences.

The table points to the SOC hierarchy doing a reasonable job of grouping similar jobs together. It also shows that most job pairs sit in different major groups and score a little above 50. It’s tempting to interpret this as indicating that the scores aren’t particularly good at differentiating jobs. But, this is actually what you’d expect of a database that’s designed to cover everything, rather than provide a representative snapshot of the labour market.

#average of each measure by where the codes first differ
sum_by_level <- rlt_pairs |>
  summarise(
    nmb_pairs         = n(),
    sim_ratings       = mean(sim_ratings),
    sim_semantic      = mean(sim_semantic),
    soc_rank_distance = mean(soc_rank_distance),
    .by = differ_at
  ) |>
  arrange(differ_at) |>
  mutate(sim_ratings       = round_half_up(sim_ratings, 0),
         sim_semantic      = round_half_up(sim_semantic, 2),
         soc_rank_distance = round_half_up(soc_rank_distance, 0))

sum_by_level

But the table also points to an uncomfortable truth about all three measures: it becomes increasingly difficult to use similarity scores to differentiate jobs once they share the same major group. Whether a job pair is in the same minor group, broad occupation, or have differences in the last two digits of the SOC code, the average similarity score shows limited movement. The plot below illustrates this by presenting the entire distribution of each measure by each group. Although the average score rises as the classification narrows, the boxes overlap so much that the classification tells you little about how alike two jobs are once they share a major group. The score still varies inside each grouping, which points to similarity scores being best used as a means to rank alternative pathways, rather than to put a precise number on how similar two jobs are.

Note: the association between similarity measures and SOC groups is not new, it’s a design principle of the system. The SOC classifies occupations based upon work performed, skills, education, training, and credentials⁸ , which are many of the same characteristics the skill based scores were based on. In addition, because managers are intentionally grouped with workers they manage and/or supervise, skill similarity is likely to be deeply nested in the hierarchy. In short, a skills-based score agreeing with a system designed around skills is exactly what you’d expect. It’s also why I felt a SOC code distance measure might provide a plausible proxy for similarity scores in the first place.

dta_plt_level_dist <- rlt_pairs |>
  transmute(differ_at,
            `Skill-based ratings`   = percent_rank(sim_ratings),
            `Job descriptions`      = percent_rank(sim_semantic),
            `SOC list position`     = percent_rank(soc_closeness)) |>
  pivot_longer(-differ_at, names_to = "measure", values_to = "percentile") |>
  mutate(measure = fct_inorder(measure))

plt_level_dist <- dta_plt_level_dist |>
  ggplot(aes(x = differ_at, y = percentile)) +
  geom_boxplot(fill = ref_col_grid, colour = ref_col_text, outlier.alpha = 0.03,
               outlier.size = 0.3, width = 0.6) +
  #the level means, which are the numbers in the table
  stat_summary(fun = mean, geom = "point", colour = ref_col_desc, size = 2.5) +
  scale_x_discrete(labels = \(x) str_remove(x, "^Differ ")) +
  facet_wrap(vars(measure), nrow = 1) +
  ref_theme_post +
  theme(axis.text.x = element_text(angle = 25, hjust = 1)) +
  labs(x = "Where the two codes first differ", y = "Percentile of the measure",
       subtitle = "All three measures rise with the hierarchy on average (gold points), and all three overlap heavily below the major group")

 
plt_level_dist

The explanatory power of descriptions appears to be steady across the hierarchy

To prove that I’m not siding with SOC distance measures as a proxy for similarity paths, the plot below presents insights resulting from a long series of arguments with Claude.

Presented in the plot are rank correlation of each measure with skill-based scores based on how job pairs are grouped within the SOC. Notice that across the largest set of pairs, SOC distance measures agree more closely with the skill-based rankings. But, in every other group job descriptions do a better job of predicting skill-based rankings. Notice also that the raw SOC code doesn’t beat the list measure either, which points to the information lost by using the ranking being minimal, even when jobs are in the same narrow classification.

Aside from this supporting the idea that the distance measure is proxying characteristics embodied by the SOC, it also points to job descriptions providing a more reliable proxy for skill similarity the deeper in the hierarchy the jobs being compared sit. And since that is where most workers are likely to look when they think about changing jobs, descriptions are likely to be the more useful measure for sense checking potential transition pathways.

# Rank correlation of the ratings score with each measure, on one set of
# pairs. Written as a function because it runs once per subset.
fnc_cor_with_ratings <- function(dta) {
  summarise(
    dta,
    `Descriptions`  = cor(sim_ratings, sim_semantic,       method = "spearman"),
    `List position` = cor(sim_ratings, soc_closeness,      method = "spearman"),
    `Raw code`      = cor(sim_ratings, soc_code_closeness, method = "spearman")
  )
}

sum_cor_by_subset <- dta_pairs_subsets |>
  mutate(result = map(pairs, fnc_cor_with_ratings)) |>
  select(cut_by, subset, nmb_pairs, result) |>
  unnest(result)

sum_cor_by_subset |>
  mutate(across(where(is.numeric), \(x) round_half_up(x, 2)))

plt_cor_by_subset <- sum_cor_by_subset |>
  ggplot(aes(y = subset)) +
  geom_vline(xintercept = 0, colour = ref_col_grid) +
  #the gap between the two measures in each row is the story
  geom_segment(aes(x = `List position`, xend = Descriptions, yend = subset),
               colour = ref_col_grid, linewidth = 2) +
  geom_point(aes(x = `Raw code`, colour = "Raw code"),
             shape = 21, fill = "white", size = 3.2, stroke = 1.1) +
  geom_point(aes(x = `List position`, colour = "List position"), size = 4) +
  geom_point(aes(x = Descriptions, colour = "Descriptions"), size = 4) +
  geom_text(aes(x = Descriptions, label = round_half_up(Descriptions, 2)),
            vjust = -1.2, size = 3, colour = ref_col_desc) +
  geom_text(aes(x = `List position`, label = round_half_up(`List position`, 2)),
            vjust = 2.1, size = 3, colour = ref_col_soc_rank) +
  #the pair count, at the left edge of each row
  geom_text(aes(x = -Inf, label = paste0(format(nmb_pairs, big.mark = ",", trim = TRUE),
                                         " pairs")),
            hjust = -0.1, size = 2.8, colour = "grey50") +
  scale_colour_manual(
    values = c(Descriptions = ref_col_desc, `List position` = ref_col_soc_rank,
               `Raw code` = ref_col_soc_raw),
    name   = NULL,
    guide  = guide_legend(override.aes = list(shape = c(16, 16, 21), fill = "white"))
  ) +
  scale_y_discrete(limits = rev) +
  scale_x_continuous(expand = expansion(mult = c(0.25, 0.08))) +
  #one panel per cut_by, each with its own five rows
  facet_wrap(vars(cut_by), scales = "free_y") +
  ref_theme_post +
  theme(panel.grid.major.y = element_blank()) +
  labs(x = "Rank correlation with the ratings score", y = NULL,
       title = "How each measure tracks the ratings, within each subset of pairs")

plt_cor_by_subset

Does the agreement hold across occupation groups?

The plot below presents the rank correlation of each measure with skill-based scores by occupational group, where job pairs are in the same major group.

Once again, description-based scores outperform the SOC distance measure in most groups, and where they don’t the gap is small. The clear exception is Architecture and Engineering, which plausibly stems from the SOC hierarchy separating engineers from the technicians and drafters who support them, a technical grouping that the description-based similarity scores largely miss.

#pairs where both occupations sit in the same major group. Each pair belongs to
#exactly one group, so no need to count both directions.
sum_cor_by_group <- rlt_pairs |>
  filter(shared_depth >= 1) |>
  mutate(major = str_sub(soc_code_a, 1, 2)) |>
  summarise(
    nmb_pairs        = n(),
    `Descriptions`   = cor(sim_ratings, sim_semantic,  method = "spearman"),
    `List position`  = cor(sim_ratings, soc_closeness, method = "spearman"),
    .by = major
  ) |>
  pivot_longer(c(Descriptions, `List position`),
               names_to = "measure", values_to = "correlation") |>
  mutate(major = factor(major, levels = lkp_soc_major$soc_major,
                        labels = lkp_soc_major$soc_major_label))

plt_cor_by_group <- sum_cor_by_group |>
  ggplot(aes(y = major)) +
  geom_vline(xintercept = 0, colour = ref_col_grid) +
  geom_point(aes(x = correlation, colour = measure), size = 2.8) +
  #pair count at the left edge of each row, so the thin groups announce themselves
  geom_text(data = \(d) distinct(d, major, nmb_pairs),
            aes(x = -Inf, label = paste0(format(nmb_pairs, big.mark = ",", trim = TRUE),
                                         " pairs")),
            hjust = -0.1, size = 2.6, colour = "grey50") +
  scale_colour_manual(values = c(Descriptions    = ref_col_desc,
                                 `List position` = ref_col_soc_rank),
                      name = NULL) +
  scale_y_discrete(limits = rev) +
  scale_x_continuous(expand = expansion(mult = c(0.3, 0.05))) +
  ref_theme_post +
  labs(x = "Rank correlation with the skill-based score", y = NULL,
       title = "How well each proxy measure tracks the skill-based score for pairs within a SOC major group")

plt_cor_by_group

Comparing the predictive power of similarity measures

The code below attempts to answer the question that inspired this post in the first place: if you don’t have detailed job characteristic data like the O*NET, how far might you get using national classification structures and job descriptions?

In an attempt to answer this, the explanatory power of each measure is compared across the same five job pair groupings and a group with all job pairs. For each group, a linear regression model is used to estimate the explanatory power of measures individually and when combined together. A model with a higher adjusted R squared score explains a larger share of differences in skill-based score. The gap between the best measure and the combined model can be thought of as the marginal predictive value of adding the other measure.

The results support the same stories that were presented before:

• Neither measure has much explanatory power for skill-based scores: A model that combines both measures at once can explain about a third of variations in skill-based scores across all pairs, but only 10 to 15 percent once two jobs are from the same major group.
• Both proxies appear to work as complements to one another:

Note: p-values and standard errors have not been reported as the observations used in the model aren’t independent e.g. as the same job is included multiple times across job pairs. This isn’t expected to impact the r-squared scores, but will influence anything that relies on the number of independent observations.

#the models, defined once: names, formulas and colours in the same order
ref_models <- tribble(
  ~model,                             ~formula,
  "Descriptions only",                sim_ratings ~ sim_semantic,
  "SOC list position only",           sim_ratings ~ soc_closeness,
  "Descriptions + SOC list position", sim_ratings ~ sim_semantic + soc_closeness
)

ref_col_models <- c(
  `Descriptions only`                = ref_col_desc,
  `SOC list position only`           = ref_col_soc_rank,
  `Descriptions + SOC list position` = ref_col_combo
)

#the adjusted R-squared of each model on one set of pairs. Written as a function
#because it runs once per subset.
fnc_r_squared <- function(dta) {
  ref_models |>
    mutate(r_squared = map_dbl(formula,
                               \(ref_formula) summary(lm(ref_formula, data = dta))$adj.r.squared)) |>
    select(model, r_squared)
}

#all pairs first, as the whole-dataset benchmark, then the pairs cut by where
#their codes first differ
sum_models <- dta_pairs_subsets |>
  filter(subset == "All pairs" | cut_by == "By where the codes first differ") |>
  mutate(result = map(pairs, fnc_r_squared)) |>
  select(subset, nmb_pairs, result) |>
  unnest(result) |>
  mutate(subset = fct_relevel(subset, "All pairs"),
         model  = factor(model, levels = ref_models$model))

sum_models |>
  mutate(r_squared = round_half_up(r_squared, 3)) |>
  select(-nmb_pairs) |>
  pivot_wider(names_from = model, values_from = r_squared)

plt_models <- sum_models |>
  ggplot(aes(y = subset)) +
  geom_vline(xintercept = 0, colour = ref_col_grid) +
  #the gap between the best single measure and the combined model is what the
  #second measure adds
  geom_segment(data = \(d) d |> summarise(x = max(r_squared[model != "Descriptions + SOC list position"]),
                                          xend = r_squared[model == "Descriptions + SOC list position"],
                                          .by = subset),
               aes(x = x, xend = xend, yend = subset), colour = ref_col_grid, linewidth = 2) +
  geom_point(aes(x = r_squared, colour = model), size = 4) +
  geom_text(data = \(d) distinct(d, subset, nmb_pairs),
            aes(x = -Inf, label = paste0(format(nmb_pairs, big.mark = ",", trim = TRUE), " pairs")),
            hjust = -0.1, size = 2.8, colour = "grey50") +
  scale_colour_manual(values = ref_col_models, name = NULL) +
  scale_y_discrete(limits = rev) +
  scale_x_continuous(expand = expansion(mult = c(0.4, 0.05))) +
  ref_theme_post +
  theme(panel.grid.major.y = element_blank()) +
  labs(x = "Share of variation in the skill-based score explained (adjusted R²)", y = NULL,
       title = "How much of the skill-based score each proxy recovers, for all pairs and by where the two codes first differ")

plt_models

What each measure brings

The upshot of the regression results is that neither one of the measures can stand in as a substitute for skill-based similarity scores, but both measures appear to hold value in the absence of occupational data as detailed as the O*NET. A classification distance measure might say little about transition pathways for jobs that resemble each other, but might provide a useful metric for splitting occupations into distinct groups between which transition is less likely. From there, similarity scores based on job descriptions might be useful for ranking potential transition paths between similarly grouped job pairs.

The analysis also points to the measures being practically useful for sense-checking similarity scores, which was what I wanted to test in the first place. For instance, where a country’s occupational classification standards align with the ILO’s International Standard Classification of Occupations, the broad structure of occupation codes should behave in a similar way to the SOC. And classification systems that include text-based information of each job and occupational group provide a means for producing semantic similarity scores. Neither measure is likely to be as rich as the O*NET. But, when the O*NET can’t (or shouldn’t be) used for analysis, both provide accessible metrics for validating transition pathways estimated from non-traditional and unstructured data sources, such as online job postings, vocational curricula and survey data.⁹

Summing up

The post has its origins in a project to identify potential occupational pathways for a country with limited data and a labour force that looked nothing like most of the OECD. Neither problem is unique, but the comparability of local occupations to their US counterparts is a critically important consideration when deciding whether to use the O*NET, particularly when the pathways are intended to inform public policy.¹⁰

One solution to this is to leverage local data sources to develop a local database of occupational characteristics as a substitute to the O*NET. But, even if sufficient data and money exist to make this possible, it can be hard to know whether the identified pathways make sense. This post was meant to test whether description-based similarity measures and classification distance scores might serve as a basic sense check, which the analysis indicates they can. The classification distance measure might point to whether a pathway crosses a major occupational boundary that doesn’t make sense, the descriptions should help check if the rankings of transition pathways within a group look sane. Neither is likely to tell you that a particular pathway is correct, but together they will hopefully point to job transitions that are implausible.

One of the things I argued with Claude about while finalizing this post, was its use of the word “cheap” to describe the two proxy measures. I didn’t like the phrasing (and I still don’t), but Claude is right that both measures are cheap. One is produced by subtracting one classification code from the another. The other comes from a transformer model that runs on my laptop, takes <100MB of space and was deployed across a set of job descriptions that were never meant to be comprehensive outlines of a job. The measures are cheap, which makes it rather extraordinary that they can explain so much of a far richer dataset.

Another point that came to mind while writing this is that by making the skill-based score the thing to predict, I’ve implicitly assumed it’s the standard other measures should be judged against. However, it’s also possible that all three measures carry valuable information about an occupation that isn’t mutually shared. This is testable with the right data, but it’s a point worth keeping in mind as the skill-based score not aligning with either measure might also reflect it lacking important information. If so, the richest and most expensive dataset in the room comes with its own blind spots.

How AI was used to write this post: AI produced the first draft of the code and based on the code used in my last set of analysis of the O*NET. I then proceeded to heavily edit this until it answered the questions a human being (me) might be interested in. The bulk of the writing is my own, with AI only used when I needed inspiration for improving how some points were communicated.

1. Raimi, D. and Greenspon, J., 2025. Finding the Right Fit: What Jobs Offer a Good Match for Fossil Fuel Workers’ Skills? (No. 25-06). Resources for the Future. ︎
2. Resources for the Future, Skills Matching Explorer, https://www.rff.org/publications/data-tools/skills-matching-explorer/ ︎
3. Nor was it necessarily appropriate: Lo Bello, S., Sanchez Puerta, M.L. and Winkler, H., 2019. From Ghana to America: The skill content of jobs and economic development (No. 12259). IZA Discussion Papers. https://papers.ssrn.com/sol3/papers.cfm?abstract_id=3390249 ︎
4. As noted in the previous post, other factors matter too, such as wage differentials, the availability (and proximity) of jobs and the intrinsic benefits of occupations being compared. ︎
5. U.S. Bureau of Labor Statistics (2018), 2018 SOC User Guide: Classification Principles and Coding Guidelines, https://www.bls.gov/soc/2018/soc_2018_class_prin_cod_guide.pdf ︎
6. Claude suggested it and I checked if it made sense ︎
7. See: Saroglou, S., Diamantaras, K., Preta, F., Delianidi, M., Benisis, A. and Meyer, C.J., 2025. Enhancing job matching: occupation, skill and qualification linking with the ESCO and EQF taxonomies. arXiv preprint arXiv:2512.03195. ︎
8. U.S. Bureau of Labor Statistics, Standard Occupational Classification: User Guide, “Classification Principles”. https://www.bls.gov/soc/soc-user-guide.htm. Accessed 19 September 2026. ︎
9. For instance, see: World Economic Forum, 2018. Towards a reskilling revolution: A future of jobs for all. Report, (link); Lassébie, J., Marcolin, L., Vandeweyer, M. and Vignal, B., 2021. Speaking the same language: A machine learning approach to classify skills in Burning Glass Technologies data. OECD Social, Employment and Migration Working Papers, (link); and Granata, J., Posadas, J. and Testaverde, M., 2021. Indonesia’s Online Vacancy Outlook: From Online Job Postings to Labor Market Intelligence 2020. World Bank: Washington, DC, USA. (link). ︎
10. For instance, Lo Bello, S., Sanchez Puerta, M.L. and Winkler find large differences between non-routine and manual tasks when comparing developed and developing countries. See: Lo Bello, S., Sanchez Puerta, M.L. and Winkler, H., 2019. From Ghana to America: The skill content of jobs and economic development (No. 12259). IZA Discussion Papers. https://papers.ssrn.com/sol3/papers.cfm?abstract_id=3390249 ︎

The post O*NET ratings, job descriptions and classification codes as measures of occupational similarity appeared first on Giles.

A Reserving Problem, Solved Without Excel

Ask a reserving actuary how they run a Chain-Ladder and you’ll usually hear “Excel” or the name of a pricey specialized tool. It turns out a modern dataframe library handles it just as well — in a few lines, for hundreds of companies at once.

We use the CAS loss reserving data, specifically the other liability line of business (LoB): 233 US insurers (“GRNAME”), 10 accident years (1998–2007), paid and incurred losses at every development lag (1-10). We treat 2007 as our reporting year, i.e. we simulate a year-end closing.

After reading in the data, we prepare it a bit

import polars as pl

# Other Liability Data Set, December 2025
df = pl.read_csv(
    "https://www.casact.org/sites/default/files/2026-03/othliab_pos_98-07.csv"
)
df_triangle = (
    df
    .group_by("GRNAME", "AccidentYear", "DevelopmentLag")
    .agg(pl.sum("IncurredLosses").alias("Incurred"), pl.sum("CumPaidLoss").alias("Paid"))
    .sort("GRNAME", "AccidentYear", "DevelopmentLag")
)

This is still the “full triangle”, i.e. all 10 development lags for all 10 accident years. Note that df_triangle is not in a triangle format but in a long data format, see Tidy Data by Hadley Wickham.

Chain-Ladder in a Few Lines of Polars

The mechanics are the textbook ones: aggregate claims by accident year i and development lag j (starts at 1, not 0), keep only the upper-left (already observed) triangle, and compute CL development factors

f^{CL}_j = \frac{\sum_{i=1998}^{2008-j} C_{i,j}}{\sum_{i=1998}^{2008-j} C_{i,j-1}}

separately for every company. In polars this is just a group_by plus a couple of window (.over(...)) expressions — no manual loops over triangles required:

def df_2_cl_factors(df, group_by=None):
    g = [] if group_by is None else group_by

    cl_factors = (
        df
        .group_by(g + ["AccidentYear", "DevelopmentLag"])
        .agg(pl.sum("Paid"), pl.sum("Incurred"))
        # Filter upper left triangle
        .filter(REPORTING_YEAR >= pl.col("AccidentYear") + pl.col("DevelopmentLag") - 1)
        .sort(g + ["AccidentYear", "DevelopmentLag"])
        .with_columns(
            PreviousPaid=pl.col("Paid").shift(1).over(g + ["AccidentYear"]),
            PreviousIncurred=pl.col("Incurred").shift(1).over(g + ["AccidentYear"]),
        )
        # Calculate volume-weighted factors per lag
        .group_by(g + ["DevelopmentLag"])
        .agg(
            number_of_years=pl.len(),
            Paid=pl.col("Paid").sum(),
            PreviousPaid=pl.col("PreviousPaid").sum(),
            Incurred=pl.col("Incurred").sum(),
            PreviousIncurred=pl.col("PreviousIncurred").sum(),
        )
        .with_columns(
            f_CL_paid=pl.when(pl.col("PreviousPaid") == 0).then(1).otherwise(pl.col("Paid") / pl.col("PreviousPaid")),
            f_CL_inc=pl.when(pl.col("PreviousIncurred") == 0).then(1).otherwise(pl.col("Incurred") / pl.col("PreviousIncurred")),
        )
        .sort(g + ["DevelopmentLag"])
    )
    return cl_factors

cl_factors = df_2_cl_factors(df_triangle, group_by=["GRNAME"])

We validated the result against R’s ChainLadder package for one insurer, Grinnell Mut Grp, and the incurred factors matched.

Next, we calculate the Chain-Ladder development factors, again for both paid and incurred, this time separately for each company. Note that we take care to only account for the upper left triangle, which is what’s usually available in practice.

def df_2_cl_factors(df, group_by=None):
    g = [] if group_by is None else group_by

    cl_factors = (
        df
        .group_by(g + ["AccidentYear", "DevelopmentLag"])
        .agg(pl.sum("Paid"), pl.sum("Incurred"))
        # Filter upper left triangle
        .filter(REPORTING_YEAR >= pl.col("AccidentYear") + pl.col("DevelopmentLag") - 1)
        .sort(g + ["AccidentYear", "DevelopmentLag"])
        .with_columns(
            PreviousPaid=pl.col("Paid").shift(1).over(g + ["AccidentYear"]),
            PreviousIncurred=pl.col("Incurred").shift(1).over(g + ["AccidentYear"]),
        )
        # Calculate volume-weighted factors per lag
        .group_by(g + ["DevelopmentLag"])
        .agg(
            number_of_years=pl.len(),
            Paid=pl.col("Paid").sum(),
            PreviousPaid=pl.col("PreviousPaid").sum(),
            Incurred=pl.col("Incurred").sum(),
            PreviousIncurred=pl.col("PreviousIncurred").sum(),
        )
        .with_columns(
            f_CL_paid=pl.when(pl.col("PreviousPaid") == 0).then(1).otherwise(pl.col("Paid") / pl.col("PreviousPaid")),
            f_CL_inc=pl.when(pl.col("PreviousIncurred") == 0).then(1).otherwise(pl.col("Incurred") / pl.col("PreviousIncurred")),
        )
        .sort(g + ["DevelopmentLag"])
    )
    return cl_factors

cl_factors = df_2_cl_factors(df_triangle, group_by=["GRNAME"])
cl_factors.filter(pl.col("GRNAME") == "Grinnell Mut Grp")
┌────────────┬────────────┬────────────┬────────┬───┬──────────┬────────────┬───────────┬──────────┐
│ GRNAME     ┆ Developmen ┆ number_of_ ┆ Paid   ┆ … ┆ Incurred ┆ PreviousIn ┆ f_CL_paid ┆ f_CL_inc │
│ ---        ┆ tLag       ┆ years      ┆ ---    ┆   ┆ ---      ┆ curred     ┆ ---       ┆ ---      │
│ str        ┆ ---        ┆ ---        ┆ i64    ┆   ┆ i64      ┆ ---        ┆ f64       ┆ f64      │
│            ┆ i64        ┆ u32        ┆        ┆   ┆          ┆ i64        ┆           ┆          │
╞════════════╪════════════╪════════════╪════════╪═══╪══════════╪════════════╪═══════════╪══════════╡
│ Grinnell   ┆ 1          ┆ 10         ┆ 59563  ┆ … ┆ 191044   ┆ 0          ┆ 1.0       ┆ 1.0      │
│ Mut Grp    ┆            ┆            ┆        ┆   ┆          ┆            ┆           ┆          │
│ Grinnell   ┆ 2          ┆ 9          ┆ 89554  ┆ … ┆ 172605   ┆ 166116     ┆ 1.727141  ┆ 1.039063 │
│ Mut Grp    ┆            ┆            ┆        ┆   ┆          ┆            ┆           ┆          │
│ Grinnell   ┆ 3          ┆ 8          ┆ 107323 ┆ … ┆ 151543   ┆ 151052     ┆ 1.38351   ┆ 1.003251 │
│ Mut Grp    ┆            ┆            ┆        ┆   ┆          ┆            ┆           ┆          │
│ Grinnell   ┆ 4          ┆ 7          ┆ 106828 ┆ … ┆ 129664   ┆ 131101     ┆ 1.135997  ┆ 0.989039 │
│ Mut Grp    ┆            ┆            ┆        ┆   ┆          ┆            ┆           ┆          │
│ Grinnell   ┆ 5          ┆ 6          ┆ 96179  ┆ … ┆ 105967   ┆ 106851     ┆ 1.088687  ┆ 0.991727 │
│ Mut Grp    ┆            ┆            ┆        ┆   ┆          ┆            ┆           ┆          │
│ Grinnell   ┆ 6          ┆ 5          ┆ 82541  ┆ … ┆ 86134    ┆ 86800      ┆ 1.05192   ┆ 0.992327 │
│ Mut Grp    ┆            ┆            ┆        ┆   ┆          ┆            ┆           ┆          │
│ Grinnell   ┆ 7          ┆ 4          ┆ 64703  ┆ … ┆ 66151    ┆ 66159      ┆ 1.019475  ┆ 0.999879 │
│ Mut Grp    ┆            ┆            ┆        ┆   ┆          ┆            ┆           ┆          │
│ Grinnell   ┆ 8          ┆ 3          ┆ 48924  ┆ … ┆ 49409    ┆ 49457      ┆ 1.008015  ┆ 0.999029 │
│ Mut Grp    ┆            ┆            ┆        ┆   ┆          ┆            ┆           ┆          │
│ Grinnell   ┆ 9          ┆ 2          ┆ 32840  ┆ … ┆ 33029    ┆ 32969      ┆ 1.006343  ┆ 1.00182  │
│ Mut Grp    ┆            ┆            ┆        ┆   ┆          ┆            ┆           ┆          │
│ Grinnell   ┆ 10         ┆ 1          ┆ 15785  ┆ … ┆ 15915    ┆ 15908      ┆ 1.002413  ┆ 1.00044  │
│ Mut Grp    ┆            ┆            ┆        ┆   ┆          ┆            ┆           ┆          │
└────────────┴────────────┴────────────┴────────┴───┴──────────┴────────────┴───────────┴──────────┘

We validated the result against R’s ChainLadder package for one insurer, Grinnell Mut Grp, and the incurred factors matched, see the linked notebook.

A Pattern That Doesn’t Play Nice

Zooming into the CL factors of individual companies, two stood out: the paid pattern of Virginia Mut Ins Co and the incurred pattern of Grinnell Mut Grp. Neither is monotonic — a red flag, since most commercial reserving tools only offer parametric curve fits that are strictly monotonic (to be fair, parametric curve fits are more for tail factor estimation).

Grinnell’s incurred CL factors rise above 1 early on, dip below 1 (case reserves getting released, maybe subrogation), then climb back above 1 with some zig-zag in later years. A parametric curve simply can’t represent that saddle-shaped pattern.

Enter Whittaker-Henderson

Whittaker-Henderson (WH) smoothing is a non-parametric, penalized smoother — no functional form assumed, just a trade-off between fitting the data and penalizing roughness. That trade-off is controlled by two knobs: the penalty order (we use the standard order=2, i.e. penalizing curvature) and the penalty strength lamb, which we set by eye.

It happens to be a perfect match here: development factors are a discrete-time signal with equal time steps, exactly what WH smoothing was designed for (it dates back to Georg Bohlmann in 1899 — arguably it should be called Bohlmann-Whittaker-Henderson). As of scipy 1.18, it ships out of the box as scipy.signal.whittaker_henderson — full disclosure, I contributed that implementation, so I might be a little biased towards finding excuses to use it

from scipy.signal import whittaker_henderson

f_smooth = whittaker_henderson(
  signal=f_cl_inc, weights=previous_incurred, lamb=1e4
).x

One neat property of WH smoothing: it preserves the weighted in-sample sum of the signal, so the smoothed factors reproduce the observed incurred losses exactly on the fitted range. Out-of-sample — i.e. for the not-yet-observed lower-right triangle — the smoothed and raw factors diverge, which is exactly where it matters for reserving.

Does It Change the Reserve?

Yes, a bit:

The total reserve shifts by about 3% (+31.1) — modest for the total, though individual accident years move more (some by double digits in percentage terms), since smoothing lets a factor’s neighbors pull it away from its own noisy ratio.

If we compare against the incurred losses after all 10 years of development — the closest proxy we have to the true ultimate loss — both CL and the WH-smoothed CL turn out to overestimate it. The un-smoothed CL is only marginally closer.

So for this one company, smoothing made the development pattern easier to reason about, but it didn’t make the forecast more accurate — the difference (about 31) is small compared to the 4,582 standard error Mack’s method reports for this triangle, so neither method is clearly better here. A good reminder to check against a holdout whenever one is available.

Takeaways

• polars makes Chain-Ladder wrangling compact and fast, even across hundreds of companies at once.
• Volume-weighted CL factors computed in polars match R’s ChainLadder package exactly — always reassuring when switching tools.
• Whittaker-Henderson smoothing is a flexible alternative to parametric curve fitting whenever development factors are noisy or non-monotonic, and it’s now built into scipy. A smoother pattern isn’t automatically a more accurate one, though — always check against a holdout when you can.

Natural next steps could be to add a tail factor, estimate reserve uncertainty (à la Mack’s method), or let REML pick lamb automatically instead of choosing it by eye.

The full notebook — all the polars code, the charts, and the R comparison — is on GitHub. This post as well as the notebook was AI reviewed.

Spotted a bug, or have a favorite way to smooth development factors? Let me know in the comments!

This post is a follow-up to my previous posts on semi-parametric option pricing. A link to the study (accepted for presentation at the osQF 2026 conference) is provided at the end of this post.

In this study (R code provided), we build an empirical pricing measure for options directly from their underlying’s historical dynamics, requiring no market of option prices to calibrate against. Starting from a single observed discounted price series, we filter out linear predictability via an AR(1) fit and preserve the remaining dependence structure through a stationary block bootstrap of the residuals. The resulting empirical distribution is then given a minimal adjustment consistent with no-arbitrage: a single scalar correction, following Duan and Simonato [1998], that enforces the martingale condition required by the Fundamental Theorem of Asset Pricing (FTAP). We are explicit that this construction is not a recovery of “the” risk-neutral measure itself in the usual economic sense, but rather a pricing measure obtained by disturbing the historical dynamics as little as possible, to verify the FTAP. We validate the methodology empirically against the implied volatility surface of DAX index European options quoted July 5, 2002, comparing prices directly rather than implied volatilities. Because the construction requires no option-market input at any stage, it extends naturally to path-dependent payoffs; we illustrate this on arithmetic Asian options. 1 Motivation Standard option pricing practice requires the assumption of a parametric family for the option underlying’s dynamics (geometric Brownian motion [Black and Scholes, 1973], stochastic volatility [Heston, 1993], jump-diffusion [Merton, 1976]), and calibrating the parameters of that parametric family to observed option prices. This works well when a liquid market of option prices exists to calibrate against. It does not help when no such market exists, for examples for path-dependent or structured payoffs written on an underlying with no quoted option market (an Asian option on a private index, an embedded option inside an insurance product, a participation certificate). This note develops an alternative construction for option pricing that requires no market prices at all. It builds an empirical pricing measure directly from the historical dynamics of the underlying asset, and imposes no parametric distributional assumption beyond what the data itself exhibits. The construction is then made consistent with the Fundamental Theorem of Asset Pricing (FTAP). We first state the FTAP theorem, then show how it motivates the step of our pricing estimator, before turning to empirical validation.

https://www.researchgate.net/publication/414513095_Semi-parametric_option_pricing_based_on_underlying’s_historical_data

I’m now interested in constructive remarks and feedback on the study, that will allow to enrich, improve and robustify the methodology.

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

With Obama’s recent comments on AI, I was thinking about what once was the Open Government initiative that, while remarkable in many senses, used to confuse data and information. That was back in 2009, and a large fraction of my work now and then consisted in obtaining data, this is values might pull out of a database such as region, city and variables of interest provided in formats meant to provide interpretations and drawing conclusions. In other words, a part of my work has the complication of dealing with data provided in the wrong format for such goal, as it is the case of releasing data in PDF format.

When an economist like me goes to a government website such as the USITC, it is often the case that he/she is looking for raw data so that it can be analyzed and shaped in different ways, usually detecting inconsistencies, and interpreted into the form ot summary tables and plots inside a document. USITC uses the correct formats for that goal but other public offices tend to insist on the wrong formats even now that AI has become a salient issue and it is often perceived as a multi-purpose tool.

Sharing data in PDF format will continue to be a common practise, and often once that overly complicates scientific work and investigative journalism as it was the case of Dollars for Docs. The popular saying is that if you have a hammer, then everything looks like a nail. AI-based workflow may increase the speed at which we can create a software prototype and still work within an organizational culture that confuses data and information.

With the rise of AI, government agencies can streamline the workflows that result in a PDF for effective information distribution. The information in a PDF document can be read in its electronic form, printed and reshared easily. However, adding AI to government workflows can add a vicious instead of a virtuous element. For the economist who wants raw data, PDFs made faster as in more frequent updates is the incorrect choice. The risk is that AI can become a new hammer, a tool to try to hit everything with in the same way as we continue to use PDFs without thinking of the suitability behind the technological choice.

I am quite sure that many data issues, not just in government but also in businesses, do not require AI at all. Around five years ago, a much less salient hype was about Machine Learning (ML), and many comments gave me the impression that ML was becoming a new hammer. In particular, many data issues such as bottlenecks and organization can be solved by moving data from multiple files into a proper Structured Query Language (SQL) database engine such as PostgreSQL or MariaDB, both with a proven track record.

Raw data also requires metadata describing the time range, source, and descriptions for the individual variables, and multiple important properties. A PDF seems to be useful to provide a data dictionary, and so could be a Word document or even a plain TXT files, which had led me to think about the lack of debate on the suitability of AI and many other tools. For the SQL case, it is justified for large raw datasets where Excel spreadsheets fall short, and there will be cases where Excel spreadsheet will be just right.

There are many other raw data formats that might be more suited to particular needs such as plain spreadsheet files (CSV). There are specific needs, as it is the case of geographical data, for which we have shapefiles (SHP) for holding spatial feature information with their respective database version known as PostGIS. Like different data format serve different goals, it is up to us to see the AI hammer in the toolbox and know when to use it.

How is it relevant to the R community? I have received an increasing number of low effort Pull Requests (PRs) to the different R packages that I maintain and I see the IA hammer in action. Most of those PRs failed to answer simple questions from my side such as “is this data free to be re-distributed?”. I usually get a positive feeling when I receive an email or see an open issue on GitHub from a user asking me about a particular dataset. A few times, fortunately a few ones, I have received a PR with a simple change because I documented a function with a bad wording or I had typos, and so I am thankful for those corrections that sometimes came without much debate.

I think that I will continue to receive AI-based PRs, which it not bad by itself, and so it is not using it to add more examples, add checks or rewrite parts of a code. My same critique will apply if, instead of AI-based changes I would be receiving changes consisting in copying and pasting without much consideration about the relevance of doing it.

Today’s just a very short sequel to last week’s post, where I had a look at some very skewed distributions to test the idea that sample sizes sometimes need to be in the tens of thousands for the sample mean to have a normal distribution. Turns out they do.

I had a bit of unfinished business at the back of my mind, which was “would a bootstrap confidence interval do any better?”. Hence today’s new set of simulations.

I compared the coverage of a 95% confidence interval for the mean constructed the traditional way—like they teach it in basic stats courses—from a few heavily skewed distributions. I also constructed a 95% confidence interval using the bias-corrected and adjusted bootstrap method, which I believe is the best candidate to work in a wide variety of bias and skew situations.

To skip to the chase, here’s the results. Turns out that a) the bootstrap does indeed do considerably better than just relying on the central limit theorem, particularly with smaller sample sizes; and b) it’s still got coverage a lot less than the 95% we wanted:

No surprise here; from what I understand of the history, this is pretty much exactly what the BCa bootstrap was developed for. So we’re on it’s home ground, and it does (relatively) well. But those actual coverage numbers are still well below 95%, for both methods.

Here’s the code that did that. It’s very similar to the one from a few days ago.

library(tidyverse)
library(actuar)
library(glue)
library(scales)
library(boot)

# set the below to TRUE if running for the first time
run_sims <- FALSE

set.seed(123)

# Number of repeats for each combination of sample size and population:
today_reps <- 1000

# Population size:
N <- 1e6

# Modified version of the function we used last week, this time just looking at
# coverage of confidence intervals and using a BCa bootstrap to compare to the
# traditional asuymptotic CLT/normal assumed one:
sim_clt2 <- function(
  x,
  n = 30,
  reps = today_reps,
  replace = TRUE,
  conf = 0.95,
  boot_R = 3001,
  ...
) {
  true_mean <- mean(x)

  samples <- replicate(
    reps,
    sample(x = x, size = n, replace = replace),
    simplify = FALSE
  )
  means <- sapply(samples, mean)

  covered_clt <- sapply(samples, function(s) {
    # rely on asymptotic normality to estimate a confidence interval and check
    # for coverage for each sample
    se <- stats::sd(s) / sqrt(length(s))
    ci <- mean(s) + c(-1, 1) * qnorm((1 - conf) / 2 + conf) * se
    ci[1] <= true_mean & true_mean <= ci[2]
  })

  covered_boot <- sapply(samples, function(s) {
    b <- boot::boot(
      data = s,
      statistic = function(x, w) {
        mean(x[w])
      },
      R = boot_R
    )
    ci_boot_res <- boot::boot.ci(b, conf = conf, type = "bca")
    ci_boot <- ci_boot_res$bca[4:5]
    ci_boot[1] <= true_mean & true_mean <= ci_boot[2]
  })

  return(list(
    coverage_clt = mean(covered_clt),
    coverage_boot = mean(covered_boot)
  ))
}

# Populations we're going to use
pops <- list(
  exp(rnorm(N)),
  exp(rnorm(N, sd = 2)),
  exp(rexp(N, rate = 2))
)

# Sample sizes we're going to use
ns <- c(10, 30, 200, 1000)

# Run simulations:
if (run_sims) {
  results <- expand_grid(pop = 1:3, n = ns) |>
    mutate(coverage_clt = NA, coverage_boot = NA)

  # this - obviously when you think about what it's doing - will take a long time
  # (~2 hours) to run. It's embarassingly parallel so could consider parallelising
  # it easily enough, but there is a lot of demands on memory so for my laptop is
  # probably not going to be worth trying this as the machine wouldn't be able to
  # do multiple goes of the 3000 rep bootstrap, 1000 rep simulation from a 1e6
  # population at once.
  for (i in 1:nrow(results)) {
    cat(i)
    param <- results[i, ]
    tmp <- sim_clt2(pops[[param$pop]], n = param$n)
    results[i, ]$coverage_clt <- tmp$coverage_clt
    results[i, ]$coverage_boot <- tmp$coverage_boot
  }

  save(results, file = glue("0333-boot-results-{Sys.Date()}.rda"))
} else {
  lf <- sort(
    list.files(pattern = "0333-boot-results.*\\.rda$"),
    decreasing = TRUE
  )
  load(lf[1])
}

# labels for the populations:
pop_labs <- c("log normal(0,1)", "log normal(0,2)", "exponential(2)")

# Draw plot:
p <- results |>
  mutate(lab = pop_labs[pop]) |>
  mutate(lab = fct_reorder(lab, coverage_boot)) |>
  ggplot(aes(x = coverage_clt, y = coverage_boot, colour = lab)) +
  geom_abline(slope = 1, intercept = 0, colour = "grey50") +
  geom_point(size = 2) +
  geom_text_repel(aes(label = comma(n)), seed = 123, alpha = 0.5) +
  coord_equal() +
  scale_x_continuous(label = percent) +
  scale_y_continuous(label = percent) +
  labs(
    x = "Confidence interval from asymptotic normality includes the mean",
    y = "Confidence interval from BCa bootstrap includes the mean",
    title = "Bootstrap outperforms asymptotic normality assumption with smaller n.",
    subtitle = "Proportion of time the 95% confidence interval actually contains the true value.
Labelled numbers indicate sample sizes. Diagonal line shows equal performance.",
    colour = "Population distribution:"
  )

 print(p)

I still haven’t looked at the point—raised by Professor Harrell himself after my last post—of the assymetry of these confidence intervals, which causes a whole new set of problems. I think I’ve run out of oomph for looking at that, but it is actually an important point to remember. Maybe some time later.

That’s it for today really. I still think the bootstrap is a close to magic as you get in frequentist statistics, and I thoroughly recommend it. It’s good stuff. But when you’ve got a sample size of 10, 30, 200—sometimes even when you’ve got 1,000, 10,000 or 50,000—there’s just limits to what you can do.

You can read the original post in its original format on Rtask website by ThinkR here: From URL to Theme: How {css2r} Steals a Website’s Colours, and Where It Gives Up

You have thirty minutes before the demo, and your Shiny app still looks like Bootstrap 5 out of the box.

We have all been there. The analysis is solid, the model is good, the tables are clean. And then someone from the communication team walks past your screen and says, very politely, that “it doesn’t really look like us”. They are right. Your app is blue-ish, your company is orange. Somewhere there is a 60-page brand book you have never opened, and you have no intention of reading it before 2pm.

So here is the naive question we asked ourselves at ThinkR: the company website already carries the brand. The colours are in its CSS, right there, publicly served over HTTP. Why not just go and take them?

That is {css2r}, and the small Shiny app that wraps it, Shiny Copy. You give it a URL, it gives you a bslib::bs_theme() call. It is live at connect.thinkr.fr/css2r. Open it in another tab, paste your company’s address, and come back. This article will still be here.

It works. And the interesting part of this article is not that it works. It is everything the web refuses to give you along the way. Buckle up, we’re talking about CSS.

Thirty seconds: a URL in, a theme out

remotes::install_github("ThinkR-open/css2r")
library(css2r)
thinkr <- css2r$new(url = "https://thinkr.fr")
#> ✔ Internet ok
#> ✔ html page downloaded
#> ✔ CSS links extracted
#> ✔ CSS links filtered
#> ✔ CSS downloaded
#> ✔ Colors extracted successfully
#> ✔ Colors analyzed successfully.
#> ℹ No Google Fonts detected.
#> ✔ Shiny theme code generated.

Nine steps, a few seconds, and the object holds everything it found:

thinkr$top_colors
#> $white_black
#>      Color Count
#> 32 #FFFFFF     1
#>
#> $top_colors
#>     Color Count
#> 1 #38404C   133
#> 2 #0046C8    72
#> 3 #F05622    40
#> 4 #20B8D6    23
cat(thinkr$shiny_code)
#> fluidPage(
#>   theme = bslib::bs_theme(
#>     bg = "#FFFFFF",
#>     fg = "#38404C",
#>     primary = "#38404C",
#>     secondary = "#0046C8"
#>   ),
#>   h1("Hello World primary", class = "text-center text-secondary"),
#>   h1("Hello World secondary", class = "text-center text-primary")
#> )

Copy, paste into your app_ui.R, done.

And if you prefer clicking to typing, that is exactly what the hosted version does: swatches, a copy button, and a live Bootstrap mockup showing what your app would look like before you write a single line. Locally, css2r::run_app() launches the same thing.

Now look at that theme again. Look at it properly.

#F05622 is the ThinkR orange. The one on the logo, on the slides, on the mugs. It came third, with 40 occurrences. It is nowhere in the generated theme. What {css2r} picked as primary is #38404C, a dark slate grey: the body-text colour, which by construction appears everywhere.

We ran our own tool on our own website and it did not find our own brand colour. That deserves an explanation, and the explanation is the rest of this article.

Under the hood: nine steps and one regular expression

The css2r R6 class is deliberately small. The pipeline is linear and every step can be called by hand with on_initialize = FALSE:

1. check_internet(), via {curl}
2. download_html(), via {rvest}
3. extract_css_links(), every <link rel="stylesheet">
4. filter_css_links(), keeping only the ones on the same domain
5. download_css_files(), via {httr}
6. extract_colors(), the heart of the machine
7. analyze_colors(), splitting neutrals from brand colours
8. detect_google_fonts(), reading the Google Fonts URL parameters
9. generate_shiny_code(), assembling the bs_theme() call

Step 6, the one that does all the real work, is five lines:

pattern <- "#[0-9A-Fa-f]{6}"
matches <- gregexpr(pattern = pattern, text = self$css_content, perl = TRUE)
colors_hex <- regmatches(x = self$css_content, m = matches) |>
  unlist() |>
  toupper()
table_colors_hex <- colors_hex |> table() |> as.data.frame(stringsAsFactors = FALSE)

Find every six-digit hex code, count how many times each one appears, sort. Then analyze_colors() puts #FFFFFF and #000000 aside as neutrals, keeps the top four, and generate_shiny_code() assigns background and foreground by relative luminance:

hex_luminance = function(hex) {
  hex <- sub("^#", "", hex)
  r <- strtoi(substr(hex, 1, 2), 16L) / 255
  g <- strtoi(substr(hex, 3, 4), 16L) / 255
  b <- strtoi(substr(hex, 5, 6), 16L) / 255
  0.2126 * r + 0.7152 * g + 0.0722 * b
}

Lightest colour becomes bg, darkest becomes fg, most frequent becomes primary, second becomes secondary. That is the entire theory.

It is a heuristic, and a good one: it produces a usable theme on a large share of ordinary websites. But it rests on three assumptions about the web, and the modern web breaks all three.

Why this is harder than it looks

To make this concrete rather than theoretical, we pointed the extraction at a handful of real sites in August 2026 and counted what was actually reachable. The results are more instructive than any diagram.

Three distinct failure modes hide in that table.

1. The CSS you cannot reach

filter_css_links() keeps only stylesheets served from the same domain as the page:

keep(.p = ~ domain(.x) == self$domain)

That rule exists for a good reason. On thinkr.fr it correctly throws away a highlight.js theme hosted on cdnjs, a syntax-highlighting palette that has nothing to do with our brand and would have polluted the count with random purples.

On stripe.com, the same rule throws away everything. All five stylesheets live on b.stripecdn.com. Same company, different hostname, and {css2r} returns empty-handed:

External stylesheets only. This website loads 5 stylesheet(s), but all of them are served from another domain, so ShinyCopy cannot tell brand CSS from third-party CSS.

Any organisation that serves its assets from a dedicated CDN, which is to say most large organisations, is invisible to us.

There is a subtlety here that cost us a bug, and it is worth naming because it is the kind of mistake that hides behind a plausible error message. urltools::domain() does not return the registrable domain despite its name: it returns the full host, subdomain included. So for a site that redirects its apex to www, the page we downloaded and the stylesheets it references ended up on www.example.com while the filter was still comparing against the example.com the user had typed. Every stylesheet was discarded, and the app confidently blamed inline CSS. The filter now compares eTLD+1, and resolves everything against the URL actually reached after redirects, which keeps the useful behaviour (cdnjs still goes) without the self-inflicted wound.

Then there is the CSS that never appears in a <link> at all. thinkr.fr ships three inline <style> blocks, 8.3 kB of critical CSS injected straight into the HTML for faster first paint. It is a completely standard performance practice, and extract_css_links() walks right past it because it only looks at link[rel='stylesheet']. Shiny Copy at least has the decency to tell you:

No stylesheet found. No <link rel="stylesheet"> tag was found on thinkr.fr. This usually means the site ships its CSS inline in the HTML, or renders the page with JavaScript. ShinyCopy supports neither yet.

And finally, the pages that have no meaningful HTML at load time because everything is rendered client-side by JavaScript. {rvest} fetches the document as the server sent it, not as the browser eventually paints it. For a single-page app, there is often nothing to read.

2. The colours you cannot parse

This one is my favourite, because the tool fails silently and partially.

#[0-9A-Fa-f]{6} matches six-digit hex codes. Nothing else. On thinkr.fr that means we see 309 colour declarations and quietly ignore 234 three-digit hex codes (#fff, #333) and 70 rgb() / rgba() calls. Roughly forty percent of the colour information on our own site is dropped on the floor, and nothing in the output tells you so. Add hsl(), which appears 60 times on lemonde.fr, plus CSS named colours like rebeccapurple, and the blind spot widens further.

The fix is not conceptually hard: normalise every colour notation to hex before counting. It is just work that has not been done yet.

The deeper problem is CSS custom properties. Look at that var(--…) column again. On tailwindcss.com there are 8 371 references to CSS variables and on posit.co 1 545. A modern design system declares its brand colour exactly once:

:root { --brand-primary: #F05622; }

…and then uses it a thousand times through var(--brand-primary). Our frequency counter sees that colour once. Meanwhile a grey declared inline in fifteen legacy components scores fifteen. We are not measuring importance, we are measuring how badly the stylesheet has aged.

3. Frequency is not identity

Which brings us back to the ThinkR orange, and to the assumption underneath the whole approach: the most frequent colour is the most important colour.

It simply is not true. The most frequent colour in a stylesheet is almost always the interface plumbing: text grey, border grey, disabled grey. Brand colours are used sparingly, and that is the entire point of a brand colour. A logo orange that appears on three buttons carries far more identity than a #38404C used on every paragraph of the site.

lemonde.fr is the purest illustration, and a worse case than our own site. Every step of the pipeline succeeds, no error is raised, and this comes out:

$top_colors
    Color Count
1 #2A303C    77
2 #E8EAEE    62
3 #F5F6F8    40
4 #464F5F    25
bslib::bs_theme(bg = "#FFFFFF", fg = "#2A303C",
                primary = "#2A303C", secondary = "#E8EAEE")

All four are greys. Le Monde’s blue, #015AAD, sits at rank 10 with 17 occurrences: ranks one through nine are neutrals, every single one. And secondary lands on #E8EAEE, a contrast ratio of 1.09 against the white background. A .btn-secondary would be, quite literally, invisible.

On tailwindcss.com the top two are white (214) and black (129), and since analyze_colors() deliberately sets those aside, we fall through to #030712, which is… very slightly less black. Below that, 300 unique colours with no meaningful frequency gap between them, because a design system publishes its entire palette in one file.

There is one more consequence worth naming. On thinkr.fr, #38404C is both the most frequent colour and the darkest one, so it becomes fg and primary at the same time. A .text-primary element is then exactly the same colour as body text: technically valid, visually pointless. {css2r} performs no contrast check and no de-duplication between roles.

So we asked a human

Here is the pragmatic conclusion we reached while building Shiny Copy: if the machine cannot reliably tell a brand colour from interface plumbing, do not let it decide alone.

The app extracts the top four candidates, generates its best guess, and then displays clickable swatches for the primary and secondary roles. One click re-themes the whole preview live:

observeEvent(input$preview_primary, {
  rv$preview_primary <- input$preview_primary
  session$setCurrentTheme(
    bslib::bs_theme_update(
      rv$site$shiny_theme,
      primary   = rv$preview_primary,
      secondary = rv$preview_secondary
    )
  )
})

The frequency ranking stops being an answer and becomes a shortlist. You still get your ThinkR orange. You just have to click on it, which takes about a second, and you are the one who knows what your brand looks like.

Except that a shortlist is only ever as good as the ranking that produced it, and lemonde.fr shows exactly where that breaks: the four swatches on offer are four greys, and the brand blue at rank 10 is not among them. There is nothing left for the human to rescue. Deferring to the user fixes a ranking that is merely imperfect; it does nothing for one that is wrong from the first row.

What it would take to go further

The limits above are not mysteries, they are a roadmap. In rough order of payoff:

• Rank by chroma, not just by count. The cheapest fix with the largest visible effect. Discard near-neutrals before choosing primary, using nothing more than the spread between the RGB channels:

  chroma <- function(hex) {
    v <- c(strtoi(substr(hex, 2, 3), 16L),
           strtoi(substr(hex, 4, 5), 16L),
           strtoi(substr(hex, 6, 7), 16L))
    max(v) - min(v)
  }

  #38404C scores 18, #F05622 scores 206. The separation is not subtle. Applied to our two examples, primary / secondary would go from #38404C / #0046C8 to #0046C8 / #F05622 on thinkr.fr, and from two greys to #015AAD on lemonde.fr. It also happens to solve the primary == fg collision for free, since a dark neutral can no longer win.

• Normalise colour notations. A parse_css_color() helper turning #fff, rgb(), rgba(), hsl() and named colours into a canonical hex would immediately recover the ~40% we currently drop.

• Read inline <style> blocks. One extra html_nodes("style") |> html_text() and a lot of critical CSS stops being invisible.

• Handle asset CDNs. Comparing registrable domains fixed the apex/www case, but b.stripecdn.com is genuinely a different domain from stripe.com and no amount of string comparison will bridge that. It needs a different rule: an allowlist, or scoring stylesheets rather than excluding them outright.

• Resolve CSS custom properties. Parse --var: #hex declarations, then count var(--var) usages against them. This is the single biggest fidelity gain available, and the one that would make the tool work on modern design systems at all.

• Weight by context, not by count. A colour on a .btn background is not worth the same as a gradient stop in a hero image. That means a real CSS parser instead of a regular expression. A much bigger project, and the point where “small useful tool” becomes “product”.

• Check contrast before emitting. A quick WCAG ratio on the bg / fg pair, and a rule preventing primary from collapsing onto fg.

• Headless rendering for JavaScript-heavy pages, via {chromote}. Powerful, heavy, and a whole other class of deployment problem.

We shipped without any of these on purpose. A tool that solves the common case in five seconds and tells you honestly when it cannot is more useful than one that never leaves the lab. Every one of those error modals in Shiny Copy exists so that you find out immediately, rather than shipping a grey app and wondering why.

Try it, break it, tell us

The app is live and free, no installation required:

connect.thinkr.fr/css2r

And if you would rather have it in your own scripts:

remotes::install_github("ThinkR-open/css2r")
css2r::run_app()

Point it at your company website. If it works, you have saved yourself an afternoon in the brand book. If it fails, we would genuinely like to know which of the three failure modes above caught you. The source is on GitHub, issues and pull requests are very welcome, and the colour-notation parser is a lovely first contribution.

And if what you actually need is a Shiny application that is properly branded, properly tested and properly deployed, get in touch. That part we do by hand.

Figures measured in August 2026; websites change, so your numbers may differ.

This post is better presented on its original ThinkR website here: From URL to Theme: How {css2r} Steals a Website’s Colours, and Where It Gives Up

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

Have you ever heard someone claim that R is “too slow”?

As a Linux user who does almost everything on the command line, I used to find this idea hard to believe. R has always been quite efficient for my work. But recently, everything clicked.

While attending an R Sprint, I had to run R inside a Windows virtual machine to test that my suggested improvements to R worked there too. Linear models and some mathematical functions ran much slower compared to another Linux virtual machine setup.

Here is why: When you install Linux, something naive such as “pacman -S r” manages high-performance numerical libraries (like OpenBLAS) and makes sure everything is neat. On Windows, base R doesn’t ship with these optimizations out of the box, and setting them up requires an intricate separate installation.

A while back, Microsoft R Open solved this exact bottleneck by shipping with Intel MKL built-in. MRO is now retired and it closed a gap.

To bridge this gap again, I created R Batteries Included. It is a custom R distribution for Windows that ships pre-configured with OpenBLAS and Intel MKL. No messy setups and no extra steps. Just high-performance R on Windows, right out of the box.

G. N. Gueye suggested that my R distribution that makes R fast on Windows should be openly available with sponsorships. I followed that advice and R Batteries Included is here:

• Code
• Windows installer

See the plots showing the speed gains. I used the famous AT&T Benchmark.

Please re-post, try it, and add a star to the repository. If you can, consider sponsoring this work on:

• GitHub
• Buy me a Coffee

R-Community on Google Groups

I am in the process of creating a welcoming, collaborative space where anyone can ask and answer R questions without judgment and share R packages or their work using R. Think of it as the helpful, supportive sides of StackOverflow or Reddit, minus the negativity and elitism.

The idea of using Google Groups is to avoid making people creating accounts on different platforms, as these days the R community is spreaded between LinkedIn, X, Bluesky, Mastodon, and others I have not explored. Google Groups works with any email regardless if it is hosted by Google and you can interact with it from groups.google.com in the browser or from your email client.

No matter your programming style, all R preferences are welcome here. Whether you are working with base R, the Tidyverse, data.table, or exploring niche packages in the wider R ecosystem, your questions and insights are valuable. The idea is are all here to learn from one another.

Code of Conduct: Strictly friendly and welcoming environment. Zero-tolerance policy for bullying or hate.

Click here to join the community form: https://docs.google.com/forms/d/e/1FAIpQLSdMAj4adRAT4Gyuwt_9dPvxRvOUPml9AD59vuI7qS7XDlp48g/viewform

London R Users Group: Upcoming Social

How does meeting up during the second week of October sound?

I am planning an informal social event at a local Wetherspoons. It will be a completely relaxed evening to grab a pint (or a soft drink!), catch up in person, and chat about all things R. Whether you are a seasoned package developer or just starting your data journey, I would love to have you there.

If you are interested in coming along, please let me know:

• Email me directly
• Or join the discussion over on the R-Community Google Group

Two small R packages of mine have gone to CRAN. coinclp binds the COIN-OR Clp linear programming solver and was accepted on 15 September 2026; ROI.plugin.coinclp registers it with the R Optimization Infrastructure and is in the submission queue behind it. Between them they put back something R lost at the end of 2021, and this post is a short account of what that was, what the new packages do, and — because it would be dishonest to leave it out — when you should reach for a different solver altogether.

Where this comes from

Clp is one of the original COIN-OR codes. COIN-OR started inside IBM Research in 2000 as an open-source home for operations research software, and Clp — John Forrest’s simplex code, with a barrier method alongside — has been its linear programming workhorse ever since. It is mature, free under the Eclipse Public License, and good at the thing simplex codes are for: solving a model, changing it a little, and solving it again from the basis you already have.

R had a binding to it for a decade. Gabriel Gelius-Dietrich’s clpAPI, written at Heinrich Heine University Düsseldorf for the sybil metabolic-modelling toolkit, was on CRAN from 2011. In 2017 Benoit Thieurmel built ROI.plugin.clp on top of it, so that anyone using ROI could pick Clp with one argument. Then clpAPI’s CRAN checks started failing, nobody fixed them, and it was archived on 30 November 2021. ROI.plugin.clp had done nothing wrong, but it depended on an archived package, so it followed six weeks later.

I had been using clpAPI in my own models since 2018 and noticed the way everyone does: a fresh R install, a script that would not load. The particular annoyance is that on Windows the Clp library is already on the machine — Rtools has shipped it since version 4.3 — so the only thing missing was a few hundred lines of C++ to get at it.

What coinclp does

The bindings are written from scratch against the current Clp callable library, with registered entry points and external pointers that clean up after themselves, and they build on R 4.5 and 4.6. There are three ways in.

The first is one call. Give it an objective, a constraint matrix, directions and a right-hand side, and get back the solution, the objective, the shadow prices and the reduced costs together:

library(coinclp)

A <- rbind(material = c(120, 210),
           labour   = c(110,  30),
           capacity = c(  1,   1))

fit <- clp_solve(c(143, 60), A, "<=", c(15000, 4000, 75), max = TRUE)

fit$objval    # 6315.625
fit$solution  # 21.875 53.125
fit$duals     # 0.0000 1.0375 28.8750

The constraint matrix can be an ordinary dense matrix, a Matrix sparse matrix, a slam triplet matrix or plain i/j/v triplets. A sparse matrix goes to Clp as sparse: the package passes only the non-zero entries, and never expands the matrix into a full grid of mostly zeros first, which for a large model is the difference between fitting in memory and not.

The second is the whole callable library: build a model, keep it, change bounds or coefficients, hand back the basis and re-solve. That is what Clp is for, and it is where the one-call interfaces of most R solver packages let you down — a parametric study or a column-generation loop that rebuilds the model each iteration throws away exactly the information that makes simplex fast. In the vignette a tightened re-solve from a saved basis takes zero iterations.

The third is for old code. Every function clpAPI exported — initProbCLP(), loadProblemCLP(), solveInitialCLP() and the rest — is reproduced with the same names and arguments, so a script written against clpAPI needs only a new library() line. None of clpAPI’s code is reused; the layer is an independent implementation of its interface.

ROI.plugin.coinclp is the ROI side. The solver is called "coinclp" rather than "clp", because ROI takes the name from the package, and unlike the 2017 plugin it returns duals, reduced costs and row activities alongside the primal solution:

library(ROI)
library(ROI.plugin.coinclp)

res <- ROI_solve(op, solver = "coinclp")
solution(res, "dual")

When not to use it

Clp is not the fastest open-source LP solver any more, and I would rather say so than have you find out. The obvious comparison is HiGHS, from Julian Hall’s group at Edinburgh, which is now the default LP solver in SciPy and in MATLAB. On Hans Mittelmann’s LPopt benchmark as of September 2026, HiGHS solves 54 of the 65 test problems within the time limit and Clp solves 40, and HiGHS is about twice as fast on the scaled geometric mean. The commercial codes are further ahead again: COPT solves all 65 and is roughly 27 times faster than Clp on the same measure.

The benchmark also shows why the honest answer is “it depends on the model”. On a handful of instances Clp is the quicker of the two — Linf_520c takes Clp 36 seconds and HiGHS 872, and Clp solves datt256 where HiGHS times out — but on more of them the reverse holds, sometimes by a wide margin. If you have one large LP to solve once, try highs or ROI.plugin.highs first; with ROI the switch is one argument, so trying both costs nothing.

And Clp solves linear programs only. If your variables are integer you need a MIP solver: HiGHS again, or Rglpk, lpSolve or COIN-OR’s own Rsymphony. The ROI plugin will refuse an integer problem rather than quietly relax it, and clp_solve() has no notion of integer variables at all.

Where Clp still earns its place is the kind of work it was built for: models that are solved many times with small changes, where a warm start from the previous basis matters more than raw speed on a cold solve; anything that already speaks Clp’s API, which is more code than you might think; and Windows machines, where it is the one LP solver you get for free with Rtools and no further installation.

Installing it

On Windows, install.packages("coinclp") is the whole job — from source with Rtools installed until CRAN’s Windows binaries appear, which usually takes a few days. Elsewhere Clp is a system library and has to be there first: coinor-libclp-dev on Debian and Ubuntu, coin-or-Clp-devel on Fedora, brew install clp on macOS, coin-or-clp from conda-forge. Until the ROI plugin clears the CRAN queue it installs from GitHub:

install.packages("coinclp")
remotes::install_github("SamLovick/ROI.plugin.coinclp")

Source, issues and the vignette are at github.com/SamLovick/coinclp and github.com/SamLovick/ROI.plugin.coinclp. Both packages are under the Eclipse Public License, matching COIN-OR. If you were a clpAPI user and something in the compatibility layer does not behave as it used to, an issue with the script that broke would be very welcome.

The post coinclp: the COIN-OR Clp linear programming solver is back on CRAN appeared first on Sam Lovick Consulting.