Включение априорной информации в регрессию хребта (RAPM) в R

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

Вызов glmnet выглядит следующим образом:

glmnet(y = shot_rates, x = dummy_matrix, weights = shift_length, lambda = previously_obtained_lambda)

(Лямбда получается путем перекрестной проверки того же набора данных, что также выполняется с использованием glmnet.)

На данный момент распределение полностью гауссово, и каждая предикторная переменная смещена к среднему значению, равному нулю. Я хочу включить предварительную информацию (предварительные средства) для каждой фиктивной переменной и, возможно, установить отдельные значения лямбда для каждой из них, но я не уверен, как мне это сделать. Я считаю, что могу использовать Penal.factors для настройки лямбда-выражений для каждой переменной, чтобы мы могли пока отложить это в сторону и сосредоточиться на предыдущих средствах.

Я изучил использование пакета bayesglm и реализацию prior.means, но мои проблемы с ним двоякие: он медленный и не принимает разреженные матрицы, что значительно замедляет работу. Для справки: моя матрица фиктивных переменных содержит примерно 600 000 строк и примерно 2 000 столбцов.

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


person Topdownhockey    schedule 26.11.2020    source источник
comment
Есть ли причина, по которой вы хотите использовать гребневую регрессию вместо построения полностью байесовской модели?   -  person sjp    schedule 27.11.2020
comment
Честно говоря, я изучил этот метод с помощью гребневой регрессии, и я обнаружил, что код довольно эффективен. Теоретически я не против полностью байесовской модели, но, по общему признанию, я не эксперт по байесовской статистике, и, что более важно, мой ранний опыт работы с байесовской моделью заключался в том, что она была бы довольно медленной, что далеко не идеально для меня, поскольку я хотел бы использовать метод проб и ошибок, чтобы усовершенствовать мой процесс. Однако я открыт для предложений; что вы имеете в виду для байесовской модели?   -  person Topdownhockey    schedule 27.11.2020


Ответы (1)


Итак, основываясь на комментарии, кажется, что в принципе байесовский подход вас устраивает, и это единственный известный мне способ иметь регуляризирующие априорные значения, не центрированные на 0. Вы также упомянули, что скорость была проблемой, поэтому я бы рекомендуют подбирать модель с использованием Стэна, который, как правило, намного быстрее, чем другие байесовские методы. Кроме того, brms и Stan имеют замечательную документацию, которая намного полезнее, чем та, которую вы обычно найдете для других статистических пакетов в R.

brms — очень полезный пакет, который позволяет подгонять модели Стэна к синтаксису, подобному lme4.

В brms можно указать априорные значения для перехвата и каждой независимой переменной следующим образом:

model_priors <- c(
  prior(normal(0, 5), class = "Intercept"),
  prior(normal(0, 1), class = "b")
)

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

model_priors <- c(
  prior(normal(0, 5), class = "Intercept"),
  prior(normal(0.5, 1), class = "b", coef = "first_predictor"),
  prior(normal(-0.5, 1), class = "b", coef = "second_predictor")
)

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

Затем вы можете включить эти априоры в модель примерно так:

modelfit <- brm(
  formula = outcome_variable ~ first_predictor + second_predictor, 
  data = df,
  prior = model_priors,
)
person sjp    schedule 26.11.2020