Home  /  Products  /  Stata 17  /  Bayesian multilevel modeling
This page announced the new features in Stata 17. Please see our Stata 18 page for the new features in Stata 18.

Bayesian multilevel modeling: Nonlinear, joint, SEM-like...

Highlights

  • Outcomes: continuous, binary, ordinal, count, survival, ...
  • Convenient random-effects specification
  • Random intercepts and coefficients
  • Nested and crossed effects
  • Latent factors
  • Multiple levels of hierarchy
  • Linear and nonlinear models
  • Univariate and multivariate models
  • Joint, multivariate growth, and SEM-like models
  • Multivariate nonlinear random-effects models
  • Flexible random-effects prior distributions
  • Posterior distributions of random effects
  • MCMC diagnostics, including multiple chains
  • Full Bayesian postestimation features support

Multilevel models are used by many disciplines to model group-specific effects, which may arise at different levels of hierarchy. Think of regions, states nested within regions, and companies nested within states within regions. Or think hospitals, doctors nested within hospitals, and patients nested within doctors nested within hospitals.

In addition to the standard benefits of Bayesian analysis, Bayesian multilevel modeling offers many advantages for data with a small number of groups or panels. And it provides principled ways to compare effects across different groups by using posterior distributions of the effects. See more discussion here.

bayesmh has a new random-effects syntax that makes it easy to fit Bayesian multilevel models. And it opens the door to fitting new classes of multilevel models. You can fit univariate linear and nonlinear multilevel models more easily. And you can now fit multivariate linear and nonlinear multilevel models!

Think of mixed-effects nonlinear models as fit by menl, or some SEM models as fit by sem and gsem, or multivariate nonlinear models that contain random effects and, as of now, cannot be fit by any existing Stata command. You can now fit Bayesian counterparts of these models and more by using bayesmh.


Let's see it work

Two-level models

Let's start simple and fit a few two-level models first.

See Bayesian multilevel models using the bayes prefix. There, we show how to use bayes: mixed and other bayes multilevel commands to fit Bayesian multilevel models. Let's replicate some of the specifications here using the new random-effects syntax of bayesmh.

We consider the same data on math scores of pupils in the third and fifth years from different schools in Inner London (Mortimore et al. 1988). We want to investigate a school effect on math scores.

Let's fit a simple two-level random-intercept model to these data using bayesmh. We specify the random intercepts at the school level as U0[school] and include them in the regression specification.

bayesmh requires prior specifications for all parameters except random effects. For random effects, it assumes a normal prior distribution with mean 0 and variance component {var_U0}, where U0 is the name of the random effect. But we must specify the prior for {var_U0}. We specify normal priors with mean 0 and variance 10,000 for the regression coefficients and an inverse-gamma prior with shape and scale of 0.01 for variance components. We specify a seed for reproducibility and use a smaller MCMC size of 1,000.

. bayesmh math5 math3 U0[school], likelihood(normal({var_0}))
  prior({math5:}, normal(10000)) 
  prior({var_U0 var_0}, igamma(0.01, 0.01) split)
  rseed(17) mcmcsize(1000)

Burn-in 2500 aaaaaaaaa1000aaaaaaaaa2000aaaaa done
Simulation 1000 .........1000 done

Model summary
Likelihood:
math5 ~ normal(xb_math5,{var_0})
 
Priors:
{math5:math3 _cons} ~ normal(0,10000) (1)
{U0[school]} ~ normal(0,{var_U0}) (1)
{var_0} ~ igamma(0.01,0.01)
 
Hyperprior:
{var_U0} ~ igamma(0.01,0.01)
(1) Parameters are elements of the linear form xb_math5. Bayesian normal regression MCMC iterations = 3,500 Random-walk Metropolis–Hastings sampling Burn-in = 2,500 MCMC sample size = 1,000 Number of obs = 887 Acceptance rate = .184 Efficiency: min = .01949 avg = .02627 Log marginal-likelihood max = .03782
Equal-tailed
Mean Std. dev. MCSE Median [95% cred. interval]
math5
math3 .6106683 .0294298 .004785 .6118322 .557625 .6666446
_cons 30.33802 .286075 .054654 30.30869 29.78776 30.90507
var_0 28.04414 1.127223 .255309 28.12474 25.9292 30.12599
var_U0 4.18167 1.324194 .293471 3.791668 2.443605 8.09385

