Forty-eight answers to one question: Reproducible research with generative AI

Why the Same Request Gets Different Replies

The obvious suspect is the randomness of sampling, which the temperature setting controls. A temperature of zero should select the most probable token at every step, yet it does not make the models that providers host deterministic. When Ouyang et al. (2025) asked ChatGPT for code several times per problem, the programs often gave different outputs, and a temperature of zero reduced that variation without removing it. Atıl et al. (2025) ran five models 10 times on each of eight tasks with settings meant to be deterministic, and accuracy still varied by up to 15 percentage points across runs.

He and Thinking Machines Lab (2025) traced much of this behaviour to the providers’ servers, where the number of requests processed together varies with demand. That number changes the order, and so the rounding, of some sums, which can be enough to tip the choice between two almost equally probable tokens. With an open model, one that anyone can download and run, 1,000 replies to one request at temperature zero contained 80 distinct texts, and computations made independent of that number brought them down to one. Models that reason before answering are especially exposed, since an early rounding difference can send the reasoning down another path (Yuan et al., 2025).

Sampling is only one source of variation. Providers can change proprietary models at any time (Spirling, 2023) and retire them, after which requests to them fail (Anthropic, n.d.-b). Between its March and June 2023 versions, GPT-4 fell from 84% to 51% accuracy in classifying numbers as prime or composite (L. Chen et al., 2024). Meaning-preserving changes to the format of a prompt, the text given to a model, altered the accuracy of one open model by up to 76 percentage points (Sclar et al., 2024). In a chat window, none of these sources is under the user’s control, and few are recorded.

What the Researcher Can Still Control

Reproducibility with generative AI therefore rests on three choices left to the researcher: the tools, the requests made of them and the record of both. For the tools, McLoughlin et al. (2026) recommend open models, which can be run locally and so run again later, and, like Palmer et al. (2024), they ask that any use of a proprietary model be justified. Even an open model may not match the published one when a company serves it. In summer 2024, 11 of 31 services offering Meta’s Llama models returned outputs whose distribution differed from that of the published model (Gao et al., 2025).

The choice of tools extends to their settings, such as the effort level, which sets how many tokens a model spends on reasoning, text and tool use, and whose default differs between models (Anthropic, n.d.-a). The second round of this experiment ran at medium effort. The access route matters too, since a model may behave differently in a chat window and through an application programming interface (API), which lets other programs send it requests (Feuerriegel et al., 2026). So does the agentic tool around the model. Nadgir et al. (2026) gave OpenAI’s GPT-5.4 39 tasks from CORE-Bench, a benchmark that asks agents to reproduce published results from the studies’ code and data. At the same effort level, it solved 94.9% of them within Codex, the company’s own agentic tool, and 51.3% within CORE-Agent, an agent built for the benchmark.

A request to an agent can ask for a script in place of an answer, one that runs without manual steps, prints its result in a fixed format and sets any random seeds. It can also ask for tests of the expected behaviour and a file that pins package versions, as Bridgeford et al. (2026) recommend. What comes back still needs checking. Vangala et al. (2026) asked three coding agents for programs that would run in a clean environment, with every dependency declared. Of the 300 programs, 205 ran without intervention, and on average, they used 15 times as many packages at run time as they declared.

Asking for a script brings AI-assisted work under the standards of reproducible computational research, which ask that every result be traceable to the steps that produced it and that the versions of the software be archived (Sandve et al., 2013; see also Peng, 2011; Wilson et al., 2017). A script kept under version control, in a system such as Git that records each version of a project’s files, documents the steps but not the software environment. When Trisovic et al. (2022) re-executed R files deposited with published research data, 75% of those that did not time out stopped with an error, and 60% still did after automatic repairs. The renv package (Ushey & Wickham, 2026) reduces that risk with a lockfile, which records the versions of R and of a project’s packages, and the output of sessionInfo() shows the environment that actually ran the code.

Reporting guidelines set out the rest of the record. The GUIDE-LLM checklist asks authors in the behavioural sciences to report the provider, the exact model identifier, the date and route of access, the settings and the prompts (Feuerriegel et al., 2026). Guidelines for software engineering add the version of the agentic tool and a trace of its tool calls, the actions it takes, such as reading a file or running a command (Baltes et al., 2026).

A Data Set With a Known Answer

I simulated the data so that the population effect would be known. From the same specification and seed, pilotr produces identical files in R and in Python, at least for a specification such as this one, which rounds reaction times to whole milliseconds (Bernabeu, 2026b). In this data set, 30 participants each respond to 24 words, once after a related prime and once after an unrelated one. The specification below sets a shifted lognormal distribution for the reaction times, a word-frequency covariate, and random intercepts and priming slopes for both participants and words.

cat(readLines(file.path(experiment, 'spec.json')), sep = '\n')

#> {
#>   "spec_version": "0.3",
#>   "name": "semantic_priming_lexical_decision",
#>   "seed": 2026,
#>   "units": { "subject": { "n": 30 }, "item": { "n": 24 } },
#>   "factors": [
#>     {
#>       "name": "prime",
#>       "levels": ["related", "unrelated"],
#>       "contrasts": { "prime_c": [-0.5, 0.5] },
#>       "vary_within": ["subject", "item"]
#>     }
#>   ],
#>   "predictors": [
#>     { "name": "frequency", "varies_by": "item", "mean": 4, "sd": 0.7 }
#>   ],
#>   "fixed": {
#>     "intercept": 6.24,
#>     "coefficients": { "prime_c": 0.05, "frequency": -0.06 }
#>   },
#>   "random": {
#>     "subject": {
#>       "intercept_sd": 0.15,
#>       "slopes": { "prime_c": 0.04 },
#>       "correlations": { "intercept,prime_c": 0.2 }
#>     },
#>     "item": {
#>       "intercept_sd": 0.08,
#>       "slopes": { "prime_c": 0.03 },
#>       "correlations": { "intercept,prime_c": -0.1 }
#>     }
#>   },
#>   "response": {
#>     "family": "shifted_lognormal", "name": "rt",
#>     "sigma": 0.3, "shift": 250, "round": 0
#>   }
#> }

The code below regenerates the data in R and compares its SHA-256 checksum, a fingerprint of the file’s contents, with those of the file given to the agents and of the file generated in Python.

library(pilotr)

spec <- load_spec(file.path(experiment, 'spec.json'))
simulated <- simulate_design(spec)
regenerated <- data.frame(
  participant = sprintf('P%02d', simulated$subject),
  item = sprintf('W%02d', simulated$item),
  prime = simulated$prime,
  frequency = round(simulated$frequency, 2),
  rt = simulated$rt
)
# write.csv() ends each line with '\n', and a binary connection stops Windows
# from turning it into '\r\n', so the file has the same bytes on every system.
regenerated_file <- tempfile(fileext = '.csv')
connection <- file(regenerated_file, 'wb')
write.csv(regenerated, connection, row.names = FALSE, quote = FALSE)
close(connection)

checksums <- c(
  `regenerated in R` = digest::digest(file = regenerated_file, algo = 'sha256'),
  `given to the agents` = digest::digest(file = file.path(experiment, 'data.csv'), algo = 'sha256'),
  `recorded from Python` = readLines(file.path(experiment, 'python_sha256.txt'))[2]
)
stopifnot(n_distinct(checksums) == 1)
checksums

#>                                                   regenerated in R 
#> "2e2ad9e3cdd2d4fe3a6a5b1e4670f4b36518276e3455110600979daaec15df61" 
#>                                                given to the agents 
#> "2e2ad9e3cdd2d4fe3a6a5b1e4670f4b36518276e3455110600979daaec15df61" 
#>                                               recorded from Python 
#> "2e2ad9e3cdd2d4fe3a6a5b1e4670f4b36518276e3455110600979daaec15df61"

All three checksums match, so the three files are identical. The next calculation draws a million participants and a million words to estimate how much the priming effect varies across each in the population.

fixed <- spec$fixed; random <- spec$random; frequency <- spec$predictors[[1]]
# effect_sd() below moves frequency with the word, so it must vary by word
stopifnot(frequency$name == 'frequency', frequency$varies_by == 'item')
mu <- fixed$intercept + fixed$coefficients$frequency * frequency$mean

# The variance, at a given prime code, of the terms that belong to one unit
# (participant or item): its random intercept and slope and their covariance.
unit_variance <- function(unit, prime_code) {
  unit$intercept_sd^2 + prime_code^2 * unit$slopes$prime_c^2 + 2 * prime_code *
    unit$intercept_sd * unit$slopes$prime_c * unit$correlations[['intercept,prime_c']]
}
residual_variance <- (fixed$coefficients$frequency * frequency$sd)^2 + spec$response$sigma^2

# The mean of a lognormal variable is exp(mu + variance / 2), and the shift
# cancels out of the difference between the two types of prime.
population_mean <- function(prime_code) {
  spec$response$shift + exp(mu + fixed$coefficients$prime_c * prime_code +
    (unit_variance(random$subject, prime_code) + unit_variance(random$item, prime_code) +
       residual_variance) / 2)
}
population_effect <- population_mean(0.5) - population_mean(-0.5)

