如何使用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

为了报告,我会说"我们拟合了一个线性混合模型,其中负面影响作为结果变量,性别作为预测因子,研究水平作为随机效应输入.在此模型中,男性水平导致负面影响显着减少(β= -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

以下是心理学期刊可能接受或不接受的个人意见.

为了报告,我会说"我们拟合了一个线性混合模型,其中负面影响作为结果变量,性别作为预测因子,研究水平作为随机效应输入.在此模型中,男性水平导致负面影响显着减少(β= -0.47,t(1321)= - 7.06,p <.001).

那是对的吗？

从频率论的角度来看,这被认为是正确的.

贝叶斯视角的关键概念是(当然,以模型为条件)

1. 真实效应小于后中位数的概率为0.5,真实效应大于后中位数的概率为0.5.频繁的人倾向于将后验中位数视为数值最优.
2. 该posterior_interval函数在中值附近产生可信区间,默认概率为0.9(尽管较低的数字产生更准确的边界估计).所以,你可以合理地说,在这些界限之间存在真实效果的概率为0.9.频繁的人往往认为置信区间就像可信区间一样.
3. 该as.data.frame功能可以让您访问原始抽奖,因此mean(as.data.frame(fitB)$male > 0)可以得出同一研究中男女结果的预期差异为正的概率.频繁的人倾向于将这些概率视为p值.

我会说,对于贝叶斯方法

我们使用马尔可夫链蒙特卡罗拟合线性模型,其中负面影响作为结果变量,性别作为预测因子,截距允许根据研究水平而变化.

然后使用上面的三个概念来讨论估算.