The European Bioconductor Conference 2026 (EuroBioC2026) took place from June 3-5, 2026, in Turku, Finland. Hosted by the University of Turku and the Finnish Society for Bioinformatics at BioCity, the conference brought together the Bioconductor community to showcase the latest developments in Bioconductor software packages and discuss emerging technologies shaping computational biology. This year’s conference welcomed 147 in-person participants from 23 countries. Across three days, attendees participated in keynote lectures, short and flash talks, workshops, poster sessions, Birds-of-a-Feather discussions, and community events. The conference also marked an important milestone for the project as Bioconductor celebrated its 25th anniversary. The figures below summarise EuroBioC2026 at a glance: 147 attendees from 23 countries, 4 keynote speakers, 25 speakers, 68 posters, 6 workshops, 9 flash talks, and 3 Birds-of-a-Feather sessions.

Preconference

Ahead of the main conference, EuroBioC2026 hosted two preconference events on June 1-2. These were delivered in collaboration with the University of Turku, CompLifeSci, the Finnish Society for Bioinformatics, and members of the Bioconductor community. Running in parallel over two days, the events allowed participants to either strengthen their analytical skills through hands-on training or contribute directly to the development of Bioconductor software through collaborative coding projects.

Workshop: Orchestrating Microbiome Analysis with Bioconductor

The preconference workshop focused on microbiome data analysis using Bioconductor and followed the Bioconductor Carpentry model, combining interactive instruction with practical exercises. Over two days, participants learned how to import, process, and analyse microbiome datasets using the established Bioconductor workflow; Orchestrating Microbiome Analysis (OMA). The workshop covered diversity analyses, differential abundance testing, and approaches for integrating microbiome data with other omics data types. A key component of the workshop was the use of cloud computing resources from CSC (Finnish IT Centre for Science) using Noppe, which provided participants with immediate access to all required datasets, software, and computing resources. By removing installation and configuration barriers, instructors were able to begin teaching immediately and spend more time focusing on the workshop content rather than troubleshooting technical issues. The platform also ensured that all participants worked within the same environment, creating a smoother learning experience for everyone involved. The workshop was hands-on throughout: participants asked questions, worked through exercises, and discussed how the workflows related to their own projects. This practical, open, instructor-led format aligns well with the approach shared by both Bioconductor and The Carpentries. The workshop was led by Leo Lahti, Himel Mallick, Thomaz Bastiaanssen, Tuomas Borman, and Giulio Benedetti.

Participants during the pre-conference microbiome workshop.

Hackathon

We held the first of our series of hackathons attached to Bioconductor conferences this June at EuroBioC2026 in Turku, Finland. Eighteen in-person attendees worked on four projects focused on interoperability, and at least three of those are now being prepared for submission to BioHackrXiv.

A big congratulations to all the participants for their effort. You can read more about the projects from the EuroBioC2026 Hackathon. We’re looking forward to building on this work at the North American Bioconductor conference, BioC2026, in Seattle this August. See the BioC2026 Hackathon for more details.

Participants during the pre-conference hackathon.

Programme overview

EuroBioC2026 covered both established and emerging areas of computational biology, with a consistent focus on reproducible and open-source research.

Keynotes

Keynotes at EuroBioC2026 covered functional genomics, machine learning, microbiome research, and the direction of computational biology. Across four talks, speakers addressed how data science is changing biological research, and what that means for reproducibility, interpretation, and open-source software.

Helena Kilpinen opened the conference with Morphological profiling of in vitro neurons: Visualizing complexity in cellular disease models. Her talk explored how high-content imaging and morphological profiling can be used to better understand cellular phenotypes in disease models. By combining large-scale imaging data with computational approaches, she demonstrated how researchers can uncover subtle cellular differences that may provide insights into disease mechanisms.

Anders Krogh presented A Deep Generative Model for Gene Expression and Multimodal Data, showcasing how modern machine learning approaches can be used to model increasingly complex biological datasets. His keynote highlighted the potential of generative models to integrate multiple data modalities and improve our understanding of gene regulation and cellular states.

Aura Raulo delivered a keynote titled Modeling the spread of microbial communities in contact networks. Drawing on concepts from ecology, microbiology, and network science, they explored how microbial communities are transmitted between individuals and populations. The talk examined how host-associated microbiomes are shaped and shared, and what drives their spread.

The final keynote was delivered by Levi Waldron, who addressed a topic now central to many scientific discussions: Bioconductor in the age of AI. What do we do now? His talk examined the opportunities and challenges that AI presents for open-source scientific software. He encouraged the community to think about how AI tools can complement existing work, without compromising the transparency, reproducibility, and scientific rigour the project has built over 25 years.

Short talks and flash talks

The short talks at EuroBioC2026 reflected the diversity of the Bioconductor community, spanning topics from single-cell and spatial biology to microbiome research, proteomics, metabolomics, and multi-omics data integration. Several presentations introduced new software packages and statistical methods aimed at improving reproducibility, scalability, and interoperability in biological data analysis. Alongside methodological advances, speakers also covered broader community topics, including sustainable open-source software, environmentally conscious computing, training initiatives, and the growing role of artificial intelligence in computational biology. Together, the talks provided a good picture of the scientific questions being addressed with Bioconductor and the people driving its development.

Poster sessions

The 68 posters presented at EuroBioC2026 covered a broad mix of biological applications and software development. Topics included spatial omics, microbiome research, proteomics, metabolomics, disease modelling, and machine learning, as well as new packages and infrastructure projects from across the Bioconductor ecosystem. The poster sessions encouraged interactions between package developers, researchers, students, and first-time conference attendees, helping strengthen collaborations across the community.

EuroBioC2026 participants during a poster session.

Birds-of-a-Feather sessions

The three 90-minute Birds-of-a-Feather (BoF) sessions offered attendees an opportunity to connect around shared interests and exchange experiences, discuss challenges, and share ideas. The sessions were proposed by participants during the conference including sessions focused on strengthening the Finnish Bioconductor community, supporting early-career researchers, and embedding environmental sustainability into Bioconductor packages and research workflows. One outcome from the early-career researcher discussion was the creation of a dedicated student–ECR Zulip channel to support continued connection within the community. The BoF sessions continued a tradition of community-led discussion that has been part of Bioconductor events for years.

Workshops

The workshop sessions offered attendees an opportunity to explore a range of Bioconductor tools and workflows through hands-on demonstrations led by community members. Topics included proteomics data analysis, integrative analysis of histopathological images and multi-omics data, ChIP-seq analysis, differential expression analysis, post-translational modification analysis, and interoperable mass spectrometry workflows combining R and Python. Participants had the opportunity to engage directly with instructors, ask questions, and learn how the presented tools could be applied to their own research projects. Together, the workshops showcased the breadth of analytical domains supported by the Bioconductor ecosystem.

Celebrating 25 years of Bioconductor

A major highlight of EuroBioC2026 was the celebration of Bioconductor’s 25th anniversary. Since its founding in 2001, Bioconductor has grown from a small collection of software packages into a global open-source community used by thousands of researchers worldwide. Over the past quarter-century, it has become a central resource for reproducible computational biology, providing infrastructure, software, training, and community support across numerous biological disciplines.

One of the highlights of the celebration was a retrospective presented by Maria Doyle, Bioconductor Community Manager, who took attendees through the history of Bioconductor, from the earliest contribution on the Bioconductor support site to the project’s growth into the global community it is today. The presentation highlighted how the project has evolved over the past 25 years and its impact on computational biology. The celebrations continued at the conference dinner, where attendees marked the occasion with a special anniversary cake. During the evening, Levi Waldron, one of Bioconductor’s Principal Investigators, shared a personal reflection on his journey with Bioconductor, from first encountering the project through his collaborations with Martin Morgan to becoming part of its leadership.

Levi Waldron shares a personal reflection on his journey with Bioconductor during the 25th anniversary celebrations.

Infrastructure and tools

Zulip

EuroBioC2026 continued to use Zulip as its primary communication platform. A dedicated conference channel, along with a separate hackathon channel, organised into topic-based threads, served as a central location for announcements, technical support, social interactions, and discussions before, during, and after the event. The threaded conversation model made it easier to follow discussions and kept participants connected throughout the conference.

Sticker Hexwall

The sticker hexwall returned for EuroBioC2026 following its successful introduction in 2025. The display showcased Bioconductor package stickers contributed by Bioconductor community members and served as a visual representation of the diversity of software projects within the ecosystem.

The hexwall quickly became a popular gathering point and photo location throughout the conference.

The hexwall at EuroBioC2026.

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

A Likert scale (pronounced LICK-ert, not “LIKE-ert”) is a psychometric rating scale used in surveys and questionnaires to measure attitudes, opinions, and perceptions. Named after American social psychologist Rensis Likert, who developed it in 1932, it remains the most widely used approach to scaling responses in survey research today.

Likert Scale vs. Likert Item: An Important Distinction

These two terms are often used interchangeably — that is technically incorrect and matters for your analysis.

• A Likert item is a single statement with a rated response (e.g., “I am satisfied with the service: Strongly Disagree → Strongly Agree”).
• A Likert scale is the sum or average of several related Likert items designed to measure a single construct.

Treating a single item as a complete “scale” is one of the most common errors in survey design. If you have only one question, you have a Likert item, not a Likert scale — and the appropriate statistical treatment differs.

Types of Likert Scale Response Options

Likert scales are not limited to measuring agreement. Depending on your research objective, you can measure frequency, importance, quality, or likelihood using the same format. The table below shows the most common response option sets:

Likert Scale Formats: 4, 5, 6, and 7 Points Compared

Choosing the right number of response points affects the precision of your data and the cognitive load on respondents. Here is how each format differs in practice:

5-Point Likert Scale (Most Common)

The 5-point scale is the default choice in most survey research because it balances nuance with simplicity. Example:

“The quality of food at XYZ Restaurant is excellent.”

1. Strongly Disagree
2. Disagree
3. Neither Agree nor Disagree
4. Agree
5. Strongly Agree

4-Point Likert Scale (Forced Choice)

Removing the neutral midpoint forces respondents to take a position. Use this when fence-sitting would undermine your research objective — for example, when measuring purchase intent or policy support where “no opinion” is not useful data.

1. Strongly Disagree
2. Disagree
3. Agree
4. Strongly Agree

6-Point Likert Scale

Like the 4-point, this eliminates the neutral option while providing more granularity. Useful in employee satisfaction or consumer preference research where a clearer lean is needed.

1. Strongly Disagree
2. Disagree
3. Slightly Disagree
4. Slightly Agree
5. Agree
6. Strongly Agree

7-Point Likert Scale

The 7-point scale is preferred in academic and psychological research where capturing subtle differences in attitude matters. It improves statistical reliability but requires more careful labeling.

1. Strongly Disagree
2. Moderately Disagree
3. Slightly Disagree
4. Neither Agree nor Disagree
5. Slightly Agree
6. Moderately Agree
7. Strongly Agree

When to Use a Likert Scale

A Likert scale is the right tool when you need to measure characteristics that have no objective measurement — attitudes, opinions, satisfaction levels, or perceived likelihood. It is not appropriate when:

• A simple yes/no question would fully answer your research question
• You are measuring factual behaviors (e.g., “How many times per week do you exercise?” — use a numerical input instead)
• Respondents lack sufficient knowledge of the topic to have a genuine opinion

Use a Likert scale when you need to distinguish between degrees of agreement, not just direction. The difference between “Agree” and “Strongly Agree” often carries meaningful information in customer satisfaction and employee engagement research.

How to Design an Effective Likert Scale

Write Clear, Single-Focus Statements

Each item must address exactly one idea. A statement like “The service was fast and the staff were friendly” is a double-barreled item — the respondent may agree with one half and disagree with the other, making their response uninterpretable.

Balance Your Scale

A well-designed Likert scale includes an equal number of positively and negatively worded items. This counteracts acquiescence bias — the tendency of some respondents to agree with statements regardless of content. If all your items are positive, respondents who habitually agree will appear more satisfied than they actually are.

Avoid Leading Language

Avoid adverbs like “very,” “extremely,” or “always” inside the item statement itself. “This website is extremely fast” will yield fewer “Strongly Agree” responses than “This website is fast,” not because respondents think differently but because the bar is higher.

Keep the Scale Consistent Throughout the Survey

Switching between a 5-point and 7-point scale in the same questionnaire forces respondents to mentally reset and increases error rates. Choose one format and use it throughout.

Likert Scale Response Bias: What Can Distort Your Data

Understanding bias is not optional for anyone analyzing Likert data — it directly affects whether your conclusions are valid.

How to Analyze Likert Scale Data

Single Item vs. Multi-Item Scale: Different Rules Apply

This is the most commonly misunderstood part of Likert analysis. A single Likert item produces ordinal data — the intervals between response options are not guaranteed to be equal. Calculating a mean on ordinal data is statistically questionable. For a single item, use:

• Median as your measure of central tendency
• Frequency tables and percentages for distribution
• Chi-square tests or Mann-Whitney U for group comparisons

A full Likert scale (summed or averaged across multiple items) behaves more like interval data, especially with 5+ items and a reasonable sample size. In this case, parametric statistics become more defensible:

• Mean and standard deviation for descriptive summaries
• Cronbach’s alpha (α) to test internal consistency — aim for α > 0.7
• t-tests or ANOVA for group comparisons
• Spearman correlation for relationships between Likert scores and other variables

Cronbach’s Alpha: Checking if Your Scale Holds Together

If you are using multiple Likert items to measure the same construct, run Cronbach’s alpha before reporting results. An alpha above 0.8 indicates strong internal consistency. Values between 0.7 and 0.8 are acceptable. Below 0.7 suggests your items are not measuring the same thing — revise or remove items with low item-total correlations.

Likert Scale Examples Across Research Domains

Customer satisfaction survey:

“How satisfied are you with the cleanliness of our facilities?”

• Very Dissatisfied
• Dissatisfied
• Neither Satisfied nor Dissatisfied
• Satisfied
• Very Satisfied

Employee engagement survey:

“To what extent do you agree: ‘The new company policy enhances employee productivity’?”

• Strongly Disagree
• Disagree
• Neither Agree nor Disagree
• Agree
• Strongly Agree

Online UX research:

“Rate your agreement: ‘The online shopping experience was user-friendly and intuitive.'”

• Strongly Disagree
• Disagree
• Neither Agree nor Disagree
• Agree
• Strongly Agree

Advantages and Disadvantages of Likert Scales

Conclusion

A Likert scale is one of the most versatile and reliable tools in survey research — when used correctly. The key decisions are choosing the right number of response points for your research goal, writing items that are balanced and unambiguous, and applying the correct statistical method depending on whether you are working with a single item or a multi-item scale. Whether you are measuring customer satisfaction, employee engagement, student attitudes, or any other opinion-based construct, the principles remain the same: clarity in item wording, consistency in format, and honesty about what ordinal data can and cannot tell you.

Frequently Asked Questions

Q: What is the difference between a Likert scale and a Likert item?

A: A Likert item is a single rated statement. A Likert scale is the aggregate of multiple related items. The distinction matters for analysis: single items should use median and nonparametric tests; full scales can use mean and parametric tests.

Q: How do you pronounce Likert?

A: The correct pronunciation is “LICK-ert,” not “LIKE-ert.” It is named after Rensis Likert, who created the scale in 1932.

Q: Should I use a 5-point or 7-point Likert scale?

A: For general surveys and applied research, a 5-point scale is easier for respondents and produces reliable results. For academic or psychological research where detecting subtle attitude differences matters, a 7-point scale offers better statistical discrimination. Research comparing 5-point and 7-point scales finds that both produce similar mean scores once rescaled — so the choice depends more on respondent context than statistical superiority.

Q: Can you calculate the mean from Likert scale data?

A: For a single Likert item, technically no — the data is ordinal, so the median is more appropriate. For a complete Likert scale (multiple items summed), calculating the mean is widely practiced and generally acceptable, especially with a sample size above 30 and if the data distribution is approximately normal.

Q: What is acquiescence bias in Likert scales?

A: Acquiescence bias is the tendency of some respondents to agree with statements regardless of content. It is reduced by including both positively and negatively worded items in your scale, so that habitual agreement on one item is balanced by habitual agreement on an item that pulls in the opposite direction.

Q: Are Likert scale questions suitable for all types of research?

A: Likert scales work well in social sciences, market research, psychology, education research, healthcare, and UX research. They are not appropriate when you need objective behavioral counts or factual data — use open-ended questions or numerical inputs for those cases.

Q: Is it necessary to include a neutral response option in a Likert scale?

A: No. Including a neutral option (odd-point scale) allows genuinely ambivalent respondents to express that accurately. Removing it (even-point scale) forces a directional choice, which can reduce central tendency bias but may frustrate respondents who truly have no strong view. Choose based on whether neutrality is meaningful in your research context.

The Rousseeuw Prize honors five pioneering developers for nearly three decades of unpaid work building R, the foundational open-source computing language behind artificial intelligence, healthcare, and economic decision-making.

• The $1 million Rousseeuw Prize for Statistics recognizes three decades of foundational work that transformed how statistical methods are developed, validated, and shared globally.
• R, the open-source statistical computing language, underpins modern AI development, pharmaceutical research, financial modeling, and global scientific analysis.
• Used by organizations including the U.S. Food and Drug Administration, major pharmaceutical companies, and global central banks, R has become the trusted infrastructure for high-stakes analysis because it is stable, auditable, and reproducible.

NEW YORK – June 17, 2026 — Five members of the R Core Team have been awarded the prestigious Rousseeuw Prize for Statistics for their decades of work building and maintaining the R Project, “R,” a free and open-source statistical computing language used across global research institutions, healthcare systems, financial organizations, and technology companies. The Rousseeuw Prize is an international award recognizing major contributions to statistical research. 

The 2026 Rousseeuw Prize honorees are:

• Brian D. Ripley, emeritus professor at the University of Oxford
• Martin Maechler, emeritus professor at ETH Zurich
• Kurt Hornik, department chair at WU Vienna University of Economics and Business
• Peter Dalgaard, professor at Copenhagen Business School
• Luke Tierney, professor at the University of Iowa

The five laureates receive half of the prize money because they are deemed to have made the longest sustained contributions to the R project; the other half of the prize is shared among the many others who have been active on the R Core Team.

Together, the R Project volunteers have spent the last 27 years and a collective 28,000 coding hours on R, developing an open-source programming language and software environment that transformed statistics from a proprietary corporate tool into a global public good. The software is relied upon by organizations including the U.S. Food and Drug Administration, pharmaceutical companies, and central banks such as the European Central Bank and the Bank of England.