The results are similar to those from bayes: mixed in Random intercepts. We could replicate the same postestimation analysis from that section after bayesmh, including the display and graphs of random effects. In addition to the main model parameters, Bayesian models also estimate all random effects. This is in contrast with the frequentist analysis, where random effects are not estimated jointly with model parameters but are predicted after estimation conditional on model parameters.

Next we include random coefficients for math as c.math3#U1[school]. We additionally need to specify a prior for the variance component {var_U1}. We added to the variance-components list using the inverse-gamma prior.

. bayesmh math5 math3 U0[school] c.math3#U1[school], likelihood(normal({var_0}))
  prior({math5:}, normal(10000))
  prior({var_U0 var_U1 var_0}, igamma(0.01, 0.01) split)
  rseed(17) mcmcsize(1000)

Burn-in 2500 aaaaaaaaa1000aaaaa....2000..... done
Simulation 1000 .........1000 done

Model summary
Likelihood:
math5 ~ normal(xb_math5,{var_0})
 
Priors:
{math5:math3 _cons} ~ normal(0,10000) (1)
{U0[school]} ~ normal(0,{var_U0}) (1)
{U1[school]} ~ normal(0,{var_U1}) (1)
{var_0} ~ igamma(0.01,0.01)
 
Hyperprior:
{var_U0 var_U1} ~ igamma(0.01,0.01)
(1) Parameters are elements of the linear form xb_math5. Bayesian normal regression MCMC iterations = 3,500 Random-walk Metropolis–Hastings sampling Burn-in = 2,500 MCMC sample size = 1,000 Number of obs = 887 Acceptance rate = .2641 Efficiency: min = .009673 avg = .02501 Log marginal-likelihood max = .04479
Equal-tailed
Mean Std. dev. MCSE Median [95% cred. interval]
math5
math3 .6076992 .0399422 .005968 .6027789 .5351981 .6909692
_cons 30.30091 .2693482 .060468 30.30717 29.77415 30.81859
var_0 27.1167 1.143676 .1871 27.13383 24.40773 28.99094
var_U0 3.644527 .3810517 .104254 3.632184 2.865729 4.504112
var_U1 .0359823 .0132122 .004248 .0369494 .0121753 .0612364

Our results are similar to those from bayes: mixed in Random coefficients.

By default, bayesmh assumes that random effects U0[id] and U1[id] are independent a priori. But we can modify this by specifying a multivariate normal distribution for them with an unstructured covariance matrix {Sigma,m}. We additionally specify an inverse Wishart prior for the covariance {Sigma,m}.

. bayesmh math5 math3 U0[school] c.math3#U1[school], likelihood(normal({var_0}))
  prior({math5:}, normal(10000))
  prior({var_0}, igamma(0.01, 0.01))
  prior({U0[school] U1[school]}, mvn(2,0,0,{Sigma,m}))
  prior({Sigma,m}, iwishart(2,3,I(2)))
  rseed(17) mcmcsize(1000)

Burn-in 2500 aaaaaaaaa1000aaaaaaaaa2000aaaaa done
Simulation 1000 .........1000 done

Model summary
Likelihood:
math5 ~ normal(xb_math5,{var_0})
 
Priors:
{math5:math3 _cons} ~ normal(0,10000) (1)
{var_0} ~ igamma(0.01,0.01)
{U0[school] U1[school]} ~ mvnormal(2,0,0,{Sigma,m}) (1)
 
Hyperprior:
{Sigma,m} ~ iwishart(2,3,I(2))
(1) Parameters are elements of the linear form xb_math5. Bayesian normal regression MCMC iterations = 3,500 Random-walk Metropolis–Hastings sampling Burn-in = 2,500 MCMC sample size = 1,000 Number of obs = 887 Acceptance rate = .2681 Efficiency: min = .02805 avg = .03997 Log marginal-likelihood max = .05502
Equal-tailed
Mean Std. dev. MCSE Median [95% cred. interval]
math5
math3 .6358644 .055642 .010505 .6413544 .5284978 .7378708
_cons | 30.19804 .2872908 .049467 30.21301 29.62466 30.79241
var_0 26.47047 1.316738 .20655 26.47233 23.97093 28.83548
Sigma_1_1 4.355775 1.134325 .173332 4.251965 2.460442 6.865151
Sigma_2_1 -.337147 .1266224 .01707 -.3318653 -.6110905 -.1305513
Sigma_2_2 .0989839 .0211883 .003369 .0971182 .0633831 .1454557

