Нека видим как работи

• Модели на две нива
• Нелинейни модели на много нива
• Модели от типа SEM
• Съвместни модели на надлъжни данни и данни за преживяемост
• Многомерни нелинейни модели на растеж

Модели на две нива

Нека започнем просто и първо приложим няколко модела на две нива.

Вижте Байесовите модели на няколко нива, като използвате префикса bayes. Там показваме как да използваме bayes: mixed и други команди на bayes multilevel за напасване на байесови модели на няколко нива. Нека повторим някои от спецификациите тук, като използваме новия синтаксис за случайни ефекти на bayesmh.

Разглеждаме същите данни за резултатите по математика на ученици от трети и пети клас от различни училища във Вътрешен Лондон (Mortimore et al. 1988). Искаме да изследваме ефекта на училището върху резултатите по математика.

Нека приложим прост модел на две нива със случайни отклонения към тези данни, като използваме bayesmh. Уточняваме случайните интерпекти на ниво училище като U0[school] и ги включваме в регресионната спецификация.

bayesmh изисква предварителни спецификации за всички параметри с изключение на случайните ефекти. За случайните ефекти се приема нормално предварително разпределение със средна стойност 0 и компонент на дисперсията {var_U0}, където U0 е името на случайния ефект. Но ние трябва да посочим приоритета за {var_U0}. За регресионните коефициенти задаваме нормални приоритети със средна стойност 0 и дисперсия 10 000, а за дисперсионните компоненти – обратен гама-приор с форма и мащаб 0,01. Посочваме семе за възпроизводимост и използваме по-малък MCMC размер от 1000.

. 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

Резултатите са подобни на тези от bayes: mixed в Random intercepts. Можем да повторим същия анализ след оценяване от този раздел след bayesmh, включително показването и графиките на случайните ефекти. В допълнение към основните параметри на модела, байесовите модели оценяват и всички случайни ефекти. Това е в контраст с френетичния анализ, при който случайните ефекти не се оценяват съвместно с параметрите на модела, а се прогнозират след оценяването под условие на параметрите на модела.

След това включваме случайни коефициенти за math като c.math3#U1[school]. Допълнително трябва да посочим приоритет за компонента на дисперсията {var_U1}. Добавяме към списъка с компонентите на дисперсията, като използваме инверсен гама-приор.

. 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

Резултатите ни са подобни на тези от bayes: mixed в коефициентите Random.

