karolinskaexercises

Karolinska Methods Week: meta-analysis examples and exercises

This document provides code to run all the tutorial examples, plus some others for you to try.

library("metafor")

## Loading required package: Matrix

## Loading required package: metadat

## 
## Loading the 'metafor' package (version 3.4-0). For an
## introduction to the package please type: help(metafor)

Slide 20-21

Data input and summaries for the Zinc example

cat(file="zincdata.txt", "study n1 mean1 sd1 n2 mean2 sd2  \n 
\"Douglas 1987\"  33 12.10 9.80  30 7.70 9.80 \n
\"Petrus 1998\"  52  4.40 1.40  49 5.10 2.80 \n
\"Prasad 2000\"  25  4.50 1.60  23 8.10 1.80 \n
\"Prasad 2008\"  25  4.00 1.04  25 7.12 1.26 \n
\"Turner 2000b\"  68  6.89 3.35  71 7.55 3.96 \n
\"Turner 2000c\"  72  7.90 4.25  71 7.55 3.96")

zinc <- read.table("zincdata.txt", header=TRUE, sep="")
zinc[2:7] <- lapply( zinc[2:7], as.numeric)
zinc <- as.data.frame(zinc, 
    col.names=c("study","n1","mean1","sd1", "n2","mean2","sd2"))

zinc$betahat <- with(zinc, mean2-mean1)
zinc$se      <- with(zinc, sqrt( sd1^2/n1 + sd2^2/n2 ) )

# plot these, on their own
forest(rma(yi=betahat, sei=se, data=zinc), slab=zinc$study, header=TRUE, xlim=c(-27, 12), ylim=c(0.5,8.5),
    ilab=zinc[,2:7], ilab.xpos=seq(-20,-11, l=6))
text( x=seq(-20,-11, l=6), y=rep(6.7, 6), adj=c(0.5, 0), c(
    expression(n[C]), expression(hat(mu)[C]), expression("SD"[C]),
    expression(n[T]), expression(hat(mu)[T]), expression("SD"[T])))

Slide 29-30

Overall estimate (fixed effects)

zinc.meta <- with( zinc, c(
    estimate = sum( betahat/se^2 )/sum( 1/se^2 ),
    se = sqrt( 1/sum( 1/se^2 ) ),
    zval = sum( betahat/se^2 )/sqrt(sum( 1/se^2 )),
    pval = pchisq( sum( betahat/se^2 )^2/sum( 1/se^2 ), df=1, lower=FALSE),
    ci.lb = sum( betahat/se^2 )/sum( 1/se^2 ) - qnorm(0.975)*sqrt( 1/sum( 1/se^2 ) ),
    ci.ub = sum( betahat/se^2 )/sum( 1/se^2 ) + qnorm(0.975)*sqrt( 1/sum( 1/se^2 ) )
))
round(zinc.meta, 4)

## estimate       se     zval     pval    ci.lb    ci.ub 
##   2.0423   0.2067   9.8814   0.0000   1.6372   2.4474

# a rather easier way to do it
library("metafor")
rma1 <- rma( yi= betahat, sei=se, data=zinc, method="FE")

# get rid of all the junk
round(unlist(summary(rma1)[2:7]),4)

##   beta     se   zval   pval  ci.lb  ci.ub 
## 2.0423 0.2067 9.8814 0.0000 1.6372 2.4474

Slide 31-32

The genetic association example

genet <- data.frame(
    names   =c("AGES","ARIC","CHS","FHS","RS","RES"),
    betahat =c( -0.37, -0.23, -0.31, -0.79, -0.13, -0.96), 
    stderr  =c(  0.26,  0.16,  0.33,  0.16,  0.25,  0.39),
    n       =c(  3219,  8047,  3277,  8096,  4737,  1760),
    perc.fem=c(   58,    53,     61,   54,    60,    56),
    av.age  =c(    51,    54,     72,   38,    68,    64)) 

# plot these, on their own

rma2 <- rma( yi=betahat, sei=stderr, data=genet, method="FE")
rma2

## 
## Fixed-Effects Model (k = 6)
## 
## I^2 (total heterogeneity / total variability):   50.14%
## H^2 (total variability / sampling variability):  2.01
## 
## Test for Heterogeneity:
## Q(df = 5) = 10.0277, p-val = 0.0745
## 
## Model Results:
## 
## estimate      se     zval    pval    ci.lb    ci.ub     ​ 
##  -0.4536  0.0896  -5.0651  <.0001  -0.6292  -0.2781  *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# slide 32

