Join our workshop on Introduction to All Things About bpvars, the R package for Forecasting with Bayesian Panel Vector Autoregressions,  which is a part of our workshops for Ukraine series! 

Here’s some more info: 

Title: All Things About bpvars, the R package for Forecasting with Bayesian Panel Vector Autoregressions

Date: Thursday, September 10th, 18:00 – 20:00 CEST (Rome, Berlin, Paris timezone) 

Speaker: Tomasz Woźniak is a Bayesian econometrician developing new econometric methods for applied macroeconomic research. He has been a specialised R user for eighteen years, and has recently joined The R Journal as an Associate Editor. He is the author of several R packages, available at https://bsvars.org/, that combine blazingly fast algorithms written in C++ with the convenience of data analysis in R. He works as a senior lecturer at the University of Melbourne, where he has an extensive research, teaching, and engagement portfolio. Tomasz supports Ukraine. Slava Ukrainie!

Description: This session is a unique opportunity to learn to use the bpvars package for forecasting with Bayesian panel vector autoregressions in a wide range of applications spanning:

* data preparation,

* specification of panel vector autoregressions,

* estimation of the models,

* estimation with missing observations,

* forecasting labour market outcomes globally,

* plotting and reporting predictions of different qualities, including point and density forecasts, marginal, conditional, or restricted forecasts,

* performing recursive expanding window forecasting,

* generating forecasting performance reports.

A sequence of hands-on exercises that fully prepares attendees to work with the package supports all of this. You are welcome to browse existing resources, including:

* bpvars website: https://bsvars.org/bpvars/

* bpvars CRAN profile: https://cran.r-project.org/package=bpvars

* package vignette: https://doi.org/10.48550/arXiv.2606.14143

* other materials: https://bsvars.org/bpvars/#resources

See you soon!

Preparation:

This session assumes the attendees know the basics of time series analysis. Install and load the bpvars package, then check the documentation by running the code below. Run the example you will find there in R. If it works, you’re ready to go. Thanks!

