In my recent book (see below), on page 166 and earlier, I made the point that, with pairwise comparisons and, more generally, whenever simultaneous statistical tests are performed, it is necessary to provide P-values that account for the familywise error rate, i.e. the probability of committing at least one incorrect rejection within the whole family of simultaneous tests (i.e. adjusted P-values). In this respect, it may be useful to recall that, for a single non-significant test, the comparison-wise error rate (E_c) is the probability of a wrong rejection for that single test (based on a non-adjusted P-value), whereas the probability of at least one wrong rejection within a family of (k) comparisons is much higher.

With pairwise comparisons, a single test is usually based on the ratio between a difference and its standard error (a t-test), which is assumed to follow a univariate t-distribution when the null hypothesis is true. When several simultaneous t-tests are performed, the vector of all t-ratios can be assumed to follow a multivariate t-distribution under the hypothesis that the null is true for all simultaneous tests (Bretz et al., 2011). Therefore, adjusted P-values can be obtained by using the probability function of a multivariate t-distribution in place of the simple univariate t-distribution.

As an example, let us reconsider the ‘mixture’ data used in Chapter 9 of the main book. Three herbicide mixtures and an untreated control were tested for their weed-control ability against an important weed in tomato, namely Solanum nigrum. In the code below, we load the data and fit a one-way ANOVA model, using the weight of weed plants per pot as the response variable and the herbicide treatment as the explanatory factor. For the sake of simplicity, we omit the usual checks of the basic assumptions (see the main book). The ANOVA table shows that the treatment effect is significant and, therefore, we proceed to compare treatment means in a pairwise fashion. The P-values shown below do not account for the familywise error rate but only for the comparison-wise error rate; these P-values can be reproduced by using the probability function of a univariate Student’s t-distribution (pt() function in R).

library(statforbiology)
library(emmeans)
library(multcomp)
dataset <- getAgroData("mixture")
dataset$Treat <- factor(dataset$Treat)
model <- lm(Weight ~ Treat, data = dataset)
anova(model)
## Analysis of Variance Table
##
## Response: Weight
##           Df  Sum Sq Mean Sq F value    Pr(>F)
## Treat      3 1089.53  363.18  23.663 2.509e-05 ***
## Residuals 12  184.18   15.35
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
groupMeans <- emmeans(model, ~Treat)
tab <- contrast(groupMeans, method = "pairwise", adjust = "none")
tab
##  contrast                         estimate   SE df t.ratio p.value
##  Metribuzin__348 - Mixture_378        4.05 2.77 12   1.461  0.1697
##  Metribuzin__348 - Rimsulfuron_30    -7.68 2.77 12  -2.774  0.0168
##  Metribuzin__348 - Unweeded         -17.60 2.77 12  -6.352  <.0001
##  Mixture_378 - Rimsulfuron_30       -11.73 2.77 12  -4.235  0.0012
##  Mixture_378 - Unweeded             -21.64 2.77 12  -7.813  <.0001
##  Rimsulfuron_30 - Unweeded           -9.91 2.77 12  -3.578  0.0038
#
# The P-value is obtained from the univariate t distribution (two-tails test)
abst <- abs(as.data.frame(tab)$t.ratio)
2 * pt(abst, 12, lower.tail = FALSE)
## [1] 1.696785e-01 1.683167e-02 3.651239e-05 1.157189e-03 4.782986e-06
## [6] 3.794451e-03

In order to obtain familywise error rates, we should switch from the univariate to the multivariate t-distribution. For example, let’s consider the first t-ratio in the previous Code Box (t = 1.461). We should ask ourselves: “what is the probability of obtaining a t-ratio as extreme as, or more extreme than, 1.461 from a multivariate t-distribution with six dimensions (i.e., the number of simoultaneous tests)?”. In this calculation, we must also consider that the 6 tests are correlated, at least to some extent, because they share some common elements, for example, the same error term in the denominator. In the simplest case (homoscedasticity and balanced data), this correlation is equal to 0.5 for all pairwise comparisons.

In earlier times, when the computing power was limited, calculating probabilities from the multivariate t-distribution was a daunting task. However, for some specific cases (e.g., linear models with homoscedastic and balanced data), adjusted P-values could be obtained by exploiting the distribution of the Studentised Range (the so-called ‘tukey’ method), which is the default option in the contrast() function of the emmeans package, as shown in the following Code box.

tab <- contrast(groupMeans, method = "pairwise")
# tab <- contrast(groupMeans, method = "pairwise", adjust = "tukey") # same as above
tab
##  contrast                         estimate   SE df t.ratio p.value
##  Metribuzin__348 - Mixture_378        4.05 2.77 12   1.461  0.4885
##  Metribuzin__348 - Rimsulfuron_30    -7.68 2.77 12  -2.774  0.0698
##  Metribuzin__348 - Unweeded         -17.60 2.77 12  -6.352  0.0002
##  Mixture_378 - Rimsulfuron_30       -11.73 2.77 12  -4.235  0.0055
##  Mixture_378 - Unweeded             -21.64 2.77 12  -7.813  <.0001
##  Rimsulfuron_30 - Unweeded           -9.91 2.77 12  -3.578  0.0173
##
## P value adjustment: tukey method for comparing a family of 4 estimates
# The P-value is obtained from the Studentised Range Distribution (two-tails test)
abst <- abs(as.data.frame(tab)$t.ratio)
ptukey(sqrt(2) * abst, 4, 12, lower.tail = FALSE)
## [1] 4.884620e-01 6.981178e-02 1.853807e-04 5.501451e-03 2.473776e-05
## [6] 1.725725e-02

This simple method yields exact familywise error rates with balanced data—which represent the vast majority of designed field experiments in agriculture—and performs reasonably well in the presence of small degrees of imbalance. Within the framework of traditional multiple-comparison testing procedures, the approach described above leads to the same results as Tukey’s HSD for balanced data and the Tukey–Kramer test for unbalanced data.

More recently, it has become possible to directly calculate probabilities from the multivariate t-distribution, which is particularly convenient because it provides a more general approach to obtaining familywise error rates. This distribution is implemented in the ‘mvtnorm’ package through the pmvt() function. To perform the calculation, we must specify, for each dimension, the interval over which the probability is to be computed (in this case, for the first t-ratio, the interval is (pm 1.461081)), the number of degrees of freedom (12), and the correlation matrix of the linear combinations, which can be directly retrieved from the ‘emmGrid’ object. The code below illustrates these calculations. The quantity ‘plev’ represents the probability of sampling within the interval (i.e. none of the six null hypotheses is wrongly rejected), whereas the familywise error rate corresponds to the probability of sampling outside the interval (i.e. at least one null hypothesis is wrongly rejected), which is obtained by subtraction.

library(mvtnorm)
t1 <- abs(as.data.frame(tab)$t.ratio)[1]
ncontr <- 6
corMat <- cov2cor(vcov(tab))
plev <- pmvt(lower = rep(-t1, ncontr), upper=rep(t1, ncontr), df = 12,
     corr = corMat)[1]
1 - plev
## [1] 0.4883843

In R, such an approach can be obtained by using the adjust = "mvt" argument.

tab <- contrast(groupMeans, method = "pairwise", adjust = "mvt")
tab
##  contrast                         estimate   SE df t.ratio p.value
##  Metribuzin__348 - Mixture_378        4.05 2.77 12   1.461  0.4885
##  Metribuzin__348 - Rimsulfuron_30    -7.68 2.77 12  -2.774  0.0698
##  Metribuzin__348 - Unweeded         -17.60 2.77 12  -6.352  0.0002
##  Mixture_378 - Rimsulfuron_30       -11.73 2.77 12  -4.235  0.0054
##  Mixture_378 - Unweeded             -21.64 2.77 12  -7.813  <.0001
##  Rimsulfuron_30 - Unweeded           -9.91 2.77 12  -3.578  0.0172
##
## P value adjustment: mvt method for 6 tests

The above function employs numerical integration methods and is based on simulation; consequently, the results are not fully reproducible. However, it is easy to see that these results are asymptotically equivalent to those obtained with the Tukey adjustment method shown above. Owing to this intrinsic complexity, the use of the adjust = "mvt" argument is not recommended for pairwise comparisons in balanced experiments, whereas it may prove useful in other situations, for example in the presence of strongly unbalanced data.

Thanks for reading—and don’t forget to check out my new book below!

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

We’re accepting abstracts for AI in Production until 23rd January. The conference takes place on 4th–5th June 2026 in Newcastle, with talks on Friday 5th across two streams: one focused on engineering and production systems, the other on machine learning and model development.

We often hear: “My work isn’t ready to talk about yet” or “I’m not sure anyone would be interested.” We want to address that hesitation directly.

Speaking at a conference isn’t primarily about promoting yourself or your organisation.

It’s a practical tool that helps you do better work. Preparing and delivering a talk forces useful reflection, invites feedback from people facing similar challenges, and turns knowledge that lives only in your head into something your team can reuse.

If you’re wondering whether your work qualifies: internal systems count, work in progress counts, partial success counts.

Submit your abstract by 23rd January on the AI in Production website.

Preparing a Talk Clarifies Your Decisions

When you sit down to explain a technical choice to an audience, you have to answer questions you might have glossed over at the time: Why did we build it this way? What constraints shaped our approach? What would we do differently now?

This isn’t about justifying your decisions to others. It’s about understanding them yourself. The process of turning a production system into a coherent narrative forces you to see patterns you were too close to notice while building it. You identify what worked, what didn’t, and why. That clarity is valuable whether or not you ever give the talk.

Many practitioners find that writing an abstract or outline reveals gaps in their thinking. A deployment strategy that seemed obvious in context becomes harder to explain without it. A monitoring approach that felt pragmatic reveals underlying assumptions. This friction is useful. It means you’re learning something about your own work.

Speaking Invites Useful Feedback

The audience at AI in Production will broadly fall across two streams: engineering (building, shipping, maintaining, and scaling systems) and machine learning (model development, evaluation, and applied ML).

Whether you’re working on infrastructure and deployment or on training pipelines and model behaviour, you’ll be in a room with people facing similar constraints: limited resources, shifting requirements, imperfect data, and operational pressures.

When you share what you’ve tried, you get feedback from people who understand the context. Someone has solved a similar problem differently. Someone has run into the same failure mode. Someone asks a question that makes you reconsider an assumption.

This kind of peer feedback is hard to get otherwise. Your team is too close to the work. Online discussions lack context. A conference talk puts your approach in front of people who can offer informed perspectives without having to understand your entire stack or organisational structure first.

Talks Help Share Responsibility and Knowledge

In many teams, knowledge about production systems sits with one or two people. They know why certain decisions were made, where the edge cases are, and how to interpret the monitoring dashboards. That concentration of knowledge creates risk.

Preparing a talk is a forcing function for documentation. To explain your system to strangers, you have to articulate what’s currently tacit. That articulation becomes something your team can use: onboarding material, decision records, runbooks.

Speaking also distributes responsibility. When you present work publicly, it stops being just yours. Your team shares ownership of the ideas. Others can critique, extend, or maintain them. This is particularly valuable for platform teams or infrastructure work, where the people who built something may not be the ones operating it six months later.

Turning Tacit Knowledge into Reusable Material

Much of what you know about your production systems isn’t written down. You understand the failure modes, the workarounds, and the operational quirks. You know which metrics matter and which are noise. You remember why you made certain tradeoffs.

A conference talk is an excuse to capture that knowledge. The slides become a reference. The abstract becomes a design document. The Q&A reveals what wasn’t clear and needs better documentation.

Even if the talk itself is ephemeral, the process of preparing it leaves artefacts. You’ve already done the hard work of running the system. Speaking about it turns that experience into something others can learn from, and you can build on.

Your Work Is Worth Sharing

If you’re maintaining AI systems in production, you’re solving problems worth talking about. Making models reliable under load, keeping training pipelines maintainable, monitoring behaviour when ground truth is delayed or absent, and managing technical debt while shipping features.

These are the problems practitioners face every day. Your approach won’t be perfect, and that’s the point. Talks about work in progress, about things that didn’t work, about compromises made under constraint are often more useful than polished success stories.

We’re looking for honest accounts of how people are actually building and operating AI systems. That might fit the engineering stream (deployment, infrastructure, monitoring, scaling) or the machine learning stream (training, evaluation, model behaviour, responsible data use). If you’re doing work in either area, you have something to contribute.

Submit an Abstract

The deadline is 23rd January. You’ll need a title and an abstract of up to 250 words. You don’t need a perfect story or a finished project. You need a problem you’ve worked on, some approaches you’ve tried, and some lessons you’ve learned.

Think about what would be useful for someone six months behind you on a similar path. Think about what you wish someone had told you before you started. Think about the conversation you’d want to have with peers who understand the constraints you’re working under.

If you’re not sure where to start, consider writing about one decision that shaped your system, one assumption that turned out to be wrong, or one constraint that changed your design. Good abstracts often start with a specific moment or choice rather than a broad overview.

Ready to submit? The deadline is 23rd January. Share one decision, one lesson, or one constraint from your production work:
https://jumpingrivers.com/ai-production/

If you have questions about whether your work fits the conference, reach out at events@jumpingrivers.com. We’re here to help make this easier.

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

Overview

The easiest way to use the freeCount R Shiny applications online is through Posit Connect Cloud, which is an online platform that simplifies the deployment of data applications and documents.

freeCount

The freeCount analysis framework provides a modular set of tools and tutorials for a structured approach to biological count data analysis. Users are guided through common data assessment, processing and analysis approaches.

The different analysis tools currently available include:

• Differential Expression (DE) Analysis – DA
• Network Analysis – NA
• Functional Analysis – FA
• Set Operations – SO

Steps

The following steps will walk you through how to run the freeCount apps online using Posit Connect Cloud.

1. Navigate to https://connect.posit.cloud/elizabethbrooks?search=freeCount
2. Select the app that you want to run

   Wait for the app to launch.

   Done! Now you are able to perform the selected analysis.

Step 1

Navigate to https://connect.posit.cloud/elizabethbrooks?search=freeCount

Step 2

Select the app that you want to run and click the name or image to open.

Wait…

Wait for the project to deploy in your Posit Cloud workspace.

Done!

Now you are able to perform the selected analysis.

Analysis Tutorials

The freeCount apps provide a set of common tools for analyzing biological data, including differential expression and network analysis. We have tutorials available to guide users through a structured analysis approach:

• Differential expression (DE) analysis is used to identify genes driving the patterns of variation associated with groups of samples.
• Functional analysis of DE results is useful for determining the functions of DE genes. Genes can have multiple functional annotations, so we need to determine which ones are important.
• Weighted gene co-expression network analysis (WGCNA) is used to investigate the function of genes at the system-level. In a network analysis genes with similar patterns of expression are grouped together into modules.
• Functional analysis of network results is useful for determining the functions of a set of interesting genes. These gene sets can be lists of genes produced from different analysis, including WGCNA.
• Set operations using Venn diagrams are generally very useful for comparing lists of things. Set operations are also a good way to identify unique or shared genes across sets of analysis results.

After passing through rOpenSci peer review, the distionary package is now newly available on CRAN. It allows you to make probability distributions quickly – either from a few inputs or from its built-in library – and then probe them in detail.

These distributions form the building blocks that piece together advanced statistical models with the wider probaverse ecosystem, which is built to release modelers from low-level coding so production pipelines stay human-friendly. Right now, the other probaverse packages are distplyr, allowing you to morph distributions into new forms, and famish, allowing you to tune distributions to data. Developed with risk analysis use cases like climate and insurance in mind, the same tools translate smoothly to simulations, teaching, and other applied settings.

This post highlights the top 3 features of this youngest version of distionary. Let’s start by loading the package.

library(distionary)

Feature 1: more than just Base R distributions

Of course, all the Base R distributions are available in distionary. Here’s everyone’s favourite Normal distribution.

dst_norm(0, 1)

Normal distribution (continuous)
--Parameters--
mean sd
0 1

plot(dst_norm(0, 1))

And good old Poisson.

dst_pois(3)

Poisson distribution (discrete)
--Parameters--
lambda
3

plot(dst_pois(3))

But there are additional game-changing distributions included, too.

A Null distribution, which always evaluates to NA. When you’re running an algorithm that encounters an issue, you can return a Null distribution instead of throwing an error. Even downstream evaluation steps won’t error out because the code still sees a distribution rather than a bare NA or NULL.

# Make a Null distribution.
null <- dst_null()
# Null distributions always evaluate to NA.
eval_quantile(null, at = c(0.25, 0.5, 0.75))

[1] NA NA NA

mean(null)

[1] NA

Empirical distributions, where the data are the distribution. These respect observed behaviour without forcing a specific shape, and are also commonly used as a benchmark for comparison against other models. Here’s an example using the Ozone concentration from the airquality dataset that comes loaded with R.

# Empirical distribution of Ozone from the `airquality` dataset.
emp <- dst_empirical(airquality$Ozone, na_action_y = "drop")
# Inspect
print(emp, n = 5)

Finite distribution (discrete)
--Parameters--
# A tibble: 67 × 2
outcomes probs
<int> <dbl>
1 1 0.00862
2 4 0.00862
3 6 0.00862
4 7 0.0259
5 8 0.00862
# ℹ 62 more rows

Compare its cumulative distribution function (CDF) to that of a Gamma distribution fitted to the Ozone levels, borrowing the probaverse’s famish package for the fitting task.

# Fit a Gamma distribution to Ozone using the famish package.
library(famish)
gamma <- fit_dst_gamma(airquality$Ozone, na_action = "drop")

# Plot the cumulative distribution functions (CDFs) together.
plot(emp, "cdf", n = 1000, xlab = "Ozone Levels (ppb)")
plot(gamma, "cdf", add = TRUE, col = "red")
legend(
 "bottomright",
 legend = c("Empirical", "Fitted Gamma"),
 col = c("black", "red"),
 lty = 1
)

These textbook distributions become much more useful once they become building blocks for building up a system. For example, they could form predictive distributions in a machine learning context, or be related to other variables. This is what the probaverse seeks to make possible.

Feature 2: friendly towards tidy tabular workflows

First, load the tidyverse to activate tidy tabular workflows. And yes, probaverse is named after the tidyverse because it aims to be a “tidyverse for probability”.

library(tidyverse)

You can safely ignore this next chunk unless you want to see how I’m wrangling some financial data for you.

# Wrangle the stocks data frame using tidyverse.
stocks <- as_tibble(EuStockMarkets) |>
 mutate(across(everything(), (x) 100 * (1 - x / lag(x)))) |>
 drop_na()

