В среду снова радикально меняем тему; на этот раз у нас практически классическое машинное обучение, да ещё и с байесовским уклоном. Уверен, что в исполнении Максима будет интересно.

BART
Ссылка на трансляцию (среда 19 ноября, 14:00)

Максим Николаев
(МКН СПбГУ)

BART, Bayesian Additive Regression Trees — это непараметрическая модель, которая наследует выразительность от ансамблей деревьев решений и описание неопределённости от байесовского подхода. Как это часто бывает, апостериорное распределение этой модели краткой аналитической формы не имеет, поэтому для работы с ней используются методы Монте-Карло на марковских цепях (MCMC).

Мы обсудим устройство модели, а также используемые методы MCMC. В оставшееся время обсудим направления для дальнейших исследований. Материал не требует специфической подготовки, но будет полезным понимание основных понятий теории цепей Маркова вплоть до стационарного распределения.

Typst - классный
https://www.youtube.com/watch?v=jY2nCypm0Xs

Очень прикольно сделанный сайт, позволяющий критически переосмыслить такие незыблемые вещи как ROC и AUC:

https://predictionperformancediscrimination.netlify.app/#/discrimination

Это выглядит слишком хорошо, чтобы быть правдой.
Не понимаю, почему метод так не популярен в биомедицине.

suppressPackageStartupMessages({
    library(tidyverse)
    library(logistf)
})

set.seed(15321)

n <- 150
p1 <- 0.015
p2 <- 0.03

x <- rep(0:1, each = n)
probs <- rep(c(p1, p2), each = n)

sim <- function() {
    df <- data.frame(
        x = x,
        y = rbinom(
            n = n * 2, size = 1, 
            prob = probs
        )
    )
    
    fit <- logistf(
        y ~ x, data = df, 
        pl = TRUE, firth = TRUE
    )
    
    c(
        n1 = sum(df$y[x == 0]), 
        n2 = sum(df$y[x == 1]),
        beta = fit$coefficients[[2]],
        lcl = fit$ci.lower[[2]], 
        ucl = fit$ci.upper[[2]],
        lrt = 1 - pchisq(-2 * diff(fit$loglik), df = 1)
    )
}

res <- replicate(10000, sim())

eff <- log((p2/(1 - p2))/(p1/(1 - p1)))

t(res) |>
    as_tibble() |>
    mutate(eff = eff) |>
    summarise(
        coverage = mean(between(eff, lcl, ucl)),
        power = mean(lrt.null < 0.05),
        `S-type` = mean(beta < 0 & lrt.null < 0.05),
        mean = mean(beta),
        true = mean(eff),
        bias = mean - true
    )

Здесь очень сложная ситуация.
Редкие события и маленькая выборка, с учетом редкости событий в одной из групп с 10% вероятностью событий не будет вовсе.
Но я ничего не фильтровал и никаких na.rm = TRUE, все 10000 моделей сошлись и ДИ почти обеспечивают заявленную альфу (полуторапроцентное перепокрытие мне представляется даже чем-то хорошим с практической точки зрения).
Просто попробуйте, кому интересно, glm() с вальдовскими интервалами при такой ситуации запустить (вернее даже при более мягкой, здесь с glm в общем-то ловить совсем нечего).

В рассматриваемом сценарии хотя бы в одной из групп не было событий в 1105 случаях, в обеих (!) группах не было событий в 10 случаях (да, в этих случаях мы тоже получили оценки эффекта, понятно, что это не очень информативно, но все же прикольно).

Интересно будет еще байесовский бутстреп с этой моделью попробовать релизовать для RD и RR

Вдруг кому-то пригодится

Не ответ на вопрос, но в последнее время я обратил внимание, что очень часто в качестве учебника по статистике рекомендуют вот эту книгу

https://www.routledge.com/Statistical-Inference/Casella-Berger/p/book/9781032593036

Вот, буквально на днях

Это очень интересные сюжеты, кстати

В курсе "Основы байесовского вывода" сегодня поговорили о двух важных общих сюжетах:

СПбГУ — 2025.11.27 — Принцип максимума энтропии и априорные распределения Джеффриса
(слайды и доска на странице курса)

Принцип максимума энтропии показывает, какие распределения вероятностей являются "наиболее характерными". Сам принцип приходит из статистической физики, из работ Гиббса и того же Эдвина Джейнса, но и в машинном обучении тоже встречается. Так что в виде небольшого лирического отступления рассказал, откуда берётся максимизация энтропии, и привёл пару примеров.

А главный сегодняшний объект, априорные распределения Джеффриса, решает проблему, которая возникает из наивного вопроса: что происходит с априорными распределениями при репараметризации? Предположим, что мы хотим выразить незнание о параметре монетки, и выражаем его равномерным априорным распределением на [0, 1]. Но если мы перейдём от вероятности орла, скажем, к log-odds, что является даже более естественной параметризацией, равномерное распределение превратится в бог знает что... Сэр Рональд Фишер критиковал за это весь байесианизм в целом, а другой сэр Гарольд Джеффрис предложил потребовать инвариантность при репараметризации как свойство априорных распределений, и получились как раз в довольно глубоком смысле "распределения полного незнания".

