ベイジアンモデルの相互作用効果をプロットする（rstanarmを使用）

Question

rstanarm（）のベイジアンリニアモデルの別の変数の値がどのように変化するかを示したいと考えています。私はモデルを適合させることができ、各パラメータの推定値を見るために後部から引き出すことができますが、他の変化や関連する不確実性として、相互作用における1つの変数の影響のある種のプロットを与える方法は明確ではありません（すなわち、限界効果プロット）。以下は私の試みである。この（線形、ガウス、アイデンティティのリンク、無インターセプト）の場合

library(rstanarm) 

# Set Seed 
set.seed(1) 

# Generate fake data 
w1 <- rbeta(n = 50, shape1 = 2, shape2 = 1.5) 
w2 <- rbeta(n = 50, shape1 = 3, shape2 = 2.5) 

dat <- data.frame(y = log(w1/(1-w1)), 
        x = log(w2/(1-w2)), 
        z = seq(1:50)) 

# Fit linear regression without an intercept: 
m1 <- rstanarm::stan_glm(y ~ 0 + x*z, 
         data = dat, 
         family = gaussian(), 
         algorithm = "sampling", 
         chains = 4, 
         seed = 123, 
         ) 


# Create data sets with low values and high values of one of the predictors 
dat_lowx <- dat 
dat_lowx$x <- 0 

dat_highx <- dat 
dat_highx$x <- 5 

out_low <- rstanarm::posterior_predict(object = m1, 
            newdata = dat_lowx) 

out_high <- rstanarm::posterior_predict(object = m1, 
             newdata = dat_highx) 

# Calculate differences in posterior predictions 
mfx <- out_high - out_low 

# Somehow get the coefficients for the other predictor? 

Answer 1

、

mu = beta_x * x + beta_z * z + beta_xz * x * z 
    = (beta_x + beta_xz * z) * x 
    = (beta_z + beta_xz * x) * z 

ので、xまたはzの限界効果をプロットする、あなただけの必要次に、あなたが

post <- as.data.frame(m1) 

を介して取得することができるそれぞれの係数の事後分布の適切な範囲、

dmu_dx <- post[ , 1] + post[ , 3] %*% t(sort(dat$z)) 
dmu_dz <- post[ , 2] + post[ , 3] %*% t(sort(dat$x)) 

そして、あなたは、あなたのデータの各観測とmuにzの影響のためにのためにmuにxの影響を計算し、以下のようなものを、使用して、データの各観察のための単一の限界効果を推定することができます各観測。

colnames(dmu_dx) <- round(sort(dat$x), digits = 1) 
colnames(dmu_dz) <- dat$z 
bayesplot::mcmc_intervals(dmu_dz) 
bayesplot::mcmc_intervals(dmu_dx) 

この場合、列名は単に観測値です。

Answer 2

ggeffects-packageのいずれかを使用することもできます。 （マージナルエフェクトの場合はsjPlotは、ggeffects）の関数をラップします。

相互作用の限界効果をプロットするには、type = "int"でsjPlot::plot_model()を使用します。継続的なモデレータ変数にプロットする値を定義するにはmdrt.valuesを使用し、ppdを使用して、予測を線形予測子の事後分布に基づいて行うか、または後予測分布から描画するかを指定します。

library(sjPlot) 
plot_model(m1, type = "int", terms = c("x", "z"), mdrt.values = "meansd") 

plot_model(m1, type = "int", terms = c("x", "z"), mdrt.values = "meansd", ppd = TRUE) 

または、他の特定の値に限界効果をプロットtype = "pred"を使用しterms -argumentの値を指定する：

plot_model(m1, type = "pred", terms = c("x", "z [10, 20, 30, 40]")) 
# same as: 
library(ggeffects) 
dat <- ggpredict(m1, terms = c("x", "z [10, 20, 30, 40]")) 
plot(dat) 

さらに多くのオプションがあり、プロットの外観をカスタマイズするさまざまな方法があります。関連するヘルプファイルとパッケージのビネットを参照してください。