R/model.py.hom.eqcor.R

Defines functions model.py.hom.eqcor

model.py.hom.eqcor <- function(prior.type = "unif", rank.prob = TRUE){
if(prior.type == "unif" & rank.prob){
modelstring<-"
model{
 for(i in 1:len){
  y[i] ~ dpois(py[i]*lambda[i])
  lambda[i] <- exp(mu[t[i]] + vi[s[i], t[i]])
 }
 for(j in 1:nstudy){
  vi[j,1:ntrt] ~ dmnorm(zeros[1:ntrt], T[1:ntrt, 1:ntrt])
 }
 for(j in 1:ntrt){
  rate[j] <- exp(mu[j] + pow(sigma, 2)/2)
  lograte[j] <- log(rate[j])
  mu[j] ~ dnorm(0, 0.001)
 }
 sigma ~ dunif(0.0001, c)
 for(j in 1:ntrt){
  for(k in 1:ntrt){ 
   T[j,k] <- 1/sigma^2*ifelse(j == k, diag, offdiag)
  }
 }
 diag <- (1 + (ntrt - 2)*rho)/(1 + (ntrt - 2)*rho - (ntrt - 1)*rho^2)
 offdiag <- (-rho/(1 + (ntrt - 2)*rho - (ntrt - 1)*rho^2))
 rho ~ dunif(-1/(ntrt - 1), 0.9999)
 for(j in 1:ntrt){        
  for(k in 1:ntrt){
   ratio[j,k] <- rate[j]/rate[k]
   logratio[j,k] <- log(ratio[j,k])
  }
 }
 rk[1:ntrt] <- (ntrt + 1 - rank(rate[]))*ifelse(higher.better, 1, 0) + (rank(rate[]))*ifelse(higher.better, 0, 1)
 for(i in 1:ntrt){
  rank.prob[1:ntrt,i] <- equals(rk[], i)
 }
}
"
}

if(prior.type == "unif" & !rank.prob){
modelstring<-"
model{
 for(i in 1:len){
  y[i] ~ dpois(py[i]*lambda[i])
  lambda[i] <- exp(mu[t[i]] + vi[s[i], t[i]])
 }
 for(j in 1:nstudy){
  vi[j,1:ntrt] ~ dmnorm(zeros[1:ntrt], T[1:ntrt, 1:ntrt])
 }
 for(j in 1:ntrt){
  rate[j] <- exp(mu[j] + pow(sigma, 2)/2)
  lograte[j] <- log(rate[j])
  mu[j] ~ dnorm(0, 0.001)
 }
 sigma ~ dunif(0.0001, c)
 for(j in 1:ntrt){
  for(k in 1:ntrt){ 
   T[j,k] <- 1/sigma^2*ifelse(j == k, diag, offdiag)
  }
 }
 diag <- (1 + (ntrt - 2)*rho)/(1 + (ntrt - 2)*rho - (ntrt - 1)*rho^2)
 offdiag <- (-rho/(1 + (ntrt - 2)*rho - (ntrt - 1)*rho^2))
 rho ~ dunif(-1/(ntrt - 1), 0.9999)
 for(j in 1:ntrt){        
  for(k in 1:ntrt){
   ratio[j,k] <- rate[j]/rate[k]
   logratio[j,k] <- log(ratio[j,k])
  }
 }
}
"
}

if(prior.type == "invgamma" & rank.prob){
modelstring<-"
model{
 for(i in 1:len){
  y[i] ~ dpois(py[i]*lambda[i])
  lambda[i] <- exp(mu[t[i]] + vi[s[i], t[i]])
 }
 for(j in 1:nstudy){
  vi[j,1:ntrt] ~ dmnorm(zeros[1:ntrt], T[1:ntrt, 1:ntrt])
 }
 for(j in 1:ntrt){
  rate[j] <- exp(mu[j] + pow(sigma, 2)/2)
  lograte[j] <- log(rate[j])
  mu[j] ~ dnorm(0, 0.001)
 }
 sigma <- 1/sqrt(inv.sig.sq)
 inv.sig.sq ~ dgamma(a, b)
 for(j in 1:ntrt){
  for(k in 1:ntrt){ 
   T[j,k] <- 1/sigma^2*ifelse(j == k, diag, offdiag)
  }
 }
 diag <- (1 + (ntrt - 2)*rho)/(1 + (ntrt - 2)*rho - (ntrt - 1)*rho^2)
 offdiag <- (-rho/(1 + (ntrt - 2)*rho - (ntrt - 1)*rho^2))
 rho ~ dunif(-1/(ntrt - 1), 0.9999)
 for(j in 1:ntrt){        
  for(k in 1:ntrt){
   ratio[j,k] <- rate[j]/rate[k]
   logratio[j,k] <- log(ratio[j,k])
  }
 }
 rk[1:ntrt] <- (ntrt + 1 - rank(rate[]))*ifelse(higher.better, 1, 0) + (rank(rate[]))*ifelse(higher.better, 0, 1)
 for(i in 1:ntrt){
  rank.prob[1:ntrt,i] <- equals(rk[], i)
 }
}
"
}

if(prior.type == "invgamma" & !rank.prob){
modelstring<-"
model{
 for(i in 1:len){
  y[i] ~ dpois(py[i]*lambda[i])
  lambda[i] <- exp(mu[t[i]] + vi[s[i], t[i]])
 }
 for(j in 1:nstudy){
  vi[j,1:ntrt] ~ dmnorm(zeros[1:ntrt], T[1:ntrt,1:ntrt])
 }
 for(j in 1:ntrt){
  rate[j] <- exp(mu[j] + pow(sigma, 2)/2)
  lograte[j] <- log(rate[j])
  mu[j] ~ dnorm(0, 0.001)
 }
 sigma <- 1/sqrt(inv.sig.sq)
 inv.sig.sq ~ dgamma(a, b)
 for(j in 1:ntrt){
  for(k in 1:ntrt){ 
   T[j,k] <- 1/sigma^2*ifelse(j == k, diag, offdiag)
  }
 }
 diag <- (1 + (ntrt - 2)*rho)/(1 + (ntrt - 2)*rho - (ntrt - 1)*rho^2)
 offdiag <- (-rho/(1 + (ntrt - 2)*rho - (ntrt - 1)*rho^2))
 rho ~ dunif(-1/(ntrt - 1), 0.9999)
 for(j in 1:ntrt){        
  for(k in 1:ntrt){
   ratio[j,k] <- rate[j]/rate[k]
   logratio[j,k] <- log(ratio[j,k])
  }
 }
}
"
}

if(!is.element(prior.type, c("unif", "invgamma"))){
  stop("specified prior type is wrong.")
}

return(modelstring)
}

Try the pcnetmeta package in your browser

Any scripts or data that you put into this service are public.

pcnetmeta documentation built on Aug. 31, 2022, 9:08 a.m.