title: How to report with APA style a Bayesian Linear (Mixed) Models using rstanarm?
body: <p>I'm currently struggling with how to report, following APA-6 recommendations, the output of <code>rstanarm::stan_lmer()</code>. </p>
  
  <p>First, I'll fit a mixed model within the frequentist approach, then will try to do the same using the bayesian framework. </p>
  
  <p>Here's the reproducible code to get the data:</p>
  
  <pre><code>library(tidyverse)
  library(neuropsychology)
  library(rstanarm)
  library(lmerTest)
  
  df &lt;- neuropsychology::personality %&gt;% 
    select(Study_Level, Sex, Negative_Affect) %&gt;% 
    mutate(Study_Level=as.factor(Study_Level),
           Negative_Affect=scale(Negative_Affect)) # I understood that scaling variables is important
  </code></pre>
  
  <p>Now, let's fit a linear mixed model in the "traditional" way to test the impact of Sex (male/female) on Negative Affect (negative mood) with the study level (years of education) as random factor.</p>
  
  <pre><code>fit &lt;- lmer(Negative_Affect ~ Sex + (1|Study_Level), df)
  summary(fit)
  </code></pre>
  
  <p>The output is the following:</p>
  
  <pre><code>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(&gt;|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
  </code></pre>
  
  <p>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 &lt; .001).</p>
  
  <p>Is that correct?</p>
  
  <p>Then, let's try to fit the model within a bayesian framework using <code>rstanarm</code>:</p>
  
  <pre><code>fitB &lt;- 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)
  </code></pre>
  
  <p>This returns:</p>
  
  <pre><code>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').
  </code></pre>
  
  <p>I think than <code>median</code> is the median of the posterior distribution of the coefficient and <code>mad_sd</code> the equivalent of standart deviation. These parameters are close to the beta and standart error of the frequentist model, which is reassuring. However, I do not know how to formalize and put the output in words.</p>
  
  <p>Moreover, if I do the summary of the model (<code>summary(fitB, probs=c(.025, .975), digits=2)</code>), I get other features of the posterior distribution: </p>
  
  <pre><code>...
  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
  ...
  </code></pre>
  
  <p>Is something like the following good?</p>
  
  <p>"we fitted a linear mixed model within the bayesian framework with negative affect as outcome variable, sex as predictor and study level was entered as a random effect. Priors for  the coefficient and the intercept were set to normal (mean=0, sd=1). Within this model, the features of the posterior distribution of the coefficient associated with the male level suggest a decrease of negative affect (mean = -0.47, sd = 0.11, 95% CI[-0.59, -0.34]).</p>
  
  <p>Thanks for your help.</p>