I'm trying to run my below model code on OpenBUGS but I keep getting the error "Sorry something went wrong in procedure UpdaterStd.Sample in module UpdaterSCAM" when trying to update the model.
My model is supposed to be an adapted identity link to a log link that is adjusted for time (linear and quadratic term for time centered at 2 hours).
Any help is appreciated.
# Normal likelihood, log link, consistency model, with linear and quadratic terms of time – 2
# Random effects model for multi-arm trials
model{
# *** PROGRAM STARTS
for(i in 1:ns){ # LOOP THROUGH STUDIES
w[i,1] <- 0 # adjustment for multi-arm trials is zero for control arm
delta[i,1] <- 0 # treatment effect is zero for control arm
mu[i] ~ dnorm(0,.0001) # vague priors for all trial baselines
for (k in 1:na[i]) { # LOOP THROUGH ARMS
var[i,k] <- pow(se[i,k],2) # calculate variances
prec[i,k] <- 1/var[i,k] # set precisions
y[i,k] ~ dnorm(theta[i,k],prec[i,k]) # normal likelihood
log(theta[i,k]) <- mu[i] + delta[i,k] + (time[i,k] - 2) + pow((time[i,k] - 2),2) # model for linear and quadratic predictor with time effect centered at 2 hours
dev[i,k] <- (y[i,k]-theta[i,k])*(y[i,k]-theta[i,k])*prec[i,k] #Deviance contribution
}
resdev[i] <- sum(dev[i,1:na[i]]) # summed residual deviance contribution for this trial
for (k in 2:na[i]) { # LOOP THROUGH ARMS
delta[i,k] ~ dnorm(md[i,k],taud[i,k]) # trial-specific LOR distributions
md[i,k] <- d[t[i,k]] - d[t[i,1]] + sw[i,k] # log ratio of means of the kth treatment and 1st treatment
taud[i,k] <- tau*2*(k-1)/k # precision of treat effects distributions (with multi-arm trial correction)
w[i,k] <- (delta[i,k] - d[t[i,k]] + d[t[i,1]]) # adjustment for multi-arm RCTs
sw[i,k] <- sum(w[i,1:k-1])/(k-1) # cumulative adjustment for multi-arm trials
}
}
totresdev <- sum(resdev[]) #Total Residual Deviance
d[1] <- 0 # treatment effect is zero for reference treatment
for (k in 2:nt){ d[k] ~ dnorm(0,.0001) } # vague priors for treatment effects
sd ~ dunif(0,5) # vague prior for between-trial SD.
tau <- pow(sd,-2) # between-trial precision = (1/between-trial variance)
for (j in 1:ns) {
for (k in 1:na[j]) {
time[j,k] ~ dnorm(0, 0.0001) # vague prior for time effect
}
}
#Ranking of treatments#
for(k in 1:nt) {
order[k]<- rank(d[],k)
# this is when the outcome is positive - omit 'nt+1-' when the outcome is negative
most.effective[k]<-equals(order[k],1)
for(j in 1:nt) {
effectiveness[k,j]<- equals(order[k],j)
}
}
for(k in 1:nt) {
for(j in 1:nt) {
cumeffectiveness[k,j]<- sum(effectiveness[k,1:j])
}
}
#SUCRAS#
for(k in 1:nt) {
SUCRA[k]<- sum(cumeffectiveness[k,1:(nt-1)]) /(nt-1)
}
# ranking on relative scale
for (k in 1:nt) {
# rk[k] < - nt+1 -rank(d[],k) # assumes events are “good”
rk[k] < - rank(d[],k) # assumes events are “bad”
best[k] < - equals(rk[k],1) #calculate probability that treat k is best
}
} # *** PROGRAM ENDS
I think the issue might be the following line
log(theta[i,k]) <- mu[i] + delta[i,k] + (time[i,k] - 2) + pow((time[i,k] - 2),2)
When I take out the log in front of theta, my model updates correctly, but when I put it back in it gives me the error. I need this to be a log link model, so I'm not sure what to do.
