如何计算 Rstan 中 glm 程序的系数（分类变量）之间的差异

Question

我有一个关于如何在 Rstan 中计算 glm 系数（分类变量）之间差异的问题。

比如我用R中的鸢尾花数据集来判断是否可以计算系数差的后验分布

首先，我进行了如下所示的基本 glm 程序并计算了系数的显着差异。

library(tidyverse)
library(magrittr)
library(multcomp)

iris_glm <-
glm(Sepal.Length ~ Species, data = iris) 

multcomp::glht(iris_glm, linfct = mcp(Species = "Tukey")) %>% 
summary(.) %>% 
broom::tidy() 

                     lhs rhs estimate std.error statistic      p.value
1    versicolor - setosa   0    0.930 0.1029579  9.032819 0.000000e+00
2     virginica - setosa   0    1.582 0.1029579 15.365506 0.000000e+00
3 virginica - versicolor   0    0.652 0.1029579  6.332686 4.294805e-10

接下来，我使用如下代码中的 stan 执行了贝叶斯 glm 程序，并计算了生成量部分中系数之间差异的后验分布。

# Make the model matrix for Rstan
iris_mod <-
model.matrix(Sepal.Length ~ Species, data = iris) %>%
as.data.frame(.)

# Input data
stan_data <-
list(N = nrow(iris_mod),
SL = iris$Sepal.Length,
Intercept = iris_mod$`(Intercept)`,
versicolor = iris_mod$Speciesversicolor,
virginica = iris_mod$Speciesvirginica)

# Stan code

data{
int N;
real <lower = 0> SL[N];
int <lower = 1> Intercept[N];
int <lower = 0, upper = 1> versicolor[N];
int <lower = 0, upper = 1> virginica[N];
}

parameters{
real beta0;
real beta1;
real beta2;
real <lower = 0> sigma;
}

transformed parameters{
real mu[N];
for(n in 1:N) mu[n] = beta0*Intercept[n] + beta1*versicolor[n] +       
beta2*virginica[n];
}

model{
for(n in 1:N) SL[n] ~ normal(mu[n], sigma);
}

generated quantities{
real diff_beta0_beta1;
real diff_beta1_beta2;
real diff_beta0_beta2;

diff_beta0_beta1 = (beta0 + beta1) - beta0;
diff_beta1_beta2 = (beta0 + beta1) - (beta0 + beta2);
diff_beta0_beta2 = (beta0 + beta2) - beta0;
}


library(rstan)
fit_stan <- 
stan(file = "iris.stan", data = stan_data, chains = 4, 
seed = 1234)

# confirmation of posterior distribution
print(fit_stan, pars = c("diff_beta0_beta1", "diff_beta1_beta2", 
"diff_beta0_beta2"))

                  mean se_mean  sd  2.5%   25%   50%   75% 97.5% n_eff Rhat
diff_beta0_beta1 0.92       0 0.1 0.73 0.86 0.92 0.99  1.13  2041    1
diff_beta1_beta2 0.65       0 0.1 0.45 0.58 0.65 0.72  0.86  4000    1
diff_beta0_beta2 1.58       0 0.1 1.38 1.51 1.58 1.64  1.78  1851    1

最后，我可以在频率论方法和贝叶斯方法之间得到相同的结果。我认为这是正确的方法，但我不确定，因为没有信息或示例。此外，我还确认这种方式可以扩展到另一个误差分布（包括泊松分布、伽马分布、二项分布、负二项分布等）。

如果还有其他好的方法或者建议，请教教我

Answer 1

您可以根据任何适当的后验分布计算绘图（例如 Stan 生成的函数）的任何函数（包括系数的差异）。

如何计算 Rstan 中 glm 程序的系数（分类变量）之间的差异

How to calculate differences between the coefficients (categorical variables) of glm procedure in Rstan

comparison

r

bayesian

glm

rstan