Подивимося, як це працює

• Дворівневі моделі
• Нелінійні багаторівневі моделі
• Моделі типу SEM
• Об’єднані моделі поздовжніх даних та даних про виживання
• Багатовимірні нелінійні моделі зростання

Дворівневі моделі

Давайте почнемо з простого і спочатку побудуємо кілька дворівневих моделей.

Дивіться статтю про багаторівневі байєсівські моделі з префіксом bayes. Там ми показуємо, як використовувати bayes: mixed та інші багаторівневі команди bayes для підгонки багаторівневих байєсівських моделей. Давайте повторимо деякі специфікації тут, використовуючи новий синтаксис випадкових ефектів bayesmh.

Ми розглядаємо одні й ті ж дані про результати з математики учнів третього та п’ятого класів з різних шкіл Внутрішнього Лондона (Mortimore et al. 1988). Ми хочемо дослідити вплив школи на результати з математики.

Підіб’ємо до цих даних просту дворівневу модель випадкового перехоплення за допомогою bayesmh. Ми позначимо випадкові перехоплення на рівні школи як U0[school] і включимо їх у специфікацію регресії.

bayesmh вимагає попередніх специфікацій для всіх параметрів, окрім випадкових ефектів. Для випадкових ефектів передбачається нормальний попередній розподіл із середнім значенням 0 і компонентою дисперсії {var_U0}, де U0 – назва випадкового ефекту. Але для {var_U0} ми повинні вказати попередній розподіл. Для коефіцієнтів регресії ми вказуємо нормальні попередні з середнім значенням 0 і дисперсією 10 000, а для компонент дисперсії – обернений гамма-пріор з формою і масштабом 0.01. Ми вказуємо початкову величину для відтворюваності і використовуємо менший розмір 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

Результати подібні до результатів bayes: mixed у розділі Випадкові перехоплення. Ми можемо повторити той самий аналіз після оцінювання з цього розділу після 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 у розділі Випадкові коефіцієнти.

За замовчуванням bayesmh припускає, що випадкові ефекти U0[id] та U1[id] є незалежними апріорі. Але ми можемо змінити це припущення, задавши для них багатовимірний нормальний розподіл з неструктурованою коваріаційною матрицею {Sigma,m}. Ми додатково задамо обернене попереднє Вішарта для коваріації {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) у Випадкових коефіцієнтах. Як і в попередньому розділі, ми можемо використовувати 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))

Але всі ці моделі можна легко підібрати вже за prefix bayes. Далі ми продемонструємо моделі, які не можуть бути підібрані за допомогою bayes:.

Нелінійні багаторівневі моделі

Ви можете підібрати байєсівські одновимірні нелінійні багаторівневі моделі за допомогою bayesmh. Синтаксис bayesmh подібний до синтаксису команди menl, яка підганяє класичні нелінійні моделі зі змішаними ефектами, за винятком того, що bayesmh додатково підтримує перехресні ефекти, такі як UV[id1#id2], та латентні фактори, такі як L[_n].

Трирівнева нелінійна модель в [BAYES] bayesmh.

Моделі типу SEM

Ви можете використовувати bayesmh для підгонки деяких моделей структурних рівнянь (SEM) та споріднених з ними моделей. SEM аналізують декілька змінних і використовують так звані латентні змінні у своїх специфікаціях. Прихована змінна – це псевдозмінна, яка має інше, неспостережуване значення для кожного спостереження. У bayesmh ви вказуєте латентні змінні як L[_n]. Ви можете використовувати різні латентні змінні в різних рівняннях, ви можете використовувати однакові латентні змінні в різних рівняннях, ви можете накладати обмеження на коефіцієнти латентних змінних тощо.

Приклади див. у статтях “Модель латентного зростання” та “Теорія відгуку на товар” у [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] в моделі виживання. Для константи {panss:_cons} та параметра Вейбулла {lnp} ми задаємо нормальні попередні із середнім значенням 0 та дисперсією 1,000. Для інших коефіцієнтів регресії ми припускаємо, що вони є нормальними з математичним сподіванням 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 методу, використовуючи специфікації з декількома рівняннями, які тепер підтримують випадкові ефекти та латентні фактори. Рівняння можуть бути повністю лінійними, повністю нелінійними або сумішшю цих двох типів. Існують різні способи моделювання залежності між кількома результатами (рівняннями). Наприклад, ви можете включити однакові випадкові ефекти, але з різними коефіцієнтами в різні рівняння, або ви можете включити різні випадкові ефекти, але пов’язати їх через попередній розподіл.

Наслідуватимемо Джонса та ін. (2005), які застосували байєсівську двовимірну модель росту для вивчення росту пташенят крячка чорнолобого. Зростання вимірювали за довжиною крила y1 та вагою y2. Для довжини крила було використано лінійну модель росту, а для ваги – логістичну модель росту.

where d

$d$

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

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

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

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

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

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

$U$

$U$ and V

$V$

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

$y_2$

$y_2$,
we have random effects B

$B$

$B$ and C

$C$

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

$dC$

$dC$, where d

$d$

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

$(U,V,C,B)$

$(U, V, C, B)$ follow a multivariate normal distribution with
an unstructured covariance. The errors from the two equations follow a
bivariate normal distribution.

Ми використовуємо фіктивний набір даних, змодельований на основі опису в Jones та ін. (2005). Ми підібрали двовимірну нормальну модель з коваріацією помилки {Sigma0,m}. Чотирьом випадковим ефектам присвоюється багатовимірна нормальна попередня модель з відповідними середніми параметрами та коваріацією {Sigma,m}. Попередні середні мають нормальний розподіл з математичним сподіванням 0 і дисперсією 100. Матриці коваріацій є оберненими попередніми Вішарта. Параметру d присвоєно експоненціальне попереднє значення з масштабом 1. Для підвищення ефективності ми використовуємо вибірку Гіббса для коваріаційних матриць і параметрів блоків. Ми також задаємо початкові значення, використовуємо менший розмір MCMC (2,500) і задаємо випадкові числа для відтворюваності.

. 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, що є високим показником. Вимірювання довжини та ваги крила чітко корелюють між собою і повинні моделюватися спільно.

Посилання

Діггл, П. Дж. 1998. Робота з відсутніми значеннями в поздовжніх дослідженнях. Останні досягнення в статистичному аналізі медичних даних, за ред. Б. С. Еверітта та Г. Данна, 203-228. Лондон: Arnold. Хендерсон, Р., П. Діггл та А. Добсон. 2000. Спільне моделювання поздовжніх вимірювань і даних про час подій. Біостатистика 4: 465-480. Джонс, Г., Р. Д. Кідвелл, А. Д. Л. Ноубл та Д. І. Хеддерлі. 2005. Знайомство з пташенятами: Калібрування та дискримінація в нелінійній багатовимірній ієрархічній моделі росту. Журнал сільськогосподарської, біологічної та екологічної статистики 10: 306-320. Мортимор, П., П. Семмонс, Л. Столл, Д. Льюїс та Р. Екоб. 1988. Шкільні питання: Молодші класи. Сомерсет, Великобританія: Відкриті книги.