“`

install.packages(“bpvars”)

library(bpvars)

?bpvars

Minimal registration fee: 20 euro (or 20 USD or 800 UAH)

Please note that the registration confirmation is sent 1 day before the workshop to all registered participants rather than immediately after registration

How can I register?

• Go to https://bit.ly/3wvwMA6 or https://bit.ly/4aD5LMC  or  https://bit.ly/3PFxtNA and donate at least 20 euro. Feel free to donate more if you can, all proceeds go directly to support Ukraine.

• Save your donation receipt (after the donation is processed, there is an option to enter your email address on the website to which the donation receipt is sent)

• Fill in the registration form, attaching a screenshot of a donation receipt (please attach the screenshot of the donation receipt that was emailed to you rather than the page you see after donation).

If you are not personally interested in attending, you can also contribute by sponsoring a participation of a student, who will then be able to participate for free. If you choose to sponsor a student, all proceeds will also go directly to organisations working in Ukraine. You can either sponsor a particular student or you can leave it up to us so that we can allocate the sponsored place to students who have signed up for the waiting list.

How can I sponsor a student?

• Go to https://bit.ly/3wvwMA6 or https://bit.ly/4aD5LMC  or https://bit.ly/3PFxtNA and donate at least 20 euro (or 17 GBP or 20 USD or 800 UAH). Feel free to donate more if you can, all proceeds go to support Ukraine!

• Save your donation receipt (after the donation is processed, there is an option to enter your email address on the website to which the donation receipt is sent)

• Fill in the sponsorship form, attaching the screenshot of the donation receipt (please attach the screenshot of the donation receipt that was emailed to you rather than the page you see after the donation). You can indicate whether you want to sponsor a particular student or we can allocate this spot ourselves to the students from the waiting list. You can also indicate whether you prefer us to prioritize students from developing countries when assigning place(s) that you sponsored.

If you are a university student and cannot afford the registration fee, you can also sign up for the waiting list here. (Note that you are not guaranteed to participate by signing up for the waiting list).

You can also find more information about this workshop series,  a schedule of our future workshops as well as a list of our past workshops which you can get the recordings & materials here.

Looking forward to seeing you during the workshop!

I’m a big fan of Github Copilots ghost text and next edit suggestions. As a data analyst they give you much more control over the pace and direction of analysis code than a full AI agent. They also let you practice the core skill an agent needs from you: writing a clear specification.

I’ll walk through it in R, but the same ideas apply in any language. I’m assuming you’re on VScode with the Copilot extension installed. This should work with the free or paid plans (though I haven’t checked free myself).

library(dplyr)
library(ggplot2)
dat <- data.frame(
    time = rep(1:10, 3),
    species = rep(c("A", "B", "C"), each = 10),
    value = c(rnorm(10, mean = 5), rnorm(10, mean = 10), rnorm(10, mean = 15))
)

datA <- dat |> filter(species == "A")
m1 <- lm(value ~ time, data = datA)
coef(m1)["time"]

ggplot(datA, aes(x = time, y = value)) +
    geom_point() +
    geom_smooth(method = "lm", se = FALSE) +
    labs(
        title = "Species A Time Series",
        x = "Time",
        y = "Value"
    )

fit_fun <- function(species_name) {
    datA <- dat |> filter(species == species_name)
    m1 <- lm(value ~ time, data = datA)
    coef(m1)["time"]

    ggplot(datA, aes(x = time, y = value)) +
        geom_point() +
        geom_smooth(method = "lm", se = FALSE) +
        labs(
            title = paste("Species", species_name, "Time Series"),
            x = "Time",
            y = "Value"
        )
}

fit_fun("A")

# Simulate a new dataset of abundance at x-y coordinates

# plot a 2D map

# fit a model with interaction between x and y

Four hundred twelve new packages were submitted to CRAN in June. Here are my Top 40 picks in nineteen categories: Bioarchaeology, Biology, Climate Studies, Computational Methods, Ecology, Epidemiology, Finance, Functional Data Analysis, Machine Learning, Medical Statistics, Networks, Pharmacokinetics, Probability, Programming, Psychometrics, Risk Analysis, Statistics, Time Series, and Utilities.

This post describes the R version of Python’s GPopt (https://docs.techtonique.net/GPopt/GPopt.html), a package for ‘Bayesian’ optimization of black-box functions and machine learning hyperparameter tuning, using Gaussian Process Regression and other conformalized surrogates. The package is available on GitHub, and through the R universe.

Keep in mind that this package is for Machine Learning hyperparameter tuning: the global minimum won’t always be found, but this isn’t an issue, since it means you aren’t overfitting the training set.

It’s ported the same way as nnetsauce for R was: with uv to create an isolated Python virtual environment containing the Python GPopt package, and reticulate to call into it from R. Every function in this R package is a thin wrapper that returns the underlying Python object; the general rule is: object accesses with .’s in Python are replaced by $’s in R.

See this post for the technique: Finally figured out a way to port python packages to R using uv and reticulate.

Install

1. Create a Python virtual environment with uv

# pip install uv # if necessary
uv venv venv
source venv/bin/activate      # on Windows: venv\Scripts\activate
uv pip install pip GPopt

Keep track of where venv/ lives – you’ll pass its path as venv_path to every function in this package.

2. Install the R package

install.packages("remotes")
remotes::install_github("Techtonique/GPopt_r") 

reticulate will be installed automatically as a dependency.

Examples

Minimizing the Branin function

This is a standard test function for optimization algorithms. GPOpt is more suitable for expensive black-box functions, but this is a good example to illustrate the usage of the package.

library(GPopt)

branin <- function(x) {
  x1 <- x[1]; x2 <- x[2]
  term1 <- (x2 - (5.1 * x1^2) / (4 * pi^2) + (5 * x1) / pi - 6)^2
  term2 <- 10 * (1 - 1 / (8 * pi)) * cos(x1)
  term1 + term2 + 10
}

opt <- GPOpt(
  lower_bound = c(-5, 0),
  upper_bound = c(10, 15),
  objective_func = branin,
  n_init = 10,
  n_iter = 40,
  venv_path = "./venv"
)

opt$optimize(verbose = 1L)
print(opt$x_min)  # best parameters
print(opt$y_min)  # best objective value

Tuning a scikit-learn model’s hyperparameters

library(GPopt)

sklearn <- get_sklearn(venv_path = "./venv")
RandomForestClassifier <- sklearn$ensemble$RandomForestClassifier

X <- as.matrix(iris[, 1:4])
y <- as.integer(iris$Species) - 1L

mlopt <- MLOptimizer(scoring = "accuracy", cv = 5, venv_path = "./venv")

param_config <- list(
  n_estimators = list(bounds = c(10, 300), dtype = "int"),
  max_depth    = list(bounds = c(1, 20),   dtype = "int")
)

mlopt$optimize(
  X_train = X, y_train = y,
  estimator_class = RandomForestClassifier(),
  param_config = param_config,
  verbose = 1L
)

print(mlopt$get_best_parameters())
print(mlopt$get_best_score())

Bayesian optimization with early stopping

library(GPopt)

opt <- BOstopping(
  f = branin,
  bounds = rbind(c(-5, 10), c(0, 15)),
  venv_path = "./venv"
)
result <- opt$optimize(n_iter = 100L)

Using a custom (conformalized) surrogate model

library(GPopt)

sklearn <- get_sklearn(venv_path = "./venv")
ns <- get_nnetsauce(venv_path = "./venv")

opt <- GPOpt(
  lower_bound = c(-5, 0),
  upper_bound = c(10, 15),
  objective_func = branin,
  acquisition="ucb",
  method="splitconformal",
  surrogate_obj = ns$PredictionInterval(sklearn$ensemble$RandomForestRegressor()),
  venv_path = "./venv"
)
opt$optimize(verbose = 1L)

Biological phenomena are often studied by examining how a numerical variable, usually called the response (e.g., the weight or height of an organism), is affected by another variable, usually called the predictor (e.g., time or fertiliser application rate). These ‘predictor-response’ relationships are commonly described using models, which express the response as a mathematical function of the predictor. In general, a model can be written as:

where is the response, is the predictor, and is a collection of parameters, usually denoted by symbols, such as the letters of the Roman or Greek alphabet. The element is the function, which determines the shape of the relationship when plotted on an x–y graph.

Because mathematical modelling plays such a central role in biology and many other scientific disciplines, biologists need to be familiar with the most important mathematical functions. More importantly, they need to be able to “read” these functions and use their parameters to describe, interpret, and quantify biological processes. With this aim in mind, I have compiled a collection of the mathematical functions most commonly encountered in biology, explaining the meaning of their parameters, with particular emphasis on their biological interpretation rather than their mathematical properties.

Curve shapes

Functions are often classified according to the shape they exhibit when plotted on an x–y graph. This approach is adopted, for example, in Ratkowsky (1990), and the classification presented below is largely based on that work.

• Polynomials
  1. Straight line function
  2. Quadratic polynomial function
• Concave/Convex curves (no inflection)
  1. Exponential function
  2. Asymptotic function / Negative exponential function
  3. Power function
  4. Logarithmic function
  5. Rectangular hyperbola
• Sigmoidal curves
  1. Logistic function
  2. Gompertz function
  3. Modified Gompertz function
  4. Log-logistic function
  5. Weibull (type-1) function
  6. Weibull (type-2) function
• Curves with maxima/minima
  1. Peaked sigmoidal function
  2. Bragg function
  3. Lorentz function
  4. Beta function

I have chosen to include a relatively large number of functions, so this post is necessarily rather long. To make navigation easier, you can first inspect Figure 1 to identify the type of curve of interest and then use the links above to jump directly to the corresponding section.

Throughout this post, denotes the predictor, denotes the response, and the letters of the Roman alphabet (, , , …) denote the model parameters, which determine the shape of the function.

D(expr = expression(a + b*X + c*X^2), name = "X")

b + c * (2 * X)

Concave/Convex curves (no inflection points)

Exponential function

The most common parameterisation of an exponential function is:

Depending on the sign of , Eq. 3 describes a concave-up monotonically increasing shape (exponential growth; ) or a concave-up monotonically decreasing shape (exponential decay; ).

The parameter represents the value of when , whereas the meaning of can be understood by considering the first derivative:

D(expression(a * exp(k * X)), "X")

a * (exp(k * X) * k)

In algebraic terms, it is:

and thus:

Therefore, represents the relative rate of change (increase/decrease), throughout the domain, which is often known in growth analysis as the Relative Growth Rate (RGR).

Other equivalent parameterisations are:

and:

Equations 3-5 are equivalent, as can be shown by setting and :

Another, slightly different, parameterisation is commonly used in bioassay studies, mainly to describe exponential decay processes:

where corresponds to in the previous models and .

For all the exponential decay equations presented above, as . A lower asymptote, , can be incorporated into Eq. 6 to account for situations in which the response does not approach zero as the predictor tends to infinity:

Exponential functions are often used to describe the growth of populations under non-limiting environmental conditions or the degradation of xenobiotics in the environment (first-order degradation kinetics). In both cases, represents time and is therefore restricted to non-negative values, while (or, equivalently, ) must be different from zero. Check the post at this link for an example of fitting.

Asymptotic function

This function appears in several different parameterisations and is also known as the monomolecular growth function, the Mitscherlich law, or the von Bertalanffy law. Owing to its biological interpretation, the most widely used parameterisation is:

It describes a monotonically increasing, concave-down curve approaching a horizontal asymptote as tends to infinity. The parameter represents the maximum attainable value of (the plateau), while is the value of when (the initial value). The parameter is proportional to the rate at which the response approaches the plateau. Indeed, the first derivative is

D(expression(a - (a - b) * exp (- m * X)), "X")

(a - b) * (exp(-m * X) * m)

Considering Eq.8, we can write:

and then:

The previous expression shows that the relative rate of change in the response (the RGR in growth analysis) is not constant, as it is for the exponential function. Instead, it depends on the attained value of . In particular, the RGR is greatest at the beginning of the process, when is smallest, and gradually approaches zero as approaches the plateau .

Another closely related parameterisation is often encountered, in which is replaced by , by , and by . Some simple algebraic manipulations are also introduced to make the model consistent with the parameterisation commonly adopted for biological assay models:

In some applications, is reparameterised through its logarithm to facilitate model fitting:

For all these equations, setting (or equivalently ) yields a curve passing through the origin. In this case, the most common parameterisation is:

which is known as the negative exponential function.

The asymptotic function is widely used in biology, not only for growth analysis but also, for example, in weed competition studies. The negative exponential function has also been used to model the absorbed Photosynthetically Active Radiation () as a function of leaf area index (). In this context, represents the incident PAR (), while is the light extinction coefficient. Check the post at this link for an example of fitting.

Power function

The power function is also known as Freundlich function or allometric function. The most common parameterisation is:

Its shape is highly flexible and is determined by the value of the parameter (see Figure 1). When , the response increases with and the curve is concave down. When , decreases as increases and the curve is concave up. Finally, when , increases with and the curve is concave up. These three cases are illustrated in Figure 1 using different colours. The function is defined only for and has no horizontal asymptote as .

The biological interpretation of the parameters and is not straightforward. Both influence the slope of the curve, as can be seen by examining the first derivative:

D(expression(a * X^b), "X")

a * (X^(b - 1) * b)

The power function is mathematically equivalent to an exponential function of , since:

The power function (named as the Freundlich equation) is widely used in agricultural chemistry, for example to model the sorption of xenobiotics in soil. It is also used to describe the relationship between the number of plant species and the sampled area (the species–area relationship). Check the post at this link for an example of fitting.

Logarithmic function

This function is linear in :

Because of the logarithmic transformation, the function is only defined for . The parameter determines the shape of the curve. When , the response increases with and the curve is concave down. Conversely, when , decreases as increases and the curve is concave up, as shown in Figure 1.

The interpretation of the parameters is fairly straightforward. The parameter is the value of the response when . Changing while keeping constant simply shifts the curve vertically without altering its shape, so that the difference between any two curves remains constant for every value of (Fig. 2).

The parameter controls the slope of the curve and, more specifically, it is equal to the slope at , as can be seen from the expression for the first derivative:

D(expression(a + b*log(X)), "X")

b * (1/X)

Changing while keeping constant produces curves that all intersect at (Fig. 2).

In biology, logarithmic functions are used, for example, to describe species–area relationships in ecology, enzyme kinetics in biochemistry, and sensory perception in neurobiology. Check the post at this link for an example of fitting.

Rectangular hyperbola

The rectangular hyperbola is commonly known as the Michaelis-Menten function and is usually parameterised as

The curve is monotonically increasing and concave down, approaching the horizontal asymptote (plateau) . It passes through the origin of the axes ( when ). The parameter is the value of that produces a response equal to . Indeed,

which is readily solved to obtain .

The first derivative is:

D(expression( (a*X) / (b + X) ), "X")

a/(b + X) - (a * X)/(b + X)^2

From this expression, we can see that the initial slope (at ) is .

Equation 14 is not defined for and has no biological meaning for . In most biological applications, however, both and are positive.

An equivalent parameterisation is:

In this form, the parameters and are often replaced by and , respectively. Although the response is not defined for , it tends to zero as .

Another common parameterisation includes the initial slope as an explicit parameter because of its biological relevance. This is obtained by dividing both the numerator and denominator by and noting that , so that :

The Michaelis-Menten function is widely used in pesticide chemistry, enzyme kinetics (Eq. 14), biological assay models (Eq. 15), and weed competition studies, where it is used to describe crop yield loss as a function of weed density (Eq. 16). An example of fitting is shown in another post, at this link.

Sigmoidal functions

Sigmoidal functions are S-shaped, with two horizontal asymptotes (a lower and an upper asymptote) and an inflection point. They can be parameterised in countless ways, which may be confusing when selecting the most appropriate model for a particular biological process. In this post, I will use the parameterisations proposed in Ritz et al (2019), which is very consistent and puts all the curves on the same ground. Furthermore, the fitting properties are good, which contribute to good convergence and reliable estimation.

For the sake of simplicity, it is useful to recognise that most sigmoidal models belong to one of three basic families, which differ only in the position of the inflection point relative to the two asymptotes:

1. Logistic: the ordinate of the inflection point lies halfway between the lower and upper asymptotes; consequently, the S-shape is symmetric (the concave-up and the concave-down portions ‘mirror’ each other).
2. Gompertz: the ordinate of the inflection point lies closer to the lower asymptote, and, thus, the concave-up portions is ‘shorter’ than the concave-down portion.
3. Modified Gompertz: the ordinate of the inflection point lies closer to the upper asymptote, and, thus, the concave-up portion is ‘longer’ than than the concave-down portion.

Each of these three equations can describe either increasing or decreasing responses, depending on the sign of one of the model parameters (see below). Their different shapes are compared side by side in Fig. 3.

Each of these three families (logistic, Gompertz and modified Gompertz) can be transformed by replacing with , giving raise to three corresponding families:

1. Log-logistic
2. Type-I Weibull
3. Type-II Weibull

Like their linear-scale counterparts, these three families can describe either increasing or decreasing responses, depending on the sign of the same model parameter (see below). In the end, we have twelve different sigmoidal functions to choose from. The choice among them depends primarily on the expected symmetry of the response and on whether the predictor is more naturally interpreted on a linear or logarithmic scale.

Thanks to their versatility, sigmoidal functions have been used for an uncountable number of biological applications. Just to mention a few examples in agriculture, I would like to cite plant growth, dose-response curves in biological assays and the time-course of seed germination. For some examples, take a look at my post at this link.

Logistic function

The logistic curve is derived from the cumulative logistic distribution function. It is symmetric about its inflection point and can be parameterised as:

where is the upper asymptote, is the lower asymptote, is the value of at the inflection point, and is proportional to the slope at the inflection point. Because the curve is symmetric, also represents the value of that produces a response halfway between and (commonly referred to as the ED50 in biological assays). The parameter may be either positive or negative and, consequently, the response may either increase or decrease as increases.

This equation is known as the four-parameter logistic model. If appropriate, constraints can be imposed on the parameter values. For example, may be fixed at 0, yielding the three-parameter logistic model. If, in addition, is fixed at 1, the model reduces to the two-parameter logistic model.

Gompertz function

The Gompertz curve can be parameterised in many different ways. I prefer a parameterisation that closely resembles that of the logistic function:

Unlike the logistic function, the Gompertz curve is not symmetric about its inflection point. It exhibits a longer lag phase at the beginning, followed by a progressively steeper increase before gradually approaching the upper asymptote. The parameters have essentially the same interpretation as those of the logistic function, except that , the abscissa of the inflection point, does not correspond to the value of producing a response halfway between and .

As with the logistic function, four-, three-, and two-parameter Gompertz models can be obtained by constraining one or both asymptotes.

Modified Gompertz function

We have seen that the logistic curve is symmetric about its inflection point, whereas the Gompertz curve exhibits a longer lag phase at the beginning, followed by a progressively steeper increase. A different asymmetric pattern can be obtained by modifying the Gompertz function as follows:

The resulting curve increases rapidly at the beginning but gradually slows down as it approaches the upper asymptote. As with the logistic and Gompertz functions, one or both asymptotes can be constrained ( and/or ), giving rise to four-, three-, and two-parameter versions of the modified Gompertz model.

Log-logistic function

The log-logistic curve is symmetric with respect to . A log-normal curve has a very similar shape, although it is used much less frequently. In biological assays (and also in germination studies), the log-logistic function is commonly parameterised as

The parameters have the same interpretation as in the logistic function. In particular, represents the value of that produces a response halfway between and (the ED50). It is easy to show that the above equation is equivalent to

Like the logistic function, the log-logistic model can be fitted in four-, three-, or two-parameter versions by constraining one or both asymptotes. It is widely used to describe dose-response relationships in biological assays, seed germination, and crop growth.

Weibull function (type 1)

The Type I Weibull function is the logarithmic counterpart of the Gompertz function, being defined on rather than on . It is parameterised as

The parameters have essentially the same interpretation as those of the other sigmoidal functions presented above. In particular, and are the lower and upper asymptotes, respectively, while is the value of corresponding to the inflection point. Unlike the log-logistic function, however, does not correspond to the ED50.

Weibull function (type 2)

The Type II Weibull function is closely related to the Type I Weibull function, but it describes a different type of asymmetry, analogous to that of the modified Gompertz function:

The parameters have the same interpretation as those of the Type I Weibull function. In particular, and are the lower and upper asymptotes, respectively, while is the value of at the inflection point. As with the Type I Weibull function, does not correspond to the ED50. One or both asymptotes may be constrained, giving rise to four-, three-, and two-parameter versions of the model.

Another ‘flexible’ sigmoid

I would also like to mention another sigmoidal function that has been widely used in biology because of its supposed flexibility: the Richards function:

In this model, the degree of asymmetry is controlled by the parameter . When , the function reduces to the logistic model, whereas, as decreases below 1, the curve progressively resembles a Gompertz function. Figure 4 illustrates the different shapes obtained for different values of .

Despite this apparent flexibility, the Richards function has poor statistical properties for parameter estimation. In particular, the additional parameter often causes strong correlations among the estimates, making the fitting process unstable and the parameter estimates difficult to interpret. Ratkowsky (1990) described it as having “more undesirable nonlinear regression behaviour than almost any nonlinear regression model in common use.” For these reasons, the Richards function will not be considered further in this blog.

      As a kid, my baseball hero was Mickey Mantle. In his 18-year Major League career, he amassed tremendous statistics, and surely he was one of the game’s all-time greats. When Mantle retired, he was third on the all-time career home run list with 536 home runs, trailing only Babe Ruth (714) and Willie Mays (587). Unfortunately, he suffered numerous injuries, leaving fans to wonder what his numbers would look like if he had enjoyed a healthier career.

      For decades many people considered (and still consider) Babe Ruth to be the greatest all-around baseball player, because he was a star pitcher and a star hitter. The counter-argument to the Mantle “what-if” is: what if Ruth had not spent his first six years as a pitcher and instead had been a full-time hitter? (Similarly, what if Ted Williams had not spent some of his prime years in the military.)

      Nowadays Shohei Ohtani is certainly a once-in-a-generation two-way star. However, he spent his early prime years in Japan’s professional baseball, and he may not have enough time left in US baseball to amass these historic lifetime stats. Aaron Judge is another exceptional modern hitter; hopefully, he stays healthy so we can watch his career numbers climb.

      To compare Mantle and Ruth’s lifetime stats objectively, I turned to data. The Sean Lahman baseball dataset contains Major League player stats back to 1871, and is available in the R library lahman. This was a good opportunity to practice data manipulation using dplyr.

      Ruth played 102 more games than Mantle over a 22-year career (2,503 versus 2,401). His lifetime statistics eclipse Mantle’s in every offensive category except stolen bases — though Ruth did steal home 10 times and hit nearly twice as many triples, suggesting he was faster than most realize. His lifetime batting average was an impressive .342, sitting just behind Tris Speaker and Ted Williams.

      Lineup protection played a massive role for both men. Ruth generally batted directly before Lou Gehrig, and Mantle batted directly before Yogi Berra. While both Ruth and Mantle drew plenty of walks, presumably pitchers rarely chose to walk them intentionally just to face Gehrig or Berra.

      I was also curious about their defensive metrics. For a fair comparison, I filtered the data to isolate only their outfield appearances. There is no value in comparing an outfield throw to an assist Ruth made on a comeback grounder while pitching. Similarly, I excluded Mantle’s infield appearances; he primarily played first base in his final two seasons and filled in briefly at other infield spots early in his career.

      Ruth played 222 more outfield games than Mantle (2,241 versus 2,019). Their total outfield putouts were nearly identical, likely because Mantle played centerfield and covered more ground. However, Ruth recorded nearly twice as many outfield assists (204 versus 117), which aligns with the arm strength expected of a former pitcher. While Ruth’s 204 assists don’t match stars such as Roberto Clemente’s 266, it remains a highly respectable number.

      Ultimately, the data shows that even if you completely ignore his pitching career, Babe Ruth built an outstanding, standalone career as both a hitter and a fielder.

              
          MANTLE  RUTH
games       2401  2503
at_bats     8102  8398
runs        1677  2174
hits        2415  2873
doubles      344   506
triples       72   136
home_runs    536   714
BA         0.298 0.342
rbi         1509  2217
sb           153   123
bb          1733  2062
so          1710  1330
games_of    2019  2241
putouts     4438  4444
assists      117   204

      Here is my R code:

library(Lahman)
library(tidyverse)

data(People)
which(People$nameLast == "Mantle")   # 11773,  playerID = mantlmi01
which(People$nameLast == "Ruth")   # 16674 , playerID = ruthba01 
data(Batting)    
data(Fielding)   # for fielding, want games in outfield POS == 'OF'

df_fielding <- Fielding %>% 
  filter(playerID %in% c('mantlmi01', 'ruthba01'), POS == 'OF') %>% 
  group_by(playerID) %>% 
  summarize(games_of = sum(G, na.rm = TRUE), putouts = sum(PO, na.rm = TRUE),
    assists = sum(A, na.rm = TRUE)
  )

df <- Batting %>% 
  filter(playerID %in% c('mantlmi01', 'ruthba01')) %>% 
   group_by(playerID) %>% 
   summarize(games = sum(G, na.rm = TRUE), at_bats = sum(AB, na.rm = TRUE),
                runs = sum(R, na.rm = TRUE), hits = sum(H, na.rm = TRUE),
                doubles = sum(X2B, na.rm = TRUE), 
                triples = sum(X3B, na.rm = TRUE), home_runs = sum(HR, na.rm = TRUE),
                BA = round(hits / at_bats,3), 
                rbi = sum(RBI, na.rm = TRUE),  sb = sum(SB, na.rm = TRUE),
                bb = sum(BB, na.rm = TRUE),  so = sum(SO, na.rm = TRUE),
                )  %>%
  merge(df_fielding, by = "playerID", all.x = TRUE) 

df_transposed <- as.data.frame(t(df))
colnames(df_transposed) <- c("MANTLE", "RUTH")
df_transposed <- df_transposed[-1, ]
df_transposed

End

Remember the “AI rat penis” incident at Frontiers in Cell & Developmental Biology?

Well, that was two and a half years ago. That’s enough water under the bridge to now check in and see how this incident affected the journal. If we look at how many papers published, we should be able to see whether the incident affected submissions.

I will say upfront that this analysis did not reveal what I thought it would. It’s well known that Frontiers journals had huge growth some years ago and I also knew that submissions had decreased in the last few years. So, I expected that this drop off was due to authors not wanting to publish in the journal after the AI figure debacle or perhaps a combination of that and the study by Mark Hanson and colleagues showing suspiciously fast turnaround times at journals, including titles from this publisher.

The plots can be quickly generated using {PubMedLagR} which is available here.

Jump to the code, or just look at the plots.

If we pull all the articles in PubMed for this journal and see how many were published per month, we get a plot like this. The dotted line shows the date of the rat penis incident.

So there was a huge boom in publications in the COVID-19 pandemic era and then a decrease. This decrease predates the rat penis incident. If anything, the number of papers published began to stabilise and perhaps even recover following the incident!

It’s a bit hard to see if the composition of articles has changed over time from the plot above. So let’s replot these data and look at the fraction of all articles that are papers versus other types. I did some spot-checking and there were several papers that are clearly Commentary, Editorial and Review types that are not properly tagged on PubMed. So take this with a pinch of salt. The articles are split 50:50 between papers and other types and this has not really changed. It could be argued that commissioned content might be unaffected whereas as directly submitted papers might drop off if authors were concerned about a journal. But there’s no evidence of that here.

This is looking at published papers. Obviously this metric lags behind the author behaviour we’re interested in. Did the rat penis incident affect submissions to the journal? We can look at the received dates of published papers:

Pretty much the same story: month-on-moth submission are steadily rising after the incident.

Obviously we only have submission data for papers that were eventually published in the journal. We don’t know what the rejection rate was over this period. It would be possible for a beleaguered journal to reduce the rejection rate to stabilise the numbers of papers it publishes. But if we assume the rejection rate is constant, then submissions follow similar dynamics.

Conclusion: the AI debacle didn’t negatively affect submissions at this journal.

Maybe potential authors were unaware of this issue (although there was plenty of publicity), or maybe they view the incident as a one-off, or perhaps they’re not concerned about science integrity issues.

It’s no secret that the biggest drivers of journal choice for authors are 1) the journal impact factor and 2) the turnaround time/hassle to publish the work. The journals with the highest impact factor and the fastest turnaround win big. A scandal like this doesn’t seem to have affected journal choice. Of course, this is just a simple analysis and there are many factors at play here.

Isn’t it ironic though that authors want to publish in a journals with high impact factors? They want their paper to be viewed as potentially highly citable because of where it is published; but if that means that their paper appeared in a journal that has a scientific integrity issue, the cry is “judge the paper on its own merits, not where it’s published”…

The code

library(PubMedLagR)
library(ggplot2)

retrieve_journal_year_records("Front Cell Dev Biol", 2015:2026, batch_size = 200, papers_only = FALSE)
all <- pubmed_xmls_to_df(clean = FALSE)
# remove duplicate rows
all <- all[!duplicated(all),]
# clear Data/ and then
retrieve_journal_year_records("Front Cell Dev Biol", 2015:2026, batch_size = 200)
pprs <- pubmed_xmls_to_df()
# classify the papers in all
all$paper <- ifelse(all$pmid %in% unique(pprs$pmid),"paper", "other")
# convert publication date to the first of the month
all$year_month <- as.Date(paste0(substr(all$pubdate,1,7),"-01"))

# frontiers colours for fun
x <- c("242 130 91", "24 157 196")
fcolours <- sapply(strsplit(x, " "), function(x)
  rgb(x[1], x[2], x[3], maxColorValue=255))

# ggplot of number of papers per month-year
p1 <- ggplot(all, aes(x = year_month, fill = paper)) +
  geom_bar() +
  scale_fill_manual(values = fcolours) +
  geom_vline(xintercept = as.POSIXct(as.Date("2024-02-14")), linetype=2) +
  labs(x = "", y = "Papers per month") +
  theme_classic()
p1

# stacked plot version
p2 <- ggplot(all, aes(x = year_month, fill = paper)) +
  geom_bar(position = "fill") +
  scale_fill_manual(values = fcolours) +
  geom_vline(xintercept = as.POSIXct(as.Date("2024-02-14")), linetype=2) +
  labs(x = "", y = "Composition") +
  theme_classic()
p2

# what about submissions?
all$sub_year_month <- as.Date(paste0(substr(all$recdate,1,7),"-01"))
p3 <- ggplot(all, aes(x = sub_year_month, fill = paper)) +
  geom_bar() +
  scale_fill_manual(values = fcolours) +
  geom_vline(xintercept = as.POSIXct(as.Date("2024-02-14")), linetype=2) +
  labs(x = "", y = "Papers per month") +
  theme_classic()
p3

I added some styling to the plots to post them here, but the code will give the same three plots as in the post.

—

The post title comes from “New Frontier” by Donald Fagen from his 1982 “The Nightfly” LP.

controller is an R package for working with controlled vocabularies. It’s first release (v0.1.0) is now available now on CRAN.

The package addresses something I find myself doing very often in analysis code: tidying messy and inconsistent terminologies. For smaller datasets, dplyr::recode() is okay for this, but writing the mapping out as an R function call gets tedious fast when dealing with a long list of terms. It becomes very tedious when you have variants distinguished only by things like capitalisation (OxA- vs. oxa-), word boundaries (Çatalhöyük vs. Çatal Höyük) or character encoding (ʿAin Ghazal vs. ʽAyn Ghazal).

controller instead defines preferred terms and their variants in a data frame. Its control() verb is the equivalent of dplyr::recode() but using this thesaurus and with a few extra bells and whistles for fuzzy matching and reporting what was (and wasn’t) changed:

library(controller)
data("colour_thesaurus")

shades <- c("daffodil", "purple", "magenta", "azure", "navy", "violet")
control(shades, colour_thesaurus)
#> Replaced values:
#> ℹ daffodil → yellow
#> ℹ azure → blue
#> ℹ navy → blue
#> ℹ violet → purple
#> Warning: Some values of `x` were not matched in `thesaurus`:
#> ✖ magenta

Fuzzy matching means we don’t need to exhaustively list those variants from things like differences in case, word boundaries, or character encoding:

control_ci(toupper(shades), colour_thesaurus)
#> Replaced values:
#> ℹ DAFFODIL → yellow
#> ℹ PURPLE → purple
#> ℹ AZURE → blue
#> ℹ NAVY → blue
#> ℹ VIOLET → purple
#> [1] "yellow"  "purple"  "MAGENTA" "blue"    "blue"    "purple"
#> Warning message:
#> Some values of `x` were not matched in `thesaurus`:
#> ✖ MAGENTA

This package has been hanging around for a while! It started off as a helper function I used for cleaning up site names from prehistoric sites in Southwest Asia. The basic idea was inspired by similar functions that used to exist in c14bazAAR for cleaning sample metadata for radiocarbon date, that I thought were quite neat. So when the maintainers of that package decided to deprecate those, I took over the thesauri as part of c14 and spun the supporting functions off into controller as a standalone package. Then over the years it acquired some more functionality for working with controlled vocabularies (a surprising gap in the R ecosystem), like reading heritage vocabularies in Historic England’s FISH format. Five years later, I am finally getting around to releasing it on CRAN because I need to release c14 on CRAN, because that’s used in analyses I’m now publishing. It’s the research software engineering of changing a lightbulb, basically.

The first release of controller includes:

• control(), control_ci(), and control_fuzzy() for recoding values
• control_names(), control_names_ci(), and control_names_fuzzy() for recoding names
• control_matches() for inspecting how matches were made
• read_fish() for reading vocabularies in Historic England’s FISH format
• colour_thesaurus, an example dataset

You can install it from CRAN:

install.packages("controller")

Or the development version from GitHub:

remotes::install_github("joeroe/controller")

You can find the full documentation at https://controller.joeroe.io.

muttest so far only printed to the console which mutants were killed and which survived.

That’s the minimal setup you need to act on mutation testing findings, but it’s far from optimal. Clear the terminal, close the IDE, and the results are gone unless you saved the output.

muttest@0.3.0 fixes that with an HTML report.

1. Run your tests with new JSON reporter.
2. Render the report.
3. Open it in your browser, and every mutant is laid over your source code.

Click the ones that survived and see exactly what slipped past your suite. That list of survivors is your to-do list: every one is a candidate test. You decide which are worth writing and which are just noise.

Read the full changelog here.

New to mutation testing?

Coverage tells you which lines ran. It says nothing about whether your tests would notice if those lines were wrong. You can delete every assertion in your suite, run covr, and still see 100%.

Mutation testing asks the harder question: if this code were subtly broken, would a test fail? It makes small changes to your source (> becomes >=, TRUE becomes FALSE) and reruns your tests against each version. Each changed version is a mutant. If a test fails, the mutant is killed. If every test still passes, the mutant survived, and you’ve found a gap.

A surviving mutant isn’t a bug to fix. It’s a missing test that could hide an expensive bug later.

If you want to dive deeper, the 0.2.0 post walks through it with an example.

The report

Run your mutation tests with JSONMutationReporter, call report(), and you get a self-contained HTML file.

See the full example report here.

Here’s what the report does that the console can’t do:

• Your source code annotated with mutants. Every file is shown with its actual code, and each mutated line is marked inline. You read the survivor in context instead of reconstructing where a given line was.
• Per-file scores. The overall number is a starting point, not a diagnosis. The report breaks the score down by file, so you can see your parser sitting at 95% while the validation module is at 40%. That tells you which file to open first.
• Per-line mutant diffs. Click a mutant and see exactly what changed: the original line, the mutated line, and whether your tests killed it. When something survived, the diff tells you what kind of a potential bug slipped through.
• Filter by status. The report opens focused on survivors, since the killed mutants already did their job. You can bring the others back when you want the full picture, but the default is the view you actually came for.
• Keyboard navigation. Step through mutants without scrolling through the whole file.

That’s the workflow – filter to survivors, walk the list top to bottom, and each one is a decision: write the test, or call it noise and move on.

A checklist you can work through.

JSON output

The new JSONMutationReporter writes results to a file that conforms to the mutation-testing-elements schema. That schema isn’t something I invented. It’s the format shared by StrykerJS, Stryker.NET, Stryker for Scala, and others. If other teams in your organization already use Stryker and its reporting tools, you can plug muttest results straight into the same dashboards.

library(muttest)

plan <- muttest_plan(
  source_files = "R/shipping.R",
  mutators = comparison_operators()
)

muttest(plan, reporter = JSONMutationReporter$new())

report()

JSONMutationReporter writes to muttest.json by default, and report() reads that same file and writes muttest.html next to it, so with the defaults there’s nothing to wire up.

If you don’t like the built-in report, the JSON file is there to build your own. I plan to keep this schema stable so don’t worry about your report breaking with every next release.

Watch it run, save it for later

You usually want both: a live progress table while the run happens, and a saved file to open afterward. Before, you picked one reporter and got one behavior.

MultiReporter runs several reporters at once. Feed it a ProgressMutationReporter and a JSONMutationReporter together. The console fills up as mutants are tested, and the JSON file is waiting when it’s done.

muttest(
  plan,
  reporter = MultiReporter$new(
    ProgressMutationReporter$new(),
    JSONMutationReporter$new()
  )
)

It’s exactly the same idea as testthat’s MultiReporter.

Scoring changes worth knowing about

Two changes in this release will move your score:

• Errored mutants now count as killed. Some mutations don’t just fail an expectation; they make the code throw before a test can even assert anything. A mutation that produces an error is still a mutation your suite detected, so 0.3.0 scores it as killed. This matches PIT, Stryker, and mutmut. Previously muttest counted these as survived, which was too harsh: the test did catch the change, it just caught it by blowing up.
• Uncovered code gets its own category. Under FileTestStrategy, a mutant in a source file with no matching test file used to look identical to a mutant your tests ran and missed. Those are different problems. A survived mutant means your tests are weak; an uncovered mutant means there are no tests at all. 0.3.0 reports uncovered mutants as a separate no coverage category and leaves them out of the score, so untested code stops masquerading as escaped mutants and dragging your number down for the wrong reason.

If your score shifts after upgrading, this is why. It’s more honest this way.

Tell me about your experience of using muttest

I hope this release makes it easier to analyze results.

muttest is still young and the interface may shift. If the report is missing something you’d want to see, if a mutation you care about has no way to be expressed, or if the new scoring surprises you in a way that feels wrong, open an issue on GitHub.

Feature requests are just as welcome as bug reports.

We all were presented with this problem – water jug problem, which – in the time of Football world cup 2026 – can be translated to any liquid. *hint hint* But the riddle is as logical as mathematical. With mathematics finding the greatest common divisor. In general, it can be used with state search or Depth First Search (DFS).

With DFS we can solve this with Bezout’s identity; which in general is a theorem which relates two arbitraty integers with their greatests common divisor; and used in algebraic language, finding common zeros of n-polznomials in n-indeterminates. So the common zeros equals the product of the degrees of the polynomials.

And now that the imagine splitting the 16L beer Giraffe in two 8L Giraffes but you are only using 11L and 7L empty giraffes

This is the proof that with beer, Algebra is more fun And because it is fun, we can also find greated common divisors using Breadth-First Search (BFS). And here is the code:

solve_jugs <- function(caps = c(16, 11, 7), start = c(16, 0, 0), goal = c(8, 8, 0)) {
  # BFS over all (a, b, c) states
  # Each state is a named integer vector of water = amount is each jug or ?????
  
  queue   <- list(list(state = start, path = list(start)))
  visited <- list()
  key     <- function(s) paste(s, collapse = "-")
  
  while (length(queue) > 0) {
    node  <- queue[[1]]
    queue <- queue[-1]
    s     <- node$state
    
    if (isTRUE(all(s == goal))) return(node$path)
    if (!is.null(visited[[key(s)]])) next
    visited[[key(s)]] <- TRUE
    
    n <- length(s)
    for (from in 1:n) {
      for (to in 1:n) {
        if (from == to || s[from] == 0 || s[to] == caps[to]) next
        
        pour    <- min(s[from], caps[to] - s[to])
        new_s   <- s
        new_s[from] <- s[from] - pour
        new_s[to]   <- s[to]   + pour
        
        if (is.null(visited[[key(new_s)]])) {
          queue <- c(queue, list(list(
            state = new_s,
            path  = c(node$path, list(new_s))
          )))
        }
      }
    }
  }
  NULL  # no solution; add message or smht :)
}

solution <- solve_jugs()

 
for (step in solution) {
  cat(sprintf("  %-4d  %-4d  %-4d\n", step[1], step[2], step[3]))
}
  

And the final solution will reveal the steps and actions:

  16    0     0   
  5     11    0   
  5     4     7   
  12    4     0   
  12    0     4   
  1     11    4   
  1     8     7   
  8     8     0   
  

Similar steps are presented on the animation above.

As always, the complete code is available on GitHub in  Useless_R_function repository. And code to the animation is here (same Github repository).

Check the repository for future updates!

Stay healthy, hydrated and happy R-coding!

At useR! 2026 in Warsaw, I contributed to a workshop on Geocomputation with R (July 6, 2026) and gave a keynote entitled A world still to be mapped: reflections on geocomputation in R (July 8, 2026). Both were opportunities to consider what the R spatial ecosystem already makes possible, its strengths, limitations, and what remains to be done.

Keynote slides: https://jakubnowosad.com/user2026/

Workshop materials: https://github.com/geocompx/user26

@online{nowosad2026,
  author = {Nowosad, Jakub},
  title = {A World Still to Be Mapped: Reflections on Geocomputation in
    {R:} Takeaways from the Talk and Workshop at {UseR!} 2026},
  date = {2026-07-21},
  url = {https://jakubnowosad.com/posts/2026-07-21-user/},
  langid = {en}
}

TLDR: The first post established that the WGI and WJP rule of law measures agree on how they rank countries. This follow-up tests the stricter question of whether the WGI can predict the WJP’s levels, rather than just the relative ranking / relative position of countries. It can: forecasting errors appear to be minimal and don’t vary significantly over time.

This is second post in a series examining relationships between the rule of law and economic and social outcomes.

In the previous post in the series I tested whether the World Governance Indicator’s (WGI) rule of law index (RoL) was a worthy substitute for the index produced by the World Justice Project when needing to take advantage of the WGI’s wider and longer coverage.

When I asked Claude for feedback it rudely provided it, arguing that my analysis was fine for confirming general agreement as to how the two measures ranked countries, but additional analysis was needed to test whether the WGI could reliably be used to predict RoL levels.

My first response was to tell Claude that nobody should be using composite indices in this way in the first place as it implied a level of precision they don’t have, but I decided it would be a useful addendum to the first post that I could post during my holiday in the Northern Territory.

Project Setup and Data

Once again, data used in this post can be downloaded here for the WGI and here for the WJP’s RoL index. These were current as of July 2025.

#load the packages we'll probably need
library(tidyverse)
library(readxl)
library(janitor)
library(countrycode)

#import WGI data
dta_wgi_2025<-read_excel("./Data/wgidataset_with_sourcedata-2025.xlsx",
                         sheet="rl") |> 
  clean_names()

#import World Justice Project RoL data
dta_wjp_rol<-read_excel('./Data/2025_wjp_rule_of_law_index_HISTORICAL_DATA_FILE.xlsx', sheet="Historical Data")|> 
  clean_names()

Data Cleaning

Data cleaning is identical to the last post and revolves around making sure that countries names and regions are consistently applied across the two indices.

#standardize column names
dta_wgi_2025<-dta_wgi_2025 |> 
  rename(country=economy_name,
         iso3c=economy_code,
         wgi_rol=governance_estimate_approx_2_5_to_2_5)

dta_wjp_rol<-dta_wjp_rol |> 
  rename(iso3c=country_code) |> 
  rename_with(~ str_replace(., "^x", "factor_"), starts_with("x"))

#change wjp's year variable to YYYY format and convert to numeric 
#(adopts the first 4 digit year when in YYYY-YYYY format)
dta_wjp_rol <- dta_wjp_rol |> 
  mutate(year = as.numeric(str_sub(year, 1, 4)))

#cold-heartedly drop columns I'm not interested in 
dta_wgi_2025<-dta_wgi_2025 |> 
  select(iso3c, income_classification, year,wgi_rol)
#drop country and region name labels so these can be standardized 
dta_wjp_rol<-dta_wjp_rol |> 
  select(-country_year,-country,-region) |> 
  rename(wjp_rol=wjp_rule_of_law_index_overall_score)


#merge dataframes
dta_rol_unified<-left_join(dta_wgi_2025,
                           dta_wjp_rol,
                           by = join_by(year, iso3c), 
                           keep=FALSE) 

#add standardized and region names country names

#define country code assignments for legacy / ambigious codes
#(Note: matches devised by Claude) 
ref_iso3c_custom_names <- c(ADO = "Andorra",
                            ANT = "Netherlands Antilles",
                            PRI = "Puerto Rico",
                            REU = "Réunion",
                            XKX = "Kosovo")

ref_iso3c_custom_regions <- c(ADO = "Europe & Central Asia",
                              ANT = "Latin America & Caribbean",
                              PRI = "Latin America & Caribbean",
                              REU = "Sub-Saharan Africa",
                              XKX = "Europe & Central Asia")

#assign country names and regions:
dta_rol_unified<-dta_rol_unified |> 
  mutate(country_name=countrycode(iso3c, 
                                  origin='iso3c',
                                  destination = 'country.name.en',    
                                  custom_match = ref_iso3c_custom_names),
         region=countrycode(iso3c, 
                            origin='iso3c',
                            destination = 'region',    
                            custom_match =ref_iso3c_custom_regions))

Testing Cardinal Substitution

The code below basically asks whether it’s possible to estimate the level of WJP’s RoL index using the WGI’s RoL measure using a simple linear model. At the outset, the model’s fit is quite high at 97 percent, but we could have guessed this from the cross-country scatter we produced in our last post.

# Cardinal substitution: does a WGI->WJP mapping reproduce WJP levels? -----
sum_rol_country <- dta_rol_unified |>
  filter(!is.na(wgi_rol), !is.na(wjp_rol)) |>
  summarise(wgi_rol_mean = mean(wgi_rol),
            wjp_rol_mean = mean(wjp_rol),
            .by = country_name)

mod_rol_calib <- lm(wjp_rol_mean ~ wgi_rol_mean, data = sum_rol_country)

#output the R squared
summary(mod_rol_calib)$r.squared |> round(2)

The second chunk calculates the Leave-One-Out Cross-Validation Root Mean Squared Error (LOO RMSE), which estimates prediction error for each country when it’s held out of the sample: a smaller LOO RMSE implies smaller forecast errors when predicting outside the sample. This is then scaled by the standard deviation of the country-level WJP means to make it easier to interpretat.

The results are encouraging: the LOO RMSE indicates that using the WGI to estimate WJP levels carries a typical error of less than three percent, which is roughly 19% of the cross-country spread in WJP scores. Suggesting the WGI can provide a reasonable proxy when the levels are important too.

# Leave-One-Out Cross-Validation Root Mean Squared Error (LOO RMSE)
rlt_rol_loo_rmse   <- sqrt(mean((residuals(mod_rol_calib) /
                                   (1 - hatvalues(mod_rol_calib)))^2))

#divide by standard deviation of WJP index means to put the figure in the context fo the data
rlt_rol_rmse_ratio <- rlt_rol_loo_rmse / sd(sum_rol_country$wjp_rol_mean)

Systematic Bias

Of course, because the LOO RMSE measures is calculated across the entire dataset, it might hide systematic bias, such as if the forecast accuracy varies depending on the WGI score.

This is tested in the plot below, by comparing the prediction error (residual) with the predicted value for the WJP’s index (using the WGI RoL measure). If there no systematic bias, the residuals should remain relatively stable across the predicted WJP values.

Alas, this doesn’t appear to be the case, with the model under-predicting the WJP’s RoL index at both ends. Although this bias is relatively mild at 0.05 (in WJP index units), it does suggest the WGI is likely to be a poorer substitute for the WJP’s index for countries at either end of the RoL spectrum.

dta_plt_rol_calib <- sum_rol_country |>
  mutate(wjp_pred = fitted(mod_rol_calib),
         resid    = wjp_rol_mean - wjp_pred)

plt_rol_calib <- ggplot(dta_plt_rol_calib, aes(wjp_pred, resid)) +
  geom_hline(yintercept = 0, linewidth = 0.3, col = "grey70") +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "loess", se = FALSE, col = "grey40", linewidth = 0.5) +
  labs(x = "Predicted WJP (from WGI)", y = "Residual (WJP units)") +
  theme_classic()

plt_rol_calib

However, being an applied economist means I know a simple trick: add non-linear terms to everything and see if the smile disappears.

And it does, the smile is now a grimace: as is shown in the plot below adding a non-linear term seems to take care of the bias at either end. Notice the average errors no longer increase at either end, making the prediction relatively constant (and small) across the RoL spectrum. Suggesting the WGI can be trusted to predict the WJP index at either end of the spectrum provided we account for the non-linearity.

#reestimate the model with a non-linear WGI term 
mod_rol_calib_non_linear <- lm(wjp_rol_mean ~ wgi_rol_mean + I(wgi_rol_mean^2),
                               data = sum_rol_country)


dta_plt_rol_calib_nl <- sum_rol_country |>
  mutate(wjp_pred = fitted(mod_rol_calib_non_linear),
         resid    = wjp_rol_mean - wjp_pred)

plt_rol_calib_nl <- ggplot(dta_plt_rol_calib_nl, aes(wjp_pred, resid)) +
  geom_hline(yintercept = 0, linewidth = 0.3, col = "grey70") +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "loess", se = FALSE, col = "grey40", linewidth = 0.5) +
  labs(x = "Predicted WJP (from WGI, quadratic)", y = "Residual (WJP units)") +
  theme_classic()

plt_rol_calib_nl

But does it drift?

The final question Claude suggested I ask was whether the WGI’s ability to predict the WJP’s index changes over time. The code below tests this by comparing a simple linear prediction model with one where the influence of the WGI changes over time.

The anova() function then compares results of both models to test whether allowing the WGI’s influence to change over time improves the model’s forecasting accuracy. If it does (indicated by a significant p-value), the WGI’s ability to predict the WJP’s index might not be stable over time.

The result: a p value of 9%, which points to the time invariant model being good enough for predicting the WJP’s index.

#estimate a simple linear model with WGI's influence kept constant over time
mod_rol_calib_flat <- lm(wjp_rol ~ wgi_rol, data = dta_rol_unified)

#estimate a model where WGI's influence varies over time
mod_rol_calib_yr   <- lm(wjp_rol ~ wgi_rol * factor(year), data = dta_rol_unified)

#compare the prediction capacity of both models 
anova(mod_rol_calib_flat, mod_rol_calib_yr)

Look what Claude made me do…

So there you have it: the WGI can serve as a useful substitute for the WJP’s rule of law index whether you’re interested in RoL rankings or levels, provided you account for the non-linearity at either end of the RoL scale.

How AI was used for this post: As this additional analysis was Claude’s suggestion and I was hours away from a holiday I leaned heavily on Claude for the first draft of the code. Much of the text is mine, but I did ask AI for suggestions to make sure I was explaining the analysis accurately.

The post Level With Me: Can the WGI Predict the WJP’s Rule of Law Scores? appeared first on Giles.

Usually, the first step in every nonlinear regression analysis is to select the function that best describes the phenomenon under study. The next step is to fit this function to the observed data, possibly by using some sort of nonlinear least squares algorithm.

We have already devoted a post to the first task, which you can find at this link.

As for the second step, the main problem is that nonlinear least squares algorithms are iterative, in the sense that they start from some initial guesses for the model parameters, which are continuously improved until the least squares solution is approximately reached. Quite often, providing such initial guesses for all model parameters becomes a problem: if our guesses are not close enough to the least squares estimates, the algorithm may stall and fail to converge. Or, even worse, it may converge to the wrong solution. How do we obtain good initial guesses for the model parameters? This is not easily accomplished, especially for students and practitioners. This is where self-starters come in handy.

Self-starter functions can automatically calculate initial values for any given dataset and, therefore, they can make nonlinear regression almost as straightforward as linear regression. From a teaching perspective, this means that the transition from linear to nonlinear models is immediate and hassle-free.

In another post, at this link, I explained how self-starters can be built for both the nls() function in the ‘stats’ package and the drm() function in the ‘drc’ package (Ritz et al., 2019). In this post, I would like to provide an overview of the self-starting functions that already exist in R, either in the ‘stats’, ‘drc’, or ‘statforbiology’ packages. I do not aim for completeness here, as other packages also contain self-starters, such as the ‘nlraa’ package (Miguez, 2025), with which I am not sufficiently familiar. The exemplary datasets are included in the ‘statforbiology’ package and they all come from real field or greenhouse experiments; further information and citations to the original works can be found in my book (Onofri, 2026), or in the additional meterial at this link.

Functions and curve shapes

As in the other post, I have classified linear/nonlinear regression functions according to the shape they show when they are plotted in an x-y graph, following the approach taken in Ratkowsky (1990):

• Polynomials
  1. Straight line function
  2. Quadratic polynomial function
• Concave/Convex curves (no inflection)
  1. Exponential function
  2. Asymptotic function / Negative exponential function
  3. Power function
  4. Logarithmic function
  5. Rectangular hyperbola
• Sigmoidal curves
  1. Logistic function
  2. Gompertz function
  3. Modified Gompertz function
  4. Log-logistic function
  5. Weibull (type-1) function
  6. Weibull (type-2) function
• Curves with maxima/minima
  1. Peaked sigmoidal function
  2. Bragg function
  3. Lorentz function
  4. Beta function

To make navigation easier, you can first inspect Figure 1 to identify the type of curve of interest and then use the links above to jump directly to the corresponding section.

First of all, we need to install (if necessary) and load these packages, by using the code below.

# installing package, if not yet available
# install.packages("statforbiology")

# loading package
library(statforbiology)

# Example of fitting four straight lines at once
library(statforbiology)
dataset <- getAgroData("metamitron")
model <- drm(log(Conc) ~ Time, fct = DRC.linear(),
             curveid = Herbicide, data = dataset)
summary(model)

Model fitted: Straight line (2 parms)

Parameter estimates:

        Estimate Std. Error t-value   p-value    
a:M    4.5159111  0.0577603  78.184 < 2.2e-16 ***
a:MP   4.6304853  0.0577603  80.167 < 2.2e-16 ***
a:MC   4.5401707  0.0577603  78.604 < 2.2e-16 ***
a:MPC  4.7546974  0.0577603  82.318 < 2.2e-16 ***
b:M   -0.0386165  0.0015585 -24.777 < 2.2e-16 ***
b:MP  -0.0318908  0.0015585 -20.462 < 2.2e-16 ***
b:MC  -0.0231879  0.0015585 -14.878 < 2.2e-16 ***
b:MPC -0.0234014  0.0015585 -15.015 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error:

 0.1687429 (88 degrees of freedom)

library(statforbiology)
dataset <- getAgroData("waterProductivity")

# nls fit
model <- nls(WP ~ NLS.poly2(Regime, a, b, c),
             data = dataset)

#drc fit
model <- drm(WP ~ Regime, fct = DRC.poly2(),
             data = dataset)

Model fitted: Second Order Polynomial (3 parms)

Parameter estimates:

                 Estimate  Std. Error  t-value   p-value    
a:(Intercept)  1.9187e+00  2.3171e-01   8.2806 1.679e-05 ***
b:(Intercept)  4.2977e-01  8.4554e-03  50.8287 2.217e-12 ***
c:(Intercept) -3.3618e-03  6.6586e-05 -50.4884 2.355e-12 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error:

 0.1441625 (9 degrees of freedom)

library(statforbiology)
data(degradation)

# nls fit
model <- nls(Conc ~ NLS.expoDecay(Time, a, k),
             data = degradation)

# drm fit
model <- drm(Conc ~ Time, fct = DRC.expoDecay(),
             data = degradation)
summary(model)

Model fitted: Exponential Decay Model (2 parms)

Parameter estimates:

                 Estimate Std. Error t-value   p-value    
C0:(Intercept) 99.6349312  1.4646680  68.026 < 2.2e-16 ***
k:(Intercept)   0.0670391  0.0019089  35.120 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error:

 2.621386 (22 degrees of freedom)

Time <- c(1, 3, 5, 7, 9, 11, 13, 20)
Weight <- c(8.22, 14.0, 17.2, 16.9, 19.2, 19.6, 19.4, 19.6)

# nls fit
model <- nls(Weight ~ NLS.asymReg(Time, b, m, a) )

# drm fit
model <- drm(Weight ~ Time, fct = DRC.asymReg(names = c("b", "m", "a")))
summary(model)

Model fitted: Asymptotic Regression Model (3 parms)

Parameter estimates:

               Estimate Std. Error t-value   p-value    
b:(Intercept)  3.756123   1.367835  2.7460  0.040500 *  
m:(Intercept)  0.337078   0.052704  6.3957  0.001385 ** 
a:(Intercept) 19.629817   0.420557 46.6758 8.525e-08 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error:

 0.639656 (5 degrees of freedom)

library(statforbiology)
speciesArea <- getAgroData("speciesArea")

#nls fit
model <- nls(numSpecies ~ NLS.powerCurve(Area, a, b),
             data = speciesArea)

# drm fit
model <- drm(numSpecies ~ Area, fct = DRC.powerCurve(),
             data = speciesArea)
summary(model)

Model fitted: Power curve (Freundlich equation) (2 parms)

Parameter estimates:

              Estimate Std. Error t-value   p-value    
a:(Intercept) 4.348404   0.337197  12.896 3.917e-06 ***
b:(Intercept) 0.329770   0.016723  19.719 2.155e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error:

 0.9588598 (7 degrees of freedom)

library(statforbiology)
speciesArea <- getAgroData("speciesArea")

#nls fit
model2 <- nls(numSpecies ~ NLS.logCurve(Area, a, b),
             data = speciesArea)

# drm fit
model2 <- drm(numSpecies ~ Area, fct = DRC.logCurve(),
             data = speciesArea)
summary(model2)

Model fitted: Linear regression on log-transformed x (2 parms)

Parameter estimates:

              Estimate Std. Error t-value   p-value    
a:(Intercept)  1.51111    1.17401  1.2871     0.239    
b:(Intercept)  4.06359    0.35575 11.4224 8.847e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error:

 1.910082 (7 degrees of freedom)

       df      AIC
model   3 28.52288
model2  3 40.92769

library(statforbiology)
dataset <- getAgroData("Ammi94_YL")

#drm fit
model <- drm(YL ~ Density, fct = DRC.YL(), data = dataset)
summary(model)

Model fitted: Yield-Loss function (Cousens, 1985) (2 parms)

Parameter estimates:

               Estimate Std. Error t-value p-value  
i:(Intercept) 0.0118359  0.0045637  2.5935 0.02228 *
A:(Intercept) 0.4841796  0.1863969  2.5976 0.02211 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error:

 0.05453432 (13 degrees of freedom)

Sigmoidal curves

In the previous post (at this link) we examined six different sigmoidal shapes based on , to be selected according to the expected symmetry of the response. The equations are the logistic, Gompertz, and modified Gompertz functions, with the following three equations, respectively:

The same curves can be based on the logarithm of , to obtain the corresponding log-logistic, Type-1 Weibull, and Type-2 Weibull equations, which are, respectively:

The meaning of the model parameters and details about the shapes can be found in another post at this link.

For all these four-parameter curves, it is possible to set to obtain the three-parameter curves (17a), (18a), (19a), (20a), (21a), and (22a). Likewise, we can constrain both and to obtain the two-parameter curves (17b), (18b), (19b), (20b), (21b), and (22b). The collection of self-starters is huge; therefore, I will provide a summary table.

List of sigmoidal functions, together with self-starters for fitting with drm() and nls(). (1) These self-starters are included in the ‘drc’ package and the corresponding curves have an increasing shape when and a decreasing shape when . For all other self-starters, which are available in the ‘statforbiology’ package, the reverse is true (which is the most common definition, apart from the case of herbicide bioassays){#tbl-models}

Name	Eq. No.	for drm()	for nls()
Logistic (4-parameters)	(17)	L.4() (1)	NLS.L4()
Logistic (3-parameters)	(17a)	L.3() (1)	NLS.L3()
Logistic (2-parameters)	(17b)	L.2() (1)	NLS.L2()
Gompertz (4-parameters)	(18)	G.4() (1)	NLS.G4()
Gompertz (3-parameters)	(18a)	G.3() (1)	NLS.G3()
Gompertz (2-parameters)	(18b)	G.2() (1)	NLS.G2()
mod-Gompertz (4-parameters)	(19)	DRC.E4()	NLS.E4()
mod-Gompertz (3-parameters)	(19a)	DRC.E3()	NLS.E3()
mod-Gompertz (2-parameters)	(19b)	DRC.E2()	NLS.E2()
Log-logistic (4-parameters)	(20)	LL.4() (1)	NLS.LL4()
Log-logistic (3-parameters)	(20a)	LL.3() (1)	NLS.LL3()
Log-logistic (2-parameters)	(20b)	LL.2() (1)	NLS.LL2()
Type-1 Weibull (4-parameters)	(21)	W1.4() (1)	NLS.W1.4()
Type-1 Weibull (3-parameters)	(21a)	W1.3() (1)	NLS.W1.3()
Type-1 Weibull (2-parameters)	(21b)	W1.2() (1)	NLS.W1.2()
Type-2 Weibull (4-parameters)	(22)	W2.4() (1)	NLS.W2.4()
Type-2 Weibull (3-parameters)	(22a)	W2.3() (1)	NLS.W2.3()
Type-2 Weibull (2-parameters)	(22b)	W2.2() (1)	NLS.W2.2()

Example 8

The dataset ‘beetGrowth’ in the ‘statforbiology’ package reports the results of an experiment in which the growth of sugar beet was evaluated under either weed-free or weed-infested conditions. In the following example, we consider the time course of crop weight under infested conditions and fit the logistic, Gompertz, and modified Gompertz functions to the observed data. We see that, based on the AIC, the three fits are equivalent in practice. We fit the same models with ‘drm()’ to show that the two fitting functions give the same results, apart from the sign of the parameter, which is reversed in the L.3() function.

library(statforbiology)
beetGrowth <- getAgroData("beetGrowth")

# nls fit
model.1 <- nls(Infested ~ NLS.L3(DAE, b, d, e), data = beetGrowth)
model.2 <- nls(Infested ~ NLS.G3(DAE, b, d, e), data = beetGrowth)
model.3 <- nls(Infested ~ NLS.E3(DAE, b, d, e), data = beetGrowth)
AIC(model.1, model.2, model.3)

        df      AIC
model.1  4 250.2906
model.2  4 250.7596
model.3  4 250.8456

summary(model.1)

Formula: Infested ~ NLS.L3(DAE, b, d, e)

Parameters:
   Estimate Std. Error t value Pr(>|t|)    
b 1.188e-01  1.874e-02   6.338 1.33e-05 ***
d 2.512e+03  1.277e+02  19.671 4.01e-12 ***
e 5.803e+01  1.823e+00  31.836 3.46e-15 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 222 on 15 degrees of freedom

Number of iterations to convergence: 6 
Achieved convergence tolerance: 4.299e-06

# drm fit
model <- drm(Infested ~ DAE, fct = L.3(), data = beetGrowth)
summary(model)

Model fitted: Logistic (ED50 as parameter) with lower limit fixed at 0 (3 parms)

Parameter estimates:

                Estimate Std. Error t-value   p-value    
b:(Intercept)   -0.11876    0.01832  -6.483 1.033e-05 ***
d:(Intercept) 2511.95664  127.95600  19.631 4.131e-12 ***
e:(Intercept)   58.03099    1.83458  31.632 3.788e-15 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error:

 221.9644 (15 degrees of freedom)

Example 9

The dataset ‘brassica’ in the ‘statforbiology’ package reports the results of a short-term pot bioassay to determine the phytotoxicity of imazethapyr in hydroponic solution (Onofri, 1994). The response variable is the weight of the test plant (Brassica rapa), as affected by the herbicide concentration in the nutrient solution, which requires a decreasing sigmoidal curve. In the code below, we fit both a 4-parameter log-logistic and a 3-parameter log-logistic (with ) function and graphically compare the two fits.

library(statforbiology)
brassica <- getAgroData("brassica")
model.1 <- drm(FW ~ Dose, fct = LL.4(), data = brassica)
model.2 <- drm(FW ~ Dose, fct = LL.3(), data = brassica)

library(statforbiology)
dataset <- getAgroData("SOLNI_hormesis")
model <- drm(FW ~ Dose, fct = BC.5(),
            data = dataset)
summary(model)

Model fitted: Brain-Cousens (hormesis) (5 parms)

Parameter estimates:

              Estimate Std. Error t-value   p-value    
b:(Intercept)  3.57034    1.25309  2.8492 0.0082868 ** 
c:(Intercept)  3.76335    0.85964  4.3778 0.0001619 ***
d:(Intercept) 14.14991    1.36933 10.3334 6.992e-11 ***
e:(Intercept)  2.85405    0.76931  3.7099 0.0009485 ***
f:(Intercept)  5.26894    3.38238  1.5578 0.1309350    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error:

 2.97547 (27 degrees of freedom)

library(statforbiology)
dataset <- getAgroData("Hordeum_GR")
model <- drm(GR50 ~ Temp, fct = DRC.bragg.4(),
            data = dataset)
summary(model)

Model fitted: Bragg equation with four parameters

Parameter estimates:

                Estimate Std. Error t-value   p-value    
b:(Intercept)  0.0052849  0.0036063  1.4654 0.2804269    
c:(Intercept)  0.0039006  0.0131491  0.2966 0.7947090    
d:(Intercept)  0.0290870  0.0010552 27.5646 0.0013135 ** 
e:(Intercept) 18.0013282  0.5220972 34.4789 0.0008401 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error:

 0.001544712 (2 degrees of freedom)

© 2025 tidyomics. Content is published under Creative Commons CC-BY-4.0 License for the text and BSD 3-Clause License for any code. | R-Bloggers

I recently read John Kay and Mervyn King’s book Radical Uncertainty as part of an economics book club and I thought I would share some thoughts on the book here.

I’ve never read anything by John Kay before, but this is the second book by King which I have read, after The End of Alchemy, which is about banking regulation. Radical Uncertainty is about statistics, and especially about the problems with assuming that we can make probabilistic calculations while neglecting the possibility of things which we haven’t even contemplated happening. I was pretty fascinated by the first half of the book, but it does become quite long-winded in the second half, and repeats many of its points over and over again.

The book makes a lot of good points, but I would also like to point out some mistakes which it makes when talking about statistical modelling.

In the language of the book, Radical Uncertainty means Knightian Uncertainty. Knightian Uncertainty is named after Frank Knight, who made a distinction between uncertainty and risk. The basic idea is that you have a probability space $(X, \mu, \Sigma)$ where $X$ is a set, $\mu$ is a probability measure on $X$ and $\Sigma$ is the sigma-algebra of events. There’s absolutely no reason not to assume that $X$ is finite, so really we have a finite set $X$ and a probability $\mu(x)$ attached to each element $x \in X$. Suppose we make a draw from $X$

According to Knight’s definition, risk is “we don’t know $x$” and uncertainty is “we don’t know $X$”. One point which the book likes to hammer on is that in a lot of real-life situations, we really can’t possibly know $X$. So probability calculations (such as those made by bankers prior to 2008) are invariably wrong since they ignore “off-model events”. In 2008, this led to disastrous consequences.

I would like to mention at this point that there is a third kind of uncertainty. What if you know $X$ but don’t know $\mu$? This kind of uncertainty arises in a lot of gambling problems, for example when a bookmaker wants to make odds on the World Cup. You know exactly which outcomes are possible (all possible tournament brackets) but you need to assign a probability to each one. I was not able to find a name for this kind of uncertainty, but since it is exactly the situation faced in Bayesian statistics, it could be called Bayesian Uncertainty.

Kay and King make a pretty convincing case that Knightian uncertainty is important. However, I do think they slip up when talking about probability. For example, describing a calculation by Nate Silver of the probability of the 9/11 attacks, they say (Chapter 11)

So Silver’s calculation of probability was meaningless.

It doesn’t make sense to say that a probability calculation is meaningless. Every probability calculation is conditional on a probability model $(X, \mu)$. Given a model, you can calculate the probability of any subset of $X$ (or any event in $\Sigma$ in the continuous case). So you could say that Silver’s model was meaningless. But, in the context of the model, any probability he calculated makes perfect sense.

There is a temptation (to which Kay and King appear to succumb) to assume that there are some sort of “true” probabilities of things which modellers are attempting to calculate. But there are no such true probabilities, because there are no true probability models. A probability model is just an expression of what you don’t know. Different people can have different models because they are ignorant of different things. Thus, they can assign different probabilities to the same event, and this is fine. In some toy examples, such as rolling a die, everybody probably agrees more or less on the same probability model. But not so in the real world.

A similar error crops up when describing Keynes’ objection to the Principle of Indifference in Chapter 1

If, to take an example, we have no information whatever as to the area or population of the countries of the world, a man is as likely to be an inhabitant of Great Britain as of France, there being no reason to prefer one alternative to the other. He is also as likely to be an inhabitant of Ireland as of France. And on the same principle he is as likely to be an inhabitant of the British Isles as of France. And yet these conclusions are plainly inconsistent. For our first two propositions together yield the conclusion that he is twice as likely to be an inhabitant of the British Isles as of France. Unless we argue, as I do not think we can, that the knowledge that the British Isles composed of Great Britain and Ireland is a ground for supposing that a man is more likely to inhabit them than France, there is no way out of the contradiction.

Here the original probability model is

If you add the knowledge that the British Isles consists of GB + Ireland, the of course $p(\text{British Isles}) = 2p(\text{France})$. There is no contradiction here, but merely a failure to imagine that you could be ignorant of the fact that BI = GB + Ireland in the first place. New knowledge leads you to update your model, as it should.

It also seems that Kay and King may not be familiar with Bayesian statistics. For example, right at the end of the book they quote Savage:

There is some temptation to introduce probabilities of a second order so that the person would find himself saying such things as ‘the probability that B is more probable than C is greater than the probability that F is more probable than G.’ But such a program seems to meet insurmountable difficulties […] once second order probabilities are introduced, the introduction of an endless hierarchy seems inescapable. Such a hierarchy seems very difficult to interpret, and it seems at best to make the theory less realistic, not more.

Well, the prior distribution used in every Bayesian statistics problem is a probability distribution over probability distributions, so you definitely can consider probabilities of probabilities!

Anyway, these are minor points. Kay and King make a convincing argument, at great length, that in general you don’t know $X$ and therefore can’t construct a very good probability model. I don’t have any problem with this. But what should you do if you can’t make a probability model? To their credit, they propose an alternative approach, which they describe as choosing a “reference narrative” and essentially acting as though it is going to happen, while also taking some time to think about what you will do if the reference narrative doesn’t happen.

This process of choosing a reference narrative could also be described in terms of probability. In this case, instead of trying to construct the whole probability space $(X, \mu)$, you simply consider some set of plausible outcomes $X$, choose whichever $x \in X$ is the most probable and assign a probability weight of $1$ to it. In other words, you just choose the likeliest outcome and assume that it’s going to happen. This is much easier than constructing a probability space, and I think that Kay and King make a convincing case/ that this is generally how people act in real life. As Stigler says in The Epic Story of Maximum Likelihood

At a superficial level, the idea of maximum likelihood must be prehistoric[…]

I wondered whether I could understand why this approach works. Why would it be better, in a world of radical uncertainty, to take the likeliest outcome and run with it rather than trying to be a Bayesian?

To make a calculation, we need some sort of measure of how well a model performs versus reality. I found a very interesting paper by Süskind which makes a distinction between probability (the degree to which we believe something) and verisimilitude (the degree to which something is actually true). Probability is conditional on a probability model; verisimilitude is a property of the model itself. One of Süskind’s examples concerns black swans. Suppose there are 99 white swans for every black swan. The model “all swans are white” is wrong, but it makes a correct prediction 99 out of 100 times. The model “all swans are green” never makes a correct prediction. Therefore, “all swans are white” is not a great model, but it does have high verisimilitude.

Süskind is at pains to point out that verisimilitude is not the same as predictive accuracy. For example, assuming that Paul Daniels has magical powers predicts the outcome of his tricks very well, but it is not a model with high verisimilitude because it is not close to the truth. This is sort-of similar to the idea in machine learning where accuracy often doesn’t measure the thing which you’re actually interested in, and it’s better to use other measures of predictive performance instead.

Let’s suppose the true world has $4$ possible states with true probabilities $(p_1, p_2, p_3, 0)$ (the reason for adding a fourth state is to have a “junk” state which a model is allowed to predict but which actually never happens).¹ Suppose $p_2 > p_1$. Suppose that observers have only observed states $1$ and $2$ in the past (due to non-stationarity, aka data drift, state $3$ has never happened yet, but in fact is quite likely). A Bayesian would get the probability distribution

A maximum-likelihoodist would get the distribution $(0, 1, 0, 0)$. And a conservative maximum-likelihoodist who assumes that their reference narrative (state $2$) is likely to happen but they should also make allowances for the unknowable might choose the distribution $(0, 1/2, 0, 1/2)$.

Now suppose we decide to measure the verisimilitude of a our probability model $(x_1, x_2, x_3, x_4)$ by

You can check that the Bayesian always does better than the maximum-likelihoodist. But the conservative maximum-likelihoodist who follows the Kay-King recipe (follow your reference narrative while trying to hedge against the unknowable) can actually get a higher verisimilitude than the Bayesian who fails to take radical uncertainty into account, for example if $p_1=1/30, p_2=2/30$.

p1 <- 1/30
p2 <- 2/30
p_true <- c(p1, p2, 1-p1-p2, 0)
1-sum((c(p1/(p1+p2), p2/(p1+p2), 0, 0) - p_true)^2) # -0.26
1-sum((c(0, 1, 0, 0) - p_true)^2) # -0.68 (worse)
1-sum((c(0, 0.5, 0, 0.5) - p_true)^2) # -0.25 (better!)

Although this example is rather contrived, hopefully it does show that in the presence of Knightian Uncertainty, it can indeed be better just to take the likeliest outcome and run with it (while looking over your shoulder at the same time).

Economic Forecasts

As you would expect, the book is particularly focussed on economics. One of its conclusions is that you can’t expect to make reliable economic forecasts due to radical uncertainty. Your forecast has to be “I don’t know” and your actions have to be in accordance with whatever scenario you think is likeliest, given whatever set of scenarios you can imagine.

This all seems very reasonable. Countless studies have shown that long-range economic forecasts are no better than random guessing. So why do we believe them? What are forecasts for?

An interesting answer to this problem is suggested in a book by Jens Beckert. The key idea is that economic activity depends on coordination. People can’t coordinate their actions unless they can agree on the future. It doesn’t matter what number they agree on; they just have to agree on something. Economic forecasts, even when they are meaningless, give people something to agree about. Ancient priests predicted the future using divination. They probably got it wrong most of the time, which leads to the question of why people didn’t simply get rid of the priests or switch to a new religion? Why didn’t incidents like the sacred chickens being thrown overboard happen more often?

One possible answer is that it didn’t matter whether their predictions were true. It also didn’t matter whether people really believed their predictions. All that mattered was that people acted as if the predictions were true. The purpose of a long-range economic forecast is not, for example, to give an accurate guess at China’s benchmark interest rate in 2036. Everyone knows that the guess is wrong.² But lots of people making deals and plans involving China in 2036 have to agree on a number, and it might as well be the one obtained by some arcane economic modelling exercise with relatively high academic stature. In this sense, it would not be right to say that economic forecasts serve a quasi-religious role in our society. Rather, it seems that forecasting literally is a religion. It works because people have faith in it. Or rather, not because people have faith in it, but because people agree to act as if they do. If people could be convinced to have faith in a computer which spits out random numbers, that would probably work just as well (and, intriguingly, be more economically efficient).³

1: I know above I said there were no true probability models, but this is just a theoretical world.

2: Most predictive models would result in extremely wide confidence/credible intervals if uncertainty was properly taken into account. For example, when I worked in public health, models were based on “contact matrices” which had been compiled “using surveys”. There was no uncertainty in these matrices; they were just taken as-is as a model input. Even so, the modellers insisted on using 80% (or sometimes even 50%) credible intervals for their predictions because the 95% intervals were so wide that they were just silly. In practice, the end user don’t care about intervals anyway. The intervals are only for academic papers. All the public wants is a point prediction.

3: If you take this point of view, astrology begins to make a lot of sense. Why did astrology used to be such a popular form of divination? Perhaps because everybody could agree on the raw data from which the predictions were being made (namely, the positions of the stars and planets). At the same time, the predictions themselves were constructed by following arcane mathematical rules such as those described in the early editions of W. W. Rouse Ball’s Mathematical Recreations. These rules were difficult to follow, and variances in predictions from the same raw material could therefore be explained by errors in following the recipe. So astrology itself behaved as a sort of seedable random number generator which produced predictions upon which everyone could agree, which also leaving open the possibility that the predictions could be wrong. Perhaps this is why economic forecasting has to be complicated. If it wasn’t, then it would be too easy for everyone to realise that it doesn’t work?

A comment on comments: I recently removed the Disqus-powered comments section from the blog because it started injecting ads. I may replace it by an alternative such as Giscus.

Tabular foundation models are a new category of machine learning model that can perform zero-shot (i.e., without any gradient updates) prediction on tabular datasets. They are pretrained on a distribution of synthetic tables and can transfer to new datasets using in-context learning. TabICL is one of the strongest open-source tabular foundation models, and this post explores how it works and compares its performance to XGBoost on a credit-risk dataset, all from R using the reticulate package to call Python.

library(reticulate)
py_install(c("tabicl", "torch", "pandas", "scikit-learn"), pip = TRUE, pip_options = "--force-reinstall --no-cache-dir")

library(readr)
library(dplyr)
library(rsample)

# Load the credit sample (path is relative to the project root, not this post)
credit <- read.csv("https://bit.ly/42ypcnJ")

# Keep a focused subset — mostly numeric features plus a few categoricals
keep_cols <- c(
  "loan_amnt", "int_rate", "installment", "grade", "sub_grade",
  "annual_inc", "dti", "revol_util", "total_acc", "open_acc",
  "delinq_2yrs", "pub_rec", "fico_range_low", "inq_last_6mths",
  "home_ownership", "verification_status", "purpose", "term",
  "bad_flag"
)

df <- credit |>
  select(all_of(keep_cols)) |>
  mutate(across(where(is.character), as.factor)) |>
  mutate(bad_flag = factor(bad_flag))

# Train / test split (stratified on the target since it's imbalanced)
set.seed(42)
split <- initial_split(df, prop = 0.8, strata = "bad_flag")
train <- training(split)
test  <- testing(split)

# Quick class balance sanity check
prop.table(table(train$bad_flag))

        0         1 
0.8838605 0.1161395 

library(reticulate)

# Import the Python classifier
tabicl <- import("tabicl")
clf <- tabicl$TabICLClassifier(random_state = 42L)

# Keep X as a data frame — reticulate auto-converts to a pandas DataFrame,
# preserving dtypes so TabICL can detect categoricals.
X_train <- train[, setdiff(names(train), "bad_flag")]
y_train <- as.character(train$bad_flag)
X_test  <- test[,  setdiff(names(test),  "bad_flag")]

# TabICL's "fit" is cheap — it just stashes the training data.
# The actual prediction happens in the forward pass at predict() time.
clf$fit(X_train, y_train)

TabICLClassifier()

# Predicted probabilities for the positive class
proba <- clf$predict_proba(X_test)
pred_prob <- as.numeric(proba[, 2])

library(tidymodels)
library(tictoc)

# XGBoost with a modest tuning grid — a realistic "I spent 10 minutes" baseline
xgb_spec <- boost_tree(
  trees = 500,
  tree_depth = tune(),
  min_n = tune(),
  learn_rate = 0.01
) |>
  set_engine("xgboost") |>
  set_mode("classification")

xgb_grid <- grid_regular(
  tree_depth(range = c(4, 10)),
  min_n(range = c(2, 20)),
  levels = 3
)

xgb_wf <- workflow() |>
  add_model(xgb_spec) |>
  add_formula(bad_flag ~ .)

folds <- vfold_cv(train, v = 5, strata = "bad_flag")

tic("XGBoost tuning")
xgb_res <- tune_grid(
  xgb_wf,
  resamples = folds,
  grid = xgb_grid,
  metrics = metric_set(roc_auc)
)
toc(log = TRUE)

XGBoost tuning: 81.94 sec elapsed

best_xgb <- select_best(xgb_res, metric = "roc_auc")
xgb_final <- finalize_workflow(xgb_wf, best_xgb) |>
  fit(train)

xgb_prob <- predict(xgb_final, test, type = "prob")$.pred_1

tic("TabICL fit + predict")
clf$fit(X_train, y_train)

TabICLClassifier()

tabicl_prob <- as.numeric(clf$predict_proba(X_test)[, 2])
toc(log = TRUE)

TabICL fit + predict: 110.86 sec elapsed

library(pROC)

xgb_auc <- auc(test$bad_flag, xgb_prob)
tabicl_auc <- auc(test$bad_flag, tabicl_prob)

results <- tibble(
  model = c("XGBoost (tuned, 5-fold CV)", "TabICL (zero-shot)"),
  roc_auc = c(xgb_auc, tabicl_auc)
)
results

# A tibble: 2 × 2
  model                      roc_auc
  <chr>                        <dbl>
1 XGBoost (tuned, 5-fold CV)   0.648
2 TabICL (zero-shot)           0.701

The U.S. Census Bureau publishes the rent that a typical household pays in every county in the country, updated every year, and gives it away through a free API. With the tidycensus package you can pull that data — numbers and the map polygons to draw it — in a single function call, then turn it into a national map in about twenty lines. Here we map median gross rent by county to see where renting is expensive and where it is cheap.

Getting the data

The data comes from the American Community Survey (ACS), the Census Bureau’s rolling survey of about 3.5 million addresses a year. We use the 5-year estimates (here 2019–2023), which pool five years of responses so that even small, rural counties get a reliable number. The variable we want is B25064_001: median gross rent in dollars — contract rent plus utilities — for renter-occupied housing.

tidycensus talks to the ACS API, and the API needs a free key. Request one at api.census.gov/data/key_signup.html; it arrives by email in a minute. Activate it, then install it into R once — after that tidycensus finds it automatically in every session, so you never put it in a script you share:

library(tidycensus)
census_api_key("YOUR_KEY_HERE", install = TRUE)  # run once; writes to ~/.Renviron

With the key set, load the packages we need:

library(tidycensus)
library(tigris)     # shift_geometry(): repositions Alaska & Hawaii
library(dplyr)
library(sf)
library(ggplot2)
library(scales)     # dollar labels
library(patchwork)  # combine the small-multiple maps

options(tigris_use_cache = TRUE)  # cache boundary files after first download

Now one call does everything. We ask for the rent variable at geography = "county", and — this is what makes tidycensus special — set geometry = TRUE so it also downloads the county boundaries and returns them joined to the data as an sf object. Leaving state unset gives us every county in the country; shift_geometry() then tucks Alaska and Hawaii under the lower 48 so the whole nation fits one frame.

us <- get_acs(
  geography  = "county",
  variables  = "B25064_001",   # median gross rent (dollars)
  year       = 2023,
  survey     = "acs5",
  geometry   = TRUE,
  resolution = "20m"           # generalized boundaries: smaller, faster, fine for a national map
)

us <- shift_geometry(us)

The result is a tidy sf data frame — one row per county, the value in estimate, its 90% margin of error in moe, and a geometry column. Always look before drawing:

us                       # an sf object: data + polygons together
## Simple feature collection with 3222 features and 5 fields
## Geometry type: GEOMETRY
## Dimension:     XY
## Bounding box:  xmin: -3112200 ymin: -1697728 xmax: 2258154 ymax: 1558935
## Projected CRS: USA_Contiguous_Albers_Equal_Area_Conic
## First 10 features:
##    GEOID                         NAME   variable estimate moe
## 1  13027       Brooks County, Georgia B25064_001      752  79
## 2  31095   Jefferson County, Nebraska B25064_001      659  50
## 3  51683      Manassas city, Virginia B25064_001     1835  62
## 4  56021      Laramie County, Wyoming B25064_001     1080  29
## 5  13135     Gwinnett County, Georgia B25064_001     1713  16
## 6  20001         Allen County, Kansas B25064_001      685  60
## 7  27065    Kanabec County, Minnesota B25064_001     1003  83
## 8  28107   Panola County, Mississippi B25064_001      859  44
## 9  31185        York County, Nebraska B25064_001      885  49
## 10 42063 Indiana County, Pennsylvania B25064_001      786  28
##                          geometry
## 1  MULTIPOLYGON (((1163909 -64...
## 2  MULTIPOLYGON (((-115252.6 3...
## 3  MULTIPOLYGON (((1580860 292...
## 4  MULTIPOLYGON (((-765818.2 5...
## 5  MULTIPOLYGON (((1073286 -32...
## 6  MULTIPOLYGON (((41912.94 35...
## 7  MULTIPOLYGON (((192562.6 95...
## 8  MULTIPOLYGON (((527838.1 -3...
## 9  MULTIPOLYGON (((-152283 398...
## 10 MULTIPOLYGON (((1383131 460...

summary(us$estimate)     # the spread of county rents
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.     NAs 
##   253.0   742.0   848.5   928.4  1021.0  2893.0      10

sum(is.na(us$estimate))  # counties with no published estimate
## [1] 10

Two things to note. A handful of counties come back NA — mostly tiny populations the Bureau suppresses for privacy; we’ll let them render in grey. And that moe column is a reminder these are survey estimates, not a census: a county’s rent is “$1,200 ± $80,” not exactly $1,200. For a national map the uncertainty is small relative to the range, but it’s there.

Map it

Rent is a sequential quantity — low to high — so we fill with a sequential palette (viridis’s plasma), which is colorblind-safe and reads clearly in print. geom_sf() draws the polygons; everything else is labels and a clean, map-friendly theme.

ggplot(us) +
  geom_sf(aes(fill = estimate), color = NA) +
  scale_fill_viridis_c(
    option   = "plasma",
    labels   = label_dollar(),
    na.value = "grey85",
    name     = "Medianngross rent"
  ) +
  labs(
    title    = "Median gross rent by county, United States",
    subtitle = "American Community Survey, 2019–2023 5-year estimates",
    caption  = "Source: U.S. Census Bureau ACS, via the R tidycensus package"
  ) +
  theme_void(base_size = 13) +
  theme(plot.title = element_text(face = "bold"),
        legend.position = c(0.92, 0.3))

Reading the map

The typical county’s median rent is about $848 — but the map is a story of a few bright clusters against a wide, darker interior. Rent runs from $253 in Issaquena County up to $2,893 in San Mateo County, an elevenfold spread. Only about 6% of counties top $1,500, and they are not scattered randomly: they concentrate in coastal California and the Bay Area, the Washington–Boston corridor, and pockets around Seattle, Denver, and the mountain-resort West. Because every number in this paragraph is computed from the same object we mapped, the text always matches the picture — re-knit it next year and both update together.

The priciest counties make the coastal concentration concrete:

d |>
  arrange(desc(estimate)) |>
  transmute(County = NAME, `Median rent` = scales::dollar(estimate)) |>
  head(8) |>
  knitr::kable()

Zooming into the states

A national map flattens what happens inside each state. The state = argument fixes that — it accepts a vector, so one call pulls several states at once. Here we grab the four most populous, then draw each on its own colour scale to bring out where rent concentrates within each one.

states <- c("California", "Texas", "Florida", "New York")

sc <- get_acs(geography = "county", variables = "B25064_001", state = states,
              year = 2023, survey = "acs5", geometry = TRUE, resolution = "20m")
sc$state <- sub(".*, ", "", sc$NAME)   # pull the state name out of "County, State"

one_state <- function(st) {
  ggplot(filter(sc, state == st)) +
    geom_sf(aes(fill = estimate), color = "grey92", linewidth = 0.05) +
    scale_fill_viridis_c(option = "plasma", labels = label_dollar(),
                         na.value = "grey85", name = NULL) +
    labs(title = st) +
    theme_void(base_size = 12) +
    theme(plot.title = element_text(face = "bold", hjust = 0.5),
          legend.key.width = unit(0.35, "cm"), legend.text = element_text(size = 8))
}

(one_state("California") | one_state("Texas")) /
(one_state("Florida")   | one_state("New York")) +
  plot_annotation(
    title   = "Median gross rent by county — each state on its own colour scale",
    caption = "Source: U.S. Census Bureau ACS, via the R tidycensus package",
    theme   = theme(plot.title = element_text(face = "bold", size = 15)))

The pattern repeats with variations. California is expensive along almost its entire coast, with the Bay Area at the top of its range. New York is the sharpest split — downstate (New York City, Long Island, Westchester) against a uniformly inexpensive upstate. Florida’s money is on the water: the southeast around Miami, the southwest around Naples, the interior far cheaper. Texas is mostly its metros — Austin and Houston — but note the bright cluster out west, the Permian Basin oil counties, where a housing crunch has nothing to do with a big city.

One caution about reading these together: because each panel has its own scale, the same colour means different rents in different states — a “bright” Texas county (around $1,500) rents for far less than a bright California one (near (2,800). Per-state scales reveal internal *pattern*; for cross-state *levels*, use one shared scale (`limits = range(sc)estimate, na.rm = TRUE)`) instead.

Make it your own

Nothing above is specific to rent. Change one variable code and you map something else entirely; change geography and you change the resolution:

• B19013_001 — median household income
• B25077_001 — median home value
• DP02_0154PE — households with a broadband subscription (percent)
• geography = "tract" with a state and county — a neighborhood-level map of a single metro

To find a variable, browse the full ACS table with load_variables(2023, "acs5") and search it — there are tens of thousands. The pattern stays the same: get_acs() with geometry = TRUE, then geom_sf().

That’s all. You now have a script that pulls an authoritative national dataset, joins it to its own geography, and maps it — and re-points at any variable or place with a one-line change.

NOAA analysts trace wildfire smoke plumes from satellite imagery every day and publish them as open GIS shapefiles. The files live at a predictable, date-based URL, so you can pull the current day’s smoke map straight into R and plot it over the U.S. in about 30 lines. We use sf to read the spatial data and ggplot2 to draw it.

The data source

The Hazard Mapping System (HMS), run by NOAA/NESDIS, is a daily analysis where operators trace smoke over North America from GOES satellite imagery. Each plume is a polygon with four attributes: the satellite, a start and end time, and a density class — Light, Medium, or Heavy. The daily shapefile is posted here:

.../Smoke_Polygons/Shapefile/YYYY/MM/hms_smokeYYYYMMDD.zip

Because the filename is just the date, we can build the URL programmatically and always get the latest map. Note that HMS finalizes a given day’s analysis the following morning (Eastern time), so today’s file may not exist until then — the robust helper at the end handles that.

Related posts on DataScience+:

• Building Heatmaps in R with ggplot2 package
• Visualising Thefts using Heatmaps in ggplot2
• How Happy is Your Country? — Happy Planet Index Visualized

Packages

library(sf)
library(dplyr)
library(ggplot2)
library(maps)

Get the data

We format the date into the URL, download and unzip the shapefile into a temp folder, and read it with st_read(). HMS names files by calendar date, so Sys.Date() gives us today’s map.

day <- Sys.Date()
ymd <- format(day, "%Y%m%d")
url <- sprintf(
  "https://satepsanone.nesdis.noaa.gov/pub/FIRE/web/HMS/Smoke_Polygons/Shapefile/%s/%s/hms_smoke%s.zip",
  format(day, "%Y"), format(day, "%m"), ymd)

zip <- file.path(tempdir(), basename(url))
dir <- file.path(tempdir(), paste0("hms_", ymd))
download.file(url, zip, mode = "wb", quiet = TRUE)
unzip(zip, exdir = dir)

smoke <- st_read(dir, quiet = TRUE)

Inspect what we got

Always look at the data before plotting it. st_read() returns an sf object — a data frame with a geometry column — so the usual tools work.

smoke                    # note the CRS line: WGS 84 / EPSG:4326
## Simple feature collection with 98 features and 4 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -144.505 ymin: 12.23946 xmax: -11.87108 ymax: 85.46261
## Geodetic CRS:  WGS 84
## First 10 features:
##    Satellite        Start          End Density                       geometry
## 1  GOES-WEST 2026198 1200 2026198 1500   Light POLYGON ((-71.21715 19.0450...
## 2  GOES-WEST 2026198 1200 2026198 1500   Light POLYGON ((-71.99762 18.4664...
## 3  GOES-WEST 2026198 1200 2026198 1500   Light POLYGON ((-71.50895 18.6244...
## 4  GOES-WEST 2026198 1200 2026198 1500   Light POLYGON ((-76.61889 20.6635...
## 5  GOES-WEST 2026198 1200 2026198 1500   Light POLYGON ((-77.00836 21.2569...
## 6  GOES-WEST 2026198 1200 2026198 1500   Light POLYGON ((-80.49727 26.2337...
## 7  GOES-WEST 2026198 1200 2026198 1500   Light POLYGON ((-80.88491 26.4339...
## 8  GOES-WEST 2026198 1200 2026198 1500   Light POLYGON ((-82.60996 32.8358...
## 9  GOES-WEST 2026198 1200 2026198 1500   Light POLYGON ((-123.1322 42.0559...
## 10 GOES-WEST 2026198 1200 2026198 1500   Light POLYGON ((-103.8145 37.1463...

table(smoke$Density)     # how many plumes of each density today
## 
##  Heavy  Light Medium 
##     35     29     34

Two things to notice. First, the CRS is already EPSG:4326 (plain longitude/latitude), so no reprojection is needed. Second, Density takes three values — Light, Medium, Heavy — which is the variable we’ll color by. The Start/End fields use a YYYYDDD HHMM UTC format (day-of-year), but we won’t need them here.

Tidy the density classes

We coerce Density to an ordered factor and arrange() on it so heavy plumes are drawn last (on top of lighter ones), and call st_make_valid() because hand-drawn polygons occasionally self-intersect and would otherwise break the crop.

smoke <- smoke %>%
  filter(Density %in% c("Light", "Medium", "Heavy")) %>%
  mutate(Density = factor(Density, levels = c("Light", "Medium", "Heavy"))) %>%
  st_make_valid() %>%
  arrange(Density)

Map it

State outlines come from the maps package (pure CRAN, no API key). We crop the smoke to a lower-48 bounding box so the plot isn’t dominated by plumes over Canada and the oceans, then layer states underneath and smoke on top, colored by density.

states   <- st_as_sf(map("state", plot = FALSE, fill = TRUE))
bbox     <- st_bbox(c(xmin = -125, xmax = -66, ymin = 24, ymax = 50), crs = 4326)
smoke_us <- st_crop(smoke, bbox)

ggplot() +
  geom_sf(data = states, fill = "grey97", color = "grey80", linewidth = 0.2) +
  geom_sf(data = smoke_us, aes(fill = Density), color = NA, alpha = 0.6) +
  scale_fill_manual(values = c(Light = "#FFD24D", Medium = "#FB8C00", Heavy = "#C62828")) +
  coord_sf(xlim = c(-125, -66), ylim = c(24, 50), expand = FALSE) +
  labs(
    title    = "Wildfire smoke over the U.S.",
    subtitle = format(day, "NOAA HMS smoke plumes, %B %d, %Y"),
    fill     = "Smoke density",
    caption  = "Source: NOAA/NESDIS Hazard Mapping System (HMS)"
  ) +
  theme_minimal(base_size = 13) +
  theme(axis.text = element_blank(), panel.grid = element_blank(),
        plot.title = element_text(face = "bold"))

Reading the map

Today’s analysis has 98 smoke plumes over North America, 35 of them classed as heavy, and 57 intersecting the lower 48. The map turns that table into geography: where the polygons stack up and darken, smoke is thicker, and the plumes trace the path the smoke has travelled from its source fires — often hundreds of miles downwind. Re-run the code on a different day and both the numbers and the picture change; that is the point of building the URL from Sys.Date().

Because the counts and the map are generated from the same object, the paragraph above always describes the map you see — nothing is hard-coded.

A robust download

Run this early in the morning and today’s file may not be posted yet. This helper walks back day by day until it finds one that exists, so the script never dies on a missing date:

get_hms_smoke <- function(day = Sys.Date(), max_back = 3) {
  for (d in seq(0, max_back)) {
    date <- day - d
    ymd  <- format(date, "%Y%m%d")
    url  <- sprintf(
      "https://satepsanone.nesdis.noaa.gov/pub/FIRE/web/HMS/Smoke_Polygons/Shapefile/%s/%s/hms_smoke%s.zip",
      format(date, "%Y"), format(date, "%m"), ymd)
    zip <- file.path(tempdir(), basename(url))
    ok  <- tryCatch({ download.file(url, zip, mode = "wb", quiet = TRUE); TRUE },
                    error = function(e) FALSE)
    if (ok && file.exists(zip) && file.size(zip) > 0) {
      dir <- file.path(tempdir(), paste0("hms_", ymd))
      unzip(zip, exdir = dir)
      message("Using HMS smoke for ", date)
      return(st_read(dir, quiet = TRUE))
    }
  }
  stop("No HMS smoke file found in the last ", max_back, " days.")
}

Caveat: smoke aloft vs. smoke you breathe

HMS is a satellite product — it maps smoke seen from above, at any altitude. A thick plume on the map can sit five kilometers up and leave ground-level air clear, so HMS is a measure of smoke transport, not of what people are breathing. For surface air quality, pair it with ground PM2.5 from the EPA/USFS AirNow Fire and Smoke Map, which also has an API — a natural follow-up analysis.

That’s all. You now have a script that pulls a hand-analyzed satellite product and turns it into a national map, and re-runs itself for any day you like.

This article was first published on DataScience+, a community of R and Python tutorial authors. Have a data-science technique worth sharing? Write for us — no pitch required.

So you’ve built your data pipeline, you’ve designed your dashboard and you’ve connected the two through a real app. The screenshots match the designs and you can show it off and wow the stakeholders. Ready to launch, right? Well, maybe. In this post we’ll cover five things that are easy to overlook: some worth a quick check, some you might not have thought to add.

Sample Dashboard

Before we dive into those five things, however, I’ll introduce the sample dashboard. It uses data from the {nycflights13} R package developed by Hadley Wickham but is otherwise built entirely without R. Instead it’s built using standard frontend web technologies: HTML, CSS and JavaScript, with data coming from a single monolithic gzipped JSON file. The JavaScript then adds some noise to the data, so it looks more like live data (and so the dashboard doesn’t appear blessed with 20/20 foresight).

The actual visual presentation is that of online departure boards for the three New York/New Jersey airports — John F. Kennedy Intl, La Guardia and Newark Liberty Intl — along with some weather information and some assumed KPIs.

Each table row and each KPI card is clickable, bringing up a modal with more detailed information. Please note that, spoiler alert, the dashboard does not (yet) function well on mobile. We will get to that shortly.

Full disclosure: for rapid prototyping, the dashboard was built through an extended prompting session with Claude Code. Most of the issues I’ll cover are genuine problems that came up in the session, but in one case I did ask Claude to leave out a feature it clearly wanted to add, purely for illustrative purposes. This example, in the state it’s in, is meant to represent a dashboard that an intelligence — human, artificial or Martian — might reasonably judge ready to ship. The rest of this post exposes some remaining flaws.

Five Things to Check

Does the design work at different screen dimensions?

This is the easiest one to check and, perhaps, the most important. While dashboards are still primarily built on desktop machines, making them viewable on tablets and mobile phones should be a part of most development plans. Many, probably most (data are inconsistent and highly topic- and region-dependent), humans use mobile devices to consume content on the World Wide Web, and you should have a very good reason for excluding them from your dashboard. It’s not just improvements in mobile hardware pushing people that way; the expansion of fifth-generation mobile network technology (5G) can make the experience much nicer than it used to be, especially when the consumers themselves are mobile.

This screenshot of our sample dashboard quickly highlights a big problem with it: it is literally impossible to bring up the data for La Guardia airport because its “button” (it’s not really a button, which is another issue we’ll come on to later) is off the right side of the page, with no possibility of scrolling to it.

Less critical, but still worth fixing, is the two-way scrolling of the table. This is avoidable by stacking table columns on narrow screens. Both of these fixes are shown in the next screenshot. (The content is much taller now so we show the full page, not just the visible part. On my phone, for example, the top of the table can just be seen without scrolling.)

There’s actually a problem with the layout in both dimensions: on a shorter screen — perhaps a small browser window on a cramped notebook — the table can become unreachable. The problem here is that, in designing for (larger) desktop browsers, the design assumed a minimum browser height and then scrolling only on the table and not on the browser window itself. This works well when the browser window is tall enough, since you get fixed information at the top with no nested scrolling. But when this (arbitrary) minimum height condition isn’t met, the design fails. So that’s another thing that should be fixed prior to launch.

Does the app support touchscreens?

Supporting mobiles and tablets isn’t just about layout. Interactions need to work when the input comes from a finger or stylus rather than a mouse and keyboard. There are actually two things to consider here:

1. Can all interactions be done through touch?
2. Is it obvious what interactions are actually possible?

The base app allows the user to “change time” through a dialog that only opens when pressing the T key. This, obviously, fails 1). We can fix it by making the clock itself into a tappable object that brings up the same dialog as pressing the T key (this also gives the user access to the pause option that was previously only possible by pressing the P key). This brings us neatly on to 2): it’s not inherently obvious that the clock is tappable. Nor is it obvious that the cards and rows can be tapped for detailed information. With a mouse at least, these change colour (slightly) when hovered over. We need something more: we can go with a common icon and some explanatory text.

Is the app accessible?

Beyond use on mobile, there are a number of additional accessibility issues with the base version of the app. At least three of these relate to keyboard usage:

1. Button-like interactive elements are not HTML button elements and therefore don’t have appropriate interactions and semantics to meet WCAG requirements;
2. Not everything that can be clicked/tapped can be interacted with by a keyboard user;
3. The modal that appears when selecting interactive elements does not implement focus trapping, meaning a keyboard user can lose focus somewhere in the background behind the modal window.

We can fix all these issues with better use of HTML elements: buttons for airport selection, KPI-modal launch, and flight modal launch, and the dialog element for the modal windows. On top of this, some text in the base app does not meet colour-contrast requirements, so we also want to fix this.

The keyboard focus indicator uses the yellow accent colour found elsewhere in the app. In hindsight this was probably a mistake. A different colour focus indicator would make it clearer what had keyboard focus as opposed to what had actually been selected (i.e. which airport was currently being shown).

Is the app performant?

The app is meant to simulate real data but allows the user to select, pause and speed up time. Consequently, the data isn’t passed around in the format it would be if it were a real app tracking live events. But we can still try to optimise how we pass around data in our simulated-data app. The biggest win we found here was converting the data from row-based to column-based. This reduced the size of the data (gzipped) by a third, from 4.8 MB to 3.2 MB, and JSON parsing times from ~550 ms to ~250 ms.

It is, of course, important not to waste time and resources on premature and unnecessary optimisations. So consider the data change above as an example of where one might want to look rather than a change that was strictly necessary.

Another place to look is at the size of images. For dashboards that often means data visualisation. Large, frequently changing raster images sent from the server can lead to performance hits. This is what you get, by default, in Shiny apps when you use renderPlot()/plotOutput(). Packages that render SVG or HTML canvas directly in the browser can greatly reduce the overhead.

Does the app feel like it responds immediately to interaction?

Actual performance is one thing, but perceived performance can also be significant. If a user clicks on a button and nothing happens for a second or two, things “feel” broken even if, in reality, the browser is just processing things silently in the background. In our particular case, there is a simulated delay of ~2 seconds when changing airport. Click or tap on an airport and nothing happens. (Disclosure: this is where I actually had to force Claude not to put something in in order to make my point.) There are a number of options here for how to deal with this. The simplest might be to set the cursor to a value of “wait” in CSS, but this would obviously not help users of touchscreen devices that lack a cursor. A better alternative is to show a loading indicator to everyone. These typically rotate, so we could use an animated GIF or a static image file and CSS animations.

For longer waits, a progress indicator would be a better option if it were possible to actually estimate progress. This would ensure it didn’t look like the page had got stuck. For a wait of a second or two, this is not necessary; the requirement is only that the user sees that their action of clicking on an airport had an effect immediately.

Summary

You can see all the fixes discussed in this post (and some more we had to skip over for brevity) in the final version of the dashboard. Hopefully this article has illustrated that there can be quite a few issues with a newly designed dashboard that may be subtle or hidden, particularly if you’re only testing with a mouse on a desktop computer. The goal here wasn’t to scare you away from publishing entirely, so hopefully I haven’t done that, but to highlight that the devices via which a dashboard is developed and consumed are often quite different: remember to design for the latter. Similarly, just because you know something works and will load eventually does not mean you can expect your users to be as patient and forgiving.

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

      The Consumer Price Index (CPI) is a widely used measure of the prices of goods and services purchased by households. It’s the primary tool for tracking inflation and changes in the cost of living over time. The index is built from monthly price collections on a “basket” of goods and services from a sample of retail and service establishments. Historical CPI data is easy to download

    &nbsp A key feature of the CPI is that prices are adjusted for quality changes. If the price of a car’s side mirror rises by $200, but $120 of that increase reflects the mirror becoming “smart” rather than “dumb,” only the remaining $80 is counted as inflation. Similarly, if the price of a medical procedure rises because new equipment improves the quality of care, the portion attributable to improved quality is removed from the inflation calculation. Consumers still pay for these quality improvements, whether they want them or not, so in many cases the CPI understates pure cost increases.

    &nbsp There are eight major categories of the CPI, and each category has its own index: Food & Beverages, Housing, Apparel, Transportation, Medical Care, Recreation, Education & Communication, and Other. These are weighted to form the overall CPI, with the largest weights as Housing at about 44% of the total, Transportation at 17%, Food & Beverages at 14%, and Medical Care at 8%. (Each of these is further sub-divided into its own index; for example, Other includes Personal Care, and Personal Care has seprate indices for Cosmetics, Perfume, Bath, and Nail Preparations.) Of course the weights will not reflect your percentages of what you buy.

    &nbsp Downloading historical CPI data from FRED (Federal Reserve Bank of St. Louis) was easier than I expected. You need an API key which you can get from https://fredaccount.stlouisfed.org/login

    &nbsp Here is a line graph showing cumulative CPI growth for the overall CPI and each of the eight major categories through June 2026, indexed to December 2016 = 1.00. The overall CPI has risen 37.1% since December 2016. Housing (rent, insurance, energy, etc.) has increased the most at 44.8%. Transportation (vehicle purchases, fuel, maintenance, insurance, public transit fares) is next at 43.5%. Food (groceries and restaurants) is up 39.9%. Medical Care is lower at 25.9%. Medical Care includes out‑of‑pocket spending on providers, hospitals, and insurance, but excludes employer‑paid and government‑paid health insurance premiums.

      The CPIs exclude a lot of things like the quality changes I mentioned above, and other items that you may pay but that the government does not classify as personal expenses. We all feel the cost of our groceries going up – I discussed this previously in groceries . Each category has its own reasons why it is increasing; I will leave that discussion to the economists, except to say that the cost of energy affects a lot of items in the cost of production and delivery.

      Like many broad measures, the CPI is an attempt to estimate the overall cost of goods and services. But what ultimately matters to you is the actual cost of the things you buy.

      Here is my R code:

library(tidyquant)
library(dplyr)
library(tidyr)
library(lubridate)
library(ggplot2)
library(ggrepel)   # repel overlapping text labels

# Set your API environment variable
Sys.setenv(FRED_API_KEY = "xxxx")


# Define the official FRED database tracking codes
cpi_series <- c(
  "CPIAUCSL", "CPIFABSL", "CPIHOSSL", "CPIAPPSL", 
  "CPITRNSL", "CPIMEDSL", "CPIRECNS", "CPIEDUNS", "CPIOGSNS"
)

# Download and clean data vectors
raw_data <- tq_get(cpi_series, get = "economic.data")   # get from FRED

# 2. Clean, Filter, and Perform Group-Indexing
cpi_processed <- raw_data %>%
  mutate(
    Year  = year(date),
    Month = month(date)
  ) %>%
  # Keep all Decembers from 2015 onward OR strictly isolate June 2026
  filter((Month == 12 & Year >= 2015) | (Year == 2026 & Month == 6)) %>%
  # Convert raw tracking codes into readable titles
  mutate(Category = case_when(
    symbol == "CPIAUCSL" ~ "Overall CPI",
    symbol == "CPIFABSL" ~ "1. Food & Bev",
    symbol == "CPIHOSSL" ~ "2. Housing",
    symbol == "CPIAPPSL" ~ "3. Apparel",
    symbol == "CPITRNSL" ~ "4. Transportation",
    symbol == "CPIMEDSL" ~ "5. Medical Care",
    symbol == "CPIRECNS" ~ "6. Recreation",
    symbol == "CPIEDUNS" ~ "7. Education & Comm",
    symbol == "CPIOGSNS" ~ "8. Other Goods"
  )) %>%
  # Chronologically sort each group, then anchor base-100 to the first row (Dec 2015)
  group_by(Category) %>%
  arrange(date, .by_group = TRUE) %>%
  mutate(Indexed_Value = (price / first(price)) * 100) %>% 
  ungroup() %>%
  # Convert Timeline to ordered categories for a clean discrete X-Axis
  mutate(Period = if_else(Month == 6, paste0(Year, " (June)"), as.character(Year))) %>%
  mutate(Period = factor(Period, levels = unique(Period[order(date)])))

# 3. Isolate final data point rows for the text tags
label_data <- cpi_processed %>%
  group_by(Category) %>%
  filter(date == max(date)) %>%
  ungroup()

# 4. Generate the Chart with the Categorical String Axis Baseline
ggplot(cpi_processed, aes(x = Period, y = Indexed_Value, color = Category, group = Category)) +
  geom_line(aes(linewidth = ifelse(Category == "Overall CPI", 1.5, 0.8))) +
  geom_point(size = 2) +
  
  # Non-overlapping direct text labels
  geom_text_repel(
    data = label_data,
    aes(label = paste0(Category, " (", round(Indexed_Value, 1), ")")),
    nudge_x = 0.5,             
    direction = "y",          
    hjust = 0,                
    segment.color = "grey50", 
    segment.size = 0.4,
    force = 2,
    fontface = "bold",
    size = 4   # millimters
  ) +
  
  # High-contrast visual color mapping matrix
  scale_color_manual(values = c(
    "Overall CPI"         = "#000000", 
    "1. Food & Bev"       = "#E64B35", 
    "2. Housing"          = "#56B4E9", 
    "3. Apparel"          = "#009E73", 
    "4. Transportation"   = "#4D8805", 
    "5. Medical Care"     = "#0072B2", 
    "6. Recreation"       = "#D55E00", 
    "7. Education & Comm" = "#CC79A7", 
    "8. Other Goods"      = "#999999"  
  )) +
  
  # Format plot margins to prevent label cropping
  scale_x_discrete(expand = expansion(mult = c(0.05, 0.35))) +
  scale_linewidth_identity() + 
  labs(
    title = "10.5-Year Cumulative CPI Growth Comparison",
    subtitle = "Base Index: December 2016 = 100",
    x = "Reporting Period",
    y = "Indexed Value (Relative to 100)"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    legend.position = "none",          
    panel.grid.minor = element_blank(),
    text = element_text(face = "bold"),
    plot.title = element_text(face = "bold", size = 14),
    axis.text = element_text(face = "bold")
  )

# 5. Save the final graphic output file
# ggsave("cpi_growth_comparison.png", width = 12, height = 7, dpi = 300, bg = "white")

# Extract and print the final 10.5-year cumulative values
final_column_summary <- cpi_processed %>%
  filter(date == max(date)) %>%
  select(Category, Indexed_Value) %>%
  mutate(Indexed_Value = round(Indexed_Value, 2)) %>%
  arrange(desc(Indexed_Value)) # Sorts from highest inflation to lowest

print(as.data.frame(final_column_summary))

End

Reposted from the original at https://blog.stephenturner.us/p/turn-new-data-into-quarto-reports-automatically

The {watcher} R package monitors your filesystem and run arbitrary code when files change. You can use this to automate things like creating parameterized Quarto reports. 

—

The watcher R package (watcher.r-lib.org) monitors your filesystem for changes, and can run code automatically when data is created or updated.

A helpful use case for this is to monitor a folder for changes, then render a Quarto report for whatever new data arrived.

Simple example here starting with an empty data directory and a Quarto template.

$ tree
.
├── data
└── report.qmd

This is a parameterized Quarto template that uses Typst to compile a simple PDF report showing a summary() of a CSV you read in.

---
title: "Automatically compiled report"
author: "Stephen Turner"
subtitle: "File: `r basename(params$csv_path)`"
date: today
format: typst
params:
    csv_path: NA
---

Compiled `r format(Sys.time(), "%Y-%m-%d %H:%M:%S %Z")` from `r params$csv_path`.

```{r}