# The priming effect of one participant (or word), averaged over the other
# unit and the residual, for a million simulated participants (or words). A
# word's frequency belongs to the word, so for words it moves with the word's
# intercept and is not averaged out.
effect_sd <- function(unit, other, word = FALSE) {
  z <- matrix(rnorm(2e6), ncol = 2)
  rho <- unit$correlations[['intercept,prime_c']]
  intercept <- unit$intercept_sd * z[, 1]
  slope <- unit$slopes$prime_c * (rho * z[, 1] + sqrt(1 - rho^2) * z[, 2])
  averaged <- residual_variance
  if (word) {
    intercept <- intercept + fixed$coefficients$frequency * frequency$sd * rnorm(nrow(z))
    averaged <- spec$response$sigma^2
  }
  unit_mean <- function(prime_code) exp(mu + intercept + prime_code * (fixed$coefficients$prime_c + slope) +
    (unit_variance(other, prime_code) + averaged) / 2)
  sd(unit_mean(0.5) - unit_mean(-0.5))
}
# The seed, that of spec.json, makes the printed standard deviations
# reproducible. Other seeds move them by a few hundredths of a millisecond,
# which leaves the one-decimal values in the text unchanged (checked below).
set.seed(2026)
effect_sds <- c(participants = effect_sd(random$subject, random$item),
                words = effect_sd(random$item, random$subject, word = TRUE))
sample_effect <- with(regenerated, mean(rt[prime == 'unrelated']) - mean(rt[prime == 'related']))
stopifnot(effect_sds[['participants']] > effect_sds[['words']])  # as the text below says

c(population_effect = population_effect, sample_effect = sample_effect,
  sd_across_participants = effect_sds[['participants']], sd_across_words = effect_sds[['words']])

#>      population_effect          sample_effect sd_across_participants        sd_across_words 
#>               21.84645               18.51944               18.28419               12.89411

In the population, responses after unrelated primes are 21.8 ms slower than after related primes, and the effect varies more across participants (standard deviation 18.3 ms) than across words (12.9 ms). In the sample of 1,440 trials, the difference between the mean reaction times is 18.52 ms.

Forty-Eight Runs of One Question

results <- read.csv(file.path(experiment, 'results.csv'), encoding = 'UTF-8',
                    colClasses = c(run_id = 'character'))
executions <- read.csv(file.path(experiment, 'executions.csv'), colClasses = c(run_id = 'character'))

# The logs' times are in UTC, which the format below assumes without reading
# the zone, so the zone is checked
utc_time <- function(x) {
  stopifnot(grepl('(Z|[+]00:00)$', x))
  as.POSIXct(x, format = '%Y-%m-%dT%H:%M:%OS', tz = 'UTC')
}
logs <- lapply(results$run_id, function(id)
  jsonlite::fromJSON(file.path(experiment, 'runs', id, 'tool_calls.json')))
results$started <- utc_time(vapply(logs, `[[`, character(1), 'started'))
results$starting_directory <- vapply(logs, `[[`, character(1), 'starting_directory')

first_date <- as.Date(min(results$started))
# Runs launched together started within a second of one another, and each
# batch at least 40 seconds after the one before, so a gap of more than 10
# seconds marks a new batch. Every batch of a round had the same size.
batch_size_of <- function(started) {
  gaps <- diff(sort(as.numeric(started)))
  size <- unique(as.vector(table(cumsum(c(TRUE, gaps > 10)))))
  stopifnot(all(gaps < 1 | gaps > 40), length(size) == 1)
  size
}
batch_size <- batch_size_of(results$started)
stopifnot(n_distinct(as.Date(results$started)) == 1)
main_runs <- results$data == 'balanced'
in_folder <- grepl('/experiment$', results$starting_directory)
in_experiment_folder <- sum(in_folder)
script_runs <- results$route == 'script'
unbalanced <- read.csv(file.path(experiment, 'data_unbalanced.csv'))

# The second round: the same 48 runs at medium effort, and two of them at
# xhigh, with what each run's transcript recorded
second <- read.csv(file.path(experiment, 'rerun', 'results.csv'), encoding = 'UTF-8',
                   colClasses = c(run_id = 'character'))
second$started <- utc_time(second$started)
second$finished <- utc_time(second$finished)
medium <- second[second$round == 'medium', ]
xhigh <- second[second$round == 'xhigh', ]

Each run received a description of the data and this question: ‘How large is the priming effect, that is, how much slower are responses after unrelated primes than after related primes? Give the estimate in milliseconds with a 95% confidence interval.’ It also had to end with a result line in a fixed format, so that code could compare the answers.

The request was either open, leaving the analysis to the agent, or specified, naming the model, the contrast coding and the function for the interval. Both came with the balanced data. A follow-up gave the specified request with an unbalanced version of the data, from which 86 of the 1,440 trials had been removed at random. Its prompt described them as incorrect responses that had been excluded. These make up the experiment’s three conditions.

Each condition had two routes. The prose agents could read the data but were told to run no code and to approximate where necessary, as a chat assistant without a code interpreter (a tool for running code) would have to. The script agents wrote an R script, could execute it to check their work and replied with a brief summary. Their instructions named lme4 (Bates et al., 2015), lmerTest (Kuznetsova et al., 2017), emmeans and dplyr among the installed packages.

The experiment folder holds the prompts, the data, the scripts, the final replies and the logs of the agents’ tool calls, with a README for repeating the work, although not the agents’ full transcripts. Below is the specified request as one script run received it. The open request lacks only the paragraph beginning ‘Analyse the data as follows’.

example_run <- design$run_id[design$route == 'script' & design$request == 'specified' &
                               design$data == 'balanced'][1]
cat(strwrap(readLines(file.path(experiment, 'runs', example_run, 'prompt.txt')), width = 90),
    sep = '\n')

#> The file data.csv holds the correct responses from a lexical decision experiment on
#> semantic priming. Each of 30 participants judged 24 target words, and each target
#> appeared once after a semantically related prime and once after an unrelated prime. The
#> columns are participant, item, prime (related or unrelated), frequency (the Zipf
#> frequency of the target word) and rt (the response time in milliseconds).
#> 
#> How large is the priming effect, that is, how much slower are responses after unrelated
#> primes than after related primes? Give the estimate in milliseconds with a 95% confidence
#> interval.
#> 
#> Analyse the data as follows. Use every row, with no exclusions and no transformation of
#> rt. Code prime as -0.5 for related and +0.5 for unrelated, in a column named prime_c. Fit
#> a linear mixed-effects model with lme4::lmer(), using REML and the default optimiser,
#> with the formula rt ~ prime_c + (1 + prime_c | participant) + (1 + prime_c | item). The
#> estimate is the fixed effect of prime_c, and the interval is its 95% Wald confidence
#> interval from confint(fit, parm = "prime_c", method = "Wald").
#> 
#> Your working directory for this task is /home/user/agent_runs/437afbd0, which contains
#> data.csv. Write an R script named analysis.R in that directory that reads data.csv and
#> carries out the analysis. You can run R there (Rscript analysis.R) to check the script. R
#> 4.3 is installed with packages that include lme4, lmerTest, emmeans and dplyr, and there
#> is no internet access. Work only inside that directory and do not look at other
#> directories. The last line the script prints must have exactly this form, with numbers in
#> milliseconds:
#> RESULT estimate=<number> lower=<number> upper=<number>
#> When the script is finished, reply with a brief summary of what it does.

Both rounds used the language model Claude Opus 5.5 (identifier claude-opus-5-5), with eight runs in each combination of condition and route. The first round ran on 29 September 2026 with Claude Code 2.1.284. Its 32 runs on the balanced data were launched eight at a time, in a random order across conditions. The 16 runs of the follow-up were launched afterwards, also in a random order.

Each agent started in a folder of its own that held only its data file, and the logs show that the prose agents read nothing else. In the first round, however, 3 of the 24 script agents also wrote exploratory files outside their folders, although their commands show no sign of executing another agent’s code. Shared temporary folders, and for some agents the experiment folder itself, were within reach, so the isolation was imperfect, and the logs and the hidden checks on this page document each exception.

Each final script was then executed three times, each time in a fresh R process and a fresh copy of its folder. A record kept the versions of R, lme4, lmerTest and Matrix, the package that lme4 uses for matrix computations.

The second round ran on 11 October 2026 with Claude Code 2.1.295 on Windows 11. Its 48 main runs used medium effort and followed the first round’s launch order, eight at a time. Their prompts changed only the folder and, for the scripts, the version of R and the path to it. I added two runs on the balanced data at xhigh effort, a prose reply to the specified request and a script for the open request. That level lies between high and the maximum, and Anthropic (n.d.-a) recommends it for long-running agentic and coding work. No agent of the second round touched a file outside its folder.

Different Words, Mostly the Same Numbers

condition_levels <- c('Balanced data, open request', 'Balanced data, specified request',
                      'Unbalanced data, specified request')
route_levels <- c('Prose reply', 'Executed script')

text_of <- function(run_id, file) {
  paste(readLines(file.path(experiment, 'runs', run_id, file), warn = FALSE, encoding = 'UTF-8'),
        collapse = '\n')
}

