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

Champions Program update

Our Champions are making great progress! They’ve now completed the training on good open source software development practices, package development, and peer review, and are moving on to explore community building and communications. Meanwhile, mentoring is underway with monthly meetings helping Champions move their projects forward. We’re already starting to see some exciting results: first versions of packages are taking shape, and some Champions will soon be sharing their work at LatinR! Stay tuned for more updates as their projects continue to grow!

R-Universe updates

R-Universe has started building and checking packages for Windows ARM64. Read more in our tech note.

R-Universe is now part of the infrastructure of the Bioconductor submission process.

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 first story features Eunseop Kim and his connection with rOpenSci. Read it on our blog: From Peer Review to Mentorship: My rOpenSci Story Stay tuned for more stories from our community as we continue celebrating 15 years of rOpenSci!

Coworking

Read all about coworking!

• Tuesday September 1st, 14:00 Europe Central (12:00 UTC) “Getting to Know SORTEE”, with Steffi LaZerte and co-host Ed Ivimey-Cook.

  • Visit SORTEE (Society for Open, Reliable, and Transparent Ecology and Evolutionary Biology).
  • Meet co-host, Ed Ivimey-Cook, and learn more about SORTEE and how you might get involved.

• 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) TBA, with Steffi LaZerte and co-host TBA.

• 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 four packages recently became a part of our software suite:

• lakefetch, developed by Jeremy Lynch Farrell: Calculates fetch (open water distance) and wave exposure metrics for lake sampling points. Downloads lake boundaries from OpenStreetMap, calculates directional fetch using a ray-casting approach, and optionally integrates National Hydrography Dataset (NHD) data https://www.usgs.gov/national-hydrography for hydrological context including outlet and inlet locations. Can estimate lake depth from surface area using empirical relationships, and integrate historical weather data for cumulative wave energy calculations. Includes an optional interactive shiny application for visualization. It has been reviewed by Jorrit Mesman and Kelly Hondula.

• galamm, developed by Øystein Sørensen: Estimates generalized additive latent and mixed models using maximum marginal likelihood, as defined in Sorensen et al. (2023) https://doi.org/10.1007/s11336-023-09910-z, which is an extension of Rabe-Hesketh and Skrondal (2004)s unifying framework for multilevel latent variable modeling https://doi.org/10.1007/BF02295939. Efficient computation is done using sparse matrix methods, Laplace approximation, and automatic differentiation. The framework includes generalized multilevel models with heteroscedastic residuals, mixed response types, factor loadings, smoothing splines, crossed random effects, and combinations thereof. Syntax for model formulation is close to lme4 (Bates et al. (2015) https://doi.org/10.18637/jss.v067.i01) and PLmixed’ (Rockwood and Jeon (2019) https://doi.org/10.1080/00273171.2018.1516541). It has been reviewed by Nicholas Clark and David Lawrence Miller.

• EpiStrainDynamics, developed by Saras Windecker together with Oliver Eales, James McCaw, and Freya Shearer: EpiStrainDynamics is a statistical framework developed for inferring temporal trends of multiple pathogens from routinely collected surveillance data. It has been reviewed by Sangeeta Bhatia and Joshua Lambert.

• RAMEN, developed by Erick I. Navarro-Delgado together with Keegan Korthauer and Michael S. Kobor: Using population data, RAMEN identifies which genetic (G), environmental (E), additive (G+E) or interaction (GxE) model better explains DNA methylation levels in genome-wide locations with high DNA methylation variability. It has been reviewed by Lluís Revilla Sancho and Ulduz Vafadarshamasbi.

Discover more packages, read more about Software Peer Review.

New versions

The following twenty-five packages have had an update since the last newsletter: DataSpaceR (v1.0.2), RAMEN (v2.1.1), frictionless (v1.3.0), cffr (v1.4.2), nodbi (v0.15.0), reviser (v0.2.0), writexl (v2.0.1), rangr (v1.0.10), lightr (v2.1.0), gert (v2.4.1), c14bazAAR (5.3.0), landscapetools (v0.6.3), comtradr (v1.0.6), lingtypology (v1.1.26), textreuse (v1.0.2), sofa (v0.4.2), pkgstats (v0.2.4), refsplitr (v1.2.3), npi (v0.3.0), rerddap (v1.3.0), stantargets (0.1.3), openalexR (v3.1.0), ernest (v1.2.5), mregions2 (v1.1.3), and spiro (v0.2.4).

Software Peer Review

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

• Five at ‘6/approved’:

  • EpiStrainDynamics, Infer temporal trends of multiple pathogens. Submitted by Saras Windecker. (Stats).

  • lakefetch, Calculate Fetch and Wave Exposure for Lake Sampling Points. Submitted by jeremylfarrell.

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

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

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

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

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

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

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

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

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

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

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

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

• Three 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.

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

• Four at ‘1/editor-checks’:

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

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

  • 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

Software Review

• From Peer Review to Mentorship: My rOpenSci Story by Eunseop Kim. From submitting a package, to reviewing one, to mentoring a Champion: my path with rOpenSci.

• FOSS Tools for Lazy Editors by Steffi LaZerte. How we streamlined the editing of our blog posts using 4 open-source tools that you could adopt too.

• The Journey of {nycOpenData}: From Classroom to Community by Christian Martinez.

• From Census Data to Demographic Analysis with ARcenso: A Reproducible Workflow in R by Andrea Gomez Vargas and Emanuel Ciardullo. How to Access and Process the 1970 and 1980 Argentine Censuses Using R. Other languages: De datos censales a análisis demográficos con ARcenso: un flujo de trabajo reproducible en R (es).

Tech Notes

• An API for Everything There Is to Know About Packages by Maëlle Salmon. Use cases of the R-Universe APIs.

• Windows ARM64 comes to R-universe by Jeroen Ooms.

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.

From Rd files to Quarto

Edgar Ruiz from Posit released pkgsite, a package for converting your package’s .Rd files to Quarto. It creates qmd files that you can integrate as you want in a Quarto website.

Mutation testing, fuzzy testing

First of all, a reminder in case you confuse the two concepts…

Mutation testing: you run tests on mutated version of the code.

Fuzzy testing: you run tests with all sorts of inputs to the code.

At useR! 2026, both topics were covered:

• mutation testing in Beyond Code Coverage: Mutation Testing in R with mutator by Assanali Amandykov, Pierre Donat-Bouillud.
• fuzzy testing in Fuzz-testing R-based research software for robustness by Marco Colombo.

roxyreqs

Also at useR!, Moritz Lang and collaborators introduced roxyreqs, a package for adding roxygen2-like documentation to testthat. Slides.

checktor, a new helper for CRAN submissions

If you want to submit your package to CRAN, you can get help through the CRAN cookbook, the collaborative list maintained by ThinkR and now a new package, checktor by James Balamuta! Read more in James’ post.

Interesting AI reads

• Recommendations When Using LLM-backed Generative AI Systems for FOSS Contributions by Software Freedom Conservancy, shared by Will Gearty.
• Who cleans up after the vibe-coding party? by Sam Learner in the Financial Times.
• AI Tooling and Open Source in which Nic Crane discusses “how AI tooling affects open source, the actions maintainers have been taking to address the less positive aspects, and emerging policies that open source projects are implementing around the topic of AI-generated pull requests”.

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!

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

Changes

• Removed the dplyr, tidyr, readr, tidyselect, tibble, rlang, and magrittr dependencies in favor of data.table (for row-binding and reshaping data) and base R.
• Removed the stringr and lubridate dependencies in favor of base R string and date handling.
• Functions now return data.table objects instead of tibbles.
• I added wb_country_coverage() to summarize the available data before downloading the series.
• I removed all the retired indicators and those that return empty tables from wb_cachelist.

Installation

You can install the latest release version from CRAN with

install.packages("wbstats")

which is the same as the latest development version from GitHub this as of 2026-08-28

remotes::install_github("pachadotdev/wbstats")

Downloading data from the World Bank

library(wbstats)

# Population for every country from 1960 until present
d <- wb_data("SP.POP.TOTL")
    
head(d)

#> # A tibble: 6 × 9
#>   iso2c iso3c country    date SP.POP.TOTL unit  obs_status footnote last_updated
#>   <chr> <chr> <chr>     <dbl>       <dbl> <chr> <chr>      <chr>    <date>      
#> 1 AF    AFG   Afghanis…  2024    42647492 <NA>  <NA>       <NA>     2025-07-01  
#> 2 AF    AFG   Afghanis…  2023    41454761 <NA>  <NA>       <NA>     2025-07-01  
#> 3 AF    AFG   Afghanis…  2022    40578842 <NA>  <NA>       <NA>     2025-07-01  
#> 4 AF    AFG   Afghanis…  2021    40000412 <NA>  <NA>       <NA>     2025-07-01  
#> 5 AF    AFG   Afghanis…  2020    39068979 <NA>  <NA>       <NA>     2025-07-01  
#> 6 AF    AFG   Afghanis…  2019    37856121 <NA>  <NA>       <NA>     2025-07-01

The current World Bank API does not provide summaries of data availability. I added the wb_country_coverage() function to supply that, which reads pre-computed summaries from GitHub.

d <- wb_country_coverage("gross domestic product", c("Mexico", "Chile"), 2010, 2020)
d

#      iso2c  iso3c country pct_complete  from    to  nobs         indicator
#     <char> <char>  <char>        <num> <int> <int> <int>            <char>
#  1:     CL    CHL   Chile         56.5  1960  2100    52     CC.EG.INTS.KW
#  2:     MX    MEX  Mexico         56.5  1960  2100    52     CC.EG.INTS.KW
#  3:     CL    CHL   Chile         34.8  1960  2025    23 EG.EGY.PRIM.PP.KD
#  4:     MX    MEX  Mexico         34.8  1960  2025    23 EG.EGY.PRIM.PP.KD
#  5:     CL    CHL   Chile         53.0  1960  2025    35 EG.GDP.PUSE.KO.PP
# ---                                                                       
# 82:     MX    MEX  Mexico         54.5  1960  2025    36    PA.NUS.PRVT.PP
# 83:     CL    CHL   Chile         53.0  1960  2025    35 SL.GDP.PCAP.EM.KD
# 84:     MX    MEX  Mexico         53.0  1960  2025    35 SL.GDP.PCAP.EM.KD
# 85:     CL    CHL   Chile        100.0  2004  2023    20   SPI.D5.2.5.HOUS
# 86:     MX    MEX  Mexico         40.0  2004  2023     8   SPI.D5.2.5.HOUS

# countries with less than 15% coverage for any variable
d[pct_complete < 15, ]

#     iso2c  iso3c country pct_complete  from    to  nobs                 indicator
#    <char> <char>  <char>        <num> <int> <int> <int>                    <char>
# 1:     MX    MEX  Mexico         12.9  1960  2100    13 UIS.XUNIT.GDPCAP.02.FSGOV

Hans Rosling’s Gapminder using wbstats

library(wbstats)
library(data.table)
library(tinyplot)

my_indicators <- c(
  life_exp = "SP.DYN.LE00.IN",
  gdp_capita ="NY.GDP.PCAP.CD",
  pop = "SP.POP.TOTL"
)

d <- wb_data(my_indicators, start_date = 2016)

d <- merge(d, wb_countries(), "iso3c")
d <- na.omit(d)

png(file="man/figures/readme-gdppc-vs-lifexp.png", width = 900, height = 600)
tinyplot(
  life_exp ~ gdp_capita | region,
  data = d,
  cex = d$pop,
  pch = 19,
  alpha = 0.7,
  palette = "tableau",
  log = "x",
  xaxl = "$",
  main = "An Example of Hans Rosling's Gapminder using wbstats",
  xlab = "GDP per Capita (log scale)",
  ylab = "Life Expectancy at Birth",
  cap = "Source: World Bank"
)
dev.off()

