Introduction

In this tutorial, you will learn how to perform some simple pair-wise meta-analyses using Bayesian methods. The platform that you will use is R with the JAGS program installed. RStudio is the interface for R which allows you to write your codes and compile them into the console. For this tutorial, you will need to have installed R, RStudio, JAGS, and the various packages associated with JAGS such as rjags, R2jags, mcmcplot, code, and boa.

Example 1 – Aspirin meta-analysis using a fixed effects model

Read-in the data

When you read your data, make sure you are using the correct working directory. You can change your working directory using the setwd command. In this tutorial am using my working directory to allow the codes to work. Make sure that you replace the directory path that I’m using with your directory path when you perform these tutorials for yourself.

I have saved the data as data1.csv. Please use this data for this example. The data can be located here: https://github.com/mbounthavong/Bayesian-meta-analysis

##### Clear data
rm(list=ls())

#### I commented out the system that I am not on. In this case, I commented out the PC system for the Mac.
##### For mac:
setwd("/Users/mbounthavong/Dropbox/UW Courses/PHARM 536 -- Beth's lecture/Beths Winbugs examples")

##### For PC:
# setwd("D:\\mboun\\Dropbox\\UW Courses\\PHARM 536 -- Beth's lecture\\Beths Winbugs examples")

##### Read the data into R.
##### The data is saved as a CSV file.
data = read.csv("data1.csv")
head(data) # Shows the first six entries
##   N          y         V
## 1 1  0.3289011 0.0388957
## 2 2  0.3845458 0.0411673
## 3 3  0.2195622 0.0204915
## 4 4  0.2222206 0.0647646
## 5 5  0.2254672 0.0351996
## 6 6 -0.1246363 0.0096167

Defining values

After you read-in your data, we will define some values that will be used in the JAGS model. These values are necessary for the JAGS model to run. These include the initial values for the priors and the headers for the data columns.

##### Values for simulation
N <- length(data$y) # There are 7 total studies. 

##### load libraries
library(rjags)
library(coda)

##### now prepare dat for JAGS
## N is the number of entries (e.g., 7)
## y is the outcome in the data (e.g., ln(OR))
## V is the variance in the data
dat <- list("N" = N, "y" = data$y, "V" = data$V)  # names list of numbers

##### Initial values
inits <- list( d = 0.0 )

Enter your model

Next, you will need to define your JAGS model using the BUGS syntax. This can be done in two ways:

  1. Using a text editor or

  2. Embed the BUGS model into R.

For this example, I will embed the BUGS model into R. We will save this as “aspirinFE.txt.” The d parameter will be given a normal distribution with mean 0 and precision 1E-5.

##### define JAGS model within R