runs <- results |>
  mutate(condition = factor(case_when(
           data == 'balanced' & request == 'open' ~ condition_levels[1],
           data == 'balanced' ~ condition_levels[2],
           TRUE ~ condition_levels[3]), levels = condition_levels),
         route_label = factor(ifelse(route == 'prose', route_levels[1], route_levels[2]),
                              levels = route_levels),
         reply = vapply(run_id, text_of, character(1), file = 'response.txt'),
         text = mapply(function(id, rte) text_of(id, if (rte == 'script') 'analysis.R' else 'response.txt'),
                       run_id, route))

# The numbers of each RESULT line as written, so that a reply's precision is
# kept (a reply's '22.0' is not the same claim as a script's '21.79')
result_numbers <- function(line) {
  regmatches(line, regexec('estimate=(-?[0-9.]+) lower=(-?[0-9.]+) upper=(-?[0-9.]+)', line))[[1]][-1]
}
# Each result line as the tables show it, 'estimate [lower, upper]'
as_written <- function(lines) {
  vapply(lines, function(line) {
    n <- result_numbers(line)
    sprintf('%s [%s, %s]', minus(n[1]), minus(n[2]), minus(n[3]))
  }, character(1), USE.NAMES = FALSE)
}
runs$result_text <- as_written(runs$result_line)

# The paired differences (unrelated minus related) in the balanced data, and
# each participant's and each word's mean difference
balanced <- read.csv(file.path(experiment, 'data.csv'))
paired <- reshape(balanced[, c('participant', 'item', 'prime', 'rt')],
                  idvar = c('participant', 'item'), timevar = 'prime', direction = 'wide')
paired$difference <- paired$rt.unrelated - paired$rt.related
# reshape() would keep only the first of any repeated trial, and a pair with a
# trial missing would leave a missing difference, so neither must occur
stopifnot(nrow(paired) * 2 == nrow(balanced), !anyNA(paired$difference))
participant_means <- tapply(paired$difference, paired$participant, mean)
word_means <- tapply(paired$difference, paired$item, mean)
# The by-participant paired t interval, which most open-request replies mention
paired_ci <- t.test(participant_means)$conf.int

summary_table <- runs |>
  group_by(condition, route_label) |>
  summarise(distinct_texts = n_distinct(text),
            distinct_results = n_distinct(result_line),
            result = first(result_text), .groups = 'drop')

All 48 texts differed, as long generated texts will, but their results agreed far more closely (Table 1). On the balanced data, the eight runs of each route reported the same result in each condition, and the prose replies agreed with the executed scripts to the precision they gave. On the unbalanced data, the eight scripts again printed a single result, whereas the prose replies gave seven different results, none of which matched the scripts. Figure 1 shows how far each of those replies lay from the executed result, in the theme and colours of my depictr package (Bernabeu, 2026a).

# The result is shown where every run in the row reported the same one. The
# checks above allow more than one only in the prose replies to the unbalanced data

# Pandoc sizes the columns of a wide table in proportion to the dashes under
# its header, which kable counts from the longest cell, markup included. Where
# that misjudges a column, the dashes are set to fixed percentages of the
# width instead. Here the six characters of each   would widen the
# column of routes, so the shares are those pandoc gave this table without them.
column_shares <- function(table, shares) {
  separator <- grep('^\\|[-:|]+\\|$', table)
  cells <- strsplit(sub('^\\|(.*)\\|$', '\\1', table[separator]), '|', fixed = TRUE)[[1]]
  stopifnot(length(separator) == 1, length(cells) == length(shares), sum(shares) == 100)
  table[separator] <- paste0('|', paste0(ifelse(startsWith(cells, ':'), ':', ''), strrep('-', shares),
                                         ifelse(endsWith(cells, ':'), ':', ''), collapse = '|'), '|')
  table
}

column_shares(knitr::kable(
  summary_table |>
    arrange(condition, desc(route_label)) |>
    # A route's two words stay together, so that both routes take one line
    transmute(Condition = condition, Route = gsub(' ', '&nbsp;', route_label, fixed = TRUE),
              `Distinct texts` = distinct_texts, `Distinct results` = distinct_results,
              `Estimate [95% interval] (ms)` = if_else(distinct_results == 1, result, 'See Table 3')),
  caption = sprintf(paste(
    'Texts and Results of the First Round in Each Condition and Route. Each row covers %s runs. The texts are the',
    'replies in the prose route and the scripts in the script route. A result is the estimate and interval',
    'in a run\'s result line, as the run wrote it, and is shown when every run in the row reported the',
    'same one.'), count_text(runs_per_cell)),
  align = c('l', 'l', 'r', 'r', 'r')), shares = c(32, 14, 13, 15, 26))

# deviations, built with the hidden checks earlier on the page, holds each of the three
# numbers of every run minus the executed scripts' number in its condition, once the
# latter is rounded to the decimals the run wrote. Every script, and every run on the
# balanced data, lies at zero, so the figure plots the prose replies to the unbalanced
# data against a line at zero for the executed result.

# Dots closer than the gap, about the width of a dot, stack upwards, so that every run
# stays visible
stack_levels <- function(x, gap) {
  levels <- integer(length(x))
  for (i in seq_along(x)) {
    taken <- levels[seq_len(i - 1)][abs(x[seq_len(i - 1)] - x[i]) < gap]
    level <- 0L
    while (level %in% taken) level <- level + 1L
    levels[i] <- level
  }
  levels
}

fig_deviations <- function(narrow = FALSE, base_size = if (narrow) 14 else 12.5) {
  dots <- deviations |>
    filter(condition == condition_levels[3], route_label == route_levels[1]) |>
    group_by(number) |> arrange(deviation, .by_group = TRUE) |>
    mutate(level = stack_levels(deviation, gap = 0.035)) |> ungroup()
  top <- max(dots$level) + 1.9
  first_panel <- factor(levels(dots$number)[1], levels = levels(dots$number))
  notes <- data.frame(number = first_panel,
                      x = c(0.012, max(dots$deviation[dots$number == first_panel]) + 0.045),
                      y = c(top - 0.35, 0.15),
                      label = c(sprintf('Executed result,\nall %s scripts', count_text(runs_per_cell)),
                                sprintf('%s prose\nreplies', stringr::str_to_sentence(count_text(nrow(unbalanced_prose))))),
                      colour = c('ink', 'muted'))
  ggplot(dots, aes(deviation, level)) +
    geom_vline(xintercept = 0, colour = fig_roles[['script']], linewidth = 1.1) +
    geom_point(colour = fig_roles[['prose']], size = 2.6, shape = 16) +
    geom_text(aes(x, y, label = label), data = notes, inherit.aes = FALSE, hjust = 0, vjust = c(1, 0),
              lineheight = 0.95, size = label_size(base_size), colour = fig_roles[notes$colour]) +
    facet_wrap(~number, nrow = if (narrow) 3 else 1) +
    scale_x_continuous(breaks = seq(-0.4, 0.6, 0.2), labels = minus_labels(1),
                       limits = c(-0.32, 0.56), expand = expansion(0)) +
    scale_y_continuous(limits = c(-0.7, top), expand = expansion(0)) +
    labs(x = 'Prose reply minus executed result (ms)', y = NULL) +
    theme_figure(base_size) +
    theme(axis.text.y = element_blank(), panel.spacing = unit(if (narrow) 0.9 else 1.4, 'lines'))
}
# The axis must hold every reply
stopifnot(range(unbalanced_prose_deviations) > -0.32, range(unbalanced_prose_deviations) < 0.56)

draw_narrow(fig_deviations(narrow = TRUE), width = 4.4, height = 5.2)
fig_deviations()

What the Agents Chose When Left to Choose

Given the open request, the eight scripts converged on one analysis, a linear mixed-effects model of the untransformed reaction times with crossed random effects for participants and words. The maximal model, with by-participant and by-word priming slopes, gave a singular fit, one in which a random-effect variance is estimated at zero or a correlation at ±1. Its by-participant slope was perfectly correlated with the intercept, and every script dropped that slope. All eight then computed the interval from a t distribution with Satterthwaite’s approximation of the degrees of freedom, which lmerTest adds to models fitted with lme4 (Kuznetsova et al., 2017). Four also included word frequency as a covariate, although in a balanced design it cannot change the estimate, since every word appears after both types of prime.

Dropping that slope follows the parsimonious strategy of Matuschek et al. (2017) but sits uneasily with the maximal structure recommended by Barr et al. (2013). The sample suggested less variation across participants than across words, the reverse of the population pattern, so the agents faced the decision as an analyst would with real data, without the population values to go by.

The prose replies to the open request reached the same interval as the scripts by decomposing the 720 paired differences into participant, word and residual components. In a balanced design, this decomposition gives the same standard error as the scripts’ model without by-participant slopes. Five of these replies equated their decomposition with the maximal model, although that model keeps a small by-participant slope and gives a slightly wider Satterthwaite interval, [−1.78, 38.82] ms. With the specified request, both routes returned the Wald interval of the maximal model, which is narrower than either Satterthwaite interval because it takes its critical value from the normal distribution.

