8  Байесовский регрессионный анализ

library(tidyverse)

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

8.1 Основы регрессионного анализа

Когда мы используем регрессионный анализ, мы пытаемся оценить два параметра:

• свободный член (intercept) — значение \(y\) при \(x = 0\);
• угловой коэффициент (slope) — изменение \(y\) при изменении \(x\) на одну единицу.

\[y_i = \hat{\beta_0} + \hat{\beta_1}\times x_i + \epsilon_i\]

Причем, иногда мы можем один или другой параметр считать равным нулю.

При этом, вне зависимости от статистической школы, у регрессии есть свои ограничения на применение:

• линейность связи между \(x\) и \(y\);
• нормальность распределение остатков \(\epsilon_i\);
• гомоскидастичность — равномерность распределения остатков на всем протяжении \(x\);
• независимость переменных;
• независимость наблюдений друг от друга.

Припомнив предыдущий раздел, мы можем представить связь и в другой нотации:

\[y_{n} \sim Normal(\mu = \hat{\beta_0} + \hat{\beta_1}\times x_{n}, \sigma = \epsilon)\]

8.2 Регрессионная модель в brms

В работе (Stepanova 2011) сообщается, что средняя скорость речи на русском языке составляет 5.31 слогов в секунду (стандартное отклонение 1.93). Скачайте данные корпуса г. Звенигород (https://lingconlab.ru/zvenigorod/) и оцените какой эффект на скорость речи носителя ea1976 оказывает разное количество гласных в высказывании, используя данные из (Stepanova 2011) в качестве априорного распределения. Визуализируйте получившееся распределение.

library(brms)
n_cores <- 15 # parallel::detectCores() - 1

read_csv("https://raw.githubusercontent.com/agricolamz/2026_HSE_b_da4l/master/data/speech_tempo_zvenigorod.csv") |> 
  filter(speaker == "ea1976") ->
  df

get_prior(speech_tempo ~ n_vowels, 
          family = gaussian(),
          data = df) 

Исходя из формулы выше, нам нужно определить априорное распределение для

• свободного члена;
• углового коэффециента;
• распределения остатков.

fit1 <- brm(speech_tempo ~ n_vowels,
            data = df,
            family = gaussian(),
            cores = n_cores,
            silent = TRUE,
            prior = c(prior(normal(5.31, 1.93), class = "Intercept", lb = 0),
                      prior(normal(0, 1), class = "sigma", lb = 0),
                      prior(normal(0, 10), class = "b", lb = 0)))

После того, как функция отработала, мы можем проанализировать то, что мы получили в результате ее работы. Мы оценили три параметра:

plot(fit1)

Стоит обратить внимание на сходимость цепей. Если все цепи накладываются друг на друга и выглядят как волосатые гусеницы — все в порядке. Вот пример плохой сходимости цепей:

Видно, что цепи исследуют разное пространство и “устаканились” каждая на своем уровне. Кроме того, может быть такое, что цепи в принципе на одном уровне, но где-нибудь по середине делает скачок в сторону, а потом возвращается.

Метрики модели можно посмотреть, просто вызвав переменную с моделью.

fit1

 Family: gaussian 
  Links: mu = identity 
Formula: speech_tempo ~ n_vowels 
   Data: df (Number of observations: 786) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Regression Coefficients:
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept     2.60      0.08     2.45     2.74 1.00     4091     3051
n_vowels      0.11      0.01     0.10     0.12 1.00     3989     2794

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     1.09      0.03     1.04     1.15 1.00     3445     2754

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

Метрики качества модели: R-hat, Bulk_ESS и Tail_ESS

R-hat – это диагностика сходимости, которая сравнивает оценку модели, которая получается внутри цепи и между цепями (введена в (Gelman and Rubin 1992)). Если цепи не перемешались (например, если нет согласия внутри цепи и между цепями) R-hat будет иметь значение больше 1. Собственно для этого рекоммендуется запускать по-крайней мере 4 цепи. Stan возвращает что-то более сложное, что называется maximum of rank normalized split-R-hat и rank normalized folded-split-R-hat, однако мы не будем в этом разбиратсья, отсылаю к статье или к посту в блоге одного из авторов с объяснением.

Самое главное: \(\hat{R} \leq 1.05\) – значит проблем со сходмиостью цепей не обнаружено.

Bulk ESS (Effective sample size) — это оценка качества сэмплирования в середине апостериорного распределения, а Tail ESS — в краях апостериорного распределения. Чем выше ESS – тем лучше.

Мы можем визуализировать полученную регрессию:

fit1 |> 
  conditional_effects() |> 
  plot(points=TRUE)

Мы можем визуализировать, насколько скорость речи в данных соотносится с полученным распределением апостериорного распределения:

fit1 |> 
  pp_check(type = "dens_overlay", ndraws = 50)

fit1 |> 
  pp_check(type = "stat")

Если хочется предсказать какие-то значения, то можно использовать функцию predict(), которая возьмет выборки из апостериорного распределения:

predict(fit1, newdata = tibble(n_vowels = c(10, 20)))

     Estimate Est.Error     Q2.5    Q97.5
[1,] 3.758827  1.108490 1.586769 5.894793
[2,] 4.870823  1.093665 2.675327 6.976102

8.3 Альтернативная модель

Мы можем использовать тот же самый код, но попробовать в качестве переменной использовать логорифм количества гласных. Для того, чтобы использовать какие-то арифемтические операции внутри регрессионной формулы, нужно использовать функцию I() или сделать какие-то модификации переменных до того, как использовать в регрессии.

get_prior(speech_tempo ~ I(log(n_vowels)), 
          data = df, 
          family = gaussian())

fit2 <- brm(speech_tempo ~ I(log(n_vowels)), 
            data = df,
            family = gaussian(),
            cores = n_cores,
            silent = TRUE,
            prior = c(prior(normal(5.31, 1.93), class = "Intercept", lb = 0), 
                      prior(normal(0, 1), class = "sigma", lb = 0),
                      prior(normal(0, 10), class = "b", lb = 0)))

plot(fit2)

fit2 |> 
  conditional_effects() |> 
  plot(points=TRUE)

fit2 |> 
  pp_check(type = "dens_overlay", ndraws = 50)

fit2 |> 
  pp_check(type = "stat")

Как интерпретировать полученную регрессию с двумя отлогорифмированными значениями?

В обычной линейной регресии мы узнаем отношения между \(x\) и \(y\): \[y_i = \beta_0+\beta_1\times x_i\]

Как изменится \(y_j\), если мы увеличем \(x_i + 1 = x_j\)? \[y_j = \beta_0+\beta_1\times x_j\]

\[y_j - y_i = \beta_0+\beta_1\times x_j - (\beta_0+\beta_1\times x_i) = \beta_1(x_j - x_i)\]

Т. е. \(y\) увеличится на \(\beta_1\) , если \(x\) увеличится на 1. Что же будет с логарифмированными переменными? Как изменится \(y_j\), если мы увеличем \(x_i + 1 = x_j\)?

\[\log(y_j) - \log(y_i) = \beta_1\times (\log(x_j) - \log(x_i))\]

\[\log\left(\frac{y_j}{y_i}\right) = \beta_1\times \log\left(\frac{x_j}{x_i}\right) = \log\left(\left(\frac{x_j}{x_i}\right) ^ {\beta_1}\right)\]

\[\frac{y_j}{y_i}= \left(\frac{x_j}{x_i}\right) ^ {\beta_1}\]

Т. е. \(y\) увеличится на \(\beta_1\) процентов, если \(x\) увеличится на 1 процент. Эта процентная логика будет работать при любой комбинации. Например, если \(y\) не логорифмирован, а \(x\) логорифмирован интерпретировать это стоит в той же логике: \(y\) увеличится на \(\beta_1\), если \(x\) увеличится на 1 процент.

В работе (Coretta 2016) собраны данные длительности исландских гласных. Используя байесовскую регрессию с априорными распределениями по умолчанию, смоделируйте длительность гласного (vowel.dur) в зависимости от придыхания (aspiration) у носителя tt01. Визуализируйте результаты модели.

8.3.1 Сравнение моделей

В работах (Vehtari et al. 2017, 2024) описан метод loo, который оценивает leave-one-out кроссвалидацию, которую можно использовать для сравнения байесовских моделей.

fit1 <- add_criterion(fit1, "loo")
fit1$criteria$loo

Computed from 4000 by 786 log-likelihood matrix.

         Estimate   SE
elpd_loo  -1187.2 22.0
p_loo         3.2  0.3
looic      2374.5 44.0
------
MCSE of elpd_loo is 0.0.
MCSE and ESS estimates assume MCMC draws (r_eff in [0.8, 1.1]).

All Pareto k estimates are good (k < 0.7).
See help('pareto-k-diagnostic') for details.

fit2 <- add_criterion(fit2, "loo")
fit2$criteria$loo

Computed from 4000 by 786 log-likelihood matrix.

         Estimate   SE
elpd_loo  -1089.3 24.1
p_loo         3.2  0.4
looic      2178.6 48.3
------
MCSE of elpd_loo is 0.0.
MCSE and ESS estimates assume MCMC draws (r_eff in [0.9, 1.2]).

All Pareto k estimates are good (k < 0.7).
See help('pareto-k-diagnostic') for details.

loo_compare(fit1, fit2, criterion = "loo")

     elpd_diff se_diff
fit2   0.0       0.0  
fit1 -97.9      10.4  

Большее значение ELPD стоит интепретировать как лучшую предсказательную способность модели. Подробнее про эту метрику можно прочитать, запустив help('pareto-k-diagnostic'). Еще больше на странице пакета loo.

8.4 Категориальные переменные

8.4.1 Default contrast coding: Treatment contrasts

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

Как этот график соотносится с оценками параметров модели?

Мы видим, что свободный член (b_Intercept) соотносится с отсутствием придыхания, а переменная b_aspirationyes с разницей средних между случаев с придыханием и без. Дело в том, что когда запускается регрессию с категориальной переменной, данная переменная автоматически разбивается на группу бинарных dummy-переменных:

aspiration	dummy_aspiration_no	dummy_aspiration_yes
без придыхания	1	0
с придыханием	0	1

Дальше для регрессионного анализа выкидывают одну из переменных, так как иначе модель не сойдется (dummy-переменных всегда n-1, где n — количество категорий в переменной).

aspiration	dummy_aspiration_yes
без придыхания	0
с придыханием	1

Если переменная dummy_aspiration_yes принимает значение 1, значит речь о случаях, где согласный придыхательный, а если принимает значение 0, то о случаях, где согласный непридыхательный. Если вставить нашу переменную в регрессионную формулу получится следующее:

\[y_i = \hat\beta_0 + \hat\beta_1 \times \text{dummy-aspiration-yes} + \epsilon_i,\]

Так как dummy_aspiration_yes принимает либо значение 1, либо значение 0, то получается, что модель предсказывает лишь два значения:

\[y_i = \left\{\begin{array}{ll}\hat\beta_0 + \hat\beta_1 \times 1 + \epsilon_i = \hat\beta_0 + \hat\beta_1 + \epsilon_i\text{, если предыдущий согласный придыхательный}\\ \hat\beta_0 + \hat\beta_1 \times 0 + \epsilon_i = \hat\beta_0 + \epsilon_i\text{, если предыдущий согласный непридыхательный} \end{array}\right.\]

Таким образом, получается, что свободный член \(\beta_0\) и угловой коэффициент \(\beta_1\) в регрессии с категориальной переменной получает другую интерпретацию. Одно из значений переменной кодируется при помощи \(\beta_0\), а сумма коэффициентов \(\beta_0+\beta_1\) дают другое значение переменной. Так что \(\beta_1\) — это разница между оценками средних двух значений переменной.

В случаях, когда значений больше — логика будет та же: будет создаваться n-1 dummy переменных, каждый коэффициент при таких переменных будет говорить о разнице средних переменной и некоторого, выделенного в интерсепт значения, которое становится таким эталоном сравнения. Автоматическая процедура отправляет в эталон сравнения уровень, название которого раньше по алфавиту. Если хочется выбрать другой эталон сравнения, то можно использовать функцию factor() и эксплицитно задать уровни или же использовать функцию fct_relevel(your_variable, "make me the baseline") из пакета forcats (часть tidyverse).

8.4.2 Sum contrasts

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

read_csv("https://raw.githubusercontent.com/agricolamz/2024_HSE_b_da4l/master/data/Coretta_2017_icelandic.csv") |> 
  filter(speaker == "tt01") |> 
  mutate(aspiration = factor(aspiration)) ->
  vowels

contrasts(vowels$aspiration)

    yes
no    0
yes   1

contrasts(vowels$aspiration) <- c(-1, 1)

get_prior(vowel.dur ~ aspiration,
          data = vowels)

vowels_fit2 <- brm(vowel.dur ~ aspiration,
                   data = vowels,
                   cores = n_cores,
                   silent = TRUE)

g)
Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 3: 
Chain 3:  Elapsed Time: 0.02 seconds (Warm-up)
Chain 3:                0.01 seconds (Sampling)
Chain 3:                0.03 seconds (Total)
Chain 3: 

plot(vowels_fit2)

vowels_fit2 |> 
  conditional_effects() |> 
  plot()

Мы видим, что графики эффектов идентичны, однако графики коэффециентов разные, потому что эффекты теперь закодированы не друг относительно друга, а как некоторое средние между группами \(\pm\) некоторое значение, отраженное в коэффициенте \(\hat\beta_1\) (в нашем случае, если выразить все при помощи средних — 89.38 \(\pm\) 11.18)

vowels_fit2

 Family: gaussian 
  Links: mu = identity 
Formula: vowel.dur ~ aspiration 
   Data: vowels (Number of observations: 175) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Regression Coefficients:
            Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept      89.36      1.68    85.88    92.64 1.00     4213     2739
aspiration1   -11.18      1.70   -14.57    -7.86 1.00     4143     3200

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma    22.43      1.21    20.22    24.99 1.00     3461     2942

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

8.4.3 One-hot encoding

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

vowels_fit3 <- brm(vowel.dur ~ -1 + aspiration,
                   data = vowels,
                   cores = n_cores,
                   silent = TRUE)

plot(vowels_fit3)

vowels_fit3 |> 
  conditional_effects() |> 
  plot()

vowels_fit3

 Family: gaussian 
  Links: mu = identity 
Formula: vowel.dur ~ -1 + aspiration 
   Data: vowels (Number of observations: 175) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Regression Coefficients:
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
aspirationno    100.64      2.41    95.83   105.36 1.00     3839     2871
aspirationyes    78.30      2.40    73.60    82.97 1.00     3163     2520

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma    22.36      1.19    20.18    24.90 1.00     3937     2839

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

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

8.4.4 Другие контрасты

• Repeated contrasts: контрасты для анализа упорядоченных переменных (например, образование, возраст и т. п.)
• Helmert contrasts: первый контрсат сравнивает между уровнями L1 и L2, следующий сравнивает уровне L3 и средние между (L1 и L2), следующий сравнивает L4 и среднее между (L1, L2 и L3) и т. д.
• и другие

Вообще есть семейство функций, позволяющее задать контрасты:

• contr.treatment()
• contr.sum()
• contr.helmert()
• contr.poly()

8.5 Иерархические модели

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

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

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

• основные эффекты – это те связи, которые нас интересуют, независимые переменные (скорость речи, наличие придыхания…);
• случайные эффекты – это те переменные, которые создают группировку в данных (корпус, носитель, стимул и т. п.).

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

Из ‘Introduction to Bayesian Data Analysis for Cognitive Science’ Bruno Nicenboim, Daniel J. Schad, and Shravan Vasishth

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

df <- read_csv("https://raw.githubusercontent.com/agricolamz/2026_HSE_b_da4l/master/data/speech_tempo_zvenigorod.csv")

Случайные эффекты записываются в круглых скобках. Почему перед ним тоже нужно писать 1 (что означает интерсепт) станет ясно из примеров ниже. Давайте сравним априорные распределения по-умолчанию, которые предлагает brms для моделей со смешанными эффектами и для моделей без них:

get_prior(speech_tempo ~ 1,
          data = df,
          family = gaussian())

get_prior(speech_tempo ~ (1|speaker),
          data = df,
          family = gaussian())

Как видно добавились много параметров sd, которыми принято обозначать параметры вариативности модели внутри групп, в то время как sigma продолжает быть параметром, определяющим вариативность между группами.

fit_hierarchical <- brm(speech_tempo ~ (1|speaker),
                        data = df,
                        family = gaussian(),
                        cores = n_cores,
                        silent = TRUE,
                        prior = c(prior(normal(5.31, 1.93), class = "Intercept", lb = 0),
                                  prior(normal(2, 5), class = "sigma", lb = 0),
                                  prior(normal(2, 10), class = "sd")),
                        iter = 4000,
                        control = list(adapt_delta = .85))

Warning: There were 3 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.

Warning: Examine the pairs() plot to diagnose sampling problems

fit_hierarchical

Warning: There were 3 divergent transitions after warmup. Increasing
adapt_delta above 0.85 may help. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup

 Family: gaussian 
  Links: mu = identity 
Formula: speech_tempo ~ (1 | speaker) 
   Data: df (Number of observations: 5623) 
  Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
         total post-warmup draws = 8000

Multilevel Hyperparameters:
~speaker (Number of levels: 9) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.88      0.28     0.52     1.61 1.00     1297     1837

Regression Coefficients:
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept     4.58      0.30     4.00     5.20 1.00     1223     1435

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     1.38      0.01     1.35     1.40 1.00     3716     3450

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

plot(fit_hierarchical)

В случаях, когда у нас нет предикторов, визуализация параметров иллюстрирует необходимые распределения.

8.6 Другие варианты моделей с группировками

8.6.1 Complete pooling model

Из ‘Introduction to Bayesian Data Analysis for Cognitive Science’ Bruno Nicenboim, Daniel J. Schad, and Shravan Vasishth

fit_complete_pooling <- brm(speech_tempo ~ 1,
                            data = df,
                            family = gaussian(),
                            cores = n_cores,
                            silent = TRUE)

fit_complete_pooling

 Family: gaussian 
  Links: mu = identity 
Formula: speech_tempo ~ 1 
   Data: df (Number of observations: 5623) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Regression Coefficients:
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept     4.64      0.02     4.60     4.68 1.00     2943     2454

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     1.59      0.02     1.56     1.62 1.00     3514     2613

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

plot(fit_complete_pooling)

8.6.2 No pooling model

Из ‘Introduction to Bayesian Data Analysis for Cognitive Science’ Bruno Nicenboim, Daniel J. Schad, and Shravan Vasishth

fit_no_pooling <- brm(speech_tempo ~ 0 + speaker,
                      data = df,
                      family = gaussian(),
                      cores = n_cores,
                      silent = TRUE,
                      prior = c(prior(normal(2, 5), class = "sigma", lb = 0),
                                prior(normal(0, 10), class = "b", lb = 0)))

fit_no_pooling

 Family: gaussian 
  Links: mu = identity 
Formula: speech_tempo ~ 0 + speaker 
   Data: df (Number of observations: 5623) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Regression Coefficients:
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
speakerea1976     4.06      0.05     3.97     4.16 1.00     8765     2506
speakerli1950     4.33      0.07     4.20     4.46 1.00     6895     2851
speakeroa1975     4.73      0.06     4.61     4.85 1.00     8704     3010
speakerom1973     4.23      0.04     4.15     4.32 1.00     7934     2554
speakertk1949     4.91      0.07     4.77     5.04 1.00     8222     2834
speakerva1940     3.50      0.06     3.39     3.61 1.00     7889     2865
speakervg1945     4.74      0.08     4.57     4.90 1.00     8704     3254
speakervi1972     6.14      0.04     6.06     6.22 1.00     7674     2863
speakerya1968     4.46      0.06     4.35     4.57 1.00     8391     3275

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     1.38      0.01     1.35     1.40 1.00     7944     2711

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

plot(fit_no_pooling)

8.7 Модель со случайным эффектом: добавляем предикторы

fit_hierarchical_intercept <- brm(speech_tempo ~ I(log(n_vowels)) + (1|speaker),
                                  data = df,
                                  family = gaussian(),
                                  cores = n_cores,
                                  silent = TRUE,
                                  prior = c(prior(normal(5.31, 1.93), class = "Intercept", lb = 0),
                                            prior(normal(2, 5), class = "sigma", lb = 0),
                                            prior(normal(2, 10), class = "sd")),
                                                              
                                  iter = 4000,
                                  control = list(adapt_delta = .85))

fit_hierarchical_intercept

 Family: gaussian 
  Links: mu = identity 
Formula: speech_tempo ~ I(log(n_vowels)) + (1 | speaker) 
   Data: df (Number of observations: 5623) 
  Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
         total post-warmup draws = 8000

Multilevel Hyperparameters:
~speaker (Number of levels: 9) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.82      0.27     0.47     1.49 1.00     1251     2176

Regression Coefficients:
             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept        1.92      0.28     1.37     2.50 1.00     1161     1896
Ilogn_vowels     1.14      0.02     1.10     1.18 1.00     4568     4384

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     1.09      0.01     1.07     1.11 1.00     4255     3545

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

plot(fit_hierarchical_intercept)

fit_hierarchical_intercept |> 
  conditional_effects() |> 
  plot(points = TRUE)

Cлучайные члены для каждого носителя

ranef(fit_hierarchical_intercept)$speaker[,,"Intercept"] |> 
  as.data.frame() |> 
  rownames_to_column("speaker") |> 
  ggplot(aes(y = speaker, x = Estimate)) +
  geom_errorbar(aes(xmin = `Q2.5`, xmax = `Q97.5`))+
  geom_vline(xintercept = 0, linetype = 2) +
  geom_point()

fit_hierarchical_intercept_slope <- brm(speech_tempo ~ I(log(n_vowels)) + (I(log(n_vowels)) + 1|speaker),
                                  data = df,
                                  family = gaussian(),
                                  cores = n_cores,
                                  silent = TRUE,
                                  prior = c(prior(normal(5.31, 1.93), class = "Intercept", lb = 0),
                                            prior(normal(2, 5), class = "sigma", lb = 0),
                                            prior(normal(2, 10), class = "sd")),
                                  iter = 4000,
                                  control = list(adapt_delta = .85))

Warning: There were 9 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.

Warning: Examine the pairs() plot to diagnose sampling problems

fit_hierarchical_intercept_slope

Warning: There were 9 divergent transitions after warmup. Increasing
adapt_delta above 0.85 may help. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup

 Family: gaussian 
  Links: mu = identity 
Formula: speech_tempo ~ I(log(n_vowels)) + (I(log(n_vowels)) + 1 | speaker) 
   Data: df (Number of observations: 5623) 
  Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
         total post-warmup draws = 8000

Multilevel Hyperparameters:
~speaker (Number of levels: 9) 
                            Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
sd(Intercept)                   0.96      0.32     0.54     1.77 1.00     2693
sd(Ilogn_vowels)                0.17      0.07     0.08     0.34 1.00     2735
cor(Intercept,Ilogn_vowels)    -0.35      0.32    -0.84     0.36 1.00     3502
                            Tail_ESS
sd(Intercept)                   3594
sd(Ilogn_vowels)                4394
cor(Intercept,Ilogn_vowels)     4622

Regression Coefficients:
             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept        1.83      0.33     1.16     2.50 1.00     1787     2635
Ilogn_vowels     1.17      0.06     1.05     1.30 1.00     3055     4217

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     1.08      0.01     1.06     1.10 1.00    10854     5436

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

plot(fit_hierarchical_intercept_slope)

fit_hierarchical_intercept_slope |> 
  conditional_effects() |> 
  plot(points = TRUE)

Cлучайные члены для каждого носителя

ranef(fit_hierarchical_intercept_slope)$speaker[,,"Intercept"] |> 
  as.data.frame() |> 
  rownames_to_column("speaker") |> 
  ggplot(aes(y = speaker, x = Estimate)) +
  geom_errorbar(aes(xmin = `Q2.5`, xmax = `Q97.5`))+
  geom_vline(xintercept = 0, linetype = 2) +
  geom_point()

Угловой кэффициент для каждого носителя

ranef(fit_hierarchical_intercept_slope)$speaker[,,"Ilogn_vowels"] |> 
  as.data.frame() |> 
  rownames_to_column("speaker") |> 
  ggplot(aes(y = speaker, x = Estimate)) +
  geom_errorbar(aes(xmin = `Q2.5`, xmax = `Q97.5`))+
  geom_vline(xintercept = 0, linetype = 2) +
  geom_point()

fit_hierarchical_slope <- brm(speech_tempo ~ I(log(n_vowels)) + (I(log(n_vowels)) + 0|speaker),
                              data = df,
                              family = gaussian(),
                              cores = n_cores,
                              silent = TRUE,
                              prior = c(prior(normal(5.31, 1.93), class = "Intercept", lb = 0),
                                        prior(normal(2, 5), class = "sigma", lb = 0),
                                        prior(normal(2, 10), class = "sd")),
                              iter = 4000,
                              control = list(adapt_delta = .85))

Warning: There were 4 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.

Warning: Examine the pairs() plot to diagnose sampling problems

fit_hierarchical_slope

Warning: There were 4 divergent transitions after warmup. Increasing
adapt_delta above 0.85 may help. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup

 Family: gaussian 
  Links: mu = identity 
Formula: speech_tempo ~ I(log(n_vowels)) + (I(log(n_vowels)) + 0 | speaker) 
   Data: df (Number of observations: 5623) 
  Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
         total post-warmup draws = 8000

Multilevel Hyperparameters:
~speaker (Number of levels: 9) 
                 Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Ilogn_vowels)     0.31      0.09     0.19     0.53 1.01      961     1621

