library(here) # Working directory management
library(fs) # File & folder manipulation
library(pipebind) # Piping goodies
library(archive) # Memory efficient object storage
library(cmdstanr) # Lightweight R interface for Stan
library(posterior) # Wrangling Stan models' output
library(data.table) # Fast data manipulation (in-RAM)
library(duckdb) # DuckDB R interface
library(dplyr) # Manipulating data.frames - core (Tidyverse)
library(tidyr) # Manipulating data.frames - extras (Tidyverse)
library(dbplyr) # DB/SQL backend for dplyr/tidyr (Tidyverse)
library(stringr) # Manipulating strings (Tidyverse)
library(purrr) # Manipulating lists (Tidyverse)
library(lubridate) # Manipulating date/time (Tidyverse)
library(ggplot2) # Best plotting library (Tidyverse)
library(plotly) # Interactive plots
library(bayesplot) # PPC/Diagnostic plots for Stan models
library(patchwork) # Combining plots
options(
mc.cores = max(1L, parallel::detectCores(logical = TRUE)),
scipen = 999L,
digits = 4L,
ggplot2.discrete.colour = \() scale_color_viridis_d(),
ggplot2.discrete.fill = \() scale_fill_viridis_d()
)
nrows_print <- 10
setDTthreads(parallel::detectCores(logical = FALSE))
Caution
By design, this post contains very few explanations.
Its goal was simply to translate Ethan’s blog post to R, and update his Stan code to use within-chain parallelization and compiler optimizations, for faster sampling.
Feel free to read the original blog post to better understand what the code is doing.
Tip
You can check the page’s source code by clicking on the </> Code button at the top-right.
Setup
NoteQuarto/knitr setup
Adding custom knitr hooks
#---------------------------#
####🔺knitr custom hooks ####
#---------------------------#
library(knitr)
## Adding the `time_it` code chunk option
knitr::knit_hooks$set(
time_it = local({
assign("TIMES", list(), .GlobalEnv)
start <- NULL
function(before, options) {
if (before) {
start <<- Sys.time()
} else {
TIMES[[options$label]] <<- difftime(Sys.time(), start)
}
}
})
)Applying a custom theme to all ggplot objects, both light and dark versions
#---------------------------#
####🔺ggplot knit_prints ####
#---------------------------#
library(knitr)
library(ggplot2)
## Inspired by: https://debruine.github.io/quarto_demo/dark_mode.html
knit_print.ggplot <- function(x, options, ...) {
if (any(grepl("patchwork", class(x)))) {
plot_dark <- x & dark_addon_mar
plot_light <- x & light_addon_mar
} else {
plot_dark <- x + dark_addon_mar
plot_light <- x + light_addon_mar
}
cat('\n<div class="light-mode">\n')
print(plot_light)
cat('</div>\n')
cat('<div class="dark-mode">\n')
print(plot_dark)
cat('</div>\n\n')
}
registerS3method("knit_print", "ggplot", knit_print.ggplot)Applying a custom theme to all gt tables
#-----------------------#
####🔺gt knit_prints ####
#-----------------------#
library(knitr)
library(gt)
knit_print.grouped_df <- function(x, options, ...) {
if ("grouped_df" %in% class(x)) {
x <- ungroup(x)
}
cl <- intersect(class(x), c("data.table", "data.frame"))[1]
nrows <- ifelse(
!is.null(options$total_rows),
as.numeric(options$total_rows),
dim(x)[1]
)
is_open <- ifelse(
!is.null(options[["details-open"]]),
as.logical(options[["details-open"]]),
FALSE
)
cat(str_glue("\n<details{ifelse(is_open, ' open', '')}>\n"))
cat("<summary>\n")
cat(str_glue("\n*{cl} [{scales::label_comma()(nrows)} x {dim(x)[2]}]*\n"))
cat("</summary>\n<br>\n")
print(gt::as_raw_html(style_table(x, nrows)))
cat("</details>\n\n")
}
registerS3method("knit_print", "grouped_df", knit_print.grouped_df)
knit_print.data.frame <- function(x, options, ...) {
cl <- intersect(class(x), c("data.table", "data.frame"))[1]
nrows <- ifelse(
!is.null(options$total_rows),
as.numeric(options$total_rows),
dim(x)[1]
)
is_open <- ifelse(
!is.null(options[["details-open"]]),
as.logical(options[["details-open"]]),
FALSE
)
cat(str_glue("\n<details{ifelse(is_open, ' open', '')}>\n"))
cat("<summary>\n")
cat(str_glue("\n*{cl} [{scales::label_comma()(nrows)} x {dim(x)[2]}]*\n"))
cat("</summary>\n<br>\n")
print(gt::as_raw_html(style_table(x, nrows)))
cat("</details>\n\n")
}
registerS3method("knit_print", "data.frame", knit_print.data.frame)Functions for interactive data presentation
#----------------------------------#
####🔺knitr interactive display ####
#----------------------------------#
library(htmltools)
library(reactable)
## Getting list to display nicely in rendered documents
make_list_reactable <- function(list_dat) {
list_name <- deparse(substitute(list_dat))
get_list_elt_dim <- function(elt) {
list_elt <- list_dat[[elt]]
list_elt_dim <- if (any(c("data.frame", "matrix") %in% class(list_elt))) {
dim(list_elt)
} else {
length(list_elt)
}
return(paste0(list_elt_dim, collapse = ", "))
}
dat <- data.frame(names(list_dat)) |>
set_names(list_name) |>
mutate(
Type = unlist(pick(list_name)) |>
map_chr(\(x) class(list_dat[[x]]) |> paste0(collapse = ", ")),
Dimensions = unlist(pick(list_name)) |> map_chr(get_list_elt_dim)
)
get_list_details <- function(dat, idx, max_print = 200, max_digits = 3) {
Element <- dat[[idx]]
style <- "padding: 0.5rem"
if (any(c("data.frame", "matrix") %in% class(Element))) {
reactable(
data.frame(Element),
outlined = TRUE,
striped = TRUE,
highlight = TRUE,
compact = TRUE
) |>
htmltools::div(style = style)
} else if ("list" %in% class(Element)) {
make_list_reactable(Element)
} else if (length(Element) > max_print) {
htmltools::div(
htmltools::p(
head(Element, max_print) |>
round(max_digits) |>
paste0(collapse = ", ") |>
paste("...", sep = ", ")
),
htmltools::p(
stringr::str_glue(
"[ omitted {length(Element) - max_print} entries ]"
),
style = "font-style: italic"
),
style = style
)
} else {
htmltools::div(
round(Element, max_digits) |> paste0(collapse = ", "),
style = style
)
}
}
reactable(
dat,
defaultColDef = colDef(vAlign = "center", headerVAlign = "center"),
details = \(idx) get_list_details(list_dat, idx),
outlined = TRUE,
striped = TRUE,
highlight = TRUE,
compact = FALSE,
fullWidth = TRUE,
defaultPageSize = 15
)
}
NoteStan setup
Installing CmdStan
## Skip this step if CmdStan is already installed
cmdstanr::check_cmdstan_toolchain(fix = TRUE, quiet = TRUE)
cpp_opts <- list(
stan_threads = TRUE,
STAN_CPP_OPTIMS = TRUE,
STAN_NO_RANGE_CHECKS = TRUE, # WARN: remove this if you haven't tested the model
PRECOMPILED_HEADERS = TRUE,
# , CXXFLAGS_OPTIM = "-march=native -mtune=native"
CXXFLAGS_OPTIM_TBB = "-mtune=native -march=native",
CXXFLAGS_OPTIM_SUNDIALS = "-mtune=native -march=native"
)
cmdstanr::install_cmdstan(cpp_options = cpp_opts, quiet = TRUE)Loading CmdStan (if already installed)
highest_cmdstan_version <- fs::dir_ls(Sys.getenv("CMDSTAN_PATH")) |>
fs::path_file() |>
purrr::keep(\(e) stringr::str_detect(e, "cmdstan-")) |>
stringr::str_split_i('-', 2) |>
purrr::reduce(\(x, y) ifelse(utils::compareVersion(x, y) == 1, x, y))
set_cmdstan_path(str_glue("{Sys.getenv('CMDSTAN_PATH')}cmdstan-{highest_cmdstan_version}"))Setting up knitr’s engine for CmdStan
## Inspired by: https://mpopov.com/blog/2020/07/30/replacing-the-knitr-engine-for-stan/
## Note: We could haved use cmdstanr::register_knitr_engine(),
## but it wouldn't include compiler optimizations & multi-threading by default
knitr::knit_engines$set(
cmdstan = function(options) {
output_var <- options$output.var
if (!is.character(output_var) || length(output_var) != 1L) {
stop(
"The chunk option output.var must be a character string ",
"providing a name for the returned `CmdStanModel` object."
)
}
if (options$eval) {
if (options$cache) {
cache_path <- options$cache.path
if (length(cache_path) == 0L || is.na(cache_path) || cache_path == "NA") {
cache_path <- ""
}
dir <- paste0(cache_path, options$label)
} else {
dir <- tempdir()
}
file <- cmdstanr::write_stan_file(
options$code,
dir = dir,
force_overwrite = TRUE
)
mod <- cmdstanr::cmdstan_model(
stan_file = file,
cpp_options = list(
stan_threads = TRUE,
STAN_CPP_OPTIMS = TRUE,
STAN_NO_RANGE_CHECKS = TRUE, # The model was already tested
PRECOMPILED_HEADERS = TRUE,
# , CXXFLAGS_OPTIM = "-march=native -mtune=native"
CXXFLAGS_OPTIM_TBB = "-mtune=native -march=native",
CXXFLAGS_OPTIM_SUNDIALS = "-mtune=native -march=native"
),
stanc_options = list("O1"),
force_recompile = TRUE
)
assign(output_var, mod, envir = knitr::knit_global())
}
options$engine <- "stan"
code <- paste(options$code, collapse = "\n")
knitr::engine_output(options, code, '')
}
)1 Data
1.1 Extracting the data
Connecting to the .sqlite DB (using DuckDB instead of SQLite):
INSTALL sqlite;
LOAD sqlite;
CALL sqlite_attach(?db_path);dbplyr automatically translates dplyr/tidyr code into SQL !
(list(
## Table 1: ascent
tbl(con, "ascent") |>
filter(country %like% "USA") |>
mutate(
route_id = str_c(
str_replace_all(crag, ' ', '_'),
"__",
str_replace_all(name, ' ', '_'),
"__",
if_else(climb_type == 1, 'boulder', 'rope')
),
ascent_date = to_timestamp(date)
) |>
select(user_id, route_id, climb_type, grade_id, method_id, ascent_date),
## Table 2: grade
tbl(con, "grade") |> select(grade_id = id, usa_routes, usa_boulders),
## Table 3: method
tbl(con, "method") |> select(method_id = id, method_name = name)
) |>
reduce(left_join) |>
select(-grade_id, -method_id) |>
compute("climbs", overwrite = TRUE))Time difference of 0.7579 secs
SELECT
ascent.user_id
, REPLACE(ascent.crag, ' ', '_')
|| '__' || REPLACE(ascent.name, ' ', '_')
|| '__' || CASE WHEN ascent.climb_type = 1 THEN 'boulder' ELSE 'rope' END
AS route_id
, ascent.climb_type as climb_type
, to_timestamp(ascent.date) AS ascent_date
, grade.usa_routes
, grade.usa_boulders
, method.name AS method_name
FROM ascent
JOIN grade ON grade.id = ascent.grade_id
JOIN method ON method.id = ascent.method_id
WHERE ascent.country = 'USA'Time difference of 0.6858 secs
data.frame [658,822 x 7]
| [ omitted 658,807 entries ] |
1.2 Processing the data
threshold_ascents_dt
climbs_clean <- climbs[
order(user_id, route_id, ascent_date, usa_routes, usa_boulders, method_name),
.SD[1],
by = .(user_id, route_id)
][, let(
usa_boulders = factor(usa_boulders, levels = bouldering_grades),
usa_routes = factor(usa_routes, levels = route_ratings),
label = as.integer(method_name %chin% c("Onsight", "Flash"))
)][,
let(route_rating = mode_cat(usa_routes), bouldering_grade = mode_cat(usa_boulders)),
by = route_id
][, threshold_ascents_dt(.SD)][, let(
route_idx = frank(route_id, ties.method = "dense"),
user_idx = frank(user_id, ties.method = "dense")
)][
order(route_idx, user_idx),
.(route_idx, route_id, user_idx, climb_type, ascent_date, route_rating, bouldering_grade, label)
]Time difference of 13.03 secs
TipAlternative
climbs_first <- climbs[
order(user_id, route_id, ascent_date, usa_routes, usa_boulders, method_name),
.SD[1],
by = .(user_id, route_id)
]
climbs_clean <- copy(climbs_first)[,
threshold_ascents_dt(.SD)
][
## Replacing all route_ratings for a given route_id by its mode
climbs_first[
usa_boulders %in% bouldering_grades
][, .(bouldering_grade = mode_cat(usa_boulders)), by = route_id],
c("route_id", "bouldering_grade") := list(i.route_id, i.bouldering_grade),
on = "route_id"
][
## Replacing all bouldering_grades for a given route_id by its mode
climbs_first[
usa_routes %in% route_ratings
][, .(route_rating = mode_cat(usa_routes)), by = route_id],
c("route_id", "route_rating") := list(i.route_id, i.route_rating),
on = "route_id"
][, let(
route_idx = frank(route_id, ties.method = "dense"),
user_idx = frank(user_id, ties.method = "dense"),
label = as.integer(method_name %chin% c("Onsight", "Flash"))
)][
order(route_idx, user_idx),
.(route_idx, route_id, user_idx, climb_type, ascent_date, route_rating, bouldering_grade, label)
]Time difference of 11.42 secsthreshold_ascents_dbp
threshold_ascents_dbp <- function(current, lim = 20) {
new <- current |> filter(n() >= lim, .by = user_id) |> filter(n() >= lim, .by = route_id) |> collect()
if (pull(count(current), n) != nrow(new)) {
threshold_ascents_dbp(new, lim)
} else {
duckdb_register(con, "ascent_temp", new)
return(tbl(con, "ascent_temp"))
}
}tbl(con, "climbs_first") |>
threshold_ascents_dbp() |>
## Replacing all route_ratings for a given route_id by its mode
left_join(
tbl(con, "climbs_first") |>
filter(usa_routes %in% route_ratings) |>
count(route_id, usa_routes) |>
slice_max(n, by = route_id) |>
summarize(route_rating = min(usa_routes), .by = route_id),
by = "route_id"
) |>
## Replacing all bouldering_grades for a given route_id by its mode
left_join(
tbl(con, "climbs_first") |>
filter(usa_boulders %in% bouldering_grades) |>
count(route_id, usa_boulders) |>
slice_max(n, by = route_id) |>
summarize(bouldering_grade = min(usa_boulders), .by = route_id),
by = "route_id"
) |>
mutate(
route_idx = dense_rank(route_id),
user_idx = dense_rank(user_id),
label = as.integer(method_name %in% c("Onsight", "Flash"))
) |>
select(route_idx, route_id, user_idx, climb_type, ascent_date, route_rating, bouldering_grade, label) |>
arrange(route_idx, user_idx) |>
compute("climbs_clean", overwrite = TRUE)Time difference of 1.796 secsdata.table [232,887 x 8]
| [ omitted 232,872 entries ] |
2 Model
2.1 Stan code
Updated Stan code using within-chain parallelization
climbing_model
functions {
array[] int sequence(int start, int end) {
array[end - start + 1] int seq;
for (n in 1 : num_elements(seq)) {
seq[n] = n + start - 1;
}
return seq;
}
// Compute partial sums of the log-likelihood
real partial_log_lik_lpmf(array[] int seq, int start, int end,
data array[] int labels, real mean_ability,
data array[] int users, vector user_ability,
data array[] int routes, vector route_difficulty) {
real ptarget = 0;
int N = end - start + 1;
vector[N] mu = mean_ability + rep_vector(0.0, N);
for (n in 1 : N) {
int nn = n + start - 1;
mu[n] += user_ability[users[nn]] - route_difficulty[routes[nn]];
}
ptarget += bernoulli_logit_lpmf(labels[start : end] | mu);
return ptarget;
}
}
data {
int<lower=1> num_ascents;
int<lower=1> num_users;
int<lower=1> num_routes;
array[num_ascents] int<lower=1, upper=num_users> users;
array[num_ascents] int<lower=1, upper=num_routes> routes;
array[num_ascents] int<lower=0, upper=1> labels;
int grainsize;
}
transformed data {
array[num_ascents] int seq = sequence(1, num_ascents);
}
parameters {
real mean_ability;
vector[num_users] user_ability;
vector[num_routes] route_difficulty;
}
model {
user_ability ~ std_normal();
route_difficulty ~ std_normal();
mean_ability ~ std_normal();
target += reduce_sum(
partial_log_lik_lpmf, seq, grainsize,
labels, mean_ability, users, user_ability, routes, route_difficulty
);
}2.2 Stan data
climbing_stan_data <- list(
num_ascents = nrow(climbs_clean),
num_users = n_distinct(climbs_clean$user_id),
num_routes = n_distinct(climbs_clean$route_id),
routes = pull(climbs_clean, route_idx),
users = pull(climbs_clean, user_idx),
labels = pull(climbs_clean, label) |> as.integer(),
grainsize = max(100, round(nrow(climbs_clean) / 50))
)2.3 Model fit
climbing_mod_fit <- climbing_model$sample(
data = climbing_stan_data,
seed = 666,
iter_warmup = 500,
iter_sampling = 1000,
refresh = 0,
chains = 6,
parallel_chains = 6,
threads_per_chain = 5
)
Note
Sampling takes ~4.81 minutes on my CPU (Ryzen 5950X, 16 Cores/32 Threads), on WSL2 (Ubuntu 22)
data.table [6 x 2]
3 Model diagnostics
Plotting random subsets of the traces:
hist_trace_plot
hist_trace_plot <- function(mod, vars) {
draws <- mod$draws(variables = vars, format = "draws_list")
wrap_plots(
mcmc_hist(draws, facet_args = list(nrow = length(vars))),
mcmc_trace(draws, facet_args = list(nrow = length(vars))),
widths = c(1, 1.5)
)
}
Everything seems good.
4 Posterior Predictions
4.1 Posterior data
Getting our Posterior Predictions (subset of 500 draws per route) into long format:
## For each player, take a subsample of 500 draws from their posterior
draws <- climbing_mod_fit$draws(variables = "route_difficulty") |>
bind(x, subset_draws(x, "route_difficulty", regex = TRUE, draw = sample.int(ndraws(x), size = 500))) |>
as.data.table() |>
_[, .(route_difficulty = list(value)), by = variable][, let(
route_idx = as.integer(str_extract(variable, "\\d{1,4}")),
variable = NULL
)]
climbs_pp <- climbs_clean[, .(route_idx, route_id, bouldering_grade, route_rating, climb_type)] |>
unique(by = "route_idx") |>
_[draws, on = "route_idx", nomatch = NULL][order(route_idx)]data.table [4,288 x 6]
| [ omitted 4,273 entries ] |
Time difference of 1.996 secs
## Getting the draws into DuckDB
climbing_mod_fit$draws(variables = "route_difficulty") |>
bind(x, subset_draws(x, "route_difficulty", regex = TRUE, draw = sample.int(ndraws(x), size = 500))) |>
as_draws_df() |>
pivot_longer(everything(), names_to = "route_idx", names_pattern = ".*\\[(\\d{1,4})\\]") |>
duckdb_register(con, "draws", df = _)
## Generating out Posterior Predictions data
tbl(con, "climbs_clean") |>
select(route_idx, route_id, bouldering_grade, route_rating, climb_type) |>
distinct(route_idx, .keep_all = TRUE) |>
inner_join(tbl(con, "draws"), by = "route_idx") |>
summarize(
.by = route_idx,
across(c(bouldering_grade, route_rating, climb_type), first),
route_difficulty = list(value)
) |>
arrange(route_idx) |>
compute("climbs_pp", overwrite = TRUE)data.frame [4,288 x 5]
| [ omitted 4,273 entries ] |
Time difference of 0.7201 secs
With dplyr, we can use the rvar format to encapsulate the samples from the model, which drastically reduces the size of the samples’ data.frame
inner_join(
as.data.frame(climbs_clean) |>
select(route_idx, route_id, bouldering_grade, route_rating, climb_type) |>
distinct(route_idx, .keep_all = TRUE),
tidybayes::spread_rvars(climbing_mod_fit, route_difficulty[route_idx], ndraws = 500),
by = "route_idx"
) |>
arrange(route_idx)data.frame [4,288 x 6]
| [ omitted 4,273 entries ] |
Time difference of 1.039 secs
4.2 Posterior plots
4.2.1 Ridgeline plots
ridgeline_plot
ridgeline_plot <- function(dat, var, title) {
## Unlisting the route_difficulties and making sure the route_ratings/bouldering_grades are ordered properly
dat <- dat[,
.(route_difficulty = unlist(route_difficulty)),
by = setdiff(names(dat), 'route_difficulty')
][, let(
bouldering_grade = factor(bouldering_grade, levels = bouldering_grades),
route_rating = factor(route_rating, levels = route_ratings)
)]
return(
ggplot(dat, aes(route_difficulty, y = {{ var }}, fill = {{ var }})) +
geom_ribbon(
aes(
fill = stage({{ var }}, after_scale = alpha(fill, 0.5)),
ymin = after_stat(group),
ymax = after_stat(group + ndensity * 1.6)
),
stat = "density",
outline.type = "upper",
color = "grey30"
) *
ggblend::blend("multiply") +
geom_vline(xintercept = 0, linetype = "dashed", color = "grey50") +
labs(title = title, x = "Route Difficulty", y = "") +
scale_y_discrete(position = "right") +
theme(
legend.position = "none",
axis.line.y = element_blank(),
plot.title = element_text(hjust = 0.5)
)
)
}Route Ratings:
climbs_pp[climb_type == 0] |>
ridgeline_plot(route_rating, "Climbing Route Posteriors")
Bouldering Grades:
climbs_pp[climb_type == 1 & bouldering_grade != "V0"] |>
ridgeline_plot(bouldering_grade, "Bouldering Problem Posteriors")
4.2.2 Strip plots
strip_plot
strip_plot <- function(dat, var, title) {
strip_plot <- (dat |>
separate_wider_delim(route_id, names = c("crag", "route_name", NA), delim = "__") |>
mutate(
route_difficulty = map_dbl(route_difficulty, mean),
bouldering_grade = factor(bouldering_grade, levels = bouldering_grades),
route_rating = factor(route_rating, levels = route_ratings),
across(c(crag, route_name), \(x) str_replace_all(x, "_", " "))
) |>
ggplot(aes(route_difficulty, y = {{ var }}, color = {{ var }})) +
geom_point(
aes(group = crag, linesize = route_name), # Adding unused aesthetics to get plotly's automated tooltips
position = position_jitter(height = 0.2),
alpha = 0.6
) +
labs(
title = title,
x = "Route Difficulty",
y = str_to_title(str_replace_all(deparse(substitute(var)), "_", " "))
) +
theme(legend.position = "none"))
return(ggplotly(strip_plot, tooltip = c("group", "linesize")))
}Route Ratings:
climbs_pp[climb_type == 0] |>
strip_plot(route_rating, "Climbing Route Difficulties")Bouldering Grades:
climbs_pp[climb_type == 1 & bouldering_grade != "V0"] |>
strip_plot(bouldering_grade, "Bouldering Problem Difficulties")
Tip💻 Expand for Session Info
─ Session info ───────────────────────────────────────────────────────────────
setting value
version R version 4.6.1 (2026-06-24)
os Ubuntu 22.04.5 LTS
system x86_64, linux-gnu
ui X11
language (EN)
collate C.UTF-8
ctype C.UTF-8
tz Europe/Berlin
date 2026-07-13
pandoc 3.8.2.1
quarto 1.10.14 @ /usr/local/bin/quarto
Quarto 1.10.14
Stan (CmdStan) 2.39.0
─ Packages ───────────────────────────────────────────────────────────────────
package * version date (UTC) lib source
archive * 1.1.13 2026-04-12 [1] RSPM
bayesplot * 1.15.0 2025-12-12 [1] RSPM
cmdstanr * 0.9.0.9001 2026-07-11 [1] Github (stan-dev/cmdstanr@2763c0f)
crayon * 1.5.3 2024-06-20 [1] RSPM
data.table * 1.18.99 2026-07-11 [1] Github (Rdatatable/data.table@806ef72)
DBI * 1.3.0 2026-02-25 [1] RSPM
dbplyr * 2.6.0.9000 2026-07-11 [1] Github (Tidyverse/dbplyr@f478e20)
dplyr * 1.2.1 2026-04-03 [1] RSPM
duckdb * 1.5.4.3 2026-07-10 [1] RSPM
fs * 2.1.0 2026-04-18 [1] RSPM
ggplot2 * 4.0.3 2026-04-22 [1] RSPM
gt * 1.3.0 2026-01-22 [1] RSPM
here * 1.0.2 2025-09-15 [1] RSPM
htmltools * 0.5.9 2025-12-04 [1] RSPM
knitr * 1.51 2025-12-20 [1] RSPM
lubridate * 1.9.5 2026-02-04 [1] RSPM
patchwork * 1.3.2 2025-08-25 [1] RSPM
pipebind * 0.1.2 2023-08-30 [1] RSPM
plotly * 4.12.0 2026-01-24 [1] RSPM
posterior * 1.7.0 2026-04-01 [1] RSPM
purrr * 1.2.2 2026-04-10 [1] RSPM
reactable * 0.4.5 2025-12-01 [1] RSPM
stringr * 1.6.0 2025-11-04 [1] RSPM
tidyr * 1.3.2 2025-12-19 [1] RSPM
[1] /home/mar/dev/projects/quarto/ma-riviere.com/renv/profiles/dev-4.6/renv/library/linux-ubuntu-jammy/R-4.6/x86_64-pc-linux-gnu
[2] /opt/R/4.6.1/lib/R/library
* ── Packages attached to the search path.
──────────────────────────────────────────────────────────────────────────────Citation
BibTeX citation:
@online{rivière2022,
author = {Rivière, Marc-Aurèle},
title = {Bayesian {Rock} {Climbing} {Rankings}},
date = {2022-04-19},
url = {https://ma-riviere.com/content/blog/posts/climbing/},
langid = {en},
abstract = {This post is a transposition to R of Ethan Rosenthal’s
{[}blog
post{]}(https://www.ethanrosenthal.com/2022/04/15/bayesian-rock-climbing/)
on modeling Rock Climbing route difficulty using a Bayesian IRT
(Item Response Theory) model. The original Stan code was updated to
use {[}within-chain
parallelization{]}(https://mc-stan.org/docs/2\_30/stan-users-guide/reduce-sum.html)
and {[}compiler
optimization{]}(https://mc-stan.org/docs/2\_30/stan-users-guide/optimization.html)
for faster CPU sampling. Several data processing solutions are
showcased, using either `data.table` or `dbplyr` (with a `DuckDB`
backend), with timings to compare their speed.}
}
For attribution, please cite this work as:
Rivière, M.-A. (2022, April 19). Bayesian Rock Climbing
Rankings. https://ma-riviere.com/content/blog/posts/climbing/