The stocks data I’ve wrangled is a table of daily percent loss for four major European stock indices. The dates don’t matter for this example, so they’ve been omitted.

stocks

# A tibble: 1,859 × 4
DAX SMI CAC FTSE
<dbl> <dbl> <dbl> <dbl>
1 0.928 -0.620 1.26 -0.679
2 0.441 0.586 1.86 0.488
3 -0.904 -0.328 0.576 -0.907
4 0.178 -0.148 -0.878 -0.579
5 0.467 0.889 0.511 0.720
6 -1.25 -0.676 -1.18 -0.855
7 -0.578 -1.23 -1.32 -0.824
8 0.287 0.358 0.193 -0.0837
9 -0.637 -1.11 -0.0171 0.522
10 -0.118 -0.437 -0.314 -1.41
# ℹ 1,849 more rows

First, let’s focus on the DAX stock index. Fit an empirical distribution like last time (notice I’m using a data mask¹ in dst_empirical() this time).

# Fit an empirical distribution to the DAX stock index.
dax <- dst_empirical(DAX, data = stocks)
# Inspect the CDF.
plot(dax, xlab = "Daily Loss (%)")

You can easily calculate some standard quantiles in tabular format so that the inputs are placed alongside the calculated outputs: just use the enframe_ prefix instead of eval_ as we did above with the Null distribution.

enframe_quantile(dax, at = c(0.25, 0.5, 0.75), arg_name = "prob")

# A tibble: 3 × 2
prob quantile
<dbl> <dbl>
1 0.25 -0.638
2 0.5 -0.0473
3 0.75 0.468

Or, more to the point here – and appealing to probaverse’s soft spot for risk-focused work – you can calculate return levels (also known as “Value at Risk” in financial applications) for specific return periods. If you don’t know what these are, they are just fancy names for quantiles.

return_periods <- c(5, 50, 100, 200, 500)
enframe_return(
 dax,
 at = return_periods,
 arg_name = "return_period",
 fn_prefix = "daily_loss_pct"
)

# A tibble: 5 × 2
return_period daily_loss_pct
<dbl> <dbl>
1 5 0.621
2 50 2.17
3 100 2.75
4 200 3.08
5 500 3.71

The tabular output becomes even more powerful when inserted into a table of models, because it facilitates comparisons and trends. To demonstrate, build a model for each stock. First, lengthen the data for this task.

# Lengthen the data using tidyverse.
stocks2 <- pivot_longer(
 stocks,
 everything(),
 names_to = "stock",
 values_to = "daily_loss_pct"
)
# Inspect
stocks2

# A tibble: 7,436 × 2
stock daily_loss_pct
<chr> <dbl>
1 DAX 0.928
2 SMI -0.620
3 CAC 1.26
4 FTSE -0.679
5 DAX 0.441
6 SMI 0.586
7 CAC 1.86
8 FTSE 0.488
9 DAX -0.904
10 SMI -0.328
# ℹ 7,426 more rows

Build a model for each stock using a group_by + summarise workflow from the tidyverse (please excuse the current need to wrap the distribution in list()). Notice that distributions become table entries, indicated here by their class <dst>.

# Create an Empirical distribution for each stock.
models <- stocks2 |>
 group_by(stock) |>
 summarise(model = list(dst_empirical(daily_loss_pct)))
# Inspect
models

# A tibble: 4 × 2
stock model
<chr> <list>
1 CAC <dst>
2 DAX <dst>
3 FTSE <dst>
4 SMI <dst>

Now you can use a tidyverse workflow to calculate tables of quantiles for each model, and expand them. In fact, this workflow is common enough that I’m considering adding a dedicated verb for it.

return_levels <- models |>
 mutate(
 df = map(
 model,
 enframe_return,
 at = return_periods,
 arg_name = "return_period",
 fn_prefix = "daily_loss_pct"
 )
 ) |>
 unnest(df) |>
 select(!model)
# Inspect
print(return_levels, n = Inf)

# A tibble: 20 × 3
stock return_period daily_loss_pct
<chr> <dbl> <dbl>
1 CAC 5 0.757
2 CAC 50 2.37
3 CAC 100 2.78
4 CAC 200 3.41
5 CAC 500 3.97
6 DAX 5 0.621
7 DAX 50 2.17
8 DAX 100 2.75
9 DAX 200 3.08
10 DAX 500 3.71
11 FTSE 5 0.542
12 FTSE 50 1.58
13 FTSE 100 2.05
14 FTSE 200 2.31
15 FTSE 500 2.87
16 SMI 5 0.552
17 SMI 50 2.03
18 SMI 100 2.52
19 SMI 200 2.91
20 SMI 500 3.55

The result is a tidy dataset that’s ready for most analyses. For example, you can easily plot a comparison of the return levels of each stock. I make these plots all the time to facilitate risk-informed decision-making.

return_levels |>
 mutate(stock = fct_reorder2(stock, return_period, daily_loss_pct)) |>
 ggplot(aes(return_period, daily_loss_pct, colour = stock)) +
 geom_point() +
 geom_line() +
 theme_bw() +
 scale_x_log10(
 "Return Period (days)",
 minor_breaks = c(1:10, 1:10 * 10, 1:10 * 100)
 ) +
 scale_y_continuous("Daily Loss", label = scales::label_number(suffix = "%")) +
 annotation_logticks(side = "b") +
 scale_colour_discrete("Stock Index")

Feature 3: make the distribution you need

You can create your own distributions with distionary by specifying only a minimal set of properties; all other properties are derived automatically and can be retrieved when needed.

Let’s say you need an Inverse Gamma distribution but it’s not available in distionary. Currently, distionary assumes you’ll at least provide the density and CDF; you could retrieve these from the extraDistr package (functions dinvgamma() and pinvgamma()). Plug them into distionary’s distribution() function and enjoy access to a variety of properties you didn’t specify, like the mean, variance, skewness, and hazard function.

# Make an Inverse Gamma distribution (minimal example).
ig <- distribution(
 density = function(x) extraDistr::dinvgamma(x, alpha = 5, beta = 20),
 cdf = function(x) extraDistr::pinvgamma(x, alpha = 5, beta = 20),
 .vtype = "continuous",
)
# Calculate anything.
mean(ig)

[1] 5

variance(ig)

[1] 8.333333

skewness(ig)

[1] 3.464085

plot(ig, "hazard", to = 20, n = 1000, xlab = "Outcome")

You might also consider giving the distribution a .name – it pays off when you’re juggling multiple models. Adding .parameters provides additional specificity to the distribution .name, but are otherwise not yet used for functional purposes.

Here is a more complete implementation of the Inverse Gamma distribution, this time implemented as a function of the two parameters. Notice I also check that the parameters are positive (cheers to the checkmate package).

dst_invgamma <- function(alpha, beta) {
 checkmate::assert_number(alpha, lower = 0)
 checkmate::assert_number(beta, lower = 0)
 distribution(
 density = (x) extraDistr::dinvgamma(x, alpha = alpha, beta = beta),
 cdf = (x) extraDistr::pinvgamma(x, alpha = alpha, beta = beta),
 quantile = (p) extraDistr::qinvgamma(p, alpha = alpha, beta = beta),
 random = (n) extraDistr::rinvgamma(n, alpha = alpha, beta = beta),
 .name = "Inverse Gamma",
 .parameters = list(alpha = alpha, beta = beta),
 .vtype = "continuous",
 )
}

Now we can make that same Inverse Gamma distribution as before:

ig2 <- dst_invgamma(5, 20)
# Inspect
ig2

Inverse Gamma distribution (continuous)
--Parameters--
alpha beta
5 20

By the way, this feature – being able to inspect other distribution properties even when they are not specified – is great for learning about probability. That’s because you can see the many ways distributions can be represented, not just by the usual density or probability mass functions seen in textbooks.

This feature also allows for extensibility of the probaverse. For example, the probaverse’s distplyr package creates mixture distributions, which do not have an explicit formula for the quantile function. However, this is not problematic – the distribution can still be defined, and distionary will figure out what the quantiles are.

What’s to come?

Currently, the distionary package provides key functionality to define and evaluate distribution objects. Future goals include:

• Broader coverage. Extending beyond univariate continuous distributions so mixed discrete/continuous and multivariate use cases feel first-class.
• Risk-native metrics. Making cost functions, tail expectations, and other decision metrics evaluate more naturally.
• Workflow ergonomics. Relaxing the requirement to supply density/CDF pairs and adding verbs that streamline common mutate–map–unnest patterns.

If this excites you, join the conversation by opening an issue or contributing.

Special thanks to the rOpenSci reviewers Katrina Brock and Christophe Dutang for insightful comments that improved this package. Also thanks to BGC Engineering Inc., the R Consortium, and the European Space Agency together with the Politecnico di Milano for supporting this project.

library(tidyverse)
library(lubridate)

event_schema <- tibble::tibble(
  match_id = character(),
  datetime = ymd_hms(character()),
  set_no = integer(),
  rally_id = integer(),
  home_team = character(),
  away_team = character(),
  team = character(),        # team performing the action
  opponent = character(),    # opponent of team
  player = character(),
  jersey = integer(),
  skill = factor(levels = c("serve","pass","set","attack","block","dig","freeball")),
  evaluation = character(),  # e.g., "error","ace","perfect","positive","negative","kill","blocked","dig"
  start_zone = integer(),    # 1..6 (or 1..9 depending system)
  end_zone = integer(),
  rotation = integer(),      # 1..6
  phase = factor(levels = c("sideout","transition")),  # derived later
  score_team = integer(),    # score for team at time of event
  score_opp  = integer(),
  point_won_by = character(), # which team won rally point
  stringsAsFactors = FALSE
)

glimpse(event_schema)

library(readr)
library(janitor)

read_events_csv <- function(path) {
  readr::read_csv(path, show_col_types = FALSE) %>%
    janitor::clean_names() %>%
    mutate(
      set_no = as.integer(set_no),
      rally_id = as.integer(rally_id),
      start_zone = as.integer(start_zone),
      end_zone = as.integer(end_zone),
      rotation = as.integer(rotation)
    )
}

normalize_names <- function(df) {
  df %>%
    mutate(
      team = str_squish(str_to_title(team)),
      opponent = str_squish(str_to_title(opponent)),
      player = str_squish(str_to_title(player)),
      evaluation = str_squish(str_to_lower(evaluation)),
      skill = factor(str_to_lower(skill),
                    levels = c("serve","pass","set","attack","block","dig","freeball"))
    )
}

# Run once inside your project
install.packages(c("renv","quarto","tidyverse","lubridate","janitor","gt","patchwork","tidymodels"))

renv::init()

# Recommended folder structure
dir.create("data/raw", recursive = TRUE, showWarnings = FALSE)
dir.create("data/processed", recursive = TRUE, showWarnings = FALSE)
dir.create("R", showWarnings = FALSE)
dir.create("reports", showWarnings = FALSE)
dir.create("figures", showWarnings = FALSE)

metric_dictionary <- tribble(
  ~metric, ~definition,
  "SO%", "Sideout percentage: points won when receiving serve / total receive opportunities",
  "BP%", "Break point percentage: points won when serving / total serving opportunities",
  "Kill%", "Kills / attack attempts",
  "Eff%", "(Kills - Errors) / attempts",
  "Ace%", "Aces / total serves",
  "Err%", "Serve errors / total serves"
)

metric_dictionary

events_raw <- read_events_csv("data/raw/events.csv")
events <- events_raw %>% normalize_names()

# Basic validation
stopifnot(all(c("match_id","set_no","rally_id","team","skill","evaluation") %in% names(events)))

# Remove obvious duplicates (same match/set/rally/team/player/skill)
events <- events %>%
  distinct(match_id, set_no, rally_id, team, player, skill, evaluation, .keep_all = TRUE)

# Ensure opponent field exists
events <- events %>%
  mutate(opponent = if_else(is.na(opponent) | opponent == "",
                            NA_character_, opponent))

# Quick data quality report
quality_report <- list(
  n_rows = nrow(events),
  n_matches = n_distinct(events$match_id),
  missing_player = mean(is.na(events$player) | events$player == ""),
  missing_zone = mean(is.na(events$start_zone)),
  skill_counts = events %>% count(skill, sort = TRUE)
)

quality_report

derive_rally_context <- function(df) {
  df %>%
    group_by(match_id, set_no, rally_id) %>%
    mutate(
      serving_team = team[which(skill == "serve")[1]],
      receiving_team = setdiff(unique(team), serving_team)[1],
      phase = case_when(
        team == receiving_team ~ "sideout",
        team == serving_team   ~ "transition",
        TRUE ~ NA_character_
      ) %>% factor(levels = c("sideout","transition"))
    ) %>%
    ungroup()
}

events <- derive_rally_context(events)

# Customize to your coding system.
eval_map <- list(
  serve = list(
    ace = c("ace"),
    error = c("error","serve_error"),
    in_play = c("in_play","good","ok","positive","negative")
  ),
  pass = list(
    perfect = c("perfect","3"),
    positive = c("positive","2","good"),
    negative = c("negative","1","poor"),
    error = c("error","0")
  ),
  attack = list(
    kill = c("kill"),
    error = c("error","attack_error"),
    blocked = c("blocked"),
    in_play = c("in_play","continuation","covered")
  )
)

is_eval <- function(x, values) tolower(x) %in% tolower(values)

serve_metrics <- events %>%
  filter(skill == "serve") %>%
  mutate(
    is_ace = is_eval(evaluation, eval_map$serve$ace),
    is_error = is_eval(evaluation, eval_map$serve$error)
  ) %>%
  group_by(match_id, team) %>%
  summarise(
    serves = n(),
    aces = sum(is_ace),
    errors = sum(is_error),
    ace_pct = aces / serves,
    err_pct = errors / serves,
    .groups = "drop"
  )

serve_metrics

