В рекомендациях всплыло.

Если вы все еще думаете, что лекции по статистике – это скучно, то обязательно посмотрите эти выступления Kristin Lennox, она просто потрясающий рассказчик

Everything wrong with statistics (and how to fix it)
All About that Bayes: Probability, Statistics, and the Quest to Quantify Uncertainty

Качество не очень, но зато очень "залипательно"
https://www.youtube.com/watch?v=z23-awbMnGs

О, на блюскай прикольный проект затеяли

WeRateDAGs

Эту тему мы уже поднимали )

WeRateDAGs запостили

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

Там их 2, наверное, нужно последовательно читать, но я про второй как раз
https://www.linkedin.com/pulse/being-wrong-efficacy-probabilities-stephen-senn/?trackingId=p4y7N%2FJHQQyfLZM7YO9O9Q%3D%3D
https://www.linkedin.com/pulse/its-harder-than-you-think-stephen-senn/?trackingId=p4y7N%2FJHQQyfLZM7YO9O9Q%3D%3D

По следам сегодняшнего обсуждения в чате medstatistic решил сделать что-нибудь полезное.
Это будет байесовская лапласовская версия t-теста Уэлча*.

Почему я решил показать именно это (в смысле t-тест и именно такую его реализацию)? Просто потому что она мне самому кажется эстетически привлекательной. Представьте, люди очень давно могли посчитать на бумажке параметры апостериорного распределения (при строго определенных семействах априорных распределений, которые называются сопряженными с нормальным) для параметров нормального распределения матожидания (μ) и точности (1/σ²), но проблема в том, что это возможно сделать только если известны σ и μ, соответственно, то есть можно сказать это был абсолютно бесполезный навык с практической точки зрения.
Но вот в конце XX века появилась идея сэмплинга Гиббса – такого способа получения выборки из совместного распределения параметров, при котором мы на самом деле производим сэмплинг только из условных распределений. В данном случае, мы можем воспользоваться нашим древним знанием о сопряженных моделях для μ и 1/σ² и сделать очень простой, но при этом крайне эффективный сэмплер (демонстрация того, что это именно так приведена ниже).

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

Итак, сэмплер Гиббса на основе сопряженных распределений:

normal_gibbs <- function(y, mu_0, tau_0, a, b, n_iter = 10^5) {
    n_obs <- length(y)
    sum_y <- sum(y)
    mu_sample <- tau_sample <- numeric(n_iter)
    
    # assigning starting values
    mu <- mean(y)
    tau <- 1 / var(y)
    
    # Gibbs sampling
    for (i in seq_len(n_iter)) {
        mu <- mu_sample[[i]] <- rnorm(
            n = 1,
            # parameters of conjugate posterior for μ given σ²
            mean = (sum_y * tau + mu_0 * tau_0) / (n_obs * tau + tau_0),
            sd = 1 / sqrt(n_obs * tau + tau_0)
        )
        tau <- tau_sample[[i]] <- rgamma(
            n = 1,
            # parameters of conjugate posterior for precision (1/σ²) given μ 
            shape = a + n_obs / 2,
            rate = b + 1 / 2 * sum((y - mu)^2)
        )
    }
    
    list("mu" = mu_sample, "sd" = sqrt(1/tau_sample))
}

Данные для примера:

set.seed(13799)

group_1 <- rnorm(10, mean = 120, sd = 15)
group_2 <- rnorm(15, mean = 130, sd = 25)

res_1 <- normal_gibbs(
    group_1,
    # parameters for prior for mu
    mu_0 = 100, tau_0 = 1 / (100^2),
    # parameters for prior for tau
    a = 0.01, b = 0.01
)

res_2 <- normal_gibbs(
    group_2,
    # parameters for prior for mu
    mu_0 = 100, tau_0 = 1 / (100^2),
    # parameters for prior for tau
    a = 0.01, b = 0.01
)

Апостериорное распределение для разности матожиданий, то есть собственно "t-тест"

Апостериорное распределение отношения стандартных отклонений, то чего в t-тесте мы точно не получим

