Диагностика линейной модели для байесовских моделей с использованием rstan

Question

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

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

# ordinary linear model diagnostics, similar to my use-case 
library(dplyr) 
library(purrr) 
library(tidyr) 
library(broom) 
# make linear models for each type of species 
xy <- 
    iris %>% 
    nest(-Species) %>% 
    mutate(model = map(data, 
        ~lm(Sepal.Length ~ Petal.Width, 
         data = .)), 
     fit = map(model, augment)) 

Здесь мы имеем фрейм данных с вложенными списками, то model столбец содержит линейную модель для каждого вида:

> xy 
# A tibble: 3 × 4 
    Species    data model     fit 
     <fctr>   <list> <list>    <list> 
1  setosa <tibble [50 × 4]> <S3: lm> <data.frame [50 × 9]> 
2 versicolor <tibble [50 × 4]> <S3: lm> <data.frame [50 × 9]> 
3 virginica <tibble [50 × 4]> <S3: lm> <data.frame [50 × 9]> 

broom::augment функция позволила нам добавить значения расстояния повару для каждой точки данных в этом кадре данных, и мы можем проверить их так:

# inspect Cook's distance values 
xy %>% 
unnest(fit) %>% 
arrange(desc(.cooksd)) 

    # A tibble: 150 × 10 
     Species Sepal.Length Petal.Width .fitted .se.fit  .resid  .hat .sigma .cooksd 
     <fctr>  <dbl>  <dbl> <dbl>  <dbl>  <dbl>  <dbl>  <dbl>  <dbl> 
1 versicolor   5.9   1.8 6.612097 0.16181001 -0.7120969 0.13725081 0.4269862 0.24507448 
2  setosa   5.0   0.6 5.335281 0.17114108 -0.3352811 0.25027563 0.3410686 0.21385214 
3 virginica   4.9   1.7 6.375829 0.13613717 -1.4758292 0.04875277 0.5826838 0.15434787 
4  setosa   5.7   0.4 5.149247 0.08625887 0.5507534 0.06357957 0.3355980 0.09396588 
5  setosa   4.3   0.1 4.870195 0.08321347 -0.5701948 0.05916942 0.3349111 0.09285408 
6 virginica   5.8   2.4 6.831411 0.14828703 -1.0314106 0.05784319 0.6035012 0.09117693 
7 virginica   7.2   1.6 6.310746 0.16207266 0.8892538 0.06909799 0.6084108 0.08293253 
8 versicolor   4.9   1.0 5.471005 0.11998077 -0.5710051 0.07546185 0.4328038 0.07544526 
9  setosa   5.8   0.2 4.963212 0.05287342 0.8367879 0.02388828 0.3228858 0.07500610 
10 versicolor   6.0   1.0 5.471005 0.11998077 0.5289949 0.07546185 0.4340307 0.06475225 
# ... with 140 more rows, and 1 more variables: .std.resid <dbl> 

И с методом autoplot, мы можем сделать информативные диагностические графики, которые показывают значение расстояния Кука, и помочь нам быстро определить точки данных с негабаритным эффектом на paramters отделки моделей:

# plot model diagnostics 
library(ggplot2) 
library(ggfortify) 
diagnostic_plots_df <- 
    xy %>% 
    mutate(diagnostic_plots = map(model, 
           ~autoplot(., 
              which = 1:6, 
              ncol = 3, 
              label.size = 3))) 

Вот только один из участков производства:

> diagnostic_plots_df[[1]] 

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

# Bayesian linear model diagnostics 
library(rstanarm) 
bayes_xy <- 
    iris %>% 
    nest(-Species) %>% 
    mutate(model = map(data, 
        ~stan_lm(Sepal.Length ~ Petal.Width, 
         data = ., 
         prior = NULL, 
         chains = 1, 
         cores = 2, 
         seed = 1)), 
     fit = map(model, augment)) 

> bayes_xy 
# A tibble: 3 × 4 
    Species    data   model     fit 
     <fctr>   <list>  <list>    <list> 
1  setosa <tibble [50 × 4]> <S3: stanreg> <data.frame [50 × 5]> 
2 versicolor <tibble [50 × 4]> <S3: stanreg> <data.frame [50 × 5]> 
3 virginica <tibble [50 × 4]> <S3: stanreg> <data.frame [50 × 5]> 

Но метод broom::augment не имеет ничего подобное расстояние значения повара:

# inspect fit diagnostics 
bayes_xy %>% unnest(fit) 

# A tibble: 150 × 6 
    Species Sepal.Length Petal.Width .fitted .se.fit  .resid 
    <fctr>  <dbl>  <dbl> <dbl>  <dbl>  <dbl> 
1 setosa   5.1   0.2 4.963968 0.06020298 0.13482025 
2 setosa   4.9   0.2 4.963968 0.06020298 -0.06517975 
3 setosa   4.7   0.2 4.963968 0.06020298 -0.26517975 
4 setosa   4.6   0.2 4.963968 0.06020298 -0.36517975 
5 setosa   5.0   0.2 4.963968 0.06020298 0.03482025 
6 setosa   5.4   0.4 5.151501 0.11299956 0.21818386 
7 setosa   4.6   0.3 5.057734 0.05951488 -0.47349794 
8 setosa   5.0   0.2 4.963968 0.06020298 0.03482025 
9 setosa   4.4   0.2 4.963968 0.06020298 -0.56517975 
10 setosa   4.9   0.1 4.870201 0.11408783 0.04313845 
# ... with 140 more rows 

И нет autoplot метода:

