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 conduct mixed effect meta-regression to attempt to explain effect size heterogeneity.
Load the data into R and then run a simple (unconditional) mixed effect meta-regression.
## Please insert your solution here. Of course, feel free to add new code cells.
In a next step we are wondering if the moderator variable ablat
can explain heterogeneity in the effect size
distribution. ablat
denotes the absolute latitude (Wikipedia article). The effect of latitude on the vaccine can be assumed as follows:
"Latitude accounts for variation in rainfall, humidity, environmental mycobacteria that may produce natural immunity, temperature, and other factors that may have biological implications for vaccine efficacy. Preparation of the live vaccine requires refrigeration; unrefrigerated, it would spoil more quickly in warmer temperatures. Furthermore, direct
exposure of the vaccine to sunlight may reduce counts of live bacteria" (Berkey et al. 1995: 396).
In order to investigate the effect of ablat
, run a mixed meta-regression model.
Check our metafor
introduction slides for help with the code for mixed effects models. Come up with a meaningful
interpretation.
Tip: rma(y, v, mods= ~ x_1 + ..., data = dataset)
## 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.
Next we were wondering what happens if transform ablat into a grouped variable named
ablat_gr. Note that we need to define a predictor as a factor
variable. There are different approaches of creating factor variables in R:
recode
command from the car package. After doing that, you have to convert the newly created variable to a factor variable by using the factor() command:table(dat.bcg$ablat)
13 18 19 27 33 42 44 52 55 2 1 1 1 2 2 2 1 1
library(car)
dat.bcg$ablat_gr <- recode(dat.bcg$ablat, "0:33 = '[0,34)'; 24:90 = '[34,90)'")
dat.bcg$ablat_gr <- factor(dat.bcg$ablat_gr)
table(dat.bcg$ablat_gr)
[0,34) [34,90) 7 6
factor()
command, e.g. factor(dat.bcg$ablat)
(note: your new factor variable will have 26 levels).This variable indicates if the study took place between 34 degrees latitude and the north/south pole (remember, it is the absolute latitude).
## Please insert your solution here. Of course, feel free to add new code cells.
If you want to change the reference level, you can use the relevel()
command. So, for instance, if you want to use level [0,34)
as new reference category:
mdl_grouped <- rma(yi ~ relevel(ablat_gr, ref = "[34,90)"), vi, data = dat.bcg, digits = 3)
summary(mdl_grouped)
Mixed-Effects Model (k = 13; tau^2 estimator: REML) logLik deviance AIC BIC AICc -7.073 14.146 20.146 21.339 23.574 tau^2 (estimated amount of residual heterogeneity): 0.084 (SE = 0.063) tau (square root of estimated tau^2 value): 0.289 I^2 (residual heterogeneity / unaccounted variability): 70.82% H^2 (unaccounted variability / sampling variability): 3.43 R^2 (amount of heterogeneity accounted for): 73.31% Test for Residual Heterogeneity: QE(df = 11) = 41.789, p-val < .001 Test of Moderators (coefficient(s) 2): QM(df = 1) = 17.145, p-val < .001 Model Results: estimate se zval pval intrcpt -1.194 0.169 -7.075 <.001 relevel(ablat_gr, ref = "[34,90)")[0,34) 0.920 0.222 4.141 <.001 ci.lb ci.ub intrcpt -1.524 -0.863 *** relevel(ablat_gr, ref = "[34,90)")[0,34) 0.485 1.356 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
If you don’t like that your output looks kind of ugly because the relevel() command is included, you could generate a new variable:
dat.bcg$ablat_gr2 <- relevel(dat.bcg$ablat_gr, ref = "[34,90)")
head(dat.bcg)
trial | author | year | tpos | tneg | cpos | cneg | ablat | alloc | yi | vi | ablat_gr | ablat_gr2 |
---|---|---|---|---|---|---|---|---|---|---|---|---|
1 | Aronson | 1948 | 4 | 119 | 11 | 128 | 44 | random | -0.8893113 | 0.325584765 | [34,90) | [34,90) |
2 | Ferguson & Simes | 1949 | 6 | 300 | 29 | 274 | 55 | random | -1.5853887 | 0.194581121 | [34,90) | [34,90) |
3 | Rosenthal et al | 1960 | 3 | 228 | 11 | 209 | 42 | random | -1.3480731 | 0.415367965 | [34,90) | [34,90) |
4 | Hart & Sutherland | 1977 | 62 | 13536 | 248 | 12619 | 52 | random | -1.4415512 | 0.020010032 | [34,90) | [34,90) |
5 | Frimodt-Moller et al | 1973 | 33 | 5036 | 47 | 5761 | 13 | alternate | -0.2175473 | 0.051210172 | [0,34) | [0,34) |
6 | Stein & Aronson | 1953 | 180 | 1361 | 372 | 1079 | 44 | alternate | -0.7861156 | 0.006905618 | [34,90) | [34,90) |
mdl_grouped <- rma(yi ~ ablat_gr2, vi, data = dat.bcg, digits = 3)
summary(mdl_grouped)
Mixed-Effects Model (k = 13; tau^2 estimator: REML) logLik deviance AIC BIC AICc -7.073 14.146 20.146 21.339 23.574 tau^2 (estimated amount of residual heterogeneity): 0.084 (SE = 0.063) tau (square root of estimated tau^2 value): 0.289 I^2 (residual heterogeneity / unaccounted variability): 70.82% H^2 (unaccounted variability / sampling variability): 3.43 R^2 (amount of heterogeneity accounted for): 73.31% Test for Residual Heterogeneity: QE(df = 11) = 41.789, p-val < .001 Test of Moderators (coefficient(s) 2): QM(df = 1) = 17.145, p-val < .001 Model Results: estimate se zval pval ci.lb ci.ub intrcpt -1.194 0.169 -7.075 <.001 -1.524 -0.863 *** ablat_gr2[0,34) 0.920 0.222 4.141 <.001 0.485 1.356 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Use rma(y, v, mods = ~x_1 + ... - 1, data = dataset)
(-1
removes the intercept) to obtain the weighted average effect sizes and standard errors for each group of ablat_gr .
IMPORTANT: Do not use the Q values!
Which group has the largest effect? And, more importantly, what does that mean?
## Please insert your solution here. Of course, feel free to add new code cells.