Category: R

grouper: An R package for Optimal Group Assignment

Introduction

Universities are increasingly using collaborative learning pedagogies, which can benefit learners through deeper understanding of course content and teamwork skills. However, the realisation of these sought-after benefits depend on how educators assign learners to groups.

Educators have formulated various mathematical models to perform this assignment. Some have developed developed models that prioritised maximising students’ project preferences. Others developed a model that prioritised students’ preferences, group sizes and group composition. Yet other models address related, but distinct, problems such as assigning students to elective courses or incorporating staff workload into student-to-project supervisor assignments.

Whichever approach is used, it is apparent that there is a need for an algorithmic solution for the assignment. This would ease the burden on the instructor, while providing an objective procedure for the assignment. Our contribution is an R package grouper that offers two flexible group allocation strategies.

Optimisation Models

grouper provides two distinct integer linear programming optimisation models.

library(grouper)
library(ompr)
library(ompr.roi)
library(ROI.plugin.glpk)

Preference-Based Assignment

The Preference-Based Assignment (PBA) model allows educators to assign student groups to topics to maximise overall student preferences for those topics. The topics can be viewed as project titles. The model allows for repetitions of each project title. This formulation also allows each project team to comprise multiple sub-groups. This is useful in cases where the project requires teams with different functionality to work together, e.g. where one team works on a front-end while the other develops a back-end model.

To execute the optimisation routine, an instructor prepares:

      1. A group composition table listing the member students within each self-formed group
      2. A preference matrix containing the preference that each self-formed group has for each topic.
      3. A YAML file defining the remaining parameters of the model.

Examples

Consider the following simple dataset with 8 students: 
pba_gc_ex002
#>   id grouping
#> 1  1        1
#> 2  2        1
#> 3  3        2
#> 4  4        2
#> 5  5        3
#> 6  6        3
#> 7  7        4
#> 8  8        4

Each student is in a self-formed group of size 2, indicated via the grouping column. Suppose that, for this set of students, the instructor wishes to assign students into two topics, with each topic having two sub-groups. This requires the preference matrix to have 4 columns – one for each topic-subgroup combination. Remember that the ordering of topics/subtopics in the preference matrix should be:

Topic1-Subtopic1, Topic2-Subtopic1, Topic1-Subtopic2, Topic2-Subtopic2

Thus there should be 4 rows in the preference matrix – one for each self-formed group.

pba_prefmat_ex002
#>      col1 col2 col3 col4
#> [1,]    4    3    2    1
#> [2,]    3    4    2    1
#> [3,]    1    2    4    3
#> [4,]    1    2    3    4

The YAML file for this model contains the following parameters:

n_topics: 2
B: 2
R: 1
nmin: 2
nmax: 2
rmin: 1
rmax: 1

B corresponds to the number of sub-topics per topic, while rmin and rmax denote the minimum and maximum number of repetitions of each topic. nmin and nmax denote the minimum and maximum number of members in each sub-topic group.

It is possible to assign each self-formed group to its optimal choice of topic-subtopic combination. In our solution, we should see that group 1 is assigned to subtopic 1 of topic 1, group 2 is assigned to sub-topic 1 of topic 2, and so on.

df_ex002_list <- extract_student_info(pba_gc_ex002, "preference", 
                                     self_formed_groups = 2, 
                                     pref_mat = pba_prefmat_ex002)
yaml_ex002_list <- extract_params_yaml(system.file("extdata", 
                                         "pba_params_ex002.yml",  
                                          package = "grouper"),
                                      "preference")
m2 <- prepare_model(df_ex002_list, yaml_ex002_list, "preference")
result2 <- solve_model(m2, with_ROI(solver="glpk"))

assign_groups(result2, assignment = "preference", 
              dframe=pba_gc_ex002, yaml_ex002_list, 
              group_names="grouping")
#>   topic2 subtopic rep group size
#> 1      1        1   1     1    2
#> 2      2        1   1     2    2
#> 3      1        2   1     3    2
#> 4      2        2   1     4    2

Diversity-Based Assignment

The Diversity-Based Assignment (DBA) model enables educators to assign students to groups and topics with the dual, but weighted, aims of maximising diversity (based on student attributes) within groups and balancing specific skill levels across different groups.

To execute the DBA optimisation routine, the instructor prepares:

      1. A group composition table containing:
        1. the member students within each self-formed group,
        2. the demographics that will be used to compute pairwise dissimilarity between students, and
        3. a numeric measure of each student’s skill.
      2. A YAML file defining the remaining parameters of the model.

Examples

Consider the following dataset, that comes with the package. There are 4 students in total. 
dba_gc_ex001
#>   id major skill groups
#> 1  1     A     1      1
#> 2  2     A     1      2
#> 3  3     B     3      3
#> 4  4     B     3      4

It is intuitive that an assignment into two groups of size two, based on the diversity of majors alone, should assign students 1 and 2 into the first group and the remaining two students into another group.

The corresponding YAML dba_gc_ex001.yml file for this exercise consists of the following lines:

n_topics:  2
R:  1
nmin: 2
nmax: 2
rmin: 1
rmax: 1

To run the assignment, we can use the following commands. We can use either the gurobi solver, or the glpk solver for this example. Both are equally fast.

# Indicate appropriate columns using integer ids.
df_ex001_list <- extract_student_info(dba_gc_ex001, "diversity",
                                      demographic_cols = 2, 
                                      skills = 3, 
                                      self_formed_groups = 4)

yaml_ex001_list <- extract_params_yaml(system.file("extdata", 
                                         "dba_params_ex001.yml",  
                                         package = "grouper"),
                                       "diversity")
m1 <- prepare_model(df_ex001_list, yaml_ex001_list,
                    assignment="diversity",w1=0.5, w2=0.5)

result3 <- solve_model(m1, with_ROI(solver="glpk"))
assign_groups(result3, assignment = "diversity", 
              dframe=dba_gc_ex001, 
              group_names="groups")
#>   topic rep group id major skill
#> 1     1   1     2  2     A     1
#> 2     1   1     3  3     B     3
#> 3     2   1     1  1     A     1
#> 4     2   1     4  4     B     3

We can see that students 2 and 3 have been assigned to topic 1, repetition 1. Students 1 and 4 have been assigned to topic 2, repetition 1. w1 and w2 both have weights 0.5, which means the skills and demographic inputs are given equal weight in the optimisation.

At present, the routines use the daisy function from the cluster package to compute a pairwise dissimilarity matrix between students. However, it is also possible to supply your own custom dissimilarity matrix. Consider the following dataset of 4 students:

dba_gc_ex003
#>   year   major self_groups id
#> 1    1    math           1  1
#> 2    2 history           2  2
#> 3    3    dsds           3  3
#> 4    4    elts           4  4

Now consider a situation where we wish to consider years 1 and 2 different from years 3 and 4, and math and dsds (STEM majors) to be different from elts and history (non-STEM majors). For each difference, we assign a score of 1. This means that students 1 and 2 would have a dissimilarity score of 1 due to their difference in majors. Students 1 and 3 would also have a score of 1, but due to their difference in years. Students 1 and 4 would have score of 2, due to their differences in majors and in years. The overall dissimilarity matrix would be:

d_mat <- matrix(c(0, 1, 1, 2,
                  1, 0, 2, 1,
                  1, 2, 0, 1,
                  2, 1, 1, 0), nrow=4, byrow = TRUE)

To run the optimisation for this model, we can execute the following code:

