1. Contabilizarea abandonului informativ în analiza datelor longitudinale.
2. Studiați efectele covariabilelor de bază asupra rezultatelor longitudinale și de supraviețuire.
3. Studiați efectele covariabilelor dependente de timp asupra rezultatului supraviețuirii.

Mai jos, demonstrăm prima aplicație, dar același concept poate fi aplicat celorlalte două.
Vom utiliza o versiune a datelor Scalei Simptomelor Pozitive și Negative (PANSS) dintr-un studiu clinic care compară diferite tratamente medicamentoase pentru schizofrenie (Diggle 1998). Datele conțin scoruri PANSS (panss) de la pacienții care au primit trei tratamente (treat): placebo, haloperidol (referință) și risperidonă (terapie nouă). Scorurile PANSS sunt măsurători pentru tulburarea psihiatrică. Am dori să studiem efectele tratamentelor asupra scorurilor PANSS în timp (week).
Un model considerat pentru scorurile PANSS este un model longitudinal liniar cu efectele tratamentelor, timpul (în săptămâni) și interacțiunea lor și interceptările aleatorii U [id].

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

Deci, cum intră în joc modelul de supraviețuire?
Mulți subiecți s-au retras din acest studiu – doar aproximativ jumătate au finalizat programul complet de măsurare. Mulți subiecți au renunțat din cauza „inadecvării pentru răspuns”, ceea ce sugerează că abandonul poate fi informativ și nu poate fi pur și simplu ignorat în analiză.
Un proces de abandon științific poate fi privit ca un proces de supraviețuire cu un abandon științific informativ (infdrop) ca eveniment de interes și cu timpul până la abandon școlar ca timp de supraviețuire. Deoarece avem date longitudinale, există mai multe observații pe subiect. Deci, timpul de abandon este împărțit în mai multe vrăji în funcție de week și este astfel reprezentat de ora de început (t0) și ora de sfârșit (t1). În punctul de timp rămas t0, o observație (sau o semnificație) este considerată trunchiată la stânga. Vom presupune un model de supraviețuire Weibull pentru procesul de abandon. Abandonul este probabil legat de tratament, deci îl includem ca predictor în modelul de supraviețuire.

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

O modalitate de a explica abandonul informativ este de a include un efect aleatoriu comun între modelele longitudinale și de supraviețuire care ar induce corelația dintre rezultatul longitudinal și procesul de abandon (Henderson, Diggle și 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)))

Am adăugat interceptări aleatorii de la modelul longitudinal la modelul de supraviețuire. De asemenea, am restrâns coeficientul pentru U[id] în prima ecuație la 1. Am făcut acest lucru pentru a sublinia că doar coeficientul pentru U[id] din a doua ecuație va fi estimat. bayesmh face de fapt această presupunere în mod implicit.
Pentru a se potrivi modelului de mai sus, trebuie să specificăm distribuții anterioare pentru parametrii modelului. Avem mulți parametri, deci este posibil să fie nevoie să specificăm priorități oarecum informative. Reamintim că modelele bayesiene, pe lângă parametrii principali ai modelului, estimează, de asemenea, toți parametrii cu efecte aleatorii U [id].
Dacă există un efect al abandonului asupra scorurilor PANSS, ar fi rezonabil să presupunem că ar fi pozitiv. Deci, specificăm o prioritate exponențială cu scara 1 pentru coeficientul lui U[id] în modelul de supraviețuire. Specificăm priorele normale cu media 0 și varianța de 1.000 pentru constanta {panss: _cons} și parametrul Weibull {lnp}. Presupunem priori normali cu medie 0 și varianță 100 pentru alți coeficienți de regresie. Și folosim priori invers-gamma cu formă și scară de 0,01 pentru componentele varianței.
Pentru a îmbunătăți eficiența eșantionării, folosim eșantionarea Gibbs pentru componentele varianței și efectuăm blocarea altor parametri. De asemenea, specificăm valorile inițiale pentru unii parametri utilizând opțiunea 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

Nu ne vom concentra pe interpretarea tuturor rezultatelor acestui model, dar vom comenta coeficientul {t1: U} pentru interceptarea aleatorie partajată. Se estimează că este 0,06, iar intervalul său credibil de 95% nu conține 0. Aceasta înseamnă că există o asociere pozitivă între scorurile PANSS și timpii de abandon: cu cât scorul PANSS este mai mare, cu atât este mai mare șansa de abandon. Deci, simpla analiză longitudinală nu ar fi adecvată pentru aceste date.

Modele de creștere neliniare multivariate

Modelele ierarhice de creștere liniară și neliniară sunt populare în multe discipline, cum ar fi știința sănătății, educația, științele sociale, ingineria și biologia. Modelele de creștere liniare și neliniare multivariate sunt deosebit de utile în științele biologice pentru a studia creșterea speciilor de animale sălbatice, unde creșterea este descrisă prin măsurători multiple care sunt adesea corelate. Astfel de modele au adesea mulți parametri și sunt dificil de adaptat folosind metode clasice. Estimarea bayesiană oferă o alternativă viabilă.
Modelele de mai sus pot fi potrivite de bayesmh folosind specificații de ecuații multiple, care acceptă acum efecte aleatorii și factori latenți. Ecuațiile pot fi toate liniare, toate neliniare sau un amestec din cele două tipuri. Există diferite moduri de a modela dependența între rezultate multiple (ecuații). De exemplu, puteți include aceleași efecte aleatorii, dar cu coeficienți diferiți în ecuații diferite, sau puteți include efecte aleatoare diferite, dar corelați-le prin distribuția anterioară.
Să-i urmărim pe Jones et al. (2005) care au aplicat un model Bayesian de creștere bivariată pentru a studia creșterea puilor cu stern negru. Creșterea a fost măsurată prin lungimea aripii y1 și greutatea y2. Se presupune un model de creștere liniară pentru lungimea aripii, iar un model de creștere logistică se presupune pentru greutate.

where d

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

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

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

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

$U$ and V

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

$y_2$,
we have random effects B

$B$ and C

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

$dC$, where d

$d$
is a constant. (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.

urmeaază o distribuție normală multi variată cu o covarianță nestructurată. Erorile din cele două ecuații urmează o distribuție normală bivariată.
Folosim un set de date fictiv simulat pe baza descrierii din Jones et al. (2005). Ne potrivim unui model normal bivariat cu covarianță de eroare {Sigma0, m}. Celor patru efecte aleatorii li se atribuie o normală multivariată anterioară cu parametrii medii și covarianța corespunzătoare {Sigma, m}. Mijloacelor anterioare li se atribuie distribuții normale cu media 0 și varianța 100. Matricilor de covarianță li se atribuie priori Wishart inverse. Parametrului d i se atribuie o prioritate exponențială cu scala 1. Folosim eșantionarea Gibbs pentru matricile de covarianță și parametrii blocului pentru a îmbunătăți eficiența. De asemenea, specificăm valorile inițiale, folosim o dimensiune MCMC mai mică de 2.500 și specificăm o semință cu număr aleatoriu pentru reproductibilitate.

. 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

Există o corelație pozitivă, 0,21, între termenii de eroare.
De asemenea, calculăm corelația dintre rata de creștere a aripilor și greutatea maximă.

. 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

Corelația estimată este de 0,73, care este ridicată. Măsurările lungimii și greutății aripii sunt clar corelate și trebuie modelate împreună.