Every prose reply to the open request, and six of the script summaries, also considered a paired t test over participants, whose interval, [5.58, 31.46] ms, excludes zero. They set it aside because it ignores the variation across words and so overstates the evidence for generalising beyond the words sampled. The effect is real here, so that test would have reached the right conclusion, but through an analysis too liberal for the purpose. The mixed-model intervals include zero, and they also contain the population effect of 21.8 ms.

Mixed Models by Hand

An agent that ran no code could still fit a mixed model, and the prose replies to the balanced data showed how. Because every participant saw every word after both types of prime, the restricted maximum likelihood (REML) fit of the mixed model reduces to a few mean squares and cross-products of the paired differences and of their sums. Cross-products are the counterparts of mean squares for covariances. Table 2 recomputes in R the intermediate quantities that recur across the replies. Where a reply does not state the residual mean square as R gives it, it rounds the value or gives the residual variance, which is half of it.

library(lme4)

decomposition <- anova(lm(difference ~ participant + item, data = paired))
# The standard error of the model without by-participant slopes, from the mean
# square of the words
reduced_se <- sqrt(decomposition['item', 'Mean Sq'] / nrow(paired))
balanced$prime_c <- ifelse(balanced$prime == 'unrelated', 0.5, -0.5)
specified_fit <- lmer(rt ~ prime_c + (1 + prime_c | participant) + (1 + prime_c | item),
                      data = balanced)

comma <- function(x) format(round(x), big.mark = ',', trim = TRUE)
quantities <- tibble::tribble(
  ~key, ~quantity, ~value,
  'sum', sprintf('Sum of the %s paired differences (ms)', comma(nrow(paired))), comma(sum(paired$difference)),
  'ms_participants', 'Mean square of the differences, participants', comma(decomposition['participant', 'Mean Sq']),
  'ms_items', 'Mean square of the differences, words', comma(decomposition['item', 'Mean Sq']),
  'ms_residual', 'Mean square of the differences, residual', comma(decomposition['Residuals', 'Mean Sq']),
  'se_reduced', 'Standard error without by-participant slopes (ms)', sprintf('%.2f', reduced_se),
  'se_specified', 'Standard error of the maximal model (ms)', sprintf('%.2f', sqrt(vcov(specified_fit)[2, 2])),
  'paired_interval', 'Paired *t* interval over participant means (ms)',
    sprintf('[%.1f, %.1f]', paired_ci[1], paired_ci[2])
)

# A reply states a value when every number in it appears in the reply with no
# digit or decimal point before it and no digit after it, ignoring thousands
# separators.
states <- function(value, texts) {
  numbers <- gsub(',', '', regmatches(value, gregexpr('[0-9][0-9,]*\\.?[0-9]*', value))[[1]])
  patterns <- sprintf('(?<![0-9.])%s(?![0-9])', gsub('.', '\\.', numbers, fixed = TRUE))
  texts <- gsub('(?<=[0-9]),(?=[0-9]{3})', '', texts, perl = TRUE)
  sum(Reduce(`&`, lapply(patterns, grepl, x = texts, perl = TRUE)))
}
prose_balanced <- runs$route == 'prose' & runs$data == 'balanced'
open_replies <- runs$reply[prose_balanced & runs$request == 'open']
specified_replies <- runs$reply[prose_balanced & runs$request == 'specified']
quantities$open <- vapply(quantities$value, states, integer(1), texts = open_replies)
quantities$specified <- vapply(quantities$value, states, integer(1), texts = specified_replies)

count_of <- function(k, n) sprintf('%d of %d', k, n)
knitr::kable(
  quantities |> transmute(Quantity = quantity, `Value in R` = value,
                          `Open request` = count_of(open, length(open_replies)),
                          `Specified request` = count_of(specified, length(specified_replies))),
  caption = paste(
    'Intermediate Quantities Recomputed in R and the Number of the First Round\'s Prose Replies to the',
    'Balanced Data Stating Each. Mean squares are in squared milliseconds. The last two columns count the replies to each request',
    'that state the value as R gives it, to the precision shown, with or without a thousands separator.'),
  align = c('l', 'r', 'r', 'r'))

From these quantities, every reply to the specified request recovered the main features of lme4’s fit. Each found a by-participant slope with a standard deviation of about 9 ms, whose correlation with the intercepts lay at the boundary (+1 in seven replies and ±1 in the other).

Where Hand Calculation Runs Out

Removing trials takes that shortcut away. Some participant-word pairs lose one observation or both, and the REML fit then has to be found iteratively. All eight scripts estimated a priming effect of 21.79 ms, which differs from the raw difference in means of 21.43 ms. The specified Wald interval now excludes zero, with a lower bound of 0.58 ms, whereas with Satterthwaite degrees of freedom, the same model’s interval still includes it.

The prose replies worked from weighted per-word effects and from variances estimated from the observed mean squares. Each put the by-participant slope variance at or near zero, whereas the numerical fit kept that slope, with a standard deviation of 6.7 ms, perfectly correlated with the intercepts. The replies therefore lay closer to the model without that slope than to the executed result (Table 3).

# unbalanced_script (the executed result, which every script printed) and
# unbalanced_prose (the prose replies, in launch order) are built with the
# hidden checks earlier on the page.
knitr::kable(
  bind_rows(
    unbalanced_script |> transmute(Source = sprintf('Executed scripts (all %d)', runs_per_cell), result_text),
    unbalanced_prose |> transmute(Source = sprintf('Prose reply %d', row_number()), result_text)
  ) |> rename(`Estimate [95% interval] (ms)` = result_text),
  caption = paste('Results for the Unbalanced Data in the First Round, as Each Run Wrote Them. The prose',
                  'replies are numbered',
                  'in launch order.'),
  align = c('l', 'r'))

The largest discrepancy in any reported number was only 0.42 ms (Figure 1), and every prose reply called its numbers approximations, with a margin wider than its discrepancy. All their lower bounds, however, lay above the executed one, which was itself only 0.58 ms above zero, close enough for the exact fit to decide whether the interval excluded zero.

How Long Each Route Took

When the analysis was specified, code was by far the quicker route (Figure 2). The median script took 0.7 minutes from request to final reply, against 13.7 for a prose reply on the balanced data and 14.6 on the unbalanced data. Once written, the scripts computed their result in 2.4 seconds or less, R’s start-up included.

Given the open request, the scripts took longer, since their agents compared models and ran checks, and their times overlapped those of the prose replies. The agents shared access to the model and ran in batches, so these durations include waiting as well as computation and are best read as orders of magnitude.

# A separate row for each route makes every median unambiguous, including in the open
# condition, where the two routes' durations overlap
fig_minutes <- function(narrow = FALSE, base_size = if (narrow) 14 else 12.5) {
  points <- runs |>
    group_by(condition, route_label) |> arrange(order, .by_group = TRUE) |>
    mutate(y = as.numeric(route_label) + seq(-0.15, 0.15, length.out = n())) |> ungroup()
  medians <- runs |>
    group_by(condition, route_label) |>
    summarise(minutes = median(minutes), .groups = 'drop') |>
    mutate(y = as.numeric(route_label), label = sprintf('%.1f', minutes))
  # The first median carries the word, so that the others can be read as numbers
  first <- which(medians$condition == condition_levels[1] & medians$route_label == route_levels[2])
  medians$label[first] <- paste('median', medians$label[first])
  ggplot(points, aes(minutes, y, colour = route_label, shape = route_label)) +
    geom_segment(aes(x = minutes, xend = minutes, y = y - 0.3, yend = y + 0.3), data = medians,
                 inherit.aes = FALSE, colour = fig_roles[['ink']], linewidth = 0.9) +
    geom_point(size = 2.4) +
    geom_text(aes(minutes, y + 0.46, label = label), data = medians, inherit.aes = FALSE,
              colour = fig_roles[['ink']], size = label_size(base_size)) +
    facet_wrap(~condition, ncol = 1) +
    scale_colour_manual(values = setNames(fig_roles[c('prose', 'script')], route_levels)) +
    scale_shape_manual(values = setNames(c(16, 15), route_levels)) +
    scale_x_log10(breaks = c(0.5, 1, 2, 5, 10, 20, 30), labels = c('0.5', '1', '2', '5', '10', '20', '30')) +
    scale_y_continuous(breaks = 1:2,
                       labels = if (narrow) c('Prose\nreplies', 'Executed\nscripts') else c('Prose replies', 'Executed scripts'),
                       limits = c(0.6, 2.75), expand = expansion(0)) +
    labs(x = if (narrow) 'Minutes to final reply (log scale)'
             else 'Minutes from request to final reply (logarithmic scale)', y = NULL) +
    theme_figure(base_size) +
    theme(axis.text.y = element_text(colour = fig_roles[['ink']], size = rel(0.92)),
          panel.spacing = unit(0.6, 'lines'))
}
stopifnot(range(runs$minutes) > 0.5, range(runs$minutes) < 30)  # within the axis's labels

draw_narrow(fig_minutes(narrow = TRUE), width = 4.4, height = 5.8)
fig_minutes()

A Second Round, at Medium Effort

