2

我对如何计算 Rstan 中 glm 的系数（分类变量）之间的差异有疑问。

例如，我在 R 中使用 iris 数据集来判断是否可以计算系数差异的后验分布。

首先，我进行了如下所示的基本 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

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

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

4

1 回答 1

0

您可以从任何适当的后验分布计算绘制（例如由 Stan 产生的那些）的任何函数（包括系数差异）。

于 2017-12-11T08:19:39.940 回答