Modelle vom Typ SEM

Sie können Bayesmh verwenden, um einige Strukturgleichungsmodelle (SEMs) und damit verwandte Modelle anzupassen. SEMs analysieren mehrere Variablen und verwenden sogenannte latente Variablen in ihren Spezifikationen. Eine latente Variable ist eine Pseudovariable, die für jede Beobachtung einen anderen, unbeobachteten, Wert hat. Mit bayesmh geben Sie latente Variablen als L[_n] an. Sie können verschiedene latente Variablen in verschiedenen Gleichungen verwenden, Sie können dieselben latenten Variablen in verschiedenen Gleichungen verwenden, Sie können Beschränkungen für die Koeffizienten der latenten Variablen festlegen und vieles mehr.
Für Beispiele siehe Latentes Wachstumsmodell und Item-Response-Theorie in [BAYES] bayesmh.

Gemeinsame Modelle von Längs- und Überlebensdaten

Sie können Bayesmh verwenden, um mehrere Outcomes unterschiedlicher Typen zu modellieren. Gemeinsame Modelle von Längs- und Überlebensergebnissen sind ein solches Beispiel. Diese Modelle sind in der Praxis aufgrund ihrer drei Hauptanwendungen beliebt:

1. 1.Berücksichtigen Sie informative Dropouts bei der Analyse von Längsschnittdaten.
2. 2.Untersuchen Sie die Auswirkungen von Baseline-Kovariaten auf Längsschnitt- und Überlebensergebnisse.
3. 3.Untersuchung der Effekte von zeitabhängigen Kovariaten auf das Überlebensergebnis.

Im Folgenden wird die erste Anwendung demonstriert, aber das gleiche Konzept kann auch auf die beiden anderen angewendet werden.
Wir verwenden eine Version der Positive and Negative Symptom Scale (PANSS), Daten aus einer klinischen Studie zum Vergleich verschiedener medikamentöser Behandlungsmethoden für Schizophrenie (Diggle 1998). Die Daten enthalten PANSS-Scores (panss) von Patienten, die drei Behandlungen (treat) erhielten: Placebo, Haloperidol (Referenz) und Risperidon (neue Therapie). PANSS-Scores sind Messwerte für psychiatrische Störungen. Wir möchten die Auswirkungen der Behandlungen auf die PANSS-Scores im Zeitverlauf (Woche) untersuchen.
Als Modell für die PANSS-Scores wird ein lineares Längsschnittmodell mit den Effekten der Behandlungen, der Zeit (in Wochen) und deren Interaktion sowie zufälligen Intercepts U[id] betrachtet.

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

Wie kommt nun das Überlebensmodell hier ins Spiel?
Viele Probanden zogen sich aus dieser Studie zurück – nur etwa die Hälfte beendete den gesamten Messplan. Viele Probanden schieden wegen „unzureichender Antwort“ aus, was darauf hindeutet, dass der Dropout informativ sein kann und in der Analyse nicht einfach ignoriert werden kann.
Ein Dropout-Prozess kann als ein Überlebensprozess mit einem informativen Dropout (infdrop) als Ereignis von Interesse und mit der Zeit bis zum Dropout als Überlebenszeit betrachtet werden. Da es sich um Längsschnittdaten handelt, gibt es mehrere Beobachtungen pro Proband. Die Dropout-Zeit wird also in mehrere Zeitabschnitte nach Wochen aufgeteilt und wird somit durch den Anfangszeitpunkt (t0) und den Endzeitpunkt (t1) dargestellt. Zum linken Zeitpunkt t0 wird eine Beobachtung (oder ein Spell) als links-abgeschnitten betrachtet. Wir werden ein Weibull-Überlebensmodell für den Dropout-Prozess annehmen. Der Dropout hängt wahrscheinlich mit der Behandlung zusammen, daher nehmen wir ihn als Prädiktor in das Überlebensmodell auf.

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

Eine Möglichkeit, den informativen Dropout zu berücksichtigen, besteht darin, einen gemeinsamen Zufallseffekt zwischen dem Längsschnitt- und dem Überlebensmodell einzuschließen, der eine Korrelation zwischen dem Längsschnitt-Ergebnis und dem Dropout-Prozess induzieren würde (Henderson, Diggle und 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)))

Wir haben zufällige Achsen aus dem Längsschnittmodell zum Überlebensmodell hinzugefügt. Außerdem haben wir den Koeffizienten für U[id] in der ersten Gleichung auf 1 beschränkt. Wir haben dies getan, um zu betonen, dass nur der Koeffizient für U[id] in der zweiten Gleichung geschätzt wird. bayesmh macht diese Annahme eigentlich standardmäßig automatisch.
Um das obige Modell anzupassen, müssen wir Prioritätsverteilungen für die Modellparameter angeben. Da wir viele Parameter haben, müssen wir möglicherweise etwas informative Prioritäten angeben. Erinnern Sie sich, dass Bayes’sche Modelle zusätzlich zu den Hauptmodellparametern auch alle Parameter der zufälligen Effekte U[id] schätzen.
Wenn es einen Effekt des Dropouts auf die PANSS-Scores gibt, wäre es vernünftig anzunehmen, dass er positiv ist. Also spezifizieren wir einen Exponentialprior mit Skala 1 für den Koeffizienten von U[id] im Überlebensmodell. Wir spezifizieren Normalprioritäten mit Mittelwert 0 und Varianz von 1.000 für die Konstante {panss:_cons} und Weibull-Parameter {lnp}. Wir nehmen Normalprioritäten mit Mittelwert 0 und Varianz 100 für andere Regressionskoeffizienten an. Und wir verwenden inverse Gamma-Prioritäten mit Form und Skala von 0,01 für die Varianzkomponenten.
Um die Sampling-Effizienz zu verbessern, verwenden wir Gibbs-Sampling für die Varianzkomponenten und führen eine Blockierung der anderen Parameter durch. Wir geben auch Anfangswerte für einige Parameter an, indem wir die Option initial() verwenden.