Specifying the analysis kept the scripts on course. Given the specified request, all 16 scripts of the second round returned the first round’s results, 2 of them to four decimals (Table 4). The open request split them. Three scripts kept the by-word priming slope and returned the first round’s interval, whereas five dropped both random slopes after singular fits and kept random intercepts only. Their Satterthwaite interval, [3.63, 33.41] ms, excludes zero. None of those five tried the model that drops only the by-participant slope. That model gives no singular fit, and lme4 puts the standard deviation of its by-word slope at 29.8 ms, well away from the boundary.

# runs holds the first round's runs and medium_runs the second's, built with
# the hidden checks above, and as_written() is defined before Table 1. A cell
# lists each result as the runs wrote it, with the number of runs that gave
# it, or, with more than max_listed results, the range of each of the three
# numbers. As in Table 1, the first round's replies to the unbalanced data,
# which gave several results, are left to Table 3. The result columns are
# left-aligned, unlike those of Tables 1 and 3, because their cells stack
# several results. Each result, and each range, is joined with non-breaking
# spaces, so that a line break can fall only between them.
max_listed <- 5
results_cell <- function(lines) {
  if (n_distinct(lines) > max_listed) {
    numbers <- t(vapply(lines, result_numbers, character(3), USE.NAMES = FALSE))
    span <- function(j) {
      value <- as.numeric(numbers[, j])
      sprintf('%s&nbsp;to&nbsp;%s', minus(numbers[which.min(value), j]), minus(numbers[which.max(value), j]))
    }
    return(sprintf('%s results, with estimates %s, lower bounds %s and upper bounds %s',
                   stringr::str_to_sentence(count_text(n_distinct(lines))), span(1), span(2), span(3)))
  }
  counts <- table(lines)
  lower <- as.numeric(vapply(names(counts), function(line) result_numbers(line)[2], character(1)))
  counts <- counts[order(-counts, lower)]
  paste(gsub(' ', '&nbsp;', sprintf('%s (%d)', as_written(names(counts)), as.vector(counts)), fixed = TRUE),
        collapse = '<br>')
}
by_cell <- function(rows) {
  rows |>
    group_by(condition, route_label) |>
    summarise(distinct = n_distinct(result_line), results = results_cell(result_line),
              minutes = sprintf('%.1f', median(minutes)), .groups = 'drop')
}
first_cells <- by_cell(runs) |> mutate(results = if_else(distinct == 1, results, 'See Table 3'))

# The stacked results would leave the first round's column too narrow for one
# interval, so this table's shares are fixed too (column_shares() is defined
# before Table 1). The headers break at a <br> and nowhere else, since the
# site keeps a header on one line, and the shares fit the table in the
# 720-pixel column of a desktop page.
column_shares(knitr::kable(
  inner_join(first_cells, by_cell(medium_runs), by = c('condition', 'route_label'),
             suffix = c('_first', '_second')) |>
    arrange(condition, desc(route_label)) |>
    select(condition, route_label, results_first, results_second, minutes_first, minutes_second),
  col.names = c('Condition', 'Route', 'Results,<br>first round', 'Results,<br>second round',
                'Minutes,<br>first round', 'Minutes,<br>second round'),
  caption = sprintf(paste(
    'Results and Times of the Runs in Each Condition and Route, by Round. The first round ran at the effort',
    'level of the session that launched it, and the second at medium effort. Each row covers %s runs in each',
    'round. A result is the estimate and 95%% interval in milliseconds as the runs wrote it, followed by the',
    'number of runs that gave it. Where a cell has more than %s results, the range of each number is',
    'given instead. Times are the median minutes from request to final reply. The %s runs at xhigh effort',
    'are reported in the text.'), count_text(runs_per_cell), count_text(max_listed), count_text(nrow(xhigh_runs))),
  align = c('l', 'l', 'l', 'l', 'r', 'r')), shares = c(17, 10, 23, 28, 11, 11))

None of the second round’s prose replies reproduced a result of the first. All but one of the replies to the open request fell back on the by-participant paired interval, which every reply of the first round had set aside. Given the specified request on the balanced data, most approximated a model without by-participant slopes, and the others left out the by-word slope variance as well. On the unbalanced data, the replies put their estimates near the mean difference over complete participant-word pairs, with lower bounds well above the executed one (Table 4). They were also quicker in every condition than the replies of the first round.

Both runs repeated at xhigh effort returned the first round’s results. The script kept the by-word slope, and the prose reply reconstructed lme4’s boundary fit after 11.0 minutes of reasoning. Two runs cannot establish an effect of the effort level, but they show that the answers at medium effort were not the only ones this model would give.

Figure 3 sets every result of both rounds side by side. Each condition came with the same data in both rounds, yet on the balanced data, no interval of the first round excluded zero, whereas 15 of the 32 intervals of the second round did. Among them is the interval of the five scripts that kept random intercepts only, which shares its estimate with that of the scripts that kept the by-word slope but is narrower.

# runs holds the first round's runs and medium_runs the second's. Each line is one
# distinct result in a condition and route, with the number of runs that reported it.
# Some scripts of the second round wrote four decimals, which round to the first
# round's numbers, so the results are compared at two decimals, as the tables show them.
round_names <- c('First round', 'Second round, medium effort')
results_both <- bind_rows(runs |> mutate(round = round_names[1]),
                          medium_runs |> mutate(round = round_names[2])) |>
  mutate(across(c(estimate, lower, upper), \(x) round(x, 2)),
         round = factor(round, levels = round_names)) |>
  count(round, condition, route_label, estimate, lower, upper, name = 'runs') |>
  group_by(round, condition) |>
  # the scripts' results above the replies', each set from the lowest lower bound
  arrange(desc(route_label), lower, .by_group = TRUE) |>
  mutate(slot = row_number() + 0.5 * (route_label == route_levels[1])) |>
  ungroup()
stopifnot(sum(results_both$runs) == nrow(runs) + nrow(medium_runs))

# The two models of the second round's open-request scripts, named beside their counts
open_two <- results_both$round == round_names[2] & results_both$condition == condition_levels[1] &
  results_both$route_label == route_levels[2]
model_name <- ifelse(results_both$lower == round(intercepts_ci[1], 2), 'random intercepts only',
                     ifelse(results_both$lower == open_script$lower, 'by-word slope kept', NA))
stopifnot(sum(open_two) == 2, !anyNA(model_name[open_two]),
          results_both$runs[open_two] == c(sum(kept_slope), sum(dropped_both))[match(model_name[open_two],
            c('by-word slope kept', 'random intercepts only'))])
# The two lines labelled with their models get a little more room between them
second_open <- results_both$round == round_names[2] & results_both$condition == condition_levels[1]
results_both$slot <- results_both$slot + 0.35 * (second_open & results_both$slot >= 2)
results_both$label <- ifelse(open_two, sprintf('×%d, %s', results_both$runs, model_name),
                             ifelse(results_both$runs > 1, sprintf('×%d', results_both$runs), ''))

fig_results <- function(narrow = FALSE, base_size = if (narrow) 14 else 12.5) {
  lines <- results_both |> mutate(y = -slot)
  if (narrow) {
    lines$panel <- factor(paste(lines$condition, lines$round, sep = '\n'),
                          levels = as.vector(t(outer(condition_levels, round_names, paste, sep = '\n'))))
    lines$label <- sub('random intercepts only', 'intercepts only', lines$label, fixed = TRUE)
    lines$label <- sub('by-word slope kept', 'slope kept', lines$label, fixed = TRUE)
  }
  # The population effect is named once, above the first panel, which gets room for it
  first <- lines$condition == condition_levels[1] & lines$round == round_names[1]
  truth <- lines[first, ][1, ] |>
    mutate(x = population_effect + 1, y = if (narrow) 2.6 else 0.25,
           label = sprintf(if (narrow) 'Population effect, %s ms' else 'Population effect,\n%s ms',
                           sprintf('%.1f', population_effect)))
  room <- lines[first, ][1, ] |> mutate(y = if (narrow) 3.3 else 1.85)
  route_names <- c('Executed scripts', 'Prose replies')
  plot <- ggplot(lines, aes(y = y, colour = route_label, shape = route_label)) +
    geom_vline(xintercept = 0, colour = fig_roles[['rule']], linewidth = 0.5, linetype = 'dashed') +
    geom_vline(xintercept = population_effect, colour = fig_roles[['ink']], linewidth = 0.5, linetype = 'dotted') +
    geom_segment(aes(x = lower, xend = upper, yend = y), linewidth = 1.1, lineend = 'butt') +
    geom_point(aes(x = estimate), size = 2.3) +
    geom_text(aes(x = upper + 1.2, label = label), hjust = 0, colour = fig_roles[['muted']],
              size = label_size(base_size), show.legend = FALSE) +
    geom_text(aes(x = x, y = y, label = label), data = truth, inherit.aes = FALSE, hjust = 0, vjust = 0,
              lineheight = 0.95, colour = fig_roles[['ink']], size = label_size(base_size)) +
    geom_blank(aes(x = population_effect, y = y), data = room, inherit.aes = FALSE) +
    scale_colour_manual(values = setNames(fig_roles[c('prose', 'script')], route_levels),
                        breaks = rev(route_levels), labels = route_names, name = NULL) +
    scale_shape_manual(values = setNames(c(16, 15), route_levels),
                       breaks = rev(route_levels), labels = route_names, name = NULL) +
    # Labels that run past the last panel draw into the margin kept for them
    coord_cartesian(xlim = c(-4, if (narrow) 64 else 50), clip = 'off') +
    scale_x_continuous(breaks = seq(0, 40, 10), expand = expansion(0)) +
    labs(x = 'Reported estimate and 95% interval (ms)', y = NULL) +
    theme_figure(base_size) +
    theme(axis.text.y = element_blank(), legend.position = 'top', legend.justification = 'left',
          legend.location = 'plot', legend.text = element_text(colour = fig_roles[['ink']], size = rel(0.92)),
          legend.margin = margin(0, 0, 4, 0), legend.key.spacing.x = unit(14, 'pt'))
  if (narrow) {
    plot + facet_grid(panel ~ ., scales = 'free_y', space = 'free_y') +
      # The condition in bold, and the round on the line below it
      geom_text(aes(x = -4, y = 1.35, label = sub('\n.*$', '', panel)), data = distinct(lines, panel),
                inherit.aes = FALSE, hjust = 0, vjust = 0, fontface = 'bold', colour = fig_roles[['ink']],
                size = label_size(base_size)) +
      geom_text(aes(x = -4, y = 0.3, label = sub('^.*\n', '', panel)), data = distinct(lines, panel),
                inherit.aes = FALSE, hjust = 0, vjust = 0, colour = fig_roles[['muted']],
                size = label_size(base_size)) +
      scale_y_continuous(expand = expansion(add = c(0.7, 2.2))) +
      theme(strip.text.y = element_blank(), panel.spacing = unit(0.3, 'lines'))
  } else {
    plot + facet_grid(condition ~ round, scales = 'free_y', space = 'free_y', switch = 'y',
                      labeller = labeller(condition = label_wrap_gen(17))) +
      scale_y_continuous(expand = expansion(add = 0.7)) +
      theme(strip.placement = 'outside',
            strip.text.y.left = element_text(angle = 0, hjust = 0, vjust = 1, face = 'bold',
                                             colour = fig_roles[['ink']], size = rel(0.92)),
            panel.spacing.x = unit(1.2, 'lines'), panel.spacing.y = unit(0.8, 'lines'),
            plot.margin = margin(6, 92, 6, 6))
  }
}
stopifnot(min(results_both$lower) > -4, max(results_both$upper) < 50)  # within the axis

