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.