Multivariate analysis using multi-level models as described by Gelman

Question

I have a data set I have been asked to provide an alternate analysis for in which 34 mothers answered a psychological questionnaire. The authors have currently conducted an exploratory analysis regressing each of the tool's four sub scales against 9 predictor variables (e.g. infant gestational age at birth), and would like me to repeat this analysis using a Bayesian regularized approach.

I originally had tried to find some WinBUGS/JAGS code that used a multivariate likelihood that I could adapt, but have been struggling. In this paper on p.207 (http://www.stat.columbia.edu/~gelman/research/published/multiple2f.pdf), Andrew Gelman describes the use of multilevel models to shrink estimates from multiple outcomes. I would like to adapt this approach as it would allow me to stick within the models provided by rstanarm which would help me finish this project and meet my deadline since it avoids fiddling with BUGS code.

Before going down the rstanarm route, I wanted to make sure I understood things correctly by running a simple analysis using lmer. These are the results from a simple regression of each sub scale against one predictor.

lm(formula = scale(B) ~ scale(GA), data = models)

Coefficients:
(Intercept)    scale(GA)  
 -1.745e-16   -1.105e-01


Call:
lm(formula = scale(C) ~ scale(GA), data = models)

Coefficients:
(Intercept)    scale(GA)  
  9.828e-17   -7.675e-02


Call:
lm(formula = scale(D) ~ scale(GA), data = models)

Coefficients:
(Intercept)    scale(GA)  
  3.099e-17   -1.479e-02


Call:
lm(formula = scale(E) ~ scale(GA), data = models)

Coefficients:
(Intercept)    scale(GA)  
 -9.262e-21    1.818e-01

From what I can gather from the notation provided in the manuscript. I re-arranged my data in long form (so each scale was one column) and specified my model as

Linear mixed model fit by REML ['lmerMod']
Formula: scale(score) ~ 1 + (scale(GA) | outcome)
   Data: t
REML criterion at convergence: 401.3155
Random effects:
 Groups   Name        Std.Dev. Corr
 outcome  (Intercept) 0.35637      
          scale(GA)   0.04067  1.00
 Residual             0.94897      
Number of obs: 144, groups:  outcome, 4
Fixed Effects:
(Intercept)  
  -0.009226  

 coef(model)

     scale(GA) (Intercept)
B -0.004745175  -0.0508061
C -0.049972606  -0.4471198
D  0.029366340   0.2481023
E  0.029562769   0.2498236

Which I read as saying that each subscale has it's own intercept and also has it's own relationship with GA, with each subscale as exchangeable. I am very new to trying to adapt multi-level models from published papers without a demo data set to replicate results on first, and this is my first random slope model so I wanted to run this by more experienced users to check if I've understood things correctly before I continue.

Answer 1

Disclaimer There are people here with more experience.

While a hierarchical model is definitely the superior approach, you can also leave out this additional difficulty, given that any Bayesian approach is "regularized" by the inclusion of priors. To say, from literally interpreting the question asked to you, I don't think you necessarily need to use multilevel analysis, and take a bit more time to understand hierarchical modelling. In that case, you can try copying the formulas used between lm and rstanarm and expect similar, but regularized estimates for the coefficients. How much regularization is applied, depends on your choice of priors. While there are defaults in rstanarm, it is important to think about this, and that requires considering the scaling of the data and also the meaning of the coefficients in the model.

To create a hierarchical model, it is not necessarily clear which variables you want to group together, or shrink. I am not sure which grouping you want to use in your model. I don't really understand grouping by outcome.

In your case, what kind of variable is the outcome? Is it an ordinal variable (Low, Medium, High)? In that case, I would suggest reading page 13 in https://cran.r-project.org/web/packages/brms/vignettes/brms_overview.pdf for a very quick example.