draw_narrow(fig_results(narrow = TRUE), width = 4.4, height = 12.5)
fig_results()

Executing the Scripts Again

Executing a script again in the same environment tests whether it is repeatable, and executing it elsewhere tests whether it travels. Every execution below used a fresh R process in the C locale, a standard set of language and character conventions. Checksums compared the output line by line, allowing for the different line endings of Windows and Linux. The scripts use dplyr, emmeans, lme4, lmerTest and tidyr, and some also call pbkrtest for the Kenward–Roger method, which adjusts the standard error as well as the degrees of freedom.

The first round’s 72 stored executions used R 4.3.3 on Ubuntu 24.04.4, with lme4 1.1-35.1, lmerTest 3.1-3 and Matrix 1.6-5, and every script printed identical output each time. A Dockerfile rebuilds that environment from a dated snapshot of the software, and on a fresh machine on GitHub, all 24 scripts reproduced their stored output, as verification.csv records.

The second round’s 25 scripts were each executed twice in the library their agents had used, under R 4.6.1 on Windows. They too printed identical output, ending with the result each agent had reported.

# run_script_with_retries(), data_file_for() and child_environment() come from
# experiment/run_script.R. The list of scripts (scripts, n_scripts) and their
# stored executions (stored_of()) are read from design.csv and executions.csv
# earlier on the page.
source(file.path(experiment, 'run_script.R'))

# Every script executed once more, and again if its R process crashes, which
# happens occasionally on Windows
rerun_attempts <- bind_rows(lapply(seq_len(n_scripts), function(i)
  run_script_with_retries(scripts$run_id[i], data_file_for(scripts$data[i]),
                          experiment_dir = normalizePath(experiment), label = 'page')))
rerun <- rerun_attempts |> group_by(run_id) |> filter(attempt == max(attempt)) |> ungroup()
rerun <- rerun[match(scripts$run_id, rerun$run_id), ]
page_environment <- child_environment()

# Rscript exits with status 1 when a script stops with an error. Any other
# nonzero status means that the R process itself ended abnormally, as when it
# crashes, which says nothing about the script.
rerun$outcome <- ifelse(rerun$status == 0, 'completed', ifelse(rerun$status == 1, 'error', 'crash'))
rerun$same_result <- rerun$outcome == 'completed' & !is.na(rerun$result) &
  rerun$result == stored_of(rerun$run_id, 'result')
rerun$same_output <- rerun$outcome == 'completed' &
  rerun$stdout_sha256 == stored_of(rerun$run_id, 'stdout_sha256')

c(stored_executions = nrow(executions), scripts = n_scripts,
  completed_now = sum(rerun$outcome == 'completed'),
  completed_after_a_crash = sum(rerun$outcome == 'completed' & rerun$attempt > 1),
  stopped_with_error_now = sum(rerun$outcome == 'error'), crashed_now = sum(rerun$outcome == 'crash'),
  same_result_now = sum(rerun$same_result), same_output_now = sum(rerun$same_output))

#>       stored_executions                 scripts           completed_now 
#>                      72                      24                      24 
#> completed_after_a_crash  stopped_with_error_now             crashed_now 
#>                       2                       0                       0 
#>         same_result_now         same_output_now 
#>                      24                       0

This page executes each of the first round’s scripts once more whenever it is knitted, under the versions of R and its packages recorded in the site’s lockfile, renv.lock, which are newer than the first round’s. All 24 scripts ran to completion and printed the same result line as the stored executions, but none printed identical output. In 2 of them, the R process crashed at first and the script completed when executed again.

To find the source of the difference, I also executed each script on one Windows 11 machine with three libraries that differed only in their versions of lme4 and lmerTest (Table 5). Matrix stayed at 1.7-6 in all three, since the stored version could not be installed under the newer R. With the earlier lme4, both versions of lmerTest reproduced the stored output, as the Docker rebuild had. The newer lme4 changed only the formatting. It added a trailing space to the lines of some random-effects tables, renamed a column from deviance to -2*log(L) and, in some tables, swapped the columns of variances and standard deviations. None of the numbers changed.

labels <- c(`renv.lock` = 'Newer (renv.lock)', `lme4 1.1-35.1` = 'Newer, with lme4 as stored',
            `lme4 1.1-35.1 and lmerTest 3.1-3` = 'Newer, with lme4 and lmerTest as stored')
# count_of() writes 'k of n', as in Table 2
row_of <- function(label, system, versions, same_output, same_result) {
  data.frame(Environment = label, `Operating system` = system, versions[, c('R', 'lme4', 'lmerTest', 'Matrix')],
             `Identical output` = count_of(sum(same_output), length(same_output)),
             `Same result line` = count_of(sum(same_result), length(same_result)), check.names = FALSE)
}
environment_table <- bind_rows(
  row_of('Stored, rebuilt with Docker', rebuild_release, versions_of(unique(verification$environment)),
         verification$identical, verification$result == stored_of(verification$run_id, 'result')),
  bind_rows(lapply(names(labels), function(library) {
    rows <- library_rows(library)
    row_of(labels[[library]], newer_system, library_versions(library), rows$same_output, rows$same_result)
  })))

knitr::kable(environment_table,
  caption = sprintf(paste(
    'Output of the First Round\'s Scripts in %s Environments, Compared With Their Stored Executions.',
    'Each script was executed once in each environment. \'Identical output\' counts the scripts whose',
    'output matched the stored output line for line, and \'Same result line\' counts those whose result',
    'line matched.'),
    stringr::str_to_title(count_text(nrow(environment_table)))),
  align = c('l', 'l', 'r', 'r', 'r', 'r', 'r', 'r'))

The R session that last knitted this page is shown below, although the scripts ran in separate R processes in the C locale.

sessionInfo()