df_ex003_list <- extract_student_info(dba_gc_ex003, "diversity",
                                       skills = NULL,
                                       self_formed_groups = 3,
                                       d_mat=d_mat)
yaml_ex003_list <- extract_params_yaml(system.file("extdata",   
                                         "dba_params_ex003.yml",
                                         package = "grouper"), 
                                       "diversity")
m3 <- prepare_model(df_ex003_list, yaml_ex003_list, w1=1.0, w2=0.0)
result <- solve_model(m3, with_ROI(solver="glpk")

assign_groups(result, "diversity", dba_gc_ex003,
              group_names="self_groups")
#>   topic rep group year   major id
#> 1     1   1     1    1    math  1
#> 2     1   1     4    4    elts  4
#> 3     2   1     2    2 history  2
#> 4     2   1     3    3    dsds  3

As you can see, the members of the two groups have maximal difference between them – they differ in terms of their year, and in terms of their major. Notice that we specified

skills = NULL

and

w2 = 0.0

This ensures that no skills columns were taken into account in this optimisation.

Gurobi Optimiser

While the routines above use the glpk optimiser, we recommend using the Gurobi optimiser. The latter is a commercial software that runs to completion much faster than glpk. For more information, please refer to this website. Note that academic licenses are available from Gurobi.

Shiny Applications

The package provides numerous options for each of the two optimisation models. However, there are also two shiny applications included with the package. They may be useful if one only needs a straightforward group assignment. 

To run the DBA shiny app, the following code will suffice: 
library(shiny)
runApp(appDir=system.file("shiny", "dbaWebApp", package="grouper"))

# Analogous code for PBA app:
# runApp(appDir=system.file("shiny", "pbaWebApp", package="grouper"))

Here is a screen shot of the diversity-based shiny application.



The system folders with the shiny apps also contain example csv files for use with the apps.

More Details

The two optimisation models are flexibly parametrised. Here are some of the features:
    • Define the number of repetitions for each topic.
    • Define the max. and min. number of group members for each topic.
The vignettes also contain the precise mathematical formulation of the optimisation models. For full details, please refer to these links:

Understanding R’s `describe()` Function: A Complete Guide to Summary Statistics

Understanding R’s describe() Function: A Complete Guide to Summary Statistics

Table of Contents

  1. Introduction to describe()
  2. Breaking Down the Output Columns
  3. Key Statistics and Their Interpretation
  4. Practical Examples
  5. When to Use Which Statistic
  6. Extending the Functionality
  7. Conclusion

Introduction to describe()

The describe() function from R’s psych package (Revelle, 2023) provides a comprehensive statistical summary of your dataset. Unlike R’s base summary() function, it includes additional metrics that are particularly useful for data exploration and assumption checking.

library(psych)
describe(your_data)

Breaking Down the Output Columns

Here’s what each column in the output represents:

Column Description Formula/Calculation Ideal Use Case
vars Variable index number – Tracking variable order
n Complete cases length(na.omit(x)) Data completeness check
mean Arithmetic average sum(x)/n Normally distributed data
sd Standard deviation sqrt(var(x)) Measuring spread
median 50th percentile quantile(x, 0.5) Skewed distributions
trimmed Mean after removing extremes mean(x, trim=0.1) Robust central tendency
mad Median absolute deviation median(abs(x-median(x))) Outlier-resistant spread
min Minimum value min(x) Range assessment
max Maximum value max(x) Range assessment
range Max – Min max(x)-min(x) Total spread
skew Distribution asymmetry sum((x-mean(x))³)/(n*sd(x)³) Detecting skew direction
kurtosis Tailedness sum((x-mean(x))⁴)/(n*sd(x)⁴)-3 Outlier propensity
se Standard error sd(x)/sqrt(n) Precision of mean estimate

Key Statistics and Their Interpretation

Central Tendency

  • Mean vs. Median: Differences indicate skewness
  • Trimmed Mean: Removes influence of outliers (default drops top/bottom 10%)

Variability

  • SD vs. MAD: Use MAD when outliers are present
  • Range: Simple but outlier-sensitive

Distribution Shape

  • Skewness:
    • >0: Right-tailed
    • <0: Left-tailed
    • 0: Symmetric
  • Kurtosis (Excess):
    • >0: Heavy-tailed (more outliers than normal)
    • <0: Light-tailed

Practical Examples

Example 1: MPG from mtcars

describe(mtcars$mpg)

Output Interpretation:

   vars  n   mean    sd median trimmed   mad min  max range skew kurtosis   se
1     1 32 20.09 6.03   19.2   19.70 5.41 10.4 33.9  23.5 0.61    -0.37 1.07
  • Right-skewed (mean > median, positive skew)
  • Light-tailed (negative kurtosis)
  • SD (6.03) > MAD (5.41): Suggests some outlier influence

When to Use Which Statistic

Scenario Recommended Statistics
Normal Distribution Mean, SD
Skewed Data Median, IQR, MAD
Outlier Detection MAD, trimmed mean, kurtosis
Parametric Testing Mean, SE
Nonparametric Analysis Median, IQR

Extending the Functionality

Adding IQR

The default describe() doesn’t show IQR, but you can add it:

library(dplyr)
describe(mtcars) %>% 
  mutate(IQR = apply(mtcars, 2, IQR, na.rm = TRUE))

Comparing Groups

Use describeBy() for grouped statistics:

describeBy(mtcars$mpg, group = mtcars$cyl)

Conclusion

R’s describe() function provides a powerful starting point for exploratory data analysis. By understanding each statistic it provides, you can:

  • Detect data quality issues
  • Choose appropriate analysis methods
  • Understand your variables’ distributions
  • Make informed decisions about data transformations

For formal reporting, consider supplementing these metrics with visualization and statistical tests.

Pro Tip: Always visualize your data alongside these statistics – numbers tell part of the story, but plots reveal the full picture!

Happy coding!

—
Reference:
Revelle, W. (2023). psych: Procedures for Psychological, Psychometric, and Personality Research. Northwestern University.

Understanding Statistical Coefficients: From Regression to Variation

The Data Analyst’s Guide to  Statistical Coefficients

Table of Contents

  1. What Are Coefficients?
  2. Regression Coefficient
  3. Coefficient of Determination (R²)
  4. Coefficient of Variation (CV)
  5. Correlation Coefficient
  6. Other Common Coefficients
  7. Implementation in R
  8. Key Takeaways

What Are Coefficients?

In statistics and data analysis, coefficients are numerical measures that quantify relationships between variables or characteristics of data distributions. They serve as fundamental indicators in statistical modeling and data interpretation.

1. Regression Coefficient

Definition

The regression coefficient measures the relationship between an independent variable (X) and a dependent variable (Y).

Formula

For linear model Y = aX + b:

  • a: Regression coefficient (change in Y per unit change in X)
  • b: Intercept

R Implementation

# Linear regression example
model <- lm(mpg ~ wt, data = mtcars)
summary(model)

# Extract coefficients
coef(model)

Interpretation

A coefficient of -5.34 for vehicle weight (wt) means each additional ton reduces mileage by 5.34 mpg on average.

2. Coefficient of Determination (R²)

Definition

R-squared represents the proportion of variance in the dependent variable explained by the model (0-1 scale).

R Code

# Get R-squared value
summary(model)$r.squared

Guidelines

  • R² = 0.75 → Model explains 75% of data variation
  • Higher values indicate better model fit

3. Coefficient of Variation (CV)

Definition

CV is a standardized measure of dispersion expressed as percentage of the mean.

Formula

CV% = (Standard Deviation / Mean) × 100%

R Function

# Calculate CV
cv <- function(x) {
  (sd(x, na.rm = TRUE)/mean(x, na.rm = TRUE)) * 100
}

# Example usage
cv(mtcars$mpg)

Interpretation Benchmarks

  • CV < 15%: Low variability
  • 15-30%: Moderate variability
  • >30%: High variability

4. Correlation Coefficient

Definition

Measures the strength and direction of linear relationship between two variables (-1 to 1).

R Implementation

# Calculate correlation
cor(mtcars$mpg, mtcars$wt)

# Correlation matrix
cor(mtcars[, c("mpg", "wt", "hp")])

Interpretation

  • 1: Perfect positive correlation
  • -1: Perfect negative correlation
  • 0: No linear correlation

Other Common Coefficients

Coefficient Description R Package/Function
Skewness Measures distribution asymmetry moments::skewness()
Kurtosis Measures tail heaviness moments::kurtosis()
Concordance Assesses agreement epiR::epi.ccc()

Implementation in R

Comprehensive Analysis

library(psych)

# Descriptive statistics (includes multiple coefficients)
describe(mtcars)

# Full regression output
summary(lm(mpg ~ ., data = mtcars))

Custom Coefficient Calculations

# Multi-coefficient function
data_analysis <- function(x) {
  list(
    mean = mean(x),
    sd = sd(x),
    cv = cv(x),
    skewness = moments::skewness(x),
    kurtosis = moments::kurtosis(x)
  )
}

lapply(mtcars[, 1:4], data_analysis)

Visualization

library(ggplot2)
ggplot(mtcars, aes(wt, mpg)) + 
  geom_point() + 
  geom_smooth(method = "lm") +
  labs(title = "MPG vs Weight with Regression Line",
       x = "Weight (tons)",
       y = "Miles per Gallon")

Key Takeaways

  1. Select coefficients based on analytical goals:
    • Variable relationships → Regression/Correlation coefficients
    • Model evaluation → R-squared
    • Variability comparison → CV
  2. R advantages:
    • Built-in functions for all major coefficients
    • Seamless integration of statistical and visual analysis
  3. Best practices:
    • Understand assumptions behind each coefficient
    • Combine statistical results with domain knowledge
    • Clearly distinguish between different coefficients
  4. Advanced applications: 
    # Robust regression (for outlier-resistant coefficients)
library(MASS)
rlm(mpg ~ wt, data = mtcars)

# Standardized coefficients
library(lm.beta)
lm.beta(model)

By mastering these statistical coefficients and their R implementations, you’ll be equipped to conduct more rigorous data analysis and communicate results effectively. Remember that coefficients are tools – their proper interpretation always depends on context and research questions.

Happy coding!

Producing Systematic Literature Reviews with Bibliometrix R and Biblioshiny

Summer School in Science Mapping (SSSM) 2025 – I International Edition

Title: Producing Systematic Literature Reviews with Bibliometrix R and Biblioshiny
Date & Location: June 9-13, 2025 – Naples, ITA

We are pleased to announce the upcoming Summer School in Science Mapping (SSSM) 2025 – I International Edition, an intensive training program focused on conducting Systematic Literature Reviews using the Bibliometrix R package and its shiny-app Biblioshiny.

Organized by the academic spin-off K-Synth in collaboration with the Department of Economics and Statistics at the University of Naples Federico II, the school will be held in Naples, Italy, from June 9 to June 13, 2025.

Aim and Scope

The SSSM 2025 is an intensive training program tailored for early-career researchers and academics seeking to enhance their expertise in bibliometric methods and scientific mapping. By integrating theoretical foundations with practical sessions, the school equips participants with robust skills in citation analysis, co-citation techniques, science mapping, and reproducible workflows for scholarly evaluation. Designed as both a learning and networking opportunity, SSSM 2025 fosters methodological development and international collaboration in a dynamic, research-oriented environment.

The school’s content covers:

– Overview of bibliometric theory and methods
– Query design and data retrieval from major scientific databases
– Descriptive, relational, and structural bibliometric analyses in R
– Practical training in Bibliometrix R package and Biblioshiny app
– Applications to real-world research review cases

Lecturers and Guest Speakers

The school will be led by Professors Massimo Aria and Corrado Cuccurullo, the developers of Bibliometrix and Biblioshiny.

Additionally, the 2025 edition will feature the following keynotes by distinguished international scholars in the field of scientometrics:

– Nicolas Robinson-Garcia (University of Granada), Scientific Director of the Computational Social Sciences and Humanities Unit (U-CHASS)
– Manuel Jesús Cobo Martín (University of Cádiz), lead developer of the SciMAT software
– Nicola De Bellis (University of Modena and Reggio Emilia), Coordinator of the Bibliometric Office and author of influential studies in the evaluation of scientific research

Target Audience and Prerequisites

This Summer School is designed for PhD students, postdoctoral researchers, and academics affiliated with universities or research institutions. Participants are expected to have a basic knowledge of R programming and be familiar with RStudio.

Registration and Fees

Registration is open on the official Bibliometrix website (check the Summer School section):

https://www.bibliometrix.org/sssm/

For any inquiries, feel free to contact the organizing committee at: [email protected]

R programming book in Greek language

For anyone interested, “Programming in R” (title translated) is a free book on R programming written in Greek.

It presents a programmer’s point of view of R, for beginners (in fact for people with absolutely no programming experience) to advanced programmers. 

Thus, the book does not delve deeply on data science, machine learning, statistics and other such topics where R is broadly used.

Its emphasis is on R as a programming language.

The book may be useful to people who are interested in R programming and can read (or are willing to translate) Greek text.

Chapter titles:

1. Getting started with R
2. The basic elements of R language
3. The essential tools of an R programmer
4. Common object types
5. Functions and functional programming
6. Classes and object-oriented programming
7. Collaboration with other programming languages
8. Package creation
9. Data and content

Link: https://repository.kallipos.gr/handle/11419/8588?&locale=en (in English, with link to download PDF)
DOI: http://dx.doi.org/10.57713/kallipos-100
ISBN: 978-618-5667-90-0

TALL – Text Analysis for ALL, a new R Shiny app for NLP and Text Mining workflows

TALL – Text Analysis for ALL is an R Shiny app that includes a wide set of methodologies specifically tailored for various text analysis tasks. It aims to address the needs of researchers without extensive programming skills, providing a versatile and general-purpose tool for analyzing textual data. With TALL, researchers can leverage a wide range of text analysis techniques without the burden of extensive programming knowledge, enabling them to extract valuable insights from textual data in a more efficient and accessible manner.

Setup

TALL can be installed in two ways, depending on whether you want the stable version or the latest development version.

Official release

You can install the official release of TALL from the Comprehensive R Archive Network CRAN and updated monthly.

if (!require("pak", quietly=TRUE)) install.packages("pak")
pak::pkg_install("tall")

Development release

If you want access to the most recent features and updates not yet available on CRAN, you can install the development version directly from our GitHub repository with:

if (!require("pak", quietly=TRUE)) install.packages("pak")
pak::pkg_install("massimoaria/tall")

Run Tall

Load the library with:

library("tall")

and then run TALL shiny app with:

tall()

Introduction

In the age of information abundance, researchers across diverse disciplines are confronted with the formidable task of analyzing voluminous textual data. Textual data, encompassing research articles, social media posts, customer reviews, and survey responses, harbors invaluable insights that can propel knowledge advancement in various fields, ranging from social sciences to healthcare and beyond. Researchers endeavor to analyze textual data to unveil patterns, discern trends, extract meaningful information, and gain deeper understandings of diverse phenomena. By leveraging sophisticated natural language processing (NLP) techniques and machine learning algorithms, researchers can delve into the semantic and syntactic structures of texts, perform topic detection, polarity detection, and text summarization, among other analyses. Additionally, the advent of digital platforms and the exponential growth of online content have generated unprecedented volumes of textual data that were previously inaccessible or challenging to acquire.

Researchers can harness the power of these textual resources to delve into novel research questions, corroborate existing theories, and generate groundbreaking insights. Through the utilization of computational tools and methodologies, researchers can efficiently process and analyze expansive volumes of text, substantially reducing the time and effort expended compared to manual analysis. Furthermore, there is a burgeoning recognition of the need for text analysis tools tailored to individuals who may not possess in-depth programming expertise. While programming languages like R and Python offer powerful capabilities for data analysis, not all researchers have the time or resources to acquire proficiency in these languages. To address this challenge, a growing number of user-friendly text analysis tools have emerged, providing researchers with a viable alternative to traditional programming-based approaches. These tools empower researchers from diverse backgrounds to effectively process and analyze textual data, fostering a more inclusive research environment and democratizing access to the transformative power of text analysis.

For researchers who lack programming skills, TALL offers a viable solution, providing an intuitive interface that allow researchers to interact with data and perform analyses without the need for extensive programming knowledge.

TALL offers a comprehensive workflow for data cleaning, pre-processing, statistical analysis, and visualization of textual data, by combining state-of-the-art text analysis techniques into an R Shiny app.

TALL workflow

First TALL seamlessly integrates the functionalities of a suite of R packages designed for NLP tasks with the user-friendly interface of web applications through the Shiny package environment.

The TALL workflow streamlines the discovery and analysis of textual data by systematically processing and exploring its content. This comprehensive framework empowers researchers with a versatile toolkit for text analysis, enabling them to efficiently navigate and extract meaningful insights from large volumes of textual data.

By leveraging the strengths of both R packages and Shiny’s interactive web interface, TALL provides a powerful and accessible platform for researchers to conduct thorough the following workflow:

  1. Import and Manipulation

  2. Pre-processing and Cleaning

  3. Statistical Text Analysis and Dynamic Visualization

Some screenshot from TALL

Import text from multiple file formats

Edit, split, and add external information

Automatic Lemmatization and PoS-Tagging through LLM

Language, Model, and Analysis Term Selection

Tagging Special Entities through multiple regex

 

Semantic Tagging

Automatic Multi-word creation

Multi-word creation by a list and Custom Term List

OVERVIEW – Descriptive statistics, concordance analysis and word frequency distributions

WORDS – Multiple methods for Topic Detection

DOCUMENTS – Main approaches for entire texts

“Introduction to R, Regression, and the rms Package”: short course by Frank Harrell

On May 11th 2026, Professor Frank Harrell with lead a workshop, Introduction to R, Regression, and the rms Package, covering foundational R and RStudio skills, linear and multiple regression concepts, and an introduction to the rms package for model fitting and diagnostics. It also introduces a reproducible workflow using Quarto, with a case study demonstration typical for empirical research. This one-day virtual course is offered through Instats.

{talib}: Technical Analysis using R

talib logo

{talib} is a new R package built on TA-Lib, which is now available on CRAN. The R-package is targeted at individuals and, perhaps, institutions who, in some form or the other, interacts with the financial markets using technical analysis.

The library is built with minimal dependencies for long-term stability and freedom in mind. All functions are built around data.frame– and matrix-classes which are portable to all other data-containers with minimal effort.

Everything in the library is built ‘bottom-up’ for maximum speed and memory efficiency. Each indicator interacts directly with R’s C API via .Call().

In this blog post I will give a brief introduction to the interface and the most important aspects of the package. The library also includes static and interactive financial charts, which will be covered in a different post.

A quick introduction to the interface

In this section I will briefly introduce the most important aspects of the function calls, formals and how <NA> are handled. Below is a simple starting point; calculating the Bollinger Bands for Bitcoin:

tail(
  talib::bollinger_bands(
    talib::BTC
  )
)
#>                     UpperBand MiddleBand LowerBand
#> 2024-12-26 01:00:00 100487.38   96698.61  92909.83
#> 2024-12-27 01:00:00 100670.65   96512.96  92355.27
#> 2024-12-28 01:00:00 100632.13   96581.91  92531.69
#> 2024-12-29 01:00:00  99628.77   95576.60  91524.43
#> 2024-12-30 01:00:00  96403.53   94231.31  92059.09
#> 2024-12-31 01:00:00  95441.13   93774.23  92107.34

In itself a simple call. Below are the formals:

str(formals(talib::bollinger_bands))
#> Dotted pair list of 8
#>  $ x        : symbol 
#>  $ cols     : symbol 
#>  $ ma       : language SMA(n = 5)
#>  $ sd       : num 2
#>  $ sd_down  : symbol 
#>  $ sd_up    : symbol 
#>  $ na.bridge: logi FALSE
#>  $ ...      : symbol

All functions share the same signature x, cols, na.bridge and ..., while everything else is indicator-specific. The cols-argument is missing, but has a default value hard-coded against the upstream TA-Lib – this is also true for the remaining indicators. The cols-argument accepts a one-sided formula as follows:

tail(
  talib::bollinger_bands(
    talib::BTC,
    cols = ~close
  )
)
#>                     UpperBand MiddleBand LowerBand
#> 2024-12-26 01:00:00 100487.38   96698.61  92909.83
#> 2024-12-27 01:00:00 100670.65   96512.96  92355.27
#> 2024-12-28 01:00:00 100632.13   96581.91  92531.69
#> 2024-12-29 01:00:00  99628.77   95576.60  91524.43
#> 2024-12-30 01:00:00  96403.53   94231.31  92059.09
#> 2024-12-31 01:00:00  95441.13   93774.23  92107.34

In this case the resulting data.frame is the same as above, as the default column for which the bands are calculated is the closing price of the asset.

All indicators are wrapped by model.frame() and its functionality can be accessed via ... as follows:

tail(
  talib::bollinger_bands(
    talib::BTC,
    cols = ~close,
    subset = 1:nrow(talib::BTC) %in% 1:100
  )
)
#>                     UpperBand MiddleBand LowerBand
#> 2024-04-04 02:00:00  72605.20   68193.82  63782.44
#> 2024-04-05 02:00:00  70646.91   67502.88  64358.84
#> 2024-04-06 02:00:00  70087.36   67344.94  64602.52
#> 2024-04-07 02:00:00  70475.41   68122.29  65769.16
#> 2024-04-08 02:00:00  71820.87   69249.88  66678.89
#> 2024-04-09 02:00:00  71848.20   69372.65  66897.09

Here we only calculate the indicator on a subset of the BTC. While this may seem like a redundant ‘wow’-feature at first glance, its primary justification is in the charting interface where only parts of an indicator is of visual interest.

The <NA>-handling in {talib} works a bit differently than na.rm. Before I demonstrate this, we add some missing values randomly to BTC

BTC <- talib::BTC
BTC$close[sample(1:100, size = 20)] <- NA

The naive approach is to calculate the indicator directly:

tail(
  talib::bollinger_bands(
    BTC
  )
)
#>                     UpperBand MiddleBand LowerBand
#> 2024-12-26 01:00:00        NA         NA        NA
#> 2024-12-27 01:00:00        NA         NA        NA
#> 2024-12-28 01:00:00        NA         NA        NA
#> 2024-12-29 01:00:00        NA         NA        NA
#> 2024-12-30 01:00:00        NA         NA        NA
#> 2024-12-31 01:00:00        NA         NA        NA

Which returns <NA> for all values. This is the default behaviour. The function faithfully returns the full object with the same number of rows – filled with <NA>-values.

To avoid this you can set na.bridge = TRUE as follows:

tail(
  object <- talib::bollinger_bands(
    BTC,
    na.bridge = TRUE
  )
)
#>                     UpperBand MiddleBand LowerBand
#> 2024-12-26 01:00:00 100487.38   96698.61  92909.83
#> 2024-12-27 01:00:00 100670.65   96512.96  92355.27
#> 2024-12-28 01:00:00 100632.13   96581.91  92531.69
#> 2024-12-29 01:00:00  99628.77   95576.60  91524.43
#> 2024-12-30 01:00:00  96403.53   94231.31  92059.09
#> 2024-12-31 01:00:00  95441.13   93774.23  92107.34

Again, with the same number of rows as BTC. This behaviour is true for all functions: N-rows in, N-rows out. What na.bridge is doing under the hood, is to extract all <NA>-values, calculate the indicator and then re-adds them in their original position.

Installation

{talib} is finally on CRAN, and can be installed as follows:

install.packages("talib")

It can also be built from source with additional CMake-flags:

install.packages(
  "talib",
  type = "source",
  configure.args = "-O3 -march=native"
)

Contributing and submitting bug-reports

{talib} is still in its early stage so contributions, even if small, bug-reports, suggestions and critiques are gratefully accepted.

Visit the repository here: https://github.com/serkor1/ta-lib-R.

Created on 2026-04-24 with reprex v2.1.1

New R Package {bdlnm} Released on CRAN: Bayesian Distributed Lag Non-Linear Models in R via INLA

CRAN, GitHub

TL;DR: {bdlnm} brings Bayesian Distributed Lag Non-Linear Models (B-DLNMs) to R using INLA, allowing to model complex DLNMs, quantify uncertainty, and produce rich visualizations.

Background

Climate change is increasing exposure to extreme environmental conditions such as heatwaves and air pollution. However, these exposures rarely have immediate effects. For example:

    • A heatwave today may increase mortality several days later
    • Air pollution can have cumulative and delayed impacts

Distributed Lag Non-Linear Models (DLNMs) are the standard framework for studying these effects. They simultaneously model:

    • How risk changes with exposure level (exposure-response)
    • How risk evolve over time (lag-response)

Usually in the presence of non-linear effects, splines are used to define these two relationships. These two basis are then combined through a cross-basis function. 

As datasets become larger and more complex (e.g., studies with different regions and longer time periods), classical approaches show limitations. Bayesian DLNMs extend this framework by:

    • Supporting more flexible model structures
    • Providing full posterior distributions
    • Enabling richer uncertainty quantification

The new {bdlnm} package extends the framework of the {dlnm} package to a Bayesian setting, using Integrated Nested Laplace Approximation (INLA), a fast alternative to MCMC for Bayesian inference.

Installing and loading the package

As of March 2026, the package is available on CRAN:

install.packages("bdlnm")
library(bdlnm)

At least the stable version of INLA 23.4.24 (or newest) must be installed beforehand. You can install the newest stable INLA version by:

install.packages(
  "INLA",
  repos = c(
    getOption("repos"),
    INLA = "https://inla.r-inla-download.org/R/stable"
  ),
  dep = TRUE
)

Now, let’s load all the libraries we will need for this short tutorial:

Load required libraries 
# DLNMs and splines
library(dlnm)
library(splines)

# Data manipulation
library(dplyr)
library(reshape2)
library(stringr)
library(lubridate)

# Visualization
library(ggplot2)
library(gganimate)
library(ggnewscale)
library(patchwork)
library(scales)
library(plotly)

# Tables
library(gt)

# Execution time
library(tictoc)

Hands-on example

We use the built-in london dataset with daily temperature and mortality (age 75+) from 2000-2012.

Before fitting any model, it is useful to explore the raw data. This plot shows daily mean temperature and mortality for the 75+ age group in London from 2000 to 2012, providing a first look at the time series we are trying to model:

col_mort <- "#2f2f2f"
col_temp <- "#8e44ad"

# Scaling parameters
a <- (max(london$mort_75plus) - min(london$mort_75plus)) /
  (max(london$tmean) - min(london$tmean))
b <- min(london$mort_75plus) - min(london$tmean) * a

p <- ggplot(london, aes(x = yday(date))) +
  geom_line(
    aes(y = a * tmean + b, color = "Mean Temperature"),
    linewidth = 0.4
  ) +
  geom_line(
    aes(y = mort_75plus, color = "Daily Mortality (+75 years)"),
    linewidth = 0.4
  ) +
  facet_wrap(~year, ncol = 3) +
  scale_y_continuous(
    name = "Daily Mortality (+75 years)",
    breaks = seq(0, 225, by = 50),
    sec.axis = sec_axis(
      name = "Mean Temperature (°C)",
      transform = ~ (. - b) / a,
      breaks = seq(-10, 30, by = 10)
    )
  ) +
  scale_x_continuous(
    breaks = yday(as.Date(paste0(
      "2000-",
      c("01", "03", "05", "07", "09", "11"),
      "-01"
    ))),
    labels = c("Jan", "Mar", "May", "Jul", "Sep", "Nov"),
    expand = c(0.01, 0)
  ) +
  scale_color_manual(
    values = c(
      "Daily Mortality (+75 years)" = col_mort,
      "Mean Temperature" = col_temp
    )
  ) +
  labs(x = NULL, color = NULL) +
  guides(color = "none") +
  theme_minimal() +
  theme(
    axis.title.y.left = element_text(
      color = col_mort,
      face = "bold",
      margin = margin(r = 8)
    ),
    axis.title.y.right = element_text(
      color = col_temp,
      face = "bold",
      margin = margin(l = 8)
    ),
    axis.text.y.left = element_text(color = col_mort),
    axis.text.y.right = element_text(color = col_temp)
  ) +
  transition_reveal(as.numeric(date))

animate(p, nframes = 300, fps = 10, end_pause = 100)



Model overview

Conceptually, DLNMs model:

    • Exposure-response: how risk changes with exposure level

    • Lag-response: how risk unfold over time

A typical model is:

Yt ∼ Poisson ( μt )
log ( μt ) = α + cb ( xt , … , x t–L ) · β + ∑ k
γk u kt

where:

    • α is the intercept
    • cb(·) is the cross-basis function, defining both the exposure-response and lag-response relationships
    • β are the coefficients associated with the cross-basis terms
    • ukt are time-varying covariates with corresponding coefficients γk

Model specification & setup

Before fitting the model, we have to define the spline-based exposure-response and lag-response functions using the {dlnm} package.

For our example, we will use common specifications in the literature in temperature-mortality studies:

    • Exposure-response: natural spline with three knots placed at the 10th, 75th, and 90th percentiles of daily mean temperature

    • Lag-response: natural spline with three knots equally spaced on the log scale up to a maximum lag of 21 days

# Exposure-response and lag-response spline parameters
dlnm_var <- list(
  var_prc = c(10, 75, 90),
  var_fun = "ns",
  lag_fun = "ns",
  max_lag = 21,
  lagnk = 3
)

# Cross-basis parameters
argvar <- list(
  fun = dlnm_var$var_fun,
  knots = quantile(london$tmean, dlnm_var$var_prc / 100, na.rm = TRUE),
  Bound = range(london$tmean, na.rm = TRUE)
)

arglag <- list(
  fun = dlnm_var$lag_fun,
  knots = logknots(dlnm_var$max_lag, nk = dlnm_var$lagnk)
)

# Create crossbasis
cb <- crossbasis(london$tmean, lag = dlnm_var$max_lag, argvar, arglag)

As it’s commonly done in these scenarios, we will also control for the seasonality of the mortality time series using a natural spline with 8 degrees of freedom per year, which flexibly controls for long-term and seasonal trends in mortality:

seas <- ns(london$date, df = round(8 * length(london$date) / 365.25))

Finally, we also have to define the temperature values for which predictions will be generated:

temp <- round(seq(min(london$tmean), max(london$tmean), by = 0.1), 1)

Fit the model

Fit the previously defined Bayesian DLNM using the function bdlnm(). We draw 1000 samples from the posterior distribution and set a seed for reproducibility:

tictoc::tic()
mod <- bdlnm(
  mort_75plus ~ cb + factor(dow) + seas,
  data = london,
  family = "poisson",
  sample.arg = list(n = 1000, seed = 5243)
)
tictoc::toc()
8.33 sec elapsed

Internally, bdlnm():

    • fits the model using INLA

    • returns posterior samples for all parameters

Minimum mortality temperature

We estimate the minimum mortality temperature (MMT), defined as the temperature at which the overall mortality risk is minimized. This optimal value will later be used to center the estimated relative risks.

tictoc::tic()
mmt <- optimal_exposure(mod, exp_at = temp)
tictoc::toc()
7.3 sec elapsed

The Bayesian framework, compared to the frequentist perspective, provides directly the full posterior distribution of the MMT. It is useful to inspect this distribution to assess whether multiple candidate optimal exposure values exist and to verify that the median provides a reasonable centering value:

ggplot(as.data.frame(mmt$est), aes(x = mmt$est)) +
  geom_histogram(
    fill = "#A8C5DA",
    bins = length(unique(mmt$est)),
    alpha = 0.6,
    color = "white"
  ) +
  geom_density(
    aes(y = after_stat(density) * length(mmt$est) / length(unique(mmt$est))),
    color = "#2E5E7E",
    linewidth = 1.2,
    adjust = 2 # <-- key change: higher = smoother
  ) +
  geom_vline(
    xintercept = mmt$summary["0.5quant"],
    color = "#2E5E7E",
    linewidth = 1.1,
    linetype = "dashed"
  ) +
  scale_x_continuous(breaks = seq(min(mmt$est), max(mmt$est), by = 0.1)) +
  labs(x = "Temperature (°C)", y = "Frequency") +
  theme_minimal()

The posterior distribution of the MMT is concentrated around 18.9ºC and is unimodal, so the median is a stable centering value for the relative risk estimates.

The posterior distribution of the MMT can also be visualized directly using the package’s plot() method: plot(mmt).

Predict exposure-lag-response effects

We predict the exposure-lag-response association between temperature and mortality from the fitted model at the supplied temperature grid:

cen <- mmt$summary[["0.5quant"]]
tictoc::tic()
cpred <- bcrosspred(mod, exp_at = temp, cen = cen)
tictoc::toc()
6.83 sec elapsed

Centering at the MMT means that relative risks (RR) are interpeted relative to this optimal temperature with minimum mortality.

Several visualizations can be produced from these predictions. While simpler visualizations can be created using the package’s plot() method, here we will use fancier ggplot2 visualizations:

3D exposure-lag-response surface

We can plot the full exposure-lag-response association using a 3-D surface:

matRRfit_median <- cpred$matRRfit.summary[,, "0.5quant"]
x <- rownames(matRRfit_median)
y <- colnames(matRRfit_median)
z <- t(matRRfit_median)

zmin <- min(z, na.rm = TRUE)
zmax <- max(z, na.rm = TRUE)
mid <- (1 - zmin) / (zmax - zmin)

plot_ly() |>
  add_surface(
    x = x,
    y = y,
    z = z,
    surfacecolor = z,
    cmin = zmin,
    cmax = zmax,
    colorscale = list(
      c(0, "#00696e"),
      c(mid * 0.5, "#80c8c8"),
      c(mid, "#f5f0e8"),
      c(mid + (1 - mid) * 0.5, "#c2714f"),
      c(1, "#6b1c1c")
    ),
    colorbar = list(title = "RR")
  ) |>
  add_surface(
    x = x,
    y = y,
    z = matrix(1, nrow = length(y), ncol = length(x)),
    colorscale = list(c(0, "black"), c(1, "black")),
    opacity = 0.4,
    showscale = FALSE
  ) |>
  layout(
    title = "Exposure-Lag-Response Surface",
    scene = list(
      xaxis = list(title = "Temperature (°C)"),
      yaxis = list(title = "Lag", tickvals = y, ticktext = gsub("lag", "", y)),
      zaxis = list(title = "RR"),
      camera = list(eye = list(x = 1.5, y = -1.8, z = 0.8))
    )
  )
Click the image to explore the interactive Plotly version

The surface reveals two distinct risk regions. Hot temperatures produce a sharp, acute risk concentrated at the first lags, peaking at lag 0 and dissipating rapidly after the first lags. Cold temperatures produce a more modest and gradual increase in the first lags that does not fully disappear at longer lags. Intermediate temperatures near the MMT sit close to the RR = 1 reference plane across all lags.

The differential lag structure observed for heat- and cold-related mortality is consistent with known physiological mechanisms. Heat-related mortality tends to occur rapidly after exposure due to acute physiological stress, whereas cold-related mortality develops more gradually through delayed cardiovascular and respiratory effects, leading to increasing risk over longer lag periods.

Lag-response curves

We can also visualizes slices of the previous surface. For example, the lag-response relationship for different temperature values:

matRRfit <- cbind(
  melt(cpred$matRRfit.summary[,, "0.5quant"], value.name = "RR"),
  RR_lci = melt(
    cpred$matRRfit.summary[,, "0.025quant"],
    value.name = "RR_lci"
  )$RR_lci,
  RR_uci = melt(
    cpred$matRRfit.summary[,, "0.975quant"],
    value.name = "RR_uci"
  )$RR_uci
) |>
  rename(temperature = Var1, lag = Var2) |>
  mutate(
    lag = as.numeric(gsub("lag", "", lag))
  )

temp_min <- min(matRRfit$temperature, na.rm = TRUE)
temp_max <- max(matRRfit$temperature, na.rm = TRUE)
mmt_pos <- (cen - temp_min) / (temp_max - temp_min)

temp_values <- c(
  seq(0, mmt_pos, length.out = 6),
  seq(mmt_pos, 1, length.out = 6)[-1]
)

temp_colours <- c(
  "#053061",
  "#2166ac",
  "#4393c3",
  "#92c5de",
  "#d1e5f0",
  "#f7f7f7",
  "#fddbc7",
  "#f4a582",
  "#d6604d",
  "#b2182b",
  "#67001f"
)

p <- ggplot() +
  geom_line(
    data = matRRfit,
    aes(x = lag, y = RR, group = temperature, color = temperature),
    alpha = 0.35,
    linewidth = 0.35
  ) +
  scale_color_gradientn(
    colours = temp_colours, # same 11-colour palette as before
    values = temp_values,
    limits = c(temp_min, temp_max),
    name = "Temperature"
  ) +
  ggnewscale::new_scale_color() +
  geom_hline(
    yintercept = 1,
    linetype = "dashed",
    color = "grey30",
    linewidth = 0.5
  ) +
  scale_x_continuous(breaks = cpred$lag_at) +
  scale_y_continuous(trans = "log10", breaks = pretty_breaks(6)) +
  labs(
    title = "Lag-response curves by temperature",
    x = "Lag (days)",
    y = "Relative Risk (RR)"
  ) +
  theme_minimal() +
  theme(legend.position = "top", panel.grid.minor.x = element_blank()) +
  transition_states(
    temperature,
    transition_length = 1,
    state_length = 0
  ) +
  shadow_mark(past = TRUE, future = FALSE, alpha = 0.6)

animate(p, nframes = 300, fps = 15, end_pause = 100)

Cold temperatures (blue) gradually increase in the initial lags and then decline gradually without fully disappearing in the longer lags. Hot temperatures (red) show a different pattern: a higher risk immediately after lag 0, which drops rapidly and largely dissipates after the first lags:



Exposure-responses curves

We can also plot the exposure-responses curves by lag and the overall cumulative curve across all the lag period:

allRRfit <- data.frame(
  temperature = as.numeric(rownames(cpred$allRRfit.summary)),
  lag = "overall",
  RR = cpred$allRRfit.summary[, "0.5quant"],
  RR_lci = cpred$allRRfit.summary[, "0.025quant"],
  RR_uci = cpred$allRRfit.summary[, "0.975quant"]
)

RRfit <- rbind(matRRfit, allRRfit)

# Split data
RRfit_lags <- RRfit |>
  filter(!lag %in% c("overall")) |>
  mutate(lag = as.numeric(lag))
RRfit_overall <- RRfit |>
  filter(lag %in% c("overall"))

temps <- cpred$exp_at
t_cold <- temps[which.min(abs(temps - quantile(temps, 0.01)))]
t_hot <- temps[which.min(abs(temps - quantile(temps, 0.99)))]

# --- Top plot ---
p_main <- ggplot() +
  # Background: all lags, fading from vivid (small) to pale (large)
  geom_line(
    data = RRfit_lags,
    aes(x = temperature, y = RR, group = lag, color = lag),
    linewidth = 0.8
  ) +
  scale_color_gradientn(
    colours = c(
      "black",
      "#2b1d2f",
      "#4a2f5e",
      "#6a4c93",
      "#8b6bb8",
      "#b39cdb",
      "#d8c9f1",
      "#f3eef5"
    ),
    values = scales::rescale(c(0, 0.5, 1, 2, 3, 4, 5, 10, 20))
  ) +
  new_scale_color() +
  new_scale_fill() +
  # Credible intervals
  geom_ribbon(
    data = RRfit_overall,
    aes(
      x = temperature,
      ymin = RR_lci,
      ymax = RR_uci,
      fill = "1"
    ),
    alpha = 0.2
  ) +
  # Highlighted curves
  geom_line(
    data = RRfit_overall,
    aes(x = temperature, y = RR, color = "1"),
    linewidth = 1.2
  ) +
  geom_hline(
    yintercept = 1,
    linetype = "dashed"
  ) +
  scale_color_manual(values = "#a6761d", labels = "Overall (CrI95%)") +
  scale_fill_manual(values = "#a6761d", labels = "Overall (CrI95%)") +
  scale_y_continuous(
    transform = "log10",
    breaks = sort(c(0.8, pretty_breaks(5)(c(0.8, 4))))
  ) +
  labs(
    x = NULL,
    y = "Relative Risk (RR)",
    color = NULL,
    fill = NULL
  ) +
  theme_minimal() +
  theme(
    legend.position = "top",
    axis.text.x = element_blank(),
    plot.margin = margin(8, 8, 0, 8)
  )

# --- Bottom plot: histogram with annotated percentile lines ---
# center the color palette to the MMT temperature
temp_min <- min(cpred$exp_at, na.rm = TRUE)
temp_max <- max(cpred$exp_at, na.rm = TRUE)
mmt_pos <- (cen - temp_min) / (temp_max - temp_min)

temp_values <- c(
  seq(0, mmt_pos, length.out = 6),
  seq(mmt_pos, 1, length.out = 6)[-1]
)

temp_colours <- c(
  "#053061",
  "#2166ac",
  "#4393c3",
  "#92c5de",
  "#d1e5f0",
  "#f7f7f7",
  "#fddbc7",
  "#f4a582",
  "#d6604d",
  "#b2182b",
  "#67001f"
)

p_hist <- ggplot(london, aes(x = tmean)) +
  geom_histogram(
    aes(y = after_stat(density), fill = after_stat(x)),
    binwidth = 0.5,
    color = "black",
    linewidth = 0.2
  ) +
  geom_vline(
    xintercept = t_cold,
    linetype = "dashed",
    color = "#053061",
    linewidth = 0.6
  ) +
  geom_vline(
    xintercept = t_hot,
    linetype = "dashed",
    color = "#67001f",
    linewidth = 0.6
  ) +
  geom_vline(
    xintercept = cen,
    linetype = "dashed",
    color = "grey20",
    linewidth = 0.6
  ) +
  annotate(
    "text",
    x = t_cold,
    y = Inf,
    label = "1st pct",
    vjust = 1.4,
    hjust = 1.1,
    size = 3.2,
    color = "#053061"
  ) +
  annotate(
    "text",
    x = t_hot,
    y = Inf,
    label = "99th pct",
    vjust = 1.4,
    hjust = -0.1,
    size = 3.2,
    color = "#67001f"
  ) +
  annotate(
    "text",
    x = cen,
    y = Inf,
    label = "MMT",
    vjust = 1.4,
    hjust = -0.1,
    size = 3.2,
    color = "grey20"
  ) +
  scale_x_continuous(limits = c(temp_min, temp_max)) +
  scale_fill_gradientn(
    colours = temp_colours,
    values = temp_values,
    limits = c(temp_min, temp_max),
    name = "Temperature"
  ) +
  labs(x = "Temperature (°C)", y = "Density") +
  theme_minimal() +
  theme(
    plot.margin = margin(20, 8, 8, 8),
    legend.position = "bottom"
  )

# --- Combine ---
p_main / p_hist + plot_layout(heights = c(3, 1))


The overall cumulative curve (mustard) is clearly asymmetric: risk increases on both sides of the MMT, but the rise is much steeper for hot temperatures than for cold temperatures. The lag-0 curve (black), which reflects the immediate effect, behaves differently for cold than heat: it is below 1 at cold temperatures (reflecting the delayed nature of cold temperature effects) and increases approximately linearly for heat. The histogram confirms that most London days fall between 5°C and 20°C, so extreme temperatures, despite their high individual risks, are relatively rare events.

Attributable risk

We can also calculate attributable numbers and fractions from a B-DLNM, which allows to quantify the impact of all the observed exposures in 75+ years mortality. We compute the number of mortality events attributable to the temperature exposures (attributable number) and the fraction of all the mortality events it constitutes (attributable fraction).

Two different perspectives can be used:

    • Backward (dir = "back"): what today’s deaths were explained by past temperature exposures?

    • Forward (dir = "forw"): what future deaths will today’s temperature exposure cause?

Let’s use the forward perspective, more commonly used:

tictoc::tic()
attr_forw <- attributable(
  mod,
  london,
  name_date = "date",
  name_exposure = "tmean",
  name_cases = "mort_75plus",
  cen = cen,
  dir = "forw"
)
tictoc::toc()
110.12 sec elapsed

Attributable fraction evolution

We can plot the time series of daily attributable fractions (AF):

col_af <- "black"

temp_colours <- c(
  "#053061",
  "#2166ac",
  "#4393c3",
  "#92c5de",
  "#d1e5f0",
  "#f7f7f7",
  "#fddbc7",
  "#f4a582",
  "#d6604d",
  "#b2182b",
  "#67001f"
)

# Compute the position of MMT within the actual temperature range
temp_min <- min(london$tmean, na.rm = TRUE)
temp_max <- max(london$tmean, na.rm = TRUE)
mmt_pos <- (cen - temp_min) / (temp_max - temp_min)

temp_values <- c(
  seq(0, mmt_pos, length.out = 6),
  seq(mmt_pos, 1, length.out = 6)[-1]
)

af_med <- attr_forw$af.summary[, "0.5quant"]

af_min <- min(af_med, na.rm = TRUE) - 0.01
af_max <- max(af_med, na.rm = TRUE) + 0.01

df <- data.frame(
  date = london$date,
  x = yday(london$date),
  year = year(london$date),
  tmean = london$tmean,
  af = af_med
)

ggplot(df, aes(x = x)) +
  geom_rect(
    aes(
      xmin = x - 0.5,
      xmax = x + 0.5,
      ymin = af_min,
      ymax = af_max,
      fill = tmean
    )
  ) +
  scale_fill_gradientn(
    colours = temp_colours,
    values = temp_values,
    limits = c(temp_min, temp_max),
    name = "Temperature (°C)"
  ) +
  geom_line(
    aes(y = af),
    color = col_af,
    linewidth = 0.7
  ) +
  scale_y_continuous(
    name = "Attributable Fraction (AF)",
    breaks = seq(0, 1, by = 0.1),
    limits = c(af_min, af_max),
    expand = c(0, 0)
  ) +
  scale_x_continuous(
    breaks = yday(as.Date(paste0(
      "2000-",
      c("01", "03", "05", "07", "09", "11"),
      "-01"
    ))),
    labels = c("Jan", "Mar", "May", "Jul", "Sep", "Nov"),
    expand = c(0, 0)
  ) +
  facet_wrap(~year, ncol = 3, axes = "all_x") +
  labs(x = NULL) +
  theme_minimal(base_size = 11) +
  theme(
    panel.spacing.x = unit(0, "pt"),
    strip.text = element_text(face = "bold", size = 10),
    legend.position = "top",
    legend.key.width = unit(2.5, "cm")
  )



Sharp spikes in AF exceeding 60% are visible in summer 2003 and 2006, coinciding with the major European heatwaves. In general, summer episodes produce higher and more abrupt peaks in AF, whereas cold winter days are associated with more sustained elevations over time, though less pronounced in magnitude.

Total attributable burden

Summing across the full study period, the table quantifies the total mortality burden attributable to non-optimal temperature exposures in the 75+ population:

rbind(
  "Attributable fraction" = attr_forw$aftotal.summary,
  "Attributable number" = attr_forw$antotal.summary
) |>
  as.data.frame() |>
  round(3) |>
  gt(rownames_to_stub = TRUE)
mean sd 0.025quant 0.5quant 0.975quant mode
Attributable fraction 0.174 0.018 0.139 0.175 0.207 0.176
Attributable number 68857.597 7131.526 55071.066 69178.391 81995.459 69842.155
Over the full 2000-2012 period, approximately 17.5% (95% CrI: 13.9%-20.7%) of all deaths in the London 75+ population were attributable to non-optimal temperatures, corresponding to roughly 69,178 deaths (95% CrI: 55,071-81,996).

Conclusions

The {bdlnm} package provides a powerful and accessible implementation of Bayesian Distributed Lag Non-Linear Models in R. By combining the flexibility of DLNMs with full Bayesian inference via INLA, it enables researchers to better quantify uncertainty and fit complex exposure-lag-response relationships. This makes it a valuable tool for studying the health impacts of climate change and other environmental risks in increasingly data-rich settings.

This framework is not limited to environmental epidemiology. In fact it can be applied to any setting involving time-varying exposures and delayed effects (e.g., market shocks may affect asset prices over several days), making it a powerful and general tool for time series analysis.

Development is ongoing. Upcoming features include:

    • Multi-location analyses: pooling exposure-lag-response curves across different cities or regions within a single model
    • Spatial B-DLNMs (SB-DLNM): explicitly modelling spatial heterogeneity in the exposure-lag-response curves of different regions

The package is on CRAN. Bug reports and contributions are welcome via GitHub.

References

    • Gasparrini A. (2011). Distributed lag linear and non-linear models in R: the package dlnm. Journal of Statistical Software, 43(8), 1-20. doi:10.18637/jss.v043.i08.

    • Quijal-Zamorano M., Martinez-Beneito M.A., Ballester J., Marí-Dell’Olmo M. (2024). Spatial Bayesian distributed lag non-linear models (SB-DLNM) for small-area exposure-lag-response epidemiological modelling. International Journal of Epidemiology, 53(3), dyae061. doi:10.1093/ije/dyae061.

    • Rue H., Martino S., Chopin N. (2009). Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations. Journal of the Royal Statistical Society: Series B, 71(2), 319-392. doi:10.1111/j.1467-9868.2008.00700.x.

    • Gasparrini A., Leone M. (2014). Attributable risk from distributed lag models. BMC Medical Research Methodology, 14, 55. doi:10.1186/1471-2288-14-55.

Agentic coding with R workshop

Join our workshop on Agentic coding with R,  which is a part of our workshops for Ukraine series! 


Here’s some more info: 


Title: Agentic coding with R 

Date: Thursday, April 2nd, 14:00 – 16:00 CET (Rome, Berlin, Paris timezone) 

Speaker: Charles Crabtree is a political scientist and Senior Lecturer in the School of Social Sciences at Monash University. His research sits at the intersection of political behavior, discrimination, and research methods, with work spanning experiments, text analysis, and large-scale observational data.

Description: This workshop introduces agentic coding for R: using AI assistants that can help you plan, write, run, and revise multi-step analysis workflows while keeping your work transparent and reproducible. Using Warp.dev as a concrete interface, we will walk through practical patterns for (1) turning messy research tasks into clear, checkable steps, (2) writing R code safely, (3) generating documentation and analysis notes as you work, and (4) developing a paper trail you can share with coauthors or future you. A key focus is adversarial agentic coding: pairing a “builder” agent with a separate “reviewer” agent that tries to break, audit, and improve the code the first agent produced—stress-testing assumptions, spotting silent failures, and proposing fixes. The emphasis is not on prompt tricks, but on reliable habits: how to constrain the agent, verify outputs, and integrate agentic help into real projects (data cleaning, modeling, tables and figures, and report generation). Participants will leave with copy-paste templates they can reuse immediately.

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?



  • 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?


  • 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!