r - 如何使用 rstanarm 以 APA 样式报告贝叶斯线性(混合)模型？

Question

我目前正在努力解决如何按照 APA-6 的建议报告 rstanarm::stan_lmer() 的输出。

首先，我将在常客方法中建立一个混合模型，然后尝试使用贝叶斯框架来做同样的事情。

这是获取数据的可重现代码:

library(tidyverse)
library(neuropsychology)
library(rstanarm)
library(lmerTest)

df <- neuropsychology::personality %>% 
  select(Study_Level, Sex, Negative_Affect) %>% 
  mutate(Study_Level=as.factor(Study_Level),
         Negative_Affect=scale(Negative_Affect)) # I understood that scaling variables is important

现在，让我们以“传统”方式拟合线性混合模型，以学习水平(受​​教育年限)作为随机因素来测试性别(男性/女性)对负面影响(消极情绪)的影响。

fit <- lmer(Negative_Affect ~ Sex + (1|Study_Level), df)
summary(fit)

输出如下:

Linear mixed model fit by REML t-tests use Satterthwaite approximations to degrees of
  freedom [lmerMod]
Formula: Negative_Affect ~ Sex + (1 | Study_Level)
   Data: df

REML criterion at convergence: 3709

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.58199 -0.72973  0.02254  0.68668  2.92841 

Random effects:
 Groups      Name        Variance Std.Dev.
 Study_Level (Intercept) 0.04096  0.2024  
 Residual                0.94555  0.9724  
Number of obs: 1327, groups:  Study_Level, 8

Fixed effects:
              Estimate Std. Error         df t value Pr(>|t|)    
(Intercept)    0.01564    0.08908    4.70000   0.176    0.868    
SexM          -0.46667    0.06607 1321.20000  -7.064 2.62e-12 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
     (Intr)
SexM -0.149

要报告它，我会说“我们拟合了一个线性混合模型，其中负面影响作为结果变量，性别作为预测变量，学习水平作为随机效应输入。在这个模型中，男性水平导致显着下降负面影响 (beta = -0.47, t(1321)=-7.06, p < .001)。

对吗？

然后，让我们尝试使用 rstanarm 在贝叶斯框架内拟合模型:

fitB <- stan_lmer(Negative_Affect ~ Sex + (1|Study_Level),
                  data=df,
                  prior=normal(location=0, scale=1), 
                  prior_intercept=normal(location=0, scale=1),
                  prior_PD=F)
print(fitB, digits=2)

返回:

stan_lmer
 family:  gaussian [identity]
 formula: Negative_Affect ~ Sex + (1 | Study_Level)
------

Estimates:
            Median MAD_SD
(Intercept)  0.02   0.10 
SexM        -0.47   0.07 
sigma        0.97   0.02 

Error terms:
 Groups      Name        Std.Dev.
 Study_Level (Intercept) 0.278   
 Residual                0.973   
Num. levels: Study_Level 8 

Sample avg. posterior predictive 
distribution of y (X = xbar):
         Median MAD_SD
mean_PPD 0.00   0.04  

------
For info on the priors used see help('prior_summary.stanreg').

我认为 median 是系数后验分布的中位数，mad_sd 相当于标准差。这些参数接近常客模型的 beta 和标准误，令人放心。但是，我不知道如何将输出形式化并用文字表示。

此外，如果我对模型进行总结 (summary(fitB, probs=c(.025, .975), digits=2))，我会得到后验分布的其他特征:

...
Estimates:
                                             mean     sd       2.5%     97.5% 
(Intercept)                                    0.02     0.11    -0.19     0.23
SexM                                          -0.47     0.07    -0.59    -0.34
...

像下面这样的东西好吗？

“我们在贝叶斯框架内拟合了一个线性混合模型，其中负面影响作为结果变量，性别作为预测变量，研究水平作为随机效应输入。系数和截距的先验设置为正常值(平均值=0， sd=1)。在此模型中，与男性水平相关的系数的后验分布特征表明负面影响减少(平均值 = -0.47，sd = 0.11，95% CI[-0.59，-0.34]) .

感谢您的帮助。

Answer 1

以下是个人观点，心理学期刊可能会接受也可能不会接受。

To report it, I would say that "we fitted a linear mixed model with negative affect as outcome variable, sex as predictor and study level was entered as a random effect. Within this model, the male level led to a significant decrease of negative affect (beta = -0.47, t(1321)=-7.06, p < .001).

Is that correct?

从频率论者的角度来看，这被认为是正确的。

从贝叶斯角度来看，关键概念是(当然，以模型为条件)

1. 真实效果小于后验中位数的概率为 0.5，真实效果大于后验中位数的概率为 0.5。频率论者倾向于将后验中位数视为数值最优值。
2. posterior_interval 函数以默认概率 0.9 生成围绕中位数的可信区间(尽管较小的数字会生成更准确的边界估计)。因此，您可以合理地说真实效果在这些界限之间的概率为 0.9。频率论者倾向于将置信区间视为可信区间。
3. as.data.frame 函数将使您能够访问原始绘图，因此 mean(as.data.frame(fitB)$male > 0) 产生在同一研究中男性和女性之间的预期结果差异为正的概率。频率论者倾向于将这些概率视为类似于 p 值。

对于贝叶斯方法，我会说

We fit a linear model using Markov Chain Monte Carlo with negative affect as the outcome variable, sex as predictor and the intercept was allowed to vary by study level.

然后用上面三个概念谈谈估算。