Modélisation bayésienne multi-niveaux

Nouveau dans

Modélisation bayésienne à plusieurs niveaux : Non-linéaire, conjointe, de type SEM…

Les modèles multiniveaux sont utilisés par de nombreuses disciplines pour modéliser les effets spécifiques à un groupe, qui peuvent apparaître à différents niveaux de la hiérarchie. Pensez aux régions, aux États nichés dans les régions et aux entreprises nichées dans les États au sein des régions. Ou pensez aux hôpitaux, aux médecins imbriqués dans les hôpitaux et aux patients imbriqués dans les médecins imbriqués dans les hôpitaux.
Outre les avantages habituels de l’analyse bayésienne, la modélisation bayésienne à plusieurs niveaux offre de nombreux avantages pour les données comportant un petit nombre de groupes ou de panels. De plus, elle fournit des moyens fondés sur des principes pour comparer les effets entre différents groupes en utilisant les distributions postérieures des effets. Voir plus de discussion ici.
bayesmh possède une nouvelle syntaxe à effets aléatoires qui facilite l’ajustement des modèles bayésiens à plusieurs niveaux. Et cela ouvre la porte à l’ajustement de nouvelles classes de modèles multiniveaux. Vous pouvez ajuster plus facilement des modèles multiniveaux linéaires et non linéaires univariés. Et vous pouvez maintenant ajuster des modèles multiniveaux linéaires et non linéaires multivariés !
Pensez aux modèles non linéaires à effets mixtes tels qu’ajustés par menl, ou à certains modèles SEM tels qu’ajustés par sem et gsem, ou encore aux modèles non linéaires multivariés qui contiennent des effets aléatoires et qui, pour l’instant, ne peuvent être ajustés par aucune commande Stata existante. Vous pouvez maintenant ajuster les contreparties bayésiennes de ces modèles et plus encore en utilisant bayesmh.

Points forts

Résultats : continu, binaire, ordinal, comptage, survie, …
Spécification pratique des effets aléatoires
Effets imbriqués et croisés
Facteurs latents
Niveaux multiples de hiérarchie
Modèles linéaires et non linéaires
Modèles univariés et multivariés
Modèles conjoints, modèles de croissance multivariés et modèles de type SEM
Modèles à effets aléatoires non linéaires multivariés
Distributions préalables flexibles des effets aléatoires
Distributions postérieures des effets aléatoires
Diagnostics MCMC, y compris les chaînes multiples
Support complet des fonctions de post-estimation bayésienne

Voyons comment cela fonctionne

Mdèles à deux niveaux
Modèles multiniveaux non linéaires
Modèles de type SEM
Modèles conjoints de données longitudinales et de survie
Modèles de croissance non linéaires multivariés

Modèles à deux niveaux

Commençons simplement et ajustons d’abord quelques modèles à deux niveaux.
Voir Modèles bayésiens à plusieurs niveaux utilisant le préfixe bayes. Nous y montrons comment utiliser bayes : mixed et d’autres commandes bayes multilevel pour ajuster des modèles bayésiens multilevel. Répétons certaines des spécifications ici en utilisant la nouvelle syntaxe à effets aléatoires de bayesmh.
Nous considérons les mêmes données sur les résultats en mathématiques des élèves de troisième et cinquième années de différentes écoles du centre de Londres (Mortimore et al. 1988). Nous voulons étudier l’effet de l’école sur les résultats en mathématiques.
Ajustons un modèle simple à deux niveaux d’intercepts aléatoires à ces données en utilisant bayesmh. Nous spécifions les intercepts aléatoires au niveau de l’école comme U0[école] et les incluons dans la spécification de la régression.
bayesmh requiert des spécifications préalables pour tous les paramètres, sauf pour les effets aléatoires. Pour les effets aléatoires, il suppose une distribution antérieure normale avec une moyenne de 0 et une composante de variance {var_U0}, où U0 est le nom de l’effet aléatoire. Mais nous devons spécifier la priorité pour {var_U0}. Nous spécifions des antériorités normales de moyenne 0 et de variance 10 000 pour les coefficients de régression et une antériorité de type gamma inverse avec une forme et une échelle de 0,01 pour les composantes de variance. Nous spécifions une graine pour la reproductibilité et utilisons une taille MCMC plus petite de 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