Если кому-то интересно погрузиться поглубже и понять почему в вызовах rnorm() и rgamma() используются именно эти странно выглядящие параметры, то можно почитать этот материал Герберта Ли

Насколько сэмплинг Гиббса, как и MCMC в целом, могут быть неэффективными? Приведу недавний пример из другого обсуждения. Отмечу только, что здесь, конечно, проблема не только в том, что мы не знаем сопряженных априорных распределений для параметров модели, но и

1) Логистическая регрессия сама по себе несопоставимо более сложная штука по сравнению с линейной регрессией (и t-тестом), особенно в байесовском лапласовском сеттинге.
2) Набор данных – очень сложный, с (почти) полным разделением.

df <- brglm2::endometrial

glm(
    NV ~ EH + HG,
    data = df,
    family = "binomial"
) |>
    broom::tidy()

#> # A tibble: 3 × 5
#>   term        estimate std.error statistic p.value
#>   <chr>          <dbl>     <dbl>     <dbl>   <dbl>
#> 1 (Intercept)   -16.5    2355.    -0.00699  0.994 
#> 2 EH             -2.36      1.28  -1.84     0.0656
#> 3 HG             18.8    2355.     0.00799  0.994

3) Современные люди использовали бы Гамильтоновское Монте-Карло (HMC)
4) Настоящие байесианцы, коим я не являюсь, знают много хакерских приемов и могут сделать на порядки лучше с помощью тонкой настройки параметров сэмплинга и репараметризации модели.

Пример:

df <- brglm2::endometrial
set.seed(525243)

fit <- MCMCglmm::MCMCglmm(
    NV ~ EH + HG,
    family = "categorical",
    data = df,
    prior = list(
        B = list(
            mu = rep(0, 3),
            V = MCMCglmm::gelman.prior(
                ~ EH + HG, 
                data = df, 
                scale = sqrt(pi^2/3+1)
            )
        ),
        R = list(V = 1, fix = 1)
    ),
    nitt = 100000,
    burnin = 1000,
    thin = 10,
    verbose = FALSE
)

summary(fit)
#> 
#>  Iterations = 1001:99991
#>  Thinning interval  = 10
#>  Sample size  = 9900 
#> 
#>  DIC: 42.77461 
#> 
#>  R-structure:  ~units
#> 
#>       post.mean l-95% CI u-95% CI eff.samp
#> units         1        1        1        0
#> 
#>  Location effects: NV ~ EH + HG 
#> 
#>             post.mean l-95% CI u-95% CI eff.samp   pMCMC   
#> (Intercept)   -1.3239  -4.6634   1.8131   1430.6 0.43071   
#> EH            -1.8184  -3.4936  -0.1290    987.7 0.03131 * 
#> HG             2.9383   0.8804   5.1687    733.0 0.00182 **
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

На что здесь следует обратить внимание: в примере с t-тестом мы сделали 2 раза сэмплинг по 10⁵ и из-за отсутствия автокорреляции между элементами можем считать, что эффективный размер выборки (ESS) примерно равен физическому размеру (то есть 100% сэмплов несут независимую информацию об апостериорном распределении).
В последнем примере мы взяли выборку размера 10⁵, при этом отбирали с лагом в 10 (то есть каждый 10-ый сэмпл из цепи попал в выборку), в итоге суммарный размер составил 10⁶. Но что мы видим в саммари: ESS для угловых параметров измеряется в сотнях.То есть эффективность сэмплинга составила примерно 10³ / 10⁶ (0.1% сэмплов несут независимую информацию об апостериорном распределении).
Я когда это увидел был в полном шоке.

Ну и тематическая шутка, как водится

В комментариях есть еще мысли про 5 этап

Это оказалось очень тонким центральным моментом методологии, не все, что похоже на сэмплинг Гиббса им является, я раньше, например, считал (и на занятиях говорил), что MICE - это разновидность сэмплинга Гиббса, оказалось, что это не так, хотя и похоже
https://tpmorris.substack.com/p/mice-is-not-a-gibbs-sampler