#> R version 4.6.1 (2026-06-24 ucrt)
#> Platform: x86_64-w64-mingw32/x64
#> Running under: Windows 11 x64 (build 26300)
#> 
#> Matrix products: default
#>   LAPACK version 3.12.1
#> 
#> locale:
#> [1] LC_COLLATE=English_United Kingdom.utf8  LC_CTYPE=English_United Kingdom.utf8   
#> [3] LC_MONETARY=English_United Kingdom.utf8 LC_NUMERIC=C                           
#> [5] LC_TIME=English_United Kingdom.utf8    
#> 
#> time zone: Europe/London
#> tzcode source: internal
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices datasets  utils     methods   base     
#> 
#> other attached packages:
#> [1] lme4_2.0-6    Matrix_1.7-6  pilotr_0.3.0  depictr_0.3.0 ggplot2_4.0.3 dplyr_1.2.1  
#> 
#> loaded via a namespace (and not attached):
#>  [1] sass_0.4.10         generics_0.1.4      renv_1.2.3          stringi_1.8.7      
#>  [5] blogdown_1.24       lattice_0.22-9      digest_0.6.39       magrittr_2.0.5     
#>  [9] evaluate_1.0.5      grid_4.6.1          RColorBrewer_1.1-3  bookdown_0.47      
#> [13] fastmap_1.2.0       jsonlite_2.0.0      processx_3.9.0      scales_1.4.0       
#> [17] numDeriv_2016.8-1.1 textshaping_1.0.5   jquerylib_0.1.4     reformulas_0.4.4   
#> [21] Rdpack_2.6.6        cli_3.6.6           rlang_1.3.0         rbibutils_2.4.1    
#> [25] splines_4.6.1       withr_3.0.3         cachem_1.1.0        yaml_2.3.12        
#> [29] otel_0.2.0          tools_4.6.1         uuid_1.2-2          nloptr_2.2.1       
#> [33] minqa_1.2.8         boot_1.3-32         vctrs_0.7.3         R6_2.6.1           
#> [37] lifecycle_1.0.5     stringr_1.6.0       MASS_7.3-66         callr_3.8.0        
#> [41] pkgconfig_2.0.3     xaringanExtra_0.8.0 pillar_1.11.1       bslib_0.11.0       
#> [45] lmerTest_3.2-1      gtable_0.3.6        Rcpp_1.1.2          glue_1.8.1         
#> [49] systemfonts_1.3.2   xfun_0.60           tibble_3.3.1        tidyselect_1.2.1   
#> [53] knitr_1.51          farver_2.1.2        htmltools_0.5.9     nlme_3.1-170       
#> [57] labeling_0.4.3      rmarkdown_2.31      svglite_2.2.2       compiler_4.6.1     
#> [61] S7_0.2.2

What These Runs Can Establish

The experiment used one language model and one simulated data set, clean and fully described in the prompt, on which the lower bound of the interval lay close to zero. Real data bring exclusions, coding errors and gaps in documentation that these agents never faced. I designed the follow-up after seeing the results on the balanced data. The specified request also supplied the model formula and the function for the interval, so the agreement among its scripts tested implementation more than analysis. The prose route imitated a chat without a code interpreter, although many chat assistants now have one. An assistant that runs code and lets the researcher keep it could close much of the gap with a coding agent.

The two rounds differed in date, in the version of Claude Code, in the operating system and in the versions of R and lme4 that the agents worked with. The first round’s effort level, that of the session that launched it, went unrecorded, and only two runs were repeated at xhigh. The differences between the rounds therefore cannot be put down to the effort level alone. With eight runs in each combination of condition and route, the counts describe these runs and do not estimate error rates. Finally, I wrote every prompt as someone who reads R, so the experiment cannot show what a researcher without R would achieve with the same agents.

What Executing the Code Added

No text in this experiment was repeated, yet every executed script printed the same output each time it ran. Because the versions of R, lme4, lmerTest and Matrix were recorded and can still be installed, the first round’s scripts printed that output on another machine too. Under newer versions, their output changed only in its formatting. The scripts also left a record of the calculation that no prose reply could supply. Even where a reply matched the scripts, only recomputing its quantities could tell a correct reply from a confident error.

A reproducible analysis can still be a poor one, however (Leek & Peng, 2015). The five scripts of the second round that dropped both random slopes printed identical output each time they ran. Their interval excludes zero and points the right way, because the effect is real in this simulation. With real data, nobody could know that. Leaving out a by-word slope when the effect varies across words makes a test too liberal (Barr et al., 2013), so the same analysis would claim an effect more readily where none existed.

Agreement between runs is a separate matter. Many-analysts studies, which give one data set and one question to many teams, find wide variation among human analysts (Breznau et al., 2022; Silberzahn et al., 2018). Autonomous AI analysts given one data set and one hypothesis also reached divergent verdicts that shifted with the language model and the prompt (Bertran et al., 2026). Alizadeh et al. (2026) ran Claude Code and Codex 20 times each on one of those many-analysts studies (Breznau et al., 2022). The agents varied their model specifications at least as much as the human teams had, yet their estimates stayed broadly in line with the human consensus. The authors located the agents’ vulnerability in interpretation more than in estimation.

The task in this experiment was far simpler, and the close agreement of the first round may reflect a shared preference as much as independent confirmation, since runs of one model are not independent analysts. When Herdağdelen et al. (2026) left Claude Opus 4.8 to make the five decisions of an analysis with no preference assigned to it, it made the same five choices in 112 of 123 runs. Forty-eight runs of one model that agree may amount to one opinion stated 48 times.

Bertran et al. (2026) propose giving AI analysts a published specification, to reveal how much of their disagreement comes from underspecified choices. Here, the specified request kept the scripts consistent across rounds, whereas the open request divided them over one random slope, which decided whether their interval excluded zero. That consistency belonged to one model, at one effort level, on one structured task, and it may not survive the model’s next update (L. Chen et al., 2024). What can be kept is the code and a record of its environment, and together they leave any such difference open to inspection.

Where Coding Literacy Comes In

By coding literacy, I mean being able to read and write code well enough to tell what a script does. The experiment did not vary it, but it shows where literacy came into the work, and at those points, much of what was needed turned out to be statistical knowledge.

Writing the specified request required knowing that the data call for crossed random effects, that the choice of random slopes changes the interval (Barr et al., 2013) and that transforming reaction times can change the effects being tested (Lo & Andrews, 2015). Coding came in only where the request named the function that fits the model, wrote the model in lme4’s formula notation and gave the call to confint() that computes the interval. I cannot say how much that helped, since no condition stated the same analysis in plain statistical language.

The same request shows the limits of a specification. It fixed a Wald interval, which can be too narrow with few participants or items (Luke, 2017), and I chose it because one function call defines it exactly and so makes the scripts comparable. On the unbalanced data, that choice decided whether the interval excluded zero. An analysis for publication should use Satterthwaite or Kenward–Roger degrees of freedom, as the open-request scripts did, yet none of the 65 runs given the specified request, in either round, mentioned either. On this point, the agents left to choose did better than my specification.

Reading the code is the second point of entry. The scripts were only 18 to 252 lines long, short enough for a researcher who reads R to confirm in a few minutes that each fits the intended model to the intended rows. In the second round, the summaries of three of the five scripts that dropped both slopes justified doing so by calling their variances zero or nearly so. A reader of R would notice that none of these scripts tried the model that drops only the by-participant slope. Printing the fits they discarded would also show that every model with the by-word slope put its standard deviation at about 30 ms. Without that reading, the researcher is left with the agent’s summary of its own work, which is more text from the same language model. Even so, what a reader would catch here is a statistical choice written in code, since no result pointed to a coding bug.

A reader of R is not the only possible remedy. Given the open request at xhigh effort, the one script I ran kept the by-word slope with nobody reading it. Three of the eight at medium effort did too, however, so this is only a hint that a setting may do part of a reader’s work. A second agent asked to review each script might do more, although runs of one model may share its preferences, and different models often err in the same way (Kim et al., 2025). The experiment did not test this.

Studies of how people use these tools complicate the picture. Anthropic, the maker of Claude Code, analysed about 400,000 sessions with it and concluded that the ability to steer the agent comes more from command of a domain than from the ability to write code (Hitzig et al., 2026). The authors call these findings preliminary, since the ratings come from a model’s reading of the transcripts and they cannot tell whether the code was used afterwards. In a randomised trial, 487 consultants at the firm BCG attempted data-science tasks outside their expertise, and those given ChatGPT outscored those given other resources. Among them, self-described competent coders came within about 5 percentage points of the firm’s data scientists on the coding and statistics tasks, whereas the group as a whole stayed 12 to 14 points below them (Wiles et al., 2026). Coding experience was not randomised, however, and the tool was a chat assistant.

Fluency in code does not guarantee independent judgement either. In another study from Anthropic, 50 professional analysts who regularly used R or Python carried out an analysis with Claude Code, whose preference on each decision had been secretly set by a coin flip. On four decisions between two reasonable options, they went with the agent’s preference 86% of the time, and in an exploratory analysis, neither experience nor expertise predicted who did so (Herdağdelen et al., 2026). The preferences were strengthened on purpose, so these figures may overstate everyday influence. The small group of less experienced analysts who chose to work without the agent made more serious errors, but because they chose for themselves, the comparison cannot isolate what the agent added.

Understanding is a further concern. In a small preprint, 54 computing students built a website, and those given a coding agent produced a better first version than those given a chatbot but understood it less well. With the agent, their background in code no longer predicted the quality of that version, but it still predicted their comprehension (Balepur et al., 2026). An agent may therefore mask differences in literacy without removing them. Developers who learned a software package with AI assistance likewise understood it less well, and they were not significantly faster (Shen & Tamkin, 2026).

Research on automation explains why this matters. Automation leaves its operators to monitor the system and to take over when it fails, while the skills needed for both decay without practice (Bainbridge, 1983). In a meta-analysis of 18 experiments, the more of a task that automation took over, the better performance was while it worked and the worse when it failed (Onnasch et al., 2014). Complacency towards automation arises in experts as well as novices (Parasuraman & Manzey, 2010). Knowledge workers more confident in generative AI also reported less critical thinking (Lee et al., 2025), and Messeri and Crockett (2024) warn that AI tools can leave scientists believing they understand more than they do. As agents improve, they give their supervisors fewer errors to catch, and so less practice at catching the ones that remain.