How I create the package data

Just in case this is useful.

What worked for me to query the data was to use the format and per_page arguments. You do not need this to work with the package.

Endpoints used for the package data:

• https://api.worldbank.org/v2/indicators?format=json&per_page=30000
• https://api.worldbank.org/v2/source?format=json&per_page=100
• https://api.worldbank.org/v2/topics?format=json&per_page=100
• https://api.worldbank.org/v2/regions?format=json&per_page=100
• https://api.worldbank.org/v2/incomelevel?format=json&per_page=10
• https://api.worldbank.org/v2/lendingtypes?format=json&per_page=10
• https://api.worldbank.org/v2/languages?format=json&per_page=100

Parts of this API endpoints description comes from https://dlthub.com/context/source/world-bank-indicators-api and other were just testing things like “language” and “languages” tp update wbstats. I did not create this package, I just assumed its maintenance as it is a very valuable resoure.

I added this because I did not find much information in the official documentation.

Naming Things by Tom Benner a tiny but neat book about the naming of identifies in code (variables, classes, methods, so not packages or libraries). It had entered my to-read list a few years ago, when I read the blog post Naming Things by Vicki Boykis. Tom Benner’s book is self published and available from Amazon (e-book and paperback) and Leanpub (e-book). I rarely buy stuff on Amazon, and prefer to read on paper, so I patiently waited for a copy to appear on my favorite second-hand book website.

Why read about naming things yet again?

At this point, I’ve been exposed to much advice on naming things, in The Art of Readable Code, A Philosophy of Software Design, The Programmer’s Brain… So why bother read yet another source of information on the topic? Well, I trusted Vicki Boykis’ recommendation, and since the book is so short – less than 100 pages, it wasn’t a dangerous bet.

The book is well organized, easy to read, and feels exhaustive. It explains why naming is important, why it is difficult, and presents 4 principles for naming: understandability, conciseness, consistency, distinguishability.

Here are some of my highlights…

Bad names, bad look

Among the numerous reasons why bad names are harmful for a project, this one caught my attention:

“[A] newcomer may develop a poor perception of the project and in the worst case, a poor perception of the team.”

What is an understandable name?

“An understandable name has high comprehension (it can be understood quickly) and high recall (it can be remembered easily).”

It reminds me of The Programmer’s Brain.

The book also recommends to avoid cleverness or irrelevant concepts: calling things based on some obscure joke or musical reference.

The ladder of abstraction

The book advises to use the “appropriate level of abstraction”.

“Do not use a name that’s so specific that you’re providing information that’s irrelevant to the audience, and do not use a name that’s so generic that it provides little or no relevant information to them.”

The book then discusses 4 names for a function that removes leading and trailing whitespace¹ from a phone number: process(), format(), trim_whitespace(), strip(). The right choice is explained to be format(): it shows the intent of the function without disclosing details that might be irrelevant or subject to change.

Booleans

The book recommends to always add is_ in the name of Booleans, e.g. is_valid.

It also states that they should be stated in the positive, with an example that I’m adapting to R below:

# Bad
if (!user_is_invalid) {
  save(user)
}

# Good

if (user_is_valid) {
  save(user)
}

This example resonated with me because it happens often to me to create a Boolean, use it with an if only to realize I should define the contrary of that Boolean instead.

And it reminds me, beyond naming, of negation-related rules in linters such as Jarl: comparison_negation, outer_negation.

The cost of renames

The book discusses the costs of a bad name (that add up over time: slow comprehension, low recall) and of a rename (one-time cost). It made me think of the renaming we did and do in igraph, including the batch renaming of functions with dots in them to snake-case equivalent (along with the correct lifecycle harness ): work for us but also for maintainers of reverse dependencies and direct users of the package.

Conclusion

Naming Things is a useful short read. After reading it, I feel I pay even more attention to names in the code I was writing of reviewing.

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

I build custom R packages and Shiny/Tabler apps, and I’m currently taking on new clients. You can find my profile on Fiverr.

Before starting my second master’s and my PhD, I ran a statistics consulting practice for nearly ten years. Most of that work fell into three buckets: designing and implementing SQL databases to streamline data analysis, building tailored R packages to simplify data access and reporting for teams, and building Shiny dashboards to summarise information through plots and KPIs.

With the rise of AI, this is a slightly harder pitch to make than it was a couple of years ago: sketching out a function or a small app is easier than ever. Where I think I still add real value is in the parts AI doesn’t do well on its own – turning a set of requirements into a maintainable package, designing something that your team will actually enjoy using, and making sure the whole thing is tested, documented, and installable by people who aren’t R experts. Besides it, I put a focus on writing everything with a minimal approach, meaning that I make an active effort in keeping the code logic simple and thinking about its long-term maintenance.

I’m also in the final year of my PhD, and between a still-unresolved pending payment from my previous university and a side project – a small guitar pedal business – that I started to cope with that financial emergency but that is not very profitable nor aligned with what I study.

I’m looking to pick up this kind of work again alongside my research. I also have a mobility disability (I get around with a cane, due to arthritis), which has made the usual part-time options like coffee shops or restaurants impractical, so consulting work I can do from a laptop is genuinely the best fit for me right now.

If your organization uses R, there are clear benefits to having an internal R package, whether you have a single R user or dozens. A package built around your organization’s specific needs opens up easier data access, shared functions for transformation and analysis, and a consistent look and feel across reports and dashboards.

Getting that first internal package off the ground can still feel daunting: what functions belong in it, how colleagues will install it and get updates, and how to keep quality consistent as more people contribute. This is exactly the kind of problem I like helping with – planning the package around your team’s actual workflow, fitting it into your existing infrastructure, and building out the core functions for data access, analysis, and reporting so the project has a solid foundation from day one.

Something similar can be said about dashboards, and I can help you to build something informative that keep quality consistent as more users help to improve it.

If any of this sounds useful, feel free to reach out through Fiverr or take a look at the packages above to get a sense of my work.

Below is a sample of the packages I’ve built over the years.

Data visualization

• d3po: A set of opinionated templates for quick data visualization using D3.js and R. It is fully compatible with RMarkdown and Shiny, and it is available under the Apache 2.0 license for use in commercial and non-commercial projects.
• tabler: A fully open-source alternative to Shiny worth considering if you need a multi-session, multi-user dashboard.

International Trade

• tradestatistics: Open trade Statistics API wrapper and utility program.
• wbstats: An R package for searching and downloading data from the World Bank API.

Econometrics

• capybara: Fast and memory efficient fitting of linear models with high-dimensional fixed effects.
• gravity: Estimation methods for gravity models.

R and C++ bindings

• cpp4r: Header-Only ‘C++’ and ‘R’ interface

Linear algebra

• armadillo4r: Provides function declarations and inline function definitions that facilitate communication between R and the Armadillo C++ library for linear algebra and scientific computing.

REDATAM format

• Open REDATAM (C++): Open Redatam is an open source software for extracting raw information from REDATAM databases. It was created to recover information of REDATAM databases for statistical analysis using standard tools such as SPSS, STATA, R, etc. It currently has both R and Python wrappers.

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

Linear Regression

The Expectation-Maximization (EM) algorithm is an iterative optimization framework used to find maximum likelihood estimates of parameters when a model depends on unobserved, latent variables.

For linear regression, the EM algorithm is unnecessary because Ordinary Least Squares (OLS) yields a closed-form analytical solution, which is \(\hat{\beta} = (X^T X)^{-1} (X^T y)\). However, framing linear regression through the EM algorithm is a quite clarifying exercise before jumping to Poisson (or binomial/logit) regression, models where data is missing (e.g., Tobit), or where the dataset is generated by a Mixture of Linear Regressions (multiple hidden lines).

To understand how EM applies, consider a dataset with \(n\) observations and \(p < n\) variables (or features). Each observation \(i\) corresponds to a vector \(x_i = (x_{i1}, x_{i2}, \ldots, x_{ip})^T \in \mathbb{R}^p\), a scalar response \(y_i\), and the model is

\[ y_i = \sum_{j = 1}^p \beta_j x_{ij} + e_i,\quad e_i \sim N(0, \sigma^2), \]

where \(\beta = (\beta_1, \beta_2, \ldots, \beta_p)^T \in \mathbb{R}^p\) are the weights to be estimated.

Stacking all \(n\) observations gives the compact matrix form \(y = X\beta + e\), where

\[ X = \begin{pmatrix} x_{11} & x_{12} & \ldots & x_{1p} \\ x_{21} & x_{22} & \ldots & x_{2p} \\ \vdots & \ddots & & \\ x_{n1} & x_{n2} & \ldots & x_{np} \end{pmatrix} \in \mathbb{R}^{n \times p}, \quad y = \begin{pmatrix} y_1 \\ y_2 \\ \vdots \\ y_n \end{pmatrix} \in \mathbb{R}^n, \quad e = \begin{pmatrix} e_1 \\ e_2 \\ \vdots \\ e_n \end{pmatrix} \sim N(0, \sigma^2 I_n). \]

In the EM framing we treat the unobserved (latent) variables \(z_i = \sum_{j = 1}^p \beta_j x_{ij}\) as the complete data, even though we can compute their distribution exactly. The point is to establish the machinery that generalises to other models.

Setup: complete-data log-likelihood

If we observed both \(y_i\) and \(z_i\), the complete-data log-likelihood for \(\beta\) and \(\sigma^2\) would be

\[ \ell_c(\beta, \sigma^2) = -\frac{n}{2}\log(2\pi\sigma^2) – \frac{1}{2\sigma^2}\sum_{i=1}^n (y_i – z_i)^2. \]

because \(y_i \mid z_i \sim N(z_i, \sigma^2)\) (the noise model), and \(z_i\) is deterministic given \(\beta\).

Expectation step

Given current parameter estimates \(\beta^{(t)}\) and \(\sigma^{2(t)}\), the E-step computes the expected complete-data log-likelihood with respect to the conditional distribution \(p(z \mid y, \beta^{(t)}, \sigma^{2(t)})\).

For linear regression the latent variable \(z_i\) is fully determined by \(\beta\). There is no uncertainty once \(\beta\) is fixed, so the conditional expectation is just the current fitted value

\[ \mathbb{E}\left[z_i \mid y_i, \beta^{(t)}\right] = \sum_{j = 1}^p \beta_j^{(t)} x_{ij}. \]

In matrix notation this is the \(i\)-th entry of \(X \beta^{(t)}\).

Substituting into the complete-data log-likelihood yields the Q-function

\[ Q(\beta, \sigma^2 \mid \beta^{(t)}, \sigma^{2(t)}) = -\frac{n}{2}\log(2\pi\sigma^2) – \frac{1}{2\sigma^2}\sum_{i=1}^n \left(y_i – \sum_{j = 1}^p \beta_j^{(t)} x_{ij}\right)^2. \]

Or equivalently, using vector notation for the residual vector \(y – X \beta^{(t)}\),

\[ Q(\beta, \sigma^2 \mid \beta^{(t)}, \sigma^{2(t)}) = -\frac{n}{2}\log(2\pi\sigma^2) – \frac{1}{2\sigma^2} \| y – X \beta^{(t)} \|^2. \]

The E-step collapses to plugging in the current fitted values. No integration is needed.

Maximization step

The M-step updates the parameters by maximising \(Q\) with respect to \(\beta\) (and \(\sigma^2\)).

Updating \(\beta\). The only term in \(Q\) that depends on \(\beta\) is the sum of squared residuals. Applying the chain rule to \(Q\) with respect to \(\beta_j\) gives

