It’s much much easier if you just post the text of your programs as code:
# Poisson likelihood, log link
# 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
r[i,k] ~ dpois(theta[i,k]) # Poisson likelihood
theta[i,k] <- lambda[i,k]*E[i,k] # failure rate * exposure
log(lambda[i,k]) <- mu[i] + delta[i,k] # model for linear predictor
#Deviance contribution
dev[i,k] <- 2*((theta[i,k]-r[i,k]) + r[i,k]*log(r[i,k]/theta[i,k])) }
# summed residual deviance contribution for this trial
resdev[i] <- sum(dev[i,1:na[i]])
for (k in 2:na[i]) { # LOOP THROUGH ARMS
# trial-specific LOR distributions
delta[i,k] ~ dnorm(md[i,k],taud[i,k])
# mean of LOR distributions (with multi-arm trial correction)
md[i,k] <- d[t[i,k]] - d[t[i,1]] + sw[i,k]
# precision of LOR distributions (with multi-arm trial correction)
taud[i,k] <- tau *2*(k-1)/k
# adjustment for multi-arm RCTs
w[i,k] <- (delta[i,k] - d[t[i,k]] + d[t[i,1]])
# cumulative adjustment for multi-arm trials
sw[i,k] <- sum(w[i,1:k-1])/(k-1)
}
}
totresdev <- sum(resdev[]) #Total Residual Deviance
d[1]<-0 # treatment effect is zero for reference treatment
# vague priors for treatment effects
for (k in 2:nt){ d[k] ~ dnorm(0,.0001) }
sd ~ dunif(0,5) # vague prior for between-trial SD
tau <- pow(sd,-2) # between-trial precision = (1/between-trial variance)
# Provide estimates of treatment effects T[k] on the natural (rate) scale
# Given a Mean Effect, meanA, for 'standard' treatment A,
# with precision (1/variance) precA
A ~ dnorm(meanA,precA)
for (k in 1:nt) { log(T[k]) <- A + d[k] }
} # *** PROGRAM ENDS
data {
int<lower=0> ns; //number of studies
int<lower=0> nt; //number of treatments
int<lower=1> no; //number of observations
int<lower=1> na[ns];
int<lower=0> s[no]; //study indicator
int<lower=0> t[no]; //treatment indicator
int<lower=0> b[no]; //baseline indicator
int<lower=0> r[no]; //number of 'success' in all arms
int<lower=0> n[no]; //number of trials in each arm
real e[no]; //number of events per person-years observed
int<lower=1> b_ndx[ns];
int<lower=1> t_ndx[no-ns];
}
parameters {
vector[nt-1] d_param; //mean for treatment effect distribution
vector[no-ns] delta_param;
vector[ns] mu; //study mean
real<lower=0, upper=5> stdev; //standard deviation for treatment effect distrution
vector[nt] A; //for absolute effects
}
transformed parameters{
vector[nt] d;
vector[no] delta;
vector[no] lambda;
vector[no] theta;
vector[no] md;
cov_matrix[max(na)] cov;
d[1] <- 0;
for(i in 2:nt) d[i] <- d_param[i-1];
for (i in 1:ns) delta[b_ndx[i]] <- 0;
for (i in 1:(no-ns)) delta[t_ndx[i]] <- delta_param[i];
lambda <- exp(mu[s] + delta); //model for linear predictor
theta <- lambda*e[no]; // event rate * exposure
md <- d[t] - d[b];
#COVARIANCE MATRIX
for(i in 1:max(na)){
cov[i,i] <- stdev^2;
for(j in i+1:max(na)){
cov[i,j]<-(stdev^2)/2;
cov[j,i]<-cov[i,j];
}
}
}
model {
A ~ normal(-3, 0.751646); // liklihood for placebo effect; precision (1.77); SD = 1/sqrt(prec) = 0.751646
//mu ~ normal(0, 10000);
r ~ poisson(theta); //poisson likelihood
{
int pos_delta;
int pos_md;
pos_delta <- 1;
pos_md <- 2;
for (i in 1:ns) {
if (na[i]==2) {
delta_param[pos_delta] ~ normal(md[pos_md],stdev);
} else {
delta_param[pos_delta:(pos_delta+na[i]-2)] ~ multi_normal(md[pos_md:(pos_md+na[i]-2)], cov[1:(na[i]-1), 1:(na[i]-1)]);
}
pos_delta <- pos_delta + na[i] - 1;
pos_md <- pos_md + na[i];
}
}
}
generated quantities {
vector[no] dev;
real totalresdev; //create scalar variable for total residual deviance
vector[nt] T; //create scalar variables for absolute effects
for(i in 1:no) dev[i] <- 2*((theta[i]-r[i]) + r[i]*log(r[i]/theta[i]));
totalresdev <- sum(dev);
T <- exp(A + d); //absolute effect of treatment
}
I see things ike this in the Stan program:
A ~ normal(-3, 0.751646);
and this in the BUGS:
A ~ dnorm(meanA,precA)
It makes it very hard to compare.
That’s a lot of code to digest and compare, so if you could perhaps isolate where you think there may be a difference, we could take a look.
This may not vary from the BUGS, but we strongly disrecommend these uniform interval priors in Stan as they don’t do what most people think they’re going to do:
real<lower=0, upper=5> stdev;
There’s a lot you can do to make that Stan program more efficient.
I’d strongly recommend indenting properly rather than having things like this that don’t scan:
for (i in 1:ns) {
if (na[i]==2) {
delta_param[pos_delta] ~ normal(md[pos_md],stdev);
} else {
delta_param[pos_delta:(pos_delta+na[i]-2)] ~ multi_normal(md[pos_md:(pos_md+na[i]-2)], cov[1:(na[i]-1), 1:(na[i]-1)]);
}
pos_delta <- pos_delta + na[i] - 1;
pos_md <- pos_md + na[i];
}
All of thos pos mods are in the loop body.
I didn’t understand this:
lambda <- exp(mu[s] + delta); //model for linear predictor
theta <- lambda*e[no]; // event rate * exposure
...
r ~ poisson(theta); //poisson likelihood
Why does only the last index of e get used?
Stan has a Poisson on the log scale, rather than
r ~ poisson(exp(...) * u);
you can use
r ~ poisson_log(... + log(u));
which will be much more artihmetically stable.