R-bloggers

maestro 1.4.0 adds custom pipeline metrics

Will Hipson — Thu, 15 Oct 2026 03:00:00 +0000

[This article was first published on data-in-flight, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

The 1.4.0 release of maestro is out, and with it comes the ability to record custom metrics from your pipelines. Execution time and success/failure were already tracked, but there was no good way to answer questions like “how many rows did this pipeline process?” or “how many API calls did it take?” Now there is.

If you haven’t heard of maestro, it’s a pipeline orchestration package. You can learn more about it here.

Get it from CRAN:

install.packages("maestro")

Recording metrics with `update_metric()`

A metric is any value that can be represented as a single number: rows processed, bytes read, API calls made, attempts on a task. To record one, call update_metric() from inside a pipeline with the name of the metric and its value. maestro automatically attributes the metric to the pipeline (and run) that is currently executing.

An important design point is that metrics are recorded on the side. They are not part of a pipeline’s return value, so they never get passed downstream in a DAG.

Here’s a two-step DAG that extracts some flight data and loads it, recording a few metrics along the way:

#' ./pipelines/metric_example.R
#' @maestroFrequency 1 hour
extract_flights <- function() {
  
  # Pretends to fetch data from an API
  df <- data.frame(
    flight_id = 1:5,
    airline_code = c("UA", "AC", "AC", "AA", "DE"),
    departure_time = as.POSIXct("2026-10-14 12:00:00") + c(100, 450, 750, 1450, 1750)
  )
  
  update_metric("rows_read", nrow(df))
  update_metric("bytes_read", as.numeric(object.size(df)))
  
  df
}

#' @maestroInputs extract_flights
load_flights <- function(.input) {
  
  write.csv(.input, tempfile(fileext = ".csv"))
  
  update_metric("rows_written", nrow(.input))
}

Retrieving metrics with `get_metrics()`

After running the schedule, get_metrics() returns a long data.frame of every metric recorded. Each row gives the pipeline it came from (pipe_name), the run it was recorded in (run_id), the metric name, and its numeric value.

library(maestro)

schedule <- build_schedule(quiet = TRUE)

status <- run_schedule(
  schedule,
  orch_frequency = "1 hour",
  quiet = TRUE
)

get_metrics(schedule)

# A tibble: 3 × 4
  pipe_name       run_id metric       value
                       
1 extract_flights vPnJoC rows_read        5
2 extract_flights vPnJoC bytes_read    1648
3 load_flights    ulLl5m rows_written     5

Notice that each pipeline gets its own run_id, even within a DAG. Because the result is a plain data.frame, it’s easy to summarise further or write out to a database or file, which makes it a natural fit for monitoring pipeline health over time.

A couple of other things to know:

Metrics are only recorded when a pipeline is executed by run_schedule(). Calling update_metric() while developing or testing a pipeline function interactively is harmless: it emits an informational message and does nothing.
Execution time and success/failure are already available via get_status() and get_schedule(), so there’s no need to record those yourself.

Metrics are per execution

Metrics belong to a single execution of a pipeline. Each time a pipeline runs it gets a new run_id, and its metrics are stored separately from those of earlier runs. Within one execution, calling update_metric() again with the same name overwrites the earlier value (the last write wins).

Running the schedule a second time shows this: the metrics from both runs are retained and distinguished by run_id.

status <- run_schedule(
  schedule,
  orch_frequency = "1 hour",
  quiet = TRUE
)

get_metrics(schedule)

# A tibble: 6 × 4
  pipe_name       run_id metric       value
                       
1 extract_flights vPnJoC rows_read        5
2 extract_flights vPnJoC bytes_read    1648
3 extract_flights bilN7s rows_read        5
4 extract_flights bilN7s bytes_read    1648
5 load_flights    ulLl5m rows_written     5
6 load_flights    R9iMP0 rows_written     5

Accumulating with `increment = TRUE`

Sometimes a metric is built up over the course of a pipeline rather than computed once, like counting API calls inside a loop or the number of attempts on a flaky task. Setting increment = TRUE adds value to the metric’s current value instead of overwriting it. If the metric hasn’t been set yet in the current run, the first call simply sets it. A negative value decrements.

#' ./pipelines/metric_increment.R
#' @maestroFrequency 1 day
fetch_paginated <- function() {
  
  results <- list()
  
  for (page in 1:3) {
    attempts <- 0
    
    repeat {
      attempts <- attempts + 1
      update_metric("attempts", 1, increment = TRUE)
      
      # Pretends to call an API that sometimes fails
      ok <- runif(1) > 0.5
      if (ok || attempts >= 3) break
    }
    
    update_metric("api_calls", 1, increment = TRUE)
    results[[page]] <- data.frame(page = page, value = rnorm(10))
  }
  
  df <- do.call(rbind, results)
  
  # Overwriting (the default) is fine for a one-off value
  update_metric("rows_read", nrow(df))
  
  df
}

In this example both api_calls and attempts accumulate. When it finishes, api_calls is 3 (one per page) while attempts is at least 3, depending on how many retries were needed. Each metric still ends up with a single value per pipeline run.

For the full details, see the metrics vignette.

Other recent changes

If you missed the 1.3.x releases, there are a couple of smaller things worth mentioning:

The new @maestroCascadeTags tag propagates selected metadata tags (label, flags, and/or loglevel) from a pipeline to all of its downstream pipelines in a DAG. With no value, it cascades all three.
get_slot_usage() gained min_datetime and max_datetime parameters so the time window can be fixed rather than anchored to the current time.

The full list of changes is in the release notes.

As always, happy orchestrating!

To leave a comment for the author, please follow the link and comment on their blog: data-in-flight.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Continue reading: maestro 1.4.0 adds custom pipeline metrics

A Trip from Covariance to Correlation and Back

Andrea Onofri — Wed, 14 Oct 2026 22:00:00 +0000

[This article was first published on Statforbiology, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

When we have a multi-trait table of data, we might be interested in investigating variability structure and the reciprocal relationships (whether causal or non-causal) among variables. For example, let’s consider the ‘mtcars’ dataset, shown below:

data(mtcars)
head(mtcars)
##                    mpg cyl disp  hp drat    wt  qsec vs am gear carb
## Mazda RX4         21.0   6  160 110 3.90 2.620 16.46  0  1    4    4
## Mazda RX4 Wag     21.0   6  160 110 3.90 2.875 17.02  0  1    4    4
## Datsun 710        22.8   4  108  93 3.85 2.320 18.61  1  1    4    1
## Hornet 4 Drive    21.4   6  258 110 3.08 3.215 19.44  1  0    3    1
## Hornet Sportabout 18.7   8  360 175 3.15 3.440 17.02  0  0    3    2
## Valiant           18.1   6  225 105 2.76 3.460 20.22  1  0    3    1

Considering the different cars (along the rows), we might be interested in investigating the structure of variability for fuel consumption (mpg), number of cylinders (cyl), displacement (disp) and horsepower (hp). How much does each variable change across the different cars? How much does each variable change when another variable changes? Is there any kind of joint variation, either positive or negative?

One suitable statistic to answer these question is the covariance, that is:

where and are the means for each variable and is the number of individuals. In words, covariance is the sum of products of the i^th deviations from the mean for both variables. In this setting, the variance is just the covariance of a variable with itself.

For a data table, variances and covariances can be calculated by using the cov() function, that returns a matrix (the variance-covariance matrix, indeed):

matr <- mtcars[,1:4]

# Covariances
Sigma <- cov(matr)
Sigma
##              mpg        cyl       disp        hp
## mpg    36.324103  -9.172379  -633.0972 -320.7321
## cyl    -9.172379   3.189516   199.6603  101.9315
## disp -633.097208 199.660282 15360.7998 6721.1587
## hp   -320.732056 101.931452  6721.1587 4700.8669

While covariances indicate direction (positive or negative association), their values depend on the scale of measurement, making it difficult to gauge the strength of the relationship at a glance. For this reason, covariances are usually ‘standardized’ by dividing by the product of the standard deviations to obtain Pearson’s correlation coefficient:

where and are the standard deviations for and . Unlike covariances, correlations go from -1 to 1 and the correlation matrix has 1s along its main diagonal (since the correlation of any variable with itself is 1). In R, the matrix of correlation is obtained by using the cor() function:

# Correlations
R <- cor(matr)
R
##             mpg        cyl       disp         hp
## mpg   1.0000000 -0.8521620 -0.8475514 -0.7761684
## cyl  -0.8521620  1.0000000  0.9020329  0.8324475
## disp -0.8475514  0.9020329  1.0000000  0.7909486
## hp   -0.7761684  0.8324475  0.7909486  1.0000000

This brings us to the main question of this post: Given a variance-covariance matrix, can we compute the correlation matrix? And can we do the reverse—moving from a correlation matrix back to a variance-covariance matrix? With reference to the matrices above, can we go from ‘Sigma’ to ‘R’ and from ‘R’ to ‘Sigma’? I often find myself needing to make these conversions, especially when working with mixed models or time series analysis.

From covariance to correlation

Considering the previous expression, it is easy to see that, if we take the covariance between ‘disp’ and ‘mpg’ above (-633.0972), the correlation can be obtained as:

-633.097208 / (sqrt(36.324103) * sqrt(15360.7998))
## [1] -0.8475514

and all the necessary information is contained in the variance-covariance matrix.

Of course, doing these calculations for all covariances one by one is tedious. The question is: is there any simple way to go from a covariance matrix to a correlation matrix with one single operation?

The answer is found in matrix algebra; let’s take a look at the variance-covariane matrix: what we have to do is to take each element in this matrix and divide it by the square roots of the diagonal elements in the same column and in the same row (see figure below). Consequently, our task is:

take the diagonal elements
take each i^th row and divide all its elements by the i^th diagonal element (row-scaling)
take each i^th column of the resulting matrix and divide all its elements by the i^th diagonal element (column-scaling)

Row-scaling can be performed by using matrix-algebra; in particular, we can take a diagonal matrix of weights and multiply it by the square matrix we want to scale, as shown below.

We see that all elements of the first row were multiplied by the first weighting element, all elements of the second row were multiplied by the second weighting element and so forth.

Column scaling is done similarly, but in reverse order: we take the target matrix and post-multiply it by the diagonal matrix of weights, as shown below.

We see that all elements of the first column were multiplied by the first weighting element, all elements of the second row were multiplied by the second weighting element, and so forth.

Consequently, we can take our variance-covariance matrix and:

Take the square roots of the diagonal elements (standard deviations) and load them into a diagonal matrix
Invert this diagonal matrix
Pre-multiply and post-multiply the variance-covariance matrix by this diagonal matrix of inverted standard deviations

# from covariance to correlation - 1
StDev <- sqrt(diag(Sigma))
StDevMat <- diag(StDev)
InvStDev <- solve(StDevMat)
InvStDev %*% Sigma %*% InvStDev
##            [,1]       [,2]       [,3]       [,4]
## [1,]  1.0000000 -0.8521620 -0.8475514 -0.7761684
## [2,] -0.8521620  1.0000000  0.9020329  0.8324475
## [3,] -0.8475514  0.9020329  1.0000000  0.7909486
## [4,] -0.7761684  0.8324475  0.7909486  1.0000000

The previous solution is correct, but, in R, it is rather inefficient. The best way is to use the cov2cor() function in the ‘stat’ package.

# From covariance to correlation - 2
cov2cor(Sigma)
##             mpg        cyl       disp         hp
## mpg   1.0000000 -0.8521620 -0.8475514 -0.7761684
## cyl  -0.8521620  1.0000000  0.9020329  0.8324475
## disp -0.8475514  0.9020329  1.0000000  0.7909486
## hp   -0.7761684  0.8324475  0.7909486  1.0000000

With base R, we can also sweep() twice:

# From covariance to correlation - 3
sweep(sweep(Sigma, 1, StDev, FUN = "/"), 2, StDev, FUN = "/")
##             mpg        cyl       disp         hp
## mpg   1.0000000 -0.8521620 -0.8475514 -0.7761684
## cyl  -0.8521620  1.0000000  0.9020329  0.8324475
## disp -0.8475514  0.9020329  1.0000000  0.7909486
## hp   -0.7761684  0.8324475  0.7909486  1.0000000

Alternatively, we can combine scale() and t() (transpose) twice, though it is noticeably less intuitive:

# From covariance to correlation - 4
scale(t(scale(t(Sigma), center = F, scale = StDev)), 
      center = F, scale = StDev)
##             mpg        cyl       disp         hp
## mpg   1.0000000 -0.8521620 -0.8475514 -0.7761684
## cyl  -0.8521620  1.0000000  0.9020329  0.8324475
## disp -0.8475514  0.9020329  1.0000000  0.7909486
## hp   -0.7761684  0.8324475  0.7909486  1.0000000
## attr(,"scaled:scale")
##        mpg        cyl       disp         hp 
##   6.026948   1.785922 123.938694  68.562868

Another very clean approach in R is based on the outer product. If we have a vectors of standard deviations ($), the outer product is:

Thus, the correlation is calculated as:

# From covariance to correlation - 5
Sigma/outer(StDev, StDev)
##             mpg        cyl       disp         hp
## mpg   1.0000000 -0.8521620 -0.8475514 -0.7761684
## cyl  -0.8521620  1.0000000  0.9020329  0.8324475
## disp -0.8475514  0.9020329  1.0000000  0.7909486
## hp   -0.7761684  0.8324475  0.7909486  1.0000000

From correlation to covariance

It is easy to see that:

Thus, if we have the correlation between ‘disp’ and ‘mpg’ (-0.8521620), we can easily calculate the covariance as:

-0.8475514 * sqrt(36.324103) * sqrt(15360.7998)
## [1] -633.0972

The only problem is that the variances are not contained in a correlation matrix, so we need to get them from elsewhere. Once we have them, we

1. Take their square roots and place them into a diagonal matrix
Pre-multiply and post-multiply the correlation matrix for this diagonal matrix of standard deviations

The solution is as follows:

# from correlation to covariance - 1
StDev <- sqrt(diag(Sigma))
StDevMat <- diag(StDev)
StDevMat %*% R %*% StDevMat
##             [,1]       [,2]       [,3]      [,4]
## [1,]   36.324103  -9.172379  -633.0972 -320.7321
## [2,]   -9.172379   3.189516   199.6603  101.9315
## [3,] -633.097208 199.660282 15360.7998 6721.1587
## [4,] -320.732056 101.931452  6721.1587 4700.8669

More efficiently, we can use the cor2cov() function in the ‘lavaan’ package.

# From correlation to covariance - 2
lavaan::cor2cov(R, StDev)
##              mpg        cyl       disp        hp
## mpg    36.324103  -9.172379  -633.0972 -320.7321
## cyl    -9.172379   3.189516   199.6603  101.9315
## disp -633.097208 199.660282 15360.7998 6721.1587
## hp   -320.732056 101.931452  6721.1587 4700.8669

Or, we can sweep() twice:

# From correlation to covariance - 3
sweep(sweep(R, 1, StDev, FUN = "*"), 2, StDev, FUN = "*")
##              mpg        cyl       disp        hp
## mpg    36.324103  -9.172379  -633.0972 -320.7321
## cyl    -9.172379   3.189516   199.6603  101.9315
## disp -633.097208 199.660282 15360.7998 6721.1587
## hp   -320.732056 101.931452  6721.1587 4700.8669

or scale() and t() twice:

# From correlation to covariance - 5
scale(t(scale(t(R), center = F, scale = 1/StDev)), 
      center = F, scale = 1/StDev)
##              mpg        cyl       disp        hp
## mpg    36.324103  -9.172379  -633.0972 -320.7321
## cyl    -9.172379   3.189516   199.6603  101.9315
## disp -633.097208 199.660282 15360.7998 6721.1587
## hp   -320.732056 101.931452  6721.1587 4700.8669
## attr(,"scaled:scale")
##         mpg         cyl        disp          hp 
## 0.165921457 0.559934979 0.008068505 0.014585154

Finally, we can use the outer product as follows:

# From correlation to covariance - 6
R * outer(StDev, StDev)
##              mpg        cyl       disp        hp
## mpg    36.324103  -9.172379  -633.0972 -320.7321
## cyl    -9.172379   3.189516   199.6603  101.9315
## disp -633.097208 199.660282 15360.7998 6721.1587
## hp   -320.732056 101.931452  6721.1587 4700.8669

Do you have a cleaner or faster trick for this in R? Drop a comment below and I’ll be happy to feature it here!

Thanks for reading, and happy coding!

(P.S. If you enjoyed this post, check out my book “Field Research Methods in Agriculture: An Introduction with R” via the link below!)

Prof. Andrea Onofri
Department of Agricultural, Food and Environmental Sciences
University of Perugia (Italy)
Send comments to: andrea.onofri@unipg.it

To leave a comment for the author, please follow the link and comment on their blog: Statforbiology.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Continue reading: A Trip from Covariance to Correlation and Back

Dogs of New York

R on kieranhealy.org — Wed, 14 Oct 2026 17:17:28 +0000

[This article was first published on R on kieranhealy.org, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

I recently updated nycdogs, a data package for R that I originally made for teaching purposes. It bundles up New York City’s dog licensing and dog bite datasets into an R package. The data is already provided by the city in a directly usable way (as CSV files), but packaging it up makes it easier to integrate with other packages and get things moving a little more quickly in class.

The nyc_license table looks like this:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15


nyc_license
#> # A tibble: 819,323 × 12
#>    animal_name animal_gender animal_birth_year breed_name       breed_rc zip_code zip   license_issued_date license_expired_date extract_year borough city 
#>                                                                                              
#>  1 Paige       F                          2014 American Pit Bu… Pit Bul… 10035    10035 2014-09-12          2017-09-12                   2016 Manhat… New …
#>  2 Yogi        M                          2010 Boxer            Boxer    10465    10465 2014-09-12          2017-10-02                   2016 Bronx   Bronx
#>  3 Ali         M                          2014 Basenji          Basenji  10013    10013 2014-09-12          2019-09-12                   2016 Manhat… New …
#>  4 Queen       F                          2013 Akita Crossbreed Akita C… 10013    10013 2014-09-12          2017-09-12                   2016 Manhat… New …
#>  5 Lola        F                          2009 Maltese          Maltese  10028    10028 2014-09-12          2017-10-09                   2016 Manhat… New …
#>  6 Ian         M                          2006 Unknown          Unknown  10013    10013 2014-09-12          2019-10-30                   2016 Manhat… New …
#>  7 Buddy       M                          2008 Unknown          Unknown  10025    10025 2014-09-12          2017-10-20                   2016 Manhat… New …
#>  8 Chewbacca   F                          2012 Labrador Retrie… Labrado… 10013    10013 2014-09-12          2019-10-01                   2016 Manhat… New …
#>  9 Heidi-bo    F                          2007 Dachshund Smoot… Dachshu… 11215    11215 2014-09-13          2017-04-16                   2016 Brookl… Broo…
#> 10 Massimo     M                          2009 Bull Dog, French French … 11201    11201 2014-09-13          2017-09-17                   2016 Brookl… Broo…
#> # ℹ 819,313 more rows

The dataset is good for teaching with for several reasons. First, people like dogs. Second, it has about 800,000 records, which is a lot in many ways but also trivial for a modern laptop to work with. Third, the data are pretty clean, but they’re not that clean. There are various things going on in there that invite us to explore it but also to be careful about what we can and can’t be confident about.

The first thing we ask about any data set like this is: what is each row in this table? In the case of nyc_license, each row is a record of one time-limited license for a dog. A row is not necessarily the record of a unique dog. When an owner renews a dog’s license, the dog will appear one time for each license period. In New York City, dog licenses can be issued for one year or five years. The table does not provide any dog_id column that ties sequences of licenses to individual dogs.

Second, the table contains several “extracts” of the license data. It seems like the way the provider updates the NYC Open Data record is to periodically pull from its internal database every year or so and append it to the existing, previously-released version. The extract_year column identifies the extract, so at least we know that. But it’s also the case that a license that was active in more than one extract year appears in each extract. As a result, the table has a large number of duplicate rows.

Third, license information is provided by dog owners. This means that there’s a fair amount of noise in the records, because people make typos and other errors. For instance, the ZIP code data has incomplete ZIPs. Some people provide invalid ZIP codes (e.g. ones that are fewer than five digits long). Others (about 3,300 records) have ZIP codes that are not in New York City. Eighty percent of these (about 2,700) are ZIP codes located elsewhere in New York state. Another 77 are located in Florida, which I kind of like. Five of them live in Boca Raton, for example. Their kids never call. You’d think it would be easy for them to just pick up the bone, I mean phone, but you never hear from them.

About 650 rows have valid ZIP codes that are in Manhattan, but these ZIPs are not ZCTAs. The trouble with ZIP codes is that people think of them as areas, as polygons, but in fact they are lines. They are defined by the Post Office as postal routes, i.e. as journeys a mail carrier takes. While they tend to form closed polygons, they are absolutely not guaranteed to. But people know their ZIP code and all kinds of analysts (from marketers to political operatives) like thinking in terms of ZIP codes. So the Census defines ZIP Code Tabulation Areas that are a kind of best-effort at answering the question “What if this ZIP code were a closed polygon?” ZCTAs exclude every P.O. box and single-address ZIP nationwide. Manhattan has several buildings that have their own ZIP code. New York City also defines MODZCTAs that merge some low- or no-population ZCTAs with their neighbors. I use these below.

Fourth, beyond the ZIP codes, owners also provide information about the age and sex of their dogs, the dog’s name, and the dog’s breed. Again, this kind of information tends to be a little noisy. For example, people vary in how they spell names. They also vary in how they classify their dog by breed. Beyond the list of “recognized” breeds (itself something that dog owners variably conform to), some people might describe their dog as, e.g., “Labrador-Sheltie mix” while others might just say “Unknown”. The dataset includes the breed_name column from the original dataset and also breed_rc, which is a recoding by me to clean up equivalent terms (“Yorkshire Terrier” and “Terrier, Yorkshire”, for example) along with a few other coding decisions.

Finally, and only because I’m choosing to stop here, the table just shows the number of licensed dogs in New York City. That is definitely not the same as the number of dogs there. Moreover, license compliance varies quite a bit in ways that are related to characteristics of the owner that we can’t observe, and also ones that we partially can, like where they live.

All of which is related to the fact that one of the lolsob-mantras of courses I teach is “Now you know the true meaning of data”. In the graphs and maps below, I’ve done the best I can, where appropriate, to uniquely identify individual dogs in the dataset (even though ultimately we cannot).

Popular Breeds

With that out of the way, it’s time to indulge your prejudices about Dogs, New Yorkers, and New York Dogs. The most common “breed” of dog in New York is “Unknown”, as befits a city of immigrants, hustlers, scrappy up-and-comers, and dogs seeking to leave their small-town past behind and remake their identities. Somewhere between eight and ten percent of records each year have “Unknown” recorded as their breed. There’s some change over time (a declining share), but I think this is an artifact of how breeds are getting recorded rather than a change in the composition of the dog population. If we leave the Unknowns out, here are the trends in annual licenses issued to the twenty most popular breeds.

Top Dogs.

Some breeds are pretty steadily popular. Smaller breeds do well in a city full of apartments. We can also see some evidence here of trends in popularity. However, these annual numbers are what we might call the prevalence of dogs, rather than the incidence of new dogs. If we want to get a better sense of fashion trends in dog breeds, we should look at puppies. (Who doesn’t want to look at puppies?) Because of the structure of the data we can’t quite do this. The best we can do, based on the dog’s age, is to examine licenses issued to dogs younger than a year old. Let’s pick out some breeds that seem to be rising, a few that seem relatively steady, and a few that might be declining. I’ve indexed these changes relative to each breed in 2016 by setting the value for that year to 100. Even though our range is narrow, we use a log scale on the y-axis because, e.g., a rise in a breed’s index from 100 to 200 reflects the same relative change as a fall from 100 to 50.

Maltipoos and Goldendoodles seem to be doing pretty well for themselves. Golden Retrievers remain steadily popular. Meanwhile, Shiba Inus, beloved fad-breed of the 2010s, seem to have been affected by their association with crypto and are strongly in decline. Much drop. Very decline. Wow. Beagles are also doing poorly. Good grief.

Where Breeds Live

We can look at the relative prevalence of breeds by their ZIP code (strictly speaking, by their MODZCTA). I hesitate a little on this just because, even with 800,000 licenses, we’ll run out of degrees of freedom fairly quickly if we break things down by breed and ZIP. Relatedly, we should ask what sort of metric we’re trying to show here. What’s the denominator? (A surprising amount of Statistics, I claim, is essentially exploring the question “What happens if I divide by this number instead of that number?”) Popular dogs are, to a first approximation, pretty popular all over the city. Most ZIP codes have a solid contingent of Yorkies. If, for each popular breed, we calculate the share of that breed as a percentage of all dogs in the ZIP code, we’ll see some ZIP-by-ZIP variation come out.

It would also be nice to adjust for the fact that different ZIP codes have quite different shares of the total number of dogs licensed in the city. (Which, again, is not the same as share of the total number of dogs in the city.) So we calculate the within-ZIP share and make a ratio of it by dividing that number by the breed’s citywide share. The ratio is 1 when the breed is just as common in that ZIP code as it is across the whole city. If it’s 2, the breed is twice as common in that ZIP as it is across the city. If it’s 0.5, the breed is half as common in that ZIP as it is in the city as a whole. The main downside of ratios like this is that they can exaggerate differences a little, especially if the base rates are small. We have enough breeds in each zip code to make me feel OK about doing this. But if you refer to the bar charts above you can see that even for disproportionately popular breeds in zips, you’re looking at changes like “twice as common with respect to a baseline of five percent of all dogs”. Here’s the map:

We see some plausible evidence that breeds are distributed across ZIP codes in a manner that tracks the demography of owners. French Bulldogs in Chelsea, Soho, and TriBeCa. Havanese on the Upper East Side. Pit Bulls in the Bronx. German Shepherds on Staten Island. Goldendoodles and King Charles Spaniels in Manhattan south of 96th.

The Diffusion of Goldendoodles

Let’s take a closer look at the Goldendoodles. In 2015, Goldendoodles weren’t that common. We can ask: where were they initially popular? Did they diffuse across New York City from focal points where people were more likely to set trends in breed fashion? Restricting ourselves to young dogs, we can try to look at each ZIP’s share of all Goldendoodle puppies licensed each year from 2015 to 2025. Here’s what that looks like:

Doodle diffusion.

In 2015, Goldendoodle puppies are licensed in a relatively small number of ZIP codes, mostly in Manhattan. Most ZIP codes are doodle-free. Over the next five years, they gradually spread across the city. By 2023 or 2024, most—though still not all—ZIP codes are licensing young Goldendoodles.

Again, the underlying numbers aren’t huge. So take all of this with a delicious treat. I mean, grain of salt. I also wonder what the next few years are going to look like in terms of dog popularity. Based on the maps above, I think Australian Shepherds are likely to continue to rise in popularity amongst owners looking to acquire a fashionable breed. And while they’re not on any of the maps above, the data also suggest that the Jindo, a Korean breed that is a relatively recent arrival in the US, has been sharply on the rise. If you ask me, though, if you’re in New York City and planning to get a dog, you should throw in your lot with Unknown and get a good, solid, dog-shaped mutt. Preferably one with strong views on the state of the Subway these days, or who knows a guy who can get you a deal on some shirts, or the special ramen, or tickets. Not tickets for anything in particular. Just, like, tickets. Works out of an upstairs bar downtown on Water St. Trust me.

To leave a comment for the author, please follow the link and comment on their blog: R on kieranhealy.org.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Continue reading: Dogs of New York

A Bayesian Network to Explore Data Center Risks

Joseph Rickert — Tue, 13 Oct 2026 00:00:00 +0000

[social4i size="small" align="align-left"] -->

[This article was first published on R Works, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Data centers required to fuel the ambitions of Artificial Intelligence companies have sparked a nationwide conversation. AI companies are moving quickly to build out infrastructure and local communities are pushing back to resist what they see as the appropriation of public resources such as electricity and water to further enrich wealthy corporations. Moreover, the sentiment shared by both the political left and right is that data centers will significantly degrade the quality of life for those living adjacent to a data center and to the extended community as well. This post implements a non-trivial Bayesian Network to study the risks involved. I hope it is of some interest to a broad audience including, people interested in clarifying their thoughts about the risks that data centers may pose to their particular community and analysts curious about Bayesian Networks. However, I think the sample model built out in this post will be especially interesting to decision analysts who know quite a bit about Bayesian Networks but may be unfamiliar with the R programming language and unaware that there is capable opensource software which can provide capabilities beyond those offered by proprietary software systems.

Other than addressing the topic of data center risks the analysis contains no original research except perhaps the contribution of ChatGPT 5.6 Luna which conjured up the innovative parts of the synthetic data generation algorithm from a series of requirements.

Data Center Risks

Rational discussion about the perceived risks associated with the construction of data centers could benefit from a disciplined approach to organizing information and opinions into risk a analysis flexible enough to the consequences of complex interactions among risk factors. The Bayesian Network presented in this post is intended to be first step towards building such a tool.

In general, the risks imposed by data centers fall into the following general categories. See Appendix C below for some references to investigative reporting supporting these risk categories.

Electricity. A data center’s power demand can affect the immediate substation, the regional transmission system, and electricity rates for customers across a utility territory. The key questions are who pays for new generation and transmission, whether the facility has a dedicated rate, and whether its demand is flexible during shortages.
Water. Water use depends heavily on cooling design, climate, server load, and the use of reclaimed water. The relevant impact area is often the watershed rather than the property boundary. During droughts, a data center may compete with households, agriculture, industry, and ecosystems.
Land and housing. The facility itself may occupy hundreds of acres when parking, substations, cooling systems, security setbacks, and transmission corridors are included. Construction can remove farmland or open space, while the associated workforce and infrastructure can increase housing demand.
Noise and air quality. The most immediate quality-of-life impacts are usually construction activity, cooling-system noise, transformer hum, generator testing, lighting, and occasional diesel emissions. These effects decline with distance but can remain significant for nearby homes, schools, and outdoor spaces.
Grid reliability and affordability. A data center can bring grid upgrades that improve regional infrastructure, but rapid demand growth can also strain generation and transmission. If costs are not assigned directly to the facility, residential and small-business customers may absorb part of them.
Governance and fairness. The distance framework should be paired with an equity analysis: who receives jobs, tax revenue, and infrastructure improvements, and who experiences noise, land loss, water competition, higher rates, or pollution? Community-benefit agreements, transparent resource-demand disclosures, enforceable noise and water standards, and clear cost-recovery rules can substantially change the outcome.

Potential Benefits

Supporters of data centers emphasize the potential benefits of building a data center. These may include increased tax revenue and increased temporary employment. However, most of the reporting referenced below indicates that they do not compensate for the risks.

Organizing the risks and benefits

Model building begins by organizing nodes into a a Directed Acyclic Graph (DAG) with arrows pointing from root nodes which are thought to drive other nodes that represent factors that intensify risk. Eventually the information in these nodes flow to nodes which quantify different kinds of risk. Finally the information from the risk nodes flow into a node that quantifies overall community risk. The following graph expands on the general categories listed above, enumerates a number of risk factors with arrows indicating the direction influence or causality. The names of the nodes are for the most part self explanatory. Note that the scale of the data center and distance from the data center are fundamental factors driving risk.

The number of nodes in this graph are probably much larger than would be required for an exploratory mode. However, I have elected to go with the larges list that seemed to me to be reasonable and coherent in order to provide an example that is more complex than a typical “toy” mode. And since it is usually much easier to remove code than add it, the model below should be useful in deriving a more parsimonious custom model.

%%| label: fig-risk-network
%%| fig-cap: "Data Center Risk Network"
%%| fig-alt: "Directed network showing data-center impacts and risks leading to overall community risk."



flowchart LR

  Scale --> LandFootprint
  Scale --> ElectricityDemand
  Scale --> WaterDemand
  Scale --> ConstructionIntensity
  Scale --> BackupGeneration
  Scale --> WorkforceDemand
  Scale --> NoiseExposure
  Scale --> AirQualityExposure
  Scale --> AccessToBenefits
  Scale --> EconomicBenefit

  Cooling --> WaterDemand
  Cooling --> NoiseExposure

  Distance --> NoiseExposure
  Distance --> AirQualityExposure
  Distance --> TrafficExposure
  Distance --> VisualLandUseExposure
  Distance --> AccessToBenefits
  Distance --> ElectricityCostRisk

  LandFootprint --> VisualLandUseExposure
  ConstructionIntensity --> TrafficExposure
  BackupGeneration --> AirQualityExposure

  ElectricityDemand --> GridStress
  GridCapacity --> GridStress
  GridStress --> ElectricityCostRisk
  GridStress --> ReliabilityRisk

  WaterDemand --> WaterCompetition
  WaterAvailability --> WaterCompetition
  WaterCompetition --> WaterSecurityRisk
  Vulnerability --> WaterSecurityRisk

  WorkforceDemand --> HousingPressure
  HousingAvailability --> HousingPressure
  WorkforceDemand --> ServicePressure
  HousingAvailability --> ServicePressure

  LandFootprint --> EnvironmentalDegradation
  WaterDemand --> EnvironmentalDegradation
  WaterAvailability --> EnvironmentalDegradation

  NoiseExposure --> QualityOfLifeRisk
  AirQualityExposure --> QualityOfLifeRisk
  TrafficExposure --> QualityOfLifeRisk
  VisualLandUseExposure --> QualityOfLifeRisk
  HousingPressure --> QualityOfLifeRisk
  Vulnerability --> QualityOfLifeRisk

  AccessToBenefits --> EconomicBenefit

  QualityOfLifeRisk --> OverallCommunityRisk
  WaterSecurityRisk --> OverallCommunityRisk
  ElectricityCostRisk --> OverallCommunityRisk
  ReliabilityRisk --> OverallCommunityRisk
  EnvironmentalDegradation --> OverallCommunityRisk
  HousingPressure --> OverallCommunityRisk
  ServicePressure --> OverallCommunityRisk
  EconomicBenefit --> OverallCommunityRisk

The DAG

The code below builds a discrete Bayesian Network based on the above directed acyclic graph using bnlearn package.

Required packages and helper functions

Creating a permutation test, two ways

Gavin Simpson — Mon, 12 Oct 2026 22:00:00 +0000

[This article was first published on From the bottom of the heap — R posts, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

In my last post on permute, I looked at ordered set partitions and the new formula interface for plots and blocks. This time I want to show two helpers for the next part of writing a permutation test: taking those permutations and doing something with them.

permutations() gives you the shuffled values directly. permutator() gives you a function that hands out the next set of observation indices each time you call it. Both use the same permutation machinery underneath, so we can use them to build the same test in two slightly different ways.

These functions are in the development version of permute. To follow along, install it from my R-universe:

install.packages("permute", repos = c(
  "https://gavinsimpson.r-universe.dev",
  "https://cloud.r-project.org"
))

Something different from jackals

Let’s use PlantGrowth, which comes with R. It records dried plant weights under a control and two treatment conditions:

library("permute")
data(PlantGrowth, package = "datasets")

aggregate(weight ~ group, data = PlantGrowth, FUN = mean)

  group weight
1  ctrl  5.032
2  trt1  4.661
3  trt2  5.526

The group means differ, but is that more variation than we’d expect if plant weights were unrelated to the treatments? We can use the overall F statistic from a linear model to measure how different the groups are. Larger values mean more variation between groups relative to the variation within them.

I’ll wrap that calculation in a small function, because we’ll be doing it quite a few times:

f_stat <- function(weight, group) {
  fit <- lm(weight ~ group)
  unname(summary(fit)$fstatistic[["value"]])
}

observed <- with(PlantGrowth, f_stat(weight, group))
observed

[1] 4.846088

This is an overall test across the groups. It doesn’t tell us which particular treatments differ.

For the permutation test, we’ll leave the group labels in place and shuffle the weights among them. Under the null hypothesis, we’re assuming those weights are exchangeable across the treatment groups.

Give me the shuffled values

Usually, shuffleSet() returns permuted indices, which we then use to rearrange our data. permutations() does that indexing step for us:

ctrl <- how(nperm = 9999)

set.seed(42)
shuffled <- permutations(PlantGrowth$weight, control = ctrl)

dim(shuffled)

[1] 9999   30

Each row contains one shuffled version of the weight vector. Columns correspond to the original observation positions, so we can pair every row with the unchanged PlantGrowth$group vector. Here are the first few values from a few rows:

as.matrix(shuffled)[seq_len(3), seq_len(5)]

     [,1] [,2] [,3] [,4] [,5]
[1,] 6.03 4.50 4.17 5.37 5.14
[2,] 4.81 5.87 5.12 5.29 4.53
[3,] 5.14 4.53 4.50 4.17 6.03

Those are plant weights, ready to use. as.matrix() just removes the extra permutation metadata for this little look.

Now fit the model for each row in shuffled and collect the F statistics:

null_values <- apply(shuffled, 1, f_stat, group = PlantGrowth$group)

That’s our permutation distribution. The -value is the proportion of results at least as large as the observed F statistic, including the observed result itself:

tol <- sqrt(.Machine$double.eps)
p_value <- (1 + sum(null_values >= observed - tol)) /
  (length(null_values) + 1)

c(F = observed, p_value = p_value)

       F  p_value 
4.846088 0.014200

The small tolerance lets us count ties despite minor numerical differences in the model fits. Adding one to the numerator and denominator includes the observed statistic in this Monte Carlo test. We count the upper tail because larger F statistics are more extreme; there’s no need to double the result for a test of differences in either direction.

We can also have a look at the data and the permutation distribution:

old_par <- par(mfrow = c(1, 2))
boxplot(weight ~ group, data = PlantGrowth,
        xlab = "Group", ylab = "Dried plant weight",
        col = "grey85")
hist(null_values, breaks = 40, main = "",
     xlab = "F statistic after shuffling weights",
     col = "grey85", border = "white")
abline(v = observed, col = "firebrick", lwd = 2)
par(old_par)

Plant weights by group (left) and the permutation distribution of the F statistic (right). The red line marks the observed F statistic.

The observed value sits out in the right tail. For this unrestricted permutation test, there’s evidence of differences among the treatment groups. The p-value is estimated from random permutations, so it will vary with the seed and the number we request.

permutations() works with an atomic vector, rather than a whole data frame. If you give it a factor, it returns the shuffled character labels. It also follows the usual shuffleSet() rules about checking and complete enumeration, so the number of returned rows can differ from the number requested for small permutation spaces. Use the size of the returned object in subsequent calculations, as we did above.

One permutation at a time

Sometimes a loop¹ is a more convenient place to fit a model, collect several results, or keep track of progress. That’s one situation where permutator() comes in.

Calling it creates a function. Calling that function returns the next permutation of observation indices. I’ll reset the seed so that we use the same set we used above:

set.seed(42)
next_perm <- permutator(nrow(PlantGrowth), control = ctrl)

null_iter <- numeric(attr(next_perm, "nperm"))

for (b in seq_along(null_iter)) {
  i <- next_perm()
  null_iter[b] <- f_stat(PlantGrowth$weight[i], PlantGrowth$group)
}

Here i contains indices, so PlantGrowth$weight[i] does the shuffling. The iterator’s "nperm" attribute tells us how many permutations are actually in the set; we use that to allocate the results vector.

Once we’ve used all of the permutations, another call returns NULL:

next_perm()

NULL

An iterator remembers how far you’ve got. If you run the loop again with this exhausted iterator, it won’t start over; create a fresh one first.

Did the two approaches give us the same F statistics?

all.equal(null_values, null_iter)

[1] TRUE

They did. We reset the seed and used the same design, so both helpers worked through the same underlying index permutations. We can calculate the -value from the loop results in exactly the same way:

p_iter <- (1 + sum(null_iter >= observed - tol)) /
  (length(null_iter) + 1)

c(values = p_value, iterator = p_iter)

  values iterator 
  0.0142   0.0142

One detail to keep in mind: permutator() generates and stores the permutation matrix when you create the iterator. It hands out rows one at a time, but it isn’t generating them lazily or removing the need to store the index matrix.

Reusing a set

You can also pass an existing index matrix to permutator(). This is useful when you want several calculations to use exactly the same permutations, or want to revisit a saved set later:

set.seed(42)
indices <- shuffleSet(nrow(PlantGrowth), control = ctrl)

next_perm <- permutator(perms = indices)
i <- next_perm()
head(PlantGrowth$weight[i])

[1] 6.03 4.50 4.17 5.37 5.14 6.11

Creating another iterator over that matrix starts from the beginning:

replay <- permutator(perms = indices)
identical(i, replay())

[1] TRUE

Supply either perms or the generation arguments such as n, nset, and control. The matrix passed to perms must contain observation indices, such as those returned by shuffleSet(). The shuffled weights from permutations() aren’t interchangeable with that matrix.

For a single vector and a short calculation, I like permutations() with apply(). When the calculation fits more naturally in a loop, or I want to reuse a saved set of indices, permutator() is handy. Both leave the choice of permutation design with how(), so the same workflow extends to restricted designs when the data call for them. The package manual has the details, and feedback is welcome on the issue tracker.

Footnotes

Yes, I know apply is just a loop…︎

To leave a comment for the author, please follow the link and comment on their blog: From the bottom of the heap — R posts.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Continue reading: Creating a permutation test, two ways

Beyond the course: Reflections from the Bioconductor–Physalia sponsorship programme

Carlo Pecoraro — Mon, 12 Oct 2026 00:00:00 +0000

[This article was first published on Bioconductor community blog, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

What difference can a sponsored place on a bioinformatics course make? For some researchers, it means access to training they would otherwise struggle to afford. For others, it is an opportunity to explore a new field, work through unfamiliar analytical methods, or find new ways to approach an existing research question.

Through the Bioconductor–Physalia collaboration, Physalia offers two free places on eligible courses to researchers based in low- and middle-income countries (LMICs). The Bioconductor Community Manager coordinates the application process and works with the Training Committee co-chairs on selecting participants, while Bioconductor helps promote the courses through its community channels. Between November 2025 and October 2026, the programme provided 18 sponsored places across nine courses, covering topics ranging from RNA-seq and proteomics to microbiome analysis, epigenomics and network biology.

The aim is to make specialist training more accessible to researchers who might otherwise miss out on these opportunities. We asked participants what they found useful about their courses and how they planned to use what they had learned. Their responses offer a glimpse into the different ways the training connects with their work.

Making training accessible

Access was a recurring theme in the feedback. For Blessing Chioma Oselu, a PhD student and Bioinformatics Research Assistant at Covenant University in Nigeria, the opportunity to attend the Single-cell RNA-seq Analysis with R/Bioconductor course meant learning directly from experts and having structured materials to return to afterward.

“The access to expert teaching and the structured materials were the most valuable parts for me,” she wrote.

Blessing plans to continue working through the recorded sessions at her own pace. For researchers balancing training with other commitments, being able to revisit course material can be just as useful as attending the live sessions.

Another participant, from Benin, came to the Epigenomics Data Analysis course without previous experience in the subject. The step-by-step walkthroughs of methods including ChIP-seq and ATAC-seq helped make the material more manageable.

“I left the course at a clearly higher level than I started, and I feel ready to work with my own data when the time comes,” the participant reflected.

They also noted that opportunities like this are rare where they are based, and that support from Bioconductor and Physalia made a real difference.

Taking new methods back to research

Participants also described how the training related to questions they were already working on, or hoped to explore.

For Amira M. F. Abozaid, a PhD student and Assistant Lecturer at Alexandria University and E-JUST in Egypt, the RNA-seq data analysis with R/Bioconductor course provided a deeper understanding of the RNA-seq workflow, from quality control and normalization to differential expression and visualization.

Amira wrote that she had already started applying the workflows to her cancer genomics projects and sharing the skills with colleagues and students at her home institutions.

Other participants were looking at different applications. Idris Ishola Adejumobi, a Research Associate at the International Institute of Tropical Agriculture in Nigeria, described how the Epigenomics Data Analysis course broadened his understanding of how epigenetic variation can influence complex traits and complement genomic information in crop improvement. He also identified further training in the integration of epigenomic, genomic, and transcriptomic data as a useful next step for agricultural and plant-breeding research.

Rabiath Adigoun, a postdoctoral researcher at the University of Abomey-Calavi in Benin, highlighted the relevance of Orchestrating Microbiome Analysis with Bioconductor to her ongoing and future research on plant microbiomes.

“The expertise and methodological skills I have acquired will substantially strengthen my ongoing and future research on the plant microbiome in Benin and across Africa,” she wrote.

Sally Karimi Kinyua, a Research Assistant based in Kenya, attended the Network Analysis in Systems Biology with R/Bioconductor course. She hopes to apply what she learned about gene networks to research into mutations associated with malaria resistance. Sally also highlighted how the sponsored place gave her access to training that would otherwise have been out of reach because of limited resources.

These examples span different fields, but the underlying theme is similar: participants are looking for ways to connect computational methods with the biological questions they want to answer.

Learning to work with unfamiliar data

Sometimes, what matters most is finding a way into a subject that initially seems difficult.

A participant in the Statistical Methods for Metatranscriptomics course described how the practical sessions helped make complex concepts easier to follow. They particularly appreciated the explanations of foundational terminology before moving on to more technical material.

The course also introduced them to ways of analysing microbial gene expression data using computational tools and reproducible workflows. Coming from a bachelor’s degree background, they described it as one of their first deeper exposures to this area of bioinformatics.

Another participant, Binod Rayamajhee, a Research Fellow at the Kathmandu Research Institute for Biological Sciences in Nepal, highlighted the course’s combination of conceptual teaching and practical work with Bioconductor tools.

The training relates directly to his work in microbiome data analysis, and he anticipates applying the methods to study host–microbe interactions in bowel disease models.

Binod also hopes to share what he learned through workshops and one-to-one mentoring with colleagues and junior researchers, particularly those working in resource-limited settings.

Continuing the conversation

The feedback collected so far gives us a better sense of what participants value in specialist bioinformatics training. Expert teaching and well-organised materials matter, as do practical exercises that make unfamiliar methods easier to approach. Just as importantly, participants want to connect what they learn to their own research.

We would like to thank everyone who shared their reflections, along with the instructors and organisers who contributed to the courses. We’d also like to thank Laurah Nyasita Ondari for helping promote these opportunities through Bioconductor’s social media channels.

We hope to hear from participants again as they put their new skills into practice. In the meantime, their experiences offer a useful reminder of why access to training matters: it gives researchers more opportunities to learn new methods, explore their research questions, and contribute to the scientific communities around them.

For the full list of the Physalia courses and workshops, please visit www.physalia-courses.org. To find out more about the collaboration, explore the Bioconductor-Physalia course listings, including upcoming courses and how to apply for a sponsored place.

Interested in developing or teaching a Bioconductor-related course? Physalia welcomes ideas from potential instructors and is happy to discuss possible topics, content and formats. Get in touch with Physalia to explore an idea.

This article is also published on the Physalia website.

To leave a comment for the author, please follow the link and comment on their blog: Bioconductor community blog.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Continue reading: Beyond the course: Reflections from the Bioconductor–Physalia sponsorship programme

Ordered set partitions and formulas in permute

Gavin Simpson — Sun, 11 Oct 2026 22:00:00 +0000

[This article was first published on From the bottom of the heap — R posts, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

I’ve been giving permute a bit of attention recently. Two additions should make it easier to set up permutation tests: an option for ordered set partitions, and formulas for specifying plots and blocks. That’s a rather grand name for the first addition, so let’s start with what it does.

These changes are in the development version. To follow along, install it from my R-universe:

install.packages("permute", repos = c(
  "https://gavinsimpson.r-universe.dev",
  "https://cloud.r-project.org"
))

Same groups, fewer permutations

The jackal data that come with permute contain mandible lengths, in millimetres, for male and female golden jackals. They’re also used in the package vignette. Let’s ask whether mean mandible length differs between the sexes.

library("permute")
data(jackal, package = "permute")

table(jackal$Sex)

  Male Female 
    10     10

We could shuffle all the observation indices and use them to rearrange the sex labels. With the default Within(type = "free") design, the number of possible index permutations is:

numPerms(nrow(jackal), control = how(within = Within(type = "free")))

[1] 2.432902e+18

But swapping a pair of males, or a pair of females, leaves the group labels unchanged. If our statistic is a difference of means, we haven’t created a different result by doing that.

What we need are the distinct ways of assigning the observations to two labelled groups, keeping the group sizes fixed. That’s the ordered set partition. Here, ordered means the groups are distinguishable — male and female have different roles — rather than that observations have to be sorted inside each group.

The new option is type = "partition" which belongs in Plots(). I’ll use the new formula interface here too; ~ Sex says to get the grouping variable from data, to which we pass jackal:

ctrl <- how(
  plots = Plots(strata = ~ Sex, type = "partition"),
  data = jackal
)

numPerms(nrow(jackal), control = ctrl)

[1] 184756

That’s quite a reduction from all the possible permutations of the observation indices! In general, if the group sizes are , the number of distinct assignments is

The group sizes come from the data; they don’t have to be equal.

We still get observation indices back from permute, so the usual indexing workflow works:

set.seed(42)
p <- shuffle(nrow(jackal), control = ctrl)

table(jackal$Sex[p])

  Male Female 
    10     10

The labels have been shuffled, but their counts haven’t changed. Within the returned index vector, observations that originally had the same label retain their relative order. This gives us one representative for each distinct arrangement of the labels, instead of all the index permutations that would produce it.

There’s an important difference from Plots(type = "free"): that option moves whole plots as intact units. The partition option rearranges group membership. how() also selects Within(type = "none") automatically for a partition design. Asking to shuffle within these groups as well would put back the redundant orderings we’re trying to leave out, so other within-plot types aren’t allowed.

A permutation test with lm()

permute isn’t just for use with the vegan package — it was designed to be agnostic to the underlying model or test statistic. As an illustration, let’s fit a linear model to the jackal data:

fit <- lm(Length ~ Sex, data = jackal)
coef(fit)

(Intercept)   SexFemale 
      113.4        -4.8

The intercept is the mean for males, the reference level. The second coefficient is the female mean minus the male mean: females have mandibles that are, on average, 4.8 mm shorter in these data.

For a permutation test, we’ll keep the lengths where they are and shuffle the sex labels. We then refit the model for each assignment and save that same coefficient. Under the null hypothesis, we’re treating the labels as exchangeable among observations; this needs to be a reasonable assumption for the data we’re analysing.

set.seed(42)
perms <- shuffleSet(nrow(jackal), nset = 999, control = ctrl)

observed <- unname(coef(fit)[2])
null <- apply(perms, 1, function(i) {
  permuted <- jackal
  permuted$Sex <- jackal$Sex[i]
  unname(coef(lm(Length ~ Sex, data = permuted))[2])
})

shuffleSet() gives us a matrix with one assignment per row. We’re sampling a manageable set here, rather than fitting a model for every possible assignment.

Here’s what those coefficients look like, with the observed difference marked in red:

hist(null, breaks = "FD", main = "",
     xlab = "Female minus male mean mandible length (mm)",
     col = "grey85", border = "white")
abline(v = observed, col = "firebrick", lwd = 2)

Permutation distribution of the female-minus-male mean mandible length. The red line marks the observed difference.

For a two-sided test, we count differences at least as large as the observed one in absolute value:

tol <- sqrt(.Machine$double.eps)
p_value <- (1 + sum(abs(null) >= abs(observed) - tol)) /
  (nrow(perms) + 1)

c(difference_mm = observed, p_value = p_value)

difference_mm       p_value 
       -4.800         0.005

The small tolerance includes ties that might otherwise be missed because of floating-point arithmetic in the model fits. Adding one to the numerator and denominator includes the observed result in this Monte Carlo test. The p-value is an estimate based on the sampled assignments, so it will vary with the random seed and the number requested.

If you don’t need to build a design explicitly, shufflePartition(jackal$Sex) and shufflePartitionSet(jackal$Sex, nset = 999) are handy shortcuts. There’s also allPartitions() for complete enumeration when the assignment space is small enough. As with allPerms(), it leaves out the observed assignment by default; use control = how(observed = TRUE) to include it.

Less typing with formulas

The formula interface works for both plot strata and blocks. To show the two together, let’s use the dune meadow data from vegan. dune contains the species data and dune.env contains the matching environmental variables.

Suppose we want to look at the association between vegetation and soil A1 horizon thickness, allowing comparisons only among sites with the same management and land use. We can express that restriction directly:

library("vegan")
data(dune, package = "vegan")
data(dune.env, package = "vegan")

formula_ctrl <- how(
  plots = Plots(strata = ~ Management / Use),
  blocks = ~ Management,
  data = dune.env
)

~ Management / Use creates a grouping factor for the observed combinations of management and land use. Plots() defaults to type = "none", while Within() defaults to "free", so sites are shuffled within those combinations. The blocks say that sites must also stay within their management category. The plot restriction already implies that here, but writing out both makes the hierarchy clear.

These restrictions assume sites are exchangeable within each management–use combination; grouping them this way doesn’t account for any remaining spatial dependence.

We can check where an actual shuffle sends the observations:

set.seed(42)
i <- shuffle(nrow(dune.env), control = formula_ctrl)

all(getBlocks(formula_ctrl)[i] == getBlocks(formula_ctrl))

[1] TRUE

all(getStrata(formula_ctrl)[i] == getStrata(formula_ctrl))

[1] TRUE

Both the blocks and plot groups are respected. We can use the design directly in a vegan analysis. Here I partial out the management–use groups in an RDA, then test the additional association with A1 horizon thickness:

env <- dune.env
env$group <- getStrata(formula_ctrl)
ord <- rda(dune ~ A1 + Condition(group), data = env)

set.seed(42)
anova(ord, permutations = formula_ctrl)

Permutation test for rda under reduced model
Blocks:  Management 
Plots: Management/Use, plot permutation: none
Permutation: free
Number of permutations: 199

Model: rda(formula = dune ~ A1 + Condition(group), data = env)
         Df Variance      F Pr(>F)
Model     1   4.3492 1.3265  0.225
Residual  8  26.2298

There’s little evidence of an additional association in this particular analysis. Some of the groups contain only one site, so those sites can’t move at all under this design. Restrictions matter, both for what question we’re asking and for how much information is available to answer it.

A couple of details are worth keeping in mind when writing these formulas. They describe grouping factors, not model terms: ~ Management + Use, ~ Management:Use, and ~ Management / Use all produce groups from the combinations of those variables. The slash doesn’t create an extra level in permute’s hierarchy.

Also, when Plots() is inside how(), supply data to how(), as above. If you’re constructing Plots() on its own, it can take its own data argument:

plot_groups <- Plots(strata = ~ Management / Use, data = dune.env)

The existing factor interface still works, so there’s no need to rewrite old code. But I find the formulas easier to read, especially when grouping variables need to be combined. You can find the details in the design documentation, and I’d welcome feedback on either addition through the issue tracker.

There are also new helpers for working with the resulting permutations: permutations() returns shuffled values directly, while permutator() lets you step through a set of permutations. I’ll show those in a follow-up post.

To leave a comment for the author, please follow the link and comment on their blog: From the bottom of the heap — R posts.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Continue reading: Ordered set partitions and formulas in permute

Going back to the basics: the correlation coefficient

Andrea Onofri — Thu, 08 Oct 2026 22:00:00 +0000

[This article was first published on Statforbiology, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

In statistics, dependence or association is any statistical relationship, whether causal or not, between two random variables or bivariate data. It is often measured by the Pearson’s correlation coefficient:

Other measures of correlation can be thought of, such as the Spearman rank correlation coefficient or Kendall rank correlation coefficient.

Assumptions for the Pearson Correlation Coefficient

The Pearson correlation coefficients makes a few assumptions, which should be carefully checked.

Interval-level measurement: Both variables should be measured on a quantitative scale.
Random sampling: Each subject in the sample should contribute one value on and one value on . The values for both variables should represent a random sample drawn from the population of interest.
Linearity: The relationship between and should be linear.
Bivariate normal distribution: This means that (i) the values of follow a normal distribution at each level of , and (ii) the values of follow a normal distribution at each level of .

Hypothesis testing

It is possible to test whether against the alternative $ r $. The test is based on the test statistic:

which follows a Student’s distribution with degrees of freedom.

Let’s take the two variables cyl and mpg from the mtcars data frame. For these variables, the Pearson’s correlation index is (please note, that the Spearman’s correlation doefficient might be preferred for these variables, although I neglected this for the sake of simplicity):

r <- cor(mtcars$cyl, mtcars$gear)
r
## [1] -0.4926866

The T statistic is:

T <- r * sqrt(32 - 2) / sqrt(1 - r^2)
T
## [1] -3.101051

The p-value for the null hypothesis is:

2 * pt(T, 30)
## [1] 0.004173297

which is clearly highly significant. The null can be rejected.

When performing this test, it is important to note that the individuals from whom the pairs of measurements were taken must be independent. If they are not, the ‘t-test’ is invalid. I discuss this aspect in another post on my blog.

Correlation in R

We have already seen that we can use the standard cor(matrix, method = …) function. To obtain significance levels as well, we can use the rcorr() function from the ‘Hmisc’ package:

# Correlations with significance levels
library(Hmisc)
corr2 <- rcorr(as.matrix(mtcars), type="pearson")
print(corr2$r, digits = 2)
##        mpg   cyl  disp    hp   drat    wt   qsec    vs     am  gear   carb
## mpg   1.00 -0.85 -0.85 -0.78  0.681 -0.87  0.419  0.66  0.600  0.48 -0.551
## cyl  -0.85  1.00  0.90  0.83 -0.700  0.78 -0.591 -0.81 -0.523 -0.49  0.527
## disp -0.85  0.90  1.00  0.79 -0.710  0.89 -0.434 -0.71 -0.591 -0.56  0.395
## hp   -0.78  0.83  0.79  1.00 -0.449  0.66 -0.708 -0.72 -0.243 -0.13  0.750
## drat  0.68 -0.70 -0.71 -0.45  1.000 -0.71  0.091  0.44  0.713  0.70 -0.091
## wt   -0.87  0.78  0.89  0.66 -0.712  1.00 -0.175 -0.55 -0.692 -0.58  0.428
## qsec  0.42 -0.59 -0.43 -0.71  0.091 -0.17  1.000  0.74 -0.230 -0.21 -0.656
## vs    0.66 -0.81 -0.71 -0.72  0.440 -0.55  0.745  1.00  0.168  0.21 -0.570
## am    0.60 -0.52 -0.59 -0.24  0.713 -0.69 -0.230  0.17  1.000  0.79  0.058
## gear  0.48 -0.49 -0.56 -0.13  0.700 -0.58 -0.213  0.21  0.794  1.00  0.274
## carb -0.55  0.53  0.39  0.75 -0.091  0.43 -0.656 -0.57  0.058  0.27  1.000
print(corr2$P, digits = 2)
##          mpg     cyl    disp      hp    drat      wt    qsec      vs      am
## mpg       NA 6.1e-10 9.4e-10 1.8e-07 1.8e-05 1.3e-10 1.7e-02 3.4e-05 2.9e-04
## cyl  6.1e-10      NA 1.8e-12 3.5e-09 8.2e-06 1.2e-07 3.7e-04 1.8e-08 2.2e-03
## disp 9.4e-10 1.8e-12      NA 7.1e-08 5.3e-06 1.2e-11 1.3e-02 5.2e-06 3.7e-04
## hp   1.8e-07 3.5e-09 7.1e-08      NA 1.0e-02 4.1e-05 5.8e-06 2.9e-06 1.8e-01
## drat 1.8e-05 8.2e-06 5.3e-06 1.0e-02      NA 4.8e-06 6.2e-01 1.2e-02 4.7e-06
## wt   1.3e-10 1.2e-07 1.2e-11 4.1e-05 4.8e-06      NA 3.4e-01 9.8e-04 1.1e-05
## qsec 1.7e-02 3.7e-04 1.3e-02 5.8e-06 6.2e-01 3.4e-01      NA 1.0e-06 2.1e-01
## vs   3.4e-05 1.8e-08 5.2e-06 2.9e-06 1.2e-02 9.8e-04 1.0e-06      NA 3.6e-01
## am   2.9e-04 2.2e-03 3.7e-04 1.8e-01 4.7e-06 1.1e-05 2.1e-01 3.6e-01      NA
## gear 5.4e-03 4.2e-03 9.6e-04 4.9e-01 8.4e-06 4.6e-04 2.4e-01 2.6e-01 5.8e-08
## carb 1.1e-03 1.9e-03 2.5e-02 7.8e-07 6.2e-01 1.5e-02 4.5e-05 6.7e-04 7.5e-01
##         gear    carb
## mpg  5.4e-03 1.1e-03
## cyl  4.2e-03 1.9e-03
## disp 9.6e-04 2.5e-02
## hp   4.9e-01 7.8e-07
## drat 8.4e-06 6.2e-01
## wt   4.6e-04 1.5e-02
## qsec 2.4e-01 4.5e-05
## vs   2.6e-01 6.7e-04
## am   5.8e-08 7.5e-01
## gear      NA 1.3e-01
## carb 1.3e-01      NA

We can also use these functions with two separate matrices to obtain the correlations between each column of one matrix and each column of the other:

# Correlation matrix from mtcars
x <- mtcars[1:3]
y <- mtcars[4:6]
cor(x, y)
##              hp       drat         wt
## mpg  -0.7761684  0.6811719 -0.8676594
## cyl   0.8324475 -0.6999381  0.7824958
## disp  0.7909486 -0.7102139  0.8879799

Relationship to slope in linear regression

The correlation coefficient and the slope in linear regression share some similarities, as both describe how changes with . However, in correlation, we deal with two random variables, whereas in simple linear regression we treat as fixed and as a random variable, with regarded as a function of (and not the other way around).

Despite their distinct interpretations, it is useful to examine the algebraic relationship between the correlation coefficient and the regression slope. Let’s simulate a dataset with two variables drawn from a multivariate normal distribution, with means of 10 and 2, respectively, and a variance-covariance matrix of:

library(MASS)
cov <- matrix(c(2.20, 0.48, 0.48, 0.20), 2, 2)
cov
##      [,1] [,2]
## [1,] 2.20 0.48
## [2,] 0.48 0.20

We use the mvrnomr() function to generate the dataset.

set.seed(1234)
dataset <- data.frame( mvrnorm(n=10, mu = c(10, 2), Sigma = cov) )
names(dataset) <- c("X", "Y")
dataset
##            X        Y
## 1  11.756647 2.547203
## 2   9.522180 2.199740
## 3   8.341254 1.862362
## 4  13.480005 2.772031
## 5   9.428296 1.573435
## 6   9.242788 1.861756
## 7  10.817449 2.343918
## 8  10.749047 2.451999
## 9  10.780400 2.436263
## 10 11.480301 1.590436

The correlation coefficient and slope are as follows:

r <- with(dataset, cor(X, Y))
b1 <- coef( lm(Y ~ X, data=dataset) )[2]
r
## [1] 0.6372927
b1
##         X 
## 0.1785312

The formula for the slope is:

From this, we see that:

and:

Indeed:

sigmaX <- with(dataset, sd(X) )
sigmaY <- with(dataset, sd(Y) )
b1 * sigmaX / sigmaY 
##         X 
## 0.6372927
r * sigmaY / sigmaX
## [1] 0.1785312

It is also easy to see that the correlation coefficient is equivalent to the regression slope of standardized against standardised :

Yst <- with(dataset, scale(Y, scale=T) )
summary( lm(Yst ~ I(scale(X, scale = T) ), data = dataset) )
## 
## Call:
## lm(formula = Yst ~ I(scale(X, scale = T)), data = dataset)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.082006 -0.067143 -0.036850  0.009214  0.237923 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)  
## (Intercept)            1.065e-18  3.478e-02   0.000   1.0000  
## I(scale(X, scale = T)) 1.785e-01  7.633e-02   2.339   0.0475 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.11 on 8 degrees of freedom
## Multiple R-squared:  0.4061, Adjusted R-squared:  0.3319 
## F-statistic: 5.471 on 1 and 8 DF,  p-value: 0.04748

Intra-class correlation (ICC)

This is a related yet distinct concept. While Pearson’s correlation coefficient is used for bivariate data—where two different variables are observed on the same independent units—intraclass correlation arises when the same variable is repeatedly observed on the same units (e.g., across different times or spaces), so that each unit acts as a grouping factor that creates clusters of measurements.

Let’s illustrate this with an example. Consider the following data table:

#|echo: false
m1 <- c(10, 12, 13, 16)
m2 <- c(21, 20, 17, 19)
dtab <- cbind(m1, m2)
row.names(dtab) <- paste0("g", 1:4)
knitr::kable(dtab)

	m1	m2
g1	10	21
g2	12	20
g3	13	17
g4	16	19

When g1, g2, g3, and g4 represent four independent units (e.g., four maize plants) on which two variables, m1 and m2 (e.g., weight and height), were measured, we are in the realm of Pearson’s correlation. On the other hand, when g1, g2, g3, and g4 are entities on which multiple measurements of the same variable (e.g., yield) are taken, we enter the realm of intraclass correlation. In this second scenario, we are not assessing joint variability between two different variables; rather, we want to determine whether belonging to the same group (‘proximity’) makes observations more similar to each other than observations from different groups.

Figure 1: Confronto concettuale tra Pearson e ICC

The most common intraclass correlation coefficient is:

where is the between-group variance component and is the within-group variance component (or, more clearly, the variance of observations within a group). For balanced data, those two components can be obtained from one-way ANOVA, considering that there are four ‘replicates’ per group:

y <- c(10, 12, 13, 16, 21, 20, 17, 19)
group <- factor( rep(c("A", "B"), each = 4) )
nmeasures <- 4
mod <- lm(y ~ group)
anova(mod)
## Analysis of Variance Table
## 
## Response: y
##           Df Sum Sq Mean Sq F value   Pr(>F)   
## group      1   84.5  84.500  18.436 0.005128 **
## Residuals  6   27.5   4.583                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
sigma2_alpha <- (84.5 - 4.583)/nmeasures
sigma2_epsilon <- 4.583
icc <- sigma2_alpha/(sigma2_alpha + sigma2_epsilon)
icc
## [1] 0.8134129

In general, if subjects in one group are very similar (small ), while different groups vary widely (high ), the ICC is high. Conversely, if measurements are just as similar across different groups as they are within the same group, the ICC is low. In theory, the ICC cannot be negative, as variance components are expected to be positive. In practice, however, the ICC calculated via the method of moments shown above can turn out negative when the group mean square is lower than the residual mean square from the ANOVA. Such negative values are usually constrained to zero or interpreted as indicating an absence of positive clustering; in both cases, intraclass correlation remains a meaningful tool.

ICC is crucial and is most often studied using linear mixed models. But… that is a story for another day…

Thanks for reading, and happy coding!

(P.S. If you enjoyed this post, check out my book “Field Research Methods in Agriculture: An Introduction with R” via the link below!)

Prof. Andrea Onofri
Department of Agricultural, Food and Environmental Sciences
University of Perugia (Italy)
Send comments to: andrea.onofri@unipg.it

This post was originally published on 2019-02-07

To leave a comment for the author, please follow the link and comment on their blog: Statforbiology.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Continue reading: Going back to the basics: the correlation coefficient

Validating data analyses written by coding agents

Seascapemodels — Thu, 08 Oct 2026 13:00:00 +0000

[This article was first published on Seascapemodels, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Claude Code (and many other agents) has been getting so good at ecological analysis that I find I’m increasingly tempted not to directly review the code it writes. In fact, it finds more mistakes in my code than I find in its code these days.

So I need a way to check the validity of the code it writes. Here’s the strategy for this I’ve been developing. Its a work in progress.

Specification

Start with a clear written specification for the data analysis or modelling. This should include equations for the modelling. This ensures I understand what I want the analysis to do (even if its not what I get initially).

Claude is very helpful at this stage to review it for logical errors.

Adversarial code review

This is a no brainer. Use the agent, or ideally another provider’s agent, to review the code for errors.

Data processing checks

Instruct Claude to check all data processing steps and fail loudly if inconsistencies are found. That’s probably sufficient these days, but if you want to be thorough you can think about what you need at each stage and tell it to check for that.

For example, make sure the final dataframe used in the model has N rows, or that the number of unique survey IDs == the number of rows in the dataframe.

There are various R packages to support this like pointblank and validate.

Testing against simulated data with known truth

Often in ecology we are fitting generative models, that is, you can simulate data right out of the model you are fitting. So invert your specification of the model to be a specification that simulates data with a known model. Then run the simulated data through your workflow to see if it recovers the known parameters.

If you are doing probabilistic inference you can also check other common verification statistics like coverage (the 95% CIs should cover the point estimate in 95% of randomized datasets).

Create the data simulation in a separate Claude and R environment, to ensure independence. Make sure the AI memory storage options are turned off, or set to incognito mode, to ensure no leakage of information through memory files too.

There are also many data simulation functions as R packages or within R packages. For instance, many advanced modelling packages also provide a simulation function. You can use these to further reduce the chances of errors in your simulation step

Other simulation tests

You can also do checks where you muddle the real data directly and run it through your workflow. For instance, you should get the same result if you reshuffle rows of your data-frame randomly.

If you reshuffle values in the covariate columns independently of the response, then the effects of covariate should be 0 on average. If your workflow is about model selection, then enforcing shared rank ordering of a covariate and the response you should (almost certainly) find that covariate pops up in every model selected by the AIC.

Test multiple methods

In the simplest case, the same model fitted with different implementations should return the same results. For example compare results from lmer lme and glmmTMB. You can get claude to fire up sub-agents to write a script for each one independently.

In more complex cases it can still be useful to do fit different models and compare results, even in the models have slightly different meanings. This may reveal errors, but it may also reveal places where your model assumptions are not robust. For instance, do different model selection methods result in the same model?

What’s next?

Part of the problem is that once you start using an agent it will write very complex R code. Claude Code has this obsession with making absolutely generalizable scripts. For instance, it commonly writes functions just to identify the root file path for a key data file.

The other part of its complexity is adding numerous print statements to the R scripts, so it can see intermediate results in the terminal output. I normally do these interactively so they stay out of my scripts.

All this extra code means it becomes even more tedious to review its work. Its like a double edged sword, the agent is powerful, but it makes more review work.

So its clear to me that we soon won’t bother reading the code they write. That means we need to adapt how we approach coding and data analysis, and especially ensure we have checks in place that we understand the code we are producing and understand its correct.

This list is my start at that.

Send a comment if you have other suggestions for ways to test code written by agents.

To leave a comment for the author, please follow the link and comment on their blog: Seascapemodels.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Continue reading: Validating data analyses written by coding agents

The Pit of Success, and What a Data Science Language Would Need to Make It

Econometrics and Free Software — Thu, 08 Oct 2026 00:00:00 +0000

[This article was first published on Econometrics and Free Software, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Rico Mariani, a long-time performance engineer at Microsoft, coined the idea of the pit of success. A system has a pit of success when the natural, lazy way of using it leads to good outcomes. You don’t need heroics or perfect memory. You fall into the right result, and climbing out to do something wrong takes deliberate effort.

I came across this framing in a talk that I was recently recommended and I further recommend it to anyone interested in the topic: Functional architecture – The pits of success – Mark Seemann. It put a name to something I had been circling around for years.

I had been arguing a version of the same point long before I knew the term. In my 2022 post titled Functional programming explains why containerization is needed for reproducibility, I argued that a pipeline is only as reproducible as the functions it is built from, and that even a pure function hides inputs such as the R version, the package versions, and the operating system. Containers are one way to make those hidden inputs explicit. In another 2022 post on learnings from functional programming, I made a related case for writing functions that fail early and never reach into the global environment. Both posts were, in a sense, attempts to build a pit of success out of habits and tooling. Neither was framed that way, and I think the framing is more useful.

Pits of failure are everywhere in data science

Most of the tools we use for analysis are pits of failure. You can be meticulous and still end up with an analysis that cannot be rerun, because the meticulousness is entirely up to you.

Think about how a typical script behaves. It reads a file by path, and nothing says that the path is an input. It uses a package, and nothing records which version. It reassigns a variable halfway through, so the value at the bottom of the file depends on what happened at the top. It converts a column from text to numbers without complaint, and quietly turns the bad entries into missing values. Each of these is convenient in the moment. Together, they make the result depend on things nobody wrote down.

Good practice is then a set of instructions: pin your dependencies, avoid globals, check your inputs, document your steps. These are road signs. They tell you where to go, and they leave the driving to you. Under deadline pressure, road signs get ignored, and not out of malice. Only a speed bump would make you drive more responsibly.

What makes a pit of success

If the goal is to make the disciplined path the easy one, we need to ask what the language itself would have to do. If reproducibility depends on remembering a checklist, then reproducibility is a convention. If the language makes the checklist unnecessary, reproducibility becomes a property of the program. From my own experience of building and using R packages, and from thinking about where analyses usually go wrong, I think a language designed for this purpose would need at least the points raised in this blog post from earlier this year, but let me rephrase them here again.

The program should be a flat text file. Notebooks have become popular in data science, but I think they are a serious mistake. They mix code, output, and state in a format that is hard to diff, review, and reuse, and they encourage the hidden-state problems that make analyses irreproducible. Notebooks suffer from exactly the same problem that gets blamed on Excel. The main argument for notebooks is that they make exploration easy. But that says more about the shortcomings of the alternatives than about any real strength of notebooks. Environments like RStudio and Positron (or Emacs) already make interactive work with data natural, and they do it while keeping the program in a plain text file that can be run, tested, and versioned like any other code.

Reproducibility should be the default, not an option. If pinning an environment requires a separate file that someone has to remember to update, it will drift. A language that declares its runtime and dependencies as part of the program, and builds them the same way every time, removes a decision from the analyst. You cannot opt out of it by forgetting.

Inputs should be declared, not discovered. A program should not be able to read a file or an environment variable that it does not name. When every dependency is visible in the code, the tooling can see the whole picture: which step depends on which, and what changes when an input changes. Undeclared inputs become a build error rather than a silent difference between two machines. This is the same point as the hidden-inputs argument from my containerization post, moved from the operating system into the language.

Hidden state should be hard to create. Immutable bindings, explicit reassignment when you really need it, and no implicit global state make a result depend only on its inputs. This is also what makes caching safe. If a step is a pure function of its declared inputs, the system can reuse its previous output with confidence, and nobody has to wonder whether a cached result is stale.

Types and schemas should fail loudly. Silent coercion is one of the most common ways data goes wrong without anyone noticing. A language built for analysis should refuse to turn text into numbers without being asked, should make missing values explicit rather than inventing them, and should check that a column a downstream step needs actually exists before anything runs.

An analysis should be laid out as a pipeline. Traditional scripts tend to morph into spaghetti very quickly: steps depend on each other in ways that are hard to see, and a change in one place can silently affect results three sections later. Writing the analysis as an explicit pipeline, with each step a named stage that takes declared inputs and produces declared outputs, makes it much easier to debug a failing step, reuse results that have not changed, and hand the work over to someone (or something) else.

The whole graph should be checkable before it runs. If the pipeline is a first-class object that the language can inspect, you can catch missing columns, broken references, and cycles in milliseconds. Waiting three hours for a simulation to discover a typo in a column name is a failure of design, not of the analyst.

Errors should be values. When a step fails, the rest of the pipeline should not necessarily stop. An error that is captured as data can be inspected, reported, and handled, and independent branches can keep going. This turns debugging from reading a stack trace from code you didn’t write into reading a result that tells you where the problem is.

Interoperability should not require glue code. Real analysis mixes R, Python, Julia, SQL, and whatever comes next. Every hand-written CSV export between two of them is a place where types get lost and assumptions go unrecorded. A shared, typed interchange format between runtimes removes a whole category of errors and a lot of boilerplate.

Feedback should be structured. Error messages that are only meant for a human reading a terminal are hard for tools to act on. If the language reports problems in a machine-readable form, the same checks serve analysts, continuous integration, and increasingly, AI coding assistants that write much of the code.

What makes this set of properties attractive is that the same features that make an analysis reproducible also make it easier for an LLM to work with. An LLM is much better at modifying a program when its inputs are explicit, its dependencies are named, its state is visible, and its errors are structured. In that sense, good language design doesn’t just put guardrails around the human analyst. It gives the machine a better representation of what the analyst is trying to do.

This suggests an interesting convergence: the language we would want to make data science safer for humans is also close to the language we would want to make data science easier for AI to operate on. A language built around these principles could be a pit of success for both the human analyst and the AI assistant helping them. Instead of pulling in different directions, the two would be working with the same constraints.

Defaults, constraints, and choice architecture

There is also a useful way to think about this from economics. We often distinguish between telling people what they ought to do and changing the rules and incentives under which they make decisions. A recommendation such as “pin your dependencies” or “don’t rely on global state” is a rule in the first sense: it asks the analyst to remember and follow a good practice. Defaults, constraints, and institutions work differently. They change what is easy to do, what is costly to do, and what is possible to do at all.

This is closely related to the idea of choice architecture in behavioural economics, and to the broader observation that the design of an environment can shape behaviour without requiring people to consciously follow a set of instructions. In software, the architecture is unusually literal: a language can make some actions easy, others difficult, and some impossible to express without making the programmer explicitly acknowledge what they are doing.

That suggests a useful question for data science. Instead of asking how we can get analysts to consistently follow reproducibility best practices, what if we designed the language so that following those practices was the easiest way to write a program?

Design tools that make falling into a pit of success easy

But have you tried asking someone who does none of this to adopt these practices? People may agree with every point, but applying them consistently takes effort and frankly brings very little immediate reward. Investing early in good practice makes life much easier later, but “later” is hard to sell to someone with a deadline yesterday.

So what is the alternative? Rather than persuading people to jump into the pit, we can build tools that make falling into it the default. I’ve been working on a domain-specific language called T that aims at exactly this. The idea is that by simply using it, you end up in a pit of success. I won’t go into the details here, since you can read more about it here and here. What I want to discuss is why LLMs make this kind of tool practical to build.

Historically, building a language or a tool was a large project that only a few people could justify. LLMs lower that cost dramatically. Someone with experience in a domain can now build a small tool that solves one specific problem well, without a team or a long development cycle. Taste, knowing which problems matter and which habits are worth enforcing, is exactly what a domain expert brings and what an LLM does not supply on its own.

I don’t know yet how widely T will be adopted. But I think it is a good example of the pit-of-success idea applied to data science, and the broader point is that the tools themselves can carry the discipline, so that individual analysts don’t have to. The future will tell whether that bet pays off.

Where the pit needs walls and an exit

A pit of success can become a trap if it is the only way to work. Exploration is messy by nature. You try a model, plot something odd, throw away three approaches before finding the one that works. A language that forbids all of that makes people go back to their old tools, usually a notebook on the side, and then the discipline lives in two places at once.

So the language needs a deliberate place for exploration that does not contaminate production. Exploration should be easy and interactive, with the results promoted into the strict pipeline once they stabilise. The pit should catch you when you go into production, not when you’re trying to think. This is the part I find hardest to get right, and I suspect any language that ignores it will be abandoned quickly. For now, the way this works in T is that you have much more freedom to try things from its REPL. But I suspect that much of this exploratory work will be LLM-driven in the future. Interactive data exploration may move from “playing around” with the data to asking an LLM to do things for you and show you the results.

If that happens, interactive exploration may become less important as a feature of the language itself. The important requirement would instead be that the language provides a good interface for an LLM to explore the data and then turn what it discovers into a reproducible pipeline.

I think much of the answer lies in the tools we already use. Editors such as Emacs and Positron or VS Code already give analysts a place to explore interactively while writing plain text files that can later be run as a pipeline. What a new language needs is an editor integration that makes this transition smooth: completions, diagnostics, and navigation working in the same window as the exploratory code. That integration is a solvable engineering task, and LLMs are now quite capable of producing the plugins, language servers, or tree-sitter grammars it requires. One question is whether the language makes the move from exploration to production feel like tidying up rather than rewriting from scratch, which is, again, the typical life-cycle of an analysis. Quick and dirty exploration in a Notebook, clean rewrite in a script (sometimes even in another language). T solves this by forcing the analysis to be a pipeline that runs in a reproducible environment by design.

What the pit cannot do

A language, even one supported by an LLM, can make reproducibility cheap. It cannot make the analysis correct. A misspecified regression run reproducibly is still misspecified, and catching that remains the work of a statistician and a reviewer. Nor can a language change the culture around it. People need to believe that rebuilding last year’s figures is worth the time, and someone needs to ask whether they still can.

The pit of success makes the right thing easy. It does not decide what the right thing is, and it does not make anyone want to jump in. Those are still human jobs.

This is roughly the bet behind T, though I think the general ideas matter more than any one implementation. If you have thought about what a data science language should refuse to let you do, I would be interested to hear what you would put on that list.

To leave a comment for the author, please follow the link and comment on their blog: Econometrics and Free Software.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Continue reading: The Pit of Success, and What a Data Science Language Would Need to Make It

Bayesian workflow [book review]

xi'an — Wed, 07 Oct 2026 22:26:44 +0000

[This article was first published on R – Xi'an's Og, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

“This original, thought-provoking, and transforming, book is much much more than an implementation manual for Bayesian Data Analysis, even though it shares almost the same perspective. (The first sentence of the book states that the authors’ `conceptions of statistical practice, and of Bayesian statistics, have changed over the years’.) By providing a modus vivendi for undertaking Bayesian modelling from scratch in realistic settings where models are not magicked out of the blue, the authors explicit and rationalise the many steps required by such a bottom-up modelling protocol (`not a checklist, not a cookbook’, and not a flowchart!) in real situations. The contents read very well and very smoothly, with a seamless conjunction of intuition, modelling advices, computational details, and comparison tools. While unsurprisingly Bayesian, the perspective adopted therein remains both open and inclusive, with a welcome humility about the limitations and challenges of Bayesian workflows. This book should thus appeal to and profit a wide variety of readers, as providing guidance through an extensive collection of highly detailed examples, with shared code and exercises.”

This book proposes a modus vivendi for Bayesian modelling in applied, realistic Bayesian analysis, where models are not magicked out of the blue. It thus emphases iterative model building, model checking, computational troubleshooting, and simulated-data experimentation, filling a gap that looks glaring in retrospect. It particularly targets users and developers of Stan, with code excerpts in R and Stan. It consists of four parts:

background on Bayesian methods and computational tools;
the Bayesian workflow proper, namely building a statistical model from its components, together with its assessment tools;
the computational aspects of fitting models, diagnosing convergence and assessing calibration;
case studies.

I was eagerly waiting for the book, as I knew Andrew, Aki, and Richard had been working on it for a few years. (The quote above is the blurb I wrote upon request from the publisher.)

The tenets of BaWoFlo—if I may resort to this acronym!—are (i) fitting multiple models, (ii) applying methods repeatedly, and (iii) resorting to simulated-data experiments, which should not come as a surprise to readers of BDA. As noted in the introduction, the protocol exposed therein can also benefit non-Bayesian experimenters. This agrees with the highly moderate, “M-open”, agnostic approach to Bayesianism adopted by the authors (“there is no safe haven”). I also welcome and share their humble perspective about the limitations and challenges of Bayesian workflows.

Examples are treated in full detail, with successive modelling and computational choices profusely commented, which is a big plus for such a practical book. This starts as early as Chapter 4, with a multiple-choice exam example. Indeed, there cannot be general principles or a generic theory that would make the approach foolproof. See, e.g., “A data model is not just a ‘likelihood’” (p.70), as when the data model is not fully generative. I very much liked the section on choosing priors (5.6), and the very rich graphs (see, e.g., Chapter 8) for assessing the impact of prior and likelihood, as well as for predictive checks. In coherent continuation of the authors’ earlier work, the book advocates LOO methods and model stacking rather than model averaging. (With a surprisingly anti-Ockham perspective in Section 9.7.)

The MCMC coverage is unsurprising, with $\hat R$ at the forefront. Chapter 12, on using fast experiments to detect fitting or computational issues, is very nice. The book builds on the immense corpus of work achieved by the authors over the decades (for the most senior ones!). By contrast, the chapter on approximate solutions (13) is way too short, and the same goes for those on calibration and software development.

The book is very US-centric, unsurprisingly given Andrew’s focus on political science. Some sections are reminiscent of Andrew’s blog entries (or the opposite). The (football) World Cup example was initiated when Andrew was in France, during the 2014 World Cup, and as a result (?) the names of the teams are in French! One chapter also reanalyses the birthdate data displayed on the cover of BDA.

Mileage varies on the applied chapters, depending on the example. A dog chapter is followed by a cat chapter! Not that the (stat)dog experiment was in any way enjoyable, especially for the dogs. Maybe the cats were running it! And then come chapters on roaches and sharks. There is also a frightening flowchart (Fig. 2.1)! And the book ends with an appendix on going through BDA to better understand BaWoFlo

[The usual disclaimer applies, namely that this review is likely to appear later in CHANCE, in my book reviews column.]

To leave a comment for the author, please follow the link and comment on their blog: R – Xi'an's Og.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Continue reading: Bayesian workflow [book review]

A Collaboration Fest for research software interoperability

Kevin Rue-Albrecht — Wed, 07 Oct 2026 00:00:00 +0000

[This article was first published on Bioconductor community blog, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

On 15-17 September 2026, members of the Bioconductor, Galaxy and nf-core communities came together at The Open University in Milton Keynes for the BioFAIR Data to Discovery and BioFAIR Collaboration Fest (CoFest) joint event.

The three-day event combined hands-on training with collaborative work on research software interoperability. It was jointly organised by BioFAIR Fellow Marisa Loach and BioFAIR Pathfinder project lead and Institute for Research Software 2026 Fellow Kevin Rue-Albrecht.

The CoFest brought together contributors based in the UK, Ireland, Germany, India and the USA, with both in-person and remote participation.

Marisa led the Data to Discovery training, introducing participants to single-cell analysis using the European instance of the Galaxy platform (https://usegalaxy.eu/), while Kevin led the hybrid CoFest focused on making methods developed as Bioconductor R packages more readily accessible through infrastructures such as Galaxy and nf-core.

A Galaxy community recap of the event also covers the training and CoFest from the Galaxy perspective.

The challenge: making methods reusable across ecosystems

There are numerous communities in the bioinformatics space focusing on different aspects like code development, making non-code GUIs, creating tutorials and documentation, and more. Some of these are Bioconductor, Galaxy and nf-core, which play complementary roles in bioinformatics. Bioconductor provides a large ecosystem of R packages implementing analysis methods; Galaxy provides an accessible environment for running and sharing reproducible analyses; and nf-core provides a community framework for portable Nextflow pipelines.

The challenge is that making the same method available across these environments often requires separate wrappers, metadata and maintenance. The BioFAIR Pathfinder project behind the event is exploring ways to reduce that duplication, including tooling for adapting Bioconductor methods as Galaxy tools and nf-core modules.

A major focus of the CoFest was BiocJobs, an approach being developed that allows analysis tasks to be described alongside Bioconductor packages and used to generate representations for different workflow systems.

Putting the idea into practice

Rather than only discussing interoperability, participants applied this approach to an exemplar single-cell analysis workflow.

The group mapped parts of that workflow to reusable package-level ‘BiocJobs’ and tested the path from descriptions maintained alongside Bioconductor packages to generated workflow components. Work covered steps including data import, quality control, normalisation, dimensionality reduction and marker analysis.

This gave contributors a concrete way to test the process, improve documentation, and identify what would be needed to make it easier for others to contribute.

On the final day, we captured the workflow we had been working towards in a single diagram: from a commit to a Bioconductor package repository, through the automatic building and publishing of tool wrappers, to their use in Galaxy, nf-core and other workflow environments.

Three days were not enough to fully implement the workflow captured in the diagram – nor was that the goal. Instead, the CoFest demonstrated the process on a small number of selected workflow steps, while documenting how others could contribute. This created a foundation for onboarding new contributors and, over time, enabling more contributors to create wrappers for other widely used Bioconductor methods.

What we learned

One of the clearest lessons was that interoperability is not just about writing software that connects one platform to another.

It also depends on deciding where information should live, who is best placed to maintain it, and how contributions can fit naturally into the practices of each community. Bringing Bioconductor package developers, Galaxy contributors, and nf-core developers into the same room made it possible to look at the process end to end rather than from the perspective of a single ecosystem.

The event also reinforced the value of working from a concrete scientific use case. Using a real single-cell workflow gave participants something shared to build, test and discuss, while exposing questions that would have been difficult to anticipate from design discussions alone. The CoFest notes capture several of those emerging questions, particularly around how reusable jobs should be structured and composed.

What comes next

The CoFest was intended as a starting point rather than a one-off sprint.

The immediate next step is to continue refining the workflow demonstrated during the event: improving the tooling and documentation, testing it on additional Bioconductor methods, and making it easier for new contributors to take part.

We are also planning follow-up events in the first quarter of 2027. These will provide an opportunity to bring the communities together again, build on the examples developed in Milton Keynes and involve more contributors in extending the approach to other widely used methods.

The longer-term aim is simple: make it easier for software developed in one research community to become useful, sustainable building blocks in another.

Interested in getting involved? Join the Bioconductor Zulip and the #biofair2026-workflow-sprint channel and follow the BiocJobs topic for ongoing discussion.

To leave a comment for the author, please follow the link and comment on their blog: Bioconductor community blog.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Continue reading: A Collaboration Fest for research software interoperability

Connecting to Eureka2 using a university’s Microsoft Authenticator based login

https://pacha.dev/blog — Tue, 06 Oct 2026 23:00:00 +0000

[This article was first published on https://pacha.dev/blog, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Before the main content: * With all the funding cuts to PhD studies, I’d appreciate if you can donate here.

I am using Omarchy, a Linux distribution based on Arch for those that want efficiency, and it is similar to Manjaro minus the desktop effects that use considerable resources.

At Surrey and probably other universities in the UK, we require to use Microsoft Authenticator to log to the VPN.

Today I had to connect to Eureka2, a powerful computing cluster, using my personal laptop. The official guide Working off campus with GlobalProtect instructs to download Global Protect from the Surrey App Store, a service that has been retired.

I installed Open Connect which is a FOSS alternative:

sudo yay -S openconnect

I had to adjust my email settings here, and verify that the default sign-in method is ‘Microsoft Authenticator – notification’.

Then I was able to connect using:

sudo openconnect --protocol=gp --usergroup=gateway vpn.surrey.ac.uk

After a few tries, what works it to use the full username (not email name) with ‘@surrey.ac.uk’. What I mean is that my university account is something like ab12345 when my email firstname.lastname, and so I used ab12345@surrey.ac.uk. For the password I had to use the password that I use for Surrey Learn (not the email password) followed by the authenticator timed PIN number.

I was finally able to see:

ESP session established with server
ESP tunnel connected; exiting HTTPS mainloop.

And finally I could type this in another terminal:

ssh mv00622@eureka2.surrey.ac.uk

Which shows:

    +++++++++++++++++++++++++++++++++++++++++
      Welcome to Surrey's Eureka2 HPC Cluster
             Based on Rocky 8
    +++++++++++++++++++++++++++++++++++++++++

Documentation for HPC users is at below site:
https://docs.pages.surrey.ac.uk/research_computing/

I hope this helps. I had to spend over half an hour figuring out the steps by reading different sources and I hope this helps other people.

To leave a comment for the author, please follow the link and comment on their blog: https://pacha.dev/blog.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Continue reading: Connecting to Eureka2 using a university’s Microsoft Authenticator based login

{renv}, {rix}, rv, uvr, ir: An explosion of R package and environment managers

Isabella Velásquez — Tue, 06 Oct 2026 00:00:00 +0000

[This article was first published on R Works, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

There has been an explosion of package/environment/platform managers in the R ecosystem recently,¹ and R users are really reaping the rewards of so many ripe and riveting tools. But, I’ve found it rough² to keep track of them all and understand what they are and how they are different!

I suspect I am not the only one. My colleague Libby Heeren and I recently held a Data Science Lab with Eric Nantz to showcase {rix} and Nix (which was awesome! Keep an eye out on the Posit YouTube channel for when the recording is published), and the first three or so questions from the audience were all about how {rix} compares to X ({renv}, rv, uvr), and so I thought it’d be helpful to write a short summary of five tools that I’ve heard about: {renv}, {rix}, rv, uvr, and ir.³

While they’re all about tracking or isolating dependencies, they operate at a different “levels” from each other. TL;DR:

{rix} snapshots the entire system environment
{renv}, rv, and uvr manage R package dependencies within a single R project
ir handles dependencies for a single R script rather than a project

Graphic design is my passion.

I don’t want to/can’t pass judgement on which is the best or which I recommend, because I honestly haven’t had the opportunity to explore all of them fully (a very diplomatic way of saying that I haven’t yet used all of them in a real-world context), but I would love to hear everybody else’s opinions on which they think is the best tool because I imagine people have them!

It’s an unsatisfactory proxy for quality/experience, but here is the GitHub star trajectory of each:

Star History Chart

(Mind the recency bias).

With that, let’s begin, going chronologically by release!

{renv}

Repo | Docs | 1.2k GitHub stars | Packages R version Environments

Use case

This is a long-established package supported by Posit (creators of RStudio) used for project-local R package management.

How to use {renv}

{renv} is iterative and retrospective: you install the packages first, and then take a snapshot.

Install the {renv} package:

Terminal

install.packages("renv")

Within a project, initialize {renv}:

Console

renv::init()

Install packages:

Install packages the regular way, i.e., install.packages().

When you install packages, they will be installed into a project-specific library.

Snapshot your library’s packages into a lock file:

Console

renv::snapshot()

This lockfile captures the exact versions of those packages:

renv.lock

{
  "R": {
    "Version": "4.4.2",
    "Repositories": [
      {
        "Name": "CRAN",
        "URL": "https://cloud.r-project.org"
      }
    ]
  },
  "Packages": {
    "dplyr": {
      "Package": "dplyr",
      "Version": "1.1.4",
      "Source": "Repository",
      "Repository": "CRAN",
      "Hash": "fedd9d00c2944ff00a0e2696ccf048ec"
    },
    "ggplot2": {
      "Package": "ggplot2",
      "Version": "3.5.1",
      "Source": "Repository",
      "Repository": "CRAN",
      "Hash": "44c6a2f8202d5b7e6c72a95ca7ea1795"
    }
  }
}

Iterate:

As you continue adding packages, continue snapshotting to update the lock file.

Some considerations

Because {renv} is iterative and retrospective, you can accidentally install packages that conflict with the packages that you had installed earlier.

While {renv} tracks your R version, it intentionally does not manage it. If you need to switch R versions, you’ll want a tool like rig. See the renv documentation for the full list of caveats.

{rix}

Repo | Docs | 381 GitHub stars | Packages R version Environments

Use case

{rix} operates one level up from the other tools on this list. Instead of managing packages within an R installation, {rix} uses Nix, a general-purpose package manager, to snapshot an entire software environment: R itself, packages, and any system-level dependencies they need.

How to use {rix}

{rix} is retrospective but in a different sense than {renv}. You declare the environment you want up front, and Nix builds it from scratch (pulling from a package index rather than your currently-installed library).

Install Nix and the {rix} package:

Terminal

install.packages("rix")

Declare your environment:

Write a small script (e.g., generate_env.R) describing the R version and packages you need:

generate_env.R

library(rix)

rix(
  r_ver = "4.4.2",
  r_pkgs = c("dplyr", "ggplot2"),
  project_path = "."
)

Running this script generates a default.nix file, a plain-text description of the environment.

Build the environment:

Running nix-build creates the environment stored in /nix/store:

Terminal

nix-build

Then, you can drop into this environment interactively by running nix-shell.

From there, you can work normally by launching R, writing scripts, rendering Quarto docs, etc.

When you want to add a new package, edit generate_env.R and rerun that script to regenerate default.nix, and then run nix-build/nix-shellagain.

Some considerations

Because {rix} depends on Nix, there’s a bigger lift required to install and learn Nix itself. In exchange, you get reproducibility that extends beyond R packages to the R version and even system dependencies.

rv

Repo | Docs | 378 GitHub stars | Packages R version Environments

Use case

This is a relatively new tool by A2-ai, written in the ever-so-popular Rust language, used for project-local R package management.

It supports multiple package repositories with per-package source control (e.g., forcing a specific package to build from source vs. binary), and claims large speed improvements (20x+) over {renv} for installs.

rv is declarative: you write your desired packages and versions into a config file, and rv resolves the entire dependency graph upfront before installing anything.

How to use rv

Install rv:

rv is a standalone command-line tool, not an R package. There are a few installation instructions; see the documentation for more. I have Homebrew on my Mac, so I can run:

Terminal

brew install rv-r

Start a new project:

Start a new project with rv init. For example, running the below creates a new folder called hello_world:

Terminal

rv init hello_world

If you open up the hello_world folder, you’ll see that rv created various files and folders in your project folder:

hello_world/
  rv/
  library/
  scripts/
    activate.R
    rvr.R
  .gitignore
  .Rprofile
  rproject.toml

Update config file:

I don’t want to just repeat the documentation because it’s quite excellent, but essentially, you would fill in your packages as dependencies in rproject.toml:

rproject.toml

r_version = "4.4.2"

# A list of packages to install and any additional configuration
# Examples:
    # "dplyr",
    # {name = "dplyr", repository = "CRAN"},
    # {name = "dplyr", git = "https://github.com/tidyverse/dplyr.git", tag = "v1.1.4"},
dependencies = [
    "dplyr",
    "ggplot2"
]

Install packages:

Then, run rv sync to install the packages in the config file, as well as all of their dependencies:

Terminal

rv sync

This resolves the dependency graph and writes the exact versions to an rv.lock file, something like:

rv.lock

[[packages]]
name = "dplyr"
version = "1.1.4"
repository = "CRAN"
source = "Repository"
hash = "fedd9d00c2944ff00a0e2696ccf048ec"

[[packages]]
name = "ggplot2"
version = "3.5.1"
repository = "CRAN"
source = "Repository"
hash = "44c6a2f8202d5b7e6c72a95ca7ea1795"

Some considerations

Because rv is declarative and resolves the whole dependency graph upfront, it’s less forgiving of “install whatever you need as you go” exploratory work that {renv} accommodates; you’re expected to know (or figure out) your dependencies before you can sync.

You may notice that rv’s rproject.toml has an r_version field, but this does not actually manage the R version. Like {renv}, rv does not download, install, or switch R versions for you (again, you’d want a tool like rig for that).

If you want a single tool that tracks both your R packages and R version, and can actually switch R versions for you, you may want to check out…

uvr

Repo | Docs | 350 GitHub stars | Packages R version Environments

Use case

uvr is built by nbafrank, an independent developer. Like rv, it’s also written in Rust, and it is explicitly modeled on Python’s uv, a.k.a., the best thing to happen to the Python ecosystem in a decade.

It bundles R package management and R version management (two things that {renv} and rv treat separately), backed by a single lockfile (uvr.lock) that is meant to be the “source of truth” for a project.

How to use uvr

uvr, like rv, is declarative: you add packages (and optionally an R version constraint) to a manifest, and uvr resolves and locks the whole dependency graph before installing anything.

Install uvr:

uvr is a standalone binary, not an R package (though a companion R package, uvr-r, lets you drive it from the R console instead of the terminal).

Terminal

curl -fsSL https://raw.githubusercontent.com/nbafrank/uvr/main/install.sh | sh

Start a new project:

Terminal

mkdir my-project && cd my-project
uvr init --r-version ">=4.4.2"

This creates a uvr.toml manifest in the project folder.

Add packages:

Terminal

uvr add dplyr ggplot2

Each add updates uvr.toml and re-resolves uvr.lock. (uvr also supports Bioconductor and GitHub sources, e.g. uvr add DESeq2 --bioc or uvr add user/repo@main.)

Install everything from the lockfile:

Terminal

uvr sync

This is what generates and keeps uvr.lock up to date, capturing both the R version and package versions, something like:

uvr.lock

r-version = "4.4.2"

[[packages]]
name = "dplyr"
version = "1.1.4"
source = "cran"
hash = "fedd9d00c2944ff00a0e2696ccf048ec"

[[packages]]
name = "ggplot2"
version = "3.5.1"
source = "cran"
hash = "44c6a2f8202d5b7e6c72a95ca7ea1795"

Run scripts inside the isolated environment:

Terminal

uvr run analysis.R

ir

Repo | Docs | 65 GitHub stars | Packages R version Environments

Use case

ir operates one level down from the project-level tools. It’s built for single scripts or Quarto documents rather than whole projects. It’s a command-line tool for running portable R scripts and rendering Quarto documents whose package requirements (and optionally an R version) live inside the file itself for single-file workflows that don’t quite need a project, but still need to be easy to share and rerun later.

How to use ir

Install ir:

Follow the installation instructions in the repo:

Terminal

curl -fsSL https://raw.githubusercontent.com/r-lib/ir/main/scripts/install.sh | sh

Declare requirements in the file itself:

Add frontmatter (for .qmd/.Rmd) or a header comment (for .R scripts) listing the packages and, optionally, an R version or a Posit Package Manager snapshot date:

script.R

#!/usr/bin/env -S ir run
#| packages:
#|   - dplyr>=1.0
#|   - ggplot2@3.5.1
#| isolated: true
#| exclude-newer: "2025-05-15"

library(dplyr)
library(ggplot2)

mtcars |> count(cyl, gear) |> ggplot(aes(x = cyl, y = gear)) + geom_col()

report.qmd

---
title: My report
ir:
  packages:
    - dplyr>=1.0
    - ggplot2@3.5.1
  isolated: true
  exclude-newer: 2025-05-15
---

Run or render the file:

Terminal

ir run script.R
ir render report.qmd

When you run or render the file, ir resolves the requirements, prepares cached package libraries, and launches R or Quarto with a runtime ready to use.

Some considerations

ir is the newest tool here (first released mid-2026) and intentionally narrow in scope: it solves reproducibility for a single file, not a whole project. So, it’s a complement to, not a replacement for, {renv}/{rix}/rv/uvr when you actually do need a shared project-level environment.

Even more package and environment managers

This is just the surface of what is available out there. The amazing Novica Nakov has created a new R tooling webpage to track new R tools as they show up. So, if you want to follow an evolving list, I recommend keeping an eye out there!

While doing my research, I ran into a few more tools that I considered adding to this post and decided to not. Here are what they are, and why I left them out:

devenv: a general-purpose Nix-based dev environment manager, not specifically for R. It shares {rix}’s core mechanism but doesn’t have {rix}’s R-specific niceties like date-pinned CRAN/Bioconductor snapshots.
rpx, which I discovered on this Reddit thread as I was writing this: I didn’t add it because I’d waffled on finishing this post for too long and it would have messed up my Spelling Bee game. If anybody wants to write about it for R Works, let me know!
rig: a command-line tool for installing, listing, and switching between multiple R versions on your machine. I didn’t add it because there’s nothing to share with others (it’s more a personal R version manager), but as noted earlier, can be used to conform to the pinned R versions tracked in tools like {renv} or rv.

Footnotes

Their names create a fun Spelling Bee game:

︎
ruff? ︎
Since I’m listing both R packages and CLI tools, I am distinguishing them by wrapping the R packages in {}.︎

To leave a comment for the author, please follow the link and comment on their blog: R Works.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Continue reading: {renv}, {rix}, rv, uvr, ir: An explosion of R package and environment managers

Goodbye R-Forge, hello Codeberg

Achim Zeileis — Mon, 05 Oct 2026 22:00:00 +0000

[This article was first published on Achim Zeileis, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

After 20 years of service the R-Forge package hosting platform will be retiring soon. Codeberg is a modern alternative that shares many of the ideas and values that R-Forge was built on. A tutorial for transitioning to Codeberg (or other Git-based platforms) is provided.

You say goodbye

R-Forge was created 20 years ago after the 2nd useR! conference in 2006. Kurt Hornik had the idea that “R should have its own SourceForge!” (note that GitHub did not yet exist at the time). However, the goal was that this code sharing platform should be created by the R community for the R community. The implementation was led by Stefan Theussl in the following months using GForge/FusionForge, a fork of the software behind SourceForge. Hosting and funding was provided by WU Wirtschaftsuniversität Wien and The R Project.

R-Forge offers a Subversion (SVN) repository for collaborative development of R packages along with various other services that were state of the art at the time: daily package builds and checks, forums, mailing lists, bug trackers, web pages. A short introduction was published in the inaugural issue of The R Journal: Theußl & Zeileis (2009). “Collaborative Software Development Using R-Forge.” The R Journal, 1(1):9-14. doi:10.32614/RJ-2009-007

As the FusionForge system is no longer maintained, R-Forge will be retired in the next months, switching to a read-only site. The system was somewhat dated anyways: The website has a certain “00s” vibe (see the screenshot of my profile page below) and it lacks services that became standard over the last two decades like issues, pull requests, and continuous integration.

I say hello

Codeberg is a non-profit, community-led effort that provides services to free and open-source projects, such as Git hosting (using Forgejo), pages, continuous integration, etc. It is run by the non-profit organization Codeberg e.V. based in Berlin, Germany, which contributors can (but do not have to) become a member of. (I joined this summer.)

Thus, Codeberg shares many of the ideas and values, that R-Forge was built on, but provides state-of-the-art open-source tools for collaborative development of R packages (and beyond). For a first impression of the look and feel, see the screenshot of my profile page below.

I say go go go

For transitioning from R-Forge to Codeberg – or some other Git-based code hosting platform – Reto Stauffer has written up a step-by-step guide. The most important steps are also listed below.

1. SVN checkout from R-Forge

The starting point for the transition is the SVN repository from R-Forge that you probably have already checked out. For illustration we are using the c403 R package that is one of the packages contained in the exams project on R-Forge.

svn checkout svn://R-Forge.R-project.org/svnroot/exams/

2. List of R-Forge contributors

To find all SVN usernames who ever contributed to the SVN on R-Forge repository (in this example, to the folder exams/pkg/c403) the command svn log -q can be used. Using the following shell command a text file with the unique names can be created. Note that exams/pkg/c403 needs to be replaced for other projects and the output file name can be adjusted too:

(cd exams/pkg/c403 && svn log -q | grep -e '^r' | awk 'BEGIN { FS = "|" } ; { print $2 }' | sort | uniq) > \
  authors_c403.txt

This creates the text file authors_c403.txt which will look something like:

  azeileis
  retos

This file needs to be modified by adding the full names and e-mail addresses of the contributors. Important: The leading spaces need to be removed in each line for the subsequent commands to work.

azeileis = Achim Zeileis 
retos = Reto Stauffer

Provided the profile pages on R-Forge are up to date for the users, it is also possible to scrape this information, see the optional step 2*. below.

3. Clone a Git version of the SVN

Using git-svn a Git version of the SVN can be easily created. For this step it is crucial that all users are listed in the authors file created in the previous step, otherwise git-svn will fail.

git svn clone svn://R-Forge.R-project.org/svnroot/exams/pkg/c403 \
    --authors-file=authors_c403.txt --no-metadata \
    git-c403

Again, adjust the URL of the SVN, the authors file name, and the name of the git repository (git-c403) as needed.

4. Optional: Filter files in the Git repository

This step is not needed for the c403 example but is particularly useful when migrating a full R-Forge repository (rather than just a pkg sub-folder). In that case it might be useful to filter out certain files or directories such as the initial README file auto-created by R-Forge or the www sub-directory (typically just containing original PHP file or a auto-generated web page from pkgdown or altdoc).

Such files can be filtered, i.e., erased completely from the repository history, using git-filter-repo:

git-filter-repo --force --path 'README' --invert-paths
git-filter-repo --force --path 'www' --invert-paths

5. Import Git repository on Codeberg

First, create a new project on Codeberg, e.g., https://codeberg.org/zeileis/c403. The project should be completely empty without any initial commit.

Then, all left to do then is to tell the Git repository where its new home is (again adjust all names correspondingly):

cd git-c403
git remote add origin ssh://git@codeberg.org/zeileis/c403.git
git push --set-upstream origin main

Analogous commands can be used for setting the origin to a project on GitHub, GitLab, etc.

All history will be immediately available in the new repository. If the authors have accounts on Codeberg, their profiles will be linked based on the e-mail address(es) associated with their Codeberg account.

6. Package builds and checks via R-universe

To automatically build and check the development versions of the R package on Codeberg, the most convenient solution is to use the R-universe service developed by Jeroen Ooms and provided by rOpenSci.

In short, to make use of it, you need a GitHub account and create a repository like zeileis.r-universe.dev. This should contain a packages.json file with the names and repository URLs of one or more R packages, e.g.,

[
    {
        "package": "c403",
        "url": "https://codeberg.org/zeileis/c403"
    }
]

As a final step, the R-universe app has to be installed on the same GitHub account. See the R-universe set up documentation for more details.

7. Codeberg pages

If you want to add a project page for your package, e.g., generated via altdoc or pkgdown, then https://codeberg.page/ can be used.

The c403 package does not have a page, yet, but for an example see https://zeileis.codeberg.page/topmodels/ which is generated from the topmodels R package using altdoc in combination with quarto.

All that needs to be done to create the page manually, is to create an empty branch pages via git switch --orphan pages, adding the static site files, and pushing them to the repository.

Alternatively, it is also possible to use Forgejo actions to do this automatically.

8. Update descriptions and links

Finally, inform users where the sources of the package can be found, i.e., update the DESCRIPTION (especially the URL and BugReports fields) and in the NEWS. Additionally, a README.md file would be expected in a project while on R-Forge this was often not used (because R-Forge would be aware of the DESCRIPTION instead). Finally, it is recommended to add a link from the old R-Forge project to the new Codeberg project (possibly also in the www folder). To assure that none of the contributors can make modifications in the old R-Forge SVN, the permissions on R-Forge can be switched to read-only (typically for the Senior Developer role).

Optional addition: 2.* Scrape author names and e-mail addresses

In step 2, after the file with the author names (authors_c403.txt in the example above) has been created, one can try to automate the addition of proper names and e-mail addresses. This is particularly useful if there are more than just a handful of users. It works when the R-Forge user pages list the proper names and the (current) e-mail addresses.

The function below can be called in R as

rforge_authors("c403")

It first reads the authors_c403.txt file, looks up all users on R-Forge (using the rvest package), and then updates the authors_c403.txt file accordingly.

rforge_authors <- function(x) {
  txt <- sprintf("authors_%s.txt", x)
  if (file.exists(txt)) {
    a <- readLines(txt)
    a <- gsub(" ", "", a)
  } else {
    stop(sprintf("Please use 'svn log' to set up an '%s' file with all user names", txt))
  }

  for (i in seq_along(a)) {  
    doc <- a[i] |>
      sprintf(fmt = "https://r-forge.r-project.org/users/%s/") |>
      rvest::read_html()
    n <- doc |>
      rvest::html_table() |>
      lapply(as.data.frame) |>
      lapply(subset, subset = X1 == "Real Name:", select = "X2", drop = TRUE) |>
      unlist() |>
      unique()
    if (length(n) < 1L) n <- "FIXME"
    e <- doc |>
      rvest::html_table() |>
      lapply(as.data.frame) |>
      lapply(subset, subset = X1 == "Email Address:", select = "X2", drop = TRUE) |>
      unlist() |>
      unique() |>
      gsub(pattern = " @nospam@ ", replacement = "@", fixed = TRUE)
    if (length(e) < 1L) e <- "FIXME"
    a[i] <- sprintf("%s = %s <%s>", a[i], n, e)
  }
  writeLines(a, txt)
  invisible(a)
}

To leave a comment for the author, please follow the link and comment on their blog: Achim Zeileis.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Continue reading: Goodbye R-Forge, hello Codeberg

World crude oil stocks by @ellis2013nz

free range statistics - R — Mon, 05 Oct 2026 16:00:00 +0000

[This article was first published on free range statistics - R, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

How much crude oil does the world have available? Inquiring minds want to know but it is surprisingly difficult to get a single number. We hear more about relative levels “less than at any point since XXXX” and changes “has dropped X million barrels in the last month” than estimates of the absolute level. The International Energy Agency (IEA) publishes a lot of information on inventory levels and changes, but consolidated global stock estimates are, as far as I can see, not publicly available.

However, the Joint Organizations Data Initiative (JODI) does publish monthly country level estimates, reported to it voluntarily by countries. JODI partners are the APEC Energy Working Group, Eurostat, GECF, IEA, the Latin American and Caribbean Energy Organization (OLACDE), OPEC and the UN Statistical Division. It’s a very reputable set of partners! And JODI’s data are available to download.

Only 38 countries have data for every month JODI collects data, but this includes the USA, Japan, UK, the key European countries, and a few other big holders like Saudi Arabia, South Korea, Taiwan and Nigeria. The key omissions include, unfortunately but unsurprisingly, China, Russia, and Venezuela.. The data are also a bit old—latest data relating to July 2026 as at the time of writing, 6 October 2026, so a lag of 2-3 months. But it’s much better than nothing.

The point of today’s post is to get me the chart below, showing the best picture we can of the stocks of crude oil of the countries that do report those stocks, compared to a rough estimate of 85 million barrels per day of global crude refinery throughput:

That chart is now part of my fuel crisis monitoring page which gets updated at least weekly.

The overall story is that we indeed have less crude oil in storage than at any time since these records began in the early 2000s; but it’s still around 20 days of refinery cover. And at the rate that it’s going down, it will take nearly two years to get down to just 10 days of refinery cover. Actually, I’d expect market panic to set in at around 15 days of cover, but I’m being conservative here because I’m interested in the point where the world really has run out of any buffer against supply disruptions.

OK, the rest of the blog is mostly about just how I got that data, made some choices about excluding some countries, and drew the final chart.

First, downloading the data. There’s a separate CSV for each year, going back to 2002. Here’s how I download all those and import them into R

library(tidyverse)
library(glue)
library(janitor)
library(countrycode)

# Set to FALSE if you don't want to download the latest 2026 data. Will need to
# update all this workflow once we get to 2027 data.
update_2026 <- TRUE

#-----------------Downloads----------------

if (update_2026) {
  dir.create("fuel-crisis", showWarnings = FALSE)
  # Download current year. This will get more complete month by month.
  df <- here("fuel-crisis/jodi-2026.csv")
  url <- "https://www.jodidata.org/_resources/files/downloads/oil-data/annual-csv/primary/primaryyear2026.csv"
  download.file(url, destfile = df)
}


# Download historical data; only needed to be done once
if (!file.exists(here("fuel-crisis/jodi-2025.csv"))) {
  for (y in 2002:2025) {
    df <- here(glue("fuel-crisis/jodi-{y}.csv"))
    url <- glue(
      "https://www.jodidata.org/_resources/files/downloads/oil-data/annual-csv/primary/{y}.csv"
    )
    download.file(url, destfile = df)
  }
}

#----------------Import data------------------

# Import all the data to R:
jodi_hist_l <- list()
for (i in 1:25) {
  df <- here(glue("fuel-crisis/jodi-{i+2001}.csv"))
  jodi_hist_l[[i]] <- read.csv(df) |> as_tibble()
}

jodi_hist <- bind_rows(jodi_hist_l) |>
  clean_names()

# see https://www.jodidata.org/_resources/files/downloads/oil-data/jodi-oil-wdb-item-names-ver2017.pdf
# for guide on what everything is
# CLOSTLV is closing stocks
# CRUDEOIL is just crude
# TOTCRUDE also includes NGL and refinery feedstocks, additives and other hydrocarbons and has a lot more missing data.

Now that I’ve got the data, I want to understand the missingness. At one point I contemplated imputing the missing values of various countries, and hence I made a big wide version of the data with the idea that I might model crude stocks based on movement in some of the other numerous variables in the data. I eventually abandoned this when I realised how little meaningful information there was on eg China and Russia to do this meaningfully. But I can still use the jodi_wide data set created below for most of my plots:

#-----------explore which countries missing---------------------
jodi_wide <- jodi_hist |>
  #  filter(energy_product == "CRUDEOIL") |>
  mutate(
    dimensions = paste(energy_product, flow_breakdown, unit_measure, sep = "|")
  ) |>
  mutate(obs_value = as.numeric(obs_value)) |>
  select(ref_area, time_period, dimensions, obs_value) |>
  drop_na() |>
  spread(dimensions, obs_value, fill = NA) |>
  rename(crude = `CRUDEOIL|CLOSTLV|KBBL`) |>
  complete(ref_area, time_period, fill = list(crude = NA)) |>
  mutate(
    country = countrycode(
      ref_area,
      origin = "iso2c",
      destination = "country.name.en"
    )
  ) |>
  mutate(
    ref_area = fct_relevel(ref_area, "US"),
    date = ym(time_period),
    date_n = as.numeric(date)
  )

last_date <- jodi_wide |>
  filter(!is.na(crude)) |>
  summarise(ld = max(date)) |>
  mutate(ld = format(ld, "%B %Y")) |>
  pull(ld)

the_caption <- glue(
  "Source: JODI-OIL World Database. Most recent data is for {last_date}."
)


jodi_wide |>
  ggplot(aes(x = date, y = crude / 1000, colour = ref_area)) +
  geom_line() +
  theme(legend.position = "none") +
  scale_y_continuous(label = comma) +
  labs(
    x = "",
    y = "Millions of barrels",
    caption = the_caption,
    title = "Stocks of crude oil",
    subtitle = "Coloured by country; legend not shown."
  )

Here’s that first exploratory plot, a line chart of all the countries. The big one here is the USA; I’ve not shown the legend because it would take up too much space.

Next I want to identify which countries have missing data, and in particular those that are potentially big contributors to world crude oil stocks but are just missing a few observations that perhaps I could impute. I start by making a country_sum summary dataset:

possible_obs <- length(unique(jodi_wide$time_period))

country_sum <- jodi_wide |>
  group_by(ref_area, country) |>
  summarise(
    crude_total = mean(crude, na.rm = TRUE),
    missing_obs = sum(is.na(crude))
  ) |>
  ungroup() |>
  # particular problem if a country has any missing observations and averages
  # 10 million barrels or more of stock:
  mutate(
    all_there = as.logical(missing_obs == 0),
    some_missing = as.logical(
      missing_obs > 0 & missing_obs < possible_obs & !is.na(crude_total)
    ),
    problem = as.logical(crude_total > 10000 & missing_obs > 0)
  ) |>
  arrange(desc(crude_total))

good_countries <- country_sum |>
  filter(all_there) |>
  pull(country)

good_countries

That gives me this list of the countries that have a complete set of observations (largest crude oil stocks listed first):

> good_countries
 [1] "United States"  "Japan"          "Saudi Arabia"  
 [4] "Germany"        "South Korea"    "Canada"        
 [7] "France"         "Turkey"         "Italy"         
[10] "Poland"         "Spain"          "United Kingdom"
[13] "Netherlands"    "Taiwan"         "Thailand"      
[16] "Norway"         "Nigeria"        "Australia"     
[19] "Sweden"         "Finland"        "Czechia"       
[22] "Hungary"        "Portugal"       "Austria"       
[25] "Slovakia"       "Belgium"        "Denmark"       
[28] "Chile"          "New Zealand"    "Azerbaijan"    
[31] "Brunei"         "Ireland"        "Switzerland"   
[34] "Estonia"        "Iceland"        "Luxembourg"    
[37] "Latvia"         "Slovenia"

It also lets me find some countries that have no observations at all:

# 20 countries that never have any observations for Crude:
never_crude <- filter(country_sum, is.na(problem))
print(never_crude$country)

> print(never_crude$country)
 [1] "Albania"             "Bangladesh"         
 [3] "Bermuda"             "Belarus"            
 [5] "China"               "Egypt"              
 [7] "Georgia"             "Hong Kong SAR China"
 [9] "Moldova"             "Malta"              
[11] "Malaysia"            "Niger"              
[13] "Nepal"               "Sudan"              
[15] "Singapore"           "Syria"              
[17] "Eswatini"            "Tajikistan"         
[19] "Vietnam"             "Yemen"

Discovering that there were 20 of these, including China and Singapore, is what made me decide there was no point trying to impute the missing values. But I did a bit more exploration of which countries had some observations but not a complete set. Firstly, a scatter plot of all such countries:

…and then a faceted time series of the 9 biggest problematic countries.

I judged that maybe I could impute values for India, Mexico , Venezuela and Brazil based on their levels going up or down on average at about the world rate but decided this just wasn’t going to be worth it. And clearly Russia dwarfs them all and is just not worth trying, given the special circumstances there.

Here’s the code for those two plots:

country_sum |>
  ggplot(aes(
    x = crude_total,
    y = missing_obs,
    label = ref_area,
    colour = problem
  )) +
  geom_text() +
  scale_x_log10(label = comma) +
  labs(
    x = "Average crude stocks (thousands of barrels)",
    y = "Number of months missing an observation",
    caption = the_caption,
    title = "Size of crude stocks by number of missing observations"
  ) +
  theme(legend.position = "none")

# The nine biggest problem countries in terms of partly missing data
jodi_wide |>
  #  filter(ref_area %in% filter(country_sum, some_missing)$ref_area) |>
  filter(ref_area %in% filter(country_sum, problem)$ref_area) |>
  mutate(country = fct_reorder(country, crude)) |>
  ggplot(aes(x = date, y = crude / 1000)) +
  facet_wrap(~country, scales = "free_y") +
  geom_line(colour = "steelblue") +
  expand_limits(y = 0) +
  scale_y_continuous(label = comma) +
  labs(
    x = "",
    y = "Millions of barrels",
    title = "Crude oil stocks of countries missing at least one data point",
    caption = the_caption
  )

One last exploratory plot to try to see the impact of including or dropping those countries with partial data, here’s an area chart:

Thinking of the Russia situation in particular, I decided best just to exclude those partially missing countries altogether.

Here’s the code for that chart (and which also makes the crude_stocks data object I will shortly be using for my real plot)

#----------------Summarise and draw chart------------

# Crude oil stocks (excludes NGL etc because often missing data):
crude_stocks <- jodi_wide |>
  left_join(select(country_sum, ref_area, some_missing), by = "ref_area") |>
  group_by(date, some_missing) |>
  summarise(
    total_crude_mbbl = sum(crude, na.rm = TRUE) / 1000,
    reporting_countries = length(unique(ref_area))
  ) |>
  ungroup()

crude_stocks |>
  ggplot(aes(x = date, y = total_crude_mbbl, fill = some_missing)) +
  geom_area() +
  labs(
    caption = the_caption,
    y = "Total crude oil (millions of barrels)",
    x = "",
    title = "Total crude oil by missingness status of countries",
    fill = "Countries missing any data:"
  )

Finally, the code to draw the actual presentation plot. A key magic number here is the world refinery usage of crude oil at 85 million barrels per day; this is the approximate “crude runs” value from a recent IEA Oil Market Report. Note that this is less than the commonly cited values of around 105 million barrels per day, which includes NGLs, biofuels, and other liquid fuels. The 85m per day value is the appropriate one for me to use when looking at the CRUDEOIL series in my original JODI data; noting that they also have a TOTCRUDE series that would relate to the 105m per day figure, but which has a lot more missing data and is arguably not really what I’m most interested in anyway.

# note the commonly used figure of 105 is for TOTCRUDE, not just CRUDEOIl (which is more like 85)
world_use_per_day <- 85

# Summary for calculation of how long until only 10 days of cover:
war_summary <- crude_stocks |>
  # restrict it to countries that have observations for all months
  filter(!some_missing) |>
  summarise(
    prewar = total_crude_mbbl[date == "2026-03-01"],
    latest = total_crude_mbbl[date == max(date)],
    days = max(date) - as.Date("2026-03-01")
  ) |>
  mutate(
    rate = (prewar - latest) / as.numeric(days),
    ten_days_left = (latest - world_use_per_day * 10) / rate
  )

# draw chart:
crude_stocks |>
  filter(!some_missing) |>
  ggplot(aes(x = date, y = total_crude_mbbl)) +
  geom_line(colour = "blue") +
  geom_hline(yintercept = world_use_per_day * 10, colour = "darkred") + #expand_limits(y = world_use_per_day * 10) +
  expand_limits(y = 0) +
  scale_y_continuous(
    label = comma_format(suffix = "m"),
    sec.axis = sec_axis(
      ~ . / world_use_per_day,
      name = glue(
        "Days of crude inventory cover,\nat {world_use_per_day} million barrels per day"
      )
    )
  ) +
  labs(
    x = glue(
      "{length(good_countries)} countries in total have data for all months in this period."
    ),
    title = "Crude oil stocks of consistently reporting countries worldwide",
    y = "Millions of barrels",
    subtitle = glue(
      "Excluding countries with any missing data (eg China, Russia, India, Venezuela).
If the rate of decline since March 2026 continued, crude inventories would fall to 10 days of cover in {round(war_summary$ten_days_left / 7)} weeks."
    ),
    caption = the_caption
  )

That gets me this graphic, which is the one I will keep up to date from now on:

That’s all for today. Take care out there, and buy yourself an electric vehicle and rooftop solar if you can!

To leave a comment for the author, please follow the link and comment on their blog: free range statistics - R.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Continue reading: World crude oil stocks by @ellis2013nz

mlreserve: machine-learning loss reserving on ‘real-world’ triangles (based on a ‘ChainLadder’ fork)

T. Moudiki — Sat, 03 Oct 2026 00:00:00 +0000

[This article was first published on T. Moudiki's Webpage - R, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

This post presents a function mlReserve() for machine-learning loss reserving on ‘real-world’ triangles, based on a fork of R’s ‘ChainLadder’ package (https://github.com/thierrymoudiki/ChainLadder/tree/master). It includes two experiments, each scored against a known true reserve: (A) Two triangles from an individual-claims simulator (Wang & Wüthrich, https://github.com/actuarial-data-science/PackageIndividualClaimsSimulator), embedded below: chain ladder vs mlReserve, in detail; (B) A mini-benchmark on freshly simulated triangles from a simpler aggregate generator, with tunable distortions.

# =============================================================================
# mlReserve: machine-learning loss reserving on "real-world" triangles
#
# Colab: Runtime > Change runtime type > R
#
# Two experiments, each scored against a known true reserve:
#   A. Two triangles from an individual-claims simulator (Wang & Wüthrich),
#      embedded below: chain ladder vs mlReserve, in detail.
#   B. A mini-benchmark on freshly simulated triangles from a simpler aggregate
#      generator (defined in step 4), with tunable distortions.
# =============================================================================

# ---- 0. Setup -------------------------------------------------------------------

options(Ncpus = 2, repos = "https://cloud.r-project.org",
        repr.plot.width = 12, repr.plot.height = 5.5)
for (pkg in c("remotes", "ggplot2", "ranger", "e1071"))
  if (!requireNamespace(pkg, quietly = TRUE)) install.packages(pkg)
if (!requireNamespace("ChainLadder", quietly = TRUE) ||
    !"mlReserve" %in% getNamespaceExports("ChainLadder"))
  remotes::install_github("thierrymoudiki/ChainLadder", upgrade = "never")

suppressPackageStartupMessages({
  library(ChainLadder)
  library(ggplot2)
})

n      <- 10                                   # triangle size (accident years)
models <- c("Chain ladder", "mlReserve: SVM", "mlReserve: random forest")
ink    <- "#0b0b0b"
colours <- c("Chain ladder"             = "#eb6834",
             "mlReserve: SVM"           = "#2a78d6",
             "mlReserve: random forest" = "#1baf7a",
             "Truth"                    = ink)
theme_set(theme_minimal(base_size = 13) +
            theme(legend.position = "top", legend.title = element_blank(),
                  panel.grid.minor = element_blank(), plot.title.position = "plot",
                  strip.text = element_text(face = "bold")))


# ---- 1. Helpers -----------------------------------------------------------------

# A full n x n square of incremental payments -> observed (cumulative) upper
# triangle + true outstanding reserve per accident year (sum of the lower triangle)
make_case <- function(full) {
  dimnames(full) <- list(origin = 1:n, dev = 1:n)
  future <- row(full) + col(full) > n + 1
  upper  <- full; upper[future] <- NA
  list(full      = full,
       upper     = incr2cum(as.triangle(upper)),
       true_ibnr = rowSums(full * future))
}

# Fit the three models; `uncertainty = FALSE` gives point estimates only
fit_models <- function(triangle, uncertainty = TRUE, nsim = 300) {
  mse <- if (uncertainty) "bootstrap" else "none"
  fits <- list(
    MackChainLadder(triangle, est.sigma = "Mack"),
    mlReserve(triangle, e1071::svm, features = "numeric", transform = "asinh",
              mse.method = mse, nsim = nsim, seed = 1),
    mlReserve(triangle, ranger::ranger, num.trees = 300, num.threads = 1,
              features = "both", transform = "asinh",
              mse.method = mse, nsim = nsim, seed = 1))
  setNames(fits, models)
}

# Reserve (IBNR) and standard error by accident year, whatever the model
reserve_by_origin <- function(fit) {
  if (inherits(fit, "MackChainLadder")) {
    s <- summary(fit)$ByOrigin
    data.frame(origin = 1:n, IBNR = s$IBNR, SE = s$Mack.S.E)
  } else {
    s <- fit$summary[rownames(fit$summary) != "total", ]
    data.frame(origin = as.integer(rownames(s)), IBNR = s$IBNR, SE = s$S.E)
  }
}
total_reserve <- function(fit) sum(reserve_by_origin(fit)$IBNR)


# ---- 2. Data for experiment A ---------------------------------------------------

# Incremental paid (thousands), aggregated from the Wang & Wuthrich individual
# claims simulator (github.com/actuarial-data-science/PackageIndividualClaimsSimulator):
#   "Speed-up"     : settlement accelerates for recent accident years (delays -45%)
#   "All combined" : speed-up + calendar inflation shock (2% -> 15% a year from
#                    calendar month 78) + a few large, slowly paid claims
squares <- list(
  "Speed-up" = c(4172.1, 16153.4, 14800.2, 7561.1, 3275, 1856.3, 726.7, 410.7, 35.8, 140, 4983.2, 18064.1, 13574.5, 8350.7, 3954.9, 1499, 510.3, 9.5, 0, 11.3, 4698.1, 19428.8, 14864, 7343, 3710.5, 1378.6, 384.7, 92.2, 6.8, 27.1, 5467.1, 18481, 14847.6, 8089.7, 2794.1, 1379.2, 549.5, 152.2, 6.8, 0, 5985.2, 21516.6, 15980.6, 8353, 2954.2, 1655.5, 636.4, 191.3, 211.3, 28.9, 8274.1, 22014.4, 15981.4, 7425.8, 2362.7, 731, 242.1, 41.6, 52.8, 0, 8148.4, 24712.9, 16329.1, 6587.5, 2804.9, 753.6, 49.5, 0, 0, 0, 9974.2, 27573.4, 14523, 5141.6, 1123.8, 350, 73.3, 0, 0, 0, 10698.8, 29877.8, 14303.3, 3200.3, 978.3, 239.1, 0, 0, 0, 0, 13669.1, 30541.7, 10984, 1837.5, 157.3, 0, 0, 0, 0, 0),
  "All combined" = c(4236.8, 16766.9, 15671.5, 8229, 3582.8, 2143.9, 838.8, 548.6, 57.6, 261.6, 5195.7, 19008.8, 15926.9, 9857.6, 4643.8, 1746.5, 712.4, 26.4, 0, 22.5, 4956.8, 20861.6, 16283.4, 8704.2, 4321.7, 1876.2, 666.5, 161.8, 14.7, 61.4, 5908.2, 20336.4, 16580.1, 9331.1, 3673.8, 2276.2, 1205.8, 332.2, 15.4, 0, 6586.1, 24045.6, 18499, 11218.3, 4576, 3174.8, 1307, 460.6, 665, 96.3, 9292, 26084.3, 21055.7, 11396.7, 4203.8, 1531.8, 585.4, 115, 163.7, 0, 9702.9, 33010, 25202.7, 12133, 6404.4, 1956.8, 385.9, 0, 0, 0, 13752.5, 42686.8, 25857.7, 10698.4, 2706.5, 3213, 4232.5, 0, 0, 0, 17094.3, 53962.4, 30132.7, 8527.6, 3087.9, 746.2, 0, 0, 0, 0, 25375.8, 63473.7, 26195.3, 5737.2, 507, 674.2, 0, 0, 0, 0))
cases <- lapply(squares, function(v) make_case(matrix(v, n, n, byrow = TRUE)))
combo <- cases[["All combined"]]

cat("Observed cumulative paid, 'All combined' (thousands):\n")
print(round(combo$upper))
fits <- suppressWarnings(fit_models(combo$upper))
cat("\nmlReserve with a random forest; true total reserve =",
    format(round(sum(combo$true_ibnr)), big.mark = ","), "\n")
print(fits[["mlReserve: random forest"]])


# ---- 3. Experiment A: one triangle in detail ------------------------------------

# 3a. Reserve by accident year vs the truth
by_origin <- do.call(rbind, lapply(models, function(m)
  data.frame(model = m, reserve_by_origin(fits[[m]]))))
by_origin <- subset(by_origin, origin >= 2)             # AY 1 is fully developed
by_origin$model <- factor(by_origin$model, levels = models)
truth_points <- data.frame(origin = 2:n, IBNR = combo$true_ibnr[2:n])

print(
  ggplot(by_origin, aes(origin, IBNR / 1000, colour = model)) +
    geom_linerange(aes(ymin = (IBNR - SE) / 1000, ymax = (IBNR + SE) / 1000),
                   position = position_dodge(0.6), linewidth = 0.9) +
    geom_point(position = position_dodge(0.6), size = 2.6) +
    geom_point(data = truth_points, aes(origin, IBNR / 1000), inherit.aes = FALSE,
               shape = 23, size = 3.4, fill = ink, colour = "white") +
    scale_colour_manual(values = colours) +
    scale_x_continuous(breaks = 2:n) +
    labs(title = "Reserve by accident year, +/- one standard error (black diamonds: truth)",
         x = "Accident year", y = "Reserve (millions)"))

# 3b. Predictive distribution of the total reserve. Each model uses its own
#     resampling scheme: ODP bootstrap (BootChainLadder) for chain ladder,
#     residual bootstrap (mlReserve) for the learners.
set.seed(1)
boot <- BootChainLadder(combo$upper, R = 2000, process.distr = "od.pois")
draws <- rbind(
  data.frame(model = models[1], total = boot$IBNR.Totals),
  data.frame(model = models[2], total = rowSums(fits[[2]]$sims.reserve.pred)),
  data.frame(model = models[3], total = rowSums(fits[[3]]$sims.reserve.pred)))
draws$model <- factor(draws$model, levels = models)

print(
  ggplot(draws, aes(total / 1000, colour = model)) +
    geom_density(linewidth = 1, adjust = 1.2) +
    geom_vline(xintercept = sum(combo$true_ibnr) / 1000, colour = ink, linewidth = 0.8) +
    annotate("text", x = sum(combo$true_ibnr) / 1000, y = Inf, label = " true reserve",
             hjust = 0, vjust = 1.5, colour = ink) +
    scale_colour_manual(values = colours) +
    labs(title = "Model-specific predictive distributions of the total reserve",
         x = "Total reserve (millions)", y = "Density"))

cat("\nPredictive quantiles of the total reserve (thousands):\n")
quants <- t(sapply(split(draws$total, draws$model), quantile, c(0.025, 0.5, 0.975)))
print(noquote(formatC(round(quants), format = "d", big.mark = ",")))

# 3c. Projected cumulative payments vs the truth, accident years 7-10
cumulative_paths <- function(case, label) {
  cl <- MackChainLadder(case$upper, est.sigma = "Mack")
  rf <- mlReserve(case$upper, ranger::ranger, num.trees = 300, num.threads = 1,
                  features = "both", transform = "asinh", mse.method = "none", seed = 1)
  projections <- list("Truth"                    = t(apply(case$full, 1, cumsum)),
                      "Chain ladder"             = unclass(cl$FullTriangle),
                      "mlReserve: random forest" = unclass(rf$FullTriangle))
  out <- expand.grid(dev = 1:n, origin = 7:n, model = names(projections),
                     stringsAsFactors = FALSE)
  out$value <- mapply(function(o, d, m) projections[[m]][o, d],
                      out$origin, out$dev, out$model)
  out$panel <- paste0(label, " · AY ", out$origin)
  subset(out, dev >= n + 1 - origin)      # from the last observed diagonal onwards
}
paths <- do.call(rbind, Map(cumulative_paths, cases, names(cases)))
paths$model <- factor(paths$model, levels = c("Truth", models[c(1, 3)]))
paths$panel <- factor(paths$panel, levels = unique(paths$panel))

options(repr.plot.height = 7)
print(
  ggplot(paths, aes(dev, value / 1000, colour = model, linetype = model)) +
    geom_line(linewidth = 1) +
    facet_wrap(~ panel, nrow = 2, scales = "free_y") +
    scale_colour_manual(values = colours) +
    scale_linetype_manual(values = c("solid", "22", "solid")) +
    scale_x_continuous(breaks = seq(2, n, 2)) +
    labs(title = "Projected cumulative payments from the last observed diagonal",
         x = "Development year", y = "Cumulative paid (millions)"))
options(repr.plot.height = 5.5)


# ---- 4. Experiment B: mini-benchmark on fresh triangles -------------------------

# Aggregate generator: gamma payment pattern with over-dispersed Poisson noise.
#   speedup   : mean payment delay shrinks by up to this share for recent AYs
#   inflation : calendar-year log-growth rate after calendar period 7 (2% before)
simulate_case <- function(seed, speedup = 0.3, inflation = 0.08, phi = 5) {
  set.seed(seed)
  ultimate <- 5e4 * (1 + 0.03 * (0:(n - 1))) * exp(rnorm(n, 0, 0.05))
  delay    <- 3 * (1 - speedup * pmax(0, (1:n - 4) / (n - 4)))
  pattern  <- t(sapply(delay, function(m) diff(pgamma(0:n, shape = 2, scale = m / 2))))
  calendar <- outer(1:n, 1:n, "+") - 1
  infl     <- exp(0.02 * pmin(calendar, 7) + inflation * pmax(calendar - 7, 0))
  mean_incr <- ultimate * pattern * infl
  make_case(phi * matrix(rpois(n * n, mean_incr / phi), n))
}

n_rep <- 20                                   # triangles per scenario
scenarios <- list(
  # "Baseline (no distortion)" = c(speedup = 0,   inflation = 0.02),  # control
  "Settlement speed-up"      = c(speedup = 0.3, inflation = 0.02),
  "Speed-up + inflation"     = c(speedup = 0.3, inflation = 0.08))

bench <- do.call(rbind, lapply(names(scenarios), function(sc) {
  do.call(rbind, lapply(seq_len(n_rep), function(seed) {
    case <- simulate_case(seed, speedup   = scenarios[[sc]][["speedup"]],
                                inflation = scenarios[[sc]][["inflation"]])
    point_fits <- suppressWarnings(fit_models(case$upper, uncertainty = FALSE))
    data.frame(scenario = sc, seed = seed, model = models,
               error = sapply(point_fits, total_reserve) - sum(case$true_ibnr))
  }))
}))

score <- aggregate(error ~ scenario + model, data = bench, FUN = function(e)
  c(bias = mean(e), RMSE = sqrt(mean(e^2)), MAE = mean(abs(e))))
score <- do.call(data.frame, score)
names(score) <- c("scenario", "model", "bias", "RMSE", "MAE")
score <- score[order(score$scenario, score$RMSE), ]

cat("\nError of the total reserve over", n_rep, "triangles per scenario (thousands):\n")
fmt <- function(x) formatC(round(x), format = "d", big.mark = ",")
print(transform(score, bias = fmt(bias), RMSE = fmt(RMSE), MAE = fmt(MAE)),
      row.names = FALSE)

long <- rbind(data.frame(score[c("scenario", "model")], metric = "RMSE", value = score$RMSE),
              data.frame(score[c("scenario", "model")], metric = "MAE",  value = score$MAE))
long$model <- factor(long$model, levels = rev(models))
print(
  ggplot(long, aes(value / 1000, model, fill = model)) +
    geom_col(width = 0.6) +
    facet_grid(metric ~ scenario, scales = "free_x") +
    scale_fill_manual(values = colours, guide = "none") +
    labs(title = "Out-of-sample error of the total reserve (lower is better)",
         x = "Millions", y = NULL))

Observed cumulative paid, 'All combined' (thousands):
      dev
origin     1     2     3     4     5     6     7     8     9    10
    1   4237 21004 36675 44904 48487 50631 51470 52018 52076 52338
    2   5196 24204 40131 49989 54633 56379 57092 57118 57118    NA
    3   4957 25818 42102 50806 55128 57004 57670 57832    NA    NA
    4   5908 26245 42825 52156 55830 58106 59312    NA    NA    NA
    5   6586 30632 49131 60349 64925 68100    NA    NA    NA    NA
    6   9292 35376 56432 67829 72032    NA    NA    NA    NA    NA
    7   9703 42713 67916 80049    NA    NA    NA    NA    NA    NA
    8  13752 56439 82297    NA    NA    NA    NA    NA    NA    NA
    9  17094 71057    NA    NA    NA    NA    NA    NA    NA    NA
    10 25376    NA    NA    NA    NA    NA    NA    NA    NA    NA

mlReserve with a random forest; true total reserve = 174,050 
mlReserve (formula interface)

        Latest Dev.To.Date Ultimate   IBNR         S.E        CV
2      57118.1   0.9994925  57147.1     29    39.76128 1.3710786
3      57832.2   0.9982222  57935.2    103    91.80065 0.8912685
4      59311.6   0.9935674  59695.6    384   244.97640 0.6379594
5      68099.8   0.9732881  69968.8   1869   888.80702 0.4755522
6      72032.5   0.9286787  77564.5   5532  1605.55443 0.2902304
7      80048.6   0.8558364  93532.6  13484  2443.15573 0.1811892
8      82297.0   0.7087359 116118.0  33821  7158.48746 0.2116581
9      71056.7   0.5017785 141609.7  70553 13547.27988 0.1920156
10     25375.8   0.1881838 134845.8 109470 23912.17707 0.2184359
total 573172.3   0.7090055 808417.3 235245 37338.60673 0.1587222

Predictive quantiles of the total reserve (thousands):
                         2.5%    50%     97.5%  
Chain ladder             316,100 356,087 397,519
mlReserve: SVM           178,587 252,035 371,639
mlReserve: random forest 162,579 217,369 301,946

Error of the total reserve over 20 triangles per scenario (thousands):
             scenario                    model    bias    RMSE     MAE
  Settlement speed-up mlReserve: random forest  67,436  67,519  67,436
  Settlement speed-up             Chain ladder  94,022  94,097  94,022
  Settlement speed-up           mlReserve: SVM  93,888  94,318  93,888
 Speed-up + inflation mlReserve: random forest  58,834  58,975  58,834
 Speed-up + inflation           mlReserve: SVM  79,344  79,800  79,344
 Speed-up + inflation             Chain ladder 100,217 100,311 100,217

The experiment is not a universal proof that random forests are better. It is a demonstration that mlReserve() can be used to explore the performance of machine-learning models on loss reserving problems, and that it can outperform the classical chain ladder method in some scenarios.

To leave a comment for the author, please follow the link and comment on their blog: T. Moudiki's Webpage - R.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Continue reading: mlreserve: machine-learning loss reserving on ‘real-world’ triangles (based on a ‘ChainLadder’ fork)

pipeflow: v0.4.0 on CRAN

R some blog — Sat, 03 Oct 2026 00:00:00 +0000

[This article was first published on R some blog, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Pipeline setup

{pipeflow} builds a pipeline by adding one R function per step with pip_add (see also the Get started article). In this post, we’ll work with the following toy pipeline:

step load produces a sequence
step clean doubles it
step fit sums it, and
step report formats the result

pip <- pip_new("my-pip") |>
    pip_add(
        step = "load",
        fun = \(n = 5) seq_len(n),
        tags = c("io", "daily")
    ) |>
    pip_add(
        step = "clean",
        fun = \(x = ~load) x * 2,
        tags = c("io", "core")
    ) |>
    pip_add(
        step = "fit",
        fun = \(x = ~clean) sum(x),
        tags = c("model", "core")
    ) |>
    pip_add(
        step = "report",
        fun = \(x = ~fit) paste("result:", x),
        tags = "report"
    )

Each step must have a unique name, and function parameters can refer to other steps — for example, x = ~load takes the output of the load step. We also assign tags to each step, which we will use to filter steps later.

Before we run it, let’s have a look:

pip
 my-pip (4 steps)
---------------------------
     step params depends state       tags
1:   load      n           new   io,daily
2:  clean      x    load   new    io,core
3:    fit      x   clean   new model,core
4: report      x     fit   new     report
---------------------------
 last run: never

We see one row per step, showing its parameters (params), the steps it depends on (depends), its current state (all new before the first run), and our tags. Once we run the pipeline, …

pip_run(pip)
info [2026-10-03 09:36:44.793 UTC]: Starting run of pipeflow 'my-pip'
info [2026-10-03 09:36:44.794 UTC]: Step 1/4 load
info [2026-10-03 09:36:44.795 UTC]: Step 2/4 clean
info [2026-10-03 09:36:44.797 UTC]: Step 3/4 fit
info [2026-10-03 09:36:44.799 UTC]: Step 4/4 report
info [2026-10-03 09:36:44.803 UTC]: Finished run of pipeflow 'my-pip'
pip
 my-pip (4 steps)
---------------------------
     step params depends state            out       tags
1:   load      n          done      1,2,3,4,5   io,daily
2:  clean      x    load  done  2, 4, 6, 8,10    io,core
3:    fit      x   clean  done             30 model,core
4: report      x     fit  done     result: 30     report
---------------------------
 last run: 2026-10-03 11:36:44

… there is a new out column with the results, which can also be accessed directly:

pip[["report", "out"]]
[1] "result: 30"

Pipeline views in action

Real pipelines get long, and often you only care about a subset of steps: a topic, a stage, or the steps that produce the desired outputs. Views (new in v0.4.0) let you work on such a subset without copying anything. They reference the underlying pipeline, so every operation applied to it — running it, updating parameters, collecting output, … — writes through to the original pipeline, restricted to the steps the view covers.

Creating views

pip_view() returns a view containing only the steps that match the given filters:

pip_view(pip, tags = "core")
 my-pip view (2 of 4 steps)
------------------------------------------
    step params depends state            out       tags
1: clean      x    load  done  2, 4, 6, 8,10    io,core
2:   fit      x   clean  done             30 model,core
------------------------------------------
 last run: 2026-10-03 11:36:44
pip_view(pip, tags = "core", state = "done")
 my-pip view (2 of 4 steps)
------------------------------------------
    step params depends state            out       tags
1: clean      x    load  done  2, 4, 6, 8,10    io,core
2:   fit      x   clean  done             30 model,core
------------------------------------------
 last run: 2026-10-03 11:36:44
pip_view(pip, step = c("clean", "fit"))
 my-pip view (2 of 4 steps)
------------------------------------------
    step params depends state            out       tags
1: clean      x    load  done  2, 4, 6, 8,10    io,core
2:   fit      x   clean  done             30 model,core
------------------------------------------
 last run: 2026-10-03 11:36:44

By default steps must match all filters (logical AND), while the values within one filter are alternatives (OR). join = "union" keeps steps that match any filter, and fixed = FALSE treats filter values as regular expressions:

pip_view(pip, tags = "report", step = "clean", join = "union")
 my-pip view (2 of 4 steps)
------------------------------------------
     step params depends state            out    tags
1:  clean      x    load  done  2, 4, 6, 8,10 io,core
2: report      x     fit  done     result: 30  report
------------------------------------------
 last run: 2026-10-03 11:36:44
pip_view(pip, step = "^f", fixed = FALSE)
 my-pip view (1 of 4 steps)
------------------------------------------
   step params depends state out       tags
1:  fit      x   clean  done  30 model,core
------------------------------------------
 last run: 2026-10-03 11:36:44

Selecting steps with `[`

The extract operator (also new in v0.4.0) provides a data.table-like way of selecting steps and returns a view by default:

pip[c("load", "fit")]
 my-pip view (2 of 4 steps)
------------------------------------------
   step params depends state       out       tags
1: load      n          done 1,2,3,4,5   io,daily
2:  fit      x   clean  done        30 model,core
------------------------------------------
 last run: 2026-10-03 11:36:44

Boolean filters are evaluated in the context of the step table.

pip[tags %like% "core"]
 my-pip view (2 of 4 steps)
------------------------------------------
    step params depends state            out       tags
1: clean      x    load  done  2, 4, 6, 8,10    io,core
2:   fit      x   clean  done             30 model,core
------------------------------------------
 last run: 2026-10-03 11:36:44
pip[step %in% c("clean", "fit") & state == "done"]
 my-pip view (2 of 4 steps)
------------------------------------------
    step params depends state            out       tags
1: clean      x    load  done  2, 4, 6, 8,10    io,core
2:   fit      x   clean  done             30 model,core
------------------------------------------
 last run: 2026-10-03 11:36:44

Running views

Running a view executes the steps it covers together with any upstream dependencies that are not up to date. The run log marks the view’s own steps as [view] and the dependencies pulled in as [upstream]:

pip_reset(pip)  # reset the pipeline to "new" state for demonstration
pip_run(pip_view(pip, step = "report"))
info [2026-10-03 09:36:44.883 UTC]: Starting run of pipeflow 'my-pip view'
info [2026-10-03 09:36:44.883 UTC]: Step 1/4 [upstream] load
info [2026-10-03 09:36:44.883 UTC]: Step 2/4 [upstream] clean
info [2026-10-03 09:36:44.884 UTC]: Step 3/4 [upstream] fit
info [2026-10-03 09:36:44.885 UTC]: Step 4/4 [view] report
info [2026-10-03 09:36:44.886 UTC]: Finished run of pipeflow 'my-pip view'

Afterwards, the original pipeline is up to date for the covered steps. Let’s collect all results related to the “core” steps:

pip[tags %like% "core"] |> pip_collect()
$clean
[1]  2  4  6  8 10

$fit
[1] 30

For a complete guide on composing views, see the Pipeline views vignette.

More new and updated features

Views are only one of the additions. Each of the following has its own article on the documentation site:

Functional and method-style API — the pipeable pip_*() functions, with every function also available as a method (p$run(), p$view(), …): Get started.
Extract and replace operators — data.table-style [ and [[ for reading and editing steps: reference.
Collect and group output — pip_collect() with grouping by tags or dependencies: Collect and group output.
Combining and modifying pipelines — combine several pipelines with rbind(), or replace, rename, remove, and reset steps: Combine pipelines · Modify existing pipelines.
Split, map, and reduce — native exec = "split", "auto", and "reduce" modes: Split, map, and reduce.
Nested pipelines — use a pipeline as a reusable building block inside a step: Nested pipelines.
Self-modifying pipelines — .self, run-time restructuring, and restart(): Self-modifying pipelines.

Under the hood, significant performance gains have been made by implementing the dependency graph in C++ (since v0.3.0) and by improving the step and parameter bookkeeping in the underlying {data.table}. Finally, by dropping lgr and jsonlite, external package dependencies have been reduced to {data.table} and Rcpp.

{pipeflow} vs {targets}

{targets} remains R’s de-facto standard for heavy-duty, reproducible pipelines: it persists results to disk, records provenance, and scales to distributed infrastructure via crew. {pipeflow} does not try to compete on that home turf but targets a different niche — the interactive session, often as the backend of a Shiny app — where a pipeline is built and modified on the fly and has to respond while a user waits. Since responsiveness is key in that use case, {pipeflow} has been optimized for low latency.

The pipeflow vs targets article compares the in-session latency of the two packages head to head. The figure below gives a flavour of the results:

It shows the scenario where one parameter has changed¹ and the pipeline is rerun, with each step doing 5 ms of work. The x-axis represents the number of steps in the pipeline (its size), and the y-axis shows the time required to re-run it.²

Across the tested sizes, {pipeflow} is consistently the faster of the two — see the article for the full benchmark and methodology.

Overall, the two are best seen as complementary: reach for {targets} for large batch jobs that must be reproducible end to end, and for {pipeflow} when the pipeline is part of an interactive application.

Wrapping up

{pipeflow} 0.4.0 is on CRAN:

install.packages("pipeflow")

The full changelog lists every change, and the documentation includes the articles linked above. If you build Shiny apps or spend a lot of time exploring parameters interactively, give it a try — feedback and issues are welcome on GitHub.

Imagine a user changing the plot axis in a Shiny app.︎
That means checking which steps are affected by the parameter change and then executing those steps.︎

To leave a comment for the author, please follow the link and comment on their blog: R some blog.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Continue reading: pipeflow: v0.4.0 on CRAN

Major update to VS Code R extension

Seascapemodels — Thu, 01 Oct 2026 14:00:00 +0000

[This article was first published on Seascapemodels, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

I was doing some analysis work with R today (don’t get to do that very often anymore) and R wasn’t playing with VS Code nicely.

Turns out there has recently been a major update to the VS Code R editor extension.

So if you are having issues, check out the Extensions page for setup instructions.

Key changes:

New jgd package to improve the plotting experience. This replaces the httpgd package that was previously used to improve the plotting experience.
The r.term.path setting is now called r.consolePath (that’s the setting where you point VS Code to the R program)
Option to install afr to use as the console for R (instead of the base console or Radian).
It was easy to fix once I realized what was wrong. See the extension page for instructions.

I’m not yet sure if its better yet, afr does seem a bit smoother and faster than radian was.

I have also been trying out Positron lately, Posit’s VS Code clone. In general that seems much easier for R based work than VS Code always has been. It plays nicely with R out of the box. Unlike VS Code, which requires lots of additional settings and extension installs.

I switched from RStudio to VS Code as my primary IDE a few years ago. I also teach my AI data analysis workshops in VS Code. Though there is a fair bit of time spent in each workshop just helping people connect VS Code to R. At the time Posit just weren’t keeping up with the AI coding tools at that time. But they seem to be catching up again with Positron.

I’ve switched my teaching to recommend Positron. So if this extension update is giving you headaches, give Positron a go.

To leave a comment for the author, please follow the link and comment on their blog: Seascapemodels.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Continue reading: Major update to VS Code R extension

BioC2026 conference recap

Laurah Ondari — Wed, 30 Sep 2026 00:00:00 +0000

[This article was first published on Bioconductor community blog, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

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

Participants by country

Participants travelled to Seattle from across North America and beyond, reflecting the international reach of the Bioconductor community. Attendees joined from the United States, Canada, Germany, Ireland, the United Kingdom, Australia, and Nigeria. The audience included researchers at different stages of their careers, with staff representing 42% of attendees, followed by students and postdoctoral researchers at 28%, faculty at 13.8%, industry at 13.8%, and government at 2.4%.

Programme overview

Keynotes

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

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

Short talks

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

Poster sessions and package demo sessions

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

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

Participants during the poster session.

Celebrating 25 years of Bioconductor

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

Bioconductor at 25 celebrations.

Community and networking

One of the strongest themes to emerge from the BioC2026 survey was the value of the Bioconductor community itself. Participants described the conference as welcoming, collaborative, and interactive, and appreciated the opportunity to reconnect with colleagues and meet new people. All respondents who answered the question said they would recommend BioC2026 to their colleagues. The smaller size of the conference was also seen as a strength. Several respondents noted that the more intimate format made it easier to engage in conversations and build connections across different roles and career stages. Networking was consistently identified as one of the most valuable aspects of the event. Participants appreciated both the formal scientific programme and the informal conversations that happened around talks, posters, and demonstrations. Looking ahead, many attendees expressed interest in even more dedicated time for networking and informal interaction.

Infrastructure and tools

Zulip

During BioC2026, Zulip provided a dedicated space for participants to connect, interact and stay engaged throughout the conference. The conference space included topic-based discussions covering: – Workshops: Dedicated topics for individual workshops, including Single Cell, Bioconductor in Workflows, microbiome, WebR, and other training sessions. – Training & education: A space for broader conversations around training activities and educational sessions. – Social: Topics for informal conversations, social activities and connecting with other participants. – Introductions: Participants could introduce themselves and connect with others attending the conference. – Conference logistics channels: For travel, rooms, schedules, Day 1 information, and general conference support. – Community engagement: Dedicated spaces for the developer survey, student and early-career researcher activities, and other community initiatives. – Help and support: Dedicated topics such as the helpdesk, lost and found, and issues accessing virtual sessions.

Sticker contest winner

Kayva Vaddadi was BioC2026’s sticker design contest winner. Kayva’s design celebrated Seattle as the host city of BioC2026, featuring well-known landmarks such as the Space Needle, the Seattle Great Wheel, Mount Rainier, and the surrounding evergreen forests. Puget Sound and a Washington State ferry reflect Seattle’s maritime landscape and the culture of connection across the region. The BioC2026 hot air balloon on the sticker was inspired by the scenic balloon views popular in the Seattle area and represents the celebratory and exploratory spirit of the Bioconductor community.

See the announcement here: https://lnkd.in/p/emAWMGX

Conference materials

All conference recordings are available on the Bioconductor YouTube channel. BioC2026 workshops will also be available in the “Archived” section of workshop.bioconductor.org, alongside workshops from previous BioC, EuroBioC, and BiocAsia events. If you’d like to add additional materials or use the platform for your own workshop or course, contact us at support.bioconductor.org.

Post-conference events

Orchestrating Large-Scale Single-Cell Analysis with Bioconductor

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

Participants during the scRNA-seq workshop.

Bioconductor North America 2026 Hackathon

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

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

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

Participants during the Hackathon.

Up next…

Later this year, the Bioconductor community will head to Melbourne, Australia, for BioCAsia2026, taking place on November 19-20, 2026. BioCAsia2026 brings together researchers, developers, and data scientists from across the Asia-Pacific region and beyond to share advances in bioinformatics, computational biology, and open-source software. The meeting will highlight practical training, community building, and the latest developments in R and Bioconductor. BioC2027 will then take place in early summer in Boston, Massachusetts, continuing the North American Bioconductor conference series and bringing the community together once again for scientific exchange, collaboration, and community building. The community will thereafter gather in Basel, Switzerland, for EuroBioC2027, taking place from September 8–10, 2027. The European Bioconductor Conference will continue to bring together researchers, developers, and educators from across Europe and beyond to share new software, methods, and applications in computational biology.

Acknowledgements

BioC2026 gratefully acknowledges the support of all sponsors and partners whose contributions and support made the conference possible.

Organising committee

We thank the local organisers, programme committee, workshop instructors, keynote speakers, volunteers, sponsors, and all participants whose contributions made BioC2026 a success.

Marc Carlson, Seattle Children’s Hospital
Nicholas Cooley, University of Limerick
Sean Davis, University of Colorado Cancer Center
Maria Doyle, University of Limerick
Mikhail Dozmorov, Virginia Commonwealth University
Jenny Drnevich, University of Illinois at Urbana-Champaign
Russell Eberts, Fred Hutchinson Cancer Center
Erica Feick, Dana-Farber Cancer Institute
Sana Hirata, Fred Hutchinson Cancer Center
Ted Laderas, Fred Hutchinson Cancer Center
Alexandru Mahmoud, University of Limerick
Matthew McCall, University of Rochester Medical Center
Lori (Shepherd) Kern, Roswell Park Comprehensive Cancer Center
Laurah Nyasita Ondari, International Institute of Tropical Agriculture (IITA)
Charlotte Soneson, Friedrich Miescher Institute for Biomedical Research
Tim Triche, Van Andel Institute
Levi Waldron, CUNY Graduate School of Public Health and Health Policy
Wes Wilson, Princess Margaret Cancer Centre

To leave a comment for the author, please follow the link and comment on their blog: Bioconductor community blog.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Continue reading: BioC2026 conference recap

My first attempt at porting a Shiny app failed, and not for a technical reason

Arthur Bréant — Tue, 29 Sep 2026 16:14:32 +0000

[This article was first published on Rtask, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

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

Diary of a rewrite, episode 2 of 4.

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

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

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

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

Here is how I went about it.

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

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

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

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

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

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

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

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

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

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

The problem is all the others.

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

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

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

Why inventorying the observers is the wrong entry point

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

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

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

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

It has been my working rule ever since.

Start from the data being served

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

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

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

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

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

60 lines with nothing left to do

Here is the example that tipped my understanding.

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

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

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

On the htmxr side, one argument:

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

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

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

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

What happens to each building block

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

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

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

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

updateSelectInput() becomes: nothing. We no longer update a control. The server renders the URL ︎ inputs sync low the push_url argument static ggplot low svglite + htmx swap DT (display) low hx_table(), hx_table_rows() long task + notifications medium external job + periodic polling interactive plotly high plotly.js, or drop the interactivity rhandsontable high JS library, or one form per cell

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

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

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

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

The most useful side effect

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

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

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

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

What you do not get back

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

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

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

Conclusion

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

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

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

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

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

Diary of a rewrite, four episodes:

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

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

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

To leave a comment for the author, please follow the link and comment on their blog: Rtask.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Continue reading: My first attempt at porting a Shiny app failed, and not for a technical reason

August 2026 Top 40 New CRAN Packages

Joseph Rickert — Tue, 29 Sep 2026 00:00:00 +0000

[This article was first published on R Works, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

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

Agriculture

agridatasets v0.1.1: Offers a rich and diverse collection of datasets focused on agriculture, agronomy, animal science, and related fields. The package includes experimental, observational, and field-trial data on crops such as rice, wheat, corn, soybean, cotton, coffee, avocado, and orange, as well as forestry species including bamboo, eucalyptus, and timber. Datasets cover plant breeding and genetics, factorial and randomized block experiments, herbicide and insecticide efficacy trials, pest and disease infestation, soil characteristics and land suitability, plant growth regulators, seed germination, and crop yield modeling. See the vignette.

Artificial Intelligence

commons v0.1.0: Implements trustworthy large language model agents. Connect raw data sources, provides a pool of trusted calculations, and a searchable context layer that demonstrates how to interpret them. Then, deploy data agents that answer questions, log interactions, and can be evaluated and improved over time. See the vignettes Introduction and Security and governance.

diffuseR v0.2.2: A native R implementation of diffusion models providing a functional interface to state-of-the-art generative AI. Inspired by the Python library diffusers from Hugging Face, functions generate and manipulate images from text prompts using models such as Stable Diffusion, with no Python dependency. Supports multiple diffusion schedulers and device acceleration. See the vignette and README for more information.

Rhobots v0.1.10: Implements the BERTopic topic modeling pipeline directly in R: Provides transformer-based sentence embedding, uniform manifold approximation and projection dimensionality reduction, hierarchical density-based spatial clustering of applications with noise clustering, and class-based term frequency-inverse document frequency topic extraction without any dependency on Python, conda, or reticulate. Every stage runs in R through torch, safetensors, tok, uwot, and dbscan. The package mirrors the accessor API of the original Python package, adds integrated quality metrics and hyperparameter search tools, and introduces part-of-speech filtered and C-value-ranked representation models. See README for more information.

Climate Studies

topocast v0.0.5: Downscales coarse-resolution raster data to a finer grid by fitting local linear regressions of a response, such as a climate variable, on one or more fine-resolution predictors, such as elevation and other terrain indices, within a moving window. Multiplicative and additive anomaly application downscale time series relative to a baseline climatology. Follows the regression-on-elevation approach used for high-resolution climate surfaces described in Karger et al. (2017). See the vignettes Getting started and How moving-window downscaling works.

Computational Methods

openfhe.R v1.5.1.1: Provides an R interface to penFHE, the open-source C++ library for fully homomorphic encryption Badawi et al. (2022), which allows computation directly on encrypted data without access to the secret key. Supports the Brakerski-Fan-Vercauteren (2012), Brakerski-Gentry-Vaikuntanathan (2014), and Cheon-Kim-Kim-Song (2017) schemes for arithmetic on encrypted numbers. There are three vignettes including Introduction and CKKS Bootstrapping.

Ecology

FINN v0.1.0: Implements a hybrid dynamic forest model that can be configured as a fully mechanistic, process-based model, like classic forest gap models, or with its demographic processes (growth, mortality, regeneration) replaced by deep neural networks, or any combination of the two. Provides functions to define a model and its mechanistic or empirical components, calibrate it to forest inventory data, and interpret the calibrated processes. FINN is implemented with the torch package but no knowledge of torch is required. The hybrid modeling approach is described in Pichler and Käber (2026). There are five vignettes including Introduction and Mortality: a binomial response and a neural-network process.

gbif.range v1.9.2: Implements an end-to-end workflow to generate ecologically informed species range maps from sparse observations using environmental clustering and convex hulls. Serves as a standalone framework or complementary approach to species distribution models. By constraining estimated ranges within authoritative or custom ecoregion boundaries, the approach prevents spurious range over-prediction common in geometric hull methods. There are five vignettes including Getting Started and Part 2: Ecoregion-Based Range Inference.

PhysMove v1.2.5: Provides tools to analyse animal movement and space-use patterns from telemetry data using methods derived from statistical physics. Methods span displacement-based approaches, distribution fitting, space-use metrics, including the influence of correlations on space-use, network-based community detection, and measures of entropy and predictability. The package enables characterization of these patterns across spatial and temporal scales. For applications of these methods in ecological studies see Rodríguez et al. (2017) and Sequeira et al. (2018). There are four vignettes including Introduction and Movement Patterns.

Epidemiology

RtForecastR v0.1.1: Provides functions for filtered (real-time/causal) and smoothed (retrospective) estimation of the time-varying effective reproduction number from case-count time series, using the EpiFilter algorithm of Parag (2021), together with a one-step-ahead in-sample prediction check, a genuine out-of-sample one-step forecast with predictive intervals, elimination probability , and forecast calibration metrics. See the vignette.

Genomics

diffwrap v0.6-3: Provides functions for differential expression analysis of read counts from messenger RNA (mRNA) sequencing (RNA-Seq) data or micro RNA (miRNA) expression values generated by the Comprehensive Analysis Pipeline for microRNA Sequencing expression_reports.sh script. The workflow follows the Bioconductor edgeR–limma expression data analysis pipeline providing options for different approaches, such as pure edgeR, voom or paired samples. Methods are described in Robinson, McCarthy and Smyth (2010), Ritchie et al. (2015), Law et al. (2014 and Sun et al. (2014). See the vignette.

Machine Learning

figsr v0.1.1: Implements a flexible, interpretable machine learning algorithm for additive tree sums. Fits a sum of shallow classification and regression trees (CART) by greedily minimizing residual impurity, growing a new tree or deepening an existing one at each step, whichever reduces the residuals most. Supports regression and two-class classification, variable importance, bootstrap ensembling and seamless integration with parsnip and tidymodels workflows. The method is described in Tan et al. (2023). See the vignette.

neuralsbi v0.3.2: Provides a native R implementation of a neural simulation-based estimator that runs on the torch back end and is focused on Neural Posterior Estimation. Given a prior over parameters and a simulator, functions train a conditional neural density estimator to approximate the Bayesian posterior, enabling amortized, likelihood-free inference. It targets applied researchers who want an approachable interface with sensible defaults and built-in posterior diagnostics. There are four vignettes including Getting started and Case study: inferring epidemic parameters (SIR).

NBvarsel v0.1.1: Performs exhaustive or groupwise (backward elimination) variable selection for binary outcome prediction models using cross-validated net benefit as the optimization criterion. It supports predictor costs, restricted cubic splines, interaction terms, permutation importance, and parallel computation. It includes visualizations for model comparison and variable importance. References include Vickers & Elkin (2006) and Van Calster et al. (2018). See the vignette.

Medical Statistics

CausalState v0.10.2: Implements sequential doubly robust and infinite-dimensional targeted maximum likelihood estimators for longitudinal modified treatment policies in settings with transitioning states, such as ICU, ward, or emergency department care episodes. Treatment is permitted in active states and becomes structurally inapplicable after a state transition (e.g. discharge or death). Methods based on Diaz et al. (2021) and Luedtke et al. (2017). There are three vignettes including Getting started and Wu-Benkeser density-ratio metalearner.

orthoMTL v0.1.0: Fits regularized multi-task learning models where relationships between tasks are controlled via orthogonality or disjoint-support constraints. Supports regression, binary classification, and censored survival data. In survival mode, time-to-event outcomes are converted into binary labels at user-defined thresholds, enabling the discovery of features with time-varying effects that standard proportional-hazards models cannot detect. Implements the penalty described in Vervier et al. (2014). See the vignette.

expoquimR v0.1.0: Provides a unified toolkit for occupational chemical exposure risk assessment, implementing three internationally recognized methods: the qualitative control-banding methods COSHH Essentials (UK Health and Safety Executive) and the method of the French National Research and Safety Institute (INRS), together with the quantitative statistical procedure of the UNE-EN 689 standard for comparing measured exposure levels against occupational exposure limits. Every step of each method is implemented as a small, independently callable, and unit-tested function, so assessments are reproducible and auditable. Optional shiny applications provide a guided, interactive workflow. References: UK Health and Safety Executive (2003) and Mallet, Pilorget and Berne (2013), ISBN:978-2-7389-2166-2. There are three vignettes including COSHH Essentials: Qualitative Chemical Risk Assessment and INRS Method: Qualitative Inhalation Risk Assessment

rdborrow v0.0.4.0: Implements causal inference methods for incorporating external control data into randomized controlled trials with longitudinal outcomes. Provides an analysis module supporting weighting-based methods such as inverse probability weighting and augmented inverse probability weighting, difference-in-differences. Methods are based on Zhou et al. (2024) and Zhou et al. (2024). There are five vignettes including Introduction and Primary Analysis Workflow.

Meta Analysis

metaselection v0.3.0: Fits a flexible class of p-value selection models for meta-analysis and meta-regression models, providing standard errors and confidence intervals based on either cluster-robust variance estimators (i.e., sandwich estimators) or cluster-level bootstrapping to handle dependent effect size estimates, as described in Pustejovsky, Citkowicz, and Joshi (2025) and Citkowicz, Pustejovsky, and Joshi (2026). Supported models include generalizations of the step-function selection model as proposed by Vevea and Hedges (1995) and the beta-function selection model as proposed by Citkowicz and Vevea (2017). See the vignette.

Networks

idiographic v0.3.4: Provides functions to make person-specific and within-person network estimation from intensive longitudinal and panel data. Estimators include ordinary vector autoregression (VAR), graphical vector autoregression, multilevel vector autoregression, rolling ordinary and graphical VAR, native Bayesian VAR, multilevel Bayesian VAR, unified Structural Equation Modeling, and Group Iterative Multiple Model Estimation. Methods are described in Saqr et al. 2025 and Epskamp et al. (2018). There are seven vignettes including an introduction and Multilevel VAR.

lame v1.3.4: Implements additive and multiplicative effects models for both cross-sectional and longitudinal network analysis and supports square and rectangular network structures. Key features include: (1) Cross-sectional network analysis with support for binary, continuous, ordinal, and count data; (2) Longitudinal network analysis with additive sender/receiver and multiplicative latent-factor effects that can evolve over time through AR(1) processes’ (Sewell and Chen (2015) and Durante and Dunson (2014)); (3) Handling of changing actor compositions across time periods in longitudinal models; and (4) Performance improvements. There are seven vignettes including Overview and Getting Started.

tabulergm v0.1.0: Creates publication-ready tables documenting exponential-family random graph models (ERGMs), a class of statistical models for social networks (Robins et al. (2007)). Tables describe model terms through their definitions, mathematical representations, and graphical representations, and can be generated from ERGM formulas or from models fitted with the ergm package (Hunter et al. (2008). Resulting tables can be integrated into quarto and rmarkdown documents.

Process Control

gpci v0.1.0: Implements a comprehensive, generalized framework for computing, estimating, and validating generalized process capability indices. Supports user-supplied probability density functions, cumulative distribution functions, survival functions, and quantile functions with uncensored data parameter estimation via Maximum Likelihood Estimation. Provides several classical and non-normal capability indices, including Cpy (Maiti, Saha and Nanda, (2010), Spmk (Dey and Saha (2019), CpTk (Saha, Dey and Maiti (2019) and others. Functions also compute parametric and non-parametric bootstrap confidence intervals, confidence levels using percentile, highest posterior density intervals and Heidelberger-Welch convergence diagnostics. See Maiti, Saha and Nanda (2010) and Saha, Dey and Maiti (2018) for background. See the vignettes Getting Started and Using Custom Distributions.

Programming

polyglotSQL v0.1.0: Provides functions to parse, tokenize, validate, format, analyze and translate SQL between more than 30 dialects (PostgreSQL, MySQL, BigQuery, Snowflake, DuckDB, T-SQL, and others) using polyglot-sql Rust crate. All processing happens locally in the R session. Includes column-level lineage, structural query analysis, query optimization, AST diffing and OpenLineage facet generation. There are four vignettes including Getting Started and Parsing, validation and lineage.

Public Transit

transittraj v1.1.0: Today’s public transit vehicles produce a large amount of automatic vehicle location (AVL) data which is very useful for planning and performance studies, but can be noisy, error-prone, and sparse. This package provides tools for cleaning AVL point data and turning it into continuous, differentiable, monotonic, and invertible vehicle trajectory functions, based on the work of Robbennolt et al. (2025) and Huang et al. (2023). See the vignettes Introduction and The AVL Cleaning Workflow.

Statistics

chaidr v0.1.0: Implements the CHAID (Chi-squared Automatic Interaction Detection) decision tree algorithm of Kass (1980) and the Exhaustive CHAID variant of Biggs, de Ville, and Suen (1991), as specified in the IBM SPSS Statistics Algorithms documentation. Supports nominal, ordinal with floating missing category, and continuous predictors, and nominal, ordinal, and continuous response variables using Pearson chi-squared, Goodman row-effects, and one-way ANOVA F tests respectively. Includes prediction, rule extraction, gains and lift analysis, validation on holdout data, and visualization via base graphics, Graphviz DOT, plotly, and conversion to partykit objects. See the Introduction and the Japanese tutorial.

citcdf v1.1.0: Enables complex hypothesis testing through conditional cumulative distribution function estimation. Method is detailed in: Gauthier et al. (2021). See the vignette.

falsifyr v1.0.0: Attacks fitted R model claims by searching for small plausible perturbations that make a target result disappear. The package focuses on claim-level fragility, smallest-kill reporting, and reproducible caveated robustness checks for ordinary fitted model objects. The methods draw on the fragility-index concept of Walsh et al. (2014), multiverse analysis of Steegen et al. (2016), specification-curve analysis of Simonsohn et al. (2020), and robust covariance estimation of Zeileis (2004). See the vignettes Attacking a Regression Claim and Interpreting Survival Scores.

FPScausal v0.1.1: Implements functional propensity score weighting for causal inference with functional treatments. The method estimates weights that balance observed confounders by removing their dependence on the functional treatment and uses a dual formulation of the weighting problem for efficient unconstrained optimization. The framework supports scalar, binary, and functional outcomes, as well as functional covariates, and can be used to estimate marginal causal effects in settings with time-varying exposures. The methodology follows Ciardulli et al. (2026). See the vignette.

LRErdd v0.1.0: Provides functions for the design and analysis of Regression Discontinuity Designs as local randomized experiments within the potential outcome approach as formalized in Li, Mattei and Mealli (2015) including functions to implement the design phase of the study, where the focus is on the selection of suitable subpopulations for which valid causal inference can be drawn. These functions provide summary statistics of pre- and post-treatment variables by treatment status and select suitable subpopulations around the threshold where pre-treatment variables are well balanced between treatment groups. See the vignette. There is also a Shiny application.

mvboxcox v0.1.4: Fits bivariate logistic Box-Cox regression models for binary outcomes and positive continuous predictors. Transformation parameters are selected by cross-validated grid search with adaptive refinement and thin-plate spline smoothing. The package also provides prediction, empirical and sampling-weighted median effects, simulation tools, and sampling-weighted model fitting. The methodology extends the logistic Box-Cox approach of Xing et al. (2021). See the vignette.

Surveys

sondage v0.9.1: Implements survey sampling algorithms for single-stage probability sampling from finite populations, written in C. Provides equal probability methods (simple random sampling, systematic, Bernoulli), unequal probability methods (conditional Poisson / maximum entropy, Sampford, Brewer, systematic PPS, Pareto, sequential Poisson, Poisson, Chromy’s minimum replacement, multinomial), balanced sampling via the cube method, and spatially balanced sampling via the local pivotal method and spatially correlated Poisson sampling. Functions compute joint inclusion probabilities, pairwise expectations, and sampling covariances for variance estimation. See Tillé (2006) for background. There are two vignettes Getting Started and Extending sondage with Custom Methods.

Time Series

scanr v0.1.1: Detects change points in long univariate time series using the SCAN framework. The implementation uses a native Rust backend exposed to R via extendr. See the vignette.

Utilities

rmoriebricklayer v0.5.1: Tools for building brick-proof, reproducible, self-contained data capsules. Resolves open-data sources through the Comprehensive Knowledge Archive Network package_showand package_search endpoints, records and verifies provenance with Secure Hash Algorithm 256 (SHA-256) digests and Internet Archive Wayback Machine snapshots, validates downloaded data against a pinned schema. Run records are captured in a manifest plus a plain-language summary so any result can be traced back to its inputs. Tests distributional drift between a pinned capsule and a fresh fetch because a re-released extract can be statistically identical yet differ byte-for-byte. Manifests can be authenticated with keyed digests (HMAC-SHA-256, RFC 2104) or post-quantum hash-based signatures. There are eight vignettes including Building reproducible data capsules and Provenance you can verify.

scimesh v0.4.0: Implements a fast, GPU-free 3D software renderer written in modern C++17 with native R bindings. Renders triangle meshes to publication-quality images entirely on the CPU, requiring no display server or graphics hardware. Features multi-light Blinn-Phong shading, screen-space ambient occlusion, anti-aliasing, depth fog, transparency, wireframe rendering, texture mapping, screen-space lines and text labels, and procedural geometry generation. Supports standard mesh file formats with PNG and PPM output. Works on high-performance computing clusters, headless servers, containers, and continuous integration pipelines, making it suitable for scientific visualization across neuro-imaging, molecular structures, and general 3D graphics. See the vignette.

Visualization

dgraphs v0.2.0: Constructs data-derived graphs from numerical observations using mutual, shared-neighbor, intersection, geodesic, radius, adaptive-radius, and minimum-spanning-tree completion methods. Provides graph conversion, pruning, diagnostics, spectral embedding, endpoint detection, and path utilities. The implemented graph constructions include methods described by Jarvis and Patrick (1973), Brito et al. (1997), Berry and Sauer (2019), and Gower and Ross (1969). See the vignette.

gghotelling v0.2.1: Calculate Hotelling’s T² ellipses and detect multivariate outliers both for base R plots and ggplot2 plots. Optionally, uses robust covariance estimation to reduce the influence of outliers on the ellipses. Also included: bagplots, kernel density plots and outlier diagnostic plots. See the vignette.

ggmultiglyph v0.1.0: Provides ggplot2 geoms for visualizing multivariate data using glyphs. The package implements several established glyph designs described in the information visualization literature, including the review by Borgo et al. (2013).

grip v0.2.0: Implements GRIP multiscale graph layout with a unified choice between hop-count and geometry-aware edge-length graph metrics in 2D and 3D. Provides layout scoring, candidate comparison, multiscale trace diagnostics, synthetic graph families, and advanced experimental geodesic-KK utilities for weighted-layout evaluation and polish. Based on Gajer and Kobourov (2002) and Gajer, Goodrich and Kobourov (2004). There are four vignettes including Getting Started and Weighted Graph Layouts.

orbis v0.1.0: Implements a layered grammar of graphics that compiles plots to a resolution-independent scene description and renders it through two back-ends: a self-contained SVG writer with embedded JavaScript for interactive figures (tooltips, hover highlighting, zoom, pan and legend toggling) and R’s own graphics devices for publication-quality output at any resolution. Geographic layers are first class. The layered grammar follows Wickham (2010); projections follow Snyder (1987) and, for Equal Earth, Savric, Patterson and Jenny (2019); line simplification uses Douglas and Peucker (1973); the default colour scales follow the guidance on perceptually uniform palettes of Crameri, Shephard and Heron (2020). See the vignette

To leave a comment for the author, please follow the link and comment on their blog: R Works.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Continue reading: August 2026 Top 40 New CRAN Packages

Celebrating the rOpenSci Community

rOpenSci — Tue, 29 Sep 2026 00:00:00 +0000

[This article was first published on rOpenSci - open tools for open science, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Read it in: Español.

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

Stretch your comfort zone

Early in 2018 I came across this post:

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

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

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

Be respectful and kind

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

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

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

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

Contribute in your own way

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

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

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

Attract great people

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

Keep the human in your human connections

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

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

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

Noam Ross, Yanina Bellini Saibene and me at the 2024 CZI meeting in Boston, USA.

Monica Gerber, Stefanie Butland, Sean Kross, Julia Stewart Lowndes, Kelly O’Briant and me at the 2024 Posit conference in Seattle, USA.

Luis D. Verde Arregoitia, Diana Garcia Cortes and me at the Georgia Aquarium after the 2025 Posit conference in Atlanta, USA.

Ronny A. Hernández Mora and me at Ronny’s family gathering near San José, Costa Rica in 2025.

Beatriz Milz and me after coffee and pastries in São Paulo, Brazil in 2025.

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

Two books I like are How to Win Friends and Influence People and Crucial Conversations (Third Edition): Tools for Talking When Stakes Are High. ︎

To leave a comment for the author, please follow the link and comment on their blog: rOpenSci - open tools for open science.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Continue reading: Celebrating the rOpenSci Community

X Reasons for policy professionals to get into programming in 202X

Giles — Mon, 28 Sep 2026 10:06:06 +0000

[This article was first published on Data Analytics and AI Archives - Giles, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

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

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

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

Background

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

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

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

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

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

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

It’s sunk costs all the way down

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

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

Source: rogierK @ Twitter (post no longer available)

Doing smart stuff quickly

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

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

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

Probability-based slop machines

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

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

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

“The majority of human production has always been slop. Mediocrity is not a bug of technology; it is the baseline of culture. The canon of art we revere today – those few thousands works in museums and textbooks – is the surviving tip of an immense iceberg of forgotten, derivative, or simply boring creations…”
D’Isa, Francesco (1 December 2025). “The Idea of ‘AI Slop’ Is Slop”. The Philosophical Salon (link)

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

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

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

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

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

The jagged frontier of AI-assisted public policy analysis

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

Source: Ethan Mollick, Centaurs and Cyborgs on the Jagged Frontier, One Useful Thing, link

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

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

“…The right analogy for AI is not humans, but an alien intelligence with a distinct set of capabilities and limitations. Just because it exceeds human ability at one task doesn’t mean it can do all related work at human level. Although AIs and humans can perform some similar tasks, the underlying “cognitive” processes are fundamentally different…”
Mollick, E. (2024, May 12). Superhuman? What does it mean for AI to be better than a human? And how can we tell? One Useful Thing (link).

There’s no accounting for taste

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

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

Vibing fast and slow

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

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

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

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

“…LLMs and AI are, I think, a massive amplifier. You can certainly use them to amplify the worst aspects of yourself, but you can also use them to amplify the best aspects…”
Wickham, H. (2026, July 25). y code when ai? A talk on AI and coding. Tidy Design. https://tidydesign.substack.com/p/y-code-when-ai

Summing up

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

To leave a comment for the author, please follow the link and comment on their blog: Data Analytics and AI Archives - Giles.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Continue reading: X Reasons for policy professionals to get into programming in 202X

Adding a free AI chatbot to a blogdown website

R-bloggers on Almost Random — Mon, 28 Sep 2026 00:00:00 +0000

[This article was first published on R-bloggers on Almost Random, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

I wanted to add a small “Ask AI” window to this website. Visitors should be able to ask about my projects, publications or background and get a short answer with a useful link. Running it should cost nothing. I described my website setup to ChatGPT, followed its suggestions and discussed the parts I wanted to change. That covered choosing a provider, writing the code and adjusting the window for desktop and mobile.

To leave a comment for the author, please follow the link and comment on their blog: R-bloggers on Almost Random.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Continue reading: Adding a free AI chatbot to a blogdown website

R-bloggers

maestro 1.4.0 adds custom pipeline metrics

Recording metrics with update_metric()

Retrieving metrics with get_metrics()

Metrics are per execution

Accumulating with increment = TRUE

Other recent changes

A Trip from Covariance to Correlation and Back

From covariance to correlation

From correlation to covariance

Dogs of New York

Popular Breeds

Where Breeds Live

The Diffusion of Goldendoodles

A Bayesian Network to Explore Data Center Risks

Data Center Risks

Potential Benefits

Organizing the risks and benefits

The DAG

Creating a permutation test, two ways

Something different from jackals

Give me the shuffled values

One permutation at a time

Reusing a set

Footnotes

Beyond the course: Reflections from the Bioconductor–Physalia sponsorship programme

Making training accessible

Taking new methods back to research

Learning to work with unfamiliar data

Sharing what comes next

Continuing the conversation

Ordered set partitions and formulas in permute

Same groups, fewer permutations

A permutation test with lm()

Less typing with formulas

Going back to the basics: the correlation coefficient

Assumptions for the Pearson Correlation Coefficient

Hypothesis testing

Correlation in R

Relationship to slope in linear regression

Intra-class correlation (ICC)

Validating data analyses written by coding agents

Specification

Adversarial code review

Data processing checks

Testing against simulated data with known truth

Other simulation tests

Test multiple methods

What’s next?

The Pit of Success, and What a Data Science Language Would Need to Make It

Pits of failure are everywhere in data science

What makes a pit of success

Defaults, constraints, and choice architecture

Design tools that make falling into a pit of success easy

Where the pit needs walls and an exit

What the pit cannot do

Bayesian workflow [book review]

A Collaboration Fest for research software interoperability

The challenge: making methods reusable across ecosystems

Putting the idea into practice

What we learned

What comes next

Connecting to Eureka2 using a university’s Microsoft Authenticator based login

{renv}, {rix}, rv, uvr, ir: An explosion of R package and environment managers

{renv}

Use case

How to use {renv}

How to share using {renv}

Some considerations

{rix}

Use case

How to use {rix}

How to share using {rix}

Some considerations

rv

Use case

How to use rv

How to share using rv

Some considerations

uvr

Recording metrics with `update_metric()`

Retrieving metrics with `get_metrics()`

Accumulating with `increment = TRUE`

Selecting steps with `[`