\[ \frac{\partial Q}{\partial \beta_j} = -\frac{1}{2\sigma^2} \sum_{i=1}^n 2\!\left(y_i – \sum_{k=1}^p \beta_k x_{ik}\right)(-x_{ij}) = \frac{1}{\sigma^2} \sum_{i=1}^n x_{ij}\!\left(y_i – \sum_{k=1}^p \beta_k x_{ik}\right). \]

Setting this to zero and multiplying through by \(\sigma^2\):

\[ \sum_{i=1}^n x_{ij} y_i = \sum_{i=1}^n x_{ij} \sum_{k=1}^p \beta_k x_{ik} = \sum_{k=1}^p \beta_k \underbrace{\sum_{i=1}^n x_{ij} x_{ik}}_{(X^T X)_{jk}}, \quad j = 1, \ldots, p. \]

The left-hand side is the \(j\)-th entry of \(X^T y\), since \((X^T y)_j = \sum_i x_{ij} y_i\). The right-hand side is the \(j\)-th entry of \(X^T X \beta\), since the \((j,k)\) entry of \(X^T X\) is exactly \(\sum_i x_{ij} x_{ik}\). Stacking all \(p\) equations (\(j = 1, \ldots, p\)) into a single matrix equation:

\[ X^T X \beta = X^T y. \]

where \(X^T X \in \mathbb{R}^{p \times p}\) is a symmetric positive-definite matrix (assuming the columns of \(X\) are linearly independent) and \(X^T y \in \mathbb{R}^p\) is a vector of inner products between each feature and the response. For illustration with \(p = 3\), these look like:

\[ \underbrace{\begin{pmatrix} \sum x_{i1}^2 & \sum x_{i1}x_{i2} & \sum x_{i1}x_{i3} \\ \sum x_{i1}x_{i2} & \sum x_{i2}^2 & \sum x_{i2}x_{i3} \\ \sum x_{i1}x_{i3} & \sum x_{i2}x_{i3} & \sum x_{i3}^2 \end{pmatrix}}_{X^T X} \begin{pmatrix} \beta_1 \\ \beta_2 \\ \beta_3 \end{pmatrix} = \underbrace{\begin{pmatrix} \sum x_{i1} y_i \\ \sum x_{i2} y_i \\ \sum x_{i3} y_i \end{pmatrix}}_{X^T y}. \]

Inverting \(X^T X\) gives the OLS formula

\[ \beta^{(t+1)} = (X^T X)^{-1} X^T y. \]

Note that this solution is independent of the current iterate \(\beta^{(t)}\), so the algorithm converges in a single M-step regardless of the initialisation. This is consistent with the fact that OLS has a closed-form solution.

Updating \(\sigma^2\). Setting \(\partial Q / \partial \sigma^2 = 0\):

\[ \sigma^{2(t+1)} = \frac{1}{n} \sum_{i=1}^n \left(y_i – \sum_{j=1}^p \beta_j^{(t+1)} x_{ij}\right)^2 = \frac{1}{n}\|y – X \beta^{(t+1)}\|^2. \]

which is the mean squared residual at the updated weights.

Convergence

Because the M-step yields the global maximum of the Q-function in closed form and that maximum does not depend on \(\beta^{(t)}\), the sequence \(\{\beta^{(t)}\}\) reaches \(\hat{\beta} = (X^T X)^{-1} X^T y\) after exactly one iteration. This is the standard EM convergence property specialised to the case where the complete-data problem is a convex problem with a unique global maximum.

Poisson Regression

Poisson regression models count data. As an aside, international trade models rely on Poisson pseudo maximum likelihood (PPML) and a continuous variable such as exports (or imports) does not follow a discrete Poisson distribution, which is why the PPML and not PML name. The PPML estimator is consistent if the conditional mean of the variate of interest is correctly specified. More on that on The Log of Gravity page.

The response \(y_i \in \{0, 1, 2, \ldots\}\) is assumed to follow a Poisson distribution whose mean depends on the covariates through a log-link:

\[ y_i \sim \text{Poisson}(\mu_i), \quad \mu_i = \exp\!\left(\sum_{j=1}^p \beta_j x_{ij}\right) = \exp(x_i^T \beta). \]

Unlike linear regression there is no closed-form solution for \(\beta\), so we need an iterative method. The EM algorithm provides one by introducing latent variables that make the complete-data problem tractable.

Latent-variable construction

Write the Poisson mean as \(\mu_i = \exp(x_i^T \beta)\) and introduce \(m\) latent binary indicators \(z_{i1}, \ldots, z_{im}\). This is one for each of \(m\) hypothetical sub-processes that together generate \(y_i\). Specifically, partition \(\mu_i\) into \(m\) equal parts \(\lambda = \mu_i / m\) and let

\[ z_{il} \sim \text{Bernoulli}(\lambda / (1 + \lambda)), \quad l = 1, \ldots, m, \]

so that \(y_i = \sum_{l=1}^m z_{il}\) in the limit \(m \to \infty\).

In practice the standard EM formulation for Poisson regression avoids this explicit construction and instead treats the complete data as the pair \((y_i, \eta_i)\), where \(\eta_i = x_i^T \beta\) is the linear predictor, and exploits the exponential-family properties that apply to Poisson log-likelihood.

For more on the exponential family and statistical sufficiency, you can check Casella and Berger. It is one of my favourite books (unlike others that put elegance over clarity).

Setup: complete-data log-likelihood

The Poisson log-likelihood for a single observation is

\[ \log p(y_i \mid \beta) = y_i \log \mu_i – \mu_i – \log(y_i!) = y_i (x_i^T \beta) – \exp(x_i^T \beta) – \log(y_i!). \]

Summing over all \(n\) observations gives the complete-data log-likelihood (dropping the constant \(\sum_i \log(y_i!)\)):

\[ \ell(\beta) \propto \sum_{i=1}^n \left[ y_i (x_i^T \beta) – \exp(x_i^T \beta) \right] = y^T X \beta – \mathbf{1}^T \exp(X\beta), \]

where \(\exp(X\beta)\) denotes element-wise exponentiation and \(\mathbf{1}\) is a vector of ones.

Expectation step

Unlike linear regression, the Poisson log-likelihood is not quadratic in \(\beta\), so the E-step does not collapse trivially. The standard approach is to construct a working quadratic surrogate (the Q-function) at the current iterate \(\beta^{(t)}\) using a second-order Taylor expansion of \(\exp(x_i^T \beta)\) around \(\eta_i^{(t)} = x_i^T \beta^{(t)}\):

\[ \exp(x_i^T \beta) \approx \exp(\eta_i^{(t)}) + \exp(\eta_i^{(t)})(x_i^T \beta – \eta_i^{(t)}) + \frac{1}{2}\exp(\eta_i^{(t)})(x_i^T \beta – \eta_i^{(t)})^2. \]

As an aside, Taylor expansions provide the foundation for Newton’s Method. Both are used heavily in industry, and a famous example is Quake’s III Fast Inverse Square Root.

Substituting into \(\ell(\beta)\) and keeping only terms that depend on \(\beta\) gives the Q-function

\[ Q(\beta \mid \beta^{(t)}) \propto -\frac{1}{2} \sum_{i=1}^n \mu_i^{(t)} \!\left(x_i^T \beta – \eta_i^{(t)} – \frac{y_i – \mu_i^{(t)}}{\mu_i^{(t)}}\right)^{\!2}, \]

where \(\mu_i^{(t)} = \exp(\eta_i^{(t)})\).

Defining the working response

\[ \tilde{y}_i^{(t)} = \eta_i^{(t)} + \frac{y_i – \mu_i^{(t)}}{\mu_i^{(t)}} \]

and the weight \(\mu_i^{(t)}\), the Q-function becomes

\[ Q(\beta \mid \beta^{(t)}) \propto -\frac{1}{2} \sum_{i=1}^n \mu_i^{(t)} \!\left(\tilde{y}_i^{(t)} – x_i^T \beta\right)^2. \]

This is exactly a weighted least-squares objective. It is the same structure as the linear regression log-likelihood with complete data, but with observation-specific weights \(\mu_i^{(t)}\).

Maximization step

Maximising the weighted-least-squares Q-function with respect to \(\beta\) is the same calculation as in the linear regression M-step, but with a diagonal weight matrix \(W^{(t)} = \text{diag}(\mu_1^{(t)}, \ldots, \mu_n^{(t)}) \in \mathbb{R}^{n \times n}\).

Updating \(\beta\). The weighted normal equations are obtained by the same entry-wise gradient argument as before. For each \(j = 1, \ldots, p\):

\[ \frac{\partial Q}{\partial \beta_j} = \sum_{i=1}^n \mu_i^{(t)} x_{ij}\!\left(\tilde{y}_i^{(t)} – x_i^T \beta\right) = 0, \]

which in matrix form is

\[ X^T W^{(t)} X \beta = X^T W^{(t)} \tilde{y}^{(t)}. \]

Here \(X^T W^{(t)} X\) is symmetric positive-definite (it is the weighted Gram matrix), and \(X^T W^{(t)} \tilde{y}^{(t)} \in \mathbb{R}^p\). Solving gives

\[ \beta^{(t+1)} = \left(X^T W^{(t)} X\right)^{-1} X^T W^{(t)} \tilde{y}^{(t)}. \]

This is the Iteratively Reweighted Least Squares (IRLS) update, which is the standard algorithm for fitting Poisson (and other GLM) models. Each M-step is a weighted OLS problem with the same structure as the linear-regression result, but the weights \(W^{(t)}\) and working responses \(\tilde{y}^{(t)}\) change at every iteration as \(\mu_i^{(t)}\) is updated.

Convergence

Because the Poisson log-likelihood \(\ell(\beta)\) is strictly concave in \(\beta\) (the Hessian \(-X^T W X\) is negative-definite if \(X\) has full column rank), the sequence \(\{\beta^{(t)}\}\) converges to the unique maximum-likelihood estimate \(\hat{\beta}\). Unlike linear regression, convergence requires multiple iterations. The algorithm terminates when \(\|\beta^{(t+1)} – \beta^{(t)}\| < \varepsilon\) for a chosen tolerance \(\varepsilon\).

More about Poisson and other models

Check McCullagh and Nelder. It covers Poisson, Logit, and Generalized Linear Models with a very detailed treatment. This is another of my favourite books.

Overview

      Do people in red states live shorter lives?

      This project examines whether state level life expectancy is statistically associated with each state’s political climate in the 2024 presidential election. Political climate is measured using the popular vote margin, defined as (Trump votes – Harris votes) / (Trump votes + Harris votes). A positive value indicates a Republican advantage, and a negative value indicates a Democratic advantage.

      I’m not predicting elections, and I’m not making a claim of causality. There are many reasons why two states differ in life expectancies – differences in average income, availability of medical care, occupations with differing job hazards, etc. I’m just asking whether these two measurable state‑level quantities move together.

      A scatterplot of life expectancy versus vote margin shows a clearly downward trend: states with higher Republican vote margins tend to have lower life expectancy, and states with higher Democratic margins tend to have higher life expectancy.

      The correlation coefficient r is -0.50. This represents a moderate negative relationship – neither weak nor strong, but unmistakenly present. A significance test yields t = -4.041, p = .00019. With 51 observations, this correlation is statistically significant at the .05 level.

Details

      I obtained 2024 presidential percentage popular vote data from the Federal Election Commission. Trump won 31 of 51 (50 plus DC) states in 2024. To visualize the distribution of states, I grouped the vote margin varaible into four bins:

• Dem Majority: margin ≤ -10%
• Dem Narrow: -10% < margin ≤ 0
• Rep Narrow: 0 < margin ≤ 10%
• Rep Majority: margin > 10%

