Let’s see it work

• Two-level models
• Nonlinear multilevel models
• SEM-type models
• Joint models of longitudinal and survival data
• Multivariate nonlinear growth models

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 y1 and weight y2. A linear growth model is assumed for wing length, and a logistic growth model is assumed for weight.

where d

$d$

$d$ is a fixed parameter, (Ui,Vi,Ci,Bi)∼N(μ,Σ)

$(U_i,V_i,C_i,B_i)\sim N(\mu ,\mathrm{\Sigma })$

$(U_i, V_i, C_i, B_i) \sim N(\mu, \Sigma)$,
and (ε1,ε2)∼N(0,Σ0)

$({\epsilon }_1,{\epsilon }_2)\sim N(0,{\mathrm{\Sigma }}_0)$

$(\varepsilon_1, \varepsilon_2) \sim N(0, \Sigma_0)$.

There are four random effects at the individual (chick) level: U

$U$

$U$ and V

$V$

$V$
are the intercept and growth rate for the wings. In the equation for y2

$y_2$

$y_2$,
we have random effects B

$B$

$B$ and C

$C$

$C$, which represent the rate and maximum
weight. The location parameter is assumed to take the form dC

$dC$

$dC$, where d

$d$

$d$
is a constant. (U,V,C,B)

$(U,V,C,B)$

$(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[id] V[id] C[id] B[id]}, 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)
 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

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.