cat("model
    {

    for ( i in 1:N ) {
    
    P[i] <- 1/V[i]
    
    y[i] ~ dnorm( d, P[i] )

    }

    ### Define the priors
    d ~ dnorm( 0, 0.00001 )

    ### Transform the ln(OR) to OR
    OR <- exp( d )

    }", file="aspirinFE.txt")

Set up the JAGS model

Use the following syntax to set up the JAGS model. Notice that the number of chains is set to “1.” We are going to only use a single chain for our MCMC simulations. The number of adaptions is set at 500. This is done to improve the calibation of the MCMC and improve the mixing of the chain.

#### Set up the JAGS model and settings
jags.m <- jags.model( file = "aspirinFE.txt", data=dat, inits=inits, n.chains=1, n.adapt=500 )
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 7
##    Unobserved stochastic nodes: 1
##    Total graph size: 33
## 
## Initializing model

Specify the parameters to be monitored

These are the parameters that the output will provide for us. Review the BUGS model and notice that the output parameters only include d and OR.

## specify parameters to be monitored
params <- c("d", "OR")

Run JAGS and save the posterior samples

We will save the coda samples as “samps” and iterate this MCMC for 10,000 iterations. The more iterations your run, the longer it will take to complete.

## run JAGS and save posterior samples
samps <- coda.samples( jags.m, params, n.iter=10000 )

Summarize the posteriors samples with a burn-in of 5000 iterations

We need to have a burn-in period to eliminate any poor mixing of the MCMC at the beginning of the simulation. The usual number of burn-ins is approximately 5,000.

## summarize posterior samples
summary(samps)
## 
## Iterations = 1:10000
## Thinning interval = 1 
## Number of chains = 1 
## Sample size per chain = 10000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##      Mean      SD  Naive SE Time-series SE
## OR 1.1153 0.03701 0.0003701      0.0003701
## d  0.1086 0.03318 0.0003318      0.0003318
## 
## 2. Quantiles for each variable:
## 
##       2.5%     25%    50%    75%  97.5%
## OR 1.04453 1.09005 1.1148 1.1398 1.1897
## d  0.04356 0.08622 0.1087 0.1309 0.1737
summary(window(samps, start=5001))  # Burn in of 5000. Start at 5001
## 
## Iterations = 5001:10000
## Thinning interval = 1 
## Number of chains = 1 
## Sample size per chain = 5000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##      Mean      SD  Naive SE Time-series SE
## OR 1.1154 0.03694 0.0005224      0.0005224
## d  0.1087 0.03312 0.0004684      0.0004684
## 
## 2. Quantiles for each variable:
## 
##       2.5%     25%    50%    75%  97.5%
## OR 1.04453 1.09018 1.1153 1.1395 1.1892
## d  0.04356 0.08634 0.1091 0.1306 0.1733

Plot distributions of the parameters

Plot the posterior distributions of the parameters. We are looking for ditributions that are stable and not erratic.

plot(samps)

Example #2 – Aspirin meta-analysis using a random effects model

In this example, we will evaluate the same meta-analysis, but we will use a random effects model instead of a fixed effects model. In random effects model, there is heterogeneity in the outcomes across the different studies. In order to address this, we need to include a new parameter into the meta-analysis called tau.

Read-in the data

I have saved the data as data1.csv. Please use this data for this example.

##### Clear the dataset
rm(list=ls())

##### For mac:
setwd("/Users/mbounthavong/Dropbox/UW Courses/PHARM 536 -- Beth's lecture/Beths Winbugs examples")

##### For PC:
# setwd("D:\\mboun\\Dropbox\\UW Courses\\PHARM 536 -- Beth's lecture\\Beths Winbugs examples")

##### Read the data into R.
##### The data is saved as a CSV file.
data = read.csv("data1.csv")
head(data) # Shows the first six entries
##   N          y         V
## 1 1  0.3289011 0.0388957
## 2 2  0.3845458 0.0411673
## 3 3  0.2195622 0.0204915
## 4 4  0.2222206 0.0647646
## 5 5  0.2254672 0.0351996
## 6 6 -0.1246363 0.0096167

Define the values for the model.

With a random effects model, the prior will need to have initial values. These are usually zero unless you have some idea of what the values should be. For this example, we will set the initial values to zero. Only tau has an initial value of 1. Since there are 7 studies in the meta-analysis, there must be 7 initial values for delta.

##### Values for simulation
N <- length(data$y) # There are 7 total studies. 

##### load libraries
library(rjags)
library(coda)

##### now prepare dat for JAGS
## N is the number of entries (e.g., 7)
## y is the outcome in the data (e.g., ln(OR))
## V is the variance in the data
dat <- list("N" = N, "y" = data$y, "V" = data$V)  # names list of numbers

##### Initial values
inits <- list( d = 0.0, tau = 1, delta = c (0, 0, 0, 0, 0, 0, 0) )

Embed the BUGS model into R.

After setting up the initial values, embed the BUGS model directly into R. We will save this as “aspirinRE2.txt.”

##### define JAGS model within R
## The BUGS model is saved as a text or jag file. I used Notepad ++ to build the model
## but you can use any text editor. Please do not use Words.

#### Alternatively, you can enter the mode directly in R:
cat("model
    {

    for (i in 1:N) {
    
    y[i] ~ dnorm( delta[i], P[i] )

    delta[i] ~ dnorm( d, prec )

    P[i] <- 1/V[i]

    }

    ### Define the priors
    d ~ dnorm( 0, 0.00001 )
    tau ~ dunif( 0, 10 )
    tau.sq <- tau * tau
    prec <- 1 / tau.sq
  
    ### Transform ln(OR) into OR
    OR <- exp ( d )

    }", file="aspirinRE2.txt")

Set up the JAGS model

## Set up the JAGS model.
jags.m <- jags.model( file = "aspirinRE2.txt", data=dat, inits=inits, n.chains=1, n.adapt=500 )
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 7
##    Unobserved stochastic nodes: 9
##    Total graph size: 48
## 
## Initializing model

Specify the parameters of interest.

Normally this is the outcome, such as the odds ratio or standardized mean difference. We also include tau and delta. Delta tells us the treatment effect of ASA. We won’t plot delta[i]; however, you can do this by including delta as part of the parameters to be monitored. Keep in mind that you may run into errors if you plot too many graphics into a single viewing pane due to the margins and size limitations. If you encounter this problem, it maay be best to plot a few parameters at a time.

## specify parameters to be monitored
params <- c("d", "OR", "tau")

Run JAGS model and save the samples

## run JAGS and save posterior samples
samps <- coda.samples( jags.m, params, n.iter=10000 )

Use a burn-in of 5,000 iterations

## summarize posterior samples
summary(samps)
## 
## Iterations = 501:10500
## Thinning interval = 1 
## Number of chains = 1 
## Sample size per chain = 10000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##       Mean      SD  Naive SE Time-series SE
## OR  1.1637 0.11715 0.0011715       0.002366
## d   0.1469 0.09571 0.0009571       0.001840
## tau 0.1626 0.11296 0.0011296       0.006094
## 
## 2. Quantiles for each variable:
## 
##          2.5%     25%   50%    75%  97.5%
## OR   0.981045 1.09357 1.148 1.2160 1.4382
## d   -0.019137 0.08944 0.138 0.1956 0.3634
## tau  0.009806 0.08769 0.144 0.2118 0.4274
summary(window(samps, start=5001))  # Burn in of 5000. Start at 5001.
## 
## Iterations = 5001:10500
## Thinning interval = 1 
## Number of chains = 1 
## Sample size per chain = 5500 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##       Mean      SD Naive SE Time-series SE
## OR  1.1656 0.11952 0.001612       0.003319
## d   0.1485 0.09574 0.001291       0.002548
## tau 0.1580 0.11639 0.001569       0.009341
## 
## 2. Quantiles for each variable:
## 
##          2.5%     25%    50%    75%  97.5%
## OR   0.983968 1.09745 1.1483 1.2158 1.4343
## d   -0.016162 0.09299 0.1383 0.1954 0.3607
## tau  0.004856 0.08291 0.1424 0.2090 0.4012

Plot the posterior distribution

Plot the traces and posterior distributions.

plot(samps)

Example 3 – Meta-analysis using a random effects Binomial model

In this example, we will perform a meta-analysis with a binomial model. The data set will include the number of events of interest and the total number of the sample. The data will be data2.csv

Read in the data

##### Clear the dataset
rm(list=ls())

##### For mac:
setwd("/Users/mbounthavong/Dropbox/UW Courses/PHARM 536 -- Beth's lecture/Beths Winbugs examples")

##### For PC:
# setwd("D:\\mboun\\Dropbox\\UW Courses\\PHARM 536 -- Beth's lecture\\Beths Winbugs examples")

##### Read the data into R.
##### The data is saved as a CSV file.
data <- read.csv("data2.csv")
head(data) # Shows the first six entries
##    rA  rB   nA   nB
## 1  49  67  615  624
## 2  44  64  758  771
## 3 102 126  832  850
## 4  32  38  317  309
## 5  85  52  810  406
## 6 246 219 2267 2257

Load libraries and define the values for the simulation

##### Values for simulation
N <- length(data$rA) # There are 7 total studies. 


##### load libraries
library(rjags)
library(coda)

##### now prepare dat for JAGS
## N is the number of entries (e.g., 7)
## rA is the number of events for Treatment A
## rB is the number of events for Treatment B
## nA is the total number (denominator) for Treatment A
## nB is the total number (denominator) for Treatment B
dat <- list("N" = N, "rA" = data$rA, "rB" = data$rB, "nA" = data$nA, "nB" = data$nB)  # names list of numbers

##### Initial values
inits <- list( d = 0, tau = 1, delta = c (0,0,0,0,0, 0,0), mu = c (0,0,0,0,0, 0,0) )

Embed the BUGS model into R

#### Alternatively, you can enter the model directly into R.
cat("model
    {
    
    ## Define the likelihood (logit) model:
    for( i in 1:N ) {
    
    rA[i] ~ dbin( pA[i], nA[i] )
    rB[i] ~ dbin( pB[i], nB[i] )
    
    logit( pA[i] ) <- mu[i]
    logit( pB[i] ) <- mu[i] + delta[i]
    
    mu[i] ~ dnorm( 0.0, 0.00001 )
    
    delta[i] ~ dnorm( d, prec )
    
    }
    
    ## Define priors
    d ~ dnorm( 0.0, 0.000001 )
    tau ~ dunif( 0, 10 )
    tau.sq <- tau*tau
    prec <- 1/( tau.sq )
    
    ## Outcomes of interest
    OR <- exp( d )
    prob.OR1 <- step( d )
    
    }", file="aspirinREbin2.txt")

Specify the JAGS model

## Set up the JAGS model
jags.m <- jags.model( file = "aspirinREbin2.txt", data=dat, inits=inits, n.chains=1, n.adapt=500 )
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 14
##    Unobserved stochastic nodes: 16
##    Total graph size: 91
## 
## Initializing model

Specify the parameters of interest

## specify parameters to be monitored
params <- c("OR", "prob.OR1", "tau")

Run JAGS and save the posterior samples

## run JAGS and save posterior samples
samps <- coda.samples( jags.m, params, n.iter=10000 )

Summarize the posterior samples

## summarize posterior samples
summary(samps)
## 
## Iterations = 501:10500
## Thinning interval = 1 
## Number of chains = 1 
## Sample size per chain = 10000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##            Mean     SD Naive SE Time-series SE
## OR       1.1648 0.1147 0.001147       0.002797
## prob.OR1 0.9628 0.1893 0.001893       0.002927
## tau      0.1605 0.1095 0.001095       0.004685
## 
## 2. Quantiles for each variable:
## 
##             2.5%    25%    50%    75%  97.5%
## OR       0.97805 1.0947 1.1512 1.2191 1.4283
## prob.OR1 0.00000 1.0000 1.0000 1.0000 1.0000
## tau      0.01749 0.0833 0.1402 0.2109 0.4339
summary(window(samps, start=5001))  # Burn in of 5000. Start at 5001.
## 
## Iterations = 5001:10500
## Thinning interval = 1 
## Number of chains = 1 
## Sample size per chain = 5500 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##            Mean     SD Naive SE Time-series SE
## OR       1.1659 0.1145 0.001544       0.003167
## prob.OR1 0.9627 0.1894 0.002554       0.003773
## tau      0.1653 0.1074 0.001448       0.005788
## 
## 2. Quantiles for each variable:
## 
##             2.5%     25%    50%    75%  97.5%
## OR       0.97505 1.09501 1.1538 1.2208 1.4195
## prob.OR1 0.00000 1.00000 1.0000 1.0000 1.0000
## tau      0.02529 0.09033 0.1435 0.2151 0.4332

Plot the traces and posterior distrubition for the outcomes.

The outcomes are odds ratio and 95% confidence intervals because this is a binomial model.

plot(samps)

Example 4.1 – Sensitivity analysis using different priors (priors set 1)

In the three examples, we started with vague priors. In Bayesian analsyis, priors can be based on our knowledge or they can be “vague.”" Vague priors are uninformative and will yield results similar to a frequentist approach. Since this is a pair-wise comparison of two treatment options, we expect to use vague priors to generate OR and 95% CI similar to that of a frequentist approach.

We are going to perform a meta-analysis by pooling the studies comparing treatment A and B. In this example, there are two treatments (A and B) that result in an event. The number of events are denoted as rA and rB for treatments A and B, respectively. The total number of subjects (denominator) is denoted by nA and nB for treatments A and B, respectively.

Read-in the data

Unlike our previous examples, we will enter the data directly into R. You can still read-in the data using hte read.csv() function; however, you should also know that there are more than one way to enter data into R.

##### Clear the dataset
rm(list=ls())

##### For mac:
setwd("/Users/mbounthavong/Dropbox/UW Courses/PHARM 536 -- Beth's lecture/Beths Winbugs examples")

##### For PC:
# setwd("D:\\mboun\\Dropbox\\UW Courses\\PHARM 536 -- Beth's lecture\\Beths Winbugs examples")

dat <- list(rA = c( 49, 44, 102, 32, 85, 246, 1570 ),
            nA = c( 615, 758, 832, 317, 810, 2267, 8587 ),
            rB = c( 67, 64, 126, 38,52, 219, 1720 ),
            nB = c( 624, 771, 850, 309, 406, 2257, 8600 ),
            Nstud = 7)

Embed the BUGS model into R

Next, we will embed the BUGS model directly into R and save it as “betaMeta1.txt.”

### Embed the BUGS model
cat("model
    {
    for( i in 1 : Nstud ) {

    rA[i] ~ dbin(pA[i], nA[i])
    rB[i] ~ dbin(pB[i], nB[i])

    logit(pA[i]) <- mu[i]
    logit(pB[i]) <- mu[i] + delta[i]

    mu[i] ~ dnorm(0.0,1.0E-5)
    delta[i] ~ dnorm(d, prec)

    }

    OR <- exp(d)
    d ~ dnorm(0.0,1.0E-6)

    tau ~ dunif(0,10)
    tau.sq <- tau*tau
    prec <- 1/(tau.sq)

    }", file="betaMeta1.txt")

Define the initial values for the parameters

The model has several parameters that need to have initial values. We will use initial values that are 0 for mu and delta. Since there are 7 studies, we will need 7 initial values for mu and delta.

inits <- list("d" = 0,
              "tau"=1, 
              "mu" = c( 0, 0, 0, 0, 0,  0, 0 ),
              "delta" = c( 0, 0, 0, 0, 0,  0, 0 )
              )

params <- c("d","tau","OR")

Set up the JAGS model

We set up the JAGS model using the embeded BUGS model entitled “betaMeta1.txt.” We will use only 1 chain and 500 adaptations in order to calibrate the MCMC and optimize the model mixing.

jags.m <- jags.model( file = "betaMeta1.txt", data=dat, inits=inits, n.chains=1, n.adapt=500 )
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 14
##    Unobserved stochastic nodes: 16
##    Total graph size: 90
## 
## Initializing model

Run JAGS and save the posterior samples

## run JAGS and save posterior samples
samps <- coda.samples( jags.m, params, n.iter=10000 )

Summarize the posterior samples

## summarize posterior samples
summary(samps)
## 
## Iterations = 501:10500
## Thinning interval = 1 
## Number of chains = 1 
## Sample size per chain = 10000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##       Mean      SD  Naive SE Time-series SE
## OR  1.1634 0.11698 0.0011698       0.003500
## d   0.1466 0.09634 0.0009634       0.002850
## tau 0.1638 0.11156 0.0011156       0.005475
## 
## 2. Quantiles for each variable:
## 
##          2.5%     25%    50%    75%  97.5%
## OR   0.982228 1.08817 1.1462 1.2180 1.4454
## d   -0.017932 0.08450 0.1365 0.1972 0.3684
## tau  0.005797 0.08992 0.1473 0.2157 0.4419
summary(window(samps, start=5001))  # Burn in of 5000. Start at 5001.
## 
## Iterations = 5001:10500
## Thinning interval = 1 
## Number of chains = 1 
## Sample size per chain = 5500 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##       Mean      SD Naive SE Time-series SE
## OR  1.1640 0.11861 0.001599       0.004554
## d   0.1470 0.09772 0.001318       0.003595
## tau 0.1636 0.11103 0.001497       0.006910
## 
## 2. Quantiles for each variable:
## 
##          2.5%     25%    50%    75%  97.5%
## OR   0.984877 1.08559 1.1463 1.2190 1.4615
## d   -0.015238 0.08213 0.1365 0.1980 0.3795
## tau  0.002002 0.08810 0.1477 0.2187 0.4484

Plot the traces and posterior distributions of the parameters of interest

plot(samps)

Record the results

Record the OR and 95% CI. It should be: OR=1.17; 95%CI: 0.97, 1.44. We will compare our results using a different set of priors.

Example 4.2 – Sensitivity analysis using different priors (priors set 2)

We are going to perform the same sensitivity analysis with a different set of priors.

Read-in the data

Unlike our previous examples, we will enter the data directly into R. You can still read-in the data using hte read.csv() function; however, you should also know that there are more than one way to enter data into R.

##### Clear the dataset
rm(list=ls())

##### For mac:
setwd("/Users/mbounthavong/Dropbox/UW Courses/PHARM 536 -- Beth's lecture/Beths Winbugs examples")

##### For PC:
# setwd("D:\\mboun\\Dropbox\\UW Courses\\PHARM 536 -- Beth's lecture\\Beths Winbugs examples")

dat <- list(rA = c( 49, 44, 102, 32, 85, 246, 1570 ),
            nA = c( 615, 758, 832, 317, 810, 2267, 8587 ),
            rB = c( 67, 64, 126, 38,52, 219, 1720 ),
            nB = c( 624, 771, 850, 309, 406, 2257, 8600 ),
            Nstud = 7)

Embed the BUGS model into R.

Next, we will embed the BUGS model directly into R and save it as “betaMeta2.txt.” Here, is where we can change the priors. Notice that tau, tau.sq, and prec are different in Example 4.2 compared to 4.1.

### Embed the BUGS model
cat("model
    {
    for( i in 1 : Nstud ) {

    rA[i] ~ dbin(pA[i], nA[i])
    rB[i] ~ dbin(pB[i], nB[i])

    logit(pA[i]) <- mu[i]
    logit(pB[i]) <- mu[i] + delta[i]

    mu[i] ~ dnorm(0.0,1.0E-5)
    delta[i] ~ dnorm(d, prec)

    }
    ## Output
    OR <- exp(d)
    d ~ dnorm(0.0,1.0E-6)
    
    ## Priors
    tau ~ dnorm( 0, 10 )
    tau.sq <- tau * tau
    prec <- 1 / ( tau.sq )

    }", file="betaMeta2.txt")

Define the initial values for the parameters.

The model has several parameters that need to have initial values. We will use initial values that are 0 for mu and delta. Since there are 7 studies, we will need 7 initial values for mu and delta.

inits <- list("d" = 0,
              "tau"=1, 
              "mu" = c( 0, 0, 0, 0, 0,  0, 0 ),
              "delta" = c( 0, 0, 0, 0, 0,  0, 0 )
              )

params <- c("d","tau","OR")

Set up the JAGS model

We set up the JAGS model using the embeded BUGS model entitled “betaMeta1.txt.” We will use only 1 chain and 500 adaptation in order to calibrate the MCMC and optimize the model mixing.

jags.m <- jags.model( file = "betaMeta2.txt", data=dat, inits=inits, n.chains=1, n.adapt=500 )
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 14
##    Unobserved stochastic nodes: 16
##    Total graph size: 90
## 
## Initializing model

Run JAGS and save the posterior samples

## run JAGS and save posterior samples
samps <- coda.samples( jags.m, params, n.iter=10000 )

Summarize the posterior samples

## summarize posterior samples
summary(samps)
## 
## Iterations = 501:10500
## Thinning interval = 1 
## Number of chains = 1 
## Sample size per chain = 10000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##         Mean      SD  Naive SE Time-series SE
## OR  1.157751 0.10035 0.0010035       0.002849
## d   0.142864 0.08439 0.0008439       0.002389
## tau 0.004506 0.16410 0.0016410       0.030887
## 
## 2. Quantiles for each variable:
## 
##          2.5%      25%      50%    75%  97.5%
## OR   0.990461  1.09350 1.145809 1.2106 1.3857
## d   -0.009585  0.08939 0.136111 0.1911 0.3262
## tau -0.302553 -0.12278 0.001683 0.1320 0.3024
summary(window(samps, start=5001))  # Burn in of 5000. Start at 5001.
## 
## Iterations = 5001:10500
## Thinning interval = 1 
## Number of chains = 1 
## Sample size per chain = 5500 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##         Mean      SD Naive SE Time-series SE
## OR   1.15853 0.09943 0.001341       0.003496
## d    0.14362 0.08330 0.001123       0.002934
## tau -0.06213 0.14960 0.002017       0.031477
## 
## 2. Quantiles for each variable:
## 
##          2.5%      25%      50%    75%  97.5%
## OR   0.993679  1.09429  1.14616 1.2095 1.3845
## d   -0.006341  0.09011  0.13642 0.1902 0.3253
## tau -0.332253 -0.16106 -0.08682 0.0386 0.2454

Plot the traces and posterior distributions of the parameters of interest

plot(samps)

Record the results from Example 4.2

Record the OR and 95% CI. It should be: OR=1.16; 95% CI: 0.99, 1.40. We will compare our results using a different set of priors.

Example 4.3 – Sensitivity analysis using different priors (priors set 3)

We are going to perform the same sensitivity analysis with a different set of priors. This is the third set of priors.

Read-in the data

Unlike our previous examples, we will enter the data directly into R. You can still read-in the data using the read.csv() function; however, you should also know that there are more than one way to enter data into R.

##### Clear the dataset
rm(list=ls())

##### For mac:
setwd("/Users/mbounthavong/Dropbox/UW Courses/PHARM 536 -- Beth's lecture/Beths Winbugs examples")

##### For PC:
# setwd("D:\\mboun\\Dropbox\\UW Courses\\PHARM 536 -- Beth's lecture\\Beths Winbugs examples")

dat <- list(rA = c( 49, 44, 102, 32, 85, 246, 1570 ),
            nA = c( 615, 758, 832, 317, 810, 2267, 8587 ),
            rB = c( 67, 64, 126, 38,52, 219, 1720 ),
            nB = c( 624, 771, 850, 309, 406, 2257, 8600 ),
            Nstud = 7)

Embed the BUGS model into R

Next, we will embed the BUGS model directly into R and save it as “betaMeta3.txt.” Here, is where we can change the priors. Notice that tau, tau.sq, and prec are different in Example 4.3 compared to 4.1 and 4.2.

## Embed the BUGS model
cat("model
    {
    for( i in 1 : Nstud ) {
    
    rA[i] ~ dbin( pA[i], nA[i] )
    rB[i] ~ dbin( pB[i], nB[i] )
    
    logit( pA[i] ) <- mu[i]
    logit( pB[i] ) <- mu[i] + delta[i]
    
    mu[i] ~ dnorm( 0.0, 1.0E-5 )
    delta[i] ~ dnorm( d, prec )
    
    }
    
    ## Output
    OR <- exp( d )
    d ~ dnorm( 0.0, 1.0E-6 )
    
    ## Priors
    tau <- sqrt( tau.sq )
    tau.sq <- 1 / ( prec )
    prec ~ dgamma( 0.001, 0.001 )
    
    }", file="betaMeta3.txt")

Define the initial values for the parameters

The model has several parameters that need to have initial values. We will use initial values that are 0 for mu and delta. Since there are 7 studies, we will need 7 initial values for mu and delta.

Notice that tau is not among the initial values like in Examples 4.1 and 4.2. In Examples 4.1 and 4.2, tau was a distribution, whereas, in Example 4.3, prec is a distribution. As a result, the initial values should only include the distributions that are defined in the model. In this case, they are d, prec, mu and delta.

inits <- list("d" = 0,
              "prec"=1, 
              "mu" = c( 0, 0, 0, 0, 0,  0, 0 ),
              "delta" = c( 0, 0, 0, 0, 0,  0, 0 )
              )

params <- c("d", "prec", "OR")

Set up the JAGS model

We set up the JAGS model using the embeded BUGS model entitled “betaMeta1.txt.” We will use only 1 chain and 500 adaptation in order to calibrate the MCMC and optimize the model mixing.

jags.m <- jags.model( file = "betaMeta3.txt", data=dat, inits=inits, n.chains=1, n.adapt=500 )
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 14
##    Unobserved stochastic nodes: 16
##    Total graph size: 88
## 
## Initializing model

Run JAGS and save the posterior samples

## run JAGS and save posterior samples
samps <- coda.samples( jags.m, params, n.iter=10000 )

Summarize the posterior samples

## summarize posterior samples
summary(samps)
## 
## Iterations = 501:10500
## Thinning interval = 1 
## Number of chains = 1 
## Sample size per chain = 10000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##          Mean        SD  Naive SE Time-series SE
## OR     1.1497   0.09354 0.0009354       0.002334
## d      0.1363   0.07957 0.0007957       0.001988
## prec 209.1433 372.95637 3.7295637      11.758617
## 
## 2. Quantiles for each variable:
## 
##           2.5%      25%     50%      75%     97.5%
## OR    0.990059  1.08997  1.1398   1.1973    1.3657
## d    -0.009991  0.08615  0.1308   0.1801    0.3116
## prec  9.903849 36.96688 82.5160 206.8525 1279.0830
summary(window(samps, start=5001))  # Burn in of 5000. Start at 5001.
## 
## Iterations = 5001:10500
## Thinning interval = 1 
## Number of chains = 1 
## Sample size per chain = 5500 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##          Mean        SD Naive SE Time-series SE
## OR     1.1474   0.09285 0.001252       0.003224
## d      0.1343   0.07939 0.001071       0.002754
## prec 227.7257 414.64235 5.591036      16.633727
## 
## 2. Quantiles for each variable:
## 
##         2.5%     25%     50%      75%     97.5%
## OR    0.9879  1.0894  1.1378   1.1962    1.3584
## d    -0.0122  0.0856  0.1291   0.1792    0.3063
## prec 10.3397 38.9423 85.9871 222.5854 1397.1768

Plot the traces and posterior distributions of the parameters of interest

plot(samps)

Record the results

Record the OR and 95% CI. It should be: OR=1.15; 95% CI: 0.997, 1.37.

Compare results from the other examples

Example 1: OR=1.17; 95% CI: 0.97, 1.44

Example 2: OR=1.16; 95% CI: 0.99, 1.40

Example 3: OR=1.15; 95% CI: 0.997, 1.37

The OR and 95% CI are not that different. However, we do report slight differences in the mean OR and 95% CI. Although the conclusions do not change, the use of different priors illustrate that there are still variations in the results. Despite these variations, the result do not change significantly.

Conclusions

Bayesian analysis is useful when performing meta-analysis because they allow us to perform indirect treatment comparisons. These examples only focused on pair-wise comparisons. In the next lesson, you will be introduced to mixed treatment comparison using Bayesian methods, which will allow us to use a combination of different pair-wise studies to make direct and indirect comparisons.