24 of Trump’s 31 states were by margins greater than 10%. A bar chart and a US map show how many states fall into each group and where they are located geographically.

      The CDC (Centers for Disease Control and Prevention) publishes life expectancy tables by state. These are period life tables, showing the life expectancy of a newborn under today’s mortality rates, assuming the age‑specific death rates observed in that year (e.g., 2022) stay fixed for the newborn’s entire lifetime. Hawaii has the highest life expectancy at 80.0 years, and Mississippi has the lowest at 70.9 years.

      The following map shows life expectancies by state. I allocated the states by quartile (shortest life expectancy, shorter, longer, longest). I believe there is a relationship especially with southern states having short life expectancies in this map, compared with Republican margins in the prior map.

      The actual correlation coefficient is r = -0.50, which is moderate, but statistically significant.

      Incidentally, I had a little challenge drawing the maps with R library usmap. That library includes Puerto Rico which was not in the voter or life expectancy data, and the map would show Puerto Rico as an NA until I excluded it in the plot_usmap statement.

R code

library(readxl)
pres <- read_excel("C:/Users/Jerry/Desktop/R_files/2024presgeresults.xlsx", n_max=51)
pres$TRUMP_PERCENT <- pres$TRUMP/(pres$TRUMP + pres$HARRIS)   # ratio of popular votes
pres$HARRIS_PERCENT <- pres$HARRIS/(pres$TRUMP + pres$HARRIS)
pres$TRUMP_MARGIN <- round(pres$TRUMP_PERCENT - pres$HARRIS_PERCENT,3)
pres <- pres[, c("STATE", "TRUMP_MARGIN")]   # states are 2 letter abbrevs
colnames(pres)[1] <- "state"  # usmap requires state Column Name to be lowercase "state"

# CDC life expectancies by state:   https://www.cdc.gov/nchs/data/nvsr/nvsr74/nvsr74-12.pdf

library(pdftools)   # extract text from pdf file
library(tidyverse)
raw_text <- pdf_text("C:/Users/Jerry/Desktop/R_files/nvsr74-12.pdf")
page_text <- raw_text[3]   # page 3 only
lines <- read_lines(page_text)

# data cleaning of life expectancy file:
clean_lines <- lines %>%
  str_trim() %>%                      # Remove leading/trailing spaces
  keep(~ .x != "")                    # Drop empty rows

# Extract column headers (row 6 contains headers):
headers <- str_split(clean_lines[6], "\\s{2,}")[[1]]  # these are partial headers

# Process the data rows (Rows 7 to the end):
data_rows <- clean_lines[7:(length(clean_lines)-3)]   # delete footnotes

# Convert text lines into data frame:
life_exp <- data_rows %>%
  # Split columns whenever there are 2 or more spaces
  str_split_fixed("\\s{2,}", n = length(headers)) %>%
  as_tibble(.name_repair = "minimal")

colnames(life_exp) <- c("State", "Tot_Rank", "Tot_LE", "Tot_SE", "Male_Rank", "Male_LE", "Male_SE",
     "Fem_Rank", "Fem_LE", "Fem_SE")
life_exp$State <- gsub("\\.", "", life_exp$State)    # delete periods
life_exp$State <- sub("\\s+$", "", life_exp$State)   # delete spaces after last char 
life_exp <- subset(life_exp, State != "United States")
# convert states from names to abbreviations; District of Columbia will be NA without next line:
life_exp$State <- c(state.abb, "DC")[match(life_exp$State, c(state.name, "District of Columbia"))]   # Convert full name to 2-letter abbreviation
life_exp <- life_exp %>%
  mutate(across(where(is.character) & -1, as.numeric))   # converts all character columns in a data frame into numeric columns, except for the very first column
life_exp <- life_exp[, c("State", "Tot_LE")]   # states are 2 letter abbrevs
colnames(life_exp)[1] <- "state"  # usmap requires state Column Name to be lowercase "state"
print(life_exp)

df <- merge(pres, life_exp, by = "state")

########  Display summaries:  ########

library(ggplot2)

common_theme <- theme(
        plot.title = element_text(size=15, face="bold"),
        plot.subtitle = element_text(size=12.5, face="bold"),
        axis.title = element_text(size=15, face="bold"),
        axis.text = element_text(size=15, face="bold"),
        legend.title = element_text(size=15, face="bold"),
        legend.text = element_text(size=15, face="bold"))

df <- df %>%
  mutate(TRUMP_MARGIN_RANGE = case_when(  
    TRUMP_MARGIN  -.10 & TRUMP_MARGIN  0 & TRUMP_MARGIN  .10 ~ "Rep Majority",
    TRUE ~ NA_character_ 
  )
)

percent_colors <- c("Dem Majority" = "#883068", "Dem Narrow" = "#4292C6", 
                    "Rep Narrow" = "#FB6A4A", "Rep Majority" = "#CB181D")
df$TRUMP_MARGIN_RANGE <- factor(
  df$TRUMP_MARGIN_RANGE, 
  levels = c("Dem Majority", "Dem Narrow", "Rep Narrow", "Rep Majority")
)

ggplot(df, aes(x = TRUMP_MARGIN_RANGE)) +
  geom_bar(fill = percent_colors) +
  geom_text(
    stat = "count", 
    aes(label = after_stat(count)),
    fontface = "bold", 
    vjust = -0.5
  ) +  
  labs(title="2024 Presidential Election Results by Vote Margin",
       y = "Number of States", x = "% Popular Vote Margin") +
  guides(fill = guide_legend(title = NULL)) + 
  scale_x_discrete(
    labels = c(
      "Dem Majority" = "Dem + 10% or more",
      "Dem Narrow" = "Dem 0 - 10%",
      "Rep Narrow" = "Rep 0 - 10%",
      "Rep Majority" = "Rep + 10% or more")) +
  common_theme

colSums(is.na(df))  # 0
colnames(df)[1] <- "state"   # usmap requires state Column Name to be lowercase "state"
length(df$state)   # 51
df$state <- trimws(toupper(df$state))   #51 states including DC

library(usmap)
unique(usmap::us_map(regions = "states")$full)  # Includes Puerto Rico
plot_usmap(data = df, , regions = "states", values = "TRUMP_MARGIN_RANGE", exclude = "Puerto Rico") +
  labs(title="2024 Presidential Election Results by Vote Margin") +
  scale_fill_manual(values=percent_colors,
    guide = guide_legend(title = NULL, direction="vertical"),
    labels = c(
      "Dem Majority" = "Dem + 10% or more",
      "Dem Narrow" = "Dem 0 - 10%",
      "Rep Narrow" = "Rep 0 - 10%",
      "Rep Majority" = "Rep + 10% or more")
    ) +
  theme(
    legend.position = "bottom",
    legend.box = "vertical",
    plot.title = element_text(size=15, face="bold"),
    legend.title = element_text(size=12, face="bold"),
    legend.text = element_text(size=12, face="bold")
  )

df <- df %>%
  mutate(
    LE_RANGE = case_when(
      ntile(Tot_LE, 4) == 1 ~ "Shortest LE",
      ntile(Tot_LE, 4) == 2 ~ "Shorter LE",
      ntile(Tot_LE, 4) == 3 ~ "Longer LE",
      ntile(Tot_LE, 4) == 4 ~ "Longest LE"
    )
  )

table(df$LE_RANGE)

LE_colors <- c(
  "Shortest LE" = "#D1E5F0",
  "Shorter LE"    = "#92C5DE",
  "Longer LE"     = "#4393C3",
  "Longest LE"  = "#8B3C59"
)

df$LE_RANGE <- factor(
  df$LE_RANGE, 
  levels = c("Shortest LE", "Shorter LE", "Longer LE", "Longest LE")
)

plot_usmap(data = df, regions = "states", values = "LE_RANGE", exclude = "Puerto Rico") +
  labs(title="Life Expectancies by State") +
  scale_fill_manual(values=LE_colors,
     guide = guide_legend(title = NULL, direction="vertical")) +
  theme(
    legend.position = "bottom",
    legend.box = "vertical",
    plot.title = element_text(size=15, face="bold"),
    legend.title = element_text(size=12, face="bold"),
    legend.text = element_text(size=12, face="bold")
  )


########  Corr coeff:  ########

library(ggrepel)
ggplot(data = df, mapping = aes(x = TRUMP_MARGIN, y = Tot_LE)) +
  geom_point(color = "#4D4D4D") +
  geom_smooth(method = "lm", color = "steelblue", se = FALSE, linewidth = 1) +
  geom_text_repel(
    data = df[df$state %in% c("HI", "WV"), ],
    aes(x = TRUMP_MARGIN, y = Tot_LE, label = state),
    size = 4, color = "black", fontface = "bold", nudge_y = 0.25
  ) +
  ggtitle("Life Expectancy vs % Popular Vote Margin by State") +
  xlab("Pop Vote Margin %:  Rep +, Dem - ") +
  ylab("Life Expectancy") + 
  annotate(
    "text",
    x = min(df$TRUMP_MARGIN) + .02,
    y = max(df$Tot_LE) - .02,
    label = "r = -0.50",
    hjust = 0, vjust = 1,
    color = "gray20", fontface = "bold", size = 4) +
  common_theme

r <- round(cor(df[, c("Tot_LE", "TRUMP_MARGIN")])[1, 2], 3)  # -.50 is moderate (neither weak nor strong)

# statistical significance of correlation coefficient:
n <- nrow(df)
dof <- n - 2
t <- round(r*sqrt(n - 2) / sqrt(1 - r^2),3)   # - 4.041
alpha <- .05
p_value <- round(2 * pt(abs(t), df = dof, lower.tail = FALSE), 5)    # .00019
stat_signif <- if(p_value < alpha, "is statistically significant", "is not statistically significant")
cat("r = ", r, "p-value = ", p_value, stat_signif)

End

We’ve previously extolled the virtues of automating the repetitive chores we encounter, allowing us to focus on the tasks that matter. Three years ago, we posted an example of how to automate the deployment of a Shiny application to shinyapps.io from a GitHub Actions. When a new commit was pushed to the GitHub repository, it triggered an automated pipeline to bundle the Shiny application source code and send it to shinyapps.io. You could relax, knowing that the production deployment was always up-to-date.

But in January 2026, Posit announced that shinyapps.io would be closing to new apps at the end of 2026, with all users moving to Posit Connect Cloud. Existing content on shinyapps.io will continue to work before automatically migrating across in early 2027. But if you followed our previous methodology for automating that deployment from a GitHub Actions workflow, you’ll need to adjust your deployment strategy.

Why move to Posit Connect Cloud?

shinyapps.io has provided a faithful service to the Shiny community for a long time. You make your Shiny application, you click “Deploy” in RStudio, some magic happens, then your work appears online for others to access. You didn’t have to think too hard about R packages, build a Docker container, or set up a cloud compute instance to grant access to your app. It was perhaps the simplicity of the deployment process and availability of a free-tier that made it so popular.

Posit Connect Cloud takes what made shinyapps.io so popular, and stacks more features and convenience on top. You’re no longer restricted to just hosting Shiny applications—Posit Connect Cloud can also support Streamlit, Bokeh, Jupyter Notebooks, and all plans allow unlimited hosting for rendered Quarto and R Markdown documents¹. You also get more functionality: the ability to set secret variables, regenerate content on a schedule, higher maximum compute limits, and SSL certificates when using custom domains.

A free-tier of Posit Connect Cloud also remains. And while Posit Connect Cloud supports more types of content beyond just Shiny applications, you will find that benefit is reflected in the higher prices on paid-tiers over their nearest shinyapps.io equivalents. Perhaps the biggest winners are users who mainly just needed a custom domain: This required the highest $349/month “Professional” tier on shinyapps.io, but is now available (with SSL certificate) on the $59/month² “Enhanced” tier and above on Posit Connect Cloud.