The results are again similar to those from bayes: mixed, covariance(unstructured) in Random coefficients. Just like in that section, we could use bayesstats ic after bayesmh to compare the unstructured and independent models.

We can also use the new mvnexchangeable() prior option to specify an exchangeable random-effects covariance structure instead of an unstructured. An exchangeable covariance is characterized by two parameters, a common variance and a correlation. We specify additional priors for those parameters.

. bayesmh math5 math3 U0[school] c.math3#U1[school], likelihood(normal({var_0}))
  prior({math5:}, normal(10000))
  prior({var_0}, igamma(0.01, 0.01))
  prior({U0[school] U1[school]}, mvnexchangeable(2,0,0,{var_U},{rho}))
  prior({var_U}, igamma(0.01, 0.01)) prior({rho}, uniform(-1,1))
  rseed(17) mcmcsize(1000)

Burn-in 2500 aaaaaaaaa1000aaaaa....2000..... done
Simulation 1000 .........1000 done

Model summary
Likelihood:
math5 ~ normal(xb_math5,{var_0})
 
Priors:
{math5:math3 _cons} ~ normal(0,10000) (1)
{var_0} ~ igamma(0.01,0.01)
{U0[school] U1[school]} ~ mvnexchangeable(2,0,0,{var_U},{rho}) (1)
 
Hyperpriors:
{var_U} ~ igamma(0.01,0.01)
{rho} ~ uniform(-1,1)
(1) Parameters are elements of the linear form xb_math5. Bayesian normal regression MCMC iterations = 3,500 Random-walk Metropolis–Hastings sampling Burn-in = 2,500 MCMC sample size = 1,000 Number of obs = 887 Acceptance rate = .2192 Efficiency: min = .004793 avg = .009362 Log marginal-likelihood max = .01871
Equal-tailed
Mean Std. dev. MCSE Median [95% cred. interval]
math5
math3 .5981686 .0363804 .010405 .5997212 .525986 .6658922
_cons 30.38414 .1865243 .043121 30.36968 30.02465 30.73264
var_0 32.47519 .509181 .219875 32.45254 31.75455 33.4238
var_U .0635754 .0293628 .013412 .0574445 .0241849 .1300754
rho -.6144082 .2107051 .088129 -.6555609 -.9454211 -.2390335

We could fit other models from Bayesian multilevel models using the bayes prefix by using bayesmh. For instance, we could fit the three-level survival model by using

. bayesmh time education njobs prestige i.female U[birthyear] UU[id<birthyear],
  likelihood(stexponential, failure(failure))
  prior({time:}, normal(0,10000)) prior({var_U var_UU}, igamma(0.01,0.01) split)

and the crossed-effects logistic model by using

. bayesmh attain_gt_6 sex U[sid] V[pid], likelihood(logit)
  prior({attain_gt_6:}, normal(0,10000)) prior({var_U var_V}, igamma(0.01,0.01))

But all these models can be easily fit already by the bayes prefix. In what follows, we demonstrate models that cannot be fit with bayes:.


Nonlinear multilevel models

