Please do not touch anything in this section, otherwise this notebook might not work properly. You have been warned! Also, if you have no clue what you are staring at, please consult our Preface chapter.
source("run_me_first.R")
In this exercise, we want you to estimate the fixed-effect and the random-effects model for the dat.bcg
data (using risk ratios).
Next, we would like you to manually estimate the “common” effect, i.e. the summary effect according to the fixed-effect model, its variance/standard error and the 95% CI. You will need the formulas and the slide “Summary” in the subsection ‘Distinction between fixed- and random-effects’ might be useful to get you started. Tip: The function sum(x)
calculates the sum of a vector x
.
Note: you will:
wi <- 1/dat.bcg$vi
)ma.fem <- sum(wi * dat.bcg$yi) / sum(wi)
)var.ma.fem <- 1/sum(wi)
)se.ma.fem <- sqrt(var.ma.fem)
)ma.fem.ll <- ma.fem - 1.96 * se.ma.fem
)## Please insert your solution here. Of course, feel free to add new code cells.
## Please insert your solution here. Of course, feel free to add new code cells.
## Please insert your solution here. Of course, feel free to add new code cells.
## Please insert your solution here. Of course, feel free to add new code cells.
## Please insert your solution here. Of course, feel free to add new code cells.
Of course, as you may already know, you do not have to estimate summary effects by hand. metafor
's rma()
function can do that for you.
For this exercise, you can apply the rma()
function to replicate the results from exercise 1, i.e. estimate a fixed-effect summary. If you forgot how to apply the ?rma
function, check the metafor introduction slides and exercises. (Remember, you only need for now to know three arguments in the ?rma function: the effect-size argument yi
, the variance argument vi
, and the argument method
. With the latter, you specify which effect-size model you want to apply.) Here, we need the fixed-effect model and therefore the argument is method = "FE"
.
Note: See the slides from the morning, the command that you need were used in the example, but be sure to
use the method = "FE"
for this exercise.
## Please insert your solution here. Of course, feel free to add new code cells.
Finally, we will estimate the random-effects model. This morning, we have discussed that there exist several
between-study variance estimates and we introduced the method-of-moments estimator (also known as DerSimonian and Laird estimator). In metafor
’s rma()
function, this can be accomplished with the method = "DL"
argument.
So, estimate the random-effects model and discuss the following questions:
## Please insert your solution here. Of course, feel free to add new code cells.
## Please insert your solution here. Of course, feel free to add new code cells.
## Please insert your solution here. Of course, feel free to add new code cells.
## Please insert your solution here. Of course, feel free to add new code cells.
## Solution according to http://www.metafor-project.org/doku.php/analyses:berkey1995.
## IMPORTANT: Different metehod ("EB") applied!
library(metafor)
dat.bcg <- escalc(measure = "RR", ai = tpos, bi = tneg, ci = cpos, di = cneg, data = dat.bcg)
dat.bcg$vi <- with(dat.bcg, sum(tneg/tpos)/(13*(tneg+tpos)) +
sum(cneg/cpos)/(13*(cneg+cpos)))
ma.rem <- rma(yi = yi, vi = vi, method = "EB", data = dat.bcg)
ma.rem
Random-Effects Model (k = 13; tau^2 estimator: EB) tau^2 (estimated amount of total heterogeneity): 0.2682 (SE = 0.1801) tau (square root of estimated tau^2 value): 0.5178 I^2 (total heterogeneity / total variability): 87.49% H^2 (total variability / sampling variability): 7.99 Test for Heterogeneity: Q(df = 12) = 85.8625, p-val < .0001 Model Results: estimate se zval pval ci.lb ci.ub -0.5429 0.1842 -2.9474 0.0032 -0.9040 -0.1819 ** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
For those of you who are already bored, we suggest that you try estimating the RE model by hand.
Note: to estimate between-studies variability you will need to compute $Q_T$ and use this information to estimate the $\tau^2_{DL}$. Formulas are provided in our slides.
## Please insert your solution here. Of course, feel free to add new code cells.
## Please insert your solution here. Of course, feel free to add new code cells.
## Please insert your solution here. Of course, feel free to add new code cells.
## Please insert your solution here. Of course, feel free to add new code cells.
## Please insert your solution here. Of course, feel free to add new code cells.