Deploying content to Posit Connect Cloud

The existing deployment methods used for shinyapps.io still work with Posit Connect Cloud, but there are some new options too:

• The Quarto CLI can deploy to Posit Connect Cloud using the quarto publish command.

• You can grant Posit Connect Cloud access to your GitHub account to perform a Git-backed deployment, where it will monitor the code for changes and automatically re-deploy when the target branch is updated. The only extra step you need to do is to commit and push a manifest.json file, which is often as simple as running the following in R:

  rsconnect::writeManifest()
  

The Git-backed deployment may be a very useful replacement to those who have previously deployed content to shinyapps.io from GitHub Actions; The deployment work is now handled by Posit Connect Cloud rather than using up your GitHub Actions allowance. If that method works for you, it’s what we’d now recommend in most cases. But there are some circumstances where you may still want to automate deployment from your own CI/CD process:

1. You only want deployment to happen after earlier pipeline checks are successful.
2. You have sensitive code elsewhere in your GitHub account, and don’t feel comfortable or aren’t allowed to grant access to your GitHub repositories to a third-party tool.
3. You’re on the free-tier of Posit Connect and have code in a Private GitHub repository.
4. You want to include extra resources that aren’t stored in the GitHub repository, such as a moderately-sized read-only dataset that rarely updates. Here you might want to use the GitHub Actions workflow to pull external resources together and create a fully self-contained application bundle of source code and static data. This can reduce data export costs on busy applications.

Deployment to Posit Connect Cloud using GitHub Actions

So you might want to automate deployment, but not be able to use the standard Git-backed deployment methods. Let’s discuss how to make it work.

Obtain a Content ID

The first thing you’ll want to do is perform an initial deployment of the application so that we have a Content ID. The easiest way to do this is using the one-click deployment method from RStudio, or through the Posit Publisher extension in Positron or VS Code. Log in to your Posit Connect Cloud account and find the content in your list. In the Settings menu, go to URL and look at the “Default URL”. It should contain a UUID-like section after the https:// and before the .share.connect.posit.cloud parts—we want to make a note of this for later.

The “Default URL” in the settings menu for your content contains the Content ID. In this example, the default URL is https://019eb78b-0c21-3b55-3fe6-38ae4d03dee4.share.connect.posit.cloud, so the Content ID is 019eb78b-0c21-3b55-3fe6-38ae4d03dee4.

This manual initial deployment is also a good opportunity to ensure that the deployed application is working in the first place—if your app doesn’t work when deployed from your IDE, then it’s unlikely to work when the same stages are performed in a GitHub Actions workflow.

Add an renv.lock file

You’ll also want to maintain an {renv} lockfile to record which packages were used during development. These matching packages will be used in the deployed version of the application for maximum compatibility. We’ll ask {renv} to use {pak} when restoring these R packages during the GitHub Actions workflow—{pak} is generally faster at package installation and can automatically install all the system dependencies needed for the packages. Remember to commit and push the renv.lock and other relevant {renv}-related files to the Git remote.

Create a Posit Connect Cloud token

Your Posit Connect Cloud account is part of your larger Posit Cloud account. In Posit Cloud, you can access a list of your “Credentials”, which are access tokens. These can be found at https://login.posit.cloud/identity/credentials.

You’ll have the option to create “New Credentials”.

Provide a name for the new token that helps identify where it will be used, then for the “Use with” option, select “Connect Cloud”. Click “OK” to generate a token.

You’ll be presented with a block of R code containing the credentials you can use to log in. These should be kept secret; Anyone with these details is able to impersonate you. It should resemble this:

rsconnect::connectCloudClientCredentials(
 clientId="01234567-89a1-b2c3-d4e5-f60123456789",
 clientSecret="SuPeR/SeCrEt/VeRy/LoNg/CoDe",
 account="<YOUR_ACCOUNT_HERE>"
)

We’ll need these when we come to set GitHub Actions variables and secrets later.

Write a GitHub Actions Workflow

We’ll be creating a GitHub Actions workflow with a number of stages. In the root of our Git repository, we’ll make a file at .github/workflows/deploy.yml.

# .github/workflows/deploy.yml
name: Deploy to Posit Connect Cloud

on:
 push:
 branches:
 - main
 - master
 workflow_dispatch:

jobs:
 deploy:
 runs-on: ubuntu-latest

 steps:
 - name: Checkout repository
 uses: actions/checkout@v4

 - name: Setup R
 uses: r-lib/actions/setup-r@v2
 with:
 r-version: "4.6.1"

 - name: Install pak
 run: |
 UBUNTU_CODENAME=$(lsb_release -cs)
 Rscript -e "
 install.packages('pak', repos = 'https://packagemanager.posit.co/cran/__linux__/${UBUNTU_CODENAME}/latest');
 "

 - name: Instruct renv to use pak in .Rprofile
 run: |
 echo 'options(renv.config.pak.enabled = TRUE)' >> .Rprofile

 - name: Restore packages from renv.lock file
 uses: r-lib/actions/setup-renv@v2

 - name: Install rsconnect if not present in renv.lock
 run: |
 Rscript -e "if (!requireNamespace('rsconnect', quietly = TRUE)) pak::pak('rstudio/rsconnect')"

 - name: Authenticate with rsconnect
 run: |
 Rscript -e '
 rsconnect::connectCloudClientCredentials(
 clientId = Sys.getenv("RSCONNECT_CLIENT_ID"),
 clientSecret = Sys.getenv("RSCONNECT_CLIENT_SECRET"),
 accountName = Sys.getenv("RSCONNECT_USERNAME"),
 name = NULL
 )
 '
 env:
 RSCONNECT_CLIENT_ID: ${{ secrets.RSCONNECT_CLIENT_ID }}
 RSCONNECT_CLIENT_SECRET: ${{ secrets.RSCONNECT_CLIENT_SECRET }}
 RSCONNECT_USERNAME: ${{ vars.RSCONNECT_USERNAME }}

 - name: Add Posit Connect deployment config file
 run: |
 mkdir -p "rsconnect/${SERVER}/${RSCONNECT_USERNAME}"
 cat > "rsconnect/${SERVER}/${RSCONNECT_USERNAME}/${APP_NAME}.dcf" <<EOF
 name: ${APP_NAME}
 title: ${APP_TITLE}
 username: ${RSCONNECT_USERNAME}
 account: ${RSCONNECT_USERNAME}
 server: ${SERVER}
 hostUrl: https://api.${SERVER}/v1
 appId: ${CONNECT_CONTENT_ID}
 EOF
 env:
 SERVER: connect.posit.cloud
 APP_NAME: ${{ vars.APP_NAME }}
 APP_TITLE: ${{ vars.APP_TITLE }}
 CONNECT_CONTENT_ID: ${{ vars.CONNECT_CONTENT_ID }}
 RSCONNECT_USERNAME: ${{ vars.RSCONNECT_USERNAME }}

 - name: Deploy to Posit Connect Cloud
 run: |
 Rscript -e '
 rsconnect::deployApp(
 appDir = ".",
 appId = Sys.getenv("CONNECT_CONTENT_ID"),
 appTitle = Sys.getenv("APP_TITLE"),
 logLevel = "verbose",
 account = Sys.getenv("RSCONNECT_USERNAME"),
 forceUpdate = TRUE
 )
 '
 env:
 CONNECT_CONTENT_ID: ${{ vars.CONNECT_CONTENT_ID }}
 APP_TITLE: ${{ vars.APP_TITLE }}
 RSCONNECT_USERNAME: ${{ vars.RSCONNECT_USERNAME }}

 - name: Clean up account details
 run: |
 Rscript -e '
 rsconnect::removeAccount(
 name = Sys.getenv("RSCONNECT_USERNAME")
 )
 '
 env:
 RSCONNECT_USERNAME: ${{ vars.RSCONNECT_USERNAME }}

As an aside, you may notice there’s a stage named “Add Posit Connect deployment config file”. What’s that needed for? When you deploy content the first time using the {rsconnect} package, it will keep a record of some metadata of where it was deployed to in a .dcf file. If you re-deploy the content, {rsconnect} will try to overwrite the existing deployment, by identifying the target by the unique Content ID. Without knowing the Content ID, {rsconnect} has to assume that it’s not safe to overwrite any existing content, and new content must be created instead. Creating a .dcf file and populating it with some details on where the content was previously installed to convinces {rsconnect} that it is safe to overwrite the existing deployment.

Set GitHub Actions variables and secrets

The deploy.yml file requires a number of secrets and variables to be configured in the GitHub Actions workflow. Remember that secrets will be censored in log messages, while variables will be visible.

Head to the “Settings” page for your GitHub repository, and in the side menu go to “Secrets and variables”, then “Actions”.

There are two secrets to set:

• RSCONNECT_CLIENT_ID: The clientId value from your Posit Cloud credentials.
• RSCONNECT_CLIENT_SECRET: The clientSecret from the Posit Cloud credentials.

Followed by four variables:

• RSCONNECT_USERNAME: Your Posit Connect Cloud username, which you set when you created an account. If you have forgotten this, look at the URL once you have logged in to Posit Connect Cloud; The URL will take the format https://connect.posit.cloud/<your-username>.
• APP_TITLE: A display title for your content. This is the title that will appear in your list of deployed content when logged into Posit Connect Cloud.
• APP_NAME: An internal application name. For simplicity, you could set this to your APP_TITLE but with spaces and punctuation replaced with hyphens. For example, “My useful application” becomes “my-useful-application”.
• CONNECT_CONTENT_ID: The Content ID we obtained after the initial deployment to Posit Connect Cloud.

Extending to dev deployments

The example we’ve shown above is designed to deploy when you push to your main or master branch. But if you want a separate deployment for development applications, you can simply extend this action to deploy when changes are pushed to a dev branch, but remember that you’ll need to target a different Content ID, otherwise changes pushed to the dev branch will overwrite your deployment from the main branch.

But keep in mind that “Basic” and “Free” accounts have a limit on the number of applications and a development deployment would count as a separate application to the main deployment. The apps on both these tiers will also be public.

Migration and hosting advice

For users with content already on shinyapps.io, Posit has provided a migration tool to help move your content now, otherwise it will be automatically moved across in early 2027. Links to your content on shinyapps.io will automatically redirect to the new content when done through this tool.

It’s worth using the tool as it allows you to preview and test that the application will work on Posit Connect Cloud. Older applications that use old dependencies or private packages are most likely to encounter issues when migrating to Posit Connect Cloud.

At Jumping Rivers, we often encounter packages and applications that need bringing up-to-date. Our R and Python experts can provide advice and solutions for migrating your content to new hosting solutions that match your needs. Contact us at hello@jumpingrivers.com to see how we can help.

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

This post was originally published on our collaborative substack site. Visit the site to follow us and read more similar posts.

Jointly authored by Chris Brown, Scott Spillias, Carla Sbrocchi and Luis D. Verde Arregoitia.

Chris had putting off trying local models as the setup seemed too complex for him, but Scott convinced him to try it out. Chris did and says it was easier than he had thought. This tutorial walks through the setup.

ollama serve                  # start the server (the desktop app does this for you)
ollama pull qwen2.5-coder     # download a model
ollama run qwen2.5-coder      # chat with a model (downloads it first if you don't have it)
ollama list                   # see which models you've already downloaded
ollama ps                     # see which models are currently loaded in memory
ollama stop qwen2.5-coder     # unload a model from memory now
ollama rm qwen2.5-coder       # delete a downloaded model from disk

