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)

Usually I don’t watch much youtube (a social phenomenon I have somewhat missed out on), but here’s an interesting video about another social phenomenon I have somewhat experienced myself .

With LLM coding, I feel like I can do anything and start working on too many things simultaneously / in multiple parallel sessions (pyfixest features and refactoring, a paper on our new demeaning algo, R bindings for the new algo, presentations on this and that, blog posts, a JOSS paper draft, etc etc).

I am paying for subscriptions with Codex, Claude, and OpenCode, and when I don’t spend all of my usage, I feel like I missed an opportunity, so sometimes I skip going to the gym to keep on working because I still have usage to spend or usage will reset in 30 minutes so I can keep on working prompting my agents to keep on working… This was particularly bad the weekend that Fable was temporarily available on the 20$ plan – I almost couldn’t leave my keyboard, the desire to keep on prompting was too strong. Needless to say that – because I prompted without much thinking – the results I got were not very good.

I usually start working on OSS after work around 21:00, and suddenly it’s 23:30 and I want to stop, but there is one last prompt to write and finally I stop and it is 00:30 but I feel no accomplishment because I haven’t been doing any deep work … but at least, I used all of my usage (#tokenmaxxing)!

Long story short, I think a more optimal strategy for me going forward is to reduce the number of subscriptions (fewer token limits to hit) and topics I work on in parallel and, to some degree, stop over-relying on an auto-mode-and-review workflow and switch back to a more “socratic” style of coding with LLMs.

When dockerizing an R Shiny application hosted via Shiny Server (built from shiny server docker image https://hub.docker.com/r/rocker/shiny), a common issue developers face is that environment variables set via ENV instructions in the Dockerfile (or passed at runtime via docker run -e) are completely missing inside the R Shiny app.

For example, if you set the following environment variable in your Dockerfile:

1

ENV SHOULD_IN_SHINY="From dockerfile"

When your app launches and runs Sys.getenv("SHOULD_IN_SHINY"), it returns an empty string "" instead of "From dockerfile". This is not a problem with Docker itself, but rather a consequence of how Shiny Server spawns R worker processes.

In this post, we will look into the root cause behind this behavior in Shiny Server and demonstrate the two recommended solutions to correctly expose environment variables to your R Shiny workers.

The Root Cause: How Shiny Server Spawns R Processes

The reason environment variables do not carry over to your R session lies in how Shiny Server executes R worker processes inside the container.

Shiny Server runs as a system service (typically as root or shiny). When launching an app instance, it re-executes R as the shiny unprivileged user using su. Specifically, the execution call combines two mutually exclusive flags:

1

su shiny --login --preserve-environment -c "... R ..."

Let’s break down what these flags request from su:

1. --login (-l): Starts a login shell. This intentionally resets the environment to a minimal whitelist (HOME, PATH, USER, TERM, etc.) and sources system startup profiles (/etc/profile, ~/.profile).
2. --preserve-environment (-p): Explicitly asks su to keep the current environment inherited from the caller (which includes Docker’s ENV variables).

Because Linux’s su utility cannot honor both conflicting behaviors, it chooses --login and prints a warning:

1

su: ignoring --preserve-environment, it's mutually exclusive with --login

Note: This warning is emitted by Shiny Server’s underlying process execution call, not by the R app itself. Attempting to suppress or patch it requires modifying and rebuilding Shiny Server C++ code, which is rarely practical.

Because the login shell wins, all custom environment variables passed to the container via ENV or docker run are wiped before R ever starts.

Solutions

Since patching Shiny Server is unnecessary, we can utilize the natural extension points provided by the login shell or R itself.

Solution 1: Use ~/.profile (Recommended for User-Level Shell Vars)

Since --login causes the shell to source ~/.profile for the shiny user, we can write our environment variables to /home/shiny/.profile during the Docker build phase.

In your Dockerfile:

1
2

# Append environment variable to the shiny user's profile
RUN echo 'export SHOULD_IN_SHINY="From dockerfile"' >> /home/shiny/.profile

When Shiny Server executes su shiny --login ..., the login shell will read /home/shiny/.profile and load SHOULD_IN_SHINY into the environment right before starting the R process.

Solution 2: Use ~/.Renviron (Recommended for R-Specific Configs)

Alternatively, R automatically inspects and loads ~/.Renviron on startup, right after the shell environment is initialized.

~/.Renviron is purpose-built for R:

• Scoped strictly to R processes.
• Does not use shell export keywords or $ variable expansions.
• Uses simple KEY=value key-value pairs.

In your Dockerfile, set up /home/shiny/.Renviron:

1
2
3

# Set up .Renviron for the shiny user
RUN echo 'SHOULD_IN_SHINY="From dockerfile"' >> /home/shiny/.Renviron \
    && chown shiny:shiny /home/shiny/.Renviron

In your Shiny app (app.R or server.R), you can access it reliably:

1
2

should_in_shiny <- Sys.getenv("SHOULD_IN_SHINY")
# Returns: "From dockerfile"

Summary

• The Problem: Shiny Server calls su shiny --login --preserve-environment, causing su to discard inherited environment variables (such as Docker ENV) in favor of a clean login shell.
• Solution 1 (~/.profile): Append export KEY="value" to /home/shiny/.profile in your Dockerfile.
• Solution 2 (~/.Renviron): Append KEY="value" to /home/shiny/.Renviron in your Dockerfile.

Both approaches integrate seamlessly with Docker builds and ensure your Shiny application receives all required configuration variables cleanly.

As for last year, I’ve got a poster accepted at COPA conference.

This year, the poster is about ‘Zero-Shot Probabilistic Stock Returns Forecasting with Pretrained RVFL Networks’.

Links to the paper and the code repository are below:

• Paper: https://www.researchgate.net/publication/412115839_Zero-Shot_Probabilistic_Stock_Returns_Forecasting_with_Pretrained_RVFL_Networks
• R Code: https://github.com/thierrymoudiki/2026_05_28_Pretrain_Ridge2_Stocks_Full_Pipeline

A few days ago I had a coffee in Paris with my long-time friend (and former Statistics & Computing editor) Gilles Celeux, and he mentioned me stopping solving and posting maths puzzles like those weekly published by Le Monde. They have indeed vanished with the retirement of the authors, but Gilles added that the arrival of LLMs would have made the exercise moot. I disagreed as (i) the fun of solving the puzzle on my own  has not gone away and (ii) the pedagogical appeal of the puzzle and its resolution remains. As the next Fiddler puzzle arrived in my mailbox, my resolution was put to the test (contrariwise to the previous entry, which did not require massive computations):

The Fiddler League consists of two teams. Over a season, they play each other 162 times. Each team has an equal chance of winning each game, and the results of games are independent. Over the season, on average, how many games would you expect the team with the better record to have won?

I started on the wrong foot with E[X|X≥81] when X is Bin(162,½), equal to 85.77 (either directly or with a Normal approximation), which differs from my second thought, E[max(X,162-X)]=86.07 (either directly or with a Normal approximation), which is larger because of the reflection produced by max. While the first computation was manageable, the second one seemed to involve simulation and I caved in prompting Claude, which provided the answer along with the connection

E[max(X,162−X)]=E[X∣X≥81](1+p₈₁​)−81p₈₁​

After some expansion, the League boasts 30 teams. Over a season, each team plays each other team five times. (Each team plays a total of 145 games.) Again, each team has an equal chance of winning each game, and the results of games are independent. Over the season, on average, how many games would you expect the team with the best record to have won?

The best record is max(X_i) with each of the 30 X_i‘s a sum of 29 Y_ij and the Y_ij=5-Y_ji distributed as Bin(5,½). The X_i‘s are thus Bin (145,½) but dependent. While I could not figure out a closed form answer for the expectation, a direct Monte Carlo resolution is obviously feasible, with Claude (rather than me) running it over 400 million repetitions, but a 30 dimensional Normal approximation exploiting the correlation of 1/29 between the components leads to roughly 85 as the expected value. (Again computed by a Claudicant simulation.)

While the conclusion that the Normal approximation is pretty accurate with so many terms in the Binomial variates is quite unsurprising, Claude saves me coding time without ruining the puzzle altogether. (And Gemini made me aware that the name of the café where Gilles and I regularly meet, L’Écir, is an Auvergne noun for a local, dangerous, mountain blizzard! Thus linking the place to the foundation of the café by Auvergne expatriates…)

Last month, Vicki Boykis recommended the book “Code that Fits in Your Head” by Mark Seemann. I was intrigued, especially by her takeaway that “Writing good software should be a slow and deliberate craft”.

The book reminded of The Pragmatic Programmer in the breadth of topics it covers. However, I much preferred its tone.

Sustainability

Throughout the book, the author makes it clear that he believes one should pay attention to the code one writes, to its internal quality. The following two sentences even got their little frame:

“The goal is not to write code fast. The goal is sustainable software.”

Arrange Act Assert

What useful phrase!

You know how in a test you can have something like:

test_that("Exception Bla handled", {
  withr::local_option("bla" = 1)
  test_thing <- 23

  x <- my_function(test_thing)

  expect_equal(x, 2)
})

The first two lines within the test_that() call prepare what’s needed. That’s the arrange phrase.

The line calling my_function() does the thing I want to test. That’s the act phrase.

The line calling expect_equal() checks that thing. That’s the assert phrase.

Arrange, act, assert! The book even mentions you can separate the three phases with empty lines.

Red Green Refactor

Another neat phrase! It describes how you could iteratively work on a part of your software.

• Red: you write a failing test for e.g. a feature.
• Green: you write the code to make the test pass.
• Refactor: you improve that code (but don’t make the test fail again!).

I appreciated how the author wrote that even the red phase isn’t easy, because you can get a passing test that you thought would fail. I think the point is made again elsewhere in the book: you can imagine a test will pass or fail under certain conditions, but better to actually prove it by running the test!

The Devil’s Advocate

A technique presented in the book is The Devil’s Advocate, in which you write wrong code on purpose, to see whether your tests will detect it. I suppose it’s a less random (but more difficult?) version of mutation testing, in which you evaluate how often your test suite detects mutants of your code. In R, mutation testing is provided by e.g. the mutator package.

Code reviews that do not block

There’s a short discussion of when to do code reviews (regularly) so as to not block your team mates. I enjoyed seeing this, as it’s important to not always be a rate-limiting factor. It shows the book’s advice is pragmatic.

Bundle smaller breaking changes

In the chapter “Augmenting code” that’s about working with existing code, the author – among other things – discusses whether to bundle or separate breaking changes into one or several releases. The decision criterion is how much work you create for “client developers” (maintainers of reverse dependencies).

It’s also one of the many places in the book where the author reminds us there’s no hard rule, that it’s the “art of software engineering”.

House or not House

The first chapter of the book presents and criticizes analogies of software engineering, like comparing developing software to building a house. That reminded of the talk Maintaining the house the tidyverse built by Hadley Wickham, that compared tidyverse maintenance to house maintenance.

The chapter’s thesis is that all analogies are useful and imperfect and you shouldn’t let them limit or cloud the way you view your work.

Justify exceptions

At some point the author writes that static code analysis brings false positives, but that if you disable specific rules, you should justify why. That’s how Jarl, a linter for R code, makes you suppress rules: you have to explain a reason.

Git is easy

I like these encouraging three sentences:

“Git isn’t the most user-friendly piece of technology on the planet, but you’re a programmer. You’ve managed to learn at least one programming language. Compared to that, learning the basics of Git is easy.”

Testing against databases

As an example of slow tests, the author mentions tests against databases. I want to use this as an excuse to plug the dittodb package to mock databases in tests.

Bisection without and with Git

In the chapter about Troubleshooting, the idea of bisection without Git is introduced. The author says it’s called bisection for lack of a better word. The idea being to go from a big piece of code with a bug to the smallest piece of code with the same bug, I suppose a better word is reprex.

“Being able to produce a minimal working example is a superpower in software engineering.”

Regarding Git bisection, here’s your friendly reminder that you can try it out using the saperlipopette R package.

Time boxing

The author extols the virtues of time boxing (e.g. the Pomodoro method). I’ve been using Entracte and would recommend giving it a go, if you have no similar setup yet!

When reading the book, between chunks of reading I’d crochet a few rounds.

Performance

I enjoyed the pragmatic discussion of performance: don’t forget to check the numbers and their significance.

Behavioural code analysis

I have never used behavioural code analysis, something that uses Git data, but its presentation reminded me of Git commands to get to know a project which is related (or the same?).

Code navigation

Yay to the mention of learning how to navigate code in your IDE. Such an important skill.

Conclusion

“Code that Fits In Your Head” is an ambitious (and rather long) book. I was glad to learn new phrases and to get the opportunity to get to know or reflect on so many topics.

Because the book was published in 2021, AI only makes an appearance as a brief footnote. The only moment I really thought AI would make a point less valid was a point about some kinds of refactoring taking ages: in the igraph R package, some changes have recently become feasible in an easier way thanks to using tools like Claude.

Some aspects of the book might be a tad annoying like the super simple diagrams or illustrations (think: a hammer drawing to illustrate the fact that to someone with a hammer, everything looks like a nail) but I suppose those play their role of breaking up pages.

Bonus: a reading list

If someone fairly new to software engineering were to ask me book recs, from looking at the stack near my desk I would recommend the following books, more specialized than “Code that Fits In Your Head”:

• The Art of Readable Code by Dustin Boswell and Trevor Foucher
• The Programmer’s Brain by Felienne Hermans (that I plan to re-read soon!)
• A Philosophy of Software Design by John Ousterhout
• A Git book, either Git in Practice by Mike McQuaid or Pro Git by Scott Chacon
• A book about team work. Maybe The Psychology of Software Teams by Cat Hicks, The Collective Edge by Colin M. Fisher, Dare to Lead by Brené Brown…
• A book about productivity/organization, but not a sanctimonious one about profound labor or the like. I enjoyed Make Time by Jake Knapp and John Zeratsky and Tranquility by Tuesday by Laura Vanderkam.
• A book critical of AI such as The AI Con by Emily M. Bender and Alex Hanna or Empire of AI by Karen Hao.

I also believe there are other ways to learn, like watching talks, but this post is about books.