Les résultats sont similaires à ceux de bayes : mixed dans Random intercepts. Nous avons pu reproduire la même analyse de post-estimation de cette section après bayesmh, y compris l’affichage et les graphiques des effets aléatoires. En plus des paramètres principaux du modèle, les modèles bayésiens estiment également tous les effets aléatoires. Ceci est en contraste avec l’analyse fréquentiste, où les effets aléatoires ne sont pas estimés conjointement avec les paramètres du modèle mais sont prédits après l’estimation conditionnellement aux paramètres du modèle.
Ensuite, nous incluons des coefficients aléatoires pour les mathématiques comme c.math3#U1[school]. Nous devons en outre spécifier une antériorité pour la composante de variance {var_U1}. Nous avons ajouté à la liste des composantes de variance en utilisant la priorisation inverse-gamma.

. 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

Nos résultats sont similaires à ceux de bayes : mixed in Random coefficients.
Par défaut, bayesmh suppose que les effets aléatoires U0[id] et U1[id] sont indépendants a priori. Mais nous pouvons modifier cela en spécifiant une distribution normale multivariée pour eux avec une matrice de covariance non structurée {Sigma,m}. Nous spécifions en outre une priorité Wishart inverse pour la 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

Les résultats sont à nouveau similaires à ceux de bayes : mixed, covariance(unstructured) dans Random coefficients. Tout comme dans cette section, nous pourrions utiliser bayesstats ic après bayesmh pour comparer les modèles non structurés et indépendants.
Nous pouvons également utiliser la nouvelle option préalable mvnexchangeable() pour spécifier une structure de covariance à effets aléatoires échangeable au lieu d’une structure non structurée. Une covariance échangeable est caractérisée par deux paramètres, une variance commune et une corrélation. Nous spécifions des prieurs supplémentaires pour ces paramètres.

. 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

Nous pouvons ajuster d’autres modèles à partir de modèles bayésiens à plusieurs niveaux en utilisant le préfixe bayes en utilisant bayesmh. Par exemple, nous pouvons ajuster le modèle de survie à trois niveaux en utilisant

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

et le modèle logistique à effets croisés en utilisant

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

Mais tous ces modèles peuvent déjà être facilement ajustés par le préfixe de Bayes. Dans ce qui suit, nous démontrons des modèles qui ne peuvent pas être ajustés avec bayes :.

Modèles multiniveaux non linéaires