In anticipation of the second edition of Data Visualization, which is coming later this year from Princeton University Press, version 2.0.0 of my socviz package is now on CRAN. The update removes some functions that aren’t needed anymore and adds a couple of sample datasets.

TLDR: There’s no official crosswalk between the ONET’s occupational taxonomy and ANZSCO, which is awkward, because Australian researchers use ONET data all the time. This post builds one, explains why you should be suspicious of relying on it and then suggests a better methodology is probably what was applied by the European Commission.

Background

In 2025 I served as an adviser for a project to map occupational transition pathways in India. Although the work leveraged locally-sourced data and occupational profiles, the methodology leaned heavily on studies that use occupational profile data from the Occupational Information Network (O*NET). With the basic idea being that the viability of a transition pathway was in-part determined by how similar two jobs are (conditional on geography, wage rate differentials, education etc).

The Boring Part

I plan to write more about the interesting parts of this project in the future, but I need to start with the boring parts first: developing a crosswalk table between the O*NET-SOC taxonomy and national standards.

I’ll focus on developing a crosswalk / correspondence table between the O*NET-SOC and the Australian and New Zealand Standard Classification of Occupations (ANZSCO). This is both because I’ll use this crosswalk in a future post that uses this standard and as an “official” crosswalk doesn’t exist (despite Australian researchers frequently using the O*NET database).

I suspect one reason official correspondence tables don’t exist already is that the OSCA is a new standard and that the O*NET SOC taxonomy doesn’t cleanly match to the ANZSCO or intermediate correspondence tables. In practice, this means the judgement of the analyst will be required to decide how to match one standard with the other so that it suits their use-case. For instance, if the research is on occupations in the Trucking industry it will be sensible to confirm the data you’re drawing on is being sensibly assigned.

This problem isn’t exclusive to occupational correspondences. The same problems will often rear their head when trying to connect datasets that use different definitions for industries, administrative boundaries and/or products. In most cases the difficulty stems from each standard using a different approach for defining groups, which results in occupations being weirdly assigned at each step of creating a map between standards.

In the case of this post, the problems reared their head at every step from SOC to ANZSCO:

• From ANZSCO to OSCA: In most cases ANZSCO occupations have been assigned to one or more OSCA group, but there are cases where the opposite occurs too, such as Production Manager (Manufacturing), which has been assigned more than one ANZSCO grouping.
• From ISCO-08 (or ESCO) to OSCA: The bridge between the OSCA to ISCO suffers similar problems. However, because ISCO-08 groupings are only provided at the unit group level, in most cases OSCA occupations are bundled into large groups. However, the opposite also occurs too, with Engineering Technologist being assigned to several ISCO-08 groupings at the same time.
• From SOC to ISCO-08 (or ESCO): one of the best crosswalks maps 3,349 occupations to 958 occupations in the O*NET. Once again, these aren’t 1:1 matches,. For example, the ISCO/ESCO job Sports, recreation and cultural centre managers is assigned two jobs from the O*NET. While the O*NET occupation legislators is assigned to more than one distinct ISCO/ESCO occupation category.

The reason I mention this upfront is to make it clear the process is messy. And if I wasn’t intending to replicate research that uses ANZSCO in a future post I wouldn’t bother. But, I am, so I thought I should share my process (and pain) so other people can learn from my mistakes and re-purpose the approach in a way that suits their analysis.

Note: Because the Occupation Standard Classification for Australia (OSCA) is the modern successor of the ANZSCO, this crosswalk (and post) will have a short shelf-life.

Data: the correspondence tables used in this post are available here for the O*NET and here for the ABS. These were originally sourced from the O*NET and ABS on 21/8/2026.

How I used AI in this post: Because developing the correspondence table mainly requires data cleaning and joining occupational definitions, Claude was heavily leaned on to write the code for this post. The write up is more or less untouched by AI.

Project Setup

The code below sets the assumptions for importing correspondence tables and saving outputs.

library(tidyverse)
library(readxl)
library(janitor)

ref_dir_data <- file.path(".", "Data")
ref_dir_out  <- file.path(".", "Outputs")

ref_file_esco_to_soc <- file.path(ref_dir_data, "ESCO_to_ONET-SOC.xlsx")
ref_file_osca_tables <- file.path(ref_dir_data, "OSCA_correspondence_tables_v2.xlsx")

# Keep a SOC link only if this share of the ANZSCO unit group's detailed
# occupations backs it. Set to 0 to keep everything.
ref_min_pct_unit_support <- 50

#import the O*NET / ESCO crosstab
dta_esco_soc_raw <- read_excel(ref_file_esco_to_soc, sheet = 1, skip = 3,
                               col_types = "text") |>
  clean_names() |>
  filter(!is.na(esco_isco_code), !is.na(o_net_soc_2019_code)) |>
  transmute(
    esco_code   = str_trim(esco_isco_code),          # e.g. "8332.5" or "8332"
    esco_name   = str_trim(esco_isco_title),
    isco_code   = str_trim(str_extract(esco_isco_code, "^[^.]+")),
    isco_digits = str_length(isco_code),
    soc_code    = str_trim(o_net_soc_2019_code),
    soc_name    = str_trim(o_net_soc_2019_title),
    link_level  = if_else(str_detect(esco_isco_code, "\\."),
                          "esco_occupation", "isco_unit_group")
  )

# A few rows sit at ISCO MINOR group level (3 digits, e.g. "213"), which is
# coarser than a unit group, so they cannot be placed and are dropped.
dta_esco_soc_unitlevel <- dta_esco_soc_raw |>
  filter(isco_digits == 4)

# ESCO maps some ISCO unit groups directly, and others only via the narrow
# occupations inside them. tHEPrefer the direct mapping; fall back to the
# occupation rows for the 84 groups that have none. 
lkp_isco_with_own_row <- dta_esco_soc_unitlevel |>
  filter(link_level == "isco_unit_group") |>
  pull(isco_code) |>
  unique()

dta_esco_soc_kept <- dta_esco_soc_unitlevel |>
  filter(link_level == "isco_unit_group" | !(isco_code %in% lkp_isco_with_own_row))

# Table 8 is written OSCA -> ISCO; we travel it ISCO -> OSCA. Same pairs.
dta_osca_to_isco <- read_excel(
  ref_file_osca_tables,
  sheet     = "Table 8",
  col_names = c("osca_code", "osca_name", "isco_code", "match_flag", "isco_name"),
  col_types = "text",
  range     = cell_limits(ul = c(6L, 1L), lr = c(NA_integer_, 5L))
) |>
  # The last row is an ABS copyright line, not data.
  filter(!str_detect(coalesce(osca_code, ""), "Commonwealth")) |>
  fill(osca_code, osca_name, .direction = "down") |>
  filter(!is.na(isco_code)) |>
  # "xxxxxx" is the ABS marker for "no counterpart exists".
  filter(osca_code != "xxxxxx", isco_code != "xxxxxx") |>
  mutate(across(c(osca_code, isco_code), str_trim)) |>
  distinct(osca_code, osca_name, isco_code, isco_name)

# Table 1 is written ANZSCO -> OSCA.
dta_anzsco_to_osca <- read_excel(
  ref_file_osca_tables,
  sheet     = "Table 1",
  col_names = c("anzsco_code", "anzsco_name", "osca_code", "match_flag", "osca_name"),
  col_types = "text",
  range     = cell_limits(ul = c(6L, 1L), lr = c(NA_integer_, 5L))
) |>
  filter(!str_detect(coalesce(anzsco_code, ""), "Commonwealth")) |>
  fill(anzsco_code, anzsco_name, .direction = "down") |>
  filter(!is.na(osca_code)) |>
  filter(anzsco_code != "xxxxxx", osca_code != "xxxxxx") |>
  mutate(across(c(anzsco_code, osca_code), str_trim)) |>
  distinct(anzsco_code, anzsco_name, osca_code) |>
  # ANZSCO codes are 6 digits (a detailed occupation). The first 4 are the unit
  # group, which is the level the O*NET analysis reports at.
  mutate(anzsco_unit_code = str_sub(anzsco_code, 1, 4))

Joins

The code below joins each correspondence pair sequentially. Because each mapping splits and merges occupational classifications differently, the unified crosswalk results isn’t a clean 1:1 correspondence. For this reason a better approach is likely to be matching occupational descriptions from either standard, such as was done by the European Commission for mapping the O*NET to ISCO-08, but I’ve already written the code so here we are…

# Start at the ESCO level of detail: one row per ESCO occupation and SOC code.
dta_esco_x_soc <- dta_esco_soc_kept

# COLLAPSE to the ISCO unit group. This is where the ESCO occupation codes and
# names get dropped -- they cannot be carried further, because the ABS tables
# are keyed on the ISCO unit group and not on ESCO.
dta_isco_x_soc <- dta_esco_x_soc |>
  distinct(isco_code, soc_code, soc_name)

dta_isco_x_soc_x_osca <- dta_isco_x_soc |>
  inner_join(dta_osca_to_isco, by = join_by(isco_code),
             relationship = "many-to-many")

dta_isco_x_soc_x_osca_x_anzsco <- dta_isco_x_soc_x_osca |>
  inner_join(dta_anzsco_to_osca, by = join_by(osca_code),
             relationship = "many-to-many")

dta_crosswalk_all_levels <- dta_isco_x_soc_x_osca_x_anzsco |>
  select(soc_code, soc_name, isco_code, isco_name, osca_code, osca_name,
         anzsco_code, anzsco_name, anzsco_unit_code) |>
  arrange(soc_code, isco_code, anzsco_code)

Collapsing Occupations to the Unit-Group Level

Because the OSCA occupations are only mapped to the ISCO-08 unit-group (the first 4 digits of the code), the code below collapses the correspondence table to provide a listing of major SOC and ESCO occupations by unit group. As you’d expect, this results in a lot of granularity being lost.

lkp_anzsco_unit_detail <- dta_crosswalk_all_levels |>
  distinct(anzsco_unit_code, anzsco_code, anzsco_name) |>
  summarise(
    nmb_unit_occs = n_distinct(anzsco_code),
    # These files carry no ANZSCO unit group titles, only occupation titles, so
    # the lowest-numbered occupation stands in as the label.
    anzsco_unit_name = anzsco_name[order(anzsco_code)][1],
    .by = anzsco_unit_code
  )

# Advisory sanity check: do the ISCO major group (1st digit) and the SOC major
# group (1st 2 digits) sit in compatible broad families? Some ESCO mappings are
# simply poor, and this catches them. It has false positives, so it is reported
# as a column and never filtered on.
lkp_valid_major_group_pairs <- tribble(
  ~isco_major, ~soc_majors,
  "0",         "55,33",                              # armed forces
  "1",         "11",                                 # managers
  "2",         "13,15,17,19,21,23,25,27,29",         # professionals
  "3",         "13,15,17,19,21,25,29,31,33,49",      # technicians
  "4",         "41,43",                              # clerical
  "5",         "31,33,35,37,39,41",                  # service and sales
  "6",         "45",                                 # agriculture
  "7",         "47,49,51",                           # trades
  "8",         "51,53",                              # plant and machine
  "9",         "35,37,41,45,47,53"                   # elementary
) |>
  separate_longer_delim(soc_majors, delim = ",") |>
  rename(soc_major = soc_majors) |>
  mutate(broad_group_match = TRUE)

rlt_crosswalk_by_unit_group <- dta_crosswalk_all_levels |>
  summarise(nmb_occ_support = n_distinct(anzsco_code),
            .by = c(soc_code, soc_name, isco_code, anzsco_unit_code)) |>
  left_join(lkp_anzsco_unit_detail, by = join_by(anzsco_unit_code)) |>
  mutate(pct_unit_support = round(100 * nmb_occ_support / nmb_unit_occs, 1),
         isco_major = str_sub(isco_code, 1, 1),
         soc_major  = str_sub(soc_code, 1, 2)) |>
  left_join(lkp_valid_major_group_pairs, by = join_by(isco_major, soc_major)) |>
  mutate(broad_group_match = coalesce(broad_group_match, FALSE)) |>
  select(-isco_major, -soc_major) |>
  arrange(soc_code, anzsco_unit_code)

