Interpreting output from a stan_mvmer

Interfaces rstanarm
1

I am fitting Bayesian multivariate generalized linear models with correlated group-specific terms via Stan to model the binary outcomes for three outcomes; y1, y1 and y3 with a single predictor. In addition I have two grouping variables; hhid and years. The model is specified as below

model <- stan_mvmer(
formula = list(
y1bin ~ wealthindex + (1 | years) + (1 | hhid)
, y2bin ~ wealthindex + (1 | years) + (1 | hhid)
, y3bin ~ wealthindex + (1 | years) + (1 | hhid)
)
, data = sim_dflist[[1]]
, family = list(binomial, binomial, binomial)
, …
)

Below are some of the outputs (estimates from the model fit)

Estimates:
mean sd
y1|wealthindex 0.452 0.266
y1|(Intercept) 0.337 0.142
y1|mean_PPD 0.578 0.012
y2|(Intercept) 0.322 0.092
y2|wealthindex 0.803 0.203
y2|mean_PPD 0.581 0.012
y3|(Intercept) 0.423 0.174
y3|wealthindex 0.413 0.303
y3|mean_PPD 0.595 0.012
Sigma[years:y1|(Intercept),y1|(Intercept)] 0.352 0.154
Sigma[years:y2|(Intercept),y1|(Intercept)] 0.035 0.063
Sigma[years:y3|(Intercept),y1|(Intercept)] 0.036 0.116
Sigma[years:y2|(Intercept),y2|(Intercept)] 0.137 0.065
Sigma[years:y3|(Intercept),y2|(Intercept)] 0.103 0.082
Sigma[years:y3|(Intercept),y3|(Intercept)] 0.550 0.231
Sigma[hhid:y1|(Intercept),y1|(Intercept)] 0.311 0.069
Sigma[hhid:y2|(Intercept),y1|(Intercept)] 0.060 0.032
Sigma[hhid:y3|(Intercept),y1|(Intercept)] 0.094 0.050
Sigma[hhid:y2|(Intercept),y2|(Intercept)] 0.069 0.033
Sigma[hhid:y3|(Intercept),y2|(Intercept)] 0.106 0.038
Sigma[hhid:y3|(Intercept),y3|(Intercept)] 0.452 0.089
log-posterior -6394.996 22.838

Groups
Groups Name Std.Dev. Corr
years y1|(Intercept) 0.59344
y2|(Intercept) 0.36975 0.158
y3|(Intercept) 0.74185 0.083 0.376
hhid y1|(Intercept) 0.55775
y2|(Intercept) 0.26280 0.411
y3|(Intercept) 0.67248 0.250 0.602

Question?

  1. What’s the interpretation difference between the Sigma in Estimates (first table), e.g., Sigma[years:y1|(Intercept),y1|(Intercept)] = 0.352 and Std.Dev. in Group (second table), e.g.,
    years y1|(Intercept) = 0.59344.

  2. If I simulated my data with a pre-specified covariance structure (Sigma) of the y1, y2 and y3 intercepts and then sampled from a MVN(intercept_betas, Sigma), how do I compare my output to the pre-specified Sigma? Or rather, what represents my Sigma in the model fit output?

I will appreciate any help from you. Thanks

  • Operating System: 2.18.2
  • rstanarm Version: Ubuntu 18.10
2

Sigma in Estimates (first table – returned by summary()) is the estimates for the variance-covariance matrix of the group-specific terms. Note that the summary() table will show you both posterior medians and posterior means.

The Std.Dev. in Group (second table – returned by print()) is the estimates for the standard deviations and correlations of the groups-specific terms. From memory I think they are posterior medians, not posterior means, but perhaps @jonah or @bgoodri can confirm (it is the same as for stan_glmer).

So in other words… The first table is variances and covariances. The second table is standard deviations and correlations.

So for e.g. (0.59344) ^ 2 = 0.352

If the Sigma you used to simulate the data is a variance-covariance matrix then you can look at the posterior medians for Sigma in the summary() output. That is, compare your pre-specified Sigma to the posterior medians for Sigma in Estimates (first table).

p.s. I should really just write a vignette for stan_mvmer (although your question is more about labelling of the output and not actually the methodology!). But FWIW, the methodology is basically described in the longitudinal submodel section of the stan_jm vignette (if you just ignore everything else in that vignette about the survival models). Although, in your case you have two grouping factors – so you have two group-specific variance-covariance matrices, and each one is assumed independent of the other. This should hopefully be evident from the number of correlation parameters you saw returned in the second table of your output.

3

@sambrilleman Thank you so much for the detailed response. This is very clear and very helpful.

Vignette for stan_mvmer will really be helpful for people like us who switching to rstanarm.

Thanks,