R/ModelIPDnetmeta.R

In htx-r/GenericModelNMA:

modelIPDnetmeta <-function(){

  for (i in 1:np) { # loop through individuals

    # binomial likelihood of 0/1 bianry outcome y for each individual i of study j in arm k
    y[i]~dbern(p[i])

    # logistic transformation - to estimate Odds Ratio (OR)
    logit(p[i]) <- u[studyid[i]]+theta[studyid[i],treat[i]]*(1-equals(treat[i],baseline[i]))
  }

  for(j in 1:nIPD) { # loop through IPD studies
     w[j,1]<- 0
     theta[j,t[j,1]]<- 0
    for (k in 2:(na[j])) { # loop through non-referent IPD arms
      # distribution of random effects
      theta[j,t[j,k]] ~ dnorm(md[j,t[j,k]],precd[j,t[j,k]])

      # accounting for correlation between effect sizes estimated in multi-arm trials
      md[j,t[j,k]]<- mean[j,k] + sw[j,k]
      w[j,k]<- theta[j,t[j,k]]-mean[j,k]
      sw[j,k]<- sum(w[j,1:(k-1)])/(k-1)
      precd[j,t[j,k]]<- prec*2*(k-1)/k

      # consistency equations
      mean[j,k] <-d[t[j,k]] - d[t[j,1]]

    }}

  #** PRIORS

  # prior distribution for log-odds in baseline arm of study j
  for (j in 1:nIPD) {u[j] ~ dnorm(0,0.1)}

  # prior distribution for heterogeneity
  tau ~ dnorm(0,1)%_%T(0,)
  prec<- 1/pow(tau,2)
  tau.sq<- pow(tau,2)

  # prior distribution for basic parameters
  d[ref] <- 0
  for(k in 1:(ref-1)) {d[k] ~ dnorm(0,0.001)}
  for(k in (ref+1):nt) {d[k] ~ dnorm(0,.01)}
}