Блин, это гениально
https://www.youtube.com/watch?v=uqwC41RDPyg

Блин, вспомнил с этим Харди-Вайнбергом про свою медицинскую юность. Ностальгия прям. Полез посмотреть, что там сейчас с проектами, где я принимал участие.

Вот наткнулся, аж слезы на глазах )

2015 год, госпиталь Мулаго, Кампала, Уганда
Во втором ряду крайний слева – это я если что )

Здесь я этого почему-то не постил еще, исправляю.
Из курса МакЭлрита

И еще из закромов. Рубрика "Andrew Gelman says".

"...like Sander Greenland who disagrees with everybody"

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

Если вы все еще думаете, что лекции по статистике – это скучно, то обязательно посмотрите эти выступления 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) Настоящие байесианцы, коим я не являюсь, знают много хакерских приемов и могут сделать на порядки лучше с помощью тонкой настройки параметров сэмплинга и репараметризации модели.