Hierarchical log odds model example

I am working through Bayesian Approaches to Clinical Trials and Health-Care Evaluation (David J. Spiegelhalter, Keith R. Abrams, Jonathan P. Myles) (referred to as BACTHCE from here on). In chapter three I saw an example (3.13) where I wanted to do the calculations myself. The example concerns a hierarchical model of death rates which is calculated via a normal approximation of odds ratios.

Data 

According to BACTHCE: 'Epidemiology, animal models and biochemical studies suggested intravenous magnesium sulphate may have a protective effect after acute myocardial infarction (AMI), articularly through preventing serious arrhythmias'. A small set of randomized trials was performed. We will analyze those data. In the data below MgD and ControlD refers to deaths for Mg treated and control patients, while MgN and ControlN refers to the number of patients enrolled in each group.
r1 <- read.csv(text="
        Trial      MgD  MgN ControlD ControlN
        Morton       1   40   2   36 
        Rasmussen    9  135  23  135
        Smith        2  200   7  200
        Abraham      1   48   1   46
        Feldstedt   10  150   8  148
        Shechter     1   59   9   56
        Ceremuzynski 1   25   3   23
        LIMIT-2     90 1159 118 1157
        ",skip=1,sep='')

Analysis

Odds rations

The first analysis will be to reproduce all the derived variables, such as log odds ratios. This will reproduce the essentials of table 3.8. All formulas are in BACTHCE. Equivalent N is the N if this were normal distributed data with SD=2.

logoddrat <- function(

    a1,b1=N1-a1,N1=a1+b1,

    a2,b2=N2-a2,N2=a2+b2) {

  log(((a1+.5)/(b1+.5))/

          ((a2+.5)/(b2+.5)))

}

sdlogoddrat <- function(

    a1,b1=N1-a1,N1=a1+b1,

    a2,b2=N2-a2,N2=a2+b2) {

  sqrt(1/(a1+.5)+1/(b1+.5) +1

          /(a2+.5)+1/(b2+.5))

}

r1$logoddsratio <- with(r1,

    logoddrat(a1=MgD,N1=MgN,a2=ControlD,N2=ControlN))

r1$sdodds <- with(r1,

    sdlogoddrat(a1=MgD,N1=MgN,a2=ControlD,N2=ControlN))

r1$equivN <- (2/r1$sdodds)^2

r1

         Trial MgD  MgN ControlD ControlN logoddsratio    sdodds     equivN

1       Morton   1   40        2       36  -0.64616697 1.0587581   3.568342

2    Rasmussen   9  135       23      135  -1.02299771 0.4057220  24.299805

3        Smith   2  200        7      200  -1.12412388 0.7372510   7.359177

4      Abraham   1   48        1       46  -0.04301739 1.1731854   2.906208

5    Feldstedt  10  150        8      148   0.21130909 0.4765711  17.611833

6     Shechter   1   59        9       56  -2.05412373 0.9000425   4.937805

7 Ceremuzynski   1   25        3       23  -1.02554609 1.0207731   3.838854

8      LIMIT-2  90 1159      118     1157  -0.29801453 0.1462380 187.042101

Odds ratios in plot

The plot shows the odds ratios as function of standard deviation of the hyper priors (tau). Mu is the overall mean. Note that I could not find formula for the overall mean in BACTHCE, and just assumed that the value for By in shrinkage was applied. 

mut <- function(tau,r1) {  

  yy <- r1$logoddsratio

  By <- r1$sdodds^2/(r1$sdodds^2+tau^2)

  mus <- sum(yy*r1$equivN*By)/sum(r1$equivN*By)

  yy <- r1$logoddsratio*(1-By)+mus*By

  data.frame(Trial=c(as.character(r1$Trial),'Mu'),tau=tau,mu=c(yy,mus)) 

}

la <- lapply(seq(0,2,.01),function(x) mut(x,r1))

lac <- do.call(rbind,la) 

lac$Trial <- factor(lac$Trial,levels=c(levels(r1$Trial),'Mu'))

ggplot(lac,aes(x=tau,y=mu,col=Trial)) +

    geom_line() +

    theme(legend.position="bottom")  +

    guides(col=guide_legend(ncol=3)) +
    scale_y_continuous('Log odds ratio')

Beta binomial hierarchical model

I wanted to know how that would compare if I did that calculation as a beta-binomial, for which I used JAGS. Lor is the estimated log odds ratio, which with -0.7 is close to what is found in BACTHCE using a method of moments estimator for tau. It can be seen that tau in the neighborhood of 0 is possible. The theta parameter represents the probability of death in both groups. What I noticed here is that the control group has a slightly larger standard deviation in theta. Is that a reason for heterogeneity in tau which was mentioned in BACTHCE?

library(R2jags)

datain <- as.list(r1[,2:5])

datain$n <- nrow(r1)

model1 <- function() {

  for (i in 1:2) {

    kappa_log[i] ~  dpar(1,1.5)

    kappa[i] <- exp(kappa_log[i])

    theta[i] ~ dnorm(thetam[i],thetatau[i])

    thetam[i] ~ dunif(0,1)

    thetatau[i] <- pow(thetasd[i],-2) 

    thetasd[i] ~ dnorm(0,.3) %_% I(0,) 

    a[i] <- kappa[i] * theta[i]

    b[i] <- kappa[i] * (1 - theta[i]) 

  } 

  for (i in 1:n) {

    MgD[i] ~ dbetabin(a[1],b[1],MgN[i])

    ControlD[i] ~ dbetabin(a[2],b[2],ControlN[i])

  }

  lor <- log((theta[1]/(1-theta[1]))/

          (theta[2]/(1-theta[2])))

}

parameters <- c('theta','lor')

jj <- jags(model.file=model1,data=datain,parameters=parameters,DIC=FALSE)

jj 

Inference for Bugs model at "C:/Users/Kees/AppData/Local/Temp/Rtmp6h4pKI/model1814317e3684.txt", fit using jags,

 3 chains, each with 2000 iterations (first 1000 discarded)

 n.sims = 3000 iterations saved

         mu.vect sd.vect   2.5%    25%    50%    75% 97.5%  Rhat n.eff

lor       -0.701   0.485 -1.637 -1.019 -0.706 -0.388 0.298 1.001  2500

theta[1]   0.057   0.021  0.029  0.043  0.053  0.065 0.116 1.003  1200

theta[2]   0.106   0.031  0.060  0.085  0.101  0.123 0.183 1.001  3000

For each parameter, n.eff is a crude measure of effective sample size,

and Rhat is the potential scale reduction factor (at convergence, Rhat=1).

Beta binomial in SAS

In clinical trials, SAS is often the tool of choice. Hence the beta binomial in PROC MCMC. The results are a bit different from JAGS, but then I had some difficulty creating the model in PROC MCMC and reruns gave somewhat varying results. I already did increase the number of samples, but there might be better tricks to get stable results.

data r1;

input  Trial $ MgD  MgN ControlD ControlN;

cards;

        Morton       1   40   2   36 

        Rasmussen    9  135  23  135

        Smith        2  200   7  200

        Abraham      1   48   1   46

        Feldstedt   10  150   8  148

        Shechter     1   59   9   56

        Ceremuzynski 1   25   3   23

        LIMIT-2     90 1159 118 1157

;;;;

run;

proc mcmc data=r1 monitor=(theta1 theta2 lor) 

     ntu=2000 nbi=4000 nmc=8000 nthin=8;

      parms (kappa_log1 kappa_log2) 1.8 (theta1 theta2) .1;

      prior kappa_log1 ~ uniform(0,5);

      prior kappa_log2 ~ uniform(0,5);

      prior theta1 ~ uniform(0,1);

      prior theta2 ~ uniform(0,1);

      kappa1 = exp(kappa_log1);

      kappa2 = exp(kappa_log2);

      a1 = kappa1*theta1;

      b1 = kappa1*(1-theta1);

      a2 = kappa2*theta2;

      b2 = kappa2*(1-theta2);

      lor = log((a1/b1)/(a2/b2));

      random Mgp ~beta(a1,b1) subject=trial;

      random Controlp ~beta(a2,b2) subject=trial;

      model MgD ~  binomial(MgN,Mgp);

      model ControlD ~ binomial(ControlN,Controlp); 

 run;  

Posterior Summaries and Intervals
Parameter	N	Mean	Standard Deviation	95% HPD Interval
theta1	1000	0.0519	0.0157	0.0244	0.0841
theta2	1000	0.0959	0.0224	0.0596	0.1404
lor	1000	-0.6790	0.3880	-1.4862	0.0872

Effective Sample Sizes
Parameter	ESS	Autocorrelation Time	Efficiency
theta1	591.5	1.6907	0.5915
theta2	247.0	4.0483	0.2470
lor	343.7	2.9097	0.3437