The evidence on learning is not uniform, since a meta-analysis of programming education found no reliable effect on exam scores, in either direction, when students could not use the tools in the exam (Maier et al., 2026). What the website and package studies show is that the quality of the output and the understanding behind it need assessing separately. For researchers in training, I would still advise learning to read and write R or Python, which for now is the most direct way to check an agent’s work. Study design, statistics and the habit of asking how a result could be wrong deserve at least as much of their time.

What Experts Expect From Coding Agents

This experiment captures one model on 29 September 2026 and 11 October 2026, and each of its runs wrote one short script. Longer projects add a problem that these scripts did not test, the upkeep of a growing codebase, from its documentation to its stale and duplicated files. Carlini (2026), who built a 100,000-line C compiler with 16 parallel agents, noted that ‘LLM-written code frequently re-implements existing functionality’ and assigned one agent to coalesce duplicate code and others to improve code quality and documentation. In SlopCodeBench, a benchmark in which agents repeatedly extend their own code, the code’s structure eroded in 77% of the runs, and quality guidance repeated at every step improved the starting point without slowing the decline (Orlanski et al., 2026). In open-source projects that adopted Cursor, an AI code editor, development sped up sharply but briefly, whereas code complexity and the warnings of static analysis, which inspects code without running it, rose persistently (He et al., 2026).

The upkeep is improving, if unevenly. In SWE-CI, a benchmark in which agents maintain real repositories over many rounds of change, each provider’s newer models scored higher, although most models still broke previously passing tests in over three quarters of the tasks (J. Chen et al., 2026). In a research note from METR (Model Evaluation and Threat Research), a nonprofit that studies the capabilities and risks of advanced AI systems, maintainers reviewing patches from successive models saw the step from Claude 4 Opus to Claude 4.5 Sonnet improve mostly code quality (Whitfill et al., 2026).

Experts have not forecast how long researchers will need to read code to check an agent’s work, but their forecasts speak to a narrower milestone (Figure 4). In spring 2026, the Forecasting Research Institute asked its Longitudinal Expert AI Panel when a model would first succeed 80% of the time at tasks that take a human expert at least eight hours. The panel’s experts are computer scientists, economists, technical staff at AI companies and policy researchers. Its superforecasters are people who forecast with high accuracy in earlier geopolitical forecasting tournaments (Forecasting Research Institute, 2025). The experts’ median forecast gave that milestone an even chance by 2030, and the superforecasters’ by 2028 (Forecasting Research Institute, 2026). Those dates concern completing tasks, however, and say nothing directly about when agents will keep a codebase healthy without supervision, or when researchers can stop reading code.

# forecast_dates, read before the section's text, holds each panel's year for each
# probability. The bar joins a panel's years for a 25% and a 75% chance, which are
# group summaries, not a range of respondents' forecasts. Every dot carries its year, so
# the figure needs no axis for them.
fig_forecasts <- function(narrow = FALSE, base_size = if (narrow) 14 else 12.5) {
  dates <- forecast_dates |>
    mutate(panel = factor(panel, levels = c('Superforecasters', 'Experts')),
           probability_label = sprintf('%.0f%%', 100 * probability),
           year_label = sprintf('%.0f', year))
  spans <- dates |> summarise(first = min(year), last = max(year), .by = panel)
  ggplot(dates, aes(year, panel)) +
    geom_segment(aes(x = first, xend = last, yend = panel), data = spans, colour = fig_roles[['grid']],
                 linewidth = 4.5, lineend = 'round') +
    geom_point(aes(size = probability == 0.5), colour = fig_roles[['ink']]) +
    # On a phone, a year is narrower than a label, so neighbouring labels lean apart
    geom_text(aes(label = year_label,
                  hjust = if (narrow) c(0.85, 0.35, 0.2)[match(probability, c(0.25, 0.5, 0.75))] else 0.5),
              nudge_y = -0.3, colour = fig_roles[['muted']], size = label_size(base_size)) +
    geom_text(aes(label = probability_label), data = filter(dates, panel == 'Experts'), nudge_y = 0.3,
              colour = fig_roles[['ink']], size = label_size(base_size)) +
    scale_size_manual(values = c(`FALSE` = 2.2, `TRUE` = 3.6)) +
    scale_x_continuous(breaks = NULL, limits = c(if (narrow) 2025.6 else 2026, 2036.4),
                       expand = expansion(0)) +
    coord_cartesian(clip = 'off') +
    scale_y_discrete(expand = expansion(add = c(0.48, 0.5))) +
    labs(x = NULL, y = NULL) +
    theme_figure(base_size) +
    theme(axis.text.y = element_text(colour = fig_roles[['ink']], size = rel(0.92)))
}
stopifnot(range(forecast_dates$year) >= 2026, range(forecast_dates$year) <= 2036)

draw_narrow(fig_forecasts(narrow = TRUE), width = 4.4, height = 1.75)
fig_forecasts()

Measurements of past progress point the same way. By METR’s measure, the length of software task that frontier models, the most capable of their day, can complete half the time doubled about every seven months between 2019 and 2025 (Kwa et al., 2025), and its update of January 2026 puts the doubling time since 2023 at about four months (METR, 2026a). That measure has limits of its own. It rests on self-contained tasks ‘with clear success criteria that can be automatically evaluated’ (METR, 2026b), and by METR’s account, relatively few of its tasks remain beyond the latest models (METR, 2026a). The health of a codebase, by contrast, can be scored only in part. Agents resolved fewer tasks when they built on code written by other agents than on code written by people, and many traditional measures of maintainability did not explain the difference (Patel et al., 2026).

Capability need not turn into use, either. The computer scientists Arvind Narayanan and Sayash Kapoor (2025) distinguish gains in capability from the development and adoption of reliable applications. They expect diffusion to stay slow in consequential settings, held back by institutions, tacit knowledge and the need for validation, and a growing share of people’s work to go into controlling AI. Upkeep illustrates their point, since whether a file is stale may depend on what its owner intends. Asked to ‘clean up old branches’, one agent issued a command to delete the branches on the remote server, although the user may have meant only those on their own computer (Hughes, 2026).

Whichever way these forecasts turn out, the case for saved and executed code would stand, since a script can be executed again whether or not anyone reads it. A well-kept codebase can still implement a poor analysis, and whether an effect should be allowed to vary across words depends on what a study means to generalise to, which no amount of tidying settles. Someone also has to answer for the result, which is why Nature decided in 2023 to accept no language model as a credited author (“Tools Such as ChatGPT,” 2023). In the meantime, the practical question is how to check the work that agents produce.

A Workflow for Research With Coding Agents

The experiment and the wider evidence point to a workflow built on the three choices reviewed at the start, the tools, the requests and the record, with checking and upkeep between them.

The tools come first, and their choice deserves a reason. I used a proprietary model because its behaviour was the object of study, a case in which McLoughlin et al. (2026) suggest that the choice may be reasonable (see also Palmer et al., 2024). Where a model’s replies are themselves the data, as in text annotation, a model whose version can be fixed, such as an open model run locally, gives a far stronger record (Spirling, 2023). With real data, the terms of consent and data protection also need checking before an agent sends anything it reads to a hosted model.

An analysis plan written before any code, ideally in a preregistration, allows every request to be a specified one, and the request can also ask for a runnable script, pinned dependencies and tests. Where a decision remains open, it may help to form a view before seeing the agent’s preference, or to ask it for the best case for each option, two approaches that Herdağdelen et al. (2026) suggest exploring.

Checking falls, for now, to someone who reads each change before accepting it and executes each script again in a fresh R process. That includes what the code does to files, since 3 of the 24 script agents of the first round wrote outside their folders. A growing project soon outruns any reader, and the checking then has to rest on independent reanalysis and on tests with many simulated data sets whose answers are known (Morris et al., 2019). One data set would not have been enough here, since both intervals of the second round’s open-request scripts contain the population effect (Figure 3), and only many data sets would show how often each misses it. Tests alone are not enough either. Maintainers of three open-source projects would not have merged roughly half of the patches by agents of 2024 and 2025 that had passed a benchmark’s tests, although those agents had no chance to respond to feedback (Whitfill et al., 2026).

The upkeep belongs in that review. After each substantial change, the documentation, stale files, inconsistencies and duplicated code all need a look. An agent can take on much of this work, as in Carlini’s project (Carlini, 2026), but asking it to may not be enough, since repeated guidance did not slow the decline in SlopCodeBench (Orlanski et al., 2026).

Finally, the record should include the effort level, alongside the versions of the model and the agentic tool, the settings, the prompts, the replies, the tool calls, the scripts and their outputs, and the software environment (Baltes et al., 2026; Feuerriegel et al., 2026). Having left the first round’s effort level unrecorded, I can only hint at its part in the differences between the rounds.

Used in this way, coding agents can make an AI-assisted analysis far more reproducible than a reply written without running code. Their summaries of their own work, however, are more text from the same model, so for now, the checking falls to a researcher who can specify, read and execute the code. How long that will hold is the least certain part of this post. That someone will have to answer for the result is among the most certain.