#spsu #lectures #bayes2025

И вот что еще удивительно: априорные распределения Джеффриса имеют непосредственное отношение к фриквентистскому методу, который недавно рассматривали.

Точечные оценки ОШ получаются одинаковыми при использовании поправки Холдейна, регрессии Фирта и апостериорного ОШ при использовании априорного распределения Джеффриса [θ ~ Beta(1/2, 1/2)] в бета-биномиальной модели.

suppressPackageStartupMessages({
    library(tidyverse)
})

df <- tibble(
    group = 0:1,
    N = c(25, 25),
    events = c(0, 3)
)

# Haldane correction for OR
df |>
    mutate(
        odds = (events + 0.5) / (N - events + 0.5)
    ) |>
    pull(odds) |>
    (\(x) x[2] / x[1])()
#> [1] 7.933333

# Firth's logistic regression
df |>
    mutate(N = N - events) |>
    pivot_longer(
        -group, 
        names_to = "outcome",
        values_to = "N"
    ) |>
    mutate(
        outcome = as.integer(outcome == "events")
    ) |>
    logistf::logistf(
        outcome ~ group, 
        data = _,
        weights = N,
        firth = TRUE
    ) |>
    coef() |>
    _[[2]] |>
    exp()
#> [1] 7.933333

# Jeffreys prior for beta-binomial model
df |>
    mutate(
        # posterior parameters
        alpha = events + 0.5,
        beta = N - events + 0.5,
        # posterior mean
        mean = alpha / (alpha + beta),
        # posterior odds
        odds = mean / (1 - mean)
    ) |>
    pull(odds) |>
    (\(x) x[2] / x[1])()
#> [1] 7.933333

Это, кстати, основа эпистемологического обоснования использования нормального распределения (точнее нормального правдоподобия), т.е. без ЦПТ. Очень прикольная тема, сейчас промотал лекцию, Николенко в ней как раз это обоснование в качестве примера приводит.

Сейчас подумал, что MCMC (вернее анализ того, что происходит внутри цепи) очень ярко отражают контринтуитивность теорвера и матстата.

Пост Вехтари в блоге Гелмана про effective sample size

То есть, сложно же себе представить такое: есть случайный процесс, в котором каждое последующее наблюдение отрицательно коррелирует с предыдущим и в итоге вы получаете больше информации о маргинальном распределении по сравнению с информацией, что несет выборка с полностью незваисимыми наблюдениями, потом просто трансформируете элементы выборки и информации о маргинальном распределении трансформированой величины становится меньше.

Пришла в голову следующая задачка:
Думаю, что многие знают, что t-тест Стьюдента является просто частным случаем МНК регрессии.

set.seed(42)
x <- rep(0:1, times = c(10, 15))
y <- rnorm(length(x), mean = x)

t.test(y ~ x, var.equal = TRUE) |>
    broom::tidy() |>
    dplyr::select(
        estimate, statistic, p.value
    )
#> # A tibble: 1 × 3
#>   estimate statistic p.value
#>      <dbl>     <dbl>   <dbl>
#> 1   -0.400    -0.755   0.458

lm(y ~ x) |>
    broom::tidy() |>
    dplyr::filter(term != "(Intercept)") |>
    dplyr::select(
        estimate, statistic, p.value
    )
#> # A tibble: 1 × 3
#>   estimate statistic p.value
#>      <dbl>     <dbl>   <dbl>
#> 1    0.400     0.755   0.458

Вопрос состоит в следующем: как можно с помощью регрессионной модели, которая потенциально может включать не только один предиктор, воспроизвести результаты t-теста Уэлча, который не делает допущения о гомоскедастичности? 🙃

t.test(y ~ x, var.equal = FALSE) |>
    broom::tidy() |>
    dplyr::select(
        estimate, statistic, p.value
    )
#> # A tibble: 1 × 3
#>   estimate statistic p.value
#>      <dbl>     <dbl>   <dbl>
#> 1   -0.400    -0.845   0.407

Я нашел такое.
В нем для воспроизведения важны не только структура дисперсии, но и использование REML и аппроксимация df по Саттервайту.

set.seed(42)
x <- rep(0:1, times = c(10, 15))
y <- rnorm(length(x), mean = x)

t.test(y ~ x, var.equal = FALSE) |>
    broom::tidy() |>
    dplyr::select(
        estimate, statistic, df = parameter, p.value
    )
#> # A tibble: 1 × 4
#>   estimate statistic    df p.value
#>      <dbl>     <dbl> <dbl>   <dbl>
#> 1   -0.400    -0.845  22.4   0.407

nlme::gls(
    y ~ x, 
    weights = nlme::varIdent(form = ~ 1 | x),
    method = "REML"
) |>
    emmeans::emmeans(~ x, mode = "satterthwaite") |>
    emmeans::contrast("revpairwise") |>
    tibble::as_tibble() |>
    dplyr::select(
        estimate, 
        statistic = t.ratio,
        df, p.value
    )
#> # A tibble: 1 × 4
#>   estimate statistic    df p.value
#>      <dbl>     <dbl> <dbl>   <dbl>
#> 1    0.400     0.845  22.4   0.407

Закину сюда, довольно интересные материалы