forest(rma1, slab=zinc$study, main="Fixed effects meta-analysis", header=TRUE)

Exercise 1

Implement fixed-effects meta-analysis for the following three examples

The first analysis includes four studies where persons who donated a kidney were compared with persons in a control group. The outcome was Systolic Blood Pressure (mmHg). The effect size was the difference in mean Systolic Blood Pressure in kidney donaters, versus control.

cat(file="kiddonor.txt", "study,    donor.mean, donor.sd,   donor.n, ctrl.mean, ctrl.sd,    ctrl.n 
Najarian 1992,  134, 15,    57, 130,    21, 50 
Undurraga 1998, 125, 18,    30, 118,    13, 30 
Talselth 1986,  140, 23,    32, 132,    29, 32 
Williams 1886,  136, 25,    38, 129,    16, 16
")
kiddonor <- read.table("kiddonor.txt", sep=",", header=TRUE)

# make the estimate and corresponding standard errors
kiddonor$estimate <- with(kiddonor, donor.mean-ctrl.mean)
kiddonor$se       <- with(kiddonor, sqrt(donor.sd^2/donor.n + ctrl.sd^2/ctrl.n) )

# Add your code here

The second analysis includes nine studies where patients were randomized to be treated for cardiac arrest based on either old guidelines or new guidelines. Outcome was survival (being alive), and we focus on the logs odds ratio as the effect size.

cat(file="cardiac.txt", 
"study, new.alive, new.n, old.alive, old.n
Steinmetz 2008, 41, 226,    21, 193
Sayre 2009, 96, 1021,   40, 660
Olasveengen 2009,   63, 482,    46, 435
Robinson 2010,  20, 170,    18, 162
Hinchey 2010,   47, 410,    18, 425
Hung 2010,  30, 430,    47, 463
Aufderheide 2010,   211,    1605,   166,    1641
Bigham 2011,    177,    2725,   294,    5054
Lick 2011,  48, 247,    9,  106")
file.show("cardiac.txt")
cardiac <- read.table("cardiac.txt", header=TRUE, sep=",")

# computing study-specific log odds ratios and corresponding standard errors

cardiac$odds.new <- with(cardiac, new.alive/(new.n-new.alive))
cardiac$odds.old <- with(cardiac, old.alive/(old.n-old.alive))
cardiac$logOR <- with(cardiac, log( odds.new/odds.old))
cardiac$se <- with(cardiac, 1/new.alive + 1/(new.n-new.alive) + 1/old.alive + 1/(old.n-old.alive))

# Add your code here

The third analysis includes eleven studies where mothers whose children suffered with chronic illnesses were evaluated for post-traumatic stress disorder (PTSD). The outcome was the proportion of mothers who showed symptoms of PTSD.

cat(file="ptsd.txt", 
"study, events, n
Pelcovitz 1996, 8,  37
Barakat 1997,   36, 400
Manne 1998, 5,  100
Fuemmeler 2001, 12, 30
Libov 2002, 15, 80
Brown 2003, 20, 85
Kazak 2004, 26, 200
Manne 2004, 13, 120
Kazak 2005, 48, 200
Landolt 2005,   15, 70
Glover, 31, 100")
file.show("ptsd.txt")
ptsd <- read.table("ptsd.txt", header=TRUE, sep=",")

# simple binomial estimates of rates, and estimated standard errors
ptsd$rate <- with(ptsd, events/n)
ptsd$rate.se <- with(ptsd, sqrt(rate*(1-rate)/n))

# estimate of log odds, with estimated standard errors
# - computed using glm()

for(i in 1:11){
  glmi <- glm(matrix(c(ptsd$events[i], ptsd$n[i]-ptsd$events[i]),1,2)~1, family="binomial")
  ptsd$logodds[i] <- coef(glmi)[1]
  ptsd$logodds.se[i] <- coef(summary(glmi))[1,2]
}

#Add your code here

Slide 36-38

Random effects analyses, and comparing it to fixed-effects

zinc.tau2 <- with(zinc, 
    max(0, (sum( (betahat-zinc.meta["estimate"])^2/se^2 ) - nrow(zinc)+ 1)/(
        sum(se^-2) - sum(se^-4)/sum(se^-2) ) ) )

zinc.meta.r <- with( zinc, c(
    estimate = sum( betahat/(se^2+zinc.tau2) )/sum( 1/(se^2+zinc.tau2) ),
    se = sqrt( 1/sum( 1/(se^2+zinc.tau2) ) ),
    zval = sum( betahat/(se^2+zinc.tau2) )/sqrt(sum( 1/(se^2+zinc.tau2) )),
    pval = pchisq( sum( betahat/(se^2+zinc.tau2) )^2/sum( 1/(se^2+zinc.tau2) ), df=1, lower=FALSE),
    ci.lb = sum( betahat/(se^2+zinc.tau2) )/sum( 1/(se^2+zinc.tau2) ) - qnorm(0.975)*sqrt( 1/sum( 1/(se^2+zinc.tau2) ) ),
    ci.ub = sum( betahat/(se^2+zinc.tau2) )/sum( 1/(se^2+zinc.tau2) ) + qnorm(0.975)*sqrt( 1/sum( 1/(se^2+zinc.tau2) ) )
))

round(zinc.meta.r, 4)

## estimate       se     zval     pval    ci.lb    ci.ub 
##   1.2051   0.7600   1.5856   0.1128  -0.2845   2.6946

# or the quick way

rma1.r <- rma( yi= betahat, sei=se, data=zinc, method="DL")
rma1.r

## 
## Random-Effects Model (k = 6; tau^2 estimator: DL)
## 
## tau^2 (estimated amount of total heterogeneity): 2.8121 (SE = 2.3063)
## tau (square root of estimated tau^2 value):      1.6769
## I^2 (total heterogeneity / total variability):   90.71%
## H^2 (total variability / sampling variability):  10.77
## 
## Test for Heterogeneity:
## Q(df = 5) = 53.8406, p-val < .0001
## 
## Model Results:
## 
## estimate      se    zval    pval    ci.lb   ci.ub   ​ 
##   1.2051  0.7600  1.5856  0.1128  -0.2845  2.6946    
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

rma2   <- rma( yi=betahat, sei=stderr, data=genet, method="FE")
rma2.r <- rma( yi=betahat, sei=stderr, data=genet, method="DL")
rma2

## 
## Fixed-Effects Model (k = 6)
## 
## I^2 (total heterogeneity / total variability):   50.14%
## H^2 (total variability / sampling variability):  2.01
## 
## Test for Heterogeneity:
## Q(df = 5) = 10.0277, p-val = 0.0745
## 
## Model Results:
## 
## estimate      se     zval    pval    ci.lb    ci.ub     ​ 
##  -0.4536  0.0896  -5.0651  <.0001  -0.6292  -0.2781  *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

rma2.r

## 
## Random-Effects Model (k = 6; tau^2 estimator: DL)
## 
## tau^2 (estimated amount of total heterogeneity): 0.0527 (SE = 0.0694)
## tau (square root of estimated tau^2 value):      0.2296
## I^2 (total heterogeneity / total variability):   50.14%
## H^2 (total variability / sampling variability):  2.01
## 
## Test for Heterogeneity:
## Q(df = 5) = 10.0277, p-val = 0.0745
## 
## Model Results:
## 
## estimate      se     zval    pval    ci.lb    ci.ub    ​ 
##  -0.4447  0.1366  -3.2553  0.0011  -0.7124  -0.1769  ** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

rma2.r$se^2 / rma2$se^2

## [1] 2.326205

Exercise 2

Compute random effects meta-analysis, using the DerSimonian Laird estimate, for one of the earlier examples. Compare the results

Slides 40-43

Plots to assess heterogeneity

forest(rma2, slab=genet$names, showweights=TRUE, header=TRUE, main="Using forest() defaulrs")

forest(rma2, slab=genet$names, header=TRUE, xlim=c(-4.0,2), xlab="Observed association",
    ilab=genet[,c("n", "betahat","stderr")], ilab.xpos=c(-3.0,-2.5,-2))
text( x=c(-3,-2.5,-2), y=rep(7.7, 3), adj=c(0.5, 0), main="With extra arguments",
c(expression(n[i]), expression(hat(beta)[i]), expression(sigma[i])))

rma2.r <- rma( yi=betahat, sei=stderr, data=genet, method="DL")
forest(rma2, slab=genet$names, showweights=TRUE, header=TRUE, ylim=c(-2.5, 8.5), 
       main="adding RE with addpoly()")
addpoly(rma2.r)

rma2.r$se^2/rma2$se^2 # eek! A factor of 2.3 in sample size is real money!

## [1] 2.326205

forest(rma1.r, slab=zinc$study, main="Fixed & random effects meta-analysis with prediction", 
header=TRUE, ylim=c(-2.5,9), addpred=TRUE)
addpoly(rma1)

funnel(rma2, main="Funnel plot, metafor")

# other versions (not packaged up)
plot( betahat~I(1/stderr^2), data=genet, 
    xlab=expression(1/sigma[i]^2), ylab=expression(hat(beta)[i]),
    main="Alt Funnel plot, hand-coded", ylim=c(-1.5,0.5), xlim=c(2,45), las=1)
abline(h=rma2["b"][[1]][1,1], lty=2, col="grey50")
curve( rma2["b"][[1]][1,1] + qnorm(0.975)/sqrt(x),0,45, add=TRUE, col=2)
curve( rma2["b"][[1]][1,1] - qnorm(0.975)/sqrt(x),0,45, add=TRUE, col=2)
legend("topright", lty=1, bty="n", col=2, "pointwise 95% bounds\n under homogeneity")
text(betahat~I(1/stderr^2), data=genet, genet$names, pos=c(1,1,1,1,3,1))

# preferable because x-axis is information

# the raincloud plot
plot( I(1/stderr^2/sum(1/stderr^2))~betahat, data=genet,
    ylab="inverse-variance weight", xlab=expression(hat(beta)[i]),
    main="Raincloud plot, genetics example", ylim=c(0,1), xlim=c(-1.75,0.5))
with(genet, segments(betahat-1.96*stderr, 1/stderr^2/sum(1/stderr^2), betahat+1.96*stderr, 1/stderr^2/sum(1/stderr^2)))
points( x= rma2["b"][[1]][1,1] , y=1, pch=19, cex=2)

# add the RE version
points( I(1/(stderr^2+rma2.r["tau2"][[1]][1])/sum(1/(stderr^2+rma2.r["tau2"][[1]][1])))~betahat, 
data=genet, col=2)
points( x= rma2.r["b"][[1]][1,1] , y=1, pch=1, cex=2, col=2)
mtext(side=3, "RE weights in red", col=2)

# the Galbraith plot
par(mar=c(4,4,0,0))
radial(rma2, zlab=expression(frac(hat(beta)[i],sigma[i])), las=1, xlab=expression("1/"*sigma[i]))

Exercise 3

Implement two plots showing heterogeneity for one of the earlier examples.

Slide 46-52

Calculating measures of heteogeneity

zinc.meta.het <- with( zinc, c(
    Q = sum( (betahat-zinc.meta["estimate"])^2/se^2 ),
    I2 = max(0,sum( (betahat-zinc.meta["estimate"])^2/se^2 )-(nrow(zinc)-1))/sum( (betahat-zinc.meta["estimate"])^2/se^2 ),
    H2 = sum( (betahat-zinc.meta["estimate"])^2/se^2 )/(nrow(zinc)-1),
    pval = pchisq( sum( (betahat-zinc.meta["estimate"])^2/se^2 ), df=nrow(zinc)-1, lower=FALSE)
    ))

round(zinc.meta.het,4)

##       Q      I2      H2    pval 
## 53.8406  0.9071 10.7681  0.0000

# compare these with the quick way
rma1

## 
## Fixed-Effects Model (k = 6)
## 
## I^2 (total heterogeneity / total variability):   90.71%
## H^2 (total variability / sampling variability):  10.77
## 
## Test for Heterogeneity:
## Q(df = 5) = 53.8406, p-val < .0001
## 
## Model Results:
## 
## estimate      se    zval    pval   ci.lb   ci.ub     ​ 
##   2.0423  0.2067  9.8814  <.0001  1.6372  2.4474  *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

confint(rma1.r)

## 
##        estimate   ci.lb    ci.ub 
## tau^2    2.8121  1.1862  46.0152 
## tau      1.6769  1.0891   6.7835 
## I^2(%)  90.7133 80.4702  99.3783 
## H^2     10.7681  5.1204 160.8371

# tackling the interval under fixed effects
zinc.meta.het["zeta2"] <- (zinc.meta.het["Q"]-nrow(zinc)+1)/sum(zinc$se^-2)
sqrt(zinc.meta.het["zeta2"])

##    zeta2 
## 1.444411

f1 <- function(x){ qchisq( 0.025, df=nrow(zinc)-1, ncp=x)}; f1 <- Vectorize(f1)
f2 <- function(x){ qchisq( 0.975, df=nrow(zinc)-1, ncp=x)}; f2 <- Vectorize(f2)
hici <- uniroot( function(ncp){ 
    qchisq( 0.025, df=nrow(zinc)-1, ncp=ncp)-zinc.meta.het["Q"]},
    c(0,100) )$root
loci <- uniroot( function(ncp){ 
    qchisq( 0.975, df=nrow(zinc)-1, ncp=ncp)-zinc.meta.het["Q"]},
    c(0,100) )$root
# giving intervals for zeta^2, and zeta:
c(est=zinc.meta.het["zeta2"], loci=loci/sum(zinc$se^-2), hici=hici/sum(zinc$se^-2))

## est.zeta2      loci      hici 
##  2.086323  1.089316  3.500136

sqrt(c(est=zinc.meta.het["zeta2"], loci=loci/sum(zinc$se^-2), hici=hici/sum(zinc$se^-2)))

## est.zeta2      loci      hici 
##  1.444411  1.043703  1.870865

curve(f2 , 0,100, ylab="Upper, lower 2.5% quantiles", xlab=expression("noncentrality parameter = "*zeta^2*sum(sigma[i]^"-2"))
, ylim=c(0,150),
main=expression("Quantiles of noncental "*chi["k-1"]^2))
curve(f1, 0,100, add=TRUE) 
abline(h=zinc.meta.het["Q"], col=2)

confint(rma1.r)

## 
##        estimate   ci.lb    ci.ub 
## tau^2    2.8121  1.1862  46.0152 
## tau      1.6769  1.0891   6.7835 
## I^2(%)  90.7133 80.4702  99.3783 
## H^2     10.7681  5.1204 160.8371

Slide 56-62

rma.regn <- rma(yi=betahat, sei=stderr, mods= ~I(n<4000), data=genet, method="FE")
rma.regn

## 
## Fixed-Effects with Moderators Model (k = 6)
## 
## I^2 (residual heterogeneity / unaccounted variability): 60.01%
## H^2 (unaccounted variability / sampling variability):   2.50
## R^2 (amount of heterogeneity accounted for):            0.00%
## 
## Test for Residual Heterogeneity:
## QE(df = 4) = 10.0018, p-val = 0.0404
## 
## Test of Moderators (coefficient 2):
## QM(df = 1) = 0.0259, p-val = 0.8721
## 
## Model Results:
## 
##                  estimate      se     zval    pval    ci.lb    ci.ub     ​ 
## intrcpt           -0.4454  0.1031  -4.3212  <.0001  -0.6474  -0.2434  *** 
## I(n < 4000)TRUE   -0.0335  0.2082  -0.1610  0.8721  -0.4416   0.3746      
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

with(genet, sum(betahat^2/stderr^2))

## [1] 35.68249

rma.regn$QM + rma.regn$QE

## [1] 10.02773

rma2$QE

## [1] 10.02773

# an example with a multiple categories, i.e. ANOVA
genet$cutn <- with(genet, cut(n,c(0,3250,8000,Inf), labels=c("S","M","L")))
rma(yi=betahat, sei=stderr, mods= ~cutn, data=genet, method="FE")

## 
## Fixed-Effects with Moderators Model (k = 6)
## 
## I^2 (residual heterogeneity / unaccounted variability): 62.02%
## H^2 (unaccounted variability / sampling variability):   2.63
## R^2 (amount of heterogeneity accounted for):            0.00%
## 
## Test for Residual Heterogeneity:
## QE(df = 3) = 7.8985, p-val = 0.0482
## 
## Test of Moderators (coefficients 2:3):
## QM(df = 2) = 2.1293, p-val = 0.3449
## 
## Model Results:
## 
##          estimate      se     zval    pval    ci.lb    ci.ub   ​ 
## intrcpt   -0.5515  0.2163  -2.5495  0.0108  -0.9755  -0.1275  * 
## cutnM      0.3559  0.2941   1.2100  0.2263  -0.2206   0.9324    
## cutnL      0.0415  0.2441   0.1701  0.8649  -0.4369   0.5200    
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# slide 47 and 49, full meta-regression

reduce.fact <- 10 # how much to shrink the boxes
rma.regn2 <- rma(yi=betahat, sei=stderr, mods= ~av.age, data=genet, method="FE")

plot(betahat~av.age, data=genet, ylim=c(-1.75,0.5), main="Meta-regression of genetic effect on study's average age", col=0)
with(genet, symbols(y=betahat, x=av.age, squares=1/stderr/reduce.fact , inches=FALSE, add=TRUE, bg=1))
with(genet, segments(av.age, betahat-1.96*stderr, av.age, betahat+1.96*stderr))
abline(rma.regn2$beta)
mtext(side=3, line=0.5, col=4, "Blue: plain fixed-effects estimate & CI, at x=overall average age")

xF <- with(genet, sum(av.age/stderr^2)/sum(1/stderr^2))
betaFhat <- rma2$beta[[1,1]]
sebetaFhat <- rma2$se
segments(xF, betaFhat-1.96*sebetaFhat, xF, betaFhat+1.96*sebetaFhat, col=4)
symbols(y=betaFhat, x=xF, squares=1/sebetaFhat/reduce.fact , inches=FALSE, add=TRUE, bg=4, col=4)

Slide 67-70

With \(k\)=6, the Zinc and genetic association candidates are not good candidates for assessment of publication bias. So we’ll use the package’s example instead

res <- rma(yi, vi, data=dat.hackshaw1998, measure="OR", method="FE")

funnel(res)

### carry out trim-and-fill analysis
taf <- trimfill(res)
taf 

## 
## Estimated number of missing studies on the left side: 7 (SE = 4.0405)
## 
## Fixed-Effects Model (k = 44)
## 
## I^2 (total heterogeneity / total variability):   30.44%
## H^2 (total variability / sampling variability):  1.44
## 
## Test for Heterogeneity:
## Q(df = 43) = 61.8207, p-val = 0.0313
## 
## Model Results:
## 
## estimate      se    zval    pval   ci.lb   ci.ub     ​ 
##   0.1556  0.0364  4.2706  <.0001  0.0842  0.2270  *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

### draw funnel plot with missing studies filled in
funnel(taf, legend=TRUE)

### Comparing Egger's test to our earlier meta-regression
rma2.full <- rma(yi=betahat, sei=stderr, ni=n, data=genet, method="FE")
#Egger's test
regtest(x=rma2.full, predictor="ni")

## 
## Regression Test for Funnel Plot Asymmetry
## 
## Model:     fixed-effects meta-regression model
## Predictor: sample size
## 
## Test for Funnel Plot Asymmetry: z = -0.3594, p = 0.7193

rma(yi=betahat, sei=stderr, mod=~n, data=genet, method="FE")

## 
## Fixed-Effects with Moderators Model (k = 6)
## 
## I^2 (residual heterogeneity / unaccounted variability): 59.59%
## H^2 (unaccounted variability / sampling variability):   2.47
## R^2 (amount of heterogeneity accounted for):            0.00%
## 
## Test for Residual Heterogeneity:
## QE(df = 4) = 9.8986, p-val = 0.0422
## 
## Test of Moderators (coefficient 2):
## QM(df = 1) = 0.1291, p-val = 0.7193
## 
## Model Results:
## 
##          estimate      se     zval    pval    ci.lb   ci.ub   ​ 
## intrcpt   -0.3632  0.2671  -1.3599  0.1739  -0.8867  0.1603    
## n         -0.0000  0.0000  -0.3594  0.7193  -0.0001  0.0001    
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Slide 86

Includes some code not shown on slides just computing posterior mean and variance of various parameters, assuming Normal likelihoods for the \(\hat\beta_i\) and Normal priors for the \(\beta_i\).

PostBeta.H <- function(tau2,psi2,nu=0,betahat,sigma2)
 { 
   k <- length(sigma2) 
   Pos.Cov <- solve(diag(sigma2^(-1)) + solve(tau2*diag(k) + psi2*matrix(rep(1,k^2),ncol=k)))
   Pos.E  <-   Pos.Cov%*%( diag(sigma2^(-1))%*%betahat +  solve(tau2*diag(k) + psi2*matrix(rep(1,k^2),ncol=k)) %*% rep(1,k)*nu ) 
   return(list(E.Beta=Pos.E,Cov.Beta=Pos.Cov))
  }

## Function to obtain posterior mean and variance of BetaF
PostBetaF <- function(tau2,psi2,nu=0,betahat,sigma2)
 { 
   k <- length(betahat)  
   post <- PostBeta.H(tau2,psi2,nu,betahat,sigma2)
   Var.BetaF <-  t(sigma2^(-1))%*% post$Cov.Beta  %*%(sigma2^(-1))/(sum(sigma2^(-1)))^2
   E.BetaF  <-   t(sigma2^(-1))%*% post$E.Beta / sum(sigma2^(-1))
   return(c(E.BetaF=E.BetaF, Var.BetaF=Var.BetaF))
  }

# non-Bayesian analysis
round( c(zinc.meta["estimate"],zinc.meta["se"]^2), 3)

## estimate       se 
##    2.042    0.043

# homogeneous, tight prior on location
round( PostBetaF(tau2=0.1, psi2=0.1,nu=0,zinc$betahat,zinc$se^2) , 3)

##   E.BetaF Var.BetaF 
##     1.566     0.032

# heterogeneous, tight prior on location
round( PostBetaF(tau2=100, psi2=0.1,nu=0,zinc$betahat,zinc$se^2) , 3)

##   E.BetaF Var.BetaF 
##     2.041     0.043

# homogeneous, vague prior on location
round( PostBetaF(tau2=0.1, psi2=100,nu=0,zinc$betahat,zinc$se^2) , 3)

##   E.BetaF Var.BetaF 
##     2.042     0.043

# heterogeneous, vague prior on location
round( PostBetaF(tau2=100, psi2=100,nu=0,zinc$betahat,zinc$se^2) , 3)

##   E.BetaF Var.BetaF 
##     2.042     0.043

Slide 96-97

Bayes using JAGS: note you will likely need to first install JAGS (https://sourceforge.net/projects/mcmc-jags/)

Sys.setenv(JAGS_HOME="C:\\Program Files\\JAGS\\JAGS-4.3.1")
library("rjags") # use install.packages("rjags") first if needed

## Loading required package: coda

## Linked to JAGS 4.3.1

## Loaded modules: basemod,bugs

# set up an external model file, using cat()
cat(file="metaprog.txt", "model{ 
for(i in 1:6){
    prec[i] <- 1/(se[i]*se[i])             # precision not variance
    betahat[i] ~ dnorm( beta[i], prec[i] ) # 
    beta[i] ~ dnorm(mu, 10)                # tau^2=0.1, homogeneous
    w[i] <- prec[i]/sum(prec[1:6])
}
mu ~ dnorm(0, 10)   # psi^2=10, tight prior on mean
betaF <- inprod(w[1:6], beta[1:6])           # inner product, sum(w[i]*beta[i])
}")

# have a look at it
file.show("metaprog.txt")

# do the calculations
jags1 <- jags.model("metaprog.txt", data=zinc[,c("betahat","se")])

## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 6
##    Unobserved stochastic nodes: 7
##    Total graph size: 45
## 
## Initializing model

update(jags1, 100000)
# summarize the posteriors for two parameters:
summary( coda.samples(jags1, c("mu","betaF"), n.iter=100000) )

## 
## Iterations = 100001:2e+05
## Thinning interval = 1 
## Number of chains = 1 
## Sample size per chain = 1e+05 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##        Mean     SD  Naive SE Time-series SE
## betaF 1.567 0.1781 0.0005631       0.001118
## mu    1.118 0.1994 0.0006307       0.001361
## 
## 2. Quantiles for each variable:
## 
##         2.5%    25%   50%   75% 97.5%
## betaF 1.2180 1.4468 1.569 1.687 1.915
## mu    0.7258 0.9832 1.118 1.253 1.507

plot( coda.samples(jags1, c("mu","betaF"), n.iter=10000) )

Exercise 4

Implement Bayesian analysis for one of the earlier fixed-effects or random-effects meta-analyses, using priors of your own choice.