# plot model diagnostics 
bayes_diagnostic_plots_df <- 
    bayes_xy %>% 
    mutate(diagnostic_plots = map(model, 
           ~autoplot(., 
              which = 1:6, 
              ncol = 3, 
              label.size = 3))) 

# Error, there doesn't seem to be an autoplot method for stan_lm output 

shinystan::launch_shinystan(bayes_xy$model[[1]]) 

# This is quite interesting, but nothing at the level of specific data points 

Некоторые из научной литературы говорит о методах, таких как model perturbations, phi-divergence, Cook's posterior mode distance and Cook's posterior mean distance, Kullback-Leibler divergence, etc., etc.. Но я не вижу нигде, что это было изучено с помощью R-кода, и я застрял.

Существует an unanswered question по этой теме в Cross-validated. Я размещаю здесь, потому что я ищу идеи о написании кода для расчета статистики влияния (а не советы о статистической теории и подходах и т. Д., Которые должны идти по этому другому вопросу)

Как я могу получить что-то как измерение расстояния Кука от выхода rstanarm::stan_lm?

Answer 1

Осмотрев часть обсуждения списка адресов электронной почты stan-users, я увидел, что вывод из пакета loo содержит «точечные вклады для каждого наблюдения».Так вот попытка работать с теми:

# Bayesian linear model diagnostics 
library(rstanarm) 
library(loo) 
bayes_xy <- 
    iris %>% 
    nest(-Species) %>% 
    mutate(model = map(data, 
        ~stan_lm(Sepal.Length ~ Petal.Width, 
         data = ., 
         prior = NULL, 
         chains = 1, 
         cores = 2, 
         seed = 1))) 


bayes_xy_loo <- 
bayes_xy %>% 
    mutate(loo_out = map(model, ~loo(.))) 

library(ggplot2) 
library(ggrepel) 
n <- 5 # how many points to label 


my_plots <- 
bayes_xy_loo %>% 
    select(loo_out) %>% 
    mutate(loo_pointwise = map(.$loo_out, ~data.frame(.$pointwise))) %>% 
    mutate(plots = map(.$loo_pointwise, 
     ~ggplot(., 
     aes(elpd_loo, 
      looic)) + 
     geom_point(aes(size = p_loo)) + 
     geom_text_repel(data = .[tail(order(.$looic), n),] , 
        aes(label = row.names(.[tail(order(.$looic), n),])), 
        nudge_x = -0.1, 
        nudge_y = -0.3) + 
     geom_label_repel(data = .[tail(order(.$elpd_loo), n),] , 
         aes(label = row.names(.[tail(order(.$elpd_loo), n),])), 
         nudge_x = 0.1, 
         nudge_y = 0.1) + 
     xlab("Expected log pointwise \npredictive density") + 
     ylab("LOO information \ncriterion") + 
     scale_size_area(name = "Estimated \neffective\nnumber \nof parameters") + 
     theme_minimal(base_size = 10))) 

do.call(gridExtra::grid.arrange, my_plots$plots) 

Однако точки, предложенные быть влиятельным не хороший матч. Например, в вопросе мы имеем обс. 7, 15 и 30 с высокими значениями расстояния Кука. В выходе loo кажется, что только общ. 15 обозначается как обычно. Поэтому, возможно, это неправильный способ сделать это.

Answer 2

Это post сказал Аки Vehtari сказал это лучше всего:

Разница между lppd_i и loo_i был использован в качестве меры чувствительности (смотри, например, Гельфанд и др, 1992). Вероятность оценки параметра формы Парето k, вероятно, будет велика, если разница между lppd_i и loo_i велика. Пока еще нет ясно, будет ли оценка параметра формы Парето k лучше, чем lppd_i-loo_i, но по крайней мере мы знаем, что оценка для lppd_i-loo_i слишком мала , если k близок к 1 или более, поэтому он может лучше посмотреть на k. В примере с обычной моделью k для одного наблюдения велико, но с ученик-t модель k меньше. Обычная модель такая же, как у модели student-t, но с очень сильная до степени свободы. Поэтому речь идет не только о наличии сильной предшествующей или большей усадки, но имеющей модель, которая может описывать наблюдения . С увеличенной усадкой и ненадежной моделью наблюдения одно замечание все еще может быть удивительным. Естественно, это не всегда лучшее решение для перехода на более надежную модель наблюдения , позволяющую «выбросы». Вместо этого было бы лучше сделать функцию регрессии более нелинейной (которая имеет менее сильную предысторию), или преобразовать ковариаты, или добавить больше ковариатов. Так что я рекомендую посмотреть параметры параметра формы Парето, но я не рекомендую увеличить усадку, если значения велики.

Вы можете получить парето оценку параметра формы к из $pareto_k элемента списка, созданного с помощью loo функции в пакете Лоо, который реэкспорт в пакете rstanarm. Если это значение выше 0,7 (по умолчанию), функция loo будет рекомендовать вам переопределить модель, оставляя это наблюдение, потому что заднее распределение, вероятно, будет слишком чувствительным к этому наблюдению, чтобы удовлетворить предположению LOOIC, что каждый наблюдение оказывает незначительное влияние на последующее распределение.

В случае OP только седьмое наблюдение имеет оценку параметра формы Парето, которая немного больше 0,5, так что наблюдение, вероятно, не оказывает чрезвычайного влияния на заднюю. Но вы определенно хотите исследовать наблюдения, которые имеют значение больше 1.0

Вы можете также вызвать метод plot для объекта Loo, особенно с опцией label_points = TRUE, отличной от значения, для визуализации параметров параметров формы Парето.