# Multivariate analysis using multi-level models as described by Gelman

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.

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.

• My variable is a subscale score that has been centered and divided by the sd. In the article I attached Gelman (seems) to use this approach in a typical eight schools example where each participant was administered a battery of tests. Gelman scales all of the scores, and then allows for random intercepts for school and outcome, and then a random slope that is allowed to vary by outcome and site. The idea (which I think I understand) is that you're shrinking the outcomes together to help guard against issues with multiple comparisons. Does that make sense? – Tim Sep 17 '17 at 10:49
• It's more like you are shrinking parameters of the model (intercepts and coefficients), not the outcomes. In the case of the schools examples, you group the treatment effect variable by schools. It does help against multiple comparison, yes. – Gijs Sep 17 '17 at 11:34
• I think the schools example is a bit confusing to start out with actually, as you are indeed combining the outcomes of more experiments with it. – Gijs Sep 17 '17 at 11:38
• You're right, I am thinking about shrinking parameters of the model. The schools example I'm referring to is the one in the linked article where each school had multiple outcomes as well (several different questionnaires). It sounds like you're saying I have parameterised this correctly to achieve the goal of shrinking estimates of the effect of covariates, is that right? – Tim Sep 17 '17 at 18:38
• Actually the more I think of this the more I think it just needs to wait until I have the time to get more into the reading and for now just use rstanarm to run the four regressions separately. I think I should be scaling each outcome before gathering rows, but even then when I try and replicate no pooling, complete pooling, and partial pooling the results aren't making sense i.e. coefficients for outcomes that are below the mean are being shrunk down instead of brought in to the overall mean. – Tim Sep 17 '17 at 19:30