Filtering out Poor Matches

The final step drops matches that are supported by a minority of occupations after joins. Each ANZSCO group holds several occupations and the joins match each of them to a SOC group separately. So, when more ANZSCO occupations are matched to the same SOC code within a unit group it’s assumed the match is stronger, while weaker matches are dropped and assumed to reflect the many quirks of trying to match definitions like this.

rlt_crosswalk_for_onet <- rlt_crosswalk_by_unit_group |>
  filter(pct_unit_support >= ref_min_pct_unit_support) |>
  summarise(isco_codes        = paste(sort(unique(isco_code)), collapse = "; "),
            pct_unit_support  = max(pct_unit_support),
            broad_group_match = any(broad_group_match),
            .by = c(anzsco_unit_code, anzsco_unit_name, soc_code, soc_name)) |>
  select(anzsco_code  = anzsco_unit_code,
         label_4digit = anzsco_unit_name,
         soc          = soc_code,
         soc_label    = soc_name,
         isco_codes, pct_unit_support, broad_group_match) |>
  arrange(anzsco_code, soc)

stopifnot(
  "Duplicate anzsco_code x soc rows would double-weight a SOC in the O*NET average" =
    nrow(rlt_crosswalk_for_onet) ==
      nrow(distinct(rlt_crosswalk_for_onet, anzsco_code, soc))
)

Exploratory analysis

Claude produced a lot of exploratory analysis and checks after I harangued it and quizzed its analysis, but I think the summary below is the most useful. It essentially checks how many occupations from the source correspondence files remained in the unit-group mapping.

The TLDR: most of the O*NET occupations made it, despite the inconsistencies encountered along the way. On one level that’s a surprisingly good outcome, particularly given how naive the final mapping approach is. But it hides a few things that I’d worry about if using the table in my analysis: occupations being more likely to drop out in some categories more than others and SOC occupations landing in the wrong groups (which could result in drawing on incorrect data from the O*NET).

# The links behind the delivered table, taken at the level where isco_code is
# still a single column. This is the same set of links as rlt_crosswalk_for_onet.
dta_final_links <- rlt_crosswalk_by_unit_group |>
  filter(pct_unit_support >= ref_min_pct_unit_support)

# Codes are lost at two different points, so both are reported. The joins are
# inner joins, so a code with no counterpart drops out silently. The support
# filter then removes more. OSCA is an intermediate hop and is not carried into
# the final table at all, so it has no final count.
rlt_codes_kept_vs_dropped <- tibble(
  classification = c("O*NET-SOC", "ISCO-08", "OSCA", "ANZSCO (unit group)"),
  nmb_in_source = c(
    n_distinct(dta_esco_soc_raw$soc_code),
    n_distinct(c(dta_esco_soc_unitlevel$isco_code, dta_osca_to_isco$isco_code)),
    n_distinct(c(dta_osca_to_isco$osca_code, dta_anzsco_to_osca$osca_code)),
    n_distinct(dta_anzsco_to_osca$anzsco_unit_code)),
  nmb_after_joins = c(
    n_distinct(dta_crosswalk_all_levels$soc_code),
    n_distinct(dta_crosswalk_all_levels$isco_code),
    n_distinct(dta_crosswalk_all_levels$osca_code),
    n_distinct(dta_crosswalk_all_levels$anzsco_unit_code)),
  nmb_in_final = c(
    n_distinct(dta_final_links$soc_code),
    n_distinct(dta_final_links$isco_code),
    NA_integer_,
    n_distinct(dta_final_links$anzsco_unit_code))
) |>
  mutate(nmb_dropped = nmb_in_source - nmb_in_final,
         pct_kept    = round(100 * nmb_in_final / nmb_in_source, 1))

rlt_codes_kept_vs_dropped

Output files

The final code block outputs the crosswalks.

ref_stamp <- format(Sys.Date(), "%y%m%d")
if (!dir.exists(ref_dir_out)) dir.create(ref_dir_out)

write_csv(dta_crosswalk_all_levels,
          file.path(ref_dir_out, paste0(ref_stamp, " - crosswalk_full.csv")))
write_csv(rlt_crosswalk_by_unit_group,
          file.path(ref_dir_out, paste0(ref_stamp, " - crosswalk_unit.csv")))
write_csv(rlt_crosswalk_for_onet,
          file.path(ref_dir_out, paste0(ref_stamp, " - crosswalk_for_onet.csv")))

Summing up

Because the main point of this post is to produce a correspondence table for a future post, I’m not going to split hairs about the assignment. However, if you’re intending to use the crosswalks for something more rigorous, I hope this post demonstrates clearly why you should.

A few things I’d check:

• Confirm assignments for occupations your analysis relies on. For me, it’s truck drivers as the next post replicates analysis that looks at occupational paths for that group.
• Validate assignments / groupings. Particularly for above, but the approach I’ve used is naive insofar as I’ve just joined each table together and used a simple filtering mechanism for dropping obvious cases. However, having done this makes me even more convinced that a simple comparison of job titles or descriptions would be easier to justify (like the European Commission’s methodology).
• What was dropped and is it valuable to your analysis.

Which is perhaps something I’ll do in a future post, but for now I’d just say that if you have a better approach, or want to point out the problems with mine, please get in touch. I wrote this up in an afternoon: this isn’t mean to be a master work. So, I’ll happily link to better approaches and make corrections to avoid leading others astray.

The post The boring part first: Building a crosswalk from O*NET-SOC to ANZSCO appeared first on Giles.

This week I took time to re-read The Programmer’s Brain by Felienne Hermans. Among the many gems one idea that stuck with me is that not learning how to do something and looking it up every time will make you less efficient. Therefore I starting feeling bad about one particular thing I never get done on my own: extracting parts of a string, even in simple cases.

Avoidance tactics

For instance, how do you extract names from the following templated sentences?

sentences <- c(
  "My name is Moomin.",
  "My name is Little My.",
  "My name is Snork Maiden."
)

I would use either one of these two tactics:

1. Removing the final period, and the first words.

sub("My name is ", "", sub(".$", "", sentences))
#> [1] "Moomin"       "Little My"    "Snork Maiden"

1. Looking up the regex syntax online (maybe even dutifully reading the stringr cheatsheet) or via a LLM. This could get me code calling stringr:

stringr::str_extract(sentences, "My name is (.*).", group = 1)
#> [1] "Moomin"       "Little My"    "Snork Maiden"
# Look arounds
stringr::str_extract(sentences, pattern = '(?<=My name is ).+(?=\\.)')
#> [1] "Moomin"       "Little My"    "Snork Maiden"

Or some base R code:

regmatches(
  sentences,
  regexpr("(?<=My name is ).+(?=\\.)", sentences, perl = TRUE)
)
#> [1] "Moomin"       "Little My"    "Snork Maiden"

My problems

Really I had two problems preventing me from being really autonomous:

• Not knowing enough regex.
• Not knowing where to put the regex, for whatever reason I felt I had to choose between adding a dependency on stringr or using the complicated two-step regexpr/regmatches syntax.

Solutions

To solve the first problem (lack of regex knowledge), I need to be more intentional about remembering the look-arounds syntax for instance, or what a group is.

What solved my second problem (thinking I had to choose between a dependency or code distateful to me) was a very simple tip by my rOpenSci colleague Jeroen Ooms: using sub()! The code below replaces the sentences with the names (capture groups) in them.

sub("My name is (.*).", "\\1", sentences)
#> [1] "Moomin"       "Little My"    "Snork Maiden"

This is code he seems to use quite often¹.

What was especially great about this tip, beside its timing when I was reading the book, is that it made me “get” groups more easily. The group is what’s between parentheses, it’s not more complicated than that (at least I don’t need to know more right now).

Conclusion

I will try to keep mindful of not being too lazy to learn some things when I can actually learn them. And now I know that I can use something else than stringr or the not so easy base R syntax with regmatches(): a simple call to sub()! Watch me win seconds every time I have to extract parts of a string.

biocjobs: "1.0"
name: deseq2-differential-expression
package: DESeq2
title: DESeq2 differential expression
description: >
  Runs the canonical DESeq2 workflow on a raw count matrix: size factor and
  dispersion estimation, negative-binomial GLM fitting and a Wald test for
  one pairwise contrast.
version: "1.0.0"
script: scripts/deseq2-differential-expression.R
depends: [apeglm, ashr]

inputs:
  - name: counts
    format: tsv
    label: Raw count matrix
    help: >
      Tab-separated matrix of raw (un-normalized) integer read counts.

outputs:
  - name: results
    format: tsv
    label: Differential expression results

options:
  - name: shrinkage
    type: choice
    choices: [apeglm, ashr, normal, none]
    default: apeglm
    label: Log2 fold change shrinkage
  - name: alpha
    type: float
    default: 0.1
    min: 0
    max: 1
    label: FDR threshold

resources:
  cpus: 1
  memory_gb: 4

citations:
  - doi: 10.1186/s13059-014-0550-8

params <- BiocJobs::jobParams("DESeq2", "deseq2-differential-expression")

Rscript -e 'BiocJobs::execJob("DESeq2", "deseq2-differential-expression")' \
    --counts counts.tsv --coldata coldata.tsv \
    --contrast_factor condition --contrast_numerator treated \
    --contrast_denominator control --alpha 0.05

© 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

Prompted () by some AI dread, I decided to go back to some basics and re-read The Programmer’s Brain by Felienne Hermans. Felienne Hermans’ work caught my attention when she gave a keynote talk at a Posit conference years ago. The book was a highlight of my week: extremely interesting, and easy to follow. Here are some notes, thanks to tiny bookmarks I added as a I read.

The main characters

The main characters in the book are

• the long-term memory (knowledge);
• the short-term memory (information right now);
• the working memory (processing power).

Everything is brought back to them.

Making an effort pays off

I am fascinated by the fact that a schoolteacher called Ballard found out that “when you actively try to recall information without additional study, you will remember more of what you learned”.

You can also strengthen your memories by actively thinking, elaborating around something.

Cognitive loads

Felienne Hermans summarizes the different types of cognitive load:

Intrinsic load: how complex the problem is in itself. Extraneous load: what outside distractions add to the problem. Germane load: cognitive load created by having to store your thought to long-term memory.

Regarding extraneous load, an example that’s given is a poorly formulated math problem.

Cognitive refactoring

I remembered this idea from my first read: you can refactor code to understand it better, a refactor you do only for yourself.

One example that’s given is replacing unfamiliar language constructs such as anonymous functions. It made me think of my overcomplicating a PR by both changing something crucial and replacing for loops with reduce, that were unfamiliar to my collaborator. I should have split the two changes in two PRs.

A related quote from the book:

‘“readable” is really in the eye of the beholder’

Help your working memory

When mentioning strategies for helping your working memory, such as creating state tables or diagrams, the author mentioned PythonTutor by Philip Guo, which visualizes the execution of a program. It reminded me of the boomer R package by my cynkra colleague Antoine Fabri, that lets you inspect the intermediate steps of a call.