```

```{r}
df <- read.csv(params$csv_path)
summary(df)
```

Now let’s set up the watcher. The watcher monitors a directory for new files, then runs quarto_render passing in the new file as a parameter.1

library(watcher)

render <- function(paths) {
  message(format(Sys.time()), ": ", length(paths), " file(s) changed")
  message(paths)
  quarto::quarto_render(
    "report.qmd",
    output_file = basename(paths),
    execute_params = list(csv_path = paths),
    quiet = TRUE
  )
}

w <- watcher(path = "data", callback = render, latency = 1)
w$start()

Now whenever new files land in data/ the watcher will automatically render the parameterized Quarto document, which just prints a summary of the data. The w$start() doesn’t tie up my R console. It’s running in the background.

Now, when I create new CSV files in the data directory, the watcher finds these and renders the reports.

> iris |> write.csv("data/iris.csv")
2026-07-15 05:45:28: 1 file(s) changed
/Users/sdt5z/Downloads/watcher-test/data/iris.csv

> penguins |> write.csv("data/penguins.csv")
2026-07-15 05:45:33: 1 file(s) changed
/Users/sdt5z/Downloads/watcher-test/data/penguins.csv

With this you could start the watcher in a background R process (e.g., running under tmux or something), monitoring a shared drive or some cloud location. Whenever new data arrives, a report gets compiled without you having to do anything.

More on the watcher package:

• Posit blog: https://opensource.posit.co/blog/2026-06-29_watcher-0-2-0/

• watcher docs: https://watcher.r-lib.org/

This is a short post to announce the release of an R package, VolcanoPlotR.

• Code
• Documentation

Background

Using proteomics, we often want to compare two experimental groups. A popular way to visualise this comparison this is via a volcano plot, where the enrichment of proteins in one condition is towards the right and their de-enrichment (or their enrichment in the other group) is towards the left. The y-axis denotes the p-value of the comparison. Years ago, I wrote an IGOR package called VolcanoPlot, that we have used in our work. For various reasons (see below), I have now ported the package to R.

Features

The current version takes an output (or multiple outputs) from MaxQuant, processes the data and can generate:

• Volcano Plots
• Mean vs mean plots
• PCA
• Text outputs of the analysis

It uses {ggplot2} for graphics and therefore they can be combined easily with other plots.

Some examples

There is full documentation available here. But in brief:

The package can be installed in the usual way.

# install.packages("pak")
pak::pak("quantixed/VolcanoPlotR")

There is an example file in the package that can be used for testing.

library(VolcanoPlotR)
# get the path to the proteinGroups.txt file included in the package
filepath <- system.file("extdata", "proteinGroups.txt", package = "VolcanoPlotR")
# get the filename fromt the path
filename <- basename(filepath)
# get the directory name
filedir <- dirname(filepath)
# run the automated procedure we will also tell it which groups to compare so
# we don't have to select them interactively
workflow_maxquant(file = filename, datadir = filedir,
                  group1 = "WT", group2 = "Control")

Which gives the following volcano plot:

This plot can be styled in several ways by the package, but also restyled using ggplot.

Rather than this automated workflow, the data can be loaded in, processed and used in other functions:

df <- load_maxquant(file = filename, datadir = filedir)
df <- process_maxquant(df, group1 = "WT", group2 = "Control")
# same output as the automated workflow
volcano_plot_maxquant(df)
# mean-mean plot with no labelling
mean_plot_maxquant(df)
# PCA plots to show the structure of the data
pca_plot_maxquant(df)
pca_plot_maxquant(df, by_protein = TRUE)

Why write an R package for volcano plots?

The best tool for analysing MaxQuant data is the Windows program Perseus. As wonderful as Perseus is, its graphics are terrible. This prompted me to write VolcanoPlot for IGOR years ago. The other motivation was that I wanted to make sure the analysis was correct (I felt back then that Perseus was too “black box” but I no longer hold that view). R is the main language we use in the lab these days and because IGOR is closed source and is no longer developed for mac, I am ported VolcanoPlot to R.

I wrote a while back about how I could recreate the look of VolcanoPlot in R, and the logical next step was to make a package, so that we no longer relied on VolcanoPlot. Potentially VolcanoPlotR could be expanded into a full-blown Perseus substitute, especially since PerseusR has been abandoned.

It’s true that there are several R packages out there, so why write anothere. There’s {EnhancedVolcano} and proteomics data analysis available in Bioconductor, as well as Joachim Goedhart’s shiny app. None of these suit our needs, for example VolcanoPlotR will combine multiple MaxQuant datasets, and I was also keen to keep the look of our graphics the same for consistency between papers.

—

The post title comes from “Eruption” by Van Halen from the Van Halen album. It’s an Eddie Van Halen instrumental which has broken the spirit of many a young guitar-slinger who’s tried to play it.

Getting Involved with rOpenSci

I first came to rOpenSci in 2022, though at the time I barely knew what it was. I was getting a statistical package of mine ready to submit to the Journal of Statistical Software, and that is how I was pointed toward rOpenSci review: the journal directs authors to rOpenSci’s statistical software standards, so going through the review looked like a convenient step along the way. At the time, my focus was on polishing the software for the journal submission, not on rOpenSci itself.

What I did not expect was how much the process itself would change my perspective. Coming from an academic background, the open review on GitHub felt very different from the closed, anonymous process I was used to. It was rigorous without being adversarial. An editor and two reviewers carefully examined the package, and their feedback was constructive and grounded in rOpenSci’s well-defined standards and guidelines. The review improved the package, and it also gave me a greater appreciation for the collaborative approach behind open-source software review.

I was later invited to review a package myself. The experience gave me a different perspective on the review process. As a reviewer, I saw that the goal was not simply to determine whether a package met some standards or merits, but to help authors improve their software through constructive feedback. When applications opened for the Champions Program, mentoring felt like a natural next step. Having experienced rOpenSci as both a software author and a reviewer, it seemed like a meaningful way to contribute to the community.

Mentoring in the Champions Program

I was matched with Sunny Tseng as her mentor. Over the program, Sunny built bbsTaiwan, an R package that makes Taiwan’s Breeding Bird Survey data much easier to access and analyze. It was a real package solving a real problem for people who study Taiwan’s birds, which made it a pleasure to work on together.

Mostly, what I gave was time and attention. We worked through package scope, unit testing, version control, and the other practical aspects of building an R package. Many of our conversations were not about solving a particular technical problem, but about discussing trade-offs, identifying useful resources, and thinking through the next steps. Those conversations ended up being one of my favorite parts of the program.

What I valued most, though, was seeing how Sunny’s work was used after the project. It is easy to think of a package as code made available for others to use, but in this case it became a tool that supported people working with the same data. While visiting Taiwan, she also ran a session introducing it to members of that community. This reinforced my view that open-source software is not just code shared in public, but a way of bringing people together around shared work.

Looking Back

More than anything, my time with rOpenSci has left me with an appreciation for how thoughtfully it is run. Its initiatives, from peer review to mentorship, are organized with real care, and they are built to do more than improve software. They are designed to connect people, bring contributors together, and keep community at the center of open science.

Fifteen years in, what strikes me about rOpenSci is that it has always been more about people than about packages. What began as a convenient step on the way to a journal turned into one of the more rewarding parts of my work, and I have found value in each perspective I have seen as an author, a reviewer, and a mentor.

This post is mainly (but not only) a test, because I had a broken xml feed for my R posts, and I wanted to see if it was fixed.

It’s about my last R posts from june and july, which are:

• Using scikit-learn models in R easily with the tisthemachinelearner R package
• How conformalization helps weak models
• Fast Conformal Prediction for Some Machine Learning Models (jackknife+ and no refitting)

Using scikit-learn models in R easily with the tisthemachinelearner R package

This post is about the tisthemachinelearner R package, that allows to use scikit-learn models in R. It is a wrapper around the tisthemachinelearner Python package. Prediction intervals can be computed using either split conformal prediction, surrogate methods or the bootstrap.

Read: https://thierrymoudiki.github.io/blog/2026/06/21/r/tisthemllearner

How conformalization helps weak models

In this post, we compare split conformal prediction across several predictive models, using R package mlS3.

Read: https://thierrymoudiki.github.io/blog/2026/06/07/r/conformalization-helps-weak-models

Fast Conformal Prediction for Some Machine Learning Models (jackknife+ and no refitting)

It’s surprisingly fast to obtain conformal jackknife+ prediction intervals for Machine Learning models of the form \(\hat{y} = Sy\) (including Ordinary Least Squares, Ridge Regression, Random Vector Functional Link Networks, Kernel Ridge Regression, smoothing splines, and local polynomial regression). No refitting involved, just Linear Algebra. Read https://www.researchgate.net/publication/408161842_Fast_Conformal_Prediction_for_Some_Machine_Learning_Models_via_Closed-Form_Jackknife.

All of this is a very early work.

I have been using Shiny Server for around 8 years and R for almost 11 years.

At some point I created the tabler R package, whose current CRAN version is a layer between Shiny and Tabler to offer a refreshed look for dashboards.

The experimental tabler version on GitHub completely drops Shiny usage on its codebase and it is also compatible with widgets created for Shiny such as d3po.

Besides local apps, I have been working on Tabler Server to manage tabler apps on a Linux machine.

Nothing of this is particularly novel. Most of this was done standing on the shoulders of giants to combine R and Linux tools efficiently.

Here is a simple example:

pak::pkg_install("pachadotdev/tabler")

library(tabler)

ui <- page(
  title = "Example 1",
  layout = "boxed",
  navbar = list(
    top = topbar(title = "Example 1")
  ),
  body = body(
    h2("Example 1"),
    p("Type your name and it is echoed back below."),
    textInput("name", "Your name", value = "world"),
    textOutput("greeting")
  )
)

server <- function(input, output, session) {
  output$greeting <- renderText({
    paste0("Hello, ", input$name, "!")
  })

  syncUrl(session) # syncUrl(session, exclude = c("exclude", "this"))
}

tablerApp(ui, server)

A highlight for this project is that it that it simplifies URLs. Moving the sliders, dropdowns or text inputs updates the URL and changing the URL also modified the state app. For example, http://127.0.0.1:3000/?name=John can be changed to http://127.0.0.1:3000/?name=George, and then the text entry can be changed to “Paul” to print “Hello, Paul!” to then set it to “Ringo” or another value. Unlike Shiny, the URL does not use quotes and allows for cleaner inputs like ?year=2000&country=gbr&lang=en instead of ?year=2000&country=%22gbr%22&lang=%22en%22.

Both tabler and tabler server are released under the Apache License 2.0, a permissive license for commercial and non-commercial projects.

I hope you like this. Please feel free to test tabler, open issues or contribute to the codebase. If you find this useful, please consider donating on Buy Me A Coffee.