Regression Coefficients:
             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept        1.92      0.05     1.82     2.01 1.00    10509     6068
Ilogn_vowels     1.13      0.10     0.93     1.33 1.01      848     1266

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     1.11      0.01     1.09     1.13 1.00     3279     2748

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

plot(fit_hierarchical_slope)

fit_hierarchical_slope |> 
  conditional_effects() |> 
  plot(points = TRUE)

Угловой кэффициент для каждого носителя

ranef(fit_hierarchical_slope)$speaker[,,"Ilogn_vowels"] |> 
  as.data.frame() |> 
  rownames_to_column("speaker") |> 
  ggplot(aes(y = speaker, x = Estimate)) +
  geom_errorbar(aes(xmin = `Q2.5`, xmax = `Q97.5`))+
  geom_vline(xintercept = 0, linetype = 2) +
  geom_point()

В работе (Ota et al. 2009) описывается эксперимент, в котором арабским и японским билингвам показывали английские слова и просили их определить, относятся ли они к одной семантической группе. В словах были в том числе и фонетически близкие слова, которые различаются фонемами, которых нет в L1. В датасет записан фрагмент данных, в которых стимулы не связаны. Постройте регрессионную модель, предсказывающую время реакции (Words.RT) на основе того, верно ли испытуемый (Subject) выполнил задание (Words.ACC). Я использовал sum contrast.

Coretta, Stefano. 2016. “Vowel Duration and Aspiration Effects in Icelandic.” University of York.

Gelman, A., and D. B. Rubin. 1992. “Inference from Iterative Simulation Using Multiple Sequences.” Statistical Science, 457–72.

Ota, Mitsuhiko, Robert J. Hartsuiker, and Sarah L. Haywood. 2009. “The KEY to the ROCK: Near-Homophony in Nonnative Visual Word Recognition.” Cognition 111 (2): 263–69. https://doi.org/https://doi.org/10.1016/j.cognition.2008.12.007.

Stepanova, S. 2011. “Russian Spontaneous Speech Rate-Based on the Speech Corpus of Russian Everyday Interaction.” ICPhS, 1902–5.

Vehtari, Aki, Andrew Gelman, and Jonah Gabry. 2017. “Practical Bayesian Model Evaluation Using Leave-One-Out Cross-Validation and WAIC.” Statistics and Computing 27 (5): 1413–32.

Vehtari, Aki, Daniel Simpson, Andrew Gelman, Yuling Yao, and Jonah Gabry. 2024. “Pareto Smoothed Importance Sampling.” Journal of Machine Learning Research 25 (72): 1–58.