subset(head(penguins, 2), bill_len > 47) |> boomer::boom()
#> 💣 subset(head(penguins, 2), bill_len > 47) 
#> · 💣 💥 head(penguins, 2) 
#> ·   species    island bill_len bill_dep flipper_len body_mass    sex year
#> · 1  Adelie Torgersen     39.1     18.7         181      3750   male 2007
#> · 2  Adelie Torgersen     39.5     17.4         186      3800 female 2007
#> · 
#> · 💣 💥 bill_len > 47 
#> · [1] FALSE FALSE
#> · 
#> 💥 subset(head(penguins, 2), bill_len > 47) 
#> [1] species     island      bill_len    bill_dep    flipper_len body_mass   sex         year       
#> <0 rows> (or 0-length row.names)
#> 
#> [1] species     island      bill_len    bill_dep    flipper_len body_mass   sex         year       
#> <0 rows> (or 0-length row.names)

Coupling that with constructive, another package of Antoine’s, might help one represent code better.

Roles of variables

The book has a list (by Jorma Sajaniemi) of the eleven roles a variable can have, e.g. “fixed value” or “stepper” (i in a for loop). Interesting vocabulary! The book even features a flowchart to help us determine a role a variable plays.

Parallels with natural languages

The author explains a technique for understanding code by circling all variables, linking them, etc. It reminds me of how I’d handle Latin text I had to translate in high school. I had a color and shape code, it looked very pretty and worked well.

Speaking of languages, the book draws some parallels between computer and natural languages. In particular, it explains how text comprehension strategies (like questioning or summarizing) apply to code reading.

Keeping notes

I will try to do better at taking notes on a piece of paper when I work. I already do in some cases, for instance when reviewing packages for rOpenSci. But the book really insists how it can support your memory, or resume work after an interruptions.

Beside those throwaway notes, it’s important to document/comment code to prevent future contributors to fall in some traps and to facilitate onboarding of new contributors. Recent example.

Further programming languages

IDEs:

“transfer between two programming languages is more likely if you program two different languages in the same IDE, which is a strong argument for using one IDE for multiple languages.”

Language choice:

“if you set out to learn a new language to expand your way of thinking, it’s important to pick one language that’s fundamentally different from the ones you’ve already mastered.”

The book also explains how some knowledge you have in one language means you might have to “unlearn” some syntaxes. It reminded of faux amis (false friends) for French-speaking learners of English, like “actually” that doesn’t mean actuellement (currently).

Names are important…

And the book explains why, gives useful tips. There’s a whole chapter on the topic.

I liked one of the conventions by Butler: “Identifiers should consist of words and only use abbreviations when they are more commonly used than the full words”. Recently I was very stubborn about not using “comb” for “combination” in igraph. I also enjoyed another conventions from that same list: “Identifier should not combine uppercase and lowercase character in non standard ways”, with the example Page_counter. I might be selectively reading the rules that I like.

The chapter conveys the perspective by Allamanis that names should be consistent across a codebase, because that helps chunking (when you parse code into meaningful bits).

The author underlines that you should evaluate the quality of names after coding, not during code, as it might be too much cognitive load. It made me think of Git commits: you can improve them after coding, when you’re coding you might not be able to create a perfect Git history.

Another tidbit that I found interesting is that when you improve names in your codebase, the places where you find bad names might be the places with hidden bugs for various reasons (like correlation between bad names and mistakes by a novice programmer or a programmer confused by the complexity of the problem at hand).

Automatization

Some things you know so well that you can do them without thinking much, which makes you more efficient. An argument for learning and deliberate practice.

Reading about code

Sometimes if you’re writing very complex code, you don’t learn much, “your brain was so engaged it could not store the solutions”.

Therefore, worked examples can help you learn: collaborating with someone, reading code on GitHub, reading books or blog posts about code. You don’t only (and necessarily) learn by doing.

Curse of expertise

I enjoyed reading again about the “curse of expertise”, that is especially relevant when teaching:

“Once you have mastered a skill sufficiently, you will inevitably forget how hard it was to learn that skill or knowledge.”

Conclusion

I would highly recommend reading The Programmer’s Brain by Felienne Hermans! Maybe even more than once like I did since I had clearly not committed everything to long-term memory.

The epilogue mentions some further reading including A Philosophy of Software Design by John Ousterhout which solved a mystery for me: that is where I had heard of that book!

Below are the slides for my Futurize – Tearing Down Parallelization Barriers in R with Transpilers talk that I presented at the useR! 2026 conference in Warzaw, Poland.

Title: Futurize – Tearing Down Parallelization Barriers in R with Transpilers
Speaker: Henrik Bengtsson
Slides: HTML (16 slides; 18 minutes)
Video: To appear

The new futurize package makes it easier than ever before to parallelize existing map-reduce calls – just pipe the call to futurize() and you’re done!

ys <- lapply(xs, fit_model) |> futurize()
ys <- map(xs, fit_model) |> futurize()
ys <- foreach(x = xs) %do% fit_model(x) |> futurize()
ys <- llply(xs, fit_model) |> futurize()

It also works with other popular domain-specific calls, e.g.

xs_smooth <- stats::kernapply(xs, k = k) |> futurize()
b <- boot(city, ratio, R = 999) |> futurize()
model <- caret::train(Species ~ ., data = iris, method = "rf", trControl = ctrl) |> futurize()
cv <- glmnet::cv.glmnet(x, y) |> futurize()
m <- lme4::allFit(models) |> futurize()

See the futurize package site for more examples and details.

I want to thank the useR! organizers, staff, volunteers, sponsors, and everyone else who contributed to this amazing event making it possible for the R community to come together in person. Just like last year’s useR! 2025 in the US, it was fantastic to see so many first and second timers attending the useR! conference in Europe. It’s very refreshing and it clear that we are on a great track to recover from not having in-person R conferences during COVID-19 pandemic. Next year’s useR! will take place in Santiago, Chile in July 2027 – exciting!

/Henrik

Links

• useR! 2026: https://user2026.r-project.org/
• futureverse website: https://www.futureverse.org/
• futurize package CRAN, GitHub, pkgdown

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

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

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

What Data Do You Need for a Sports Betting Model?

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

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

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

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

Useful R Packages for Sports Data

The SportsDataverse ecosystem provides packages for several major sports:

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

Install the core packages with:

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

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

Getting a Sports Odds API Key

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

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

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

Add the following line to the .Renviron file:

ODDS_API_KEY=YOUR_PRIVATE_API_KEY

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

Sys.getenv("ODDS_API_KEY")

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

Download Current Sports Betting Odds in R

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

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

api_key <- Sys.getenv("ODDS_API_KEY")

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

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

odds_raw <- resp_body_json(response, simplifyVector = FALSE)

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

Convert the API Response into Tidy Data

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

tidy_odds <- function(events) {

  map_dfr(events, function(event) {

    map_dfr(event$bookmakers, function(bookmaker) {

      map_dfr(bookmaker$markets, function(market) {

        map_dfr(market$outcomes, function(outcome) {

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

odds_df <- tidy_odds(odds_raw)

glimpse(odds_df)

The resulting table can contain columns such as:

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

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

library(lubridate)

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

Understanding Moneylines, Spreads and Totals

The API uses different market identifiers:

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

You can filter the dataset to analyze a single market:

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

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

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

Convert Decimal Odds into Implied Probabilities

Decimal odds can be converted into raw implied probability using:

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

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

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

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

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

How to Collect Historical Betting Odds

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

Historical data should ideally include:

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

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

When saving a current snapshot, include the collection time:

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

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

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

readr::write_csv(odds_snapshot, file_name)

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

Combine Betting Odds with Sports Performance Data

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

install.packages("nflreadr")

library(nflreadr)
library(dplyr)

pbp <- load_pbp(2025)

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

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

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

Sports Data Sources for NFL, NBA, MLB and NHL

NFL Data

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

NBA Data

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

MLB Data

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

NHL Data

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

Build a Simple Probability Model

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

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

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

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

Identify Potential Value Bets

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

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

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

Backtest the Model by Season

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

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

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

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

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

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

Using Bayesian Models for Sports Prediction

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

A Bayesian workflow can:

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

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

Automate Daily Sports Data Updates

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

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

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

Common Sports Betting Backtesting Mistakes

Using Information That Was Not Available Before the Game

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

Ignoring Changes in the Betting Line

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

Testing Too Many Strategies

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

Using Accuracy as the Only Metric

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

Assuming a Small Positive Return Proves an Edge

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

From Raw Data to a Complete Betting System

A complete sports betting workflow can be summarized as:

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

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

Frequently Asked Questions

Can I get sports betting data for free in R?

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

What is the best R package for sports betting odds?

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

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

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

How many years of data do I need?

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

Can a sports betting model guarantee profits?

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

Conclusion

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

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

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

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

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

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

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

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

Let’s talk about the weather

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

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

library(dplyr)

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

library(ggplot2)

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

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

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

A little summer upskilling

This is also where R becomes particularly useful.

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

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

There probably is.

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

That’s exactly where structured training can help.

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

Make the summer count

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

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

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

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

No complicated spreadsheets required.

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

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

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

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

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

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

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

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

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

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

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

II. Using R on an Android Device

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

III. Using R on an iPhone

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

IV. Sample Program: Raspberry Pi

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

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

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

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

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

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

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

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

Architectural Comparison: Modeltime/GluonTS vs. Native R Torch

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

Deep Dive into the Code Architecture

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

1. Bounded Student-t Distribution Head

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

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

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

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

2. Variance-Controlled Stochastic Monte Carlo Sampling

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

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

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

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

3. Granular Interactive Plotly Visualization

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

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

Complete R Script

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

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

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

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

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

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

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

seqs <- create_sequences(train_data$close_scaled, lookback_length)

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

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

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

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

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

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

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

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

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

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

print(metrics_summary)

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


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

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

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

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

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

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

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

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

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

Conclusion

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

Yesterday I was in Berlin for the day and presented at the local PyData Meetup on our new solver for fixed effects regression. It was of course a lot of fun! The development was prompted by an issue we received via github a couple of months ago:

Hours vs minutes! Urgh. Luckily, we’ve made some progress in between, so the issue could be closed with with

If you are curious how the new solver works, you can find the slides for the talk here, and please make sure to take a look at the code or to even run a regression with it! We also have a longer write up up in preparation, but as unfortunately is too often the case, the Pareto principle bites and we’ve been refining the last 20% for quite a while now.

Yesterday I was in Berlin for the day and presented at the local PyData Meetup on our new solver for fixed effects regression. It was of course a lot of fun! The development was prompted by an issue we received via github a couple of months ago:

Hours vs minutes! Urgh. Luckily, we’ve made some progress in between, so the issue could be closed with with

If you are curious how the new solver works, you can find the slides for the talk here, and please make sure to take a look at the code or to even run a regression with it! We also have a longer write up up in preparation, but as unfortunately is too often the case, the Pareto principle bites and we’ve been refining the last 20% for quite a while now.

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

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

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

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

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

Introduction

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

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

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

Example self-contained R script using ir

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

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

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

This script can be run with

ir run render.R

Example self-contained R script using uvr

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

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

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

This script can be run with

uvr run render.R

Automation with just in a complex directory structure

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

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

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

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

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

Bonus 1 – Example self-contained Quarto document using ir

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

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

Rest of Quarto document follows ...

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

ir render tutorial.qmd

More details are given in the ir Quarto docs.

Bonus 2 – Making the R script executable

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

For ir we add

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

and for uvr we add

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

We then make the script executable

chmod +x render

and run it with

./render

Summary

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

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

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

Installation

Old version, depends on Shiny:

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

New version, does not use Shiny:

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

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

Quick Start

Single-script app

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

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

Modular R package app

Create an R package with modular components:

library(tabler)

pkg_template("mydashboard")

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

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

Loading/Progress bar

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

Login page

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

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

Does it run Doom?

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

Available Layouts

There are additional examples for each of the following layouts:

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

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

Differences with Shiny

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

License

Apache License (>= 2)