Podívejme se, jak to funguje

• Dvouúrovňové modely
• Nelineární víceúrovňové modely
• Modely typu SEM
• Společné modely longitudinálních dat a dat o přežití
• Vícerozměrné nelineární modely růstu

Dvouúrovňové modely

Začněme jednoduše a nejprve sestavme několik dvouúrovňových modelů.

Viz Bayesovské víceúrovňové modely s použitím předpony Bayes. Tam si ukážeme, jak pomocí příkazu bayes: mixed a dalších příkazů bayes multilevel fitovat bayesovské víceúrovňové modely. Zopakujme si zde některé specifikace s použitím nové syntaxe bayesmh s náhodnými efekty.

Uvažujeme stejná data o výsledcích žáků třetích a pátých ročníků v matematice z různých škol ve vnitřním Londýně (Mortimore et al. 1988). Chceme zkoumat vliv školy na výsledky z matematiky.

Napasujme na tato data jednoduchý dvouúrovňový model s náhodným průmětem pomocí bayesmh. Náhodné intercepty na úrovni školy specifikujeme jako U0[school] a zahrneme je do regresní specifikace.

bayesmh vyžaduje předchozí specifikace pro všechny parametry kromě náhodných efektů. Pro náhodné efekty předpokládá normální předběžné rozdělení se střední hodnotou 0 a složkou rozptylu {var_U0}, kde U0 je název náhodného efektu. Musíme však specifikovat prioritu pro {var_U0}. Pro regresní koeficienty zadáme normální prioritu se střední hodnotou 0 a rozptylem 10 000 a pro složky rozptylu zadáme prioritu inverzní gama s tvarem a škálou 0,01. Pro reprodukovatelnost zadáváme seed a používáme menší velikost MCMC 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

Výsledky jsou podobné výsledkům z bayesova modelu: smíšené v Náhodné zásahy. Mohli jsme zopakovat stejnou postestimulační analýzu z této části po bayesmh, včetně zobrazení a grafů náhodných efektů. Kromě hlavních parametrů modelu odhadují bayesovské modely také všechny náhodné efekty. To je rozdíl oproti frekvenční analýze, kde se náhodné efekty neodhadují společně s parametry modelu, ale předpovídají se po odhadu podmíněném parametry modelu.

Dále zahrneme náhodné koeficienty pro math jako c.math3#U1[school]. Dále musíme určit prioritu pro složku rozptylu {var_U1}. Do seznamu komponent rozptylu jsme přidali prioritu inverzní gama.

. 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

Naše výsledky jsou podobné výsledkům z bayes: mixed analýzy: smíšené v Náhodných koeficientech.

Ve výchozím nastavení bayesmh předpokládá, že náhodné efekty U0[id] a U1[id] jsou a priori nezávislé. Můžeme to však upravit tak, že pro ně zadáme vícerozměrné normální rozdělení s nestrukturovanou kovarianční maticí {Sigma,m}. Pro kovarianci {Sigma,m} navíc určíme inverzní Wishartovu prioritu.

. 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

Výsledky jsou opět podobné výsledkům z bayes: mixed, covariance(unstructured) v náhodných koeficientech. Stejně jako v této části bychom mohli použít bayesstats ic po bayesmh k porovnání nestrukturovaného a nezávislého modelu.

Můžeme také použít novou prior volbu mvnexchangeable() a zadat výměnnou kovarianční strukturu náhodných efektů místo nestrukturované. Výměnná kovariance je charakterizována dvěma parametry, společným rozptylem a korelací. Pro tyto parametry zadáváme další priory.

. 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

Další modely z bayesovských víceúrovňových modelů jsme mohli fitovat pomocí předpony bayes pomocí bayesmh. Například bychom mohli fitovat tříúrovňový model přežití pomocí přípony

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

a logistický model se zkříženými efekty s použitím

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

Všechny tyto modely však lze snadno napasovat již pomocí bayes prefixu. V následujícím textu si ukážeme modely, které nelze fitovat pomocí bayesovské předpony:.

Nelineární víceúrovňové modely