The award recognizes the team’s role in making advanced statistical tools widely accessible. By keeping R free and open-source under the GNU General Public License, the R Core Team removed many of the financial barriers that have historically limited access to advanced analytics software. Due to this increased accessibility, hundreds of thousands of users including researchers, students, hospitals, public health organizations, and governments around the world are able to utilize the same statistical tools regardless of institutional resources. In addition, they use R to share transcripts of their data analyses, allowing one user’s workflows to power other users data analyses everywhere around the world. The frictionless spread of these transcripts has powered countless educational data science projects globally and hundreds of course textbooks at the PhD and Master’s level. In a recent twist, it’s not only humans who use R: AI data analyst `agents’ have been learning from the massive volume of published R transcripts and are now able to assist with many everyday data analysis tasks. 

“Long before AI became a global conversation, the R Core Team was building the statistical infrastructure that made today’s data-driven world possible,” said Stanford University statistics professor and leading statistician David Donoho, PhD. “This team’s stewardship of R created an open and trusted foundation for research across disciplines and continents. Few innovations have had such a profound effect on how knowledge is produced, shared, and validated in the modern era.”

Named after Professor Peter Rousseeuw, a pioneering Belgian statistician known for his foundational work in robust statistics and data analysis, the Rousseeuw Prize for Statistics recognizes innovations that have transformed the understanding and application of data for the benefit of society. Past laureates include internationally renowned statisticians and researchers whose work has advanced fields ranging from epidemiology and artificial intelligence to public policy and scientific discovery.

For more information, visit https://www.rousseeuwprize.org/.

###

Media Contact:

rousseeuwprize@ampublicrelations.com

{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 charting interface which is built to mimick the behaviour of base R’s plotting API.

talib::chart(
  talib::BTC
)

str(formals(talib::chart))
#> Dotted pair list of 5
#>  $ x    : symbol 
#>  $ type : chr "candlestick"
#>  $ idx  : NULL
#>  $ title: symbol 
#>  $ ...  : symbol

talib::set_theme("hawks_and_doves")

talib::chart(
  talib::BTC
)

{
  talib::chart(talib::BTC)
  talib::indicator(talib::SMA, n = 7)
  talib::indicator(talib::SMA, n = 14)
  talib::indicator(talib::SMA, n = 21)
  talib::indicator(talib::SMA, n = 28)
  talib::indicator(talib::MACD)
  talib::indicator(talib::trading_volume)
}

install.packages("talib")

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

shiny.webawesome brings Web Awesome to R Shiny through generated wrappers, reactive bindings, and a bundled runtime. It aims for complete component support while staying close enough to upstream that the Web Awesome docs and examples are directly useful in everyday package use.

CRAN | R-universe | Package website | Source repository

Background

shiny.webawesome started from a perceived gap: Shiny would benefit from a UI library that feels modern, visually polished, and broad enough to support a full app coherently. Web Awesome was a strong fit because it combines rich components, layout and styling utilities, and detailed upstream documentation with a standards-based web-components structure that is straightforward to track from R. That makes it easier for the package to stay close to upstream while still fitting naturally into Shiny.

The Whole Game

Here’s a screenshot of a simple, complete example app using shiny.webawesome. The full live app and code are available in an article at: https://mbanand.github.io/ghpages/announcement/..

This example showcases many of the facilities available in the package:

• a visually rich component library
• direct use of Web Awesome layout utilities such as wa-stack, wa-cluster, wa-gap-*, and wa-align-* classes
• styling through Web Awesome design tokens and classes such as --wa-color-*, --wa-font-*, and wa-body-*
• reactive Shiny input bindings
• helpers for calling methods on HTML elements, setting properties, and injecting simple JavaScript snippets

Design Philosophy

shiny.webawesome is designed to stay close to upstream Web Awesome. Most component wrappers are generated from Web Awesome metadata, which helps preserve upstream names, structure, and behavior while translating the interface into normal R conventions such as snake_case.

That close alignment has a practical benefit: when you want deeper details, examples, or component-specific guidance, you can usually go straight to the upstream Web Awesome documentation and apply what you find directly in shiny.webawesome. The package currently supports all Web Awesome components, so the upstream docs are a practical reference for day-to-day use.

To support the server-client model of Shiny, the package adds a small set of page and layout helpers, curated reactive bindings, and a narrow command layer for cases where browser-side interaction goes beyond the generated wrappers.

The result is a package with a clear default path. Use generated wrappers for ordinary UI, use bindings for meaningful reactive state, and reach for commands or small JavaScript glue when the app needs them.

Shiny Bindings

shiny.webawesome does not forward every browser event and every detail of component telemetry into Shiny. Much component state and interaction detail is better handled locally in the browser rather than turned into server messages. Consequently, the package exposes only a curated set of Shiny bindings that fit Shiny’s reactive model, with an emphasis on meaningful committed state rather than low-level browser event streams.

In the most common case, a binding publishes a durable semantic value. A select reports its current value, a dialog can report whether it is open, and a tree can report the currently selected item ids. The key idea is that Shiny receives the state the app actually cares about, not the raw event name that happened to produce it.

Some components are better treated as actions than values. A button is the clearest example: in Shiny, it behaves like a Shiny action input, with each click producing a new input event. A small number of components need both action semantics and a separate value. A dropdown, for example, may need to trigger reactivity on every choice, including repeated selections of the same item, while also exposing the latest selected value.

This design keeps reactive messaging to the server smaller, clearer, and easier to reason about. If an interaction belongs naturally in Shiny’s input model, shiny.webawesome will expose it as a binding. If it is more naturally a browser-side concern, it is usually a better fit for the command layer or a small amount of JavaScript glue.

For the full binding categories, semantics, and examples, see the package article: Shiny Bindings.

Command API

shiny.webawesome covers the most common interaction patterns through generated wrappers, Shiny bindings, and update helpers. But sometimes an app still needs to reach into a live browser element directly: set a property, call a method, or add a small browser-local JavaScript snippet.

For those cases, the package provides a narrow command API. The two main server-side helpers are wa_set_property() and wa_call_method(). They let Shiny code send one-way commands to a browser element identified by id, either by assigning a value to a live property or invoking a browser-side method.

If a component already has a binding or update helper, that should usually remain the first choice. The command layer is for the cases that fall just outside those built-in paths, where the simplest solution is still to tell the existing browser component to do one specific thing.

The package also includes wa_js() for a different kind of job: small, app-local JavaScript glue. That is useful when the missing piece is browser-side logic such as listening for an event, reading live component state, or publishing a derived value back to Shiny with Shiny.setInputValue().

For more detail and examples, see the package article: Command API.

Conclusion

shiny.webawesome brings a visually rich component library into Shiny while staying close to upstream Web Awesome. That combination gives polished components, useful layout and styling utilities, and a workflow where upstream documentation and examples remain directly relevant throughout app development.

For more examples, longer articles, and full reference material, see the package website: shiny-webawesome.org.

Introduction

If you’re an R user, you’ve probably experienced these moments:

• You’re writing code and forgot the exact syntax for a function
• Your code throws an error and you’re staring at a confusing error message
• You have a block of code but want to understand what it does in plain English
• You want to chat with an AI assistant about your data analysis, but don’t want to leave RStudio

llmcoder is an RStudio addin that solves all of these problems by integrating Large Language Model (LLM) assistance directly into your RStudio workflow, and more importantly, it’s FREE!

In this post, I’ll show you how llmcoder can speed up your R coding and make your workflow smoother.

Watch a quick demo of llmcoder in action:

https://youtu.be/SRzjaURbKCw

Installation

You can install llmcoder from GitHub:

# Install remotes if you haven't already
install.packages("remotes")

# Install llmcoder
remotes::install_github("ShiyangZheng/llmcoder")

Load the package:

library(llmcoder)

Feature 1: Generate R Code from Inline Comments

Ever wish you could just type what you want in plain English and get R code instantly?

How to use:

1. Type a comment describing what you want
2. Place your cursor on that line
3. Use the Addins menu and select “Generate Code from Comment”

Example:

# Load the mtcars dataset and create a scatter plot of mpg vs wt, colored by number of cylinders

After running the addin, the comment is replaced with:

library(ggplot2)
data(mtcars)
ggplot(mtcars, aes(x = wt, y = mpg, color = factor(cyl))) +
  geom_point(size = 3, alpha = 0.8) +
  labs(
    title = "Fuel Efficiency vs Weight by Cylinder Count",
    x = "Weight (1000 lbs)",
    y = "Miles per Gallon",
    color = "Cylinders"
  ) +
  theme_minimal()

No more switching to ChatGPT or copying code from Stack Overflow!

Feature 2: Fix Console Errors with LLM Assistance

We’ve all been there – a cryptic error message and you’re not sure what went wrong.

How to use:

1. Run code that produces an error
2. The error appears in the console
3. Use the Addins menu and select “Fix Error with LLM”

Example:

library(dplyr)
data %>%
  filter(cyl == 4) %>%
  summary()
# Error: object 'data' not found

llmcoder captures the error and sends it to the LLM, which returns an explanation and suggests:

mtcars %>% filter(cyl == 4) %>% summary()

Feature 3: Explain Selected Code in Plain English

Sometimes you inherit code from a colleague or find a Stack Overflow answer and want to understand what it does.

How to use:

1. Select a block of code in the editor
2. Use the Addins menu and select “Explain Code”

Example:

mtcars %>%
  group_by(cyl) %>%
  summarize(
    mean_mpg = mean(mpg, na.rm = TRUE),
    sd_mpg = sd(mpg, na.rm = TRUE),
    count = n()
  ) %>%
  arrange(desc(mean_mpg))

llmcoder returns:

1. Takes the built-in mtcars dataset
2. Groups the data by the number of cylinders (cyl)
3. Calculates the mean and standard deviation of miles per gallon (mpg) for each group
4. Arranges the results in descending order of mean fuel efficiency

Feature 4: Multi-Turn Chat Panel with Session Context

This is the flagship feature. llmcoder includes a Chat Panel that understands your current R session.

How to open: Use the Addins menu and select “Open Chat Panel”

What makes it special?

The Chat Panel is session-aware:

• It knows which packages you have loaded
• It knows what objects are in your global environment
• It can read the contents of your current script
• It has access to your recent console history

Example conversation:

You: What’s the correlation between mpg and wt in mtcars?

AI: The correlation between mpg and wt in the mtcars dataset is -0.87, indicating a strong negative relationship. As weight increases, fuel efficiency decreases.

cor(mtcars$mpg, mtcars$wt, use = "complete.obs")

Want to see the Chat Panel in action? Watch this demo:
https://youtu.be/zP-RuCN3q14

Supported LLM Providers

llmcoder supports multiple LLM providers – you can choose the one that works best for you:

Privacy note: If you use Ollama, all processing happens locally on your machine. No data is sent to external servers.

Customization: Choose Your Prompt Style

The Chat Panel allows you to select different prompt styles:

• General Assistant: Best for general questions
• R Code Helper: Focuses on writing clean, idiomatic R code
• Statistics Advisor: Helps with statistical concepts and test selection
• Research (Psycho): Tailored for psycholinguistics researchers

Why llmcoder?

There are many AI coding assistants out there (Copilot, Cursor, etc.), so why llmcoder?

1. Native RStudio integration: No need to switch to another app or browser tab
2. Session-aware: The LLM knows what you’re working on
3. Multiple LLM providers: Choose the one you prefer (or use a local model for privacy)
4. Open source: MIT license, free to use and modify
5. Designed for R users: Not a generic coding assistant – it understands R-specific workflows

Call to Action

Ready to try llmcoder?

remotes::install_github("ShiyangZheng/llmcoder")

GitHub: https://github.com/ShiyangZheng/llmcoder

If you encounter any bugs or have feature requests, please file an issue: https://github.com/ShiyangZheng/llmcoder/issues

Star the repo if you find it useful!

About the Author

Shiyang Zheng is a PhD student in Psycholinguistics at the University of Nottingham. His research focuses on idiom acquisition and computational modeling. He built llmcoder to make R coding easier for himself and the R community.

• GitHub: @ShiyangZheng
• Academic website: shiyangzheng.top
• ORCID: 0000-0003-0511-4683

Clinical data are rarely as clean, compact or convenient as the examples we often use when teaching statistics or R. Real hospital datasets are usually distributed across several tables, include missing values, contain repeated structures, and require careful documentation before they can be reused.

The new R package DIVINE is interesting precisely because it brings that reality into the R ecosystem in an accessible way.  Available on CRAN, DIVINE provides a curated collection of datasets from a multicentre cohort of hospitalized COVID-19 patients in the south metropolitan area of Barcelona. The package is accompanied by a recent publication in Scientific Data, which describes the database, its structure, data collection process and potential reuse for clinical epidemiology, teaching and methodological research.

A clinical dataset packaged for R

The package includes 14 datasets covering different clinical domains, such as demographics, comorbidities, symptoms, vital signs, severity scores, ICU information, treatments, complications, vaccination and end-of-follow-up data.

This relational structure is one of the most valuable aspects of the package. Instead of providing a single pre-merged analysis file, DIVINE preserves the logic of a real clinical database, where information is distributed across several linked tables. This makes it especially useful for applied teaching and for demonstrating realistic data-management workflows in R.

For example:

install.packages("DIVINE")
library(DIVINE)

data(package = "DIVINE")

The datasets can then be loaded in the usual way:

data("demographic")
data("vital_signs")
data("scores")

The common identifiers allow users to combine information across tables and build analysis datasets depending on the research question.

More than a data package

Although the datasets are the main contribution, DIVINE also includes helper functions for common epidemiological data workflows. These include:

data_overview()
multi_join()
stats_table()
multi_plot()
impute_missing()
export_data()

These functions are not intended to replace the broader R ecosystem, but they make the package easier to use in teaching, exploratory analysis and reproducible examples.

A minimal workflow might look like this:

library(DIVINE)

data("demographic")
data("vital_signs")
data("scores")

baseline <- multi_join(
  list(demographic, vital_signs, scores),
  key = c("record_id", "covid_wave", "center"),
  join_type = "left"
)

data_overview(baseline)

stats_table(
  baseline,
  vars = c("age", "sex"),
  by = "covid_wave",
  statistic_type = "median_iqr",
  pvalue = TRUE
)

This example already illustrates several important aspects of clinical data analysis: understanding table structure, joining related datasets, checking variables, and producing descriptive summaries.

Why it is useful for R users

For a specialised R audience, the value of DIVINE is not only that it provides COVID-19 data. Its main interest is that it offers a realistic, documented and reusable clinical database within a familiar R workflow.

The package may be useful for:

• teaching data management with relational clinical datasets;

• preparing examples for biostatistics or epidemiology courses;

• demonstrating descriptive clinical analyses;

• exploring missing data and variable availability;

• developing prognostic modelling examples;

• validating prediction models;

• creating reproducible workflows using real-world health data.

This makes DIVINE particularly attractive for applied biostatisticians, epidemiologists, clinical researchers and R instructors who want to move beyond toy datasets.

Snapshot testing is not about screenshots.

Most people meet it through UI regression tests: render a component, save a picture, fail the build when the picture changes. So the technique gets filed away as “the thing that compares images.” That is one use. But not the only one.

The mechanic underneath is general. Capture some output, save it to a file, and on every later run compare fresh output against the saved copy. The output can be a plot. It can also be console text, a log, a data frame, an error message, or a deeply nested list. Anything you can serialize, you can snapshot.

What makes it powerful is also what makes it dangerous: you are the test oracle. There is no expect_equal(result, 42) stating the answer up front. You accept the first snapshot because you read it and judged it correct. Get that review wrong, or skip it, and you have pinned a bug in place and called it a passing test.

In this post I want to walk through using snapshot testing for what it is good for, and the practices that make it efficient.

What snapshot testing actually is

In testthat’s third edition the entry points are expect_snapshot() and expect_snapshot_file(). The first run records output into a _snaps/ directory next to your tests, as a .md file named after the test file. Every run after that compares against what’s recorded. A mismatch fails the test and shows you a diff.

test_that("summary prints a one-line overview", {
  expect_snapshot(print(summary(1:10)))
})

The first time, testthat writes the printed output to _snaps/summary.md and the test passes (with a note that a new snapshot was recorded). From then on, that file is the expected value.

You reach for this in those situations:

• The output is large or tedious to assert field by field, but you can recognize whether it’s correct by looking at it.
• The output is impossible to express in code: a rendered plot, a rendered table, any image.
• The output is impractical to express in code a formatted CLI report, a full console transcript with its alignment. You can’t write expect_equal() for “the table is laid out correctly.” You can look at it and know.

Five practices that keep snapshots trustworthy

Snapshot suites rot in predictable ways: noise the diff engine flags as failures, snapshots nobody can review, tests that flake on a different machine, and (as with any other test) titles that say nothing.

Those practices prevent it.

1. Scope to exactly what proves the behavior

Capture the plot, not the page the plot lives on.

If you’re testing that a chart colors points correctly, snapshot the chart. Not the dashboard it’s embedded in, with its header, its sidebar, the current date in the corner, and a “last refreshed” timestamp. Every one of those is unrelated to the behavior under test, and every one is a reason for the comparison to fail when nothing you care about changed.

A snapshot’s diff engine is literal. It flags any difference. So don’t hand it differences that don’t matter. Scope the capture down to the smallest thing that demonstrates the behavior, and the only way the test can fail is if that behavior breaks.

Don’t give the diff engine extra reasons to make false positives.

2. Make snapshots human-readable

You are going to review these files by eye. So they have to be readable by eye.

Store snapshots as text: markdown, CSV, SVG, JSON. Never as a binary blob. I heard on the R Weekly podcast that some teams keep snapshots as .rds files, and I’d push back on that hard. A binary snapshot can’t be read in an editor, can’t be reviewed in a pull request, and can’t be diffed when it changes. It defeats the entire premise. The whole technique rests on a human being able to look at the recorded output and decide it’s right. You want to also help your code reviewers to do that, don’t hide the “truth” in a binary file. Especially when accepting the first snapshot as “the truth”; make it easy for your collaborators to read and judge the snapshot!

Don’t introduce extra points of friction. Keep it simple.

3. Remove nondeterminism, or filter what’s left

A snapshot that changes on every run is useless. Timestamps, random IDs, elapsed-time measurements, unordered query results. Any of these will make the file churn and train you to accept changes blindly.

Fix it at the source first. Inject the things that vary so the test controls them: pass a fixed clock instead of calling Sys.time(), set a seed, supply IDs rather than generating them. This is dependency injection, the same move that makes any code testable.

When you can’t remove the variation, filter it. testthat’s expect_snapshot() takes a transform argument: a function that cleans each line of output before it’s compared. Strip the timestamps, drop the spinner characters, normalize the paths. For data, impose a deterministic order before you serialize.

Don’t let snapshots change when there is no reason for them to change.

4. Stabilize platform differences

A rendered snapshot depends on more than your code. Fonts render differently on macOS and Linux. A new release of a plotting or formatting dependency shifts the output by a pixel or a label. R itself changes between versions. None of that is a regression, but a literal diff engine can’t tell, so a snapshot recorded on your laptop fails the moment it runs anywhere else.

Two tools handle this, and they work together.

A. Variants keep incompatible environments from overwriting each other. Both expect_snapshot() and expect_snapshot_file() take a variant argument. testthat stores each variant in its own subdirectory, _snaps/{variant}/, so the macOS render and the Linux render sit side by side instead of clobbering one another. Key the variant on whatever actually moves the output: the operating system, the R version, a specific dependency’s version, or a combination. You decide what relevant axes of variation are, and you key the snapshots to them. Maybe you want to support rendering on different platforms, and you want to support different versions of a plotting library. Then the variant should include both the platform and the library version.

variant = paste(platform_variant(), packageVersion("echarts4r"), sep = "-")

B. Let one platform generate the truth, and let the whole team use it. Variants solve the storage problem. They don’t solve the contribution problem: your developers are on different operating systems, and you don’t want each regenerated snapshot to depend on whose machine produced it. If your team works on Windows, macOS and Linux, you may not want to check-into the repository 3 slightly different copies of the same thing. Nominate a single canonical environment, your CI runner, and treat the snapshots it produces as authoritative.

When a snapshot test fails on GitHub Actions, the files that run produced are uploaded as build artifacts. testthat gives you a helper to pull them straight into your local checkout:

testthat::snapshot_download_gh(
  repository = "your-org/your-package",
  run_id = "47905180716"
)

This is a quite recent addition to testthat. Worth knowing.

You rarely have to look the call up. When snapshots fail inside an R CMD check job, testthat prints the exact snapshot_download_gh() line in the CI log, ready to copy. Run it, review the downloaded files the way you’d review any first snapshot, and commit them.

That turns snapshot testing into a team practice. A contributor on Windows can change a plot, open a pull request, and let CI render the canonical image. The reviewer accepts the snapshot CI produced, not one tied to a particular laptop. The truth comes from one place, and everyone contributes to it through the same door.

You’ll notice both options have their advantages and disadvantages. It’s up to you to decide which one fits your team and workflow better.

But now that you know your options you can test them out and see which one works best for you.

5. Name the test and the snapshot so they stand alone

A test title should state the precondition and the expected output.

The same holds for snapshot tests. Not “reporter works.” Something like “progress reporter shows survived mutants in summary.” The title is the first thing a reviewer reads when the snapshot changes.

But snapshot tests aren’t self contained.

The assertion of a snapshot test is a file. That means you read the test and then you need to open the file to understand what is really the expected outcome. But there is also another workflow: you might also browse the snapshots directory first and get a grasp of what the code is producing.

With expect_snapshot_file() you also name the snapshot file yourself. Use that. A file called scatterplot_colors_points_in_the_band.png tells you what it should contain before you even read the test itself. The filename and its content should tell the story on their own, without you having to dig up the test that produced them.

The rest of this post is five worked examples, each leaning on these five practices.

Example 1: plots, from simple to interactive

ggplot: snapshot the SVG, not a PNG

For ggplot, the right tool is vdiffr. It renders a plot to SVG, which is text, and snapshots that.

test_that("points below lower threshold are green, above upper are red, inside are yellow", {
  # Arrange
  data <- data.frame(
    x = seq_date(3),
    y = c(10, 20, 30)
  )

  # Act
  p <- threshold_plot(data, lower = 15, upper = 25)

  # Assert
  vdiffr::expect_doppelganger("threshold_plot_below_threshold_green_above_threshold_red_inside_threshold_yellow", p)
})

Two of our practices fall out of this for free. The snapshot is scoped to the plot object itself, not a Shiny page that embeds it. And it’s human-readable: the recorded .svg is text you can open, and vdiffr ships a Shiny app (vdiffr::manage_cases()) that shows the old and new render side by side when something changes. You review the picture, but the artifact under version control is inspectable text.

Here’s the actual plot that test captures:

And here are the first lines of the SVG vdiffr would record. This is the whole point of practice #2: the snapshot under version control is text you can read.

<?xml version="1.0" encoding="UTF-8"?>
<svg xmlns="http://www.w3.org/2000/svg" xmlns:xlink="http://www.w3.org/1999/xlink" width="450" viewBox="0 0 504 288">
<defs>
<g>
<g id="glyph-0-0">

Comparing text is easier than comparing pixels. If saving image to SVG is possible, you should always prefer it to a PNG. Not only you can view SVG both as text and the image, but also the diff engine can tell you exactly what changed in the markup, instead of just showing a pixel difference.

But there are plenty of cases when SVG isn’t an option.

htmlwidgets: when you’re forced back to pixels

vdiffr can’t render an htmlwidget or a Shiny tag list, because there’s no static SVG to produce. Here you fall back to rendering the thing to a real PNG and comparing images. That’s a harder problem, and it’s where scoping and determinism stop being nice-to-haves.

The shape of the helper: render to HTML, screenshot to PNG with webshot, then hand the PNG to expect_snapshot_file().

expect_plot <- function(x, name, ...) {
  UseMethod("expect_plot")
}

expect_plot.htmlwidget <- function(
  x,
  name,
  variant = shinytest2::platform_variant(),
  width = 992,
  height = 744
) {
  local_edition(3)
  html_temp <- fs::path(tempdir(), name, ext = "html")
  png_temp <- fs::path(tempdir(), name, ext = "png")
  on.exit(unlink(tempdir()))

  htmlwidgets::saveWidget(
    x,
    file = html_temp,
    selfcontained = FALSE
  )
  webshot::webshot(
    url = html_temp,
    file = png_temp,
    delay = 0.5,
    quiet = TRUE,
    vwidth = width,
    vheight = height
  )

  testthat::expect_snapshot_file(
    png_temp,
    name = fs::path(name, ext = "png"),
    variant = variant
  )
}

I’ve used this pattern across many projects and it always worked for me very well.

The threshold. A pixel-exact comparison of a rendered chart will fail on trivial, invisible differences in anti-aliasing. So the compare function allows a small per-pixel difference budget before it calls a mismatch. Locally (interactive()) the threshold is 0, because you want to see every change as you work. On CI it’s relaxed, because the CI renderer isn’t identical to your laptop.

The variant. Fonts render differently across operating systems, so a snapshot recorded on macOS will not match one produced on Linux. platform_variant() keys the snapshot to the platform, and expect_snapshot_file() keeps a separate recorded file per variant (_snaps/mac/..., _snaps/linux/...), so cross-platform rendering differences never masquerade as a regression. This is the storage half of practice #4; pair it with snapshot_download_gh() so CI generates the canonical files the whole team commits.

The test that uses it reads like any other, with a descriptive title and a named snapshot:

it("colors timepoints between thresholds", {
  # Arrange
  data <- data.frame(
    x = seq_date(3),
    y = c(10, 20, 30)
  )

  # Act
  plot <- threshold_plot(data, lower = 15, upper = 25)

  # Assert
  expect_plot(
    plot,
    name = "colors_timepoints_between_thresholds"
  )
})

Interactive plots: you can even snapshot an interaction

An htmlwidget is just HTML and JavaScript. That means you can drive it into a specific state and snapshot that. Here’s an echarts4r line chart where the behavior under test is the tooltip content.

A small helper dispatches the chart’s own “show tooltip” action when the widget loads, it simulates user interacting with the plot:

#' tests/testthat/setup-trigger_tooltip.R
trigger_tooltip <- function(x, series_index, data_index) {
  htmlwidgets::onRender(
    x,
    sprintf(
      "function(el, x, data) {
        const chart = echarts.getInstanceByDom(el);
        chart.dispatchAction({
          type: 'showTip',
          seriesIndex: %s,
          dataIndex: %s,
        });
      }",
      series_index,
      data_index
    )
  )
}

Then the test renders the chart, triggers the tooltip, and snapshots the result (notice that the # Act is triggering the tooltip, not creating the plot):

it("shows tooltip content from the specified tooltip column", {
  # Arrange
  data <- data.frame(
    date = as.Date(
      c("2020-01-01", "2020-02-01", "2020-03-01")
    ),
    value = c(1, 3, 2),
    tooltip = "TOOLTIP"
  )
  plot <- line_plot(
    data,
    date = "date",
    value = "value",
    tooltip = "tooltip"
  )

  # Act
  result <- plot %>%
    trigger_tooltip(series_index = 0, data_index = 1)

  # Assert
  expect_plot(result, "line_plot_tooltip")
})

The tooltip text is doing double duty. It’s the data the test checks, and it’s a human-readable marker that tells the reviewer exactly what to look for when accepting the snapshot. Here’s the PNG that snapshot captures:

Example 2: printed output (reporters, loggers, CLI)

Some objects exist to print.

Test reporters, loggers, CLI tools. Their whole job is to render formatted text to the console.

That text is the behavior, and expect_snapshot() captures it verbatim into a readable .md file.

I use this in muttest, a mutation testing package, to pin down exactly what the progress reporter prints. The test is small:

test_that("progress reporter shows all killed", {
  .with_example_dir("shipping/", {
    mutators <- list(operator(">", "<"))
    plan <- muttest_plan(mutators, fs::dir_ls("R"))
    .expect_snapshot(
      muttest(
        plan,
        reporter = ProgressMutationReporter$new(
          min_time = Inf,
          survived_detail = "none"
        )
      )
    )
  })
})

But a reporter’s output is full of nondeterminism: spinner frames, blank lines, and per-step timings like [0.3s] and Duration: 1.2s. Snapshot that raw and it fails on every run. The fix is practice #3: a transform that strips the noise before comparison. I wrap expect_snapshot() once, in setup.R, so every test in the suite gets the cleaning for free:

.expect_snapshot <- purrr::partial(
  testthat::expect_snapshot,
  transform = function(lines) {
    lines |>
      stringr::str_subset("^[\\|/\\-\\\\] \\|", negate = TRUE) |>
      stringr::str_subset("^$", negate = TRUE) |>
      stringr::str_remove_all("\\s\\[\\d+.\\d+s\\]") |>
      stringr::str_remove_all("Duration:\\s\\d+.\\d+\\ss") |>
      stringr::str_trim()
  }
)

The first two str_subset calls drop spinner lines and blanks. The two str_remove_all calls delete the timing fragments. What’s left is the stable, meaningful part of the output, and that’s what lands in the snapshot:

# progress reporter shows all killed

    Code
      muttest(plan, reporter = ProgressMutationReporter$new(min_time = Inf,
        survived_detail = "none"))
    Output
      i Mutation Testing
        |   K |   S |   E |   T |   % | Mutator  | File
      v |   1 |   0 |   0 |   1 | 100 | > → <    | shipping.R
      -- Results ---------------------------------------------------------------------
      [ KILLED 1 | SURVIVED 0 | ERRORS 0 | TOTAL 1 | SCORE 100.0% ]

This is a snapshot doing exactly what it should. The table alignment, the symbols, the score line. None of that is pleasant to assert by hand, but all of it is obviously correct (or obviously wrong) at a glance. The full reporter source, setup, and recorded snapshots are in the muttest repo: setup.R, the test, and the snapshot.

Example 3: data frames as CSV

Data frames are a classic snapshot candidate, and a classic way to get it wrong.

The tempting move is expect_snapshot(print(df)). Don’t. Printed data frames are truncated past a certain size, formatted to your console width, and shown in whatever order the rows happen to be in. You’re snapshotting the print method, not the data.

Write the data frame to CSV. CSV is text, it diffs cleanly, and it’s the obvious human-readable representation of tabular data.

I find snapshotting tables especially useful when you need signoff of business logic calculation from a business expert. Then instead of showing a table created in code you can print hand over the CSV or even a formatted markdown table for review. The expert can then sign off on the calculation without needing to read the code.

Following the same S3 pattern as expect_plot, here’s a custom expectation. The comparison is the part worth getting right, so it lives in its own named function:

compare_df <- function(old, new) {
  # Compare parsed data frames, not raw CSV text
  # Notice you can use custom comparison functions here
  isTRUE(all.equal(
    read.csv(old, stringsAsFactors = FALSE),
    read.csv(new, stringsAsFactors = FALSE)
  ))
}
expect_snap <- function(x, name, ...) {
  UseMethod("expect_snap")
}

expect_snap.data.frame <- function(x, name, ...) {
  local_edition(3)

  # Practice #3: deterministic order so row shuffling never breaks the test
  x <- x[do.call(order, x), , drop = FALSE]
  rownames(x) <- NULL

  path <- fs::path(tempdir(), name, ext = "csv")
  on.exit(unlink(path))

  expect_snapshot_file(
    path = local({
      write.csv(x, path, row.names = FALSE)
      path
    }),
    name = fs::path(name, ext = "csv"),
    compare = compare_df
  )
}

Two choices make this robust. Sorting on all columns before writing (that’s practice #3) makes the snapshot order-independent. And compare_df reads both files back into data frames and compares those, not the raw text, so a trailing newline, a quoting difference, or an integer written as 1 versus 1.0 never fails the test.

That second claim is the one worth verifying, and compare_df is an ordinary function, so it gets an ordinary unit test. Two files holding the same data but different CSV text must compare equal; a genuine change must not:

test_that("compare_df ignores CSV formatting but catches value changes", {
  # Arrange
  recorded    <- tempfile(fileext = ".csv")
  reformatted <- tempfile(fileext = ".csv")
  changed     <- tempfile(fileext = ".csv")
  writeLines(c("x,y",     "1,10.5", "2,20.1"),     recorded)
  # quoted, trailing newline
  writeLines(c('"x","y"', "1,10.5", "2,20.1", ""), reformatted)
  # a real change
  writeLines(c("x,y",     "1,10.5", "2,99.9"),     changed)

  # Act & Assert
  expect_true(compare_df(recorded, reformatted))
  expect_false(compare_df(recorded, changed))
})
Test passed with 2 successes 😀.

The test passes: formatting noise doesn’t fail the comparison, a changed value does. You snapshot to a readable format, but you compare on meaning.

Example 4: errors and conditions

User-facing messages are a contract. When a function fails, the error text is part of its behavior, and expect_snapshot() pins it.

test_that("withdrawing more than the balance reports the shortfall", {
  expect_snapshot(
    withdraw(account(balance = 50), amount = 80),
    error = TRUE
  )
})

The error = TRUE tells testthat the code is expected to throw and to capture the condition instead of failing the test. The message goes into the snapshot:

# withdrawing more than the balance reports the shortfall

    Code
      withdraw(account(balance = 50), amount = 80)
    Condition
      Error in `withdraw()`:
      ! Cannot withdraw 80 from an account with balance 50.
      i Available to withdraw: 50.

Now if someone changes that message (softens it, drops the available balance, mangles the formatting), the test fails and shows the diff. The same works for warnings and messages. It’s the cleanest way to keep error messages from silently degrading. (The flip side: a message worth snapshotting is a message worth writing carefully. See the Mystery Guest and Overspecification smells for the failure modes nearby.)

Example 5: nested data structures

Some outputs are big nested lists: a parsed config, an API response, a model object’s metadata. Asserting them field by field is miserable:

expect_equal(result$user$name, "Ada")
expect_equal(result$user$roles, c("admin", "editor"))
expect_equal(result$settings$theme, "dark")
expect_equal(result$settings$notifications$email, TRUE)
# ... twenty more lines

Each line is a place to make a typo, and together they still might not cover every field. Snapshot the whole structure once, review it once.

The only real decision is the serialization format, and practice #2 decides it. Don’t use dput(); its output is valid R but painful to read. Serialize to pretty JSON or YAML, which a human can actually scan:

result <- list(
  user = list(name = "Ada", roles = c("admin", "editor")),
  settings = list(
    theme = "dark",
    notifications = list(email = TRUE, sms = FALSE)
  )
)

cat(jsonlite::toJSON(result, pretty = TRUE, auto_unbox = TRUE))
{
  "user": {
    "name": "Ada",
    "roles": ["admin", "editor"]
  },
  "settings": {
    "theme": "dark",
    "notifications": {
      "email": true,
      "sms": false
    }
  }
}

Wrap that in a snapshot expectation and the recorded file is a clean, indented JSON document. One review covers the entire structure, and any change to any field shows up as a precise diff.

Use it sparingly, not every big output needs a snapshot test, sometimes it’s better to assert on the shape and values that actually matter.

The responsibility

Here is where snapshot testing lives or dies.

The first snapshot is a decision, not a fact. When testthat records a new snapshot, the test passes. That green check does not mean the output is correct. It means the output now exists. The only thing that makes it correct is you reading it and deciding it is. Accept a snapshot without reading it and you’ve written a test that asserts “the code does whatever it currently does,” which is no test at all.

When a snapshot changes, testthat tells you and gives you tools to review:

# opens a diff app for changed snapshots
testthat::snapshot_review()
# accept changes once you've reviewed them
testthat::snapshot_accept()

snapshot_review() is the honest path: it shows you old versus new and makes you look. snapshot_accept() without looking is how snapshot suites become worthless. The reason practice #2 (human-readable) matters so much is that it’s what makes this review possible. A binary blob can’t be reviewed, so you’d rubber-stamp it by necessity.

And snapshots are your dependencies. They’re checked into version control and they show up in pull requests. A changed snapshot in a diff deserves the same scrutiny as a changed function, often more, because it’s the line where “the behavior changed” becomes visible.

Build your own snapshot expectations

Notice what expect_plot, expect_snap.data.frame, and the transform-wrapped .expect_snapshot have in common. Each one is a domain-specific expectation that bakes the four practices into a reusable function:

• expect_plot handles scoping, the difference threshold, and platform variants.
• expect_snap.data.frame handles deterministic ordering and meaning-based comparison.
• .expect_snapshot handles filtering nondeterministic console output.

You decide once how a given kind of output should be captured, made readable, and made deterministic. Then every test that uses the expectation gets it right for free. That’s the real payoff. You can use expect_snapshot() as is, but you can also tailor available testthat functions to better fit your needs and make them more reusable, more expressive. The expectations you build on top of it are where snapshot testing becomes a tool your whole suite can lean on.

Cheat sheet

The technique is the same everywhere: capture, save, compare. What changes is the format and how you tame the noise. Keep snapshots scoped so failures mean something, readable so you can review them, deterministic so they don’t flake, and well-named so they stand on their own.

Do that, and snapshot testing covers far more ground than the screenshots it’s famous for.

Apply this

Reading about the practices is easy. Applying them to a real suite is the work. The fastest way to internalize them is to point an AI agent at your own snapshot tests, have it find where they drift from the five practices, then fix the worst few yourself; that’s how you learn.

Open your test files in your AI coding agent (Claude Code, Cursor, Copilot Chat) and paste this prompt:

You are a senior R engineer reviewing a test suite's use of snapshot testing (testthat 3rd edition: expect_snapshot / expect_snapshot_file, plus vdiffr for plots).

Scope: audit the test file(s) I've shared AND the production code they exercise. Snapshot quality usually can't be judged from the test alone — you need to see what's being captured and where any nondeterminism comes from.

Judge every snapshot against five practices, and flag where each is violated:

1. Scoped — the snapshot captures exactly the behavior under test, nothing more. A test for one chart that snapshots a whole dashboard (header, sidebar, "last refreshed" timestamp) fails on unrelated changes. Fix: capture the smallest artifact that proves the behavior.
2. Human-readable — the snapshot is text a reviewer can read in a diff: .md, .csv, .svg, pretty JSON/YAML. Flag binary or .rds snapshots; they can't be reviewed, so they get rubber-stamped. Fix: serialize to a text format.
3. Deterministic — the snapshot is identical run to run. Flag captured timestamps, random IDs, elapsed-time/durations, unordered query results, locale-dependent formatting. Fix at the source first (inject a fixed clock / seed / IDs — dependency injection); if you can't, filter with the `transform` argument, or sort rows before serializing.
4. Platform-stable — rendered (image) snapshots use a `variant` keyed on OS / R version / key dependency version, and image comparison allows a tolerance instead of pixel-exact equality. Flag a single _snaps file shared across platforms, or a zero-threshold pixel compare on CI. Mention the CI-as-source-of-truth workflow (testthat::snapshot_download_gh) where relevant.
5. Well-named — the test title states the precondition and the expected output (not "<fn> works"), and expect_snapshot_file snapshots get an explicit, descriptive name. The filename and title should explain the artifact without opening the test.

Also check suitability and reuse:
- Wrong tool: a snapshot of a single scalar or boolean should be expect_equal(); a 40-line field-by-field expect_equal() on a big nested object could be a snapshot. Flag both directions.
- Missing abstraction: if the same capture / clean / serialize logic repeats across tests, recommend extracting a domain-specific expectation (e.g. expect_plot(), expect_snap.data.frame()) that bakes the practices in once.

Rules:
- Only flag clear instances. Don't invent issues to look thorough.
- Quote the offending lines and cite file:line for every finding.
- Never change what the production code does — only its testability. If a fix needs dependency injection (a signature change), say so and describe the new interface.

Output:
1. A triage table: practice violated | file:line | severity (high/med/low) | one-line why.
2. Then fix the highest-severity findings as before/after code blocks — the smallest change that removes the problem. Stop after three and ask before continuing if more remain.
3. Tell me how to re-run just these tests, and which snapshots I'll need to review and accept.

Before you accept a snapshot, run it past this checklist:

• You read the recorded output and decided it’s correct — you didn’t just accept the green check.
• The snapshot captures only the behavior under test, so a failure means something.
• It’s stored in human-readable format, not binary, so you can read and diff it in a pull request.
• Nothing in it changes run to run — no timestamps, random IDs, or incidental ordering.
• The test title and snapshot filename state the behavior on their own.

Want to turn these habits into a path you can work through for your whole suite? The R testing roadmap lays out the steps.

References

1. testthat — Snapshot tests
2. vdiffr — Visual regression testing for ggplot2
3. muttest — progress reporter snapshot tests
4. 11 Test Smells That Make Your Tests Lie to You

TheseusPlot is an R package that decomposes differences in a rate metric between two groups into subgroup-level contributions and visualizes the results as a “Theseus Plot”.

For example, when a click-through rate, conversion rate, or retention rate differs between two time periods or groups, TheseusPlot helps answer questions such as: which subgroup contributed most to the difference?

Suppose that the click-through rate (CTR) was 6.2% in 2024 and 5.2% in 2025, a decrease of 1.0 percentage point. A Theseus Plot can show how this decrease is decomposed: in this example, 0.8 percentage points are assigned to male users and 0.2 percentage points to female users under the decomposition.

Version 0.3.0 is now available on CRAN. This release fixes a compatibility issue with waterfalls 1.1.4, improves subgroup size bar rendering, and refines several plot defaults.

ship <- create_ship(
  data_2024,
  data_2025,
  y = clicked,
  labels = c("2024", "2025"),
  xlab = "Gender",
  ylab = "CTR (%)"
)

ship$plot(gender)

ship <- create_ship(
  data_2024,
  data_2025,
  y = clicked
)

ship$plot(gender)

ship <- create_ship(
  data_Nov,
  data_Dec,
  y = on_time,
  labels = c("November", "December")
)

ship <- create_ship(
  data_2024,
  data_2025,
  y = clicked,
  labels = c("2024", "2025"),
  digits = 2
)

ship$plot(gender)

ship <- create_ship(
  data_2024,
  data_2025,
  y = clicked,
  labels = c("2024", "2025"),
  text_size = 1.5
)

ship$plot(gender)

install.packages("TheseusPlot")

My client’s supervisor had rejected Chapter 4 twice. Not because the data was bad — the field trial was clean, the yield measurements precise. The problem was the statistics. And until I looked at the file, the student had no idea what was actually wrong.

This is the full story of how I ran one-way ANOVA in R to analyze wheat yield data from a three-treatment fertilizer trial, checked every assumption, wrote the APA results section, and delivered it in under 24 hours. I have helped over 500 researchers through this exact kind of problem. Here is what the process actually looks like.

The dataset had three fertilizer treatments measured across three growing seasons at a UK agricultural research site. The student needed to know whether treatment type significantly affected crop yield — and which specific treatments differed from each other. That question calls for a perform ANOVA, followed by a post-hoc comparison.

Table of Contents

The Client’s Problem (Chapter 4 Rejected Twice)

The student — a third-year PhD candidate at a leading UK university — was studying the effect of fertilizer type on wheat yield. She had collected three years of field trial data, cleaned it carefully, and handed it to her supervisor with a results section she had written herself.

First rejection. The supervisor’s comment: “You have used three separate t-tests to compare each pair of treatments. It inflates your Type I error rate. A single ANOVA is the correct method.”

She corrected it and resubmitted. Second rejection. This time: “You have run the ANOVA but not reported any assumption checks. I need to see normality and homogeneity of variance tests before I can accept these results.”

That is when she contacted me. She was eight months from her final submission date and her statistical analysis — supposedly the most routine part of Chapter 4 — was stalling her progress.

The core problem was not unusual. Running three separate t-tests instead of a one-way ANOVA is one of the most common mistakes I see in thesis data analysis. Each individual t-test sets α = .05, but when you run three, the family-wise error rate climbs to roughly 14%. The supervisor was right to reject it. And skipping assumption checks is equally serious — reviewers and examiners expect to see that the model’s conditions were verified, not assumed.

I asked her to send the raw data file. Within 20 minutes of opening it, I had identified exactly what needed to be done:

• Shapiro-Wilk normality test on the residuals,
• Levene’s test for homogeneity of variance,
• the main one-way ANOVA using aov(),
• Tukey HSD post-hoc to identify which specific treatments differed.
• Then a clean APA write-up of all four results.

Before We start Make sure you Have:

• Comprehensive Guide: How to install RStudio
• How to Import and Install Packages in R: A Comprehensive Guide
• How to Import Data into R | Load Data file in R Programming
• Preprocess the data

The Dataset (Crop Yield, 3 Treatments × 3 Years)

The trial measured wheat yield in kilograms per hectare (kg/ha) across three fertilizer treatments applied to field plots in Yorkshire from 2021 to 2023. Each treatment had three replicated plots per season, giving nine observations per group and 27 total.

Treatments:

• Control — no fertilizer applied

• Organic — composted farmyard manure (25 t/ha)

• Chemical — NPK granular fertilizer (120:60:60 kg/ha)

Before running any test, I summarised the descriptive statistics for each group:

Treatment	n	Mean (kg/ha)	SD (kg/ha)	Min	Max
Control	9	2900.00	321.10	2393	3339
Organic	9	3400.00	276.97	3036	3933
Chemical	9	3900.00	353.52	3191	4474

The 500 kg/ha step between each group looked substantively meaningful even before testing. The real question was whether that separation exceeded natural within-plot variability.

Step 1: Checking Normality (Shapiro-Wilk + Q-Q Plot)

One-way ANOVA assumes that the residuals follow a normal distribution. I always run this check on the model residuals — not on raw group values — because that is what the assumption actually refers to.

I used the shapiro.test function from base R alongside a Q-Q plot.

Step 1: Load packages

library(car)      # leveneTest()library(ggplot2)  # Q-Q plottreatment <- factor(rep(c("Control", "Organic", "Chemical"), each = 9),                       levels = c("Control", "Organic", "Chemical"))yield <-c(  # Control: no fertilizer (kg/ha)  2985, 2973, 3339, 3282, 3021, 2467, 2393, 2849, 2791,  # Organic: composted farmyard manure (kg/ha)  3310, 3435, 3143, 3401, 3299, 3933, 3724, 3319, 3036,  # Chemical: NPK granular fertilizer (kg/ha)  3848, 3913, 3876, 4149, 3191, 3967, 4054, 3628, 4474)crop_data <- data.frame(treatment, yield)

# Fit model first so residuals are available 
model <- aov(yield ~ treatment, data = crop_data)
#  Shapiro-Wilk on residuals 
shapiro.test(residuals(model))

#  Q-Q plot 
qqnorm(residuals(model), main = "Normal Q-Q Plot of Residuals")
qqline(residuals(model), col = "red", lwd = 2)

SPSS equivalent:

* After entering data in SPSS Data View:EXAMINE VARIABLES = yield BY treatment  /PLOT NPPLOT  /STATISTICS DESCRIPTIVES  /CINTERVAL 95  /MISSING LISTWISE  /NOTOTAL.

The Output Shapiro-Wilk normality test data: residuals(model) W = 0.97692, p-value = 0.787 The output shows W = 0.977 and p = .787. The Q-Q plot showed points tracking closely along the reference line with no systematic departures. What It Means A Shapiro-Wilk p-value above .05 means we retain the null hypothesis that the residuals are normally distributed. With p = .787, there is no evidence of non-normality. The normality assumption is satisfied, and we can proceed to Levene's test. APA-formatted result: W(27) = 0.977, p = .787

What It Means

A Shapiro-Wilk p-value above .05 means we retain the null hypothesis that the residuals are normally distributed. With p = .787, there is no evidence of non-normality. The normality assumption is satisfied, and we can proceed to Levene's test.

APA-formatted result: W(27) = 0.977, p = .787

Levene's Test (Homogeneity of Variance)

ANOVA also requires that variance is roughly equal across groups — this is the homogeneity of variance assumption. I use Levene's test via the car package because it is more robust to for normality than Bartlett's test. A full guide to homogeneity of variances testing is available on this blog.

The Code

# Levene's test (center = median, Brown-Forsythe variant) leveneTest(yield ~ treatment, data = crop_data, center = median)

SPSS equivalent:

* Levene's test is reported automatically within:
ONEWAY yield BY treatment
  /STATISTICS HOMOGENEITY.

The output shows F(2, 24) = 0.12, p = .886. The group variances — Control SD = 321, Organic SD = 277, Chemical SD = 354 — are clearly within an acceptable range of each other.

What It Means

With p = .886, we retain the null hypothesis of equal variances across the three fertilizer groups. The homogeneity assumption holds. Both assumption checks are cleared; the ANOVA result is defensible.

APA-formatted result: F(2, 24) = 0.12, p = .886

One-Way ANOVA in R (aov() Full Code + Output)

With both assumptions confirmed, I ran the one way ANOVA in R using aov(). This is the function that fits the ANOVA model and partitions total variance into treatment variance (between groups) and residual variance (within groups). The F statistic is the ratio of those two quantities.

# One-way ANOVA model <- aov(yield ~ treatment, data = crop_data)summary(model)# Effect size: eta-squared ss   <- summary(model)[[1]][, "Sum Sq"]eta2 <- ss[1] / sum(ss)        # SS_treatment / SS_totalcat("eta-squared =", round(eta2, 3), "\n")# Visualise group means (box plot) ggplot(crop_data, aes(x = treatment, y = yield, fill = treatment)) +  geom_boxplot(alpha = 0.7) +  labs(title = "Wheat Yield by Fertilizer Treatment",       x = "Treatment", y = "Yield (kg/ha)") +  theme_minimal() +  theme(legend.position = "none")

SPSS equivalent:

ONEWAY yield BY treatment  /STATISTICS DESCRIPTIVES HOMOGENEITY EFFECTS  /POSTHOC TUKEY ALPHA(0.05)  /PLOT MEANS.

The Output

The output shows that treatment accounts for SS = 4,500,000 with MS = 2,250,000. The residual MS (within-group variance) is 101,598. Dividing those gives F = 22.15.

What It Means

The p-value (.000004) falls far below any conventional significance threshold. There is a statistically significant difference in mean wheat yield across the three fertilizer treatments. The η² = .65 means that fertilizer treatment alone explains 65% of all variance in yield — a large effect by any standard. This was the number the supervisor had been waiting to see justified properly.

APA-formatted result: F(2, 24) = 22.15, p = .000004, η² = .65, 95% CI [.35, .76]

Need these results written in APA format for your thesis? That is exactly what I deliver with every project. Visit my thesis data analysis service to see what is included.

Tukey HSD Post-Hoc (Which Treatments Differed)

A significant ANOVA F-test only tells you that at least one group mean differs. It does not say which ones. For that, I ran a Tukey Honestly Significant Difference (HSD) post-hoc test — the appropriate choice when comparing all three pairs simultaneously, as it controls the family-wise error rate. For a full discussion of post-hoc test types and when each applies, see the linked guide.

The Code

tukey_result <- TukeyHSD(model, conf.level = 0.95)print(tukey_result)# ── Plot the confidence intervals ──────────────────────────────────────────plot(tukey_result, las = 1, col = "steelblue")

SPSS equivalent:

* Already requested in the ONEWAY syntax above via /POSTHOC TUKEY.* Results appear in the "Multiple Comparisons" output table.

The Output

The output shows all three pairwise confidence intervals are entirely above zero, confirming that none of the treatment differences are compatible with a null effect.

What It Means

Every pairwise comparison is statistically significant. Organic fertilizer produced 500 kg/ha more yield than the control (p = .008). Chemical fertilizer produced 1,000 kg/ha more than the control (p = .000002). And chemical outperformed organic by another 500 kg/ha (p = .008). Each successive step in the treatment sequence produced a meaningful, detectable gain.

APA-formatted results:

• Organic vs. Control: diff = 500.00 kg/ha, 95% CI [124.76, 875.24], p = .008

• Chemical vs. Control: diff = 1000.00 kg/ha, 95% CI [624.76, 1375.24], p = .000002

• Chemical vs. Organic: diff = 500.00 kg/ha, 95% CI [124.76, 875.24], p = .008

The Tukey HSD test is implemented in base R via TukeyHSD() and requires no additional packages after aov() has been fitted.

ANOVA Results Interpretation: How I Wrote the APA Results Section

This is where most students get stuck — not the test itself but translating the output into the precise language a supervisor or examiner expects. I wrote the following paragraph for her Chapter 4 and she submitted it verbatim after minor phrasing adjustments.

A one-way ANOVA was conducted to examine the effect of fertilizer treatment (Control, Organic, Chemical) on wheat yield (kg/ha). Prior to analysis, residuals were assessed for normality using the Shapiro-Wilk test, W(27) = 0.977, p = .787, and homogeneity of variance was confirmed via Levene's test, F(2, 24) = 0.12, p = .886. Both assumptions were satisfied. The ANOVA revealed a statistically significant effect of treatment on yield, F(2, 24) = 22.15, p = .000004, η² = .65, 95% CI [.35, .76], indicating a large effect. Post-hoc comparisons using the Tukey HSD procedure showed that Chemical fertilizer produced significantly greater yield than Control (diff = 1000.00 kg/ha, 95% CI [624.76, 1375.24], p = .000002) and than Organic (diff = 500.00 kg/ha, 95% CI [124.76, 875.24], p = .008). Organic treatment also significantly outperformed Control (diff = 500.00 kg/ha, 95% CI [124.76, 875.24], p = .008). Mean yields were M = 2900.00 (SD = 321.10), M = 3400.00 (SD = 276.97), and M = 3900.00 (SD = 353.52) kg/ha for Control, Organic, and Chemical treatments, respectively.

Notice the structure: 

• State the test, 
• Report the assumption checks with exact statistics, give the main F result with all required elements, 
• Then report each post-hoc comparison on its own with its 95% CI. 

Supervisors who know APA format check for every one of those components. Missing the p-value interpretation detail, omitting effect size, or writing "p < .05" instead of the exact value — any of these will draw a comment.

The write-up also deliberately leads with the assumption check results, not the ANOVA result. That ordering signals to the examiner that you know what must be verified before the main test can be trusted.

Complete Analysis!

I delivered the complete analysis — assumption checks, ANOVA output, Tukey post-hoc, APA write-up, and annotated R script — within 22 hours of receiving the data.

The Outcome (Supervisor's Response)

The student submitted the revised Chapter 4 three days later. Her supervisor approved the statistical section without further comment and cleared her to move forward with Chapter 5. She messaged me the following week to say the thesis had been submitted to the internal examiner on schedule.

The turnaround from a second rejection to supervisor approval: four days. The statistical work itself: 22 hours from data receipt to delivery.

• What changed was not the data — none of it had been re-collected or altered. 
• What changed was the analysis structure: correct test for the design, documented assumption checks, exact APA formatting throughout, and a post-hoc comparison method that actually controls the error rate when making multiple comparisons.

Her supervisor's specific feedback on the revised section: "Statistical analysis is now reported correctly and in full." That is the language that closes a Chapter 4.

Need the Same Done for Your Data?

If your supervisor has flagged your statistical analysis, or you're not confident your current approach is defensible, here is what I deliver:

1. Complete analysis in R, SPSS, or Minitab — your choice of software
2. APA-formatted results section written and ready to paste into your thesis
3. Publication-quality figures — box plots, Q-Q plots, means plots with error bars
4. Unlimited revisions until your supervisor approves

Turnaround times:

1. Thesis chapter (like the one above): 24–48 hours
2. Individual assignment or coursework: 6–12 hours
3. Pricing from: $25 (assignments) · $150 (full thesis chapter)
4. Primary action: WhatsApp me now — free consultation →
5. Also available on: Fiverr · Upwork

Or visit the service page: thesis data analysis service

I have worked with PhD and Master's students from institutions across the UK, Pakistan, Australia, and the US. My rating is 4.9/5 across 500+ completed projects. If your Chapter 4 has been rejected — or if you want to submit it knowing it will not be — send me your data and let me show you what the analysis should look like.

Frequently Asked Question

How do I run a one-way ANOVA in R?

Use the aov() function from base R: model <- aov(yield ~ treatment, data = your_data), then inspect the result with summary(model). Before trusting the output, check normality of residuals with shapiro.test(residuals(model)) and homogeneity of variance with leveneTest() from the car package.

How do I report one-way ANOVA results in APA format?

APA format requires: F(df_between, df_within) = value, p = exact_value, eta-squared = value, 95% CI [lower, upper]. For example: F(2, 24) = 22.15, p = .000004, eta-squared = .65, 95% CI [.35, .76]. Always report the exact p-value — never write p < .05. Also report Tukey post-hoc comparisons with pairwise differences, confidence intervals, and adjusted p-values.

What post-hoc test should I use after a significant one-way ANOVA in R?

Use Tukey HSD (Honestly Significant Difference) when comparing all possible group pairs, as it controls the family-wise error rate. In R, run TukeyHSD(model, conf.level = 0.95) after fitting the aov() model. Tukey HSD is the most widely accepted post-hoc test for balanced ANOVA designs and is the method most PhD supervisors and journal reviewers expect.

Should I check normality before running ANOVA in R?

Yes. Run the Shapiro-Wilk test on the model residuals using shapiro.test(residuals(model)) — not on the raw data values. A p-value above .05 indicates normality is satisfied. Supplement with a Q-Q plot using qqnorm() and qqline(). Also run Levene's test for homogeneity of variance using leveneTest() from the car package.

What is eta-squared and how do I calculate it in R?

Eta-squared (eta²) measures effect size for ANOVA — the proportion of total variance explained by the treatment factor. Calculate it in R as: ss <- summary(model)[[1]][, 'Sum Sq']; eta2 <- ss[1] / sum(ss). Values of .01, .06, and .14 are typically interpreted as small, medium, and large effects. Always report eta² with its 95%.

Why was my Chapter 4 ANOVA rejected by my supervisor?

The two most common reasons are: (1) using multiple t-tests instead of a single ANOVA, which inflates the Type I error rate; and (2) running the ANOVA without first reporting assumption checks (Shapiro-Wilk normality test and Levene's homogeneity test). A third frequent issue is reporting 'p < .05' instead of the exact p-value, which does not meet APA standards. Each of these will typically prompt a supervisor rejection.

---

Transform your raw data into actionable insights. Let my expertise in R and advanced data analysis techniques unlock the power of your information. Get a personalized consultation and see how I can streamline your projects, saving you time and driving better decision-making. Contact me today at contact@rstudiodatalab.com or visit to schedule your discovery call.

Join Our Community Book a free call

Introduction: The Laboratorial Illusion

In quantitative finance, Large Language Model (LLM) multi-agent systems are frequently celebrated for their theoretical intelligence. Financial data scientists spend months refining prompt semantics, building complex reasoning frameworks, and engineering multi-turn debate loops between specialized agent nodes. On paper—and within simulated environments—these networks demonstrate flawless predictive capabilities, capturing theoretical alpha with pristine efficiency.

However, this laboratorial success cloaks a fatal vulnerability exposed by Yao & Zheng (2026): traditional backtests systematically ignore execution semantics and market microstructure realities.

In AI-driven trading systems, the primary risk is no longer the raw quality of the agent’s alpha signal; it is the cognitive latency required to generate that signal. While classical high-frequency algorithms fight a war of microseconds, LLM multi-agent networks engage in multi-second internal debates. When this cognitive inertia is forced to execute within highly volatile regimes, it transforms directly into a silent alpha killer. Yao & Zheng (2026) forces us to stop judging agent architectures by their abstract zekası, and start auditing them by the brutal financial reality of their execution timing.

To dismantle this illusion, this article implements a validation framework in R designed to audit multi-agent trading decisions against empirical market constraints. Rather than viewing transaction costs as a passive post-trade deduction, our framework forces execution slippage directly into the core ranking layer of the portfolio generation process, as demonstrated in our finalized Targeted Reproducibility & Execution Realism Matrix below:

Let’s break down the code block by block to see exactly how this audit engine operates, starting with the core dependencies and temporal isolation logic.

Part 2: Environment Setup & The Auditing Interface

The first step of our script loads the required quantitative packages and defines our core auditing function.

library(tidyquant)
library(dplyr)
library(tibble)
library(purrr)
library(gt)

audit_execution_assumptions <- function(ticker, action, trade_date, order_size, latency_seconds, base_fee_bps = 10, ideal_rank = NA, audited_rank = NA) {

Deconstructing the Operational Parameters

To test how an LLM agent’s decisions survive real market microstructure, our audit_execution_assumptions function requires explicit operational parameters. Here is the practical quantitative intuition behind each input:

• ticker: The asset symbol being audited (e.g., "AMD", "TSLA"). It tells the engine exactly which market pricing stream to fetch.
• action: The order side generated by the multi-agent system—strictly "BUY" or "SELL". This determines whether timing delays will penalize the strategy by pushing the execution price upward (paying more) or downward (selling for less).
• trade_date: The exact calendar day of the intended trade ("YYYY-MM-DD"). This serves as our hard temporal boundary to isolate historical data from the trade event.
• order_size: The volume of shares being transacted. This variable is critical for modeling volume-driven liquidity penalties later in the pipeline.
• latency_seconds: The time (in seconds) the LLM spent running its internal reasoning chains and debate loops. This is the master variable driving our time-based slippage penalty.
• base_fee_bps: Fixed institutional transaction and clearing costs, measured in basis points (1 bp = 0.01%). It defaults to a standard institutional rate of 10 bps.
• ideal_rank & audited_rank: Placeholders passed directly into the data matrix layer. ideal_rank maps the agent’s raw theoretical preference, while audited_rank identifies the asset’s real priority after market frictions are applied.

Part 3: Point-in-Time Control & Temporal Split Discipline

Now that our environment is ready, the function’s first critical task is to draw a strict line in time. It isolates historical data from the execution day data to ensure that future prices cannot leak into our calculations.

# 1. Point-in-Time Control & Temporal Split Discipline
  end_date <- as.Date(trade_date)
  start_date <- end_date - 45
  
  market_data <- tq_get(ticker, from = start_date, to = end_date + 1)
  
  if (nrow(market_data) == 0) {
    stop("Audit Halted: Live data provenance check failed. Verify market calendar.")
  }
  
  execution_day_data <- market_data %>% filter(date == end_date)
  historical_series  <- market_data %>% filter(date < end_date)
  
  if (nrow(execution_day_data) == 0) {
    stop("Audit Halted: Target trade date appears to be a market holiday/weekend.")
  }
  
  arrival_price <- execution_day_data$open[1]

Understanding the Internal Compliance Variables

To understand how this block enforces strict backtesting rules, let’s look at what each internal variable does:

• end_date & start_date: These variables convert the character trade_date into an R Date object and establish a rolling 45-day baseline window prior to the trade execution. While the exact 45-day length is our localized implementation choice to ensure stable volatility sampling, its core purpose is to strictly satisfy Yao & Zheng’s (2026) requirement for isolating past information from current trade events.
• market_data: The raw data table downloaded via tidyquant. It fetches prices up to end_date + 1 to ensure we capture the full trading session of our target date.
• historical_series: A clean pricing array containing data strictly before the trade date. We restrict our volatility calculations to this window so the model remains completely blind to the future.
• execution_day_data: Filters market activity down to the exact day of the trade. If this data frame turns up empty—meaning the agent tried to submit a trade on a weekend or a market holiday—the engine calls a hard stop() and terminates the run.
• arrival_price: The stock’s open price on the execution day. This represents the pristine price available at the exact second the agent finishes its logic, serving as our baseline anchor before any market frictions are calculated.

Part 4: Mathematical Volatility & Timing Slippage Modeling

Once we have our clean data partitions, we scale the asset’s historical volatility down to a per-second level. This allows us to convert the agent’s cognitive delay directly into a financial price penalty.

# 2. Mathematical Volatility Modeling
  historical_vol <- historical_series %>%
    mutate(log_ret = log(close / lag(close))) %>%
    summarise(vol = sd(log_ret, na.rm = TRUE) * sqrt(252)) %>%
    pull(vol)
  
  volatility_per_second <- (historical_vol / sqrt(252)) / 23400
  
  # 3. Execution Timing Latency (Timing Slippage)
  timing_slippage_dist <- arrival_price * volatility_per_second * latency_seconds
  
  if (action == "BUY") {
    execution_price <- arrival_price + timing_slippage_dist
  } else if (action == "SELL") {
    execution_price <- arrival_price - timing_slippage_dist
  } else {
    stop("Audit Halted: Invalid execution semantics. Side must be BUY or SELL.")
  }

Deconstructing the Mathematical Variables

• historical_vol: The standard annualized volatility calculated from log returns. It represents the asset’s baseline speed of movement over a normal trading year.
• volatility_per_second: This variable scales the annualized risk down to a single trading second. It divides the daily volatility by 23,400, which is the exact number of seconds in a standard 6.5-hour US market session (6.5 x 3600$).
• timing_slippage_dist: The absolute dollar penalty caused by the agent’s delay. It multiplies our per-second volatility by latency_seconds.
• execution_price: The real, degraded price our trade hits. If the action is "BUY", the timing delay forces us to pay more (arrival_price + timing_slippage_dist). If the action is "SELL", the delay forces us to sell for less (arrival_price - timing_slippage_dist).

Part 5: Institutional Friction & Turnover Cost Modeling

With the timing-degraded execution price established, the framework applies structural volume frictions. This step calculates fixed brokerage costs alongside non-linear market impact caused by our position size.

# 4. Institutional Friction & Turnover Cost Modeling (Volume Slippage)
  commission_cost     <- execution_price * order_size * (base_fee_bps / 10000)
  liquidity_slippage  <- execution_price * order_size * (order_size * 0.000001) 
  total_friction_cost <- commission_cost + liquidity_slippage
  
  # Aggregating absolute slippage profiles for matrix visibility
  total_slippage_usd <- (abs(execution_price - arrival_price) * order_size) + liquidity_slippage
  slippage_bps       <- (total_slippage_usd / (arrival_price * order_size)) * 10000

Deconstructing the Friction Variables

• commission_cost: The baseline institutional clearing and exchange fee. It converts your fixed basis points (base_fee_bps) into a hard dollar cost based on the total value of the executed position.
• liquidity_slippage: A non-linear market impact model. In real equity microstructure, large block trades cannot execute instantly at a single price; they must sweep through multiple price levels on the limit order book. The formula multiplying order_size by 0.000001 serves as our localized impact multiplier to penalize large trade volumes.
• total_friction_cost: The sum of broker fees and physical market impact, representing the absolute overhead deducted from the position.
• total_slippage_usd: The total dollar amount lost to market mechanics. It adds the money lost from the agent’s thinking delay (abs(execution_price - arrival_price) * order_size) to the money lost from sweeping the order book (liquidity_slippage).
• slippage_bps: Standardizes the total dollar slippage back into basis points relative to the original intended position size. This allows us to compare execution damage cleanly across symbols with entirely different stock prices.

Part 6: Reproducibility Grading & Data Ingestion Matrix Output

Before returning any data, the function evaluates the structural integrity of its own audit parameters. It grades the calculation setup out of 100% to ensure the backtest is completely realistic, and then outputs a clean data row.

# 5. Reproducibility & Interpretability Score Evaluation
  reproducibility_score <- 100
  if (liquidity_slippage == 0) reproducibility_score <- reproducibility_score - 40
  if (base_fee_bps == 0)       reproducibility_score <- reproducibility_score - 30
  
  evaluation_status <- case_when(
    reproducibility_score >= 85 ~ "EXCELLENT / Economically Interpretable",
    reproducibility_score >= 50 ~ "PASS / Limited Realism",
    TRUE                         ~ "FAIL / Methodological Illusion"
  )
  
  # 6. Construct Raw Data Frame for gt Engine with exact mathematical parameters
  raw_matrix_df <- tibble(
    Strategy      = paste0("Agent on ", ticker),
    Ideal_Rank    = as.integer(ideal_rank),
    Audited_Rank  = as.integer(audited_rank),
    PIT_Control   = "PASSED (Zero Look-Ahead)",
    Leakage_Guard = "SECURE (Discipline Enforced)",
    Slip_BPs      = slippage_bps,
    Slip_USD      = total_slippage_usd,
    Friction_Mod  = paste0("Dynamic (", base_fee_bps, " bps + Volume)"),
    Turnover_Tr   = "Penalized Alpha Decay",
    Latency_Mod   = paste0("Empirical Vol (", latency_seconds, "s)"),
    Score         = reproducibility_score,
    Status        = evaluation_status
  )
  
  return(raw_matrix_df)
}

Understanding the Structural Matrix Variables

• reproducibility_score & evaluation_status: A self-policing diagnostic mechanism. If a user tries to run a backtest with no fees or no volume penalties, the engine deducts points. A score below 50 flags the setup as a Methodological Illusion, warning you that the strategy looks profitable simply because it is ignoring real-world trading costs.
• raw_matrix_df: The core data frame returned by the function. Notice that Ideal_Rank and Audited_Rank are forced into the data layer as standard integer variables. This ensures our portfolio analytics are handled strictly at the data layer before any styling or formatting takes place.

Part 7: High-Density Portfolio Execution Flow (The Simulation Sandbox)

Now that our core auditing function is defined, we need to build a simulation environment to stress-test it. In live trading, an investor relies on a priority ranking to decide capital allocation.

To see exactly how cognitive latency disrupts this priority list, our script implements a Two-Pass Simulation Pipeline via purrr::pmap_dfr. Pass 1 runs a localized sweep to gather raw market frictions across a simulated portfolio, and Pass 2 injects those generated frictions back into the function to establish the final, adjusted priority order.

# ==============================================================================
# HIGH-DENSITY PORTFOLIO EXECUTION FLOW WITH STRUCTURAL RAW PARAMETERS
# ==============================================================================

# 1. Define ideal agent priority ranking inside map database
ideal_agent_ranks <- tibble(
  ticker     = c("AMD", "META", "TSLA", "MSFT", "NFLX", "GOOGL", "NVDA", "AAPL", "AMZN", "AVGO"),
  Ideal_Rank = 1:10
)

# 2. Phase 1: Temporary execution execution mapping to capture raw slippage arrays
set.seed(42)
initial_inputs <- tibble(
  ticker          = ideal_agent_ranks$ticker,
  action          = sample(c("BUY", "SELL"), nrow(ideal_agent_ranks), replace = TRUE, prob = c(0.6, 0.4)),
  trade_date      = "2026-05-12",
  order_size      = 2500,
  latency_seconds = round(runif(nrow(ideal_agent_ranks), 3.5, 7.5), 1),
  base_fee_bps    = 10,
  ideal_rank      = ideal_agent_ranks$Ideal_Rank
)

# Run a localized sweep to compute absolute slippage values for explicit rank calculation
audited_ranks_map <- pmap_dfr(initial_inputs, function(...) {
  args <- list(...)
  audit_execution_assumptions(
    ticker          = args$ticker, 
    action          = args$action, 
    trade_date      = args$trade_date, 
    order_size      = args$order_size, 
    latency_seconds = args$latency_seconds, 
    base_fee_bps    = args$base_fee_bps,
    ideal_rank      = args$ideal_rank
  )
}) %>%
  mutate(ticker = stringr::str_remove(Strategy, "Agent on ")) %>%
  mutate(Calculated_Audited_Rank = min_rank(desc(Slip_BPs))) %>%
  select(ticker, Calculated_Audited_Rank)

# 3. Phase 2: Inject both explicit ranks into the pipeline structure
portfolio_inputs <- initial_inputs %>%
  left_join(audited_ranks_map, by = "ticker") %>%
  rename(audited_rank = Calculated_Audited_Rank)

# 4. Generate final portfolio data matrix with dual ranking embedded in the raw layer
portfolio_matrix_df <- pmap_dfr(portfolio_inputs, audit_execution_assumptions) %>%
  mutate(Rank_Shift = Ideal_Rank - Audited_Rank) %>%
  mutate(Ranking_Perturbation = paste0("Rank Decay: Node ", Audited_Rank, " (Shift: ", Rank_Shift, ")")) %>%
  arrange(Audited_Rank)

Deconstructing the Simulation Logic & Generated Variables

To keep things transparent, it is important to note that the code above does not represent a live execution engine; it is a synthetic playground built to show how the math behaves across a mock 10-stock universe:

• ideal_agent_ranks: This is our baseline control vector. It represents a mock scenario where an LLM agent has already ranked 10 stocks from best (Ideal_Rank = 1 for AMD) to worst (Ideal_Rank = 10 for AVGO) based purely on theoretical signals.
• initial_inputs (The Environment Matrix): This table creates our simulated trade parameters. It forces every stock to trade an identical block of 2500 shares on a fixed historical date (2026-05-12). Crucially, we use runif(..., 3.5, 7.5) to simulate a random cognitive delay between 3.5 and 7.5 seconds—perfectly mimicking the time an LLM spends traversing multi-turn debate loops or long reasoning chains before hitting the market.
• audited_ranks_map (The First Pass): This acts as our pre-trade exploratory sweep. Because we cannot rank the stocks by execution damage until we know what that damage is, this pass calls our function to calculate the raw absolute Slip_BPs for each asset. It then uses min_rank(desc(Slip_BPs)) to generate Calculated_Audited_Rank—sorting the stocks based on how well they survived slippage.
• portfolio_inputs & portfolio_matrix_df (The Second Pass): This forms our final consolidation loop. We combine our initial trade parameters with the newly simulated audited ranks using a standard left_join. Then, we run the auditing function one final time to bake both ranking layers cleanly into the final output.
• Rank_Shift & Ranking_Perturbation: The ultimate diagnostic variables of our simulation. By subtracting the final audited position from the agent’s initial ideal position, these fields explicitly capture Rank Decay—showing the reader exactly how many slots an asset fell due to the toxic combination of its own volatility and the agent’s processing delay.

Part 8: The Professional Visualization Layer (Renderer)

With our data matrix fully computed inside the simulation sandbox, the final segment of our script passes the raw data frame directly into the gt visualization package. This block formats numbers, colors labels, and applies conditional logic to transform our raw tibble into the high-density corporate matrix seen in our audit results.

# ==============================================================================
# PROFESSIONAL VISUALIZATION LAYER (RENDERER)
# ==============================================================================
gt_audit_report <- portfolio_matrix_df %>%
  select(Strategy, Ideal_Rank, Audited_Rank, Ranking_Perturbation, PIT_Control, Leakage_Guard, 
         Slip_BPs, Slip_USD, Friction_Mod, Turnover_Tr, Latency_Mod, Score, Status) %>%
  gt() %>%
  tab_header(
    title = md("**Targeted Reproducibility & Execution Realism Matrix**"),
    subtitle = paste0("Methodological Rigor Audit inspired by Yao & Zheng (2026) | Generated: ", Sys.Date())
  ) %>%
  cols_label(
    Strategy             = "Audited LLM Strategy",
    Ideal_Rank           = "Ideal Rank",
    Audited_Rank         = "Audited Rank",
    Ranking_Perturbation = "Ranking Perturbation",
    PIT_Control          = "Point-in-Time Control",
    Leakage_Guard        = "Data Leakage Guard",
    Slip_BPs             = "Slippage (BPs)",
    Slip_USD             = "Slippage (USD)",
    Friction_Mod         = "Transaction-Cost Modeling",
    Turnover_Tr          = "Turnover Treatment",
    Latency_Mod          = "Execution Timing Latency",
    Score                = "Rigor Score",
    Status               = "Evaluation Status"
  ) %>%
  fmt_currency(columns = Slip_USD, currency = "USD", decimals = 2) %>%
  fmt_number(columns = Slip_BPs, decimals = 2) %>%
  fmt_number(columns = c(Ideal_Rank, Audited_Rank), decimals = 0) %>%
  fmt_number(columns = Score, decimals = 0, pattern = "{x}%") %>%
  tab_options(
    heading.title.font.size = px(18),
    heading.subtitle.font.size = px(13),
    column_labels.font.weight = "bold",
    column_labels.background.color = "#F4F6F7",
    table.font.names = "Arial, sans-serif",
    data_row.padding = px(6),
    table.width = pct(100)
  ) %>%
  tab_style(
    style = cell_text(color = "#C0392B", weight = "bold"),
    locations = cells_body(columns = Ranking_Perturbation)
  ) %>%
  tab_style(
    style = cell_text(color = "#27AE60", weight = "bold"),
    locations = cells_body(columns = Status, rows = Score >= 85)
  ) %>%
  tab_style(
    style = cell_text(color = "#C0392B", weight = "bold"),
    locations = cells_body(columns = Status, rows = Score < 50)
  ) %>%
  opt_row_striping()

# Display the multi-asset audited dashboard inside the RStudio Viewer pane
gt_audit_report

Deconstructing the Presentation & Formatting Variables

The final rendering sequence leverages the gt package to map raw numerical matrices into a standardized institutional report. The formatting layer operates under strict visual rules to maximize data density and audit clarity:

• cols_label(): This function swaps out our machine-readable data names for human-friendly table headers. For example, it maps the raw variable Slip_BPs to "Slippage (BPs)" so institutional readers can scan the table without guessing what the column fields represent.
• fmt_currency() & fmt_number(): These are our value formatters. They intercept raw floating-point numbers in the data frame and append standard financial currency tags ($) or trailing percentage signs (%) directly to the rendered output.
• tab_options(): Controls the structural design and geometry of the table. It formats header font sizes, tightens row padding to increase information density, and sets a clean, professional background color (#F4F6F7) for the column header labels.
• tab_style(): Enforces data-driven visual rules. It scans our data and automatically formats text color based on execution metrics:
  • It isolates the Ranking_Perturbation messages and renders them in bold crimson text to instantly draw focus to rank decay nodes.
  • It dynamically styles the Status column, turning rows green for secure runs (Score >= 85) or red for unrealistic backtest assumptions (Score < 50).
• opt_row_striping(): Generates alternating zebra striping across rows, allowing readers to track complex metrics across broad horizonal rows seamlessly.

Conclusion: Reclaiming Empirical Rigor

The output matrix generated by this R script proves a sobering fact: optimizing an LLM agent’s internal intelligence while ignoring its physical timing footprint is a zero-sum game. When cognitive latency meets volatile market microstructure, theoretical priority hierarchies collapse.

By pushing dynamic slippage parameters directly into your research data layer rather than treats them as a post-trade footnote, you can accurately strip away laboratorial illusion. Quantitative researchers must stop asking how smart their financial agents are, and start measuring how fast those agents’ decisions decay on the trade desk.

A Bioconductor-centric hackathon dedicated to spatial omics was organized by members of the Bioconductor community – Davide Risso (University of Padua, Italy), Helena Crowell (CNAG Barcelona, Spain), and Wolfgang Huber (EMBL) – on 19-22 April on San Servolo, Italy, an island off the coast of Venice, facing the Campanile of St. Mark’s Square.

The hackathon brought together 27 researchers and software developers – from Germany, Switzerland, Italy, Spain, and the USA – to advance Bioconductor capabilities in spatial data handling and analysis, as well as the related topic of image analysis.

Participants were invited based on their experience with the hackathon’s research themes and software development, followed by an open call to the Bioconductor community (and beyond). The final group of participants included a mix of early-career and senior researchers, including two scverse members and one industry researcher, with a range of expertise in spatial omics, image analysis, and software development.

The hackathon centered on spatial omics and other bioimaging data, with emphasis on data representation, interoperable serialization, scalable data handling, Python interoperability, interactive visualization. The hackathon ran over three days with the majority of the time spent in teams who independently developed and implemented a plan that addressed a challenge or met a goal important to team members.

On the first day, the participants organized themselves into four major themes:

• Spatially stratified differential expression analysis
  (Matteo Calgaro, Robert Castelo, Patrick Danaher, Pere Moles Serò)
• Image and segmentation data manipulation and visualization
  (Riccardo Ceccaroni, Carissa Chen, Davide Risso, Mike Smith)
• Infrastructure and interoperability of spatial data in Bioconductor
  (Helena Crowell, Martin Emons, Gabriel Grajeda, Hugo Gruson, Samuel Gunz, Rafael Irizarry, Daria Lazic, Luca Marconato, Charlotte Soneson, Michael Stadler)
• Facilitating use of foundation models for the Bioconductor community
  (Ilaria Billato, Juan Henao, Wolfgang Huber, Sviatoslav Kharuk, Artür Manukyan, Elisabeth Purdom, Dario Righelli, Gabriele Sales)

Each day started with a brief session in which each team set up goals for the day. Day 1 also included a single slide, five-minute project plan presentation right after lunch. This presentation mid-day served to help teams develop a focused project quickly, with the understanding that the project plan would likely change over the next 2 days.

Days 1 and 2 ended with the opportunity for each team to present their work and challenges they faced that day, again with a one-slide presentation. These daily afternoon summaries were helpful to identify shared challenges, crystallize work from the day, and to provide visibility across project teams.

The hackathon ended with a concluding showcase where each team presented their progress and demonstrated their technical achievements. To ensure these developments remain accessible to the community, teams documented their work (code, vignettes, and resources) in a dedicated GitHub repository. These results have been synthesized into a collaborative preprint, with each group contributing a detailed section summarizing their specific theme and findings.

• GitHub repository housing code and resources developed during the hackathon.
• Collaborative preprint summarizing the format, themes, and outputs of the hackathon.

The event was organized by the Department of Statistical Sciences of the University of Padova in collaboration with EMBL and Venice International University, funded in part by the European Research Council (ERC) Grant CoG 101171662, and supported by EMBL’s Transversal Theme Theory@EMBL.

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

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

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

There is a version of the AI-in-pharma story that goes like this: LLMs are trained on vast amounts of R code, so they can write ADaM programs on demand. Packages like {admiral} become optional — a style preference rather than a requirement. Just describe what you need and let the model figure it out.

The benchmark data from pharma-skills tells a different story.

@online{dickinson2026,
  author = {Dickinson, Jeff},
  title = {Why We Still Need \{Admiral\} in an Age of {AI}},
  date = {2026-06-14},
  url = {https://pharmaverse.github.io/blog/posts/2026-06-14_admiral_in_age_of_ai/admiral_in_age_of_ai.html},
  langid = {en}
}

A reasonably straightforward post today. I wanted to look at monthly tourism numbers in Samoa. In fact I started to do this for Pacific islands in general, but the data wrangling challenges were sufficient that I only got as far as Samoa for now. There’s interest in these at the moment because they would be a relatively timely indicator of possible economic damage from the fuel crisis related to the Iran war. Visitor numbers, inflation, and merchandise trade are part of the very select number of monthly published economic statistics in this part of the world (for a select subset of countries).

There’s two things happening in this post:

• getting hold of the Samoa’s visitor arrivals numbers; and
• seasonally adjusting them.

The latter is more interesting to me than the former, but unfortunately the former took most of the time.

Data wrangling

Here’s what the visitor numbers look like, once we’ve got them all into one data frame:

The Samoa Bureau of Statistics has a nice Excel workbook up to May 2023, but from that date onwards the data is only available as far as I can see in PDF reports. Luckily these are all available as links from a single page, but there seems to be either a skill issue on my part or some kind of block on the website that stops any systematic download of them all, so I had to download all the PDFs by hand one at a time.

Then Claude helped me write a parser to find the visitor arrivals number in each PDF. Actually, Claude wasn’t any good at actually finding the right number, but it did give me a pattern I could adapt, much as in the old days I would have used Stack Overflow. There are a lot of numbers in each PDF and we need the right one—visitor arrivals, not total arrivals (yes, that’s an em dash, but I write them all by hand). In this case it turns out that the trick is that the PDFs all include the sentence “Overall visitor numbers for the month under review stood at [number]”, and luckily there is no other use of the word “stood” in the document.

The other fiddly thing was extracting the actual date each PDF referred to. Then everything needed to be tested. In the end, it would have been quicker to just manually type the 35 numbers I needed. But here’s the code that does all the data wrangling.

library(tidyverse)
library(rvest)
library(httr2)
library(lubridate)
library(pdftools)   
library(readxl)
library(seasonal)
library(tsibble)
library(fable)
library(feasts)
library(ggtext)
library(scales)

#' Extract date from messy filenames
extract_date <- function(messy_date){
  tibble(path = messy_date) |>
  mutate(
    stem      = tools::file_path_sans_ext(basename(path)),
    # Remove any leading word followed by _ or - (e.g. "Migration_")
    stem      = str_remove(stem, "^([A-Za-z]+[_-])+(?=[A-Z])"),
    month_str = str_extract(stem, "[A-Za-z]{3,}"),
    # Normalise non-standard abbreviations
    month_str = str_replace(month_str, "^Sept$", "Sep"),
    year_str  = str_extract(stem, "\\d{4}|\\d{2}"),
    year      = if_else(nchar(year_str) == 2,
                        as.integer(paste0("20", year_str)),
                        as.integer(year_str)),
    date      = as.Date(paste(year, month_str, "01"), format = "%Y %B %d") |>
                  coalesce(as.Date(paste(year, month_str, "01"), format = "%Y %b %d"))
  ) |>
  pull(date)
}

test <- c("samoa_pdfs/April_25.pdf", "samoa_pdfs/Feb_25.pdf", "samoa_pdfs/Feb_26.pdf", 
"samoa_pdfs/Jan_26.pdf", "samoa_pdfs/January_25.pdf", "samoa_pdfs/July_25.pdf", 
"samoa_pdfs/June-25.pdf", "samoa_pdfs/March_2025.pdf", "samoa_pdfs/March_2026.pdf", 
"samoa_pdfs/May_25.pdf", "samoa_pdfs/Migration_April-2026.pdf",
"samoa_pdfs/Sept-24.pdf", "samoa_pdfs/Migration_Rep_June_2023.pdf")

extract_date(test)

#-------------------PDFs for recent data------------------
# For the more recent years no Excel tables are published, so need
# to use the PDFs and extract total from there
# These had to be downloaded by hand - nothing I tried was able to automate
# that. Download from https://www.sbs.gov.ws/migration/ and save
# in a subfolder /samoa_pdfs/.

pdf_dir <- "samoa_pdfs"
tbl <- tibble(local_path = list.files(pdf_dir, pattern = ".pdf$", full.names = TRUE))

parse_pdf_visitors <- function(path) {
  txt <- tryCatch(pdf_text(path), error = function(e) {
    message("  [WARN] pdftools could not read: ", basename(path))
    NULL
  })
  if (is.null(txt)) return(NA_integer_)
  full_text <- paste(txt, collapse = "\n")
  patterns <- c(
    "stood at [^0-9(]{0,10}\\(?([0-9,]+)\\)?"
  )
  for (pat in patterns) {
    m <- str_match(full_text, pat)[, 2]
    if (!is.na(m)) return(as.integer(str_remove_all(m, ",")))
  }
  message("  [WARN] Could not parse visitor count from: ", basename(path))
  NA_integer_
}

found <- tbl |> 
  filter(file.exists(local_path))

message("  Parsing ", nrow(found), " local PDFs...")
pdf_tbl <- found |>
  mutate(visitors = map_int(local_path, \(p) {
    message("  ", basename(p))
    parse_pdf_visitors(p)
  })) |>
  filter(!is.na(visitors)) |>
  mutate(date = extract_date(local_path))

#-----------------Excel versions for older data------------
# For May 2023 and earlier we can get the data for multiple months
# at a time from Table 1 of the Excel tables. The May 2023 Excel
# file goes back to 2017 January (although the rows are hidden)

fn <- "May_23.xlsx"
if(!file.exists(fn)){
  download.file("https://sbs.gov.ws/images/sbs-documents/social/Arrival/2023/May_23.xlsx",
                destfile = fn, mode = "wb")
}
x <- read_excel(fn, sheet = "Table 1", 
                 range = "D48:D130",
                 col_names = "visitors") |> 
  drop_na() |> 
  pull(visitors)

historical <- tibble(visitors = x, 
               date = seq(as.Date("2017-01-01"), as.Date("2023-05-01"), by = "month"))

#-----------combine and test-----------------------
samoa_visitors <- pdf_tbl |> 
  select(date, visitors) |>
  bind_rows(historical) |> 
  arrange(date) |> 
  mutate(date_month = yearmonth(date)) |> 
  as_tsibble(index = date_month)


# Test - some hand picked test cases, 4 from PDFs and 3 from the Excel
samoa_test <- tribble(~date, ~correct_visitors,
                      "2023-08-01", 16471,
                      "2024-04-01", 12644,
                      "2024-08-01", 17248,
                      "2026-04-01", 14188,
                      "2018-02-01", 7413,
                      "2020-12-01", 195,
                      "2022-06-01", 866) |> 
  mutate(date = as.Date(date))

stopifnot(
  samoa_test |> 
    anti_join(samoa_visitors, by = c("date", "correct_visitors" = "visitors")) |> 
    nrow() == 0
)

And here’s the code that draws the basic time series chart I used earlier:

the_caption = "Source: Samoa Bureau of Statistics"

ggplot(samoa_visitors, aes(x = date, y = visitors)) +
  geom_line() +
  scale_y_continuous(label = comma) +
  labs(x = "",
       y = "Visitor arrivals",
       title = "Visitor arrivals per month to Samoa",
       subtitle = "Unadjusted originals",
       caption = the_caption)

Modelling seasonal decomposition

Phew. Ok, on to the more fun part of actually modelling. I intend to use the X-13ARIMA-SEATS tool. X‑13ARIMA‑SEATS is a program developed and maintained by the US Census Bureau for seasonal adjustment and time‑series decomposition. It will automatically fit a SARIMA (seasonal autoregressive integrated moving average) time series model, has built in methods for identifying and dealing with outliers, by default adjusts for number of trading days and the moving Easter holiday, and allows the user to specify additional regression explanatory variables if you want. It’s the go-to across the world for seasonal adjustment of official statistics.

X-13ARIMA-SEATS is available in R through the seasonal package (by Christoph Sax and Dirk Eddelbuettel). Downstream of that, the feasts and fable packages by (Mitchell O’Hara-Wild, Rob Hyndman and Earo Wangmake) make it easier to work with in a tabular, tidyverse approach.

The Covid period is an obvious dominating factor in tourism anywhere in recent decades, and shows up in the first chart I showed above. I’m also interested in the period since the USA-Israel-Iran war began to see if there is an impact from that. There’s only two months of data (March and April 2026) since the war started, so an impact would have to be dramatic to show up, but it’s worth checking.

X-13ARIMA-SEATS by default will fit forecasts when you ask it to model so you need to provide values of x regressor variables to cover not only the period of the data but a few periods out ahead. I’m doing this as simple time series vectors—I haven’t yet worked out the easy way to do this in the tabular world of fable. Luckily, it seems I can use these vectors later even in fable. In a moment I’ll use both seasonal directly and via fable to make sure I get the same results. In fact, updating my dated knowledge of time series modelling to use fable was one of my main motivations for this whole exercise.

Here’s the code that makes the x regressors I’ll be using for the presence of the Covid pandemic and for the Iran war:

#----------------x regressor variables-----------------

# Covid time series indicator to use as a regressor
covid_reg <- ts(
  as.numeric(seq(as.Date("2017-01-01"), as.Date("2029-04-01"), by = "month") %in%
    seq(as.Date("2020-04-01"), as.Date("2022-07-01"), by = "month")),
  start     = c(2017, 1),
  frequency = 12
)

# Iran war time series indicator
war_reg <- ts(
  as.numeric(seq(as.Date("2017-01-01"), as.Date("2029-04-01"), by = "month") %in%
    seq(as.Date("2026-03-01"), as.Date("2026-07-01"), by = "month")),
  start     = c(2017, 1),
  frequency = 12
)

Directly with seasonal

Ok, it’s modelling time. First, using just an old fashioned time series vector of Samoa’s visitor arrivals and the seasonal package directly, here’s fitting the X-13ARIMA-SEATS model with defaults using both the Covid and Iran war regressors:

sa_ts <- ts(samoa_visitors$visitors, frequency = 12, start = c(2017, 1))

fit_ts_war <- seas(sa_ts, xreg = cbind(covid_reg, war_reg))
summary(fit_ts_war)

Super simple. That gets us this output:

Call:
seas(x = sa_ts, xreg = cbind(covid_reg, war_reg))

Coefficients:
                  Estimate Std. Error z value Pr(>|z|)    
xreg1             -3.10069    0.14567 -21.286  < 2e-16 ***
xreg2             -0.15937    0.13098  -1.217    0.224    
LS2020.Mar        -0.85036    0.14567  -5.838 5.29e-09 ***
AO2020.Dec        -0.75529    0.12350  -6.115 9.63e-10 ***
LS2021.May        -0.79012    0.13233  -5.971 2.36e-09 ***
AO2021.Jul        -1.36973    0.12378 -11.066  < 2e-16 ***
LS2022.May         0.89068    0.13233   6.731 1.68e-11 ***
LS2022.Aug        -1.07689    0.19608  -5.492 3.97e-08 ***
AR-Nonseasonal-01 -0.59006    0.07041  -8.380  < 2e-16 ***
MA-Seasonal-12     0.99961    0.08044  12.427  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

SEATS adj.  ARIMA: (1 1 0)(0 1 1)  Obs.: 112  Transform: log
AICc:  1608, BIC:  1634  QS (no seasonality in final):3.295  
Box-Ljung (no autocorr.): 26.42   Shapiro (normality): 0.9721 *

Some key things to note here.

The data was log-transformed, which is good—I would certainly have chosen to do this, or at least a square root transform, given the way visitor arrivals variance increases as its mean does, all around the world.

The model eventually adopted is described as ARIMA (1 1 0)(0 1 1). This means the main series as a autoregression term on one lag, after one round of differencing; and the seasonal part has a moving average term on on lag, after one round of differencing. This is a very normal model for tourism numbers. it indicates a general trend/drive (the differencing in the main series), a tendency for the values in one month to be related one way or another to those in the month before (in this case with a negative correlation of -0.59) and a strong annual seasonality effect that changes slowly over time.

The Covid dummy variable (xreg1) is strongly negatively significant. With a coefficient of -3.1 and the log transform of the response variable, this means that during the Covid period the actual arrivals were on average exp(-3.1) = 0.045 (ie 4.5% or down 95.5%) of during the non-Covid periods.

In contrast, we don’t have a statistically significant effect for the Iran war (xreg2). For subsequent analysis I will take out that x regressor as I don’t want it complicating the recent trend and seasonally adjusted figures.

Six months’ data have been singled out as outliers and controlled for appropriately. These are all in the difficult-to-model Covid period of 2020 to 2022.

One point to note is that we don’t have an Easter effect. I am 100% sure there is really an Easter effect in Samoa’s visitor arrivals, but 9 years of data isn’t enough to show it. Easter is sometimes in April and sometimes in March. But since 2017, it happens to have been in April every year apart from 2024. That’s just not enough variation to distinguish it from regular monthly seasonal impacts.

To check my recollection that Easter is indeed checked for by default in X-13ARIMA-SEATS, I fit a model to the well known Box and Jenkins airline data:

> # Another examplefor comparison
+ m <- seas(AirPassengers)
+ summary(m)

Call:
seas(x = AirPassengers)

Coefficients:
                    Estimate Std. Error z value Pr(>|z|)    
Weekday           -0.0029497  0.0005232  -5.638 1.72e-08 ***
Easter[1]          0.0177674  0.0071580   2.482   0.0131 *  
AO1951.May         0.1001558  0.0204387   4.900 9.57e-07 ***
MA-Nonseasonal-01  0.1156204  0.0858588   1.347   0.1781    
MA-Seasonal-12     0.4973600  0.0774677   6.420 1.36e-10 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

SEATS adj.  ARIMA: (0 1 1)(0 1 1)  Obs.: 144  Transform: log
AICc: 947.3, BIC: 963.9  QS (no seasonality in final):    0  
Box-Ljung (no autocorr.): 26.65   Shapiro (normality): 0.9908 

Here we see that the number of week days in a month and the moving Easter holiday are indeed in the final model, without me having to have asked for them to be checked. Easter has a positive impact on this air passengers series (1949 to 1960), and the number of week days in a month has a smaller negative impact.

So after that Easter-checking digression, I refit the model for my Samoa visitor arrivals without the war regressor and get an essentially identical result:

> fit_ts <- seas(sa_ts, xreg = covid_reg)
> summary(fit_ts)

Call:
seas(x = sa_ts, xreg = covid_reg)

Coefficients:
                  Estimate Std. Error z value Pr(>|z|)    
xreg              -3.10045    0.14690 -21.106  < 2e-16 ***
LS2020.Mar        -0.83292    0.14617  -5.698 1.21e-08 ***
AO2020.Dec        -0.75463    0.12449  -6.062 1.35e-09 ***
LS2021.May        -0.78975    0.13356  -5.913 3.36e-09 ***
AO2021.Jul        -1.36737    0.12476 -10.960  < 2e-16 ***
LS2022.May         0.89113    0.13356   6.672 2.52e-11 ***
LS2022.Aug        -1.07733    0.19782  -5.446 5.15e-08 ***
AR-Nonseasonal-01 -0.58577    0.07071  -8.284  < 2e-16 ***
MA-Seasonal-12     0.99912    0.07827  12.765  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

SEATS adj.  ARIMA: (1 1 0)(0 1 1)  Obs.: 112  Transform: log
AICc:  1607, BIC:  1631  QS (no seasonality in final):2.501  
Box-Ljung (no autocorr.): 26.63   Shapiro (normality): 0.9706 *

With fable and feasts

tsibble, fable and feasts are a brilliant set of packages that let you work with time series in R in a more tabular and tidyverse-friendly way than the various older time series data structures let you. Most of my work with time series was before they were available, though, so I’m lacking confidence in how they work. Luckily it seems to be pretty straightforward.

I’d already done the critical step earlier with as_tsibble(index = date_month), specifying my main samoa_visitors tibble is actually a time series tibble. Now the modelling using that tsibble is pretty simple:

#----------fable version-------
fit_fb <- samoa_visitors |> 
  model(X_13ARIMA_SEATS(
    visitors ~ xreg(covid_reg)
  ))

report(fit_fb)

Note we use report() rather than summary() to get the end result. I’m not going to print it here because it is literally identical to what we got earier with summary(fit_ts).

The main appeal for me in using the fable/feasts approach is that it fits better with both my data wrangling workflow and my approach to ggplot2 graphics. so here is a nice decomposition of hte original timeseries, produced with autoplot(

fit_fb |> 
  components() |> 
  autoplot() 

Note that in this decomposition, the original ‘visitors’ series and the trend are on the original scale, but the ‘seasonal’ and ‘irregular’ components are expressed as multipliers. So a seasonal value of 1.2 means in a given month the value is 20% higher as a result of the seasonality than otherwise.

And here is my final, presentation version of the data:

Produced with this code:

comp_data <- fit_fb |> 
  components() |> 
  # the 'trend' that comes straight from the decomposition does
  # not adjust ofr the Covid coefficient and looks pretty weird
  # so it is more intuitive to present it after adjustment for
  # Covid, which we need to calculate by multiplying (because of
  # the log transform that SEATS used autoamtically):
  mutate(covid = as.numeric(covid_reg)[1:nrow(samoa_visitors)],
         trend_adj = trend * exp(covid * coef(fit_ts)[1])) 

comp_data|> 
  ggplot(aes(x = date_month, y = season_adjust)) +
  geom_line(linewidth = 1.3) +
  geom_line(aes(y = trend_adj), colour = "steelblue", alpha = 0.9, linewidth = 1.2) + 
  geom_point(aes(y = visitors), colour = "grey70") +
  scale_y_continuous(label = comma) +
  labs(x = "",
       y ="Visitor arrivals",
       title = "Visitor arrivals per month to Samoa",
       subtitle = "<span style='color:grey60'>Original</span>, seasonally adjusted <span style='color:steelblue'>and trend (adjusted for Covid period).</span>",
       caption = "Source: data from Samoa Bureau of Statistics. Seasonal adjustment by freerangestats.info.") +
  theme(plot.subtitle = element_markdown())

After all that, what insight do we have? Well, not a lot really, other than the blindingly obvious trends of the devastation to the industry of Covid and the slow and noisy growth trend since then. We are at least well placed to make commentary on the impact of the fuel crisis on tourism, and can say that there isn’t any evident yet. If and when we do see the impact, we’ll be able to talk about in terms of trends and random variation, after having removed the seasonal element. So that’s useful.

Well, that’s all. Perhaps I’ll get around in a later post to adding the other Pacific island countries with monthly tourism data—Fiji, Vanuatu, Cook Islands, French Polynesia being the key ones I’m aware of.

If you are reading this post on R-bloggers, you will probably know that I have been publishing my selection of the “Top 40” new R packages on CRAN for quite some time. I did this first as part of my work at Revolution Analytics, then on R Views for RStudio and Posit, and now here on R Works. It used to take about a day’s worth of pleasurable work spread out over a month to select forty interesting packages. For a hundred or so packages, I could look at all of the package webpages, download and play with a small number of them. Now, the “Top 40” has become a real hamster-on-the-wheel project. The following plot shows my count of the number of new packages to make it to CRAN since I began publishing on R Works.

library(tidyverse)
file_path <- "new-cran-pkgs.csv"

if (!file.exists(file_path)) {
  stop(paste("File not found! Please check the path:", file_path))
}

# Read text and numbers safely
raw_data <- read.csv(file_path, colClasses = c("character", "numeric"), stringsAsFactors = FALSE)

plot_data <- raw_data |> mutate(Date = my(Month)) |>
  arrange(Date)

new_pkg <- ggplot(plot_data, aes(x = Date, y = Num_pkgs, group = 1)) +
  geom_line(color = "#1f77b4", size = 1) +
  geom_point(color = "#d62728", size = 1.2) +

  labs(
    title = "Monthly Volume of New CRAN Packages",
    x = "Date",
    y = "Number of Packages",
    caption = "Source: R Works monthly Top 40 posts"
  ) +

  theme(
    plot.title = element_text(face = "bold"),
    panel.grid.minor = element_blank(),
    axis.text.x = element_text(angle = 45, hjust = 1)
  )

new_pkg

Why the sharp increase? What could possibly be going on? Well, I have a guess. It is apparently now just too easy to package up some code and ship it to CRAN. That’s understandable: it’s just too easy to write and deploy any kind of software. The following plot that John Burn-Murdoch recently published in the Financial Times based on NBER research shows the explosion of apps in our Agentic AI era.

The plot also indicates that the new apps can’t be making much of a positive contribution to people’s lives or corporate profits. They are apparently not being used, or reviewed, or even discovered.

So let’s ask the same kinds of questions about new R packages. Are most of them really making a contribution to R and the R Community? Are they contributing new statistical methods, extending the reach of R into new application areas, offering efficient high-performance code, or doing something else that would arguably benefit the R Community?

My impression as an engaged dilettante is no: most new R packages are not making a contribution. One obvious indicator of quality is documentation. A significant number of new R packages don’t provide sufficient documentation to explain what they offer. For example, in May, 40 of the 323 new CRAN packages had no README file, no vignettes, and no URL linking to a repository. To my way of thinking, with the possible exceptions of packages that have some sort of discoverable out-of-band documentation (e.g., a journal publication) or are not meant to be called by end users because they are infrastructure for some suite of packages, packages that don’t describe what, why, and how they work are not contributions.

As a hamster on the wheel, I would be very pleased to hear what you have to say. If you are motivated, please leave a comment at Issue #68 on the R Works GitHub repository.

You might call them all “mock”.

Mock the database. Mock the API. Mock the function. The word becomes a catch-all for any test double, any object you substitute for a real dependency in a test. Lumping them together makes it harder to choose the right tool, and the wrong choice leads to brittle, misleading tests.

There are five distinct types, each with a specific job. Knowing which is which is how you stop writing tests that do the wrong thing.

The code under test

All five examples use a single function: process_payment. It charges a card, logs the attempt, and optionally notifies the customer.

process_payment <- function(order, payment_gateway, logger, notifier = NULL) {
  logger$log(paste("Processing order", order$id))

  result <- payment_gateway$charge(order$amount, order$card_token)

  if (!result$success) stop("Payment failed: ", result$error)

  if (!is.null(notifier)) {
    notifier$send(order$customer_id, result$transaction_id)
  }

  result$transaction_id
}

It has three dependencies: payment_gateway, logger, and notifier. Each one will be replaced with a different kind of double depending on what we’re trying to test.

1. Dummy

Definition: an object passed to satisfy a required parameter but never actually used by the test.

process_payment always calls logger$log. The logger is required. But for a test that’s only checking whether the correct transaction ID is returned, we don’t care what gets logged. We just need something that won’t blow up when called.

test_that("returns the transaction ID on successful payment", {
  # Arrange
  order <- list(
    id = "ord-1",
    amount = 100,
    card_token = "tok_visa",
    customer_id = "cust-42"
  )
  dummy_logger <- list(log = function(...) invisible(NULL))
  stub_gateway <- list(
    charge = function(amount, token) {
      list(success = TRUE, transaction_id = "txn-abc")
    }
  )

  # Act
  result <- process_payment(
    order,
    payment_gateway = stub_gateway,
    logger = dummy_logger
  )

  # Assert
  expect_equal(result, "txn-abc")
})
Test passed with 1 success 🥇.

dummy_logger accepts any call and does nothing. The test doesn’t assert on it at all. Its only job is to satisfy the function signature.

A dummy should be the simplest thing that compiles. Recording calls or setting expectations would make it something else. If you find yourself writing a dummy that crashes, or does something unexpected when called, the code path you’re testing actually does use the dependency.

Worth knowing.

2. Stub

Definition: a replacement that returns pre-programmed responses, used to control what the code under test receives.

A stub lets you put the system in a specific state without involving real infrastructure. If you want to test what process_payment does when a card is declined, you don’t need a real payment API. You just return the response you want.

test_that("throws an error when payment is declined", {
  # Arrange
  order <- list(
    id = "ord-2",
    amount = 200,
    card_token = "tok_declined",
    customer_id = "cust-7"
  )
  dummy_logger <- list(log = function(...) invisible(NULL))
  stub_gateway <- list(
    charge = function(amount, token) {
      list(success = FALSE, error = "insufficient funds")
    }
  )

  # Act & Assert
  expect_error(
    process_payment(
      order,
      payment_gateway = stub_gateway,
      logger = dummy_logger
    ),
    "insufficient funds"
  )
})
Test passed with 1 success 🥇.

The stub provides inputs to the system under test. You assert on what the code did with those inputs (in this case, that it threw the right error).

Notice that process_payment accepts payment_gateway as an argument. That’s dependency injection: the function doesn’t create or import its own gateway, so the test can pass in anything with the same interface. Without it you’d need a patching library to intercept the real dependency mid-call. With it, a plain list with a charge function is enough. Stubs work best when the code is designed this way: dependencies accepted as arguments, not hardwired inside.

If you practice test-first development, you’ll notice that you use this pattern all the time. You can’t write the test without it. You don’t know what to patch in a function that doesn’t exist yet! It’s only natural to inject all dependencies as you write the interface of your code.

When the dependency isn’t declared in the interface, when the function calls another function directly by name, mockery::stub() can patch it for the duration of a test:

# A function that calls charge_card() internally, with no way to inject it
process_payment_legacy <- function(order) {
  result <- charge_card(order$amount, order$card_token)
  if (!result$success) {
    stop("Payment failed: ", result$error)
  }
  result$transaction_id
}

charge_card <- function(amount, token) {
  stop("would call real payment API")
}

test_that("returns transaction ID when charge succeeds", {
  # Arrange
  order <- list(amount = 100, card_token = "tok_visa")
  mockery::stub(
    process_payment_legacy,
    "charge_card",
    function(amount, token) {
      list(success = TRUE, transaction_id = "txn-stub")
    }
  )

  # Act
  result <- process_payment_legacy(order)

  # Assert
  expect_equal(result, "txn-stub")
})
Test passed with 1 success 🥳.

mockery::stub() replaces charge_card inside the scope of process_payment_legacy for that one test call, without touching the real function anywhere else.

mockery::stub() has a catch. The stub is targeted by function name as a string, so if you rename charge_card, the stub silently stops working and the test passes against the real function with no warning. The test is also coupled to an implementation detail: if you refactor process_payment_legacy to call payment_gateway$charge() instead, the stub breaks even if the behavior is unchanged. That’s the Overspecification smell.

Use mockery::stub() when you’re working with legacy code that wasn’t built with testability in mind and you can’t refactor the interface right now. It lets you get tests in place quickly. Treat it as a stepping stone: once the characterization tests are green, refactor toward dependency injection and replace the patch with a plain stub passed as an argument.

To sum up: when you need to control what a dependency returns and don’t care how it was called, reach for a stub.

3. Spy

Definition: a stub that also records calls made to it, so you can assert on them afterward.

Sometimes the behavior you’re testing is a side effect. A notification that should have been sent, a message that should have been logged. The code doesn’t return a value you can assert on. It calls something. A spy captures those calls.

make_notifier_spy <- function() {
  calls <- list()
  list(
    send = function(customer_id, transaction_id) {
      calls[[length(calls) + 1]] <<- list(
        customer_id    = customer_id,
        transaction_id = transaction_id
      )
    },
    calls = function() calls
  )
}
test_that("notifies the customer after successful payment", {
  # Arrange
  order <- list(
    id = "ord-3",
    amount = 50,
    card_token = "tok_visa",
    customer_id = "cust-99"
  )
  dummy_logger <- list(log = function(...) invisible(NULL))
  stub_gateway <- list(
    charge = function(amount, token) {
      list(success = TRUE, transaction_id = "txn-xyz")
    }
  )
  spy_notifier <- make_notifier_spy()

  # Act
  process_payment(
    order,
    payment_gateway = stub_gateway,
    logger = dummy_logger,
    notifier = spy_notifier
  )

  # Assert
  expect_length(spy_notifier$calls(), 1)
  expect_equal(spy_notifier$calls()[[1]]$customer_id, "cust-99")
  expect_equal(spy_notifier$calls()[[1]]$transaction_id, "txn-xyz")
})
Test passed with 3 successes 🌈.

The spy is a stub with memory. You call the code, then interrogate the spy to see what happened.

You don’t always need to build a spy by hand. mockery::mock() also collects calls, so it can serve as a spy when you want the recording behaviour without writing the closure yourself:

test_that("notifies the customer after successful payment (mockery spy)", {
  # Arrange
  order <- list(
    id = "ord-3b",
    amount = 50,
    card_token = "tok_visa",
    customer_id = "cust-99"
  )
  dummy_logger <- list(log = function(...) invisible(NULL))
  stub_gateway <- list(
    charge = function(amount, token) {
      list(success = TRUE, transaction_id = "txn-xyz")
    }
  )
  spy_send <- mockery::mock()

  # Act
  process_payment(
    order,
    payment_gateway = stub_gateway,
    logger = dummy_logger,
    notifier = list(send = spy_send)
  )

  # Assert
  mockery::expect_called(spy_send, 1)
  expect_equal(mockery::mock_args(spy_send)[[1]][[1]], "cust-99")
  expect_equal(mockery::mock_args(spy_send)[[1]][[2]], "txn-xyz")
})
Test passed with 3 successes 😀.

The handwritten version is clearer when you want the recording mechanism visible to readers, useful in a codebase where not everyone knows mockery. mockery::mock() is more concise once the team is familiar with the library.

The difference from a mock comes down to return values. A spy records calls and nothing else. A mock records calls and can also return pre-programmed values, which makes it useful when you need the dependency to behave a specific way and want to assert on how it was used.

4. Mock

Definition: “a double pre-programmed with expectations that form a specification of the calls it should receive. A true mock can throw if it receives a call it doesn’t expect, and is checked during verification to confirm it got all the calls it was expecting.”^{[1, Fowler]}

mockery::mock() is looser than that definition. It accepts any call without complaining and doesn’t enforce expectations upfront. It records every call it receives (the arguments, the order, the count) and returns pre-programmed values you supply. Verification is your responsibility in the Assert step.

test_that("sends exactly one notification with correct arguments", {
  # Arrange
  order <- list(
    id = "ord-4",
    amount = 75,
    card_token = "tok_visa",
    customer_id = "cust-11"
  )
  dummy_logger <- list(log = function(...) invisible(NULL))
  stub_gateway <- list(
    charge = function(amount, token) {
      list(success = TRUE, transaction_id = "txn-def")
    }
  )
  mock_notifier <- list(send = mockery::mock())

  # Act
  process_payment(
    order,
    payment_gateway = stub_gateway,
    logger = dummy_logger,
    notifier = mock_notifier
  )

  # Assert
  mockery::expect_called(mock_notifier$send, 1)
  mockery::expect_args(mock_notifier$send, 1, "cust-11", "txn-def")
})
Test passed with 5 successes 🥳.

Use a mock when the interaction itself is what you’re testing: whether the code called the dependency in the right way, with the right arguments.

They’re also the easiest double to overuse. Assert on every call to every dependency and you’ve written an overspecified test, one that breaks whenever the implementation changes even when the behavior stays the same.

Prefer a spy when you only need to record calls. A plain list with a function that appends to a vector is often enough. Reach for a mock when you also need to control what the dependency returns. The risk is the same with any interaction-based assertion: check every call to every dependency and you end up with a test that mirrors the implementation rather than the behaviour, breaking whenever the internals change even when the outcome doesn’t.

5. Fake

Definition: a working implementation that’s simpler than the real thing, suitable for tests but not production.

A fake isn’t just a pre-programmed response. It has real behavior. An in-memory database is a fake: it stores and retrieves data like the real thing, just without persistence or network overhead. It behaves correctly across multiple calls, which a stub can’t do.

make_fake_payment_gateway <- function() {
  transactions <- list()

  list(
    charge = function(amount, token) {
      if (amount <= 0) {
        return(list(success = FALSE, error = "invalid amount"))
      }
      if (token == "tok_declined") {
        return(list(success = FALSE, error = "card declined"))
      }

      id <- paste0("txn-", length(transactions) + 1)
      transactions[[id]] <<- list(
        amount = amount,
        token = token
      )
      list(success = TRUE, transaction_id = id)
    },
    find = function(transaction_id) {
      transactions[[transaction_id]]
    }
  )
}
test_that("successful charges are recorded in the gateway", {
  # Arrange
  order <- list(
    id = "ord-5",
    amount = 120,
    card_token = "tok_visa",
    customer_id = "cust-3"
  )
  dummy_logger <- list(log = function(...) invisible(NULL))
  fake_gateway <- make_fake_payment_gateway()

  # Act
  txn_id <- process_payment(
    order,
    payment_gateway = fake_gateway,
    logger = dummy_logger
  )

  # Assert
  recorded <- fake_gateway$find(txn_id)
  expect_equal(recorded$amount, 120)
  expect_equal(recorded$token, "tok_visa")
})
Test passed with 2 successes 🎊.

Fakes work well when you need to test behaviour across multiple operations: place an order, query its status, refund it. A stub would need to be reprogrammed for each call. A fake just handles it.

They’re also a good fit for acceptance tests and manual inspection. An acceptance test exercises a full user-facing behaviour end-to-end, several layers of the application working together. At that level you don’t want stubs reprogrammed for individual calls; you want a dependency that behaves realistically across the whole flow. A fake payment gateway, a fake email sender, a fake file store: these let your acceptance test suite run in CI without connecting to external services, needing credentials, or leaving side effects behind. You can also wire the same fakes into a development mode of the app. Spin up the Shiny app pointing at the in-memory gateway and you can click through every payment scenario without touching a real API.

The cost is that fakes take time to build and maintain. They need to be kept in sync with the real interface they’re replacing. For a small, stable interface that’s used heavily across your test suite and in manual workflows, the investment pays off. For a dependency you only use in one unit test, a stub is simpler.

When to reach for each one

Double	Has behaviour	Records calls	Returns programmed values	Rejects unexpected calls	When to use
Dummy					Fill a required parameter you won’t touch
Stub	Pre-programmed only				Control what the code receives
Spy	Pre-programmed only				Assert on side effects after the fact
Mock (mockery)	Pre-programmed only				Assert on calls and control what the code receives
Mock	Pre-programmed only				Pin an exact interaction as a hard contract
Fake					Replace stateful or multi-call dependencies

The key difference is between stub and mock. A stub returns values. You assert on the outcome. A mock records calls and can return pre-programmed values. Using a mock where a stub would do couples your test to implementation details. Using a stub where a mock is needed means missing the interaction you were trying to verify.

When in doubt: if you’re asserting on a return value or a state change, use a stub. If you’re asserting that a specific call was made, use a spy or a mock. If the dependency has real state that needs to survive across calls, build a fake.

Appendix: implementing an eager mock by hand

mockery::mock() is sufficient for everyday use. Skip this if you’re not curious about mock that throws failures during execution of code under test.

This is what a mock matching Fowler’s definition looks like in R. It takes a list of expected calls in the Arrange step, fails immediately on anything unexpected, and exposes a verify() function to confirm every expected call was made.

make_mock_notifier <- function(expected_calls) {
  received <- list()

  list(
    send = function(customer_id, transaction_id) {
      call <- list(
        customer_id = customer_id,
        transaction_id = transaction_id
      )
      match <- any(sapply(expected_calls, identical, call))
      if (!match) {
        testthat::fail(sprintf(
          "Unexpected call: send('%s', '%s')",
          customer_id,
          transaction_id
        ))
      }
      received[[length(received) + 1]] <<- call
    },
    verify = function() {
      for (exp in expected_calls) {
        found <- any(sapply(received, identical, exp))
        if (!found) {
          testthat::fail(sprintf(
            "Expected call never made: send('%s', '%s')",
            exp$customer_id, exp$transaction_id
          ))
        }
      }
      testthat::succeed()
    }
  )
}

The mock rejects unexpected calls on the spot:

test_that("throws immediately when called with unexpected arguments", {
  # Arrange
  order <- list(
    id = "ord-4b",
    amount = 75,
    card_token = "tok_visa",
    customer_id = "cust-11"
  )
  dummy_logger <- list(log = function(...) invisible(NULL))
  stub_gateway <- list(
    charge = function(amount, token) {
      list(success = TRUE, transaction_id = "txn-def")
    }
  )
  mock_notifier <- make_mock_notifier(
    expected_calls = list(list(
      customer_id = "cust-WRONG",
      transaction_id = "txn-def"
    ))
  )

  # Act — throws before we even reach Assert
  process_payment(
    order,
    payment_gateway = stub_gateway,
    logger = dummy_logger,
    notifier = mock_notifier
  )
})
── Failure: throws immediately when called with unexpected arguments ───────────
Unexpected call: send('cust-11', 'txn-def')
Backtrace:
    ▆
 1. └─global process_payment(...)
 2.   └─notifier$send(order$customer_id, result$transaction_id)

Error:
! Test failed with 1 failure and 0 successes.

And verify() catches expected calls that were never made:

test_that("fails verification when an expected call was never made", {
  # Arrange
  order <- list(id = "ord-4c", amount = 75, card_token = "tok_visa", customer_id = "cust-11")
  dummy_logger <- list(log = function(...) invisible(NULL))
  stub_gateway <- list(charge = function(amount, token) list(success = TRUE, transaction_id = "txn-def"))
  mock_notifier <- make_mock_notifier(
    expected_calls = list(
      list(customer_id = "cust-11", transaction_id = "txn-def"),
      list(customer_id = "cust-99", transaction_id = "txn-xyz")  # will never be called
    )
  )

  # Act
  process_payment(order, payment_gateway = stub_gateway, logger = dummy_logger, notifier = mock_notifier)

  # Assert
  mock_notifier$verify()
})
── Failure: fails verification when an expected call was never made ────────────
Expected call never made: send('cust-99', 'txn-xyz')

Error:
! Test failed with 1 failure and 1 success.

The happy path passes both checks:

test_that("passes when all expected calls are made and no unexpected ones occur", {
  # Arrange
  order <- list(
    id = "ord-4d",
    amount = 75,
    card_token = "tok_visa",
    customer_id = "cust-11"
  )
  dummy_logger <- list(log = function(...) invisible(NULL))
  stub_gateway <- list(
    charge = function(amount, token) {
      list(success = TRUE, transaction_id = "txn-def")
    }
  )
  mock_notifier <- make_mock_notifier(
    expected_calls = list(list(
      customer_id = "cust-11",
      transaction_id = "txn-def"
    ))
  )

  # Act
  process_payment(
    order,
    payment_gateway = stub_gateway,
    logger = dummy_logger,
    notifier = mock_notifier
  )

  # Assert
  mock_notifier$verify()
})
Test passed with 1 success 🥇.

References

1. Martin Fowler — TestDouble
2. Gerard Meszaros — Test Double Patterns

What makes rOpenSci’s software peer review work is people who care deeply about the quality and usability of scientific software, and who give their time and expertise to help others build it better. Today, we’re pleased to announce three new members of our editorial team.

We welcome Joel H. Nitta, and Nicholas Tierney as new editors, and formally introduce Ronny Hernández Mora, who joined the editorial team in August 2025. Each brings a distinct perspective shaped by their work, their communities, and their experience with open-source R development. Together, they strengthen our capacity to serve the growing number of package authors who submit their work for review, and to uphold the collaborative, friendly and transparent standards that rOpenSci software peer review is known for.

Ronny Hernández Mora

Ronny is a PhD student at the University of Alberta and a research software developer with a background in data analysis and remote sensing. His current research explores how perception-driven drone systems can support detection, monitoring, and decision making in real world settings.

Before starting his PhD, Ronny worked as a data developer at ixpantia, developing data tools, automation pipelines, APIs, and production analytical applications for organizations across Latin America and the United States. Alongside his PhD, he works with Openscapes and contributes to opensource communities through teaching, mentorship, talks, and collaborative software projects.

Ronny on GitHub, Website.

I first heard about rOpenSci through a Data Latam podcast episode, and from there I started following the organization and its work. Over time, I got to know people connected to the community, but it was not until 2024 that I became directly involved by volunteering as a package reviewer.

I have always liked the idea of building software collaboratively: creating tools that can be understood, reviewed, reused, and improved by others. That is one of the things I value most about rOpenSci: the way it combines technical review with openness, care, and constructive feedback. I am grateful for the opportunity to contribute as an editor and to support authors and reviewers through that process and continue learning through this community effort.

Joel H. Nitta

Joel is currently an associate professor at Chiba University, Japan, where he studies the evolution and ecology of ferns. Throughout his career as a botanist, he has also cultivated a keen interest in reproducible data analysis and has authored several R packages. Two of these are currently part of rOpenSci, canaper for spatial phylogenetic analysis, and dwctaxon for validation and maintenance of taxonomic databases. He is also an official maintainer of two rOpenSci packages, rgnparser for parsing taxonomic names, and restez for querying the GenBank DNA database. Outside of rOpenSci, Joel is active in the Bio”Pack”athon community in Japan, serves on the organizing committee of the Pteridophyte Phylogeny Group, and is a certified instructor for The Carpentries. Joel’s hobbies include various forms of being active outside (hiking, running, backpacking, botanizing, cycling), playing the euphonium, and tabletop games.

Joel on GitHub, Website.

rOpenSci has been incredibly important in my journey with data science and R. Above all, the extremely knowledgeable and helpful rOpenSci community has provided invaluable support as I transitioned from “R package user” to “R package developer”. While it may be possible to do this on your own, it is much more enjoyable and efficient when you are in the company of like-minded people. I especially appreciate this because I did not study computer science in undergrad and I am largely self-taught when it comes to R, having first learned it out of necessity during my graduate studies. I am very excited to officially join the editorial team and have the opportunity to give back to the organization that has supported me so much.

Nicholas Tierney

Nicholas (Nick) Tierney is a statistician, Research Software Engineer, and freelance consultant with a PhD in Statistics who specialises in data analytics, R package development, and teaching. Previously, he worked with Professor Nick Golding at The Kids Research Institute Australia and was a Research Fellow at Monash University with Professor Dianne Cook, where he developed tools for exploratory data analysis including visdat, naniar, and brolgar. Nick actively writes about R related projects at his blog, “credibly curious”. When not coding, Nick enjoys outdoor adventures and hiked the entire Pacific Crest Trail in 2023, documenting his journey at njt.micro.blog.

Nick on GitHub, Website.

I remember seeing rOpenSci online in 2014, at their first rOpenSci hackathon in San Francisco, USA. I was inspired by not just the projects they worked on, but the collective of people who were kind and generous as much as they were brilliant. I helped run an offshoot of the rOpenSci Unconf in Australia, which ran in 2016, 2017, 2018, and 2019. I received amazing mentorship and support from Karthik Ram and Stefanie Butland in learning how to run these unconferences. I think these have had a lasting impact on connecting people from Australia and New Zealand.

I had my R package, visdat, reviewed by rOpenSci in 2017, and this experience was formative to my understanding of software review, and greatly improved my own day to day practice. In 2024, with Eric Scott, and Andrew Brown, we submitted the geotargets package to rOpenSci.

rOpenSci has always been at the forefront of building a community of practice, and their model of peer review is a gold standard. Science would be better if more journals practised such transparency and kindness. It is an enormous honour to be able to give back to a community that has been so formative in my career and my life.

About the Software Peer Review Program

rOpenSci’s software peer review program brings together volunteers to collaboratively review scientific and statistical software according to transparent, constructive, and open standards. We have two different type of review: one for general research software packages and another for package that implement statistical methods. Editors like Ronny, Joel, and Nick are central to that process: they handle initial submission checks, identify and coordinate reviewers, and guide authors through the review until their package is ready.

Get Involved

Thinking about submitting a package? Start here:

• rOpenSci Software Peer Review: scope, process, and guidelines;
• rOpenSci Packages: Development, Maintenance, and Peer Review: the full guide for package authors;
• rOpenSci Statistical Software Peer Review: the full guide for statistical software submissions;
• Public software review threads on GitHub: see software peer review in action.

Would you like to contribute as a reviewer? We’d love to have you. Fill out the rOpenSci Reviewer Sign-Up Form and we’ll match you with packages that fit your expertise.

A warm welcome to Ronny, Joel, and Nick! We’re very glad to have you with us.

Read it in: .

We are very happy to introduce the new rOpenSci Champions. This group will experience the program and work in Spanish, allowing us to continue to strengthen the open science and research software development community in this language. We are excited about the projects they will develop, which address real challenges from different disciplines and territories in Latin America.

We invite you to meet each of these people and the projects they will be working on throughout the program.

Bastián Olea Herrera

My name is Bastián, I live in Chile, I am a sociologist by training and I have a Master’s degree in sociology. I like to create content and tutorials about R, and recently we are organizing a community of R users in Santiago, where I live. I am dedicated to data analysis in the public sector, where I work mainly with social data at different territorial levels. This kind of data always needs the same types of cleaning and corrections, but at the same time they usually vary a lot in small details between each source. That is why I applied to the program to develop a package that facilitates working with data at the communal and regional level in Chile. I hope to meet other R users and learn together, as well as receive support from people with more experience and knowledge, so that we can create useful tools for many people!

Denisse Fierro Arcos

My name is Denisse Fierro Arcos, I am a marine scientist originally from Ecuador but currently living in Australia. At the moment I am close to finishing my PhD program which focuses on developing best practices for the use of oceanographic models in marine ecosystem studies. I am also currently working as a researcher at the University of Tasmania where I am developing a marine ecosystem model that will allow us to project the possible impacts of climate change on marine species and fisheries. It is precisely this model that is the focus of my project with the Champions Program. We want to publish a package in R that will allow other marine researchers to easily run our model.

Durga Valentina Linares Herrera

Hello! My name is Durga Valentina Linares Herrera and I live in Lima, Peru. I am a social scientist, with a degree in Political Science from the Universidad Antonio Ruiz de Montoya and a diploma in Data Science for Social Sciences and Public Management from the Pontificia Universidad Católica del Perú. I work as a research assistant at the Research Center of the Universidad del Pacífico, where I explore the intersection between technology and work from a sociological perspective and with mixed methodologies. My current projects revolve around the Peruvian labor market and its exposure to artificial intelligence, and the evolution of the strike as a social phenomenon in Peru in the last three decades.

For the program I presented as a project {epen}, an R package to download, process, and analyze microdata from the Permanent National Employment Survey and its predecessor, the Lima EPE. This idea was born from my experience working directly with these databases. I could see how difficult it can be to access public labor data in a fast, clear, and reproducible way, especially for those of us who come from the social sciences and do not always come to these tools with a previous technical background, as was also my case. With this package, I seek to facilitate that path for researchers, students, and public sector analysts who need to build reliable labor indicators without having to start from scratch each time, something I find especially relevant in this time of so many changes.

In recent years, with the intention of developing new skills in a context where technology is quickly advancing, I moved towards a more quantitative approach and learned to program in a self-taught way with free resources on the internet. That process has been important for my professional development, but it also made me want to turn that learning into something useful for other people: a tool to help reduce the entry barriers that I myself encountered when I started. That’s how the motivation behind this package came about, and this Champions Program appeared just at the right time. I am excited to participate in an initiative like this and to be part of a network like rOpenSci, whose commitment to a more diverse and inclusive open science seems very valuable to me. I hope that this experience will allow me to consolidate in community a path that I once started on my own and pave the way for future packages aimed at improving access to Peruvian public data in R.

Evelia Lorena Coss Navarrete

I am a Biotechnology Engineer from the Polytechnic University of Sinaloa (UPSIN), originally from Mazatlan, and I did my Master’s and PhD in Plant Biotechnology at Cinvestav, where I studied the conservation of lncRNAs in plants.

My participation in the rOpenSci Champions Program seeks to strengthen communities such as VieRnes de Bioinformatics, R-Ladies Morelia and RSG-Mexico, promoting the creation of R packages with international standards, reproducible documentation and open review. My vision is to bridge the gap between the global rOpenSci community and local initiatives in Latin America, promoting a more open, inclusive and sustainable science.

Gladys Choque Ulloa

Hello! My name is Gladys Choque Ulloa, I am originally from Peru and currently reside in Brazil. I have a degree in Statistics, a Master’s in Statistics, and I am currently pursuing my PhD in Computer Science at the ICMC-USP of the University of São Paulo, Brazil. There, I am developing research in computational neuroscience focused on the automatic diagnosis of mental disorders, using Machine Learning models, Neural Networks, LLMs, Causal Inference, and Time Series. In addition to my academic work, I am founder of the organization Women in DataLab, where I work to reduce the gender gap in technology and data science.

I applied to the rOpenSci Champions Program because I strongly believe in the power of open science and reproducible software to democratize knowledge. During the program, my goal is to hone my skills in developing R tools and scientific software under global standards. With this experience I want to strengthen my technical profile and act as a bridge for more researchers in Latin America to adopt collaborative and open practices, enhancing the impact of our scientific community internationally.

José Daniel Conejeros Pavez

I am José Daniel Conejeros, MSc in Statistics from the Pontificia Universidad Católica de Chile. I am originally from Chile and currently develop my work between Chile and Italy. In Italy I work as a Lagrange Fellow at the ISI Foundation and as a young researcher at the SENTINET Center (Surveillance, epidemiology and new technologies for emerging infectious threats), working on issues of data science, complex systems, and public health.

My work is situated at the intersection of statistics, epidemiology, and computational social science to understand infectious disease dynamics and health outcomes for different populations. Currently, I am developing spatmask, an R package oriented to the masking and anonymization of spatial data (geomasking), with the goal of enabling reproducible analyses without compromising the privacy of individuals. This project aims to bridge the gap between the use of sensitive data for research and the ethical and regulatory restrictions that limit its access and use.

I decided to apply to the Champions Program because many of the challenges I face are not only technical, but also organizational and cultural. I want to understand how to develop scientific software collaboratively, how to document it correctly and how to foster reproducible practices in contexts where open science is not yet fully installed. I am interested in learning how to build tools that not only work, but that are understandable, auditable, and useful for research teams and Public Policy decision makers.

Through this program, I hope to strengthen my skills in scientific software development, open review and collaborative work. My goal is to translate this learning into concrete transfer: training, collaboration, and community building around open science in the region.

Linda Cabrera Orellana

Hello! My name is Linda, I am Ecuadorian and currently reside in Granada, Spain. I work as a Senior Data Analyst in a digital marketing agency and I share my experience as a teacher in business schools, where I teach classes on AI applied to analytics and data visualization.

My project consists of developing an R package with an educational approach to detect, explain, and help correct common structural errors in manually created datasets, specially designed for people with no technical background in data. I am applying to the program because I want to turn my practical and teaching experience into a real contribution to the open scientific software ecosystem, and to strengthen my skills in R package design, documentation, and maintenance.

I also hope that the project will serve as an educational resource in universities in Ecuador and in the R-Ladies community, and that it will build a bridge between those who collect data and those who analyze it, promoting data literacy in Spanish.

María Florencia Tames

My name is María Florencia Tames, I am a professor at the National University of Córdoba (Argentina) and I work in research in the area of air quality, exposure to air pollutants and environmental inequalities in urban contexts in Latin America.

With a team I developed the AirExposure R package, a tool for estimating daily exposure to air pollutants by integrating information on ambient concentrations, mobility, and daily routines. Through the rOpenSci Champions Program, my goal is to improve and consolidate this package as an open, reproducible and accessible tool for the community.

I am especially interested in strengthening my skills in open software development, incorporating best practices such as documentation, testing, and peer review, and learning to build tools that can be used beyond their original research context.

I am also motivated by the program’s focus on community and working in Spanish, as access to programming resources and training is often limited by language. I hope to be able to share what I have learned through training activities and contribute to the development of an open scientific software community in Latin America.

Marina Cecilia Cock

I am from Santa Rosa, La Pampa, Argentina. I have a degree in Natural Resources and Environment Engineering from the National University of La Pampa (UNLPam) and a PhD in Agrarian Sciences with orientation in Ecology from the University of Buenos Aires (UBA).

I am currently working as a research assistant at the National Council of Scientific and Technical Research (CONICET) and as a teaching assistant in Biogeography and Statistics at the National University of La Pampa.

I applied to the program to learn about the R package review process. I am interested in participating because I find it especially valuable to join a community where learning is shared and collaborative. I am looking to strengthen my knowledge in R and then be able to transmit it both in my teaching role and within the R community.

Patricia Andrea Loto

I have a degree in Information Systems and a Diploma in Data Science, Machine Learning and its Applications. I am currently pursuing a Master’s Degree in Information Technology at the Universidad Nacional del Nordeste. I live and work in Argentina.

I work as a software developer and data analyst in the public sector, and as a university teacher in the area of systems and programming. I am a co-founder of RSE Argentina, member of the organizing committee of LatinR, and co-organizer of R-Ladies Resistencia-Corrientes. From these spaces I work to promote open science, reproducible research software and technology inclusion in Latin America.

In the 2026 cohort I am participating as mentee with the guidance of Guadalupe Pascal. My project is an R package to systematize the creation of Software and Data Management Plans through standardized templates developed with Quarto, which facilitate the documentation, preservation, and reuse of data and scientific code under international open science standards.

I applied to the Champions Program because I see it as a bridge: between where I am technically and where I want to be. I hope to learn about the rOpenSci peer review process, deepen my understanding of testing, technical documentation, and everything that makes a package really useful to others. I’m also looking to connect with a network of developers and researchers who share open science values. Finally, I want research software developed from the Global South to have greater visibility and for our contributions to be recognized, and I think rOpenSci is the ideal platform for that.

Estefania Torrejón

I am Peruvian and a biologist graduate from the Universidad Nacional Mayor de San Marcos (UNMSM), with a Master’s degree in Biomedical Sciences from the Institute of Hygiene and Tropical Medicine (IHMT) in Lisbon, Portugal. I currently reside in Portugal, where I work as a predoctoral researcher in bioinformatics at the Metabolic Diseases Research Lab of NOVA Medical School. I am also director of the International Relations Department of the Peruvian Society of Bioinformatics and Computational Biology.

The project I am developing in the framework of the rOpenSci Champions Program consists of the preparation and submission of my R package, EV-Net, to the CRAN repository. EV-Net is a bioinformatics tool that identifies and prioritizes molecules present in the cargo of extracellular vesicles (EVs) with high regulatory potential on a receptor tissue of interest. EVs are structures surrounded by a lipid bilayer that carry a great diversity of active molecules. These vesicles can move through the organism and reach specific tissues, which is why they are recognized as key mediators of cell-to-cell communication (CCC). However, most of the bioinformatics tools available to study CCC do not consider EV-mediated communication. To address this limitation, my collaborators, supervisors, and I developed EV-Net.

Being an rOpenSci Champion is an invaluable opportunity to receive the necessary training, coaching and mentoring to bring EV-Net up to the required quality standards and to be successfully incorporated into CRAN.

Next steps

With the presentation of this new group of Champeons, we begin the fourth edition of the program, the second in Spanish. This group has already started the training stage and met their mentors. They will be working for 12 months developing new packages, preparing existing packages to submit to the peer review process, and reviewing other people’s packages.

If you want to follow the development of their projects and where and when their dissemination and communication activities will take place, don’t miss our blog articles, news in our newsletter and social networks.

Stanislaw Ulam, Los Alamos, 1963 was bored in a meeting and he started dooddling integers in a spiral and circled the primes. Diagonal lines appeared. He later showed it to Martin Gardner, to Ulam surprise, Gardner published his findings in Scientific American. We are still confused to this day.

But Ulam was not doodling some little houses or boats, or cars like a normal person. No. He wrote numbers in a spiral. Then circled all the prime numbers. Then stared at what he had created with the dawning horror of a man who has seen too much.

So where does the spiral comes from? Start with 1 in the middle. Write the number in a spiral outward. And then highlight all the prime numbers. If you draw long enough diagonal lines will appear.

Primary function has initial position and directions calculated, in order to have the symmetry of the plot.

ulam_prime_spiral <- function(
    n         = 51,    
    theme     = c("Cosmic","Blueish","Classy","Psycho"),
    show_nums = FALSE,   # print values (n ≤ 21 only)
    animate   = FALSE,       
    speed     = 1,          
    verbose   = TRUE
) {
  
  theme <- match.arg(theme)
  
  #n must be odd so integers have a unique centre cell 
  if (n %% 2 == 0) { n <- n + 1L }
  if (n < 5) stop("Size muste be > 5.")
  total <- n^2
  mid   <- (n + 1L) / 2L            
  
  dr <- c( 0L, -1L,  0L,  1L)   # row deltas:  E  N  W  S
  dc <- c( 1L,  0L, -1L,  0L)   # col deltas:  E  N  W  S
  
  mat   <- matrix(0L, n, n)     
  ord_r <- integer(total)        
  ord_c <- integer(total)        
  
  r <- mid;  
  cc <- mid  
  
  mat[r, cc] <- 1L
  ord_r[1]   <- r
  ord_c[1]   <- cc
  
  d    <- 1L  # current direction index (1=E 2=N 3=W 4=S)
  step <- 1L  # current arm length
  num  <- 2L                     
  
  while (num <= total) {
    for (half in 1:2) {          # each dir is twice 
      for (i in seq_len(step)) {
        r  <- r  + dr[d]
        cc <- cc + dc[d]
        mat[r, cc] <- num
        ord_r[num] <- r
        ord_c[num] <- cc
        num <- num + 1L
        if (num > total) break  
      }
      d <- (d %% 4L) + 1L       # turn left: E→N→W→S→E
      if (num > total) break    
    }
    step <- step + 1L            
  }
  

and we determine the prime:

# Sieve of Eratosthenes   
  is_prime    <- rep(TRUE, total)
  is_prime[1] <- FALSE
  p <- 2L
  while (p * p <= total) {
    if (is_prime[p])
      is_prime[seq.int(p * p, total, p)] <- FALSE
    p <- p + 1L
  }

  prime_mat <- matrix(is_prime[mat], n, n)
  n_primes  <- sum(prime_mat)
  density   <- 100 * n_primes / total

With for type of visuals (Classy, Psycho, Blueish and Cosmic), if you decide to create a grid smaller than 21×21, you can have also the numbers displayed. And of course, the diagonals are visible as well:

As always, the complete code is available on GitHub in  Useless_R_function repository. The complete version of code is here: https://github.com/tomaztk/Useless_R_functions/blob/main/functions/ulam_spiral.R

Check the repository for future updates!

Stay healthy and happy R-coding!

Back to basics! Learning non-polar amino acids, what zwitterions actually are, and dipping into the applied math — Rodrigues rotation and Lennard-Jones potential. Slowly building toward optimal phi/psi!

Motivations

We’ve explored quite a bit lately in molecular dynamic simulation and then protein-protein docking as well the last time. There is still so much to learn. I’ve decided to go back to basics, revisiting our old friends amino acids and try to understand the natural properties behind each one and see if that will make more sense in the future when we’re exploring more. While making notes for myself of all the amino acids, I’ll also try to understand some of the basic math behind the structures. Are you ready !? Lol, I’m not, but let’s go anyway!

Objectives:

• Amino Acids
  • Non-polar
• Rodriguez Rotation Formula
• Lennard-Jone Potential Enegery
  • Bond Stretch
  • Bond Angle
  • Proper dihedral
  • Non-bonded
  • Calculating LJ
• Opportunities For Improvement
• Lessons Learnt

Amino Acids

Amino acids are the building blocks of proteins, each sharing a common backbone: a central α-carbon bonded to an amino group (–NH₂), a carboxyl group (–COOH), a hydrogen atom, and a variable side chain (R group) that defines each amino acid’s identity and chemistry.

Non-polar amino acids

Non-polar amino acids have hydrophobic side chains — they avoid water and tend to cluster in the interior of folded proteins, forming the hydrophobic core that drives protein stability. Understanding each one’s shape and bulk is directly relevant to how they pack, how they constrain backbone flexibility, and how substitutions affect enzyme active sites.

library(tibble)
library(kableExtra)

aa_nonpolar <- tribble(
  ~aa,  ~aa3,  ~name,           ~functional_group,  ~smiles_sidechain,    ~charge_ph7, ~mw_da, ~pka,     ~md_note,                                                                                ~main_function,
  "G",  "Gly", "Glycine",       "H (none)",         "[H]",                "Neutral",   75.03,  NA_real_, "Minimal VDW radius; unrestricted phi/psi; near-zero excluded volume",                  "Conformational flexibility; tight turns; active site geometry",
  "A",  "Ala", "Alanine",       "Methyl",           "C",                  "Neutral",   89.09,  NA_real_, "Low steric perturbation; high alpha-helix propensity in force fields",                 "Helix former; hydrophobic core; alanine-scanning mutagenesis",
  "V",  "Val", "Valine",        "Isopropyl",        "CC(C)",              "Neutral",   117.15, NA_real_, "Beta-branching restricts psi; favors extended beta-sheet; large gamma-carbons",        "Beta-sheet core; hydrophobic packing; sickle-cell HbS Glu6Val",
  "L",  "Leu", "Leucine",       "Isobutyl",         "CCC(C)C",            "Neutral",   131.17, NA_real_, "Flexible chi2; common rotamers at -65/-65 and -65/175; high hydrophobic SASA",         "Hydrophobic core; leucine zippers; most abundant non-polar in proteomes",
  "I",  "Ile", "Isoleucine",    "sec-Butyl",        "CCC(C)",             "Neutral",   131.17, NA_real_, "Beta-branching + gamma-branch; most restricted chi1/chi2; large buried SASA",          "Hydrophobic core; beta-barrel interiors; transmembrane helices",
  "P",  "Pro", "Proline",       "Pyrrolidine ring", "C1CCNC1",             "Neutral",   115.13, NA_real_, "Fixed phi ~-60; no backbone NH donor; cis/trans isomerism at Xaa-Pro bond",            "Helix breaker; beta-turns; collagen Gly-Pro-X repeats",
  "F",  "Phe", "Phenylalanine", "Benzyl",           "Cc1ccccc1",          "Neutral",   165.19, NA_real_, "Rigid aromatic ring; pi-pi stacking and cation-pi in MD energy decomposition",         "Hydrophobic core; aromatic clusters; ligand binding pockets",
  "W",  "Trp", "Tryptophan",    "Indolylmethyl",    "Cc1c[nH]c2ccccc12", "Neutral",   204.23, NA_real_, "Indole NH can H-bond; amphipathic at membrane interface; strong 280nm absorbance",     "Membrane anchoring; fluorescence probe; ligand binding; rarest standard AA",
  "M",  "Met", "Methionine",    "Thioether",        "CCSC",               "Neutral",   149.20, NA_real_, "Flexible sulfur geometry; oxidizable to sulfoxide in long MD runs; check reactive FF", "Translation initiation; hydrophobic core; redox sensing"
)

aa_nonpolar |>
  dplyr::select(aa:mw_da) |>
  kbl()

aa_nonpolar |>
  dplyr::select(aa,aa3,md_note,main_function) |>
  kbl()

Claude generated most of the above information. We’ll add onto the md_note section as we encounter certain things during our MD sims.

What’s Zwitterion?

A zwitterion is a molecule that has both positive and negative charges but is overall electrically neutral. In amino acids, the amino group (–NH₂) can accept a proton to become positively charged (–NH₃⁺), while the carboxyl group (–COOH) can lose a proton to become negatively charged (–COO⁻). At physiological pH (~7.4), most amino acids exist as zwitterions, with the amino group protonated and the carboxyl group deprotonated. This dual charge allows amino acids to interact with both polar and non-polar environments, contributing to their solubility in water and their ability to form various interactions in proteins.

What Does Non-polar Actually Mean?

It is worth clarifying what “non-polar” actually refers exclusively to the side chain (R group) — specifically that it consists largely of carbon and hydrogen bonds with no net dipole and no ionizable groups, making it hydrophobic and largely indifferent to water. It says nothing about the backbone, which is the same for all amino acids and always carries polar bonds (C=O, N–H). In fact, as mentioned above, all amino acids including non-polar ones exist as zwitterions at physiological pH — a property that comes entirely from the backbone, not the side chain.

Note to self: All amino acids’ backbones are zwitterions; the R-side chain determines polarity and hydrophobicity. Also, net charge neutral == overall charges equals zero, does not mean the molecule is non-polar.

Rodriguez Rotation Formula

Rodrigues’ rotation formula is a method for rotating a 3D vector in space around a specified axis by a given angle. The formula is expressed as:

\(v_{rotation} = v.\cos(\theta) + \sin(\theta)(k \times v) + (1 - \cos(\theta))(k(k \cdot v))\)

where v is the original vector, k is the unit vector along the axis of rotation, and θ is the angle of rotation in radians.

This formula apparently is very popular in computer graphics and robotics, but I can see how it can be useful in molecular dynamics as well when we want to rotate a molecule or a part of it around an axis. Especially when we want to estimate the least energy conformation of a molecule. The direction application of this formula in amino acid sequence would be in rearranging the atoms based on phi and psi which are the angles of rotations around the N-Cα and Cα-C bonds of the amino acid backbone, respectively. By applying Rodrigues’ rotation formula, we can calculate the new positions of the atoms in the amino acid after rotating them by the specified angles, allowing us to explore different conformations of the molecule. How I remember which angle is which is Nancy Phi (sounds like some detective show and also N->C) and C C Psi (All with S sound, also Carbon to carbon). We’ll leave the hand calculation until next time, but let’s learn how to rotate a coordinate based on an axis with Rodriguez!

Below I’ll write the code first, then explain. Please feel free to use your mouse to hover over the plotly object and check out the coordinates.

library(plotly)
library(pracma)

#### Let's start simple
x1 <- c(0,0,0)
x2 <- c(1,1,1)
x3 <- c(1,2,1)

rodrigues <- function(v, k, theta) {
  k <- k / sqrt(sum(k^2))
  cos(theta)*v + sin(theta)*pracma::cross(k, v) + (1 - cos(theta))*sum(k*v)*k
}

k <- x2 - x1
v <- x3 - x1

result <- rodrigues(v,k,pi/2) #notice this, pi/2 == 90 degrees

pts <- data.frame(
  x = c(x1[1],x2[1],x3[1]),
  y = c(x1[2],x2[2],x3[2]),
  z = c(x1[3],x2[3],x3[3]),
  label = c("x1","x2","x3")
)

pts_rs <- data.frame(
  x = c(x1[1],x2[1],result[1]),
  y = c(x1[2],x2[2],result[2]),
  z = c(x1[3],x2[3],result[3]),
  label = c("x1","x2","x3_new")
)

plot_ly() |>
  add_trace(data=pts, x=~x, y=~y, z=~z,
            type="scatter3d", mode="lines+markers+text",
            text=~label,
            marker=list(size=8, color="blue", opacity=0.5),
            line=list(width=4, color="blue", dash="solid")) |>
  add_trace(data=pts_rs, x=~x, y=~y, z=~z,
            type="scatter3d", mode="lines+markers+text",
            text=~label,
            marker=list(size=8, color="red", opacity=0.5),
            line=list(width=4, color="red", dash="dash"))

So with the above, we want to start off with 3 points, x1, x2, x3. they all represent their xyz coordinates.

Funny thing is, xyz coordinate here is different from what I learnt xyz as. I’ve always thought x is horizontal, y is vertical and z is depth. But in this case, x is depth, y is horizontal and z is vertical. I guess it depends on how you look at it.

Then we want to rotate x3 around the axis defined by x1 and x2 by 90 degrees (pi/2 radians). The rodrigues function takes in the vector v (which is the vector from x1 to x3), the axis k (which is the vector from x1 to x2), and the angle theta (which is pi/2). It returns the new coordinates of x3 after rotation.

Finally, we plot the original points and the rotated point using plotly. The original points are in blue, and the rotated point is in red. You can hover over the points to see their coordinates.

Now if we were to maneuver the 3d plot and align both x1 and x2 into a dot, we can clearly see that it moved 90 degrees anti-clockwise!

Now there is a pretty cool rule to know where the rotation should occur anti-clockwise vs clockwise is by using your hand !!! Remember this from high school?

It would be really cool to derive the above formula. There are a lot of videos that have done this. I’m still trying to conceptualize it, let’s leave that for another blog! It sounds interesting and may be a good exercise, especially when we’re venturing into 3d spaces.

Lennard-Jones Potential Energy

The Lennard-Jones potential energy formula is a mathematical model used to describe the interaction between a pair of neutral atoms or molecules. It is given by the equation:

\(V(r) = 4\epsilon \left[ \left( \frac{\sigma}{r} \right)^{12} - \left( \frac{\sigma}{r} \right)^6 \right]\)

Where:

• \(V(r)\) is the potential energy as a function of the distance \(r\) between the two particles.
• \(\epsilon\) is the depth of the potential well, representing the strength of the attractive interaction.
• \(\sigma\) is the finite distance at which the inter-particle potential is zero, representing the effective diameter of the particles.
• The term \(\left( \frac{\sigma}{r} \right)^{12}\) represents the repulsive part of the potential, which dominates at short distances due to the Pauli exclusion principle.
• The term \(\left( \frac{\sigma}{r} \right)^6\) represents the attractive part of the potential, which dominates at longer distances due to van der Waals forces.

The Lennard-Jones potential is widely used in molecular dynamics simulations to model the interactions between non-bonded atoms or molecules, particularly in the context of van der Waals forces. It helps to predict the behavior of particles in a system, such as their equilibrium positions and the energy landscape of molecular interactions.

Wow there are a bunch of terms and word above! I’m getting dizzy just to keep track of what is what. Let’s push through this. To use the above formula, we’d have to have some understanding of the parameters epsilon and sigma. These parameters are typically derived from experimental data or quantum mechanical calculations and are specific to the types of atoms or molecules involved in the interaction. Where to get these parameters? Here you go – openbabel:: gaff.dat

When you opened gaff.dat there are bunch of numbers! Let’s find the numbers that are meaningful for us. All the below are separated by new lines as you scroll down.

Bond Stretch

The column names should be: type mass(g/mol) polarizability(ų) source

Bond Angle

The column names should be: types K r0 source count rmsd

Proper Dihedral

The column names should be: types div barrier phase periodicity

Non-bonded

The column names should be: type R*(Å) ε(kcal/mol)

We are purely interested in the non-bonded section where R* is our sigma = R* × 2 / 2^(1/6) and ε is our epsilon. With the above parameters, we can then calculate the Lennard-Jones potential energy between any two atoms in a molecule. Let’s do a simple calculation for ethanol.

Calculating LJ

library(tidyverse)
library(igraph)
# Coordinates (x, y, z) in Angstroms, according to pubchem ethanol molecule
coords <- rbind(
  O   = c( -1.1712,   0.2997,   0.0000),
  C2  = c( -0.0463,  -0.5665,   0.0000),
  C1  = c(  1.2175,   0.2668,   0.0000),
  H4  = c( -0.0958,  -1.2120,   0.8819),
  H5  = c( -0.0952,  -1.1938,  -0.8946),
  H1  = c(  2.1050,  -0.3720,  -0.0177),
  H2  = c(  1.2426,   0.9307,  -0.8704),
  H3  = c(  1.2616,   0.9052,   0.8886),
  H6  = c( -1.1291,   0.8364,   0.8099)
)

# AMBER GAFF parameters (sigma Å, epsilon kcal/mol)
Rstar_to_sigma <- function(Rstar) 2 * Rstar / 2^(1/6)

sigma <- c(
  C1=Rstar_to_sigma(1.9080), C2=Rstar_to_sigma(1.9080),
  O=Rstar_to_sigma(1.7210),
  H1=Rstar_to_sigma(1.4870), H2=Rstar_to_sigma(1.4870),
  H3=Rstar_to_sigma(1.4870), H4=Rstar_to_sigma(1.4870),
  H5=Rstar_to_sigma(1.4870), H6=0.0000
)

epsilon <- c(
  C1=0.1094, C2=0.1094, O=0.2104,
  H1=0.0157, H2=0.0157, H3=0.0157,
  H4=0.0157, H5=0.0157, H6=0.0000
)

# Bonds 
bonds <- tribble(
  ~from, ~to,
  "C1", "C2",
  "C2", "O",
  "O", "H6",
  "C1", "H1",
  "C1", "H2",
  "C1", "H3",
  "C2", "H4",
  "C2", "H5"
)

# count bonds between two atoms
g <- graph_from_data_frame(bonds, directed = FALSE)
g_dist <- distances(g)

# LJ function
lj <- function(r, eps, sig) 4 * eps * ((sig/r)^12 - (sig/r)^6)

# Loop all pairs
atoms <- rownames(coords)
pairs <- combn(atoms, 2, simplify=FALSE) # combn so we don't repeat
total_V <- vector(mode = "numeric", length = length(pairs)) 
 
for (i in 1:length(pairs)) {
  
  # each pair
  p <- pairs[[i]]
  from <- p[1]
  to <- p[2]
  num_bond <- g_dist[from, to]
  
  if (num_bond <= 2) next                   
  
  # params needed for LJ
  r   <- sqrt(sum((coords[from,] - coords[to,])^2))
  sig <- (sigma[from] + sigma[to]) / 2
  eps <- sqrt(epsilon[from] * epsilon[to])
  
  # scale if num bond is 3 (4 atoms)
  scale <- if (num_bond == 3) 0.5 else 1.0
  
  # LJ calc
  V <- scale * lj(r, eps, sig)

  cat(from, "-", to, " num of bonds=", num_bond, " r=", r, " V=", V, "\n")
  total_V[i] <- V
}

## O - H1  num of bonds= 3  r= 3.344395  V= -0.02733269 
## O - H2  num of bonds= 3  r= 2.642383  V= 0.110613 
## O - H3  num of bonds= 3  r= 2.659841  V= 0.09535471 
## C1 - H6  num of bonds= 3  r= 2.546942  V= 0 
## H4 - H1  num of bonds= 3  r= 2.521587  V= 0.01461147 
## H4 - H2  num of bonds= 3  r= 3.074579  V= -0.007593085 
## H4 - H3  num of bonds= 3  r= 2.514978  V= 0.01576022 
## H4 - H6  num of bonds= 3  r= 2.295394  V= 0 
## H5 - H1  num of bonds= 3  r= 2.507028  V= 0.01720963 
## H5 - H2  num of bonds= 3  r= 2.510736  V= 0.01652426 
## H5 - H3  num of bonds= 3  r= 3.070262  V= -0.007612401 
## H5 - H6  num of bonds= 3  r= 2.845344  V= 0 
## H1 - H6  num of bonds= 4  r= 3.550289  V= 0 
## H2 - H6  num of bonds= 4  r= 2.908137  V= 0 
## H3 - H6  num of bonds= 4  r= 2.392984  V= 0

cat("\nTotal V_LJ:", sum(total_V), "kcal/mol\n")

## 
## Total V_LJ: 0.2275351 kcal/mol

Alright, the above we basically were trying to calculate all pairwise atoms that is more than 3 bonds (if it’s exactly 3 bonds, we scale it by half). Alright, now that we know how to do that on simple molecular, next time we can use this to minimize on as we’re seeking optimal phi and psi!

Opportunities For Improvement

• Derive Rodriguez Rotation formula
• put both Rodriguez rotation formula and LJ calculation into action to find the optimal phi and psi of amino acid sequence of a protein!
• need to include secondary/tertiary structure interactions too

Lessons learnt

• refreshed on rotation and vectors
• learnt rodriguez rotation formula
• learnt LJ formula
• learnt what gaff.dat actually contains.
• learnt about phi and psi and what they actually mean.
• learnt about net charge == total charge; polar vs non-polar is R-chain dependent.

If you like this article:

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

One of the challenges of working on R-universe is that there is never really a release day.

Unlike software projects that accumulate changes for a big launch, R-universe evolves continuously. New features, infrastructure improvements, UI tweaks, and build system enhancements are silently deployed all the time without most people noticing.

Every now and then, however, a few features emerge that are worth highlighting. In this technote we look at five recent additions that make R-universe a little nicer, faster, or more convenient to use.

1. Social media cards that actually look good

Sharing package links on social media used to be a somewhat underwhelming experience, but not anymore! We provide beautiful preview images every package, article, and universe, for example:

• https://ropensci.r-universe.dev/card.png
• https://ropensci.r-universe.dev/card.svg
• https://ropensci.r-universe.dev/targets/card.png
• https://ropensci.r-universe.dev/targets/card.svg

Each card includes package or universe stats and is automatically exposed through the appropriate HTML headers (og:image, og:title, etc). Whenever somebody shares a package link, the preview should look a bit more polished without requiring any work from package maintainers.

R-universe now generates social media preview cards for each package, like this one: ropensci.r-universe.dev/targets. You can also get the card manually from the /{package}/card.png API (or svg).

[image or embed]

— Jeroen Ooms (@jeroenooms.bsky.social) 12:01 · May 2, 2026

When a link to a vignette article is shared, R-universe automatically extracts the title and section headings from the document to generate a more informative description. For example this one.

All this won’t guarantee your package goes viral, but at least it looks cool

2. PACKAGES.rds support (or: implementing R internals in JavaScript)

This feature is mostly invisible, but improves performance of installing packages in R, and therefore also the workflow of building packages in R-universe:

Every CRAN-like repository needs an index file which lists all the content from that repo. This file may be provided in a text-based PACKAGES format and/or a binary PACKAGES.rds format (rds is R’s internal binary serialization format, see ?saveRDS).

Historically R-universe supported only the former text-based format, because all repository metadata is generated on-request in JavaScript on the server side, and emitting DCF text streams from our database is fast and easy. However, on the R side, loading RDS is a bit faster than parsing a text, which becomes noticeable for large repositories like https://bioc.r-universe.dev/.

So therefore we now also serve the PACKAGES.rds files. The implementation exists in this NPM package which reverse engineers a subset of the R RDS serializer, so that we can easily run it in our express stack. On MacOS and Windows it defaults to the new zstd compression, which makes it even faster than CRAN.

3. Fancy sort/filter bars in the WebUI

The styling of universe-level pages that list packages, articles, and datasets have been improved, gaining some nice interactive filter and sort capabilities. For example the /packages page now allows you to do a (fuzzy) search looking for keywords that appear in package descriptions/tags/authors/etc, and sort the packages based on their of stars / downloads / dependents / etc.

A similar filter bar is available on the /articles and /datasets pages to help you search those as well.

4. For the impatient: trigger a sync manually

R-universe automatically checks for updates in upstream git repositories and package registries approximately once per hour. Occasionally, however, you have just pushed a commit, fixed a build issue, updated a vignette, or corrected a typo, and waiting an hour suddenly feels like a very long time.

To accommodate the impatient among us, a new sync button has been added to the universe sidebar. Clicking the button immediately triggers a sync to check for any updates.

5. Making check results easier to find and share

For some organizations, package checks are among the most important parts of R-universe. We’ve made several improvements to make check results easier to access and easier to share with collaborators.

First, package pages now support direct links to the check table using the #checktable anchor, for example: https://jeroen.r-universe.dev/curl#checktable. Second, build logs and build artifacts linked in this table can now be downloaded without requiring GitHub authentication. The underlying files still live on GitHub Actions, but R-universe now proxies the download links. This means users can access logs and build artifacts directly from the package page, even if they do not have a GitHub account. So “I don’t have GitHub” is no longer available as an excuse for ignoring check failures.

Finally, build runs now include a deployment summary generated via the GitHub Actions Job Summaries feature: navigate to any build run and scroll down, there you find a summary table showing exactly the data deployed to R-universe during that run, including package checks and deployment results. This makes it easier to inspect builds directly from GitHub without having to cross-reference multiple pages.