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方法で、我々はクックの距離値を示して有益な診断プロットを作り、かつ迅速にモデルの偶然に上の特大影響でデータポイントを識別するために私たちを助けることができます。

# 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 

ななどのメソッドに関する学術文献の交渉の一部,phi-divergence, Cook's posterior mode distance and Cook's posterior mean distance, Kullback-Leibler divergence,etc.,etc.。しかし、これはRのコードで探検された場所はどこにも見えず、私は立ち往生しています。

このトピックのクロスバリデーションはan unanswered questionです。私は影響統計を計算するためのコードを書くアイデアを探しているからです（統計的な理論やアプローチなどに関するアドバイスではなく、他の質問に行くべきです）

rstanarm::stan_lmの出力からCook's Distance measureのように？

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) 

は、しかし、影響力のあることが示唆ポイントは良い試合ではありません。たとえば、質問にはobsがあります。クックの距離値が高い7,15、および30。 loo出力では、obsのみと思われます。 15は通常どおりに識別されます。だから、これは正しい方法ではないかもしれません。

Answer 2

このpostはアキVehtariによって前記の最良それを言った：

をlppd_iとloo_i差が感度指標 として使用されている（例えば、Gelfand他1992参照します）。 lppd_iとloo_iの差が大きい場合、パレート形状パラメータの推定値kは になる可能性があります。まだわかりません パレート形状パラメータの推定kが lppd_i-loo_iよりも良いかどうかはわかりませんが、少なくともlppd_i-loo_iの見積もりは小さすぎることがわかります kが1に近いかそれ以上であればkを見る方が良いです。 ノーマルモデルのスタック例では、1回の観測ではkが大きくなりますが、student-tモデルkは で小さくなります。通常のモデルはstudent-tモデルと同じですが、自由度に先立って が非常に強いモデルです。だから、それ以前のものよりも強いものを持っているだけではありませんが、観察を記述できるモデルを持っているのは、 です。収縮が大きくなり、非堅牢な観測モデルでは、 の観測は依然として驚くべきことです。 当然「異常値」を許容するより頑強な 観測モデルに変更するのが最良の解決策ではありません。その代わりに、より非線形である（以前はあまり強くない）か、 共変を変換するか、より多くの共変量を追加する方がよいでしょう。 パレート形状のパラメータ値を調べることをおすすめしますが、値が大きい場合は縮みを増やすことをお勧めしません。

あなたはrstanarmパッケージによって再輸出さLOOパッケージのloo関数によって生成されるリストの$pareto_k要素からパレート形状パラメータ推定値kを得ることができます。この値が0.7（デフォルト）より大きい場合、loo関数は、この観測値を残してモデルを再フィットすることを推奨します。なぜなら、事後分布がこの観測値に対してあまりにも敏感で、LOOICの仮定観察は事後分布にはごくわずかな影響しか及ぼさない。

OPのケースでは、第7の観測値だけが0.5よりわずかに大きいパレート形状パラメータ推定値を有するため、観察はおそらく後部に対して極端な影響を与えない。しかし、あなたは間違いなくあなたはまた、特にパレート形状パラメータ推定値を視覚化するために、デフォルト以外のオプションlabel_points = TRUEで、トイレオブジェクトのplotメソッドを呼び出すことができます1.0

より大きい値を持つ観測を調査したいです。