»  Home »  Products »  Stata 17 »  Bayesian multilevel modeling

# 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.

## 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(dropinf) 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(dropinf) 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.

See Bayesian longitudinal/panel-data models.

See Bayesian multilevel models using the bayes prefix.

The Stata Blog: Comparing transmissibility of Omicron lineages