Vous pouvez ajuster des modèles multiniveaux bayésiens univariés non linéaires à l’aide de bayesmh. La syntaxe de bayesmh est la même que celle de la commande menl qui ajuste les modèles non linéaires classiques à effets mixtes, sauf que bayesmh supporte en plus les effets croisés tels que UV[id1#id2] et les facteurs latents tels que L[_n].
Voir le modèle non linéaire à trois niveaux dans [BAYES] bayesmh.

Modèles de type SEM

Vous pouvez utiliser bayesmh pour ajuster certains modèles d’équations structurelles (SEM) et les modèles qui leur sont liés. Les SEM analysent plusieurs variables et utilisent des variables dites latentes dans leurs spécifications. Une variable latente est une pseudo-variable qui a une valeur différente, non observée, pour chaque observation. Avec bayesmh, vous spécifiez les variables latentes comme L[_n]. Vous pouvez utiliser différentes variables latentes dans différentes équations, vous pouvez partager les mêmes variables latentes entre plusieurs équations, vous pouvez placer des contraintes sur les coefficients des variables latentes, et plus encore.
Pour des exemples, voir Latent growth model et Item response theory dans [BAYES] bayesmh.

Modèles conjoints de données longitudinales et de survie

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:

Tenir compte des abandons informatifs dans l’analyse des données longitudinales.
Étudier les effets des covariables de base sur les résultats longitudinaux et de survie.
Étudier les effets des covariables dépendantes du temps sur les résultats de survie.

Nous présentons ci-dessous la première application, mais le même concept peut être appliqué aux deux autres.
Nous utiliserons une version de l’échelle des symptômes positifs et négatifs (PANSS) provenant d’un essai clinique comparant différents traitements médicamenteux de la schizophrénie (Diggle 1998). Les données contiennent les scores PANSS (panss) de patients qui ont reçu trois traitements (treat) : placebo, haloperidol (référence), et risperidone (nouvelle thérapie). Les scores PANSS sont des mesures des troubles psychiatriques. Nous aimerions étudier les effets des traitements sur les scores PANSS au fil du temps (semaine).
Le modèle considéré pour les scores PANSS est un modèle linéaire longitudinal avec les effets des traitements, du temps (en semaines) et de leur interaction et des intercepts aléatoires U[id].

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

Alors comment le modèle de survie entre-t-il en jeu ici ?
De nombreux sujets se sont retirés de l’étude – seule la moitié environ a terminé le programme de mesure complet. De nombreux sujets ont abandonné pour cause de « réponse inadéquate », ce qui suggère que l’abandon peut être informatif et ne peut être simplement ignoré dans l’analyse.
Un processus d’abandon peut être considéré comme un processus de survie avec un abandon informatif (infdrop) comme événement d’intérêt et avec le temps jusqu’à l’abandon comme temps de survie. Comme nous avons des données longitudinales, il y a plusieurs observations par sujet. Ainsi, le temps d’abandon est divisé en plusieurs périodes selon la semaine et est donc représenté par le temps de début (t0) et le temps de fin (t1). Au point de temps gauche t0, une observation (ou une période) est considérée comme tronquée à gauche. Nous supposerons un modèle de survie de Weibull pour le processus d’abandon. L’abandon est susceptible d’être lié au traitement, nous l’incluons donc comme prédicteur dans le modèle de survie.

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

Une façon de tenir compte de l’abandon informatif est d’inclure un effet aléatoire partagé entre les modèles longitudinal et de survie qui induirait une corrélation entre le résultat longitudinal et le processus d’abandon (Henderson, Diggle et 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)))

Nous avons ajouté les intercepts aléatoires du modèle longitudinal au modèle de survie. Nous avons également contraint le coefficient de U[id] dans la première équation à 1. Nous avons fait cela pour souligner que seul le coefficient de U[id] dans la deuxième équation sera estimé. bayesmh fait en fait cette hypothèse automatiquement par défaut.
Pour ajuster le modèle ci-dessus, nous devons spécifier des distributions préalables pour les paramètres du modèle. Nous avons beaucoup de paramètres, donc nous pouvons avoir besoin de spécifier des prieurs quelque peu informatifs. Rappelons que les modèles bayésiens, en plus des paramètres du modèle principal, estiment également tous les paramètres d’effets aléatoires U[id].
S’il existe un effet de l’abandon sur les scores PANSS, il serait raisonnable de supposer qu’il soit positif. Nous spécifions donc une priorité exponentielle avec une échelle de 1 pour le coefficient de U[id] dans le modèle de survie. Nous spécifions des prieurs normaux avec une moyenne de 0 et une variance de 1 000 pour la constante {panss:_cons} et le paramètre de Weibull {lnp}. Nous spécifions des prieurs normaux avec une moyenne de 0 et une variance de 100 pour les autres coefficients de régression. Et nous utilisons des prieurs de type gamma inverse avec une forme et une échelle de 0,01 pour les composantes de la variance.
Pour améliorer l’efficacité de l’échantillonnage, nous utilisons l’échantillonnage de Gibbs pour les composantes de la variance et effectuons le blocage des autres paramètres. Nous spécifions également des valeurs initiales pour certains paramètres en utilisant l’option initial().

. 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

Nous ne nous attarderons pas sur l’interprétation de tous les résultats de ce modèle, mais nous commenterons le coefficient {t1:U} pour l’intercept aléatoire partagé. Il est estimé à 0,06, et son intervalle crédible à 95% ne contient pas 0. Cela signifie qu’il existe une association positive entre les scores PANSS et les délais d’abandon : plus le score PANSS est élevé, plus le risque d’abandon est grand. Ainsi, une simple analyse longitudinale ne serait pas appropriée pour ces données.

Modèles de croissance non linéaires multivariés