Bayesovské jednorozměrné nelineární víceúrovňové modely můžete fitovat pomocí bayesmh. Syntaxe bayesmh je stejná jako u příkazu menl, který fituje klasické nelineární modely se smíšenými efekty, s tím rozdílem, že bayesmh navíc podporuje zkřížené efekty, jako je UV[id1#id2], a latentní faktory, jako je L[_n].

Viz Tříúrovňový nelineární model v [BAYES] bayesmh.

Modely typu SEM

Pomocí bayesmh můžete fitovat některé modely strukturálních rovnic (SEM) a modely s nimi související. SEM analyzují více proměnných a ve svých specifikacích používají tzv. latentní proměnné. Latentní proměnná je pseudoproměnná, která má pro každé pozorování jinou, nepozorovanou hodnotu. V programu bayesmh se latentní proměnné zadávají jako L[_n]. V různých rovnicích můžete používat různé latentní proměnné, stejné latentní proměnné můžete sdílet napříč rovnicemi, na koeficienty latentních proměnných můžete klást omezení a další.

Příklady viz Latentní růstový model a Teorie odezvy na položku v [BAYES] bayesmh.

Společné modely longitudinálních dat a dat o přežití

Pomocí bayesmh můžete modelovat více výsledků různých typů. Jedním z takových příkladů jsou společné modely longitudinálních výsledků a výsledků přežití. Tyto modely jsou v praxi oblíbené díky svým třem hlavním aplikacím:

1. Zohlednění informativního výpadku při analýze longitudinálních dat.
2. Studie vlivu výchozích kovariát na longitudinální výsledky a výsledky přežití.
3. Studium vlivu kovariát závislých na čase na výsledek přežití.

Níže si ukážeme první aplikaci, ale stejný koncept lze použít i na další dvě.

Použijeme verzi škály pozitivních a negativních příznaků (PANSS), která pochází z klinické studie porovnávající různé způsoby léčby schizofrenie (Diggle 1998). Data obsahují skóre PANSS (panss) od pacientů, kteří dostávali tři léčebné postupy (treat): placebo, haloperidol (referenční) a risperidon (nová léčba). Skóre PANSS je měřením psychické poruchy. Chtěli bychom studovat účinky léčby na skóre PANSS v průběhu času (week).

Uvažovaný model pro skóre PANSS je longitudinální lineární model s účinky léčby, času (v týdnech) a jejich interakcí a náhodnými intercepty U[id].

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

Jak se zde uplatňuje model přežití?

Mnoho subjektů z této studie odstoupilo – jen asi polovina dokončila celý plán měření. Mnoho subjektů odstoupilo z důvodu „nedostatečné odpovědi“, což naznačuje, že odstoupení může být informativní a nelze je v analýze jednoduše ignorovat.

Na proces vyřazení lze nahlížet jako na proces přežívání, přičemž informativní vyřazení (infdrop) je událost, která nás zajímá, a doba do vyřazení je doba přežívání. Protože máme longitudinální data, na jeden subjekt připadá více pozorování. Doba vysazení je tedy rozdělena na více období podle týdnů a je tedy reprezentována časem začátku (t0) a časem konce (t1). V levém časovém bodě t0 je pozorování (nebo kouzlo) považováno za levostranné. Pro proces vyřazení budeme předpokládat Weibullův model přežití. Vysazení pravděpodobně souvisí s léčbou, proto jej zahrneme jako prediktor do modelu přežití.

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

Jedním ze způsobů, jak zohlednit informativní předčasné ukončení studia, je zahrnout společný náhodný efekt mezi longitudinální model a model přežití, který by vyvolal korelaci mezi longitudinálním výsledkem a procesem předčasného ukončení studia (Henderson, Diggle a 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)))

Do modelu přežití jsme přidali náhodné intercepty z longitudinálního modelu. Omezili jsme také koeficient U[id] v první rovnici na hodnotu 1. Udělali jsme to proto, abychom zdůraznili, že se bude odhadovat pouze koeficient pro U[id] ve druhé rovnici. bayesmh ve skutečnosti tento předpoklad automaticky provádí ve výchozím nastavení.

Abychom mohli výše uvedený model přizpůsobit, musíme zadat prioritní rozdělení parametrů modelu. Máme mnoho parametrů, takže možná budeme muset zadat poněkud informativní priory. Připomeňme, že bayesovské modely kromě hlavních parametrů modelu odhadují také všechny parametry náhodných efektů U[id].

Pokud existuje vliv předčasného ukončení studia na skóre PANSS, bylo by rozumné předpokládat, že bude pozitivní. Proto v modelu přežití určíme exponenciální prioritu se škálou 1 pro koeficient U[id]. Pro konstantu {panss:_cons} a Weibullův parametr {lnp} zadáme normální priory se střední hodnotou 0 a rozptylem 1 000. Pro ostatní regresní koeficienty předpokládáme normální priory se střední hodnotou 0 a rozptylem 100. A pro složky rozptylu použijeme inverzní gama priory s tvarem a škálou 0,01.

