Байесовское многоуровневое моделирование

Новое в

Байесовское многоуровневое моделирование: Нелинейное, совместное, SEM-подобное…

Многоуровневые модели используются во многих дисциплинах для моделирования групповых эффектов, которые могут возникать на разных уровнях иерархии. Подумайте о регионах, штатах, вложенных в регионы, и компаниях, вложенных в штаты внутри регионов. Или подумайте о больницах, врачах, вложенных в больницы, и пациентах, вложенных в врачей, вложенных в больницы.

В дополнение к стандартным преимуществам байесовского анализа, байесовское многоуровневое моделирование предлагает множество преимуществ для данных с небольшим количеством групп или панелей. Кроме того, оно предоставляет принципиальные способы сравнения эффектов в различных группах с помощью апостериорных распределений эффектов. Более подробное обсуждение смотрите здесь.

bayesmh имеет новый синтаксис случайных эффектов, который упрощает подгонку байесовских многоуровневых моделей. Это открывает возможности для подбора новых классов многоуровневых моделей. Вы можете легче подгонять одномерные линейные и нелинейные многоуровневые модели. А теперь вы можете использовать многомерные линейные и нелинейные многоуровневые модели!

Подумайте о нелинейных моделях со смешанными эффектами, для которых подходит menl, или о некоторых моделях SEM, для которых подходят sem и gsem, или о многомерных нелинейных моделях, содержащих случайные эффекты, которые на данный момент не могут быть подогнаны ни одной из существующих команд Stata. Теперь вы можете подобрать байесовские аналоги этих и других моделей с помощью bayesmh.

Основные моменты

Исходы: непрерывные, бинарные, порядковые, счетные, выживаемость, …
Удобная спецификация случайных эффектов
Случайные интерцепты и коэффициенты
Вложенные и перекрестные эффекты
Латентные факторы
Многочисленные уровни иерархии
Линейные и нелинейные модели
Одномерные и многомерные модели
Совместные, многомерные модели роста и SEM-подобные модели
Многомерные нелинейные модели случайных эффектов
Гибкие предварительные распределения случайных эффектов
Постериорные распределения случайных эффектов
Диагностика MCMC, включая множественные цепи
Полная поддержка функций байесовской постоценки

Давайте посмотрим, как это работает

Двухуровневые модели
Нелинейные многоуровневые модели
Модели типа 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 — 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 in 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 in Random coefficients.

По умолчанию 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) в Random coefficients. Как и в этом разделе, мы можем использовать bayesstats ic после bayesmh для сравнения неструктурированной и независимой моделей.

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

. 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]. Вы можете использовать разные латентные переменные в разных уравнениях, вы можете использовать одни и те же латентные переменные в разных уравнениях, вы можете накладывать ограничения на коэффициенты латентных переменных и многое другое.

Примеры см. в разделах Латентная модель роста и Item response theory в [BAYES] bayesmh.

Совместные модели продольных данных и данных о выживании

Вы можете использовать bayesmh для моделирования нескольких исходов разных типов. Одним из таких примеров являются совместные модели продольных исходов и исходов выживания. Эти модели популярны на практике благодаря трем основным областям применения:

Учет информативного отсева при анализе продольных данных.
Исследование влияния исходных ковариаций на продольные результаты и результаты выживания.
Изучение влияния зависящих от времени ковариаций на результаты выживания.

Ниже мы продемонстрируем первое применение, но ту же концепцию можно применить и к двум другим.

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

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

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

Как же здесь работает модель выживания?

Многие испытуемые выбыли из исследования — только около половины прошли полный цикл измерений. Многие испытуемые выбыли по причине «неадекватности ответа», что говорит о том, что отсев может быть информативным и не может быть просто проигнорирован в анализе.

Процесс выбывания можно рассматривать как процесс выживания с информативным выбыванием (infdrop) в качестве интересующего события и временем до выбывания в качестве времени выживания. Поскольку у нас продольные данные, на одного испытуемого приходится несколько наблюдений. Поэтому время отсева разбивается на несколько периодов по week и, таким образом, представлено временем начала (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, разумно предположить, что он будет положительным. Поэтому для коэффициента U[id] в модели выживания мы задаем экспоненциальный приоритет с размахом 1. Для константы {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. Для длины крыла предполагается линейная модель роста, а для веса — логистическая модель роста.

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

где $d$

d

$d$ является фиксированным параметром, $(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)$ ,
и $(ε_{1}, ε_{2}) \sim N (0, Σ_{0})$

(ε_{1}, ε_{2}) \sim N (0, Σ_{0})

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

Существует четыре случайных эффекта на уровне индивидуума (цыпленка): $U$

U

$U$ и $V$

V

$V$
являются перехватом и темпом роста для крыльев. В уравнении для $y_{2}$

y_{2}

$y_2$ ,
у нас есть случайные эффекты $B$

B

$B$ и $C$

C

$C$ , которые представляют собой норму и максимальный
вес. Предполагается, что параметр местоположения имеет вид $d C$

d C

$dC$ , где $d$

d

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

(U, V, C, B)

$(U, V, C, B)$ следуют многомерному нормальному распределению с
неструктурированной ковариацией. Ошибки от двух уравнений следуют
двумерному нормальному распределению.

Мы используем вымышленный набор данных, смоделированный на основе описания в Jones et al. (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, что является высоким показателем. Измерения длины крыла и веса явно коррелируют и должны моделироваться совместно.

Ссылки

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.

Просмотреть все новые функции

Заказать Stata 17

Обновление