. 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

Wir werden uns nicht auf die Interpretation aller Ergebnisse dieses Modells konzentrieren, aber wir werden den Koeffizienten {t1:U} für den gemeinsamen zufälligen Intercept kommentieren. Er wird auf 0,06 geschätzt, und sein 95%-Glaubwürdigkeitsintervall enthält nicht 0. Das bedeutet, dass es eine positive Assoziation zwischen PANSS-Scores und Dropout-Zeiten gibt: je höher der PANSS-Score, desto höher die Dropout-Wahrscheinlichkeit. Eine einfache Längsschnittanalyse wäre also für diese Daten nicht geeignet.

Multivariate nichtlineare Wachstumsmodelle

Hierarchische lineare und nichtlineare Wachstumsmodelle sind in vielen Disziplinen beliebt, z. B. in den Gesundheitswissenschaften, im Bildungswesen, in den Sozialwissenschaften, im Ingenieurwesen und in der Biologie. Multivariate lineare und nichtlineare Wachstumsmodelle sind besonders in den Biowissenschaften nützlich, um das Wachstum von Wildtierarten zu untersuchen, bei denen das Wachstum durch mehrere Messungen beschrieben wird, die oft korreliert sind. Solche Modelle haben oft viele Parameter und sind mit klassischen Methoden nur schwer zu erfassen. Die Bayes’sche Schätzung bietet eine praktikable Alternative.
Die obigen Modelle können von Bayesmh unter Verwendung von Mehrfachgleichungsspezifikationen angepasst werden, die jetzt zufällige Effekte und latente Faktoren unterstützen. Die Gleichungen können alle linear, alle nichtlinear oder eine Mischung aus beiden Typen sein. Es gibt verschiedene Möglichkeiten, die Abhängigkeit zwischen mehreren Ergebnissen (Gleichungen) zu modellieren. Sie können z. B. dieselben zufälligen Effekte einbeziehen, aber mit unterschiedlichen Koeffizienten in verschiedenen Gleichungen, oder Sie können verschiedene zufällige Effekte einbeziehen, diese aber durch die Prioritätsverteilung korrelieren.
Folgen wir Jones et al. (2005), die ein bivariates Bayes’sches Wachstumsmodell anwandten, um das Wachstum von Schwarzstirnseeschwalbenküken zu untersuchen. Das Wachstum wurde durch die Flügellänge y1 und das Gewicht y2 gemessen. Für die Flügellänge wird ein lineares Wachstumsmodell angenommen, für das Gewicht ein logistisches Wachstumsmodell.

wobeid

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

Es gibt vier Zufallseffekte auf der individuellen (Küken-)Ebene: U U

$U$ and V

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

$y_2$,
wir haben zufällige Effekte  B

$B$ and C

$C$, , die die Rate und das maximale
Gewicht darstellen. Der Standortparameter wird in der Form angenommen dC

$dC$, where d

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

$(U, V, C, B)$ folgen einer multivariaten Normalverteilung mit
einer unstrukturierten Kovarianz. Die Fehler aus den beiden Gleichungen folgen einer
bivariaten Normalverteilung.

Wir verwenden einen fiktiven Datensatz, der basierend auf der Beschreibung in Jones et al. (2005) simuliert wurde. Wir passen ein bivariates Normalmodell mit der Fehlerkovarianz {Sigma0,m} an. Den vier zufälligen Effekten wird ein multivariater Normalprior mit den entsprechenden Mittelwertparametern und der Kovarianz {Sigma,m} zugewiesen. Den prioren Mittelwerten werden Normalverteilungen mit Mittelwert 0 und Varianz 100 zugeordnet. Den Kovarianzmatrizen werden inverse Wishart-Prioren zugewiesen. Dem Parameter d wird ein Exponentialprior mit der Skalierung 1 zugewiesen. Wir verwenden Gibbs-Sampling für Kovarianzmatrizen und Blockparameter, um die Effizienz zu verbessern. Wir geben auch Anfangswerte an, verwenden eine kleinere MCMC-Größe von 2.500 und geben einen Zufallszahlen-Seed für die Reproduzierbarkeit an.

. 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

Es besteht eine positive Korrelation, 0,21, zwischen den Fehlertermen.
Wir berechnen auch die Korrelation zwischen der Geschwindigkeit des Flügelwachstums und dem Maximalgewicht.

. 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

Die geschätzte Korrelation beträgt 0,73, was hoch ist. Die Messungen von Flügellänge und Gewicht sind eindeutig korreliert und sollten gemeinsam modelliert werden.