По подразбиране bayesmh приема, че случайните ефекти U0[id] и U1[id] са априори независими. Но ние можем да променим това, като зададем многомерното нормално разпределение за тях с неструктурирана ковариационна матрица {Sigma,m}. Допълнително задаваме обратен Wishart приоритет за ковариацията {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

Резултатите отново са подобни на тези от bayes: mixed, covariance(unstructured) в Random coefficients. Точно както в този раздел, бихме могли да използваме bayesstats ic след bayesmh, за да сравним неструктурираните и независимите модели.

Можем също така да използваме новата предварителна опция mvnexchangeable(), за да зададем обменна ковариационна структура със случайни ефекти вместо неструктурирана. Обменната ковариация се характеризира с два параметъра – обща дисперсия и корелация. За тези параметри се посочват допълнителни прийоми.

. 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

Можем да пригодим други модели от многостепенните модели на Бейс, използвайки префикса bayes, като използваме bayesmh. Например можем да пригодим модела за оцеляване на три нива, като използваме

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

и логистичния модел с кръстосани ефекти, като се използва

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

Но всички тези модели могат лесно да се пригодят чрез префикса на Бейс. По-нататък ще демонстрираме модели, които не могат да бъдат подбрани с префикса на Бейс:.

Нелинейни модели на много нива

Можете да приложите едномерни нелинейни многостепенни модели на Бейс с помощта на bayesmh. Синтаксисът на bayesmh е същият като този на командата menl, която подбира класически нелинейни модели със смесени ефекти, с изключение на това, че bayesmh допълнително поддържа кръстосани ефекти като UV[id1#id2] и латентни фактори като L[_n].

Вижте тристепенния нелинеен модел в [BAYES] bayesmh.

Модели от типа SEM

Можете да използвате bayesmh, за да подберете някои модели със структурни уравнения (SEM) и свързани с тях модели. SEM анализират множество променливи и използват т.нар. латентни променливи в спецификациите си. Латентната променлива е псевдопроменлива, която има различна, ненаблюдавана стойност за всяко наблюдение. С bayesmh задавате латентни променливи като L[_n]. Можете да използвате различни латентни променливи в различни уравнения, можете да споделяте едни и същи латентни променливи в различни уравнения, можете да поставяте ограничения върху коефициентите на латентните променливи и др.

За примери вижте Latent growth model (модел на латентния растеж) и Item response theory (теория на отговора на елемента) в [BAYES] bayesmh.

Съвместни модели на надлъжни данни и данни за преживяемост

Можете да използвате bayesmh, за да моделирате множество резултати от различни типове. Един такъв пример са съвместните модели на резултати от надлъжни изследвания и преживяемост. Тези модели са популярни в практиката поради трите си основни приложения:

1. Отчитане на информационното отпадане при анализа на надлъжни данни.
2. Проучване на ефектите на изходните ковариати върху резултатите от надлъжните изследвания и преживяемостта.
3. Проучване на ефектите на зависимите от времето ковариативите върху резултатите от преживяемостта.

По-долу ще демонстрираме първото приложение, но същата концепция може да се приложи и за другите две.

Ще използваме версия на скалата за положителни и отрицателни симптоми (PANSS), данни от клинично изпитване, сравняващо различни начини на лечение с лекарства за шизофрения (Diggle 1998). Данните съдържат оценки по PANSS (panss) от пациенти, които са получили три вида лечение (treat): плацебо, халоперидол (референтна терапия) и рисперидон (нова терапия). Оценките PANSS са измервания за психиатрично разстройство. Бихме искали да проучим ефекта на лечението върху резултатите от PANSS за определен период от време (week).

Моделът, който се разглежда за оценките по PANSS, е надлъжен линеен модел с ефектите на леченията, времето (в седмици) и тяхното взаимодействие и случайни интерцепции U[id].

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

И така, как моделът на оцеляване се прилага тук?

Много участници се оттеглиха от това проучване – само около половината завършиха пълния график на измерванията. Много субекти са се отказали поради „неадекватност на отговора“, което предполага, че отпадането може да има информативен характер и не може просто да се пренебрегне в анализа.

Процесът на отпадане може да се разглежда като процес на оцеляване с информативно отпадане (infdrop) като събитие от интерес и с време до отпадането като време на оцеляване. Тъй като разполагаме с надлъжни данни, за всеки субект има множество наблюдения. Така че времето на отпадане е разделено на няколко периода според седмицата и по този начин е представено от началния момент (t0) и крайния момент (t1). В лявата времева точка t0 наблюдението (или заклинанието) се счита за ляво прекъснато. За процеса на отпадане ще приемем модел на преживяемост на Вейбул. Отпадането вероятно е свързано с лечението, така че го включваме като предиктор в модела на оцеляването.

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

Един от начините за отчитане на информативното отпадане е включването на общ случаен ефект между моделите на надлъжното наблюдение и оцеляването, който би предизвикал корелация между резултата от надлъжното наблюдение и процеса на отпадане (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)))

Добавихме случайни интерпекти от модела на надлъжното наблюдение към модела на оцеляването. Също така ограничихме коефициента за U[id] в първото уравнение до 1. Направихме това, за да подчертаем, че ще бъде оценен само коефициентът за U[id] във второто уравнение. bayesmh всъщност прави това допускане автоматично по подразбиране.

За да приложим горния модел, трябва да посочим предварителни разпределения за параметрите на модела. Имаме много параметри, така че може да се наложи да посочим донякъде информативни приоритети. Припомнете си, че байесовите модели, освен основните параметри на модела, оценяват и всички параметри на случайните ефекти U[id].

Ако има ефект на отпадането върху резултатите от PANSS, би било разумно да се предположи, че той ще бъде положителен. Затова определяме експоненциален приоритет със скала 1 за коефициента U[id] в модела на оцеляването. Определяме нормални приоритети със средна стойност 0 и дисперсия 1 000 за константата {panss:_cons} и параметъра на Weibull {lnp}. Приемаме нормални прийоми със средна стойност 0 и дисперсия 100 за другите регресионни коефициенти. А за компонентите на дисперсията използваме инверсно-гама приори с форма и мащаб 0,01.