Pro zvýšení efektivity vzorkování použijeme Gibbsovo vzorkování pro složky rozptylu a provedeme blokování ostatních parametrů. Pro některé parametry také určujeme počáteční hodnoty pomocí volby 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

Nebudeme se věnovat interpretaci všech výsledků tohoto modelu, ale okomentujeme koeficient {t1:U} pro společný náhodný intercept. Jeho hodnota je odhadnuta na 0,06 a jeho 95% věrohodný interval neobsahuje 0. To znamená, že existuje pozitivní souvislost mezi skóre PANSS a dobou vysazení: čím vyšší je skóre PANSS, tím vyšší je pravděpodobnost vysazení. Jednoduchá longitudinální analýza by tedy pro tyto údaje nebyla vhodná.

Vícerozměrné nelineární modely růstu

Hierarchické lineární a nelineární modely růstu jsou oblíbené v mnoha oborech, například ve zdravotnictví, vzdělávání, sociálních vědách, inženýrství a biologii. Vícerozměrné lineární a nelineární modely růstu jsou zvláště užitečné v biologických vědách při studiu růstu volně žijících druhů, kde je růst popsán více měřeními, která jsou často korelovaná. Takové modely mají často mnoho parametrů a je obtížné je přizpůsobit klasickým metodám. Bayesovské odhady představují životaschopnou alternativu.

Výše uvedené modely lze fitovat pomocí bayesmh s použitím specifikací více rovnic, které nyní podporují náhodné efekty a latentní faktory. Rovnice mohou být všechny lineární, všechny nelineární nebo směs obou typů. Existují různé způsoby modelování závislosti mezi více výsledky (rovnicemi). Můžete například zahrnout stejné náhodné efekty, ale s různými koeficienty v různých rovnicích, nebo můžete zahrnout různé náhodné efekty, ale korelovat je prostřednictvím rozdělení prior.

Postupujme podle Jonese et al. (2005), kteří použili bayesovský bivariační model růstu ke studiu růstu mláďat rybáka černohlavého. Růst byl měřen pomocí délky křídel y1 a hmotnosti y2. Pro délku křídel se předpokládá lineární model růstu a pro hmotnost logistický model růstu.

kde 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)$,
a (ε1,ε2)∼N(0,Σ0)

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

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

Na úrovni jednotlivce (mláděte) jsou čtyři náhodné efekty: U

$U$

$U$ a V

$V$

$V$
jsou intercepce a míra růstu křídel. V rovnici pro y2

$y_2$

$y_2$,
máme náhodné efekty B

$B$

$B$ a C

$C$

$C$, které představují rychlost a maximální
hmotnost. Předpokládá se, že parametr umístění má tvar dC

$dC$

$dC$, kde d

$d$

$d$
je konstanta. (U,V,C,B)

$(U,V,C,B)$

$(U, V, C, B)$ se řídí vícerozměrným normálním rozdělením s
nestrukturovaná kovariance. Chyby z obou rovnic se řídí
dvourozměrným normálním rozdělením.

Používáme fiktivní soubor dat simulovaný na základě popisu v Jones et al. (2005). Vhodíme dvourozměrný normální model s kovariancí chyb {Sigma0,m}. Čtyřem náhodným efektům je přiřazena vícerozměrná normální priorita s odpovídajícími parametry střední hodnoty a kovariance {Sigma,m}. Prioritním prostředkům je přiřazeno normální rozdělení se střední hodnotou 0 a rozptylem 100. Kovarianční matice mají přiřazeny inverzní Wishartovy priority. Parametru d je přiřazena exponenciální priorita se škálou 1. Pro kovarianční matice a blokové parametry používáme Gibbsovo vzorkování, abychom zvýšili efektivitu. Rovněž určujeme počáteční hodnoty, používáme menší velikost MCMC 2 500 a pro reprodukovatelnost určujeme semeno náhodného čísla.

. 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

Mezi chybovými členy existuje kladná korelace 0,21.

Vypočítáme také korelaci mezi rychlostí růstu křídel a maximální hmotností.

. 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

Odhadovaná korelace je 0,73, což je vysoká hodnota. Měření délky křídla a hmotnosti jsou jednoznačně korelovaná a měla by být modelována společně.

Odkazy

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.