pass_metrics <- events %>%
  filter(skill == "pass") %>%
  mutate(
    perfect = is_eval(evaluation, eval_map$pass$perfect),
    positive = is_eval(evaluation, eval_map$pass$positive),
    negative = is_eval(evaluation, eval_map$pass$negative),
    error = is_eval(evaluation, eval_map$pass$error),
    # A common numeric scale (0..3)
    pass_score = case_when(
      perfect ~ 3,
      positive ~ 2,
      negative ~ 1,
      error ~ 0,
      TRUE ~ NA_real_
    )
  ) %>%
  group_by(match_id, team, player) %>%
  summarise(
    passes = n(),
    perfect_pct = mean(perfect, na.rm = TRUE),
    positive_pct = mean(positive, na.rm = TRUE),
    error_pct = mean(error, na.rm = TRUE),
    avg_pass = mean(pass_score, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  arrange(desc(avg_pass), desc(passes))

pass_metrics %>% slice_head(n = 20)

attack_metrics <- events %>%
  filter(skill == "attack") %>%
  mutate(
    kill = is_eval(evaluation, eval_map$attack$kill),
    error = is_eval(evaluation, eval_map$attack$error),
    blocked = is_eval(evaluation, eval_map$attack$blocked)
  ) %>%
  group_by(match_id, team, player) %>%
  summarise(
    attempts = n(),
    kills = sum(kill),
    errors = sum(error),
    blocks = sum(blocked),
    kill_pct = kills / attempts,
    error_pct = errors / attempts,
    blocked_pct = blocks / attempts,
    eff = (kills - errors) / attempts,
    .groups = "drop"
  ) %>%
  arrange(desc(eff), desc(attempts))

attack_metrics %>% slice_head(n = 20)

defense_metrics <- events %>%
  filter(skill %in% c("block","dig")) %>%
  mutate(
    point = evaluation %in% c("stuff","kill_block","point"),
    error = evaluation %in% c("error","net","out")
  ) %>%
  group_by(match_id, team, player, skill) %>%
  summarise(
    actions = n(),
    points = sum(point),
    errors = sum(error),
    point_rate = points / actions,
    .groups = "drop"
  )

defense_metrics

rallies <- events %>%
  group_by(match_id, set_no, rally_id) %>%
  summarise(
    serving_team = team[which(skill == "serve")[1]],
    point_won_by = first(na.omit(point_won_by)),
    .groups = "drop"
  ) %>%
  mutate(
    receiving_team = if_else(point_won_by == serving_team, NA_character_, NA_character_)
  )

# Derive receiving team robustly by looking at teams in the rally
rallies <- events %>%
  group_by(match_id, set_no, rally_id) %>%
  summarise(
    teams_in_rally = list(unique(team)),
    serving_team = team[which(skill == "serve")[1]],
    point_won_by = first(na.omit(point_won_by)),
    .groups = "drop"
  ) %>%
  mutate(
    receiving_team = map2_chr(teams_in_rally, serving_team, ~ setdiff(.x, .y)[1]),
    sideout_success = point_won_by == receiving_team,
    break_point_success = point_won_by == serving_team
  )

so_bp <- rallies %>%
  pivot_longer(cols = c(serving_team, receiving_team),
               names_to = "role", values_to = "team") %>%
  group_by(match_id, team, role) %>%
  summarise(
    opps = n(),
    points = sum(if_else(role == "receiving_team", sideout_success, break_point_success)),
    pct = points / opps,
    .groups = "drop"
  ) %>%
  mutate(metric = if_else(role == "receiving_team", "SO%", "BP%")) %>%
  select(match_id, team, metric, opps, points, pct)

so_bp

first_ball_sideout <- function(df) {
  # Identify: for each rally receiving team, find the first pass and first attack.
  df %>%
    group_by(match_id, set_no, rally_id) %>%
    mutate(
      serving_team = team[which(skill == "serve")[1]],
      receiving_team = setdiff(unique(team), serving_team)[1]
    ) %>%
    ungroup() %>%
    group_by(match_id, set_no, rally_id, receiving_team) %>%
    summarise(
      pass_eval = evaluation[which(skill == "pass" & team == receiving_team)[1]],
      first_attack_eval = evaluation[which(skill == "attack" & team == receiving_team)[1]],
      point_won_by = first(na.omit(point_won_by)),
      fbso = point_won_by == receiving_team & first_attack_eval %in% c("kill"),
      .groups = "drop"
    )
}

fbso <- first_ball_sideout(events) %>%
  mutate(
    pass_bucket = case_when(
      tolower(pass_eval) %in% eval_map$pass$perfect ~ "perfect",
      tolower(pass_eval) %in% eval_map$pass$positive ~ "positive",
      tolower(pass_eval) %in% eval_map$pass$negative ~ "negative",
      tolower(pass_eval) %in% eval_map$pass$error ~ "error",
      TRUE ~ "unknown"
    )
  ) %>%
  group_by(match_id, receiving_team, pass_bucket) %>%
  summarise(
    opps = n(),
    fbso_points = sum(fbso, na.rm = TRUE),
    fbso_pct = fbso_points / opps,
    .groups = "drop"
  ) %>%
  arrange(desc(fbso_pct))

fbso

rotation_efficiency <- events %>%
  group_by(match_id, set_no, rally_id) %>%
  summarise(
    serving_team = team[which(skill == "serve")[1]],
    point_won_by = first(na.omit(point_won_by)),
    # rotation of the receiving team at first pass (common reference)
    receiving_team = setdiff(unique(team), serving_team)[1],
    receive_rotation = rotation[which(skill == "pass" & team == receiving_team)[1]],
    .groups = "drop"
  ) %>%
  group_by(match_id, receiving_team, receive_rotation) %>%
  summarise(
    opps = n(),
    so_points = sum(point_won_by == receiving_team, na.rm = TRUE),
    so_pct = so_points / opps,
    .groups = "drop"
  ) %>%
  arrange(desc(so_pct))

rotation_efficiency

# We assume "set" rows include target_zone or target_player info; if not, join from your tagging.
# This example uses end_zone as a proxy for set location (e.g., 4/2/3/back).
setter_distribution <- events %>%
  group_by(match_id, set_no, rally_id) %>%
  mutate(
    serving_team = team[which(skill == "serve")[1]],
    receiving_team = setdiff(unique(team), serving_team)[1],
    receive_pass_score = case_when(
      skill == "pass" & team == receiving_team & tolower(evaluation) %in% eval_map$pass$perfect ~ 3,
      skill == "pass" & team == receiving_team & tolower(evaluation) %in% eval_map$pass$positive ~ 2,
      skill == "pass" & team == receiving_team & tolower(evaluation) %in% eval_map$pass$negative ~ 1,
      skill == "pass" & team == receiving_team & tolower(evaluation) %in% eval_map$pass$error ~ 0,
      TRUE ~ NA_real_
    )
  ) %>%
  ungroup() %>%
  group_by(match_id, set_no, rally_id) %>%
  summarise(
    team = first(receiving_team),
    pass_score = first(na.omit(receive_pass_score)),
    set_zone = end_zone[which(skill == "set" & team == first(receiving_team))[1]],
    score_diff = (first(na.omit(score_team)) - first(na.omit(score_opp))),
    pressure = abs(score_diff) <= 2,  # "close score" proxy
    .groups = "drop"
  ) %>%
  filter(!is.na(set_zone), !is.na(pass_score)) %>%
  mutate(pass_bucket = factor(pass_score, levels = c(0,1,2,3),
                              labels = c("error","negative","positive","perfect")))

setter_distribution_summary <- setter_distribution %>%
  group_by(team, pass_bucket, pressure, set_zone) %>%
  summarise(n = n(), .groups = "drop") %>%
  group_by(team, pass_bucket, pressure) %>%
  mutate(pct = n / sum(n)) %>%
  arrange(team, pass_bucket, pressure, desc(pct))

setter_distribution_summary

library(ggplot2)

serve_zones <- events %>%
  filter(skill == "serve") %>%
  count(team, end_zone, name = "serves") %>%
  group_by(team) %>%
  mutate(pct = serves / sum(serves)) %>%
  ungroup()

ggplot(serve_zones, aes(x = factor(end_zone), y = pct)) +
  geom_col() +
  facet_wrap(~ team) +
  labs(
    title = "Serve Target Distribution by Zone",
    x = "End Zone (Serve Target)",
    y = "Share of Serves"
  )

serve_pressure <- events %>%
  group_by(match_id, set_no, rally_id) %>%
  summarise(
    serving_team = team[which(skill == "serve")[1]],
    receiving_team = setdiff(unique(team), serving_team)[1],
    serve_end_zone = end_zone[which(skill == "serve")[1]],
    pass_eval = evaluation[which(skill == "pass" & team == receiving_team)[1]],
    point_won_by = first(na.omit(point_won_by)),
    .groups = "drop"
  ) %>%
  mutate(
    pass_score = case_when(
      tolower(pass_eval) %in% eval_map$pass$perfect ~ 3,
      tolower(pass_eval) %in% eval_map$pass$positive ~ 2,
      tolower(pass_eval) %in% eval_map$pass$negative ~ 1,
      tolower(pass_eval) %in% eval_map$pass$error ~ 0,
      TRUE ~ NA_real_
    ),
    pressure = pass_score <= 1,
    ace = FALSE # if you track aces at serve level, set it here
  )

serve_pressure_summary <- serve_pressure %>%
  group_by(serving_team, serve_end_zone) %>%
  summarise(
    serves = n(),
    avg_opp_pass = mean(pass_score, na.rm = TRUE),
    pressure_rate = mean(pressure, na.rm = TRUE),
    bp_rate = mean(point_won_by == serving_team, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  arrange(desc(bp_rate))

serve_pressure_summary

attack_tendencies <- events %>%
  filter(skill == "attack") %>%
  count(team, player, start_zone, end_zone, name = "attempts") %>%
  group_by(team, player) %>%
  mutate(pct = attempts / sum(attempts)) %>%
  ungroup() %>%
  arrange(team, player, desc(pct))

attack_tendencies %>% slice_head(n = 30)

attack_with_pass <- events %>%
  group_by(match_id, set_no, rally_id) %>%
  mutate(
    serving_team = team[which(skill == "serve")[1]],
    receiving_team = setdiff(unique(team), serving_team)[1],
    pass_eval = evaluation[which(skill == "pass" & team == receiving_team)[1]]
  ) %>%
  ungroup() %>%
  filter(skill == "attack", team == receiving_team) %>%
  mutate(
    pass_bucket = case_when(
      tolower(pass_eval) %in% eval_map$pass$perfect ~ "perfect",
      tolower(pass_eval) %in% eval_map$pass$positive ~ "positive",
      tolower(pass_eval) %in% eval_map$pass$negative ~ "negative",
      tolower(pass_eval) %in% eval_map$pass$error ~ "error",
      TRUE ~ "unknown"
    ),
    kill = tolower(evaluation) %in% eval_map$attack$kill,
    error = tolower(evaluation) %in% eval_map$attack$error
  ) %>%
  group_by(team, player, start_zone, pass_bucket) %>%
  summarise(
    attempts = n(),
    kill_pct = mean(kill, na.rm = TRUE),
    eff = (sum(kill) - sum(error)) / attempts,
    .groups = "drop"
  ) %>%
  arrange(desc(eff))

attack_with_pass

shot_chart <- events %>%
  filter(skill == "attack") %>%
  mutate(
    outcome = case_when(
      tolower(evaluation) %in% eval_map$attack$kill ~ "kill",
      tolower(evaluation) %in% eval_map$attack$error ~ "error",
      tolower(evaluation) %in% eval_map$attack$blocked ~ "blocked",
      TRUE ~ "in_play"
    )
  )

ggplot(shot_chart, aes(x = factor(end_zone), fill = outcome)) +
  geom_bar(position = "fill") +
  facet_wrap(~ player) +
  labs(
    title = "Attack Outcome Mix by Target Zone (End Zone)",
    x = "Target Zone",
    y = "Share"
  )

library(broom)

# Create a rally-level modeling table
rally_model_df <- events %>%
  group_by(match_id, set_no, rally_id) %>%
  summarise(
    serving_team = team[which(skill == "serve")[1]],
    receiving_team = setdiff(unique(team), serving_team)[1],
    pass_eval = evaluation[which(skill == "pass" & team == receiving_team)[1]],
    pass_score = case_when(
      tolower(pass_eval) %in% eval_map$pass$perfect ~ 3,
      tolower(pass_eval) %in% eval_map$pass$positive ~ 2,
      tolower(pass_eval) %in% eval_map$pass$negative ~ 1,
      tolower(pass_eval) %in% eval_map$pass$error ~ 0,
      TRUE ~ NA_real_
    ),
    serve_zone = end_zone[which(skill == "serve")[1]],
    point_won_by = first(na.omit(point_won_by)),
    .groups = "drop"
  ) %>%
  filter(!is.na(pass_score), !is.na(serve_zone)) %>%
  mutate(
    sideout_success = point_won_by == receiving_team
  )

# Baseline xSO model
xso_fit <- glm(
  sideout_success ~ pass_score + factor(serve_zone),
  data = rally_model_df,
  family = binomial()
)

tidy(xso_fit)
summary(xso_fit)

rally_model_df <- rally_model_df %>%
  mutate(xSO = predict(xso_fit, type = "response"))

rally_model_df %>%
  group_by(receiving_team) %>%
  summarise(
    actual_SO = mean(sideout_success),
    expected_SO = mean(xSO),
    delta = actual_SO - expected_SO,
    .groups = "drop"
  ) %>%
  arrange(desc(delta))

# If you have event-level score columns, you can build a win probability model.
# Here we illustrate a simple logistic model from score differential and set number.

wp_df <- events %>%
  filter(!is.na(score_team), !is.na(score_opp)) %>%
  mutate(score_diff = score_team - score_opp) %>%
  group_by(match_id, set_no, rally_id) %>%
  summarise(
    team = first(team),
    score_diff = first(score_diff),
    point_won_by = first(na.omit(point_won_by)),
    .groups = "drop"
  ) %>%
  mutate(won_point = point_won_by == team)

wp_fit <- glm(won_point ~ score_diff + factor(set_no), data = wp_df, family = binomial())
wp_df <- wp_df %>%
  mutate(win_prob_point = predict(wp_fit, type = "response"))

wp_fit %>% broom::tidy()

# Minimal Elo example (team-level). You can replace with your season match table.
matches <- tibble(
  match_id = c("m1","m2","m3"),
  date = as.Date(c("2025-09-01","2025-09-05","2025-09-10")),
  home = c("Team A","Team B","Team A"),
  away = c("Team B","Team C","Team C"),
  winner = c("Team A","Team C","Team A")
)

elo_update <- function(r_home, r_away, home_won, k = 20) {
  p_home <- 1 / (1 + 10^((r_away - r_home)/400))
  s_home <- ifelse(home_won, 1, 0)
  r_home_new <- r_home + k * (s_home - p_home)
  r_away_new <- r_away + k * ((1 - s_home) - (1 - p_home))
  list(home = r_home_new, away = r_away_new, p_home = p_home)
}

teams <- sort(unique(c(matches$home, matches$away)))
ratings <- setNames(rep(1500, length(teams)), teams)

elo_log <- vector("list", nrow(matches))

for (i in seq_len(nrow(matches))) {
  m <- matches[i,]
  rH <- ratings[[m$home]]
  rA <- ratings[[m$away]]
  upd <- elo_update(rH, rA, home_won = (m$winner == m$home))
  ratings[[m$home]] <- upd$home
  ratings[[m$away]] <- upd$away
  elo_log[[i]] <- tibble(match_id = m$match_id, p_home = upd$p_home,
                         home = m$home, away = m$away,
                         winner = m$winner,
                         r_home_pre = rH, r_away_pre = rA,
                         r_home_post = upd$home, r_away_post = upd$away)
}

bind_rows(elo_log) %>% arrange(match_id)
tibble(team = names(ratings), elo = as.numeric(ratings)) %>% arrange(desc(elo))

library(stringr)

# Build simple sequences per rally: skill chain for receiving team until point ends
rally_sequences <- events %>%
  arrange(match_id, set_no, rally_id) %>%
  group_by(match_id, set_no, rally_id) %>%
  summarise(
    serving_team = team[which(skill == "serve")[1]],
    receiving_team = setdiff(unique(team), serving_team)[1],
    seq = paste(skill, collapse = "-"),
    point_won_by = first(na.omit(point_won_by)),
    .groups = "drop"
  )

# Count bigrams (transitions) from sequences
extract_bigrams <- function(seq_str) {
  tokens <- str_split(seq_str, "-", simplify = TRUE)
  tokens <- tokens[tokens != ""]
  if (length(tokens) < 2) return(tibble(from = character(), to = character()))
  tibble(from = tokens[-length(tokens)], to = tokens[-1])
}

transitions <- rally_sequences %>%
  mutate(bigrams = map(seq, extract_bigrams)) %>%
  select(match_id, bigrams) %>%
  unnest(bigrams) %>%
  count(from, to, name = "n") %>%
  group_by(from) %>%
  mutate(p = n / sum(n)) %>%
  ungroup() %>%
  arrange(from, desc(p))

transitions

library(tidymodels)

df <- rally_model_df %>%
  mutate(
    serve_zone = factor(serve_zone),
    receiving_team = factor(receiving_team)
  )

set.seed(2026)
split <- initial_split(df, prop = 0.8, strata = sideout_success)
train <- training(split)
test  <- testing(split)

rec <- recipe(sideout_success ~ pass_score + serve_zone, data = train) %>%
  step_impute_median(all_numeric_predictors()) %>%
  step_dummy(all_nominal_predictors())

model <- logistic_reg() %>%
  set_engine("glm")

wf <- workflow() %>%
  add_recipe(rec) %>%
  add_model(model)

fit <- wf %>% fit(data = train)

pred <- predict(fit, test, type = "prob") %>%
  bind_cols(test %>% select(sideout_success))

roc_auc(pred, truth = sideout_success, .pred_TRUE)
accuracy(predict(fit, test) %>% bind_cols(test), truth = sideout_success, estimate = .pred_class)

# For player/team variation, you can use lme4 (not tidymodels-native).
install.packages("lme4")
library(lme4)

# Example: include receiving_team as a random intercept
xso_glmm <- glmer(
  sideout_success ~ pass_score + factor(serve_zone) + (1 | receiving_team),
  data = rally_model_df,
  family = binomial()
)

summary(xso_glmm)

# Bayesian logistic regression with partial pooling by receiving team
install.packages("brms")
library(brms)

bayes_fit <- brm(
  sideout_success ~ pass_score + factor(serve_zone) + (1 | receiving_team),
  data = rally_model_df,
  family = bernoulli(),
  chains = 2, cores = 2, iter = 1500,
  seed = 2026
)

summary(bayes_fit)
posterior_summary(bayes_fit)

so_bp_wide <- so_bp %>%
  select(team, metric, pct) %>%
  pivot_wider(names_from = metric, values_from = pct)

so_bp_long <- so_bp %>%
  ggplot(aes(x = team, y = pct, fill = metric)) +
  geom_col(position = "dodge") +
  coord_flip() +
  labs(title = "Sideout % and Break Point % by Team", x = NULL, y = "Rate")

so_bp_long

rot_plot_df <- rotation_efficiency %>%
  mutate(receive_rotation = factor(receive_rotation, levels = 1:6))

ggplot(rot_plot_df, aes(x = receive_rotation, y = receiving_team, fill = so_pct)) +
  geom_tile() +
  labs(title = "Rotation Sideout Heatmap", x = "Rotation (Receiving)", y = "Team")

library(gt)

attack_metrics %>%
  filter(attempts >= 10) %>%
  arrange(desc(eff)) %>%
  gt() %>%
  fmt_percent(columns = c(kill_pct, error_pct, blocked_pct), decimals = 1) %>%
  fmt_number(columns = eff, decimals = 3) %>%
  tab_header(title = "Attack Leaderboard (Min 10 Attempts)")

install.packages(c("shiny","bslib"))
library(shiny)
library(bslib)
library(tidyverse)

# Assume you already computed:
# - serve_pressure_summary
# - rotation_efficiency
# - attack_tendencies

ui <- page_sidebar(
  title = "Volleyball Analytics Dashboard (R + Shiny)",
  sidebar = sidebar(
    selectInput("team", "Select Team", choices = sort(unique(serve_pressure_summary$serving_team))),
    hr(),
    helpText("Key views: serve targets, rotation sideout, attack tendencies.")
  ),
  layout_columns(
    card(
      card_header("Serve Targets by Zone"),
      plotOutput("servePlot", height = 260)
    ),
    card(
      card_header("Rotation Sideout %"),
      plotOutput("rotPlot", height = 260)
    ),
    card(
      card_header("Top Attack Tendencies"),
      tableOutput("attackTable")
    )
  )
)

server <- function(input, output, session) {

  output$servePlot <- renderPlot({
    df <- serve_pressure_summary %>% filter(serving_team == input$team)
    ggplot(df, aes(x = factor(serve_end_zone), y = bp_rate)) +
      geom_col() +
      labs(x = "Serve End Zone", y = "Break Point Rate", title = paste("Serve Effectiveness -", input$team))
  })

  output$rotPlot <- renderPlot({
    df <- rotation_efficiency %>% filter(receiving_team == input$team) %>%
      mutate(receive_rotation = factor(receive_rotation, levels = 1:6))
    ggplot(df, aes(x = receive_rotation, y = so_pct)) +
      geom_col() +
      labs(x = "Rotation", y = "Sideout %", title = paste("Rotation Sideout -", input$team))
  })

  output$attackTable <- renderTable({
    attack_tendencies %>%
      filter(team == input$team) %>%
      group_by(player) %>%
      slice_max(order_by = pct, n = 5) %>%
      ungroup() %>%
      arrange(desc(pct)) %>%
      mutate(pct = round(pct * 100, 1))
  })
}

shinyApp(ui, server)

# Create a Quarto (.qmd) file like reports/match_report.qmd
# Then render in R:
# quarto::quarto_render("reports/match_report.qmd")

# Example render call:
quarto::quarto_render(
  input = "reports/match_report.qmd",
  execute_params = list(match_id = "match_001")
)

---
title: "Match Report"
format:
  html:
    toc: true
    code-fold: show
execute:
  echo: true
  warning: false
  message: false
params:
  match_id: "match_001"
---

set.seed(2026)

bootstrap_ci <- function(x, B = 2000, conf = 0.95) {
  n <- length(x)
  boots <- replicate(B, mean(sample(x, n, replace = TRUE)))
  alpha <- (1 - conf) / 2
  quantile(boots, probs = c(alpha, 1 - alpha), na.rm = TRUE)
}

so_ci <- rallies %>%
  mutate(sideout_success = point_won_by == receiving_team) %>%
  group_by(receiving_team) %>%
  summarise(
    so = mean(sideout_success),
    ci_low = bootstrap_ci(sideout_success)[1],
    ci_high = bootstrap_ci(sideout_success)[2],
    n = n(),
    .groups = "drop"
  )

so_ci

The post Volleyball Analytics with R: The Complete Guide to Match Data, Sideout Efficiency, Serve Pressure, Heatmaps, and Predictive Models appeared first on R Programming Books.

Enjoyed learning the process of setting up a cluster of tiny PCs for parallel computing. A note to myself on installing Ubuntu, passwordless SSH, automating package installation across nodes, distributing R simulations, and comparing CV5 vs CV10 performance. Fun project!

Motivations

Part of something I want to learn this year is getting a little more into parallel computing. How we can distribute simulation computations across different devices. Lately, we have more reasons to do this because quite a few of our simulations require long running computation and leaving my laptop running overnight or several days is just not a good use it. We have also tried cloud computing as well and without knowing how those distributed cores are, well, distributed, it’s hard for me to conceptualize how these are done and what else we could optimize. Hence, what is a better way of doing it on our own! Sit tight, this is going to be a bumpy one. Let’s go!

Objectives

• Which PCs to get?
• Install Ubuntu
• Align and fix IPs
• Passwordless ssh
• Send multiple commands via ssh
  • Install R
  • Create A Template R script For Simulation
  • Install Packages On All Nodes
  • Upload Rscript to Nodes
  • Run Rscript
  • Extract data
• Compare time
• Opportunities for improvement
• Lessons learnt

Which PCs to Get?

Preferably something functional and cheap! Something like a used Lenovo M715q Tiny PCs or something similar.

Install Ubuntu

1. Download Ubuntu Server
2. Create a bootable USB using balenaEtcher
3. When starting Lenovo up, press F12 continuously until it shows an option to boot from USB. If F12 does not work, reboot and press F1 to BIOS. Go to Startup Tab, change CSM Support to Enabled. Then set Primary Boot Priority to USB by moving priority to first. Then F10 to save configuration and exit. It will then reboot to USB.
4. Make sure it’s connected to internet via LAN for smoother installation.
5. Follow the instructions to install Ubuntu, setting username, password etc. Then reboot.
6. Make sure to remove USB drive, if you didn’t it’ll remind you. Et voila!

The installations were very quick, compared to the other OS I’ve installed in the past. Very smooth as well. I thoroughly enjoyed seeting these up.

Align and Fix IPs

For organizational purpose, make sure you go to your router setting and set your computer clusters to convenient IPs such as 192.168.1.101, 192.168.1.102, 192.168.1.103 etc. You may have to reboot your computer clusters after changing it on your router.

Passwordless SSH

Next, you want to set up passwordless SSH. This is crucial for R to work!

1. Create a key

ssh-keygen -t ed25519

2. Send Copy of Key To Your Node

ssh-copy-id -i .ssh/my_key.pub username1@192.168.1.101

it will prompt you to enter your password, then after that you won’t need a pssword to ssh in.

Passwordless Sudo

This is optional. But if you’re like me, don’t want to repeat lots of typing on installation, and see if you can use bash or R to install packages, you’d need this.

ssh -t username2@192.168.1.102 'echo "$(whoami) ALL=(ALL) NOPASSWD: ALL" | sudo tee /etc/sudoers.d/$(whoami)'

It would prompt you to enter your password. You would have to do this for all your nodes

Send Multiple Commands Via SSH

Install R

for host in username1@192.168.1.101 username2@192.168.1.102 username3@192.168.1.103; do
  ssh -t $host 'sudo apt update && sudo apt install -y r-base r-base-dev'
done

This is basically installing R on all of our clusters one after another.

Create A Template R script For Simulation

Why do we do this? We want to take advantage of the multicore of each nodes as opposed to using clusters on future as the overhead network may add on to the time and makes optimization less efficiency. Instead, we will send a script to each node so that it can fork its own cores to run the simulation. Also, if we specify packages on our script, we can automate the process of installing these packages on our nodes.

library(future)
library(future.apply)
library(dplyr)
library(SuperLearner)
library(ranger)
library(xgboost)
library(glmnet)

plan(multicore, workers = 4)

set.seed(1)

n <- 10000
W1 <- rnorm(n)
W2 <- rnorm(n)
W3 <- rbinom(n, 1, 0.5)
W4 <- rnorm(n)

# TRUE propensity score model
A <- rbinom(n, 1, plogis(-0.5 + 0.8*W1 + 0.5*W2^2 + 0.3*W3 - 0.4*W1*W2 + 0.2*W4))

# TRUE outcome model
Y <- rbinom(n, 1, plogis(-1 + 0.2*A + 0.6*W1 - 0.4*W2^2 + 0.5*W3 + 0.3*W1*W3 + 0.2*W4^2))

# Calculate TRUE ATE
logit_Y1 <- -1 + 0.2 + 0.6*W1 - 0.4*W2^2 + 0.5*W3 + 0.3*W1*W3 + 0.2*W4^2
logit_Y0 <- -1 + 0 + 0.6*W1 - 0.4*W2^2 + 0.5*W3 + 0.3*W1*W3 + 0.2*W4^2

Y1_true <- plogis(logit_Y1)
Y0_true <- plogis(logit_Y0)
true_ATE <- mean(Y1_true - Y0_true)

df <- tibble(W1 = W1, W2 = W2, W3 = W3, W4 = W4, A = A, Y = Y)

tune <- list(
  ntrees = c(500,1000),
  max_depth = c(5,7),
  shrinkage = c(0.001,0.01)
)

tune2 <- list(
  ntrees = c(250, 500, 1000),
  max_depth = c(3,5,7,9),
  shrinkage = c(0.001,0.005,0.01)
)

learners <- create.Learner("SL.xgboost", tune = tune, detailed_names = TRUE, name_prefix = "xgb")
learners2 <- create.Learner("SL.xgboost", tune = tune2, detailed_names = TRUE, name_prefix = "xgb")

# Super Learner library
SL_library <- list(
  c("SL.xgboost", "SL.ranger", "SL.glm", "SL.mean"),
  c("SL.xgboost","SL.ranger"),
  c("SL.xgboost","SL.glm"),
  list("SL.ranger", c("SL.xgboost", "screen.glmnet")),
  c("SL.glmnet","SL.glm"),
  c("SL.ranger","SL.glm"),
  c(learners$names, "SL.glm"),
  c(learners$names, "SL.glmnet"),
  c("SL.gam","SL.glm"),
  c(learners2$names, "SL.glm"))

# sample
allnum <- START:END
n_sample <- length(allnum)
n_i <- 6000

# Function to run one TMLE iteration
run_tmle_iteration <- function(seed_val, df, n_i, SL_library) {
  set.seed(seed_val)
  data <- slice_sample(df, n = n_i, replace = T) |> select(Y, A, W1:W4)

  # Prepare data
  X_outcome <- data |> select(A, W1:W4) |> as.data.frame()
  X_treatment <- data |> select(W1:W4) |> as.data.frame()
  Y_vec <- data$Y
  A_vec <- data$A

  # Outcome model
  SL_outcome <- SuperLearner(
    Y = Y_vec,
    X = X_outcome,
    family = binomial(),
    SL.library = SL_library,
    cvControl = list(V = 5)
  )

  # Initial predictions
  outcome <- predict(SL_outcome, newdata = X_outcome)$pred

  # Predict under treatment A=1
  X_outcome_1 <- X_outcome |> mutate(A=1)
  outcome_1 <- predict(SL_outcome, newdata = X_outcome_1)$pred

  # Predict under treatment A=0
  X_outcome_0 <- X_outcome |> mutate(A=0)
  outcome_0 <- predict(SL_outcome, newdata = X_outcome_0)$pred

  # Bound outcome predictions to avoid qlogis issues
  outcome <- pmax(pmin(outcome, 0.9999), 0.0001)
  outcome_1 <- pmax(pmin(outcome_1, 0.9999), 0.0001)
  outcome_0 <- pmax(pmin(outcome_0, 0.9999), 0.0001)

  # Treatment model
  SL_treatment <- SuperLearner(
    Y = A_vec,
    X = X_treatment,
    family = binomial(),
    SL.library = SL_library,
    cvControl = list(V = 5)
  )

  # Propensity scores
  ps <- predict(SL_treatment, newdata = X_treatment)$pred

  # Truncate propensity scores
  ps_final <- pmax(pmin(ps, 0.95), 0.05)

  # Calculate clever covariates
  a_1 <- 1/ps_final
  a_0 <- -1/(1 - ps_final)
  clever_covariate <- ifelse(A_vec == 1, 1/ps_final, -1/(1 - ps_final))

  epsilon_model <- glm(Y_vec ~ -1 + offset(qlogis(outcome)) + clever_covariate,
                       family = "binomial")
  epsilon <- coef(epsilon_model)

  updated_outcome_1 <- plogis(qlogis(outcome_1) + epsilon * a_1)
  updated_outcome_0 <- plogis(qlogis(outcome_0) + epsilon * a_0)

  # Calc ATE
  ate <- mean(updated_outcome_1 - updated_outcome_0)

  # Calc SE
  updated_outcome <- ifelse(A_vec == 1, updated_outcome_1, updated_outcome_0)
  se <- sqrt(var((Y_vec - updated_outcome) * clever_covariate +
                   updated_outcome_1 - updated_outcome_0 - ate) / n_i)

  return(list(ate = ate, se = se))
}

# Run iterations in parallel
for (num in 1:length(SL_library)) {
  if (num %in% c(1:9)) { next }
  cat(num)
  cat("TMLE iterations in parallel with 4 workers (multicore)...n")
  start_time <- Sys.time()

  results_list <- future_lapply(START:END, function(i) {
    result <- run_tmle_iteration(i, df, n_i, SL_library[[num]])
    if (i %% 100 == 0) cat("Completed iteration:", i, "n")
    return(result)
  }, future.seed = TRUE)

  end_time <- Sys.time()
  run_time <- end_time - start_time

  # Extract results
  predicted_ate <- sapply(results_list, function(x) x$ate)
  pred_se <- sapply(results_list, function(x) x$se)

  # Results
  results <- tibble(
    iteration = START:END,
    ate = predicted_ate,
    se = pred_se,
    ci_lower = ate - 1.96 * se,
    ci_upper = ate + 1.96 * se,
    covers_truth = true_ATE >= ci_lower & true_ATE <= ci_upper
  )

  # Summary stats
  summary_stats <- tibble(
    metric = c("true_ATE", "mean_estimated_ATE", "median_estimated_ATE",
               "sd_estimates", "mean_SE", "coverage_probability", "bias"),
    value = c(
      true_ATE,
      mean(predicted_ate),
      median(predicted_ate),
      sd(predicted_ate),
      mean(pred_se),
      mean(results$covers_truth),
      mean(predicted_ate) - true_ATE
    )
  )

  # Create output directory if it doesn't exist
  if (!dir.exists("tmle_results")) {
    dir.create("tmle_results")
  }

  # Save detailed results (all iterations)
  write.csv(results, paste0("tmle_results/tmle_iterations",num,".csv"), row.names = FALSE)

  # Save summary statistics
  write.csv(summary_stats, paste0("tmle_results/tmle_summary",num,".csv"), row.names = FALSE)

  # Save simulation parameters
  sim_params <- tibble(
    parameter = c("n_population", "n_sample_iterations", "n_bootstrap_size",
                  "SL_library", "n_workers", "runtime_seconds"),
    value = c(n, n_sample, n_i,
              paste(SL_library[[num]], collapse = ", "),
              4, as.numeric(run_time, units = "secs"))
  )
  write.csv(sim_params, paste0("tmle_results/simulation_parameters",num,".csv"), row.names = FALSE)

  # Save as RData for easy loading
  save(results, summary_stats, sim_params, true_ATE, file = paste0("tmle_results/tmle_results",num,".RData"))

}

What we did above is basically a template script (We are saving this as par_test_script.R), one where we can edit where to begin and end in terms of which iteration to start and end per node. And also instruction to save result. This is when we can put a little more effort in incorporating some instructions to inform us when task is completed (e.g., via email) and also it would also be nice to know what is the ETA of the entire task by perhaps benchmarking how long the first iteration took to complete, then multiple by total iters per node. Again, this can be sent via email, and also maybe only for the first node as opposed to all nodes, so we’re not bombarded with messages beginning and the end.

Install Packages On All Nodes

## List all of our nodes
my_clusters <- list(
  c("username1@192.168.1.101"),
  c("username2@192.168.1.102"),
  c("username3@192.168.1.103"))


## Grab all of the packages needed on our script
packages <- gsub("library(([^)]+))", "1",grep("^library",readLines("par_test_script.R"),value = T))

## Create function to run sudo
remote_r_sudo <- function(host, r_code, intern = FALSE) {
  escaped <- gsub('"', '\"', r_code)
  cmd <- sprintf("ssh %s 'sudo Rscript -e "%s"'", host, escaped)
  system(cmd, intern = intern)
}

## Loop over to install
for (cluster_i in my_clusters) {
  print(cluster_i)
  for (package in packages) {
  command <- sprintf('if (!require("%s")) install.packages("%s")', package, package)
  remote_r_sudo(cluster_i, command)
  }
}

Make sure your computer doesn’t go to sleep with this. If this is the first time your nodes are installing these extensive libraries, it will take a while. Another way we can do this is to use future_lapply for all nodes and also tmux for all installations so that we don’t need to rely on our local workstation to be turned on to continue with the installation. See below on how we used tmux to set and forget method.

Upload Rscript to Nodes

Alright, now we have installed the appropriate packages above, let’s upload scripts to our nodes.

Distribute Work

num_list <- list()
clust_num <- 3
total_loop <- 1000
div_iter <- total_loop/clust_num
final_iter <- total_loop #only use this for custom e.g., if one node did not work and it's in charge of 300:500, we can put 500 for this and set first_iter as 300
first_iter <- 1
last_iter <- round(div_iter,0) + first_iter

for (i in 1:clust_num) {
  if (i == clust_num) {
    num_list[[i]] <- paste0(first_iter,":",final_iter)
    next
  }
  num_list[[i]] <- paste0(first_iter,":",last_iter)
  first_iter <- round(first_iter + div_iter, 0)
  last_iter <- round(last_iter + div_iter, 0)
}

num_list

## [[1]]
## [1] "1:334"
##
## [[2]]
## [1] "334:667"
##
## [[3]]
## [1] "667:1000"

for (i in 1:length(my_clusters)) {
  username <- sub("@.*","",my_clusters[[i]])
  system(sprintf("sed 's/START:END/%s/g' par_test_script.R > par_test_script1.R & scp par_test_script1.R %s:/home/%s/par_test_script1.R",num_list[[i]],my_clusters[[i]],username))
}

We’ll iterate and insert the appropriate iters for each node and save it to par_test_script1.R. Then upload to each nodes with the code above.

Check set.seed in multicore

sample_df <- function(seed, df, n = 6000) {
  set.seed(seed)
  df_sample <- slice_sample(n = n, .data = df)
  return(df_sample)
}

future_lapply(100, function(x) sample_df(seed=x,df=df))

When we did the above on local computer and also in terminal with multicore. It’s still the same! Woo hoo!

The interesting thing is I didn’t have to set future.seed = T or future.seed = some_number for this. However, if we put a number on future.seed, it will return the reproducible data! This is great, next time I’ll just use this seed and I don’t have to use set.seed(i).

Run Rscript

for (i in 1:length(my_clusters)) {
  # set your tmux new session name, here we call it "test"
  cluster_name <- "test"

  # terminate any existing tmux with the existing name
  system(sprintf("ssh %s 'tmux kill-session -t %s 2>/dev/null || true'", my_clusters[[i]], cluster_name))

  # create new tmux session
  system(sprintf("ssh %s 'tmux new-session -d -s %s'", my_clusters[[i]], cluster_name))

  # run rscript in tmux
  system(sprintf("ssh %s 'tmux send-keys -t %s "Rscript par_test_script1.R > result_%d.txt"' ENTER",
                 my_clusters[[i]], cluster_name, i))
}

The code above is quite self-explanatory. Once the above code is run and completed, there we have it! it should be running in the background! You can do a spot check and see if it’s actually running. Once completed, we’ll extract the data.

Extract Data

Since we have 10 combinations we want to assess, we will set nums as 1:10 and get our data! On your template script you can set however you want to save your data, and for extraction, just look for those and download them, read and merge! Or however you want to do it.

nums <- 1:10
df <- tibble()
for (num in nums) {
  print(num)
for (i in 1:length(my_clusters)) {
  response <- system(sprintf("scp %s:tmle_results/simulation_parameters%d.csv simulation_parameters%d.csv", my_clusters[[i]], num, num), intern = F)
  if (response == 1) { next }
  df_i <- read_csv(paste0("simulation_parameters",num,".csv"), show_col_types = F)
  sl_i <- df_i |> filter(parameter == "SL_library") |> pull(value)
  df <- rbind(df, df_i |> mutate(method = sl_i, num = num))
}
}

df_sim_param <- df

df <- tibble()
for (num in nums) {
for (i in 1:length(my_clusters)) {
  response <- system(sprintf("scp %s:tmle_results/tmle_iterations%d.csv tmle_iterations%d.csv", my_clusters[[i]], num, num), intern = F)
  if (response == 1) { print(paste0(my_clusters[[i]]," is missing num", num)) ; next }
  df_i <- read_csv(paste0("tmle_iterations",num,".csv"), show_col_types = F) |>
    mutate(num = num)
  df <- rbind(df, df_i)
}
}

df_iter <- df

Take note that sometimes you may encounter issues, if for some reason the node is unable to complete the task, you can identify it then redistribute those tasks to the entire computer cluster.

Compare Time

Let’s take at our compute time for 1 cluster, 3 cluster with 5-fold cv, 3 cluster with 10-fold cv.

Looking at the time, we can definitely see the improvement of time from 1 cluster to 3 cluster. Take a look at our good old tuned xgboost and logistic regression, it took use previously for a quadcore 3.29 hours to complete, down to 1.8 hours. You’d imagine that if we use 3 pc’s as a cluster, we would notice improvement to ~1.1 hours, but I guess not for xgboost. Will have to investigate this. But if we look at xgboost + logistic regression without tuning, we went from 0.47 hours to 0.17 hours which made sense! Very interesting. Now if we up our CV to 10 fold, then we see that it took longer (makes senses), but still better than using 1 quadcore. I’ve heard people said that if you increase your K-fold CV, you reduce your bias, but increase variance. Let’s see if that’s true in our case here.

Wow, not too shabby! Indeed when we went from cv5 to cv10, we have reduced bias and slightly increased variance! How about that. Everything except gam + lr, which make sense because we don’t really tune them. Though that being said, I wonder what’s under the hood that controls the knot for gam in superlearner. Will need to check that out. With this, it looks like tuned xgboost + lr might have the best numbers. Well, now we’ve seen bias and variance, what about coverage?

5-fold CV

library(tidyverse)

num_df <- sim_param_cv5_clus5 |>
  select(num, method)

df_coverage <- df_iter_cv5_clus3 |>
  group_by(num) |>
  arrange(ate) |>
  mutate(iter = row_number()) |>
  mutate(cover = case_when(
    covers_truth == F & ate < true_ATE ~ "right_missed",
    covers_truth == F & ate > true_ATE ~ "left_missed",
    covers_truth == T ~ "covered"
  )) |>
  select(num, cover) |>
  group_by(num, cover) |>
  tally() |>
  ungroup(cover) |>
  mutate(prop = n*100/sum(n)) |>
  pivot_wider(id_cols = num, names_from = "cover", values_from = "prop") |>
  mutate(text = paste0("right missed: ",right_missed,"% covered: ",covered,"% left missed: ",left_missed,"%")) |>
  select(num, text)

method <- tibble(
  num = c(1:9),
  method = c("xgb + rf + lr + mean","xgb + rf","xgb + lr","rf + (xgb + preprocess w glmnet)","glmnet + lr","rf + lr","tuned xgb + lr","tuned xgb + glmnet","gam + lr")
)

plot <- df_iter_cv5_clus3 |>
  group_by(num) |>
  arrange(ate) |>
  mutate(iter = row_number()) |>
  mutate(cover = case_when(
    covers_truth == F & ate < true_ATE ~ "right_missed",
    covers_truth == F & ate > true_ATE ~ "left_missed",
    covers_truth == T ~ "covered"
  )) |>
  ggplot(aes(x=iter,y=ate,color=cover)) +
  geom_point(alpha=0.2) +
  geom_errorbar(aes(x=iter,ymin=ci_lower,ymax=ci_upper), alpha=0.2) +
  geom_hline(aes(yintercept=0.0373518), color = "blue") +
  geom_text(data = df_coverage,
            aes(x = 500, label = text),
            y = -0.05,
            inherit.aes = FALSE,
            size = 3,
            hjust = 0.5) +
  scale_color_manual(values = c("covered" = "#619CFF",
                                  "left_missed" = "#F8766D",
                                  "right_missed" = "#00BA38")) +
  theme_bw() +
  facet_wrap(.~num, ncol = 1,labeller = as_labeller(setNames(method$method, method$num))) +
  theme(legend.position = "bottom")

lr: logistic regression, xgb: xgboost, rf : random forest, gam : generalized additive model.

Wow, look at gam + lr’s assymetrical coverage! This is so true then when we’re assessing, a point estimate of coverage is not adequate to assess the global usefulness of a method. We can see that this method is very bias indeed with asymmetrical tails. Since CV5 and CV10 do not have significant difference in coverage, we’ll skip the visualization.

Opportunities for improvement

• plenty of opportunities to turn our personal project into a package that will help us
• Use parallel computing on local to run system (such as installation) since this takes a lot of time
• Write function to let us know when tasks are completed
• Write function to estimate time of completion
• Write function to redistribute missing iterations
• learn openMPI
• make a package for the functions above so I can reuse in the future

Lessons Learnt:

• used more sprintf with this learning experience when using with system.
• learn that in future_lapply in multicore future.seed=100 or whatever number will help reproduce the same data
• Made a few pipeline to install packages on multiple nodes
• learnt set.seed in multicore works fine
• observed reduced bias with increase variance from cv5 to cv10

If you like this article:

• please feel free to send me a comment or visit my other blogs
• please feel free to follow me on BlueSky, twitter, GitHub or Mastodon
• if you would like collaborate please feel free to contact me

I often hear the comment that LLMs/generative AI (large language models) can’t be trusted for research tasks.

Image Google’s Nano Banana tasked with “Generate an image of a male African researcher holding a balloon that is pulling them up above a tidal wave of AI generated slop that is full of errors. The balloon has a research paper inside of it. Generate the image in the style of a Simpsons cartoon.”

But this is the wrong way to think about LLMs. Humans also can’t be trusted to do scientific research accurately. They make mistakes. That’s why we have systems for review.

The more important question is: Are LLMs more accurate than humans at completing a given task?

I actually think LLMs might lead to better scientific coding and statistical analysis.

A common example of what LLMs get criticised for is writing code or performing statistical analyses. The LLM might hallucinate non-truths, or at least mislead you into thinking the analysis you have done is scientifically accurate.

The implication is that we should not be using them for particular tasks, like designing statistical models.

Its right to be skeptical of AI produced output. However, we also need to be skeptical of human produced output. Humans make mistakes as well.

As scientists peer-review is baked into our culture. But code review is much rarer. We also don’t have many systematic reviews of scientific coding that have quantified the rate of mistakes.

I suspect that mistakes in scientific coding are more common than we’d like to believe.

In one (rare) example, researchers reviewed population modelling analyses and found mathematical errors were common. One type of error occured in 62% of studies!

Now I haven’t set an LLM agent the task of doing the equivalent population models to see what its error rate is. However, my tests (which are under review) of agents at quite complicated stats and ecological modelling are showing 80-90% performance at accurately completing the tasks.

So the LLM agents are potentially doing better than the humans and making fewer mistakes.

Why I think LLMs might lead to better research is that they give us more time for code review.

As an ecological modeller I invest a ton of time into writing code, then checking that code works the way I want (and in a mathematically accurate way).

LLMs are now doing more of the code writing for me. Used effectively, this gives me more time to review the code for accuracy, as well as checking the code is an accurate representation of the scientific theory.

A human with an LLM partner could choose to: (1) produce crap work faster than pre-LLM, OR (2) produce higher quality work in a similar amount of time to what it took them pre-LLM.

I’m arguing that we should be aiming to produce the higher quality work. We can do this if we use LLMs to speed up code, then use the extra time for more quality assurance.

More generally, don’t get fooled by the argument that “genAI makes mistakes, so it can’t be trusted”.

Its the wrong way to think about the problem, and I think will lead us to being blind-sided by the oncoming flood of research slop created with genAI.

A better way to think about it is: “genAI and humans both make mistakes, how can we design workflows so that their strengths complement each other and we produce higher quality work”.

This will give us outcomes that are of higher quality than the pre-LLM world, and hopefully will rise above the huge quantity of AI generated slop that is currently happening.

1. Introduction: The Strategic Edge of Agentic Finance

In the contemporary landscape of quantitative finance, the bottleneck is no longer data availability, but the speed of insight generation. Leveraging the Microsoft AI Foundry ecosystem, we have moved beyond static scripting into the realm of Autonomous Financial Agents. This article explores how a specialized agent can navigate precious metal volatility by analyzing the Gold/Silver ratio with high-performance precision.

2. Infrastructure: Model Deployment on Microsoft AI Foundry

The intelligence behind this analysis is not a local script but a deployed model instance on Microsoft AI Foundry. We utilize the GPT-4o model, deployed as a scalable web service within the Foundry environment.

• Endpoint Security: By using the Azure OpenAI service within AI Foundry, we ensure that financial queries and data remain within a secure, enterprise-grade perimeter.
• Agentic Logic: The “Agent” is more than just a model; it is a programmed entity with a specific System Prompt that defines its persona as a Quantitative Researcher. This allows the model to “reason” through the necessary steps: from library loading to data merging and final visualization.

3. The Technical Bridge: Python-R Integration

One of the most powerful features of our AI Foundry Agent is its multi-lingual capability. It bridges the gap between Python and R using the rpy2 library, creating a high-performance research pipeline.

The R Ecosystem in Play:

• tidyquant & timetk: These packages are the engine for our time-series analysis. tidyquant handles the seamless fetching of GC=F and SI=F data, while timetk manages the complex task of plotting with built-in smoothing algorithms.
• dplyr & lubridate: Essential for the “tidy” manipulation of data, allowing the Agent to perform inner_join operations and date-based filtering with surgical precision.

4. Methodology: Taming the Noise with Visual Precision

To extract actionable trends, the Agent is instructed to apply a LOESS smoothing algorithm. By strictly setting .line_size = 1.5 and .smooth_size = 1.5, we ensure the trendline is bold enough to be the primary focus for analysts, effectively “taming” the daily price volatility.

5. Conclusion: Scaling Quantitative Research

The synergy between Microsoft AI Foundry, deployed LLMs, and specialized R packages represents the future of financial research. We have replaced manual data wrangling with an autonomous, standardized agent that can be scaled across thousands of different asset pairs with a single command.

The ABI Connection (Bridging Python to R in VS Code)

For the script to run locally in VS Code, we must establish a robust Application Binary Interface (ABI) connection. This is handled by the rpy2 library, which serves as the translation layer between Python and the R interpreter.

• Synchronization: The script uses a localconverter to transform Python data types into R objects in real-time.
• Environment Sync: Before the Agent’s code is executed, the script automatically synchronizes the working directory (setwd) so that files generated by R (like the ratio_plot.png) are immediately accessible to the Python environment for rendering.

import os

# Force rpy2 to use ABI mode to avoid the Windows CFFI conflict
os.environ['RPY2_CFFI_MODE'] = 'ABI'

import rpy2.robjects as robjects
from rpy2.robjects.packages import importr
print("Interface initialized in ABI mode.")

The Integrated Agent Script:

import os
import httpx
from openai import AzureOpenAI
import rpy2.robjects as robjects
from rpy2.robjects import pandas2ri
from rpy2.robjects.conversion import localconverter
from IPython.display import Image, display

#Microsoft AI Foundry - Azure OpenAI Connection
client = AzureOpenAI(
    api_version="2024-12-01-preview",
    azure_endpoint="AZURE_OPENAI_ENDPOINT",
    api_key="AZURE_OPENAI_KEY",
    http_client=httpx.Client(verify=False, trust_env=False)
)

def run_updated_agent(user_request):
    system_instructions = (
        "You are a Quantitative Researcher. MANDATORY: All output, comments, and labels in English. "
        "Strict Operational Guidelines:n"
        "1. Libraries: library(tidyquant), library(timetk), library(lubridate), library(dplyr), library(ggplot2).n"
        "2. Analysis: Fetch GC=F and SI=F for 3 years, merge via inner_join, and calculate 'ratio'.n"
        "3. Visualization: Use timetk::plot_time_series with .interactive = FALSE and .title = "Gold/Silver Ratio".n"
        "4. Precision: Set .line_size = 2 and ALWAYS set .smooth_size = 2 for the smoothing line.n"
        "5. Set title font face and axis texts font face to 'bold', and size to 16 with theme() function.n"
        "6. EXPORT: Save using 'ggsave("ratio_plot.png", width = 10, height = 6)'.n"
        "7. Output ONLY raw R code."
    )

    response = client.chat.completions.create(
        model="gpt-4o",
        messages=[
            {"role": "system", "content": system_instructions},
            {"role": "user", "content": user_request}
        ]
    )

    # Cleaning any markdown or headers to get raw code
    agent_code = response.choices[0].message.content.strip()
    if agent_code.startswith("```"):
        agent_code = "n".join(agent_code.split("n")[1:-1])

    print("-" * 40)
    print(agent_code)
    print("-" * 40)

    try:
        with localconverter(robjects.default_converter + pandas2ri.converter):
            # Synchronize working directory
            robjects.r(f'setwd("{os.getcwd().replace("", "/")}")')
            robjects.r(agent_code)

            if os.path.exists("ratio_plot.png"):
                display(Image(filename="ratio_plot.png"))
    except Exception as e:
        print(f"Agent Error: {e}")

# Execution
run_updated_agent("Plot the Gold/Silver ratio for the last 3 years with a smooth line.")

January 12, 2026

I’m back with the Oscars Best Picture model, albeit a little late. I had a busy holiday season, but the story of December was surprising: The Secret Agent was the favorite, followed by One Battle After Another. This was largely due to The Secret Agent’s runtime, which is right in the sweet spot for Best Picture winners.

However, remember from the last two years that my model is assuming these movies have been nominated for Best Picture. The biggest barrier facing The Secret Agent is being nominated; it is a non-English language film. While the only non-English language film to win has been Parasite (2019), the nominations are rare enough that, given that the film has been nominated, having no English dialogue isn’t a barrier to winning. The DGA and PGA nominations make me think The Secret Agent won’t be nominated, however.

It also surprised me that One Battle After Another wasn’t favored more, given Paul Thomas Anderson is a generationally phenomenal writer-director, but none of his films have won Best Picture. Looking into the data, it looks like the “career award” is not much of a thing for Best Picture (like it seems to be for the acting and directing categories). Just the opposite: If a director has had a film nominated or won Best Picture before, it actually hurts their chances of winning in my model.

That was then, this is now, though. No more awards in my models will name nominees or winners before the Oscar nominations. Where do we stand going into the announcement?

One Battle After Another is the favorite, at about 15% chance of winning. Following it closely is The Secret Agent (10%), followed by Marty Supreme (9%), Hamnet (8%), Wicked: For Good (8%), and Frankenstein (8%).

You can read more about the details of the model from my post last year and the year before. The big change I’ve made here is calibrate the probabilities so that the model isn’t too sure of itself. I will see you on the other side of the Oscar nominations.

I was recently reading Davis Vaughan’s blog post Semi-automating 200 Pull Requests with Claude Code and it really resonated with me, as I’ve been using LLMs for tedious tasks like that for some time now. Davis’s key insight: structure = success. When you can scope a task tightly and provide clear context, LLMs become genuinely useful tools.

If you’ve been following my work, you know that reproducible pipelines have been my main focus for some time now. It’s the reason I wrote {rix} for reproducible R environments, {rixpress} for declarative pipelines, and even a Python port called ryxpress. I genuinely believe these tools make data science better: more reproducible, more debuggable, more shareable.

But I also know that getting people to adopt new tools is hard. Learning a new way of structuring your code takes time and effort, and most people are busy enough already. Here’s where LLMs enter the picture: they can help translate your existing scripts into this more structured format. You provide your monolithic script, explain what you want, and the LLM does the grunt work of restructuring it.

The typical way we write analytics scripts (long chains of %>% calls in R or method-chaining in Python) works fine for interactive exploration, but quickly turns into spaghetti that’s hard to modify, test, or debug. Take my old Luxembourg Airport analysis as an example: it works, but turning that kind of script into a proper pipeline with caching, explicit dependencies, and testability is tedious work.

But we’re in 2026 where LLMs now make this trivial.

library(rixpress)

# Step 0: Load the data
avia <- rxp_r_file(
  name = avia,
  path = "avia_par_lu.tsv",
  read_function = readr::read_tsv
)

# Step 1: Select and reshape (wide → long)
avia_long <- rxp_r(
  name = avia_long,
  expr =
    avia %>%
      select("unit,tra_meas,airp_prtime", contains("20")) %>%
      gather(date, passengers, -`unit,tra_meas,airp_prtime`)
)

# Step 2: Split composite key column
avia_split <- rxp_r(
  name = avia_split,
  expr =
    avia_long %>%
      separate(
        col = `unit,tra_meas,airp_prtime`,
        into = c("unit", "tra_meas", "air_prtime"),
        sep = ","
      )
)

# Step 3: Recode transport measure
avia_recode_tra_meas <- rxp_r(
  name = avia_recode_tra_meas,
  expr =
    avia_split %>%
      mutate(
        tra_meas = fct_recode(
          tra_meas,
          `Passengers on board` = "PAS_BRD",
          `Passengers on board (arrivals)` = "PAS_BRD_ARR",
          `Passengers on board (departures)` = "PAS_BRD_DEP",
          `Passengers carried` = "PAS_CRD",
          `Passengers carried (arrival)` = "PAS_CRD_ARR",
          `Passengers carried (departures)` = "PAS_CRD_DEP",
          `Passengers seats available` = "ST_PAS",
          `Passengers seats available (arrivals)` = "ST_PAS_ARR",
          `Passengers seats available (departures)` = "ST_PAS_DEP",
          `Commercial passenger air flights` = "CAF_PAS",
          `Commercial passenger air flights (arrivals)` = "CAF_PAS_ARR",
          `Commercial passenger air flights (departures)` = "CAF_PAS_DEP"
        )
      )
)

# Step 4: Recode unit
avia_recode_unit <- rxp_r(
  name = avia_recode_unit,
  expr =
    avia_recode_tra_meas %>%
      mutate(
        unit = fct_recode(
          unit,
          Passenger = "PAS",
          Flight = "FLIGHT",
          `Seats and berths` = "SEAT"
        )
      )
)

# Step 5: Recode destination
avia_recode_destination <- rxp_r(
  name = avia_recode_destination,
  expr =
    avia_recode_unit %>%
      mutate(
        destination = fct_recode(
          `air_prtime`,
          `WIEN-SCHWECHAT` = "LU_ELLX_AT_LOWW",
          `BRUSSELS` = "LU_ELLX_BE_EBBR",
          `GENEVA` = "LU_ELLX_CH_LSGG",
          `ZURICH` = "LU_ELLX_CH_LSZH",
          `FRANKFURT/MAIN` = "LU_ELLX_DE_EDDF",
          `HAMBURG` = "LU_ELLX_DE_EDDH",
          `BERLIN-TEMPELHOF` = "LU_ELLX_DE_EDDI",
          `MUENCHEN` = "LU_ELLX_DE_EDDM",
          `SAARBRUECKEN` = "LU_ELLX_DE_EDDR",
          `BERLIN-TEGEL` = "LU_ELLX_DE_EDDT",
          `KOBENHAVN/KASTRUP` = "LU_ELLX_DK_EKCH",
          `HURGHADA / INTL` = "LU_ELLX_EG_HEGN",
          `IRAKLION/NIKOS KAZANTZAKIS` = "LU_ELLX_EL_LGIR",
          `FUERTEVENTURA` = "LU_ELLX_ES_GCFV",
          `GRAN CANARIA` = "LU_ELLX_ES_GCLP",
          `LANZAROTE` = "LU_ELLX_ES_GCRR",
          `TENERIFE SUR/REINA SOFIA` = "LU_ELLX_ES_GCTS",
          `BARCELONA/EL PRAT` = "LU_ELLX_ES_LEBL",
          `ADOLFO SUAREZ MADRID-BARAJAS` = "LU_ELLX_ES_LEMD",
          `MALAGA/COSTA DEL SOL` = "LU_ELLX_ES_LEMG",
          `PALMA DE MALLORCA` = "LU_ELLX_ES_LEPA",
          `SYSTEM - PARIS` = "LU_ELLX_FR_LF90",
          `NICE-COTE D'AZUR` = "LU_ELLX_FR_LFMN",
          `PARIS-CHARLES DE GAULLE` = "LU_ELLX_FR_LFPG",
          `STRASBOURG-ENTZHEIM` = "LU_ELLX_FR_LFST",
          `KEFLAVIK` = "LU_ELLX_IS_BIKF",
          `MILANO/MALPENSA` = "LU_ELLX_IT_LIMC",
          `BERGAMO/ORIO AL SERIO` = "LU_ELLX_IT_LIME",
          `ROMA/FIUMICINO` = "LU_ELLX_IT_LIRF",
          `AGADIR/AL MASSIRA` = "LU_ELLX_MA_GMAD",
          `AMSTERDAM/SCHIPHOL` = "LU_ELLX_NL_EHAM",
          `WARSZAWA/CHOPINA` = "LU_ELLX_PL_EPWA",
          `PORTO` = "LU_ELLX_PT_LPPR",
          `LISBOA` = "LU_ELLX_PT_LPPT",
          `STOCKHOLM/ARLANDA` = "LU_ELLX_SE_ESSA",
          `MONASTIR/HABIB BOURGUIBA` = "LU_ELLX_TN_DTMB",
          `ENFIDHA-HAMMAMET INTERNATIONAL` = "LU_ELLX_TN_DTNH",
          `ENFIDHA ZINE EL ABIDINE BEN ALI` = "LU_ELLX_TN_DTNZ",
          `DJERBA/ZARZIS` = "LU_ELLX_TN_DTTJ",
          `ANTALYA (MIL-CIV)` = "LU_ELLX_TR_LTAI",
          `ISTANBUL/ATATURK` = "LU_ELLX_TR_LTBA",
          `SYSTEM - LONDON` = "LU_ELLX_UK_EG90",
          `MANCHESTER` = "LU_ELLX_UK_EGCC",
          `LONDON GATWICK` = "LU_ELLX_UK_EGKK",
          `LONDON/CITY` = "LU_ELLX_UK_EGLC",
          `LONDON HEATHROW` = "LU_ELLX_UK_EGLL",
          `LONDON STANSTED` = "LU_ELLX_UK_EGSS",
          `NEWARK LIBERTY INTERNATIONAL, NJ.` = "LU_ELLX_US_KEWR",
          `O.R TAMBO INTERNATIONAL` = "LU_ELLX_ZA_FAJS"
        )
      )
)

# Step 6: Final cleaned dataset
avia_clean <- rxp_r(
  name = avia_clean,
  expr =
    avia_recode_destination %>%
      mutate(passengers = as.numeric(passengers)) %>%
      select(unit, tra_meas, destination, date, passengers)
)

# Step 7: Quarterly arrivals
avia_clean_quarterly <- rxp_r(
  name = avia_clean_quarterly,
  expr =
    avia_clean %>%
      filter(
        tra_meas == "Passengers on board (arrivals)",
        !is.na(passengers),
        str_detect(date, "Q")
      ) %>%
      mutate(date = yq(date))
)

# Step 8: Monthly arrivals
avia_clean_monthly <- rxp_r(
  name = avia_clean_monthly,
  expr =
    avia_clean %>%
      filter(
        tra_meas == "Passengers on board (arrivals)",
        !is.na(passengers),
        str_detect(date, "M")
      ) %>%
      mutate(date = ymd(paste0(date, "01"))) %>%
      select(destination, date, passengers)
)

# Populate and build the pipeline
rxp_populate(
  list(
    avia,
    avia_long,
    avia_split,
    avia_recode_tra_meas,
    avia_recode_unit,
    avia_recode_destination,
    avia_clean,
    avia_clean_quarterly,
    avia_clean_monthly
  )
)

rxp_make()

library(rix)

rix(
  date = "2026-01-10",
  r_pkgs = c(
    "readr",
    "dplyr",
    "tidyr",
    "forcats",
    "lubridate",
    "stringr",
    "rixpress"
  ),
  ide = "none",
  project_path = ".",
  overwrite = TRUE
)

Read it in: Español.

We are pleased to announce the opening of a new call for applications for the rOpenSci Champions Program in Spanish, which will begin in 2026. We will be accepting applications beginning in January 12, 2026 and until February 20, 2026 for both the roles of Champion as well as Mentor.

As in the previous cohort, the 2026 program will be developed entirely in Spanish and will have a regional focus on Latin America with the objective of further strengthening the research and open science software in this region.

Key dates of the 2026 call

• Opening of the call: Monday, January 12, 2026
• Community Call: Research Software in Latin America Wednesday, January 21, 2026. In Spanish
• Application Clinic: Thursday, February 5, 2026
• Closing date: Friday, February 20, 2026

The Community Call on January 21 will feature the participation of Champions and Mentors from previous cohorts. They will share their experiences and answer questions about the program, and we invite you to join us! See the recording of last year’s event.

During February we will hold one application clinic. It is an open space where you will be able to receive help to complete the application form, resolve doubts and receive direct guidance from the program team.

What is the Champions Program?

This program seeks to identify, support and recognize people who are already leading, or who want to take a step further, in building open science and sustainable research software communities.

Throughout 12 months the selected individuals will participate in:

• An initial cohort-based training of 13 2-hour workshops.
• The development of a research software project: creation of a new R package, submission of an existing R package to rOpenSci peer review, or as part of a team of reviewers of other R packages.
• 1-on-1 mentoring and cohort events.
• Leading community building activities, outreach and communications, and exchange of experiences and strengthening of regional networks.

The program also offers a stipend to recognize the time and work of participants for those who complete the program and a certificate of participation.

Who is it for?

Champions

Potential Champions are people who:

• Are from or live and work in Latin America;
• Use R and have previous experience in software development, reproducible analysis or open science;
• Are interested in developing an open research software project and strengthening their local or regional communities;
• Are able to communicate effectively in Spanish (oral and written);
• Accept our Code of Conduct;
• Commit to participating in the program activities organized for Champions during these 12 months.

Mentors

Potential Mentors are people who:

• Have experience developing, maintaining or contributing to R packages (ideally within the rOpenSci ecosystem, although this is not required);
• Have participated in software peer review processes;
• Have relevant mentoring experience;
• Have an interest in mentoring others, sharing experiences and actively contributing to the growth of the program;
• Are able to communicate effectively in Spanish (oral and written);
• Accept our Code of Conduct;
• Commit to participating in the program activities organized for mentors during these 12 months.

Why participate?

Being part of the rOpenSci Champions Program means:

• Having the opportunity to develop or strengthen a research software project with expert accompaniment.
• Being part of a regional network of people committed to open science and open software.
• Acquiring technical tools and community leadership skills.
• Receiving support to share your work and experiences with your community or institution and at regional and international events.
• Actively contributing to a more diverse, sustainable and inclusive open software ecosystem.

How to apply?

Applications will be made through online forms:

• Complete this form to apply to become a Champion
• Complete this form to apply to become a Mentor,

These forms are also available on the program website.

The forms must be completed in Spanish. On the program’s website you will also find more details regarding the requirements for Champions and Mentors, as well as answers to frequently asked questions.

Learn more about our work

On our Champion’s Program website you can find the detailed schedule of activities and the complete information of the program, as well as projects carried out by previous cohorts and the other types of activities carried out by the participants.

If you have any questions, we invite you to participate in the Community Call in January, join us in the application clinic in February, or contact our Community Manager.

We look forward to your applications with great enthusiasm! We want to continue building community in Latin America together with you.

Happy New Year from the team at Jumping Rivers!

As we move through the midpoint of the 2020s, it’s a good time to reflect on the changes that we have seen so far in this decade. In the world of data science nothing has dominated headlines quite like the rapid growth and uptake of generative artificial intelligence (GenAI).

Large language models (LLMs) such as ChatGPT, Claude and Gemini have incredible potential to streamline day-to-day tasks, whether that’s processing vast amounts of information, providing a human-like chat interface for customers or generating code. But they also come with notable risks if not harnessed responsibly.

Anyone that has interacted with these models is likely to have come across hallucination, where the model confidently presents false information as though it is factually correct. This can happen for a variety of reasons:

• LLMs often have no access to real-time information: how would a model that was trained last year know today’s date?
• The training data may be missing domain-specific information: can we really trust an off-the-shelf model to have a good understanding of pharmaceuticals and medicinal drugs?
• The model may be over-eager to come across as intelligent, so it decides to provide a confident output rather than a more nuanced, honest answer.

Often we need to give the model access to additional contextual information before we can make it “production-ready”. We can achieve this using a retrieval-augmented generation (RAG) workflow. In this blog post we will explore the steps involved and set up an example RAG workflow using free and open source packages in R.

What is RAG?

In a typical interaction with an LLM we have:

• A user prompt: the text that is submitted by the user.

• A response: the text that is returned by the LLM.

• (optional) A system prompt: additional instructions for how the LLM should respond (for example, "You respond in approximately 10 words or less").

In a RAG workflow we provide access to an external knowledge store which can include text-based documents and webpages. Additional contextual info is then retrieved from the knowledge store (hence “retrieval”) and added to the user prompt before it is sent. In doing so we can expect to receive a higher quality output.

How does it work?

Before going further, we must first introduce the concept of vectorisation.

Contrary to what you might believe, LLMs do not understand non-numerical text! They are mathematical models, meaning they can only ingest and output numerical vectors.

So how can a user interact with a model using plain English? The trick is that mappings exist which are able to convert between numerical vectors and text. These mappings are called “vector embeddings” and are used to convert the user prompt into a vector representation before it is passed to the LLM.

So, when setting up our RAG knowledge store, we have to store the information using a compatible vector representation. With this in mind, let’s introduce a typical RAG workflow:

1. Content: we decide which documents to include in the knowledge store.
2. Extraction: we extract the text from these documents in Markdown format.
3. Chunking: the Markdown content is split into contextual “chunks” (for example, each section or subsection of a document might become a chunk).
4. Vectorisation: the chunks are “vectorised” (i.e. we convert them into a numerical vector representation).
5. Index: we create an index for our knowledge store which will be used to retrieve relevant chunks of information.
6. Retrieval: we register the knowledge store with our model interface. Now, when a user submits a prompt, it will be combined with relevant chunks of information before it is ingested by the model.

At the retrieval step, a matching algorithm is typically used so that only highly relevant chunks are retrieved from the knowledge store. In this way, we are able to keep the size of the user prompts (and any incurred costs) to a minimum.

Setting up a RAG workflow in R

We will make use of two packages which are available to install via the Comprehensive R Archive Network (CRAN). Both are actively maintained by Posit (formerly RStudio) and are free to install and use.

{ragnar}

The {ragnar} package provides functions for extracting information from both text-based documents and webpages, and provides vector embeddings that are compatible with popular LLM providers including OpenAI and Google.

We will use {ragnar} to build our knowledge store.

{ellmer}

The {ellmer} package allows us to interact with a variety of LLM APIs from R. A complete list of supported model providers can be found in the package documentation.

Note that, while {ellmer} is free to install and use, you will still need to set up an API token with your preferred model provider before you can interact with any models. We will use the free Google Gemini tier for our example workflow. See the Gemini API documentation creating an API key, and the {ellmer} documentation for authenticating with your API key from R.

Example RAG workflow

We begin by loading the {ragnar} package.

library("ragnar")

The URL provided below links to the title page of the “Efficient R Programming” textbook, written by Robin Lovelace and our very own Colin Gillespie. We’re going to use a couple of chapters from the book to construct a RAG knowledge store.

url = "https://csgillespie.github.io/efficientR/"

Let’s use {ragnar} to read the contents of this page into a Markdown format.

md = read_as_markdown(url)

We could vectorise this information as it is, but first we should split it up into contextual chunks.

chunks = markdown_chunk(md)
chunks
#> # @document@origin: https://csgillespie.github.io/efficientR/
#> # A tibble: 2 × 4
#> start end context text
#> * <int> <int> <chr> <chr>
#> 1 1 1572 "" "# Efficient R programmin…
#> 2 597 2223 "# Welcome to Efficient R Programming" "## Authorsnn[Colin Gil…

The chunks are stored in a tibble format, with one row per chunk. The text column stores the chunk text (in the interests of saving space we have only included the start of each chunk in the printed output above).

The title page has been split into two chunks and we can see that there is significant overlap (chunk 1 spans characters 1 to 1572 and chunk 2 spans characters 597 to 2223). Overlapping chunks are perfectly normal and provides added context as to where each chunk sits relative to the other chunks.

Note that you can visually inspect the chunks by running ragnar_chunks_view(chunks).

It’s time to build our knowledge store with a vector embedding that is appropriate for Google Gemini models.

# Initialise a knowledge store with the Google Gemini embedding
store = ragnar_store_create(
 embed = embed_google_gemini()
)

# Insert the Markdown chunks
ragnar_store_insert(store, chunks)

The Markdown chunks are automatically converted into a vector representation at the insertion step. It is important to use the appropriate vector embedding when we create the store. A knowledge store created using an OpenAI embedding will not be compatible with Google Gemini models!

Before we can retrieve information from our store, we must create a store index.

ragnar_store_build_index(store)

We can now test the retrieval capabilities of our knowledge store using the ragnar_retreive() function. For example, to retrieve any chunks relevant to the text Who are the authors of “Efficient R Programming”? we can run:

relevant_knowledge = ragnar_retrieve(
 store,
 text = "Who are the authors of "Efficient R Programming"?"
)
relevant_knowledge
#> # A tibble: 1 × 9
#> origin doc_id chunk_id start end cosine_distance bm25 context text
#> <chr> <int> <list> <int> <int> <list> <lis> <chr> <chr>
#> 1 https://csgi… 1 <int> 1 2223 <dbl [2]> <dbl> "" "# E…

Note that the operators in "Efficient R Programming" have been used to print raw double quotes in the character string.

Without going into too much detail, the cosine_distance and bm25 columns in the returned tibble provide information relating to the matching algorithm used to identify the chunks. The other columns relate to the location and content of the chunks.

From the output tibble we see that the full content of the title page (characters 1 to 2223) has been returned. This is because the original two chunks both contained information about the authors.

Let’s add a more technical chapter from the textbook to the knowledge store. The URL provided below links to Chapter 7 (“Efficient Optimisation”). Let’s add this to the knowledge store and rebuild the index.

url = "https://csgillespie.github.io/efficientR/performance.html"

# Extract Markdown content and split into chunks
chunks = url |>
 read_as_markdown() |>
 markdown_chunk()

# Add the chunks to the knowledge store
ragnar_store_insert(store, chunks)

# Rebuild the store index
ragnar_store_build_index(store)

Now that our knowledge store includes content from both the title page and Chapter 7, let’s ask something more technical, like What are some good practices for parallel computing in R?.

relevant_knowledge = ragnar_retrieve(
 store,
 text = "What are some good practices for parallel computing in R?"
)
relevant_knowledge
#> # A tibble: 4 × 9
#> origin doc_id chunk_id start end cosine_distance bm25 context text
#> <chr> <int> <list> <int> <int> <list> <lis> <chr> <chr>
#> 1 https://csgi… 1 <int> 1 2223 <dbl [2]> <dbl> "" "# E…
#> 2 https://csgi… 2 <int> 1 1536 <dbl [1]> <dbl> "" "# 7…
#> 3 https://csgi… 2 <int> 22541 23995 <dbl [1]> <dbl> "# 7 E… "## …
#> 4 https://csgi… 2 <int> 23996 26449 <dbl [2]> <dbl> "# 7 E… "The…

Four chunks have been returned:

• One chunk from the title page of the textbook.

• One chunk from the start of Chapter 7.

• Two chunks from Section 7.5 (“Parallel Computing”).

It makes sense that we have chunks from Section 7.5, which appears to be highly relevant to the question. By including the title page and the start of Chapter 7, the LLM will also have access to useful metadata in case the user wants to find out where the model is getting its information from.

Now that we have built and tested our retrieval tool, it’s time to connect it up to a Gemini interface using {ellmer}. The code below will create a chat object allowing us to send user prompts to Gemini.

chat = ellmer::chat_google_gemini(
 system_prompt = "You answer in approximately 10 words or less."
)

A system prompt has been included here to ensure a succinct response from the model API.

We can register this chat interface with our retrieval tool.

ragnar_register_tool_retrieve(chat, store)

To check if our RAG workflow has been set up correctly, let’s chat with the model.

chat$chat("What are some good practices for parallel computing in R?")
#> Use the `parallel` package, ensure you stop clusters with `stopCluster()` (or
#> `on.exit()`), and utilize `parLapply()`, `parApply()`, or `parSapply()`.

The output looks plausible. Just to make sure, let’s check where the model found out this information.

chat$chat("Where did you get that answer from?")
#> I retrieved the information from "Efficient R programming" by Colin Gillespie
#> and Robin Lovelace.

Success! The LLM has identified the name of the textbook and if we wanted to we could even ask about the specific chapter. A user interacting with our model interface could now search online for this textbook to fact-check the responses.

In the example workflow above, we manually selected a couple of chapters from the textbook to include in our knowledge store. It’s worth noting that you can also use the ragnar_find_links(url) function to retrieve a list of links from a given webpage.

Doing so for the title page will provide the links to all chapters.

ragnar_find_links("https://csgillespie.github.io/efficientR/")
#> [1] "https://csgillespie.github.io/efficientR/"
#> [2] "https://csgillespie.github.io/efficientR/building-the-book-from-source.html"
#> [3] "https://csgillespie.github.io/efficientR/collaboration.html"
#> [4] "https://csgillespie.github.io/efficientR/data-carpentry.html"
#> [5] "https://csgillespie.github.io/efficientR/hardware.html"
#> [6] "https://csgillespie.github.io/efficientR/index.html"
#> [7] "https://csgillespie.github.io/efficientR/input-output.html"
#> [8] "https://csgillespie.github.io/efficientR/introduction.html"
#> [9] "https://csgillespie.github.io/efficientR/learning.html"
#> [10] "https://csgillespie.github.io/efficientR/performance.html"
#> [11] "https://csgillespie.github.io/efficientR/preface.html"
#> [12] "https://csgillespie.github.io/efficientR/programming.html"
#> [13] "https://csgillespie.github.io/efficientR/references.html"
#> [14] "https://csgillespie.github.io/efficientR/set-up.html"
#> [15] "https://csgillespie.github.io/efficientR/workflow.html"

You could then iterate through these links, extracting the contents from each webpage and inserting these into your RAG knowledge store. Just note, however, that including additional information in your store will likely increase the amount of text being sent to the model, which could raise costs. You should therefore think about what information is actually relevant for your LLM application.

Summary

In summary, we have introduced the concept of retrieval-augmented generation for LLM-powered workflows and built an example workflow in R using open source packages.

Before finishing, we are excited to announce that our new course “LLM-Driven Applications with R & Python” has just been added to our training portfolio. You can search for it here.

If you’re interested in practical AI-driven workflows, we would love to see you at our upcoming AI In Production 2026 conference which is running from 4-5 June in Newcastle-Upon-Tyne. If you would like to present a talk or workshop, please submit your abstracts before the deadline on 23 January.

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

In this post, I present the novelties of python package rtopy; a package allowing (whose ultimate objective is to) translate R to Python without much hassle. The intro is still available in available in https://thierrymoudiki.github.io/blog/2024/03/04/python/r/rtopyintro.

The novelties mainly concern the RBridge class and the call_r function. The RBridge class is more about persistency, while the call_r function is more about ease of use.

See for yourself in the following – hopefully comprehensive – examples (classification, regression, time series, hypothesis testing).

contents

1. Installation
2. RBridge class
3. call_r function
4. Advanced RBridge Usage Examples

 %load_ext rpy2.ipython

The rpy2.ipython extension is already loaded. To reload it, use:
  %reload_ext rpy2.ipython

%%R

install.packages("pak")
pak::pak(c("e1071", "forecast", "randomForest"))

library(jsonlite)

!pip install rtopy

"""
Advanced RBridge Usage Examples
================================

Demonstrates using R packages, statistical modeling, and data processing
through the Python-R bridge `rtopy`.
"""

import numpy as np
import pandas as pd
from rtopy import RBridge, call_r


# ============================================================================
# Example 1: Support Vector Machine with e1071
# ============================================================================
print("=" * 70)
print("Example 1: SVM Classification with e1071")
print("=" * 70)

# Generate training data
np.random.seed(42)
n_samples = 100

# Class 0: centered at (-1, -1)
X0 = np.random.randn(n_samples // 2, 2) * 0.5 + np.array([-1, -1])
# Class 1: centered at (1, 1)
X1 = np.random.randn(n_samples // 2, 2) * 0.5 + np.array([1, 1])

X_train = np.vstack([X0, X1])
y_train = np.array([0] * (n_samples // 2) + [1] * (n_samples // 2))

# Create R code for SVM training and prediction
svm_code = '''
library(e1071)

train_svm <- function(X, y, kernel_type = "radial") {
    # Convert to data frame
    df <- data.frame(
        x1 = X[, 1],
        x2 = X[, 2],
        y = as.factor(y)
    )

    # Train SVM
    model <- e1071::svm(y ~ x1 + x2, data = df, kernel = kernel_type, cost = 1)

    # Make predictions on training data
    predictions <- predict(model, df)

    # Calculate accuracy
    accuracy <- mean(predictions == df$y)

    # Return results
    list(
        predictions = as.numeric(as.character(predictions)),
        accuracy = accuracy,
        n_support = model$tot.nSV
    )
}
'''

rb = RBridge(verbose=True)
result = rb.call(
    svm_code,
    "train_svm",
    return_type="dict",
    X=X_train,
    y=y_train,
    kernel_type="radial"
)

print(f"Training Accuracy: {result['accuracy']:.2%}")
print(f"Number of Support Vectors: {result['n_support']}")
print(f"Sample Predictions: {result['predictions'][:10]}")


# ============================================================================
# Example 2: Time Series Analysis with forecast package
# ============================================================================
print("n" + "=" * 70)
print("Example 2: Time Series Forecasting with forecast")
print("=" * 70)

# Generate time series data
time_series = np.sin(np.linspace(0, 4*np.pi, 50)) + np.random.randn(50) * 0.1

ts_code = '''
library(forecast)

forecast_ts <- function(x, h = 10) {
    # Convert to time series object
    ts_data <- ts(x, frequency = 12)

    # Fit ARIMA model
    fit <- auto.arima(ts_data, seasonal = FALSE)

    # Generate forecast
    fc <- forecast(fit, h = h)

    # Return results
    list(
        forecast_mean = as.numeric(fc$mean),
        forecast_lower = as.numeric(fc$lower[, 2]),  # 95% CI
        forecast_upper = as.numeric(fc$upper[, 2]),
        model_aic = fit$aic,
        model_order = paste0("ARIMA(",
                            paste(arimaorder(fit), collapse = ","),
                            ")")
    )
}
'''

result = rb.call(
    ts_code,
    "forecast_ts",
    return_type="dict",
    x=time_series.tolist(),
    h=10
)

print(f"Model: {result['model_order']}")
print(f"AIC: {result['model_aic']:.2f}")
print(f"5-step forecast: {np.array(result['forecast_mean'])[:5]}...")


# ============================================================================
# Example 3: Random Forest with randomForest package
# ============================================================================
print("n" + "=" * 70)
print("Example 3: Random Forest Regression")
print("=" * 70)

# Generate regression data
np.random.seed(123)
X = np.random.rand(200, 3) * 10
y = 2*X[:, 0] + 3*X[:, 1] - X[:, 2] + np.random.randn(200) * 2

rf_code = '''
library(randomForest)

train_rf <- function(X, y, ntree = 500) {
    # Create data frame
    df <- data.frame(
        x1 = X[, 1],
        x2 = X[, 2],
        x3 = X[, 3],
        y = y
    )

    # Train random forest
    rf_model <- randomForest(y ~ ., data = df, ntree = ntree, importance = TRUE)

    # Get predictions
    predictions <- predict(rf_model, df)

    # Calculate R-squared
    r_squared <- 1 - sum((y - predictions)^2) / sum((y - mean(y))^2)

    # Get feature importance
    importance_scores <- importance(rf_model)[, 1]  # %IncMSE

    list(
        r_squared = r_squared,
        mse = rf_model$mse[ntree],
        predictions = predictions,
        importance = importance_scores
    )
}
'''

result = rb.call(
    rf_code,
    "train_rf",
    return_type="dict",
    X=X,
    y=y.tolist(),
    ntree=500
)

print(f"R-squared: {result['r_squared']:.3f}")
print(f"MSE: {result['mse']:.3f}")
print(f"Feature Importance: {result['importance']}")


# ============================================================================
# Example 4: Statistical Tests with stats package
# ============================================================================
print("n" + "=" * 70)
print("Example 4: Statistical Hypothesis Testing")
print("=" * 70)

# Generate two samples
group1 = np.random.normal(5, 2, 50)
group2 = np.random.normal(6, 2, 50)

stats_code = '''
perform_tests <- function(group1, group2) {
    # T-test
    t_result <- t.test(group1, group2)

    # Wilcoxon test (non-parametric alternative)
    w_result <- wilcox.test(group1, group2)

    # Kolmogorov-Smirnov test
    ks_result <- ks.test(group1, group2)

    list(
        t_test = list(
            statistic = t_result$statistic,
            p_value = t_result$p.value,
            conf_int = t_result$conf.int
        ),
        wilcox_test = list(
            statistic = w_result$statistic,
            p_value = w_result$p.value
        ),
        ks_test = list(
            statistic = ks_result$statistic,
            p_value = ks_result$p.value
        ),
        summary_stats = list(
            group1_mean = mean(group1),
            group2_mean = mean(group2),
            group1_sd = sd(group1),
            group2_sd = sd(group2)
        )
    )
}
'''

result = rb.call(
    stats_code,
    "perform_tests",
    return_type="dict",
    group1=group1.tolist(),
    group2=group2.tolist()
)

print(f"Group 1 Mean: {result['summary_stats']['group1_mean']:.2f} ± {result['summary_stats']['group1_sd']:.2f}")
print(f"Group 2 Mean: {result['summary_stats']['group2_mean']:.2f} ± {result['summary_stats']['group2_sd']:.2f}")
print(f"nT-test p-value: {result['t_test']['p_value']:.4f}")
print(f"Wilcoxon p-value: {result['wilcox_test']['p_value']:.4f}")


# ============================================================================
# Example 5: Data Transformation with dplyr
# ============================================================================
print("n" + "=" * 70)
print("Example 5: Data Wrangling with dplyr")
print("=" * 70)

# Create sample dataset
data = pd.DataFrame({
    'id': range(1, 101),
    'group': np.random.choice(['A', 'B', 'C'], 100),
    'value': np.random.randn(100) * 10 + 50,
    'score': np.random.randint(1, 101, 100)
})

dplyr_code = '''
library(dplyr)

process_data <- function(df) {
    # Convert list columns to data frame
    data <- as.data.frame(df)

    # Perform dplyr operations
    result <- data %>%
        filter(score > 50) %>%
        group_by(group) %>%
        summarise(
            n = n(),
            mean_value = mean(value),
            median_score = median(score),
            sd_value = sd(value)
        ) %>%
        arrange(desc(mean_value))

    # Convert back to list format for JSON
    as.list(result)
}
'''

result = rb.call(
    dplyr_code,
    "process_data",
    return_type="pandas",
    df=data
)

print("nGrouped Summary Statistics:")
print(result)


# ============================================================================
# Example 6: Clustering with cluster package
# ============================================================================
print("n" + "=" * 70)
print("Example 6: K-means and Hierarchical Clustering")
print("=" * 70)

# Generate clustered data
np.random.seed(42)
cluster_data = np.vstack([
    np.random.randn(30, 2) * 0.5 + np.array([0, 0]),
    np.random.randn(30, 2) * 0.5 + np.array([3, 3]),
    np.random.randn(30, 2) * 0.5 + np.array([0, 3])
])

cluster_code = '''
library(cluster)

perform_clustering <- function(X, k = 3) {
    # Convert to matrix
    data_matrix <- as.matrix(X)

    # K-means clustering
    kmeans_result <- kmeans(data_matrix, centers = k, nstart = 25)

    # Hierarchical clustering
    dist_matrix <- dist(data_matrix)
    hc <- hclust(dist_matrix, method = "ward.D2")
    hc_clusters <- cutree(hc, k = k)

    # Silhouette analysis for k-means
    sil <- silhouette(kmeans_result$cluster, dist_matrix)
    avg_silhouette <- mean(sil[, 3])

    list(
        kmeans_clusters = kmeans_result$cluster,
        kmeans_centers = kmeans_result$centers,
        kmeans_withinss = kmeans_result$tot.withinss,
        hc_clusters = hc_clusters,
        silhouette_score = avg_silhouette
    )
}
'''

result = rb.call(
    cluster_code,
    "perform_clustering",
    return_type="dict",
    X=cluster_data,
    k=3
)

print(f"K-means Within-cluster SS: {result['kmeans_withinss']:.2f}")
print(f"Average Silhouette Score: {result['silhouette_score']:.3f}")
print(f"nCluster Centers:n{np.array(result['kmeans_centers'])}")
print(f"nCluster sizes: {np.bincount(result['kmeans_clusters'])}")


print("n" + "=" * 70)
print("All examples completed successfully!")
print("=" * 70)

======================================================================
Example 1: SVM Classification with e1071
======================================================================
Training Accuracy: 100.00%
Number of Support Vectors: 9
Sample Predictions: [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]

======================================================================
Example 2: Time Series Forecasting with forecast
======================================================================
Model: ARIMA(3,1,0)
AIC: -10.21
5-step forecast: [0.29557391 0.4948255  0.64553023 0.80823028 0.93656539]...

======================================================================
Example 3: Random Forest Regression
======================================================================
R-squared: 0.972
MSE: 11.996
Feature Importance: [62.57255479535195, 86.55470841243113, 21.4933655703039]

======================================================================
Example 4: Statistical Hypothesis Testing
======================================================================
Group 1 Mean: 5.33 ± 2.06
Group 2 Mean: 5.37 ± 2.28

T-test p-value: 0.9381
Wilcoxon p-value: 0.8876

======================================================================
Example 5: Data Wrangling with dplyr
======================================================================

Grouped Summary Statistics:
  group   n  mean_value  median_score   sd_value
0     C  23   49.711861            76  11.367167
1     A  14   49.219788            74   9.744709
2     B  23   47.459312            80  10.126835

======================================================================
Example 6: K-means and Hierarchical Clustering
======================================================================
K-means Within-cluster SS: 39.38
Average Silhouette Score: 0.713

Cluster Centers:
[[-0.03545142  3.12736567]
 [ 2.9470395   3.04927708]
 [-0.07207628 -0.0825784 ]]

Cluster sizes: [ 0 30 30 30]

======================================================================
All examples completed successfully!
======================================================================

import matplotlib.pyplot as plt
import seaborn as sns

# Set a style for better aesthetics
sns.set_style("whitegrid")

# Create a scatter plot of the clustered data
plt.figure(figsize=(10, 7))
sns.scatterplot(
    x=cluster_data[:, 0],
    y=cluster_data[:, 1],
    hue=result['kmeans_clusters'],
    palette='viridis',
    s=100, # size of points
    alpha=0.8, # transparency
    legend='full'
)

# Plot the cluster centers
centers = np.array(result['kmeans_centers'])
plt.scatter(
    centers[:, 0],
    centers[:, 1],
    marker='X',
    s=200, # size of centers
    color='red',
    edgecolors='black',
    label='Cluster Centers'
)

plt.title('K-means Clustering of Generated Data')
plt.xlabel('Feature 1')
plt.ylabel('Feature 2')
plt.legend()
plt.grid(True)
plt.show()

from rtopy import RBridge, call_r

# ============================================================================
# Optional: SVM Classification (High vs Low Price)
# ============================================================================
print("n" + "=" * 70)
print("Optional: SVM Classification on Boston")
print("=" * 70)

svm_boston_class_code = '''
library(MASS)
library(e1071)

train_boston_svm_class <- function(kernel_type = "radial", cost = 1) {

    data(Boston)

    # Binary target: expensive vs cheap housing
    Boston$high_medv <- as.factor(ifelse(Boston$medv >
                                         median(Boston$medv), 1, 0))

    model <- svm(
        high_medv ~ . - medv,
        data = Boston,
        kernel = kernel_type,
        cost = cost,
        scale = TRUE
    )

    preds <- predict(model, Boston)

    accuracy <- mean(preds == Boston$high_medv)

    list(
        accuracy = accuracy,
        n_support = model$tot.nSV,
        confusion = table(
            predicted = preds,
            actual = Boston$high_medv
        )
    )
}
'''

result = rb.call(
    svm_boston_class_code,
    "train_boston_svm_class",
    return_type="dict",
    kernel_type="radial",
    cost=1
)

print(f"Classification Accuracy: {result['accuracy']:.2%}")
print(f"Number of Support Vectors: {result['n_support']}")
print("Confusion Matrix:")
print(result["confusion"])


======================================================================
Optional: SVM Classification on Boston
======================================================================
Classification Accuracy: 90.51%
Number of Support Vectors: 209
Confusion Matrix:
[[237, 29], [19, 221]]

# ============================================================================
# Optional: SVM Classification (High vs Low Price)
# ============================================================================
print("n" + "=" * 70)
print("Optional: SVM Classification on Boston")
print("=" * 70)

svm_boston_class_code = '''
library(MASS)
library(e1071)

train_boston_svm_class <- function(kernel_type = "radial", cost = 1) {

    data(Boston)

    # Binary target: expensive vs cheap housing
    Boston$high_medv <- as.factor(ifelse(Boston$medv >
                                         median(Boston$medv), 1, 0))

    model <- svm(
        high_medv ~ . - medv,
        data = Boston,
        kernel = kernel_type,
        cost = cost,
        scale = TRUE
    )

    preds <- predict(model, Boston)

    accuracy <- mean(preds == Boston$high_medv)

    list(
        accuracy = accuracy,
        n_support = model$tot.nSV,
        confusion = table(
            predicted = preds,
            actual = Boston$high_medv
        )
    )
}
'''

result = rb.call(
    svm_boston_class_code,
    "train_boston_svm_class",
    return_type="dict",
    kernel_type="radial",
    cost=1
)

print(f"Classification Accuracy: {result['accuracy']:.2%}")
print(f"Number of Support Vectors: {result['n_support']}")
print("Confusion Matrix:")
print(result["confusion"])


======================================================================
Optional: SVM Classification on Boston
======================================================================
Classification Accuracy: 90.51%
Number of Support Vectors: 209
Confusion Matrix:
[[237, 29], [19, 221]]

So you have used the excellent exiftool to extract all of the GPS-related information from a directory of photos in JPG format and write to a CSV file:

exiftool '-*GPS*' -ext jpg -csv . > outfile.csv

You’ve used R/leaflet to plot coordinates (latitude and longitude) before, but what about that tag named GPSImgDirection? It would be nice to have some kind of marker which indicates the direction in which you were facing when the photo was taken.

For me, a Google search provided hints but not one single, obvious straightforward solution to this problem (the generative AI effect? time will tell…), so here’s what I’ve put together from several sources, in particular this StackOverflow post.

The key points are:

• use awesomeIcons() to create a directional icon which can be rotated
• add the icons to your map using addAwesomeMarkers()

Here’s some example code which uses the Font Awesome icon long-arrow-up. Since “up” (north) corresponds to zero degrees, applying a rotation corresponding to GPSImgDirection should result in the correct orientation for the marker. The GPS-related tags in this case come from an iPhone 13.

library(readr)
library(leaflet)
library(sp)

outfile <- read_csv("outfile.csv")

# create dataset
# ugly but it works
dataset <- outfile %>%
  mutate(GPSLatitude = str_replace(GPSLatitude, " deg", "d"),
         GPSLatitude = GPSLatitude %>%
         char2dms() %>%
           as.numeric(),
         GPSLongitude = str_replace(GPSLongitude, " deg", "d"),
         GPSLongitude = GPSLongitude %>%
           char2dms() %>%
           as.numeric(),
         GPSHPositioningError = str_replace(GPSHPositioningError, " m", ""),
         GPSHPositioningError = GPSHPositioningError %>%
           as.numeric()) %>%
  select(latitude = GPSLatitude,
         longitude = GPSLongitude,
         GPSTimeStamp,
         GPSImgDirection,
         GPSHPositioningError)

# create the marker icons
icons <- awesomeIcons(iconRotate = dataset$GPSImgDirection,
                      icon = "long-arrow-up",
                      library = "fa",
                      markerColor = "white",
                      squareMarker = TRUE)

# create map
# can filter on positioning error if desired
leaflet(data = dataset %>%
  addProviderTiles(provider = providers$CartoDB.Positron) %>%
  addAwesomeMarkers(icon = icons, label = ~GPSTimeStamp)

Here’s a screenshot of the resulting interactive map.

# Core libraries for rugby data acquisition
library(tidyverse)
library(rvest)
library(httr)
library(jsonlite)

# Example: pulling match data from an API
response <- GET("https://api.example.com/rugby/match/9876")
raw_json <- content(response, "text")
match_data <- fromJSON(raw_json)

# Scraping a match statistics table
page <- read_html("https://example-rugby-site.com/match/9876")

team_stats <- page %>%
  html_node("table.match-stats") %>%
  html_table()

team_stats

# Standardizing and validating team statistics
team_stats_clean <- team_stats %>%
  janitor::clean_names() %>%
  mutate(across(where(is.character), str_trim)) %>%
  mutate(
    possession = as.numeric(possession),
    territory = as.numeric(territory)
  )

# Basic validation check
stopifnot(all(team_stats_clean$possession 

# Creating phase identifiers from ruck events
events <- events %>%
  arrange(match_id, event_time) %>%
  mutate(
    phase_id = cumsum(event_type == "ruck")
  )

# Summarising phase-level performance
phase_summary <- events %>%
  group_by(match_id, team, phase_id) %>%
  summarise(
    duration = max(event_time) - min(event_time),
    carries = sum(event_type == "carry"),
    meters = sum(meters_gained, na.rm = TRUE),
    turnovers = sum(event_type == "turnover"),
    .groups = "drop"
  )

# Player-level performance profile
player_profile <- events %>%
  group_by(player_id, player_name, position) %>%
  summarise(
    minutes_played = max(event_time) / 60,
    tackles = sum(event_type == "tackle"),
    missed_tackles = sum(event_type == "missed_tackle"),
    carries = sum(event_type == "carry"),
    meters = sum(meters_gained, na.rm = TRUE),
    offloads = sum(event_type == "offload"),
    penalties_conceded = sum(event_type == "penalty_conceded"),
    .groups = "drop"
  ) %>%
  mutate(
    tackles_per_min = tackles / minutes_played,
    meters_per_carry = meters / carries
  )

# Defensive performance by field channel
defense_analysis <- events %>%
  filter(event_type %in% c("tackle", "missed_tackle")) %>%
  group_by(team, field_channel) %>%
  summarise(
    tackles = sum(event_type == "tackle"),
    misses = sum(event_type == "missed_tackle"),
    success_rate = tackles / (tackles + misses),
    .groups = "drop"
  )

# Kicking distance and efficiency
kicks <- events %>%
  filter(event_type == "kick") %>%
  mutate(kick_distance = end_x - start_x)

kicking_summary <- kicks %>%
  group_by(team, kick_type) %>%
  summarise(
    avg_distance = mean(kick_distance, na.rm = TRUE),
    kicks = n(),
    .groups = "drop"
  )

# Building a basic win probability model
wp_data <- matches %>%
  mutate(
    score_diff = team_score - opponent_score,
    time_remaining = 80 - minute
  )

wp_model <- glm(
  win ~ score_diff + time_remaining + territory,
  data = wp_data,
  family = binomial()
)

summary(wp_model)

# Creating a clean match summary table
summary_table <- team_stats_clean %>%
  select(team, possession, territory, tackles, line_breaks, penalties_conceded)

knitr::kable(summary_table)

The post Rugby Analytics with R: Complete Guide to Performance Analysis in Rugby Union and League appeared first on R Programming Books.

I remember the day that I started to use R programming. I had a basic interface to write and execute the code. After that experience, R Studio emerged as a powerful IDE for R programming for me. It provided a user-friendly interface, integrated tools, and features that enhance productivity and streamline the coding process and was a huge shift for me in my R programming journey.

In July 2022, R Studio was rebranded to Posit. Apparently, a new era was about to start because the world’s needs were evolving, and R had a stronger companion in the Python programming language.

R Studio Interface (Source: biocorecrg.github.io)

To satisfy the needs of both R and Python users, Posit introduced a new product called Positron. It is a data science-oriented IDE that supports both R and Python programming languages, in contrast to R Studio. Of course, this emerging tool has tempted some R Studio users who are also using VSCode since it offers some advantages over R Studio.

Positron Interface (Source: https://positron.posit.co/)

The main difference between Positron and R Studio is their multi-language support. Positron allows users to work with both R and Python in a single environment, making it easier for data scientists who use both languages. Additionally, Positron offers better integration with Jupyter Notebooks, which are widely used in the data science community.

AI-based assistants are also integrated into Positron, providing users with suggestions and code completions based on their coding patterns. This feature can significantly enhance productivity and reduce the learning curve for new users.

If you are playing with the data, it offers more flexibility and versatility compared to R Studio. You can examine not only the data frames on your enviroment, but also .csv and parquet without importing them.

Another advantage of Positron is to offer extensions that make the IDE more customizable and adaptable to different workflows. Users can install extensions to add new features, improve functionality, and tailor the environment to their specific needs.

Package versions and R versions crashes sometimes becomes annoying if you have encounter this during your R Studio experiment. But, with Positron, you can manage different R versions at the same time on the same machine without conflicts. This is particularly useful for users who work on multiple projects with different R version requirements.

Lastly, it is being improved continuously with frequent updates and new features being added regularly. This ensures that users have access to the latest tools and technologies in the data science field.

So the question is: Should we give up on using R Studio?

Actually, no. Because it is not going away, and it still provides some advantages over Positron.

R Studio still has strong properties that tempt users use it. You can use RMarkdown and Quarto to create dynamic documents, reports, and presentations that combine code, text, and visualizations. R Studio also has a robust ecosystem of packages and extensions that enhance its functionality and provide specialized tools for various data analysis tasks.

You can save and reload your workspace. Besides, you have several panels that help you to manage your files, plots, packages, and help documents easily. You can track your codes written in the past and bring them back easily without spending long time. And, you can import your datasets without typing code!

From a developer perspective, R Studio has specific tools that makes developing packages and app easier compared to Positron.

In conclusion, both Positron and R Studio have their own strengths and weaknesses. The choice between the two ultimately depends on the user’s specific needs and preferences. If you require multi-language support, better Jupyter integration, and AI-based assistance, Positron may be the better choice. However, if you prioritize RMarkdown, a robust package ecosystem, and workspace management, R Studio may be more suitable.

Introduction: Copyright Law and Whatnot

I recently watched a clip from a debate between Timothy Nguyen of Google Deepmind and Danish author Janne Teller. The debate, entitled Technology and Freedom, took place at Hay-on-Wye in summer 2025. On the subject of copyright, Nguyen says:

The reason AI is so powerful is because it’s scraped all this data on the internet and of course that has all these issues in terms of copyright law and whatnot. But that’s also democratised knowledge and so there are ways in which it’s been good and bad. But now we have this very intelligent system which has all this knowledge from books, but then maybe there are going to be some authors who aren’t going to be very happy. So there are always going to be winners and losers.

Teller replies:

This is an undermining of any intellectual property rights we have developed up to now. Anything you have written on a Facebook post which is public will be considered by this Metaverse as something they can use to develop their AI and you might say OK, that’s a new form of sharing. Anything you contribute, everybody owns it. But then that speaks to nationalising all technology platforms. You want to have everything everyone else has created. But then we want to have your work also and have control over it.

The clip cuts off here and I haven’t seen a video of the full debate, so I don’t know how Nguyen replied. But I think Teller made a good point. It’s not just that LLMs have been trained (illegally) on masses of copyrighted material. They have also been trained on data from the internet, which is a public good, and perhaps the people who unwittingly created all the training data should be entitled to some sort of compensation. Even the slopigarchs themselves acknowledge this. For example, in 2017, Elon Musk said that the pace of change is:

a massive social challenge. And I think ultimately we will have to have some kind of universal basic income (UBI). I don’t think we’re going to have a choice.

At the moment we are facing two possible outcomes. Either AI progress grinds to a halt and the bubble bursts, or AI breakthroughs continue to happen at a rapid pace, replacing human jobs, and everyone ends up becoming unemployed until they can find other jobs to do. Every previous improvement in technology, no matter how disruptive, eventually ended up with people finding other things to do, so the economy will keep going somehow. But before we reach that point we may find ourselves facing serious social unrest. As Teller suggests, perhaps it is the AI companies themselves who should pay for this. After all, they did steal everyone else’s work to train their models. But if, in some grim future in which companies like OpenAI become profitable, we eventally get compensation, how much compensation should we get?

It seems like this question has no answer. But actually there’s a simple heuristic for evaluating the relative contributions of the model and the data, which I want to explain in this post. Not only is this heuristic relevant for musing about the future of AI, but it’s also surprisingly useful in everyday data science, too.

The Cover-Hart Theorem

Consider a classification problem in which the input is a data point $x$ contained in some metric space (i.e. a set equipped with a notion of distance) $(X, d)$, and the output is a classification into one of $M$ classes. The classifier is evaluated by the percentage of data points which it classifies correctly (the accuracy). If $A$ is the accuracy then $R = 1-A$ is called the error rate.

The Bayes Rate $R^ast$ for the problem is defined to be the best possible accuracy which any classifier could have. Why isn’t $R^ast$ just 100%? That’s because the same point might appear in more than one class! See the example below.

Suppose a data set $mathcal{X}$ is given. It consists of some points $x_i in X$ and the corresponding classes $theta_i$. We want to use the data set $mathcal{X}$ to build a classifier.

The 1-Nearest Neighbour or 1-NN classifier is the classifier $C$ which simply assigns an unseen data point $x$ to the class of the closest point to $x$ in $mathcal{X}$ (for simplicity, let’s assume that $mathcal{X}$ doesn’t contain any duplicate points). That is, if $d(x, x_i) = min_{y in mathcal{X}}d(x, y)$ then $C(x) := theta_i$. Note that to define the 1-NN classifier, we need $X$ to be a metric space, or else there is no notion of the nearest neighbour.

The theorem which Cover and Hart proved in 1967 is that the error rate $R$ of the 1-NN classifier satisfies

asymptotically as the number of data points in $mathcal{X}$ goes to $infty$, and provided that the points in $mathcal{X}$ are an iid sample from some distribution.

In other words, if you are given a data set and asked to build a predictive model, just doing the most naive thing possible and looking up the closest point in your data set to the point you want to classify already gets you halfway to the lowest possible error.

Example

Here is an example which I used to use when teaching this topic in university courses.

Let’s consider a single predictor $x$. There are two classes labelled $0$ and $1$. The distribution of $x$ for class $1$ is $N(1, 1)$ and the distribution of $x$ for class $0$ is $N(-1, 1)$. Suppose the population is equally distributed among the two classes.

The best possible classifier would classify a point $x$ into whichever class has the higher density for that particular value of $x$. The purple area represents the proprtion of points which would be misclassified. Since 50% of the population is in each class, the purple area is equal to

bayes_rate <- 1-pnorm(1)
# 0.1586553

Now suppose we are supplied with a training dataset consisting of 50 points from each class

set.seed(100)

N <- 100
df_train <- data.frame(x = c(rnorm(N/2, 1, 1), rnorm(N/2, -1, 1)), y = rep(c(1, 0), each=N/2))

The following function classifies a point using the nearest neighbour with the metric being $d(x, y) = lvert x-y rvert$.

classify_point <- function(x, df){
  df$y[which.min(abs(x-df$x))]
}

To see whether the Cover-Hart Theorem works in this example, let’s create a test data set of 1000 new points.

M <- 1000
df_test <- data.frame(x = c(rnorm(M, 1, 1), rnorm(M, -1, 1)), y = c(rep(1, M), rep(0, M)))

The error rate of the 1-NN classifier on this data set can be calculated as follows

pred <- sapply(as.list(df_test$x), function(x) classify_point(x, df_train))
1 - sum(pred == df_test$y)/length(pred)
# 0.216

As expected, $0.216$ lies between the Bayes rate and twice the Bayes rate.

Of course, many other classifiers will perform better. For example, logistic regression already almost achieves perfect accuracy.

model <- glm(y~x, data=df_train, family="binomial")
pred_logistic <- round(predict(model, df_test, type="response"))
1 - sum(pred_logistic == df_test$y)/length(pred)
# 0.16

If you run the whole script again with the same seed but with N=10000 points in the training data, you will even find that logistic regression gets an error rate which is lower than the Bayes rate! This happens because the training and test sets are finite samples from the actual data distribution, so there is some sampling error.

Practical Use

There are two ways to use this in practice. First, suppose that you are presented with a data set and build a quick and dirty classifier using 1-NN and achieve an accuracy of 80%. Then the error rate $R$ of the 1-NN classifier is 20% and the Cover-Hart Theorem tells you that the Bayes rate $R^ast ge R/2$, so the Bayes rate cannot be less than 10%, which means that you can’t expect to achieve an accuracy of better than 90% using some other algorithm. This might be a helpful guide to how much time you should spend trying to build a better classifier. In practice, the quick and dirty classifier you build will be something other than 1-NN¹, and it usually has better performance than 1-NN, so this can actually be a useful way to estimate the Bayes rate on a new data set.

Secondly, suppose that you are presented with a classification algorithm with an accuracy of 95%. Then you can estimate that the Bayes rate $R^ast$ is at most 5%, because $R^ast$ is the lowest possible error rate among all classifiers. This means that the error rate of a 1-NN classifier $R$ cannot be more than 10%. But that means that a 1-NN classifier would have given you 90% accuracy. Since the 1-NN classifier is just another name for “look at the data”, you could already achieve 90% accuracy by looking at the data alone without building your fancy model. In other words, the data is doing $90/95 = 94.7%$ of the work!²

Problems with the Cover-Hart Theorem

In practice, Cover-Hart should be used only as a heuristic and not as something which is expected to hold in all cases. This is because it makes very strong assumptions about the data.

• It’s only true if you have infinitely many data points, so it will only be approximately true in any real-life situation. How close the Cover-Hart Theorem is to being true in any real-life situation might also depend very strongly on the metric being used.
• More seriously, the data points need to be independent and identically distributed (iid). This is never true, despite the fact that textbooks and courses seem to suggest otherwise. In fact, I think it’s rare for training and test sets even to come from the same distribution.

For example, consider image classification. Cover-Hart suggests that you can classify any image correctly if you find the closest image, perhaps in the sense of Euclidean distance, in some sufficiently large reference data set. But clearly the reference data set would have to be massive, and the cost of searching for the closest image would probably be extremely high.

What does Cover-Hart say about AI?

The Cover-Hart Theorem, then, doesn’t suggest a sensible way to build an AI model. For example, suppose you want to generate the next word, given a string of text. A 1-NN classifier would be supplied with a corpus of data. It would need to search through this data and find the piece of text that was the closest match to the given string, and then extract the next word from that piece of text. For some kinds of text, like The capital of France is, this might work well, but clearly it’s not going to be a good approach in general.

This isn’t how Large Language Models work at all, so how is the Cover-Hart Theorem relevant to LLMs? Well, I think it could be used as a heuristic for measuring the relative contribution of the model and the data. For example, let’s suppose an LLM has an accuracy performance of $A$ percent on some benchmark. Then, as explained above, a 1-NN classifier could be expected to achieve an error rate of $2(1-A)$ and so you could estimate that the data by itself is contributing roughly

of the overall performance. This could be taken as a measure of how much the data is “worth” versus the model.

For example, if an AI company achieves 80% on some benchmark, then the people who contributed the data in some sense deserve $(2(0.8)-1)/0.8 = 75%$ of the credit.

Who was right?

So was Teller correct? Do the people who generated the data deserve most of the profits (if there are any) from AI? Well, that depends on what you mean by “AI”.

In the case of LLMs, assuming that they really are able to replace people in the workplace, I think the Cover-Hart Theorem could provide a first step for deciding how to regulate or tax. But the term “AI” encompasses a lot of different models, and some of those models don’t use training data at all. For example, AlphaZero reached grandmaster-level performace in chess and superhuman performance in go by playing against itself. And this isn’t a new idea; in the 1990s TD-Gammon was already able to outperform humans in backgammon by taking a similar approach. Personally I find these kind of algorithms even more impressive than LLMs, but that’s just my opinion.

One more thing. Suppose we did find ourselves in a world in which a government was choosing to tax AI companies based on the above formula. Then we could reach a bizarre scenario in which, in order to avoid tax, the AI companies would be competing to make the ratio $(2A-1)/A$ as small as possible. This would mean that, instead of boasting about the accuracy of their models on self-chosen benchmarks, these firms would find themselves in a paradoxical race to claim that their accuracy was as low as possible.

I think that would be hilarious.

[1] By the way, the 1-NN classifier is one of the very few classifiers which outputs just a class without any notion of the strength of class membership, so you can’t define an AUC for it. This is one of the classifiers which suffers from the so-called class imbalance problem, which they ask about in every data science interview. In pratice, class imbalance is never really a problem because nobody compares classifiers by using accuracy alone.

[2] Of course, this might not be the full story. For one thing, you will probably be interested in other metrics besides accuracy. For another, your algorithm might have other advantages over the 1-NN classifier, such as coming up with predictions more quickly.

Read it in: Español.

rOpenSci’s activities and spaces are supported by a Code of Conduct that applies to all people participating in the rOpenSci community, including rOpenSci staff and leadership. It applies to all modes of interaction including GitHub project repositories, the rOpenSci discussion forum, Slack, Community Calls, Co-working and social sessions, training and mentoring sessions, and in person at rOpenSci-hosted events, including affiliated social gatherings. Our Code of Conduct is developed and enforced by a committee including rOpenSci staff and an independent community member.

Here we report our annual review of rOpenSci’s Code of Conduct, reporting process, and internal guidelines for handling reports and enforcement.

Updates

1. No changes of Committee members for 2026. The committee members are Natalia Morandeira (independent member), Mark Padgham (rOpenSci Software Research Scientist) and Yanina Bellini Saibene (rOpenSci Community Manager). We are responsible for receiving, investigating, adjudicating, enforcing, and responding to all reports of potential violations of our Code.

2. No changes have been made to the text of the Code in English, Spanish or Portuguese.

3. Because the text has remained unchanged, the version number is kept at 2.5, dated January 30, 2024.

Committee Members Activities

An important aspect of our work as Committee members is to make sure that the processes we follow are transparent, consistent, and fair. To support this work we’ve developed a set of templates that guide us through different stages of incident response and reporting. We shared openly these templates and explain them in detail in the blog post “rOpenSci Code of Conduct Committee Templates”.

Yanina Bellini Saibene attended the BoF “Community safety and wellbeing: Can we do more than “good enough”? at the Open Science CZI meeting to discuss Code of Conduct best practices with other organizations in the Open Science movement.

Contact Information

We welcome your feedback by email to conduct at ropensci.org, and we thank you for continuing to work with us to ensure that rOpenSci remains a safe, enjoyable, friendly and enriching environment for everyone who participates.