За да подобрим ефективността на извадката, използваме извадка на Гибс за компонентите на дисперсията и извършваме блокиране на другите параметри. Също така задаваме начални стойности за някои параметри, като използваме опцията 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

Няма да се спираме на тълкуването на всички резултати от този модел, но ще коментираме коефициента {t1:U} за споделената случайна интерцепция. Той се оценява на 0,06, а неговият 95% надежден интервал не съдържа 0. Това означава, че има положителна връзка между резултатите по PANSS и времето на отпадане: колкото по-висок е резултатът по PANSS, толкова по-голяма е вероятността за отпадане. Така че обикновеният надлъжен анализ не би бил подходящ за тези данни.

Многомерни нелинейни модели на растеж

Йерархичните линейни и нелинейни модели на растежа са популярни в много дисциплини, като например в областта на здравеопазването, образованието, социалните науки, инженерството и биологията. Многомерните линейни и нелинейни модели на растежа са особено полезни в биологичните науки за изучаване на растежа на животинските видове, където растежът се описва с множество измервания, които често са корелирани. Такива модели често имат много параметри и е трудно да се пригодят с помощта на класическите методи. Байесовото оценяване предоставя жизнеспособна алтернатива.

Горепосочените модели могат да бъдат подбрани от bayesmh, като се използват спецификации с множество уравнения, които сега поддържат случайни ефекти и латентни фактори. Уравненията могат да бъдат изцяло линейни, изцяло нелинейни или смес от двата типа. Съществуват различни начини за моделиране на зависимостта между множество резултати (уравнения). Например можете да включите едни и същи случайни ефекти, но с различни коефициенти в различни уравнения, или можете да включите различни случайни ефекти, но да ги корелирате чрез предварително разпределение.

Нека да последваме Jones et al. (2005), които прилагат двувариантен модел на растежа по метода на Бейс, за да изследват растежа на пиленца на черночела рибарка. Растежът е измерен чрез дължината на крилата y1 и теглото y2. За дължината на крилата е приет линеен модел на растеж, а за теглото – логистичен модел на растеж.

където d

$d$

$d$ е фиксиран параметър, (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)$,
и (ε1,ε2)∼N(0,Σ0)

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

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

Има четири случайни ефекта на ниво индивид (пиле): U

$U$

$U$ и V

$V$

$V$
са пресечната точка и темпът на нарастване на крилата. В уравнението за y2

$y_2$

$y_2$,
имаме случайни ефекти B

$B$

$B$ и C

$C$

$C$, които представляват скоростта и максималната
тегло. Приема се, че параметърът на местоположението има формата dC

$dC$

$dC$, където d

$d$

$d$
е константа. (U,V,C,B)

$(U,V,C,B)$

$(U, V, C, B)$ следва многомерното нормално разпределение с
неструктурирана ковариация. Грешките от двете уравнения следват
двумерното нормално разпределение.

Използваме фиктивен набор от данни, симулиран въз основа на описанието в Jones et al. (2005). Прилагаме двумерен нормален модел с ковариация на грешката {Sigma0,m}. На четирите случайни ефекта се присвоява многомерен нормален приоритет със съответните средни параметри и ковариация {Sigma,m}. На предшестващите средства се присвояват нормални разпределения със средна стойност 0 и дисперсия 100. На ковариационните матрици са зададени обратни Wishart приоритети. На параметъра d се присвоява експоненциален приоритет със скала 1. Използваме извадка на Гибс за ковариационните матрици и блоковите параметри, за да подобрим ефективността. Посочваме също така начални стойности, използваме по-малък MCMC размер от 2500 и посочваме семена от случайни числа за възпроизводимост.

. 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

Налице е положителна корелация (0,21) между членовете на грешката.

Изчисляваме и корелацията между скоростта на растеж на крилата и максималното тегло.

. 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

Оценената корелация е 0,73, което е висока стойност. Измерванията на дължината на крилото и теглото са ясно свързани и трябва да се моделират съвместно.

Препратки

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.