You can fit Bayesian univariate nonlinear multilevel models using bayesmh. The bayesmh syntax is the same as the menl command that fits classical nonlinear mixed-effects models, except that bayesmh additionally supports crossed effects such as UV[id1#id2] and latent factors such as L[_n].

See Three-level nonlinear model in [BAYES] bayesmh.


SEM-type models

You can use bayesmh to fit some structural equation models (SEMs) and models related to them. SEMs analyze multiple variables and use so-called latent variables in their specifications. A latent variable is a pseudo variable that has a different, unobserved, value for each observation. With bayesmh, you specify latent variables as L[_n]. You can use different latent variables in different equations, you can share the same latent variables across equations, you can place constraints on coefficients of latent variables, and more.

For examples, see Latent growth model and Item response theory in [BAYES] bayesmh.


Joint models of longitudinal and survival data

You can use bayesmh to model multiple outcomes of different types. Joint models of longitudinal and survival outcomes are one such example. These models are popular in practice because of their three main applications:

1. Account for informative dropout in the analysis of longitudinal data.

2. Study effects of baseline covariates on longitudinal and survival outcomes.

3. Study effects of time-dependent covariates on the survival outcome.

Below, we demonstrate the first application, but the same concept can be applied to the other two.

We will use a version of the Positive and Negative Symptom Scale (PANSS) data from a clinical trial comparing different drug treatmeans for schizophrenia (Diggle 1998). The data contain PANSS scores (panss) from patients who received three treatments (treat): placebo, haloperidol (reference), and risperidone (novel therapy). PANSS scores are measurements for psychiatric disorder. We would like to study the effects of the treatments on PANSS scores over time (week).

A model considered for PANSS scores is a longitudinal linear model with the effects of treatments, time (in weeks), and their interaction and random intercepts U[id].

. bayesmh panss i.treat##i.week U[id], likelihood(normal({var}))

So how does the survival model come into play here?

Many subjects withdrew from this study—only about half completed the full measurement schedule. Many subjects dropped out because of "inadequate for response", which suggests that the dropout may be informative and cannot be simply ignored in the analysis.

A dropout process can be viewed as a survival process with an informative dropout (infdrop) as an event of interest and with time to dropout as a survival time. Because we have longitudinal data, there are multiple observations per subject. So the dropout time is split into multiple spells according to week and is thus represented by the beginning time (t0) and end time (t1). At the left time point t0, an observation (or a spell) is considered left-truncated. We will assume a Weibull survival model for the dropout process. The dropout is likely to be related to the treatment, so we include it as the predictor in the survival model.

. bayesmh t1 i.treat, likelihood(stweibull({lnp}) failure(infdrop) ltruncated(t0))

One way to account for informative dropout is to include a shared random effect between the longitudinal and survival models that would induce correlation between the longitudinal outcome and the dropout process (Henderson, Diggle, and Dobson 2000).

. bayesmh (panss i.treat##i.week U[id]@1, likelihood(normal({var})))
         (t1 i.treat U[id], likelihood(stweibull({lnp}) failure(infdrop) ltruncated(t0)))

We added random intercepts from the longitudinal model to the survival model. We also constrained the coefficient for U[id] in the first equation to 1. We did this to emphasize that only the coefficient for U[id] in the second equation will be estimated. bayesmh actually makes this assumption automatically by default.

To fit the above model, we need to specify prior distributions for model parameters. We have many parameters, so we may need to specify somewhat informative priors. Recall that Bayesian models, in addition to the main model parameters, also estimate all the random-effects parameters U[id].

If there is an effect of dropout on the PANSS scores, it would be reasonable to assume that it would be positive. So we specify an exponential prior with scale 1 for the coefficient of U[id] in the survival model. We specify normal priors with mean 0 and variance of 1,000 for the constant {panss:_cons} and Weibull parameter {lnp}. We assume normal priors with mean 0 and variance 100 for other regression coefficients. And we use inverse-gamma priors with shape and scale of 0.01 for the variance components.

To improve sampling efficiency, we use Gibbs sampling for variance components and perform blocking of other parameters. We also specify initial values for some parameters by using the initial() option.

. bayesmh (panss i.treat##i.week U[id]@1, likelihood(norm({var})) )
          (t1 i.treat U[id], likelihood(stweibull({lnp}), failure(dropinf) ltruncated(t0))),
  prior({panss:_cons} {lnp},     normal(0,10000))
  prior({panss:i.treat##i.week}, normal(0,100))
  prior({t1:i.treat _cons},      normal(0,100))  
  prior({t1:U},                  exponential(1))                               
  prior({var_U} {var},           igamma(.01, .01) split)
  block({panss:i.treat#i.week}) block({panss:i.week})
  block({t1:i.treat U _cons}, split)
  block({var_U} {var}, split gibbs)
  initial({t1:U} 1 {U[id]} 10 {panss:_cons} 100)
  mcmcsize(2500) rseed(17)
  
Burn-in 2500 aaaaaaaaa1000aaaaaaaaa2000aaaaa done
Simulation 2500 .........1000.........2000..... done

Model summary
Likelihood:
panss ~ normal(xb_panss,{var})
t1 ~ stweibull(xb_t1,{lnp})
 
Priors:
{panss:_cons} ~ normal(0,10000) (1)
{panss:i.treat i.week i.treat#i.week} ~ normal(0,100) (1)
{U[id]} ~ normal(0,{var_U}) (1)
{var} ~ igamma(.01,.01)
{t1:i.treat _cons} ~ normal(0,100) (2)
{t1:U} ~ exponential(1) (2)
{lnp} ~ normal(0,10000)
 
Hyperprior:
{var_U} ~ igamma(.01,.01)
(1) Parameters are elements of the linear form xb_panss. (2) Parameters are elements of the linear form xb_t1. Bayesian normal regression MCMC iterations = 5,000 Metropolis–Hastings and Gibbs sampling Burn-in = 2,500 MCMC sample size = 2,500 No. of subjects = 685 Number of obs = 685 No. of failures = 63 Time at risk =863.6239911317825 Acceptance rate = .4608 Efficiency: min = .003913 avg = .03232 Log marginal-likelihood max = .2742
Equal-tailed
Mean Std. dev. MCSE Median [95% cred. interval]
panss
treat
2 -.3822383 2.271158 .359186 -.6136486 -4.434188 4.837333
3 -2.523494 3.80129 1.21543 -2.476083 -9.917332 4.074579
week
1 -4.993821 1.896496 .305945 -5.012834 -8.444717 -2.05126
2 -6.936372 1.57161 .453628 -6.939513 -10.50464 -3.908946
4 -4.844521 1.251981 .166785 -4.877955 -7.435107 -2.411917
6 -10.18118 1.816353 .361756 -10.03946 -14.2258 -6.98241
8 -10.25389 1.943371 .215791 -10.24783 -14.31332 -6.847999
treat#week
2 1 6.310975 2.434069 .390195 6.23198 1.213006 10.90485
2 2 7.027649 1.762338 .388778 6.840791 4.187137 10.5907
2 4 5.048269 1.863907 .351182 4.95867 1.458491 8.918415
2 6 15.32668 3.174594 .558032 14.99079 9.634857 21.59519
2 8 15.06479 3.072411 .646168 15.33875 8.316151 20.73611
3 1 -4.369372 2.892201 .659758 -4.479573 -9.651484 1.705437
3 2 -3.592772 2.198812 .444487 -3.576276 -7.864927 .982366
3 4 -11.22279 2.857886 .70199 -11.46703 -16.1155 -4.78894
3 6 -6.514449 1.87237 .483044 -6.457851 -9.986309 -2.939939
3 8 -2.078064 2.111598 .334112 -2.20723 -6.045809 1.881032
U 1 0 0 1 1 1
_cons 92.73873 2.162198 .367931 92.93747 88.40111 96.73237
t1
U .0596402 .0100107 .0009 .0598603 .0399564 .0783648
treat
2 .7984438 .3316247 .043614 .8106603 .1467264 1.457444
3 -.5172767 .4007821 .070892 -.5204102 -1.312922 .2484082
_cons -4.179667 .3958759 .082368 -4.19354 -4.944542 -3.359592
var 160.0879 9.848987 .376194 159.7498 142.23 180.8697
lnp .4849039 .0896179 .019375 .4879265 .2611824 .6596941
var_U 289.2579 39.46373 1.72886 285.8929 222.4319 372.773

We will not focus on the interpretation of all the results from this model, but we will comment on the coefficient {t1:U} for the shared random intercept. It is estimated to be 0.06, and its 95% credible interval does not contain 0. This means there is a positive association between PANSS scores and dropout times: the higher the PANSS score, the higher the chance of dropout. So, simple longitudinal analysis would not be appropriate for these data.


Multivariate nonlinear growth models

Hierarchical linear and nonlinear growth models are popular in many disciplines, such as health science, education, social sciences, engineering, and biology. Multivariate linear and nonlinear growth models are particularly useful in biological sciences to study the growth of wildlife species, where the growth is described by multiple measurements that are often correlated. Such models often have many parameters and are difficult to fit using classical methods. Bayesian estimation provides a viable alternative.

The above models can be fit by bayesmh using multiple-equation specifications, which now support random effects and latent factors. The equations can be all linear, all nonlinear, or a mixture of the two types. There are various ways to model the dependence between multiple outcomes (equations). For instance, you can include the same random effects but with different coefficients in different equations, or you can include different random effects but correlate them through the prior distribution.

Let's follow Jones et al. (2005) who applied a Bayesian bivariate growth model to study the growth of black-fronted tern chicks. The growth was measured by wing length \(y_1\) and weight \(y_2\). A linear growth model is assumed for wing length, and a logistic growth model is assumed for weight.

$$ \left( \begin{array}{c} y_{1,ij} \\ y_{2,ij} \end{array} \right) = \left( \begin{array}{c} U_i + V_i\, time_{ij} \\ C_i/\{1+dC_i\exp(-B_i\, time_{ij})\} \end{array} \right) + \left( \begin{array}{c} \varepsilon_{1,ij} \\ \varepsilon_{2,ij} \end{array} \right) $$

where \(d\) is a fixed parameter, \((U_i, V_i, C_i, B_i) \sim N(\mu, \Sigma)\), and \((\varepsilon_1, \varepsilon_2) \sim N(0, \Sigma_0)\).

There are four random effects at the individual (chick) level: \(U\) and \(V\) are the intercept and growth rate for the wings. In the equation for \(y_2\), we have random effects \(B\) and \(C\), which represent the rate and maximum weight. The location parameter is assumed to take the form \(dC\), where \(d\) is a constant. \((U, V, C, B)\) follow a multivariate normal distribution with an unstructured covariance. The errors from the two equations follow a bivariate normal distribution.

We use a fictional dataset simulated based on the description in Jones et al. (2005). We fit a bivariate normal model with error covariance {Sigma0,m}. The four random effects are assigned a multivariate normal prior with the corresponding mean parameters and covariance {Sigma,m}. The prior means are assigned normal distributions with mean 0 and variance 100. The covariance matrices are assigned inverse Wishart priors. Parameter d is assigned an exponential prior with scale 1. We use Gibbs sampling for covariance matrices and block parameters to improve efficiency. We also specify initial values, use a smaller MCMC size of 2,500, and specify a random-number seed for reproducibility.

. bayesmh (y1 = ({U[id]} + time*{V[id]}))
         (y2 = ({C[id]}/(1+{d}*{C[id]}*exp(-{B[id]}*time)))),
  likelihood(mvnormal({Sigma0,m}))
  prior({U V C B}, mvnormal(4,{u},{v},{c},{b},{Sigma,m}))
  prior({u v c b},  normal(0, 100))
  prior({Sigma0,m}, iwishart(2,3,I(2)))
  prior({Sigma,m},  iwishart(4,5,I(4)))
  prior({d}, exp(1))
  block({d u v b c}, split)
  block({Sigma0,m} {Sigma,m}, gibbs split)
  init({U[id] u} -10 {V[id] v} 10 {C[id] c} 100 {d} 1)
  mcmcsize(2500) rseed(17)
 
Burn-in 2500 aaaaaaaaa1000aaaaaaaaa2000aaaaa done
Simulation 2500 .........1000.........2000..... done

Model summary
Likelihood:
y1 y2 ~ mvnormal(2,,,{Sigma0,m})
 
Priors:
{Sigma0,m} ~ iwishart(2,3,I(2))
{U[id] V[id] C[id] B[id]} ~ mvnormal(4,{u},{v},{c},{b},{Sigma,m})
{d} ~ exponential(1)
 
Hyperpriors:
{u v c b} ~ normal(0,100)
{Sigma,m} ~ iwishart(4,5,I(4))
 
Expressions:
expr1 : {U[id]} + time*{V[id]}
expr2 : {C[id]}/(1+{d}*{C[id]}*exp(-{B[id]}*time))
Bayesian multivariate normal regression MCMC iterations = 5,000 Metropolis—Hastings and Gibbs sampling Burn-in = 2,500 MCMC sample size = 2,500 Number of obs = 414 Acceptance rate = .4713 Efficiency: min = .01174 avg = .2265 Log marginal-likelihood max = .7028
Equal-tailed
Mean Std. dev. MCSE Median [95% cred. interval]
d .0634061 .0025888 .000478 .0635744 .0579154 .0680656
u -12.84796 3.011731 .255283 -12.97586 -18.25202 -6.451113
v 5.977761 .2446379 .023368 5.990374 5.422395 6.404792
c 78.42872 3.602142 .368536 78.7988 70.10973 84.34357
b .2208688 .0471093 .002637 .2229167 .1242395 .3148616
Sigma0_1_1 7.956314 .5825538 .017417 7.926544 6.871581 9.158582
Sigma0_2_1 2.625951 .6406367 .021819 2.632427 1.430312 3.875557
Sigma0_2_2 18.85203 1.342218 .038113 18.81303 16.36956 21.67296
Sigma_1_1 192.8405 67.11091 2.92639 179.5316 101.754 362.8648
Sigma_2_1 -8.029962 4.209058 .21859 -7.334189 -17.74035 -1.783869
Sigma_3_1 -108.4137 63.18093 3.39159 -97.77067 -258.3206 -18.55377
Sigma_4_1 .4582266 .6998019 .021466 .4405483 -.8234645 1.983518
Sigma_2_2 1.193545 .4200058 .025011 1.10642 .6352668 2.223882
Sigma_3_2 12.45667 5.664299 .404336 11.29209 5.259565 27.34906
Sigma_4_2 -.0023492 .0557342 .001842 -.0034794 -.1104773 .1078309
Sigma_3_3 234.2312 95.14968 6.93288 212.8518 117.8635 471.0824
Sigma_4_3 -.2949588 .829987 .032991 -.2727646 -2.063978 1.386505
Sigma_4_4 .0454308 .0136201 .000325 .0428103 .0257433 .0790052

Error covariances and random-effects covariance values are not 0, which suggests that the wing length and weight measurements are related.

We use bayesstats summary to compute correlation estimates.

. bayesstats summary (corr0: {Sigma0_1_2}/sqrt({Sigma0_1_1}*{Sigma0_2_2}))

Posterior summary statistics                      MCMC sample size =     2,500

       corr0 :  {Sigma0_1_2}/sqrt({Sigma0_1_1}*{Sigma0_2_2})

Equal-tailed
Mean Std. dev. MCSE Median [95% cred. interval]
corr0 .2141337 .048555 .001653 .2161776 .1162883 .3055137

There is a positive correlation, 0.21, between the error terms.

We also compute the correlation between the rate of wing growth and the maximum weight.

. bayesstats summary (corr24: {Sigma_2_3}/sqrt({Sigma_2_2}*{Sigma_3_3}))

Posterior summary statistics                      MCMC sample size =     2,500

      corr24 :  {Sigma_2_3}/sqrt({Sigma_2_2}*{Sigma_3_3})

Equal-tailed
Mean Std. dev. MCSE Median [95% cred. interval]
corr24 .7328452 .1065301 .005605 .7508832 .4838739 .8959725

The estimated correlation is 0.73, which is high. The wing length and weight measurements are clearly correlated and should be modeled jointly.


References

Diggle, P. J. 1998. Dealing with missing values in longitudinal studies. In Recent Advances in the Statistical Analysis of Medical Data, ed. B. S. Everitt and G. Dunn, 203–228. London: Arnold.

Henderson, R., P. Diggle, and A. Dobson. 2000. Joint modeling of longitudinal measurements and event time data. Biostatistics 4: 465–480.

Jones, G., R. J. Keedwell, A. D. L. Noble, and D. I. Hedderley. 2005. Dating chicks: Calibration and discrimination in a nonlinear multivariate hierarchical growth model. Journal of Agricultural, Biological, and Environmental Statistics 10: 306–320.

Mortimore, P., P. Sammons, L. Stoll, D. Lewis, and R. Ecob. 1988. School Matters: The Junior Years. Somerset, UK: Open Books.


Additional resources

Learn more in the Stata Bayesian Analysis Reference Manual.

See Bayesian longitudinal/panel-data models.

See Bayesian multilevel models using the bayes prefix.

The Stata Blog: Comparing transmissibility of Omicron lineages