Les modèles de croissance linéaires et non linéaires hiérarchisés sont populaires dans de nombreuses disciplines, telles que les sciences de la santé, l’éducation, les sciences sociales, l’ingénierie et la biologie. Les modèles de croissance linéaires et non linéaires multivariés sont particulièrement utiles dans les sciences biologiques pour étudier la croissance des espèces sauvages, où la croissance est décrite par de multiples mesures qui sont souvent corrélées. Ces modèles ont souvent de nombreux paramètres et sont difficiles à ajuster à l’aide des méthodes classiques. L’estimation bayésienne offre une alternative viable.
Les modèles ci-dessus peuvent être ajustés par bayesmh en utilisant des spécifications à équations multiples, qui supportent maintenant les effets aléatoires et les facteurs latents. Les équations peuvent être toutes linéaires, toutes non linéaires, ou un mélange des deux types. Il existe plusieurs façons de modéliser la dépendance entre plusieurs résultats (équations). Par exemple, vous pouvez inclure les mêmes effets aléatoires mais avec des coefficients différents dans différentes équations, ou vous pouvez inclure différents effets aléatoires mais les corréler par le biais de la distribution antérieure.
Suivons Jones et al. (2005) qui ont appliqué un modèle de croissance bayésien à deux variables pour étudier la croissance des poussins de sterne à front noir. La croissance a été mesurée par la longueur des ailes y1 et le poids y2. Un modèle de croissance linéaire est supposé pour la longueur de l’aile, et un modèle de croissance logistique est supposé pour le poids.

(y 1, i j y 2, i j) = (U i + V i t i m e i j C i / {1 + d C i exp (- B i t i m e i j)}) + (ε 1, i j ε 2, i j)

where $d$

d

$d$ is a fixed parameter, $(U_{i}, V_{i}, C_{i}, B_{i}) \sim N (μ, Σ)$

(U_{i}, V_{i}, C_{i}, B_{i}) \sim N (μ, Σ)

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

(ε_{1}, ε_{2}) \sim N (0, Σ_{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 $y_{2}$

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 $d C$

d C

$dC$ , where $d$

d

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

(U, V, C, B)

$(U, V, C, B)$ suivent une distribution normale multivariée avec
une covariance non structurée. Les erreurs des deux équations suivent une
distribution normale bivariée.

Nous utilisons un ensemble de données fictives simulées sur la base de la description de Jones et al. (2005). Nous ajustons un modèle normal bivarié avec une covariance d’erreur {Sigma0,m}. Les quatre effets aléatoires se voient attribuer un antécédent normal multivarié avec les paramètres moyens correspondants et la covariance {Sigma,m}.< On attribue aux moyennes antérieures des distributions normales de moyenne 0 et de variance 100. Les matrices de covariance sont affectées de prieurs Wishart inverses. Le paramètre d est affecté d'une priorité exponentielle d'échelle 1. Nous utilisons l'échantillonnage de Gibbs pour les matrices de covariance et les paramètres de bloc pour améliorer l'efficacité. Nous spécifions également les valeurs initiales, utilisons une taille MCMC plus petite de 2 500, et spécifions une graine de nombre aléatoire pour la reproductibilité.

. 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

Il existe une corrélation positive, 0,21, entre les termes d’erreur.
Nous calculons également la corrélation entre le taux de croissance de l’aile et le poids maximal.

. 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

La corrélation estimée est de 0,73, ce qui est élevé. Les mesures de longueur d’aile et de poids sont clairement corrélées et doivent être modélisées conjointement.

Références

Diggle, P. J. 1998. Dealing with missing values in longitudinal studies « . Dans Recent Advances in the Statistical Analysis of Medical Data, édité par B. S. Everitt et G. Dunn, 203-228. Londres : Arnold.
Henderson, R., P. Diggle et 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, et 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, et R. Ecob. 1988. School Matters : The Junior Years. Somerset, Royaume-Uni : Open Books.

Ressources supplémentaires

Pour en savoir plus, consultez le manuel de référence de l’analyse bayésienne de Stata.
Voir les modèles bayésiens de données longitudinales/panel.
Voir Modèles bayésiens multiniveaux utilisant le préfixe bayes.

VOIR TOUTES LES NOUVELLES FONCTIONNALITÉS

COMMANDEZ STATA 17

UPGRADE