Sinteza dovezilor pentru luarea deciziilor în asistență medicală

Specificați modelul JAGS

Următorul cod specifică modelul JAGS cu efecte aleatorii. Studiul de cote specifice de jurnal, se presupune că sunt probe dintr -o distribuție aleatorie comună: unde D este este media populației de cote de jurnal și este variația dintre studii. Presupunem o distribuție prealabilă uniformă pentru .

Consultați manualul JAGS pentru detalii.

# Model
model_code <- "model {
  # Likelihood
  for (j in 1:Nstud) {
    P(j) <- 1/V(j)      # Calculate precision
    Y(j) ~ dnorm(delta(j),P(j))
    delta(j) ~ dnorm(d,prec)
  }

  # Priors
  d ~ dnorm(0,1.0E-6)
  prec <- 1/tau.sq
  tau.sq <- tau*tau   # between-study variance
  tau ~ dunif(0,10)   # Uniform on SD
}" %>% textConnection

# Note that model_code is a string and the textConnection function is used 
# to pass this string to the jags() function further down in the code.

Acest cod inițializează valorile parametrilor. Rețineți că folosim o funcție pentru a specifica distribuții rezonabile din care vor fi inițializate mai multe lanțuri MCMC.

# Note that the WinBugs code in the book initializes the MCMC algorithm with a single chain as follows.
#initial_values <- list(list(
# delta=c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),d=0,tau=1))

# We use this code  which initializes the MCMC algorithm with three chains 
# whose values are drawn from probability distributions.
set.seed(9999)
inits  <- function() {
  list(
    delta = rnorm(22, 0, 0.5),
    d = rnorm(1, 0, 1),
    tau = runif(1, 0, 3)
  )
}

Acest următor bloc de cod specifică modelul MCMC. Parametrii cheie sunt:

• N.CHAINS: Umberul lanțurilor MCMC să alerge
• n.adapt: ​​numărul de iterații pentru a rula în JAGS fază adaptativă. Vezi Andrieu & Thoms (2008) pentru o discuție despre MCMC adaptiv
• n.iter: numărul de iterații pe lanț, inclusiv ardere
• N.Burnin: Numărul de iterații de aruncat ca ardere
• N.Thin: Rata de subțiere: Fiecare probă KTH va fi aruncată pentru a reduce autocorelația. A se vedea Link & Eaton (2012) pentru o discuție despre subțiere.

Examinați ieșirea

Distribuțiile posterioare pentru parametrii modelului sunt principalele rezultate ale modelului. jags Eșantioane aceste distribuții și calculează media posterioară, eroarea standard, intervalele credibile și statisticile de diagnostic pentru parametrii D, patru dintre parametri și devianța. Valorile prezentate aici sunt foarte apropiate de cele raportate în text. (Rețineți, în scopuri de prezentare, afișăm doar patru valori ale parametri.)

options(width = 999)
res <- out_initial
h_res <- head(res$summary, 4)
t_res <- tail(res$summary, 4)
signif(rbind(h_res, t_res), 3)

            mean     sd     2.5%     25%    50%    75%   97.5% Rhat n.eff overlap0     f
d         -0.247 0.0660 -0.37200 -0.2900 -0.248 -0.203 -0.1130 1.00  1000        0 1.000
tau        0.132 0.0811  0.00849  0.0674  0.125  0.184  0.3080 1.01   405        0 1.000
delta(1)  -0.239 0.1660 -0.57600 -0.3270 -0.244 -0.156  0.1220 1.00  3930        1 0.930
delta(2)  -0.288 0.1580 -0.64600 -0.3670 -0.277 -0.196  0.0112 1.00 15000        1 0.971
delta(20) -0.240 0.1300 -0.50100 -0.3200 -0.243 -0.167  0.0403 1.00 10800        1 0.959
delta(21) -0.325 0.1410 -0.65300 -0.4020 -0.309 -0.231 -0.0831 1.00 11900        0 0.994
delta(22) -0.319 0.1430 -0.65000 -0.3930 -0.302 -0.225 -0.0779 1.00 15000        0 0.993
deviance   7.810 4.5800 -1.22000  4.5100  8.050 11.300 16.0000 1.00  1950        1 0.952

Statistici de diagnostic

Statisticile de diagnostic sunt interpretate după cum urmează:

• Rhatstatistica gelman-rubină, este un diagnostic care compară variația din lanțuri cu variația dintre lanțuri. Valorile apropiate de 1 (de obicei mai puțin de 1,1) indică faptul că lanțurile au convergent.
• n.eff: Oferă o estimare a câte eșantioane independente sunt echivalente probe din lanț. Un număr mai mare sugerează estimări mai fiabile.
• overlap0 = 0 indică faptul că intervalul credibil de 95% nu include 0, ceea ce sugerează un efect semnificativ statistic.
• f este proporția posteriorului cu același semn ca și media.

În plus, The jagsFuncția returnează două statistici alternative de devianță penalizată, criteriul informațiilor despre devianță, DIC și devianța preconizată penalizată, PD, care sunt generate prin intermediul dic.samples() funcţie. Ambele statistici penalizează complexitatea modelului: mai mic este mai bun. Pentru acest model DIC = 18.3016 și PD = 10.4929

După cum s -a calculat jags()DIC este o aproximare la devianța penalizată utilizată atunci când sunt disponibili doar estimatori de puncte, adică, unde este probabilitatea modelului dat datele. Apropierea deține asimptotic atunci când numărul efectiv de parametri este mult mai mic decât dimensiunea eșantionului și când distribuțiile posterioare ale parametrilor modelului sunt distribuite în mod normal. Vezi pagina 6 din rjags PDF pentru detalii și Măsuri bayesiene de complexitate și potrivire a modeluluide Spiegelhalter și colab. (2002) pentru teorie.

Pentru fundal pe PD, vezi Funcții de pierdere penalizate pentru compararea modelului BayesianPlummer (2008).

dev <- data.frame(x = out_initial$pD, y = out_initial$DIC)
names(dev) <- c("pD", "DIC")
dev

        pD      DIC
1 10.49292 18.30155

Rețineți că toate rezultatele pe care le -am discutat pot fi generate simultan prin imprimarea obiectului de ieșire a modelului: out_initial.

jagsUI::densityplot(
  out_initial,
  parameters = c("d", "tau", "deviance", "delta(1)"),
  layout = c(2, 2)
)

colnames <- c("d", sprintf("delta(%d)", 1:22))
caterplot(out_initial$samples(, colnames)) 

Diagnostic MCMC

Un diagnostic cheie MCMC raportat de jags() function din pachetul Jagsui este suficient.ADAPT, o valoare logică care indică dacă numărul de iterații pentru adaptare a fost suficient. Pentru acest model:

cat("Sufficient Adaptation =",
    out_initial$mcmc.info$sufficient.adapt,
    "n")

Sufficient Adaptation = TRUE 

jagsUI::traceplot(
  out_initial,
  parameters = c("d", "tau", "deviance", "delta(1)"),
  layout = c(2, 2)
)