library(metaSEM)
## Warning: package 'metaSEM' was built under R version 3.4.1
## Loading required package: OpenMx
## Warning: package 'OpenMx' was built under R version 3.4.1
## OpenMx is not compiled to take advantage of computers with multiple cores.
## "SLSQP" is set as the default optimizer in OpenMx.
## mxOption(NULL, "Gradient algorithm") is set at "central".
## mxOption(NULL, "Optimality tolerance") is set at "6.3e-14".
## mxOption(NULL, "Gradient iterations") is set at "2".
data("BCG")
head(BCG)
## Trial Author Year VD VWD NVD NVWD Latitude Allocation
## 1 1 Aronson 1948 4 119 11 128 44 random
## 2 2 Ferguson & Simes 1949 6 300 29 274 55 random
## 3 3 Rosenthal et al 1960 3 228 11 209 42 random
## 4 4 Hart & Sutherland 1977 62 13536 248 12619 52 random
## 5 5 Frimodt-Moller et al 1973 33 5036 47 5761 13 alternate
## 6 6 Stein & Aronson 1953 180 1361 372 1079 44 alternate
## ln_OR v_ln_OR ln_Odd_V ln_Odd_NV v_ln_Odd_V cov_V_NV
## 1 -0.9386941 0.357124952 -3.392829 -2.454135 0.258403361 0
## 2 -1.6661907 0.208132394 -3.912023 -2.245832 0.170000000 0
## 3 -1.3862944 0.433413078 -4.330733 -2.944439 0.337719298 0
## 4 -1.4564435 0.020314413 -5.385974 -3.929530 0.016202909 0
## 5 -0.2191411 0.051951777 -5.027860 -4.808719 0.030501601 0
## 6 -0.9581220 0.009905266 -2.023018 -1.064896 0.006290309 0
## v_ln_Odd_NV
## 1 0.098721591
## 2 0.038132394
## 3 0.095693780
## 4 0.004111504
## 5 0.021450177
## 6 0.003614956
library(metafor)
## Loading required package: Matrix
##
## Attaching package: 'Matrix'
## The following objects are masked from 'package:OpenMx':
##
## %&%, expm
## Loading 'metafor' package (version 2.0-0). For an overview
## and introduction to the package please type: help(metafor).
# Calculate the ES using "escalc"
dat = escalc(measure="RR",ai=VD,bi=VWD,ci=NVD,
di = NVWD, data = BCG, append = TRUE)
# print
print(dat)
## Trial Author Year VD VWD NVD NVWD Latitude Allocation
## 1 1 Aronson 1948 4 119 11 128 44 random
## 2 2 Ferguson & Simes 1949 6 300 29 274 55 random
## 3 3 Rosenthal et al 1960 3 228 11 209 42 random
## 4 4 Hart & Sutherland 1977 62 13536 248 12619 52 random
## 5 5 Frimodt-Moller et al 1973 33 5036 47 5761 13 alternate
## 6 6 Stein & Aronson 1953 180 1361 372 1079 44 alternate
## 7 7 Vandiviere et al 1973 8 2537 10 619 19 random
## 8 8 TPT Madras 1980 505 87886 499 87892 13 random
## 9 9 Coetzee & Berjak 1968 29 7470 45 7232 27 random
## 10 10 Rosenthal et al 1961 17 1699 65 1600 42 systematic
## 11 11 Comstock et al 1974 186 50448 141 27197 18 systematic
## 12 12 Comstock & Webster 1969 5 2493 3 2338 33 systematic
## 13 13 Comstock et al 1976 27 16886 29 17825 33 systematic
## ln_OR v_ln_OR ln_Odd_V ln_Odd_NV v_ln_Odd_V cov_V_NV
## 1 -0.93869414 0.357124952 -3.392829 -2.454135 0.258403361 0
## 2 -1.66619073 0.208132394 -3.912023 -2.245832 0.170000000 0
## 3 -1.38629436 0.433413078 -4.330733 -2.944439 0.337719298 0
## 4 -1.45644355 0.020314413 -5.385974 -3.929530 0.016202909 0
## 5 -0.21914109 0.051951777 -5.027860 -4.808719 0.030501601 0
## 6 -0.95812204 0.009905266 -2.023018 -1.064896 0.006290309 0
## 7 -1.63377584 0.227009675 -5.759296 -4.125520 0.125394166 0
## 8 0.01202060 0.004006962 -5.159237 -5.171258 0.001991576 0
## 9 -0.47174604 0.056977124 -5.551354 -5.079608 0.034616627 0
## 10 -1.40121014 0.075421726 -4.604582 -3.203372 0.059412111 0
## 11 -0.34084965 0.012525134 -5.602952 -5.262102 0.005396166 0
## 12 0.44663468 0.534162172 -6.211804 -6.658439 0.200401123 0
## 13 -0.01734187 0.071635117 -6.438403 -6.421061 0.037096258 0
## v_ln_Odd_NV yi vi
## 1 0.098721591 -0.8893 0.3256
## 2 0.038132394 -1.5854 0.1946
## 3 0.095693780 -1.3481 0.4154
## 4 0.004111504 -1.4416 0.0200
## 5 0.021450177 -0.2175 0.0512
## 6 0.003614956 -0.7861 0.0069
## 7 0.101615509 -1.6209 0.2230
## 8 0.002015386 0.0120 0.0040
## 9 0.022360497 -0.4694 0.0564
## 10 0.016009615 -1.3713 0.0730
## 11 0.007128967 -0.3394 0.0124
## 12 0.333761049 0.4459 0.5325
## 13 0.034538860 -0.0173 0.0714
# Call `rma' to fit the BCG data
meta.RE = rma(yi, vi, data = dat)
# Print the summary
meta.RE
##
## Random-Effects Model (k = 13; tau^2 estimator: REML)
##
## tau^2 (estimated amount of total heterogeneity): 0.3132 (SE = 0.1664)
## tau (square root of estimated tau^2 value): 0.5597
## I^2 (total heterogeneity / total variability): 92.22%
## H^2 (total variability / sampling variability): 12.86
##
## Test for Heterogeneity:
## Q(df = 12) = 152.2330, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## -0.7145 0.1798 -3.9744 <.0001 -1.0669 -0.3622 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
forest((meta.RE))
# Call `rma' to fit the BCG data
meta.FE = rma(yi, vi, method = "FE", data = dat)
# Print the summary
meta.FE
##
## Fixed-Effects Model (k = 13)
##
## Test for Heterogeneity:
## Q(df = 12) = 152.2330, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## -0.4303 0.0405 -10.6247 <.0001 -0.5097 -0.3509 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
forest((meta.FE))
# Call `rma' to fit the BCG data
meta.DL = rma(yi, vi, method = "DL", data = dat)
# Print the summary
meta.DL
##
## Random-Effects Model (k = 13; tau^2 estimator: DL)
##
## tau^2 (estimated amount of total heterogeneity): 0.3088 (SE = 0.2299)
## tau (square root of estimated tau^2 value): 0.5557
## I^2 (total heterogeneity / total variability): 92.12%
## H^2 (total variability / sampling variability): 12.69
##
## Test for Heterogeneity:
## Q(df = 12) = 152.2330, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## -0.7141 0.1787 -3.9952 <.0001 -1.0644 -0.3638 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
forest((meta.DL))
# Call `rma' to fit the BCG data
meta.HS = rma(yi, vi, method = "HS", data = dat)
# Print the summary
meta.HS
##
## Random-Effects Model (k = 13; tau^2 estimator: HS)
##
## tau^2 (estimated amount of total heterogeneity): 0.2284 (SE = 0.1279)
## tau (square root of estimated tau^2 value): 0.4779
## I^2 (total heterogeneity / total variability): 89.63%
## H^2 (total variability / sampling variability): 9.64
##
## Test for Heterogeneity:
## Q(df = 12) = 152.2330, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## -0.7045 0.1587 -4.4408 <.0001 -1.0155 -0.3936 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
forest((meta.HS))
# Call `rma' to fit the BCG data
meta.HE = rma(yi, vi, method = "HE", data = dat)
# Print the summary
meta.HE
##
## Random-Effects Model (k = 13; tau^2 estimator: HE)
##
## tau^2 (estimated amount of total heterogeneity): 0.3286 (SE = 0.2071)
## tau (square root of estimated tau^2 value): 0.5732
## I^2 (total heterogeneity / total variability): 92.56%
## H^2 (total variability / sampling variability): 13.44
##
## Test for Heterogeneity:
## Q(df = 12) = 152.2330, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## -0.7159 0.1833 -3.9059 <.0001 -1.0751 -0.3567 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
forest((meta.HE))
# Call `rma' to fit the BCG data
meta.SJ = rma(yi, vi, method = "SJ", data = dat)
# Print the summary
meta.SJ
##
## Random-Effects Model (k = 13; tau^2 estimator: SJ)
##
## tau^2 (estimated amount of total heterogeneity): 0.3455 (SE = 0.1497)
## tau (square root of estimated tau^2 value): 0.5878
## I^2 (total heterogeneity / total variability): 92.90%
## H^2 (total variability / sampling variability): 14.08
##
## Test for Heterogeneity:
## Q(df = 12) = 152.2330, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## -0.7172 0.1871 -3.8343 0.0001 -1.0839 -0.3506 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
forest((meta.SJ))
# Call `rma' to fit the BCG data
meta.EB = rma(yi, vi, method = "EB", data = dat)
# Print the summary
meta.EB
##
## Random-Effects Model (k = 13; tau^2 estimator: EB)
##
## tau^2 (estimated amount of total heterogeneity): 0.3181 (SE = 0.1737)
## tau (square root of estimated tau^2 value): 0.5640
## I^2 (total heterogeneity / total variability): 92.33%
## H^2 (total variability / sampling variability): 13.04
##
## Test for Heterogeneity:
## Q(df = 12) = 152.2330, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## -0.7150 0.1809 -3.9525 <.0001 -1.0695 -0.3604 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
forest((meta.EB))
# Call `rma' to fit the BCG data
meta.ML = rma(yi, vi, method = "ML", data = dat)
# Print the summary
meta.ML
##
## Random-Effects Model (k = 13; tau^2 estimator: ML)
##
## tau^2 (estimated amount of total heterogeneity): 0.2800 (SE = 0.1443)
## tau (square root of estimated tau^2 value): 0.5292
## I^2 (total heterogeneity / total variability): 91.38%
## H^2 (total variability / sampling variability): 11.60
##
## Test for Heterogeneity:
## Q(df = 12) = 152.2330, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## -0.7112 0.1719 -4.1374 <.0001 -1.0481 -0.3743 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
forest((meta.ML))
metaReg = rma(yi, vi, mods = ~Latitude+Year+Allocation, data = dat)
# Print the meta-regression results
metaReg
##
## Mixed-Effects Model (k = 13; tau^2 estimator: REML)
##
## tau^2 (estimated amount of residual heterogeneity): 0.1796 (SE = 0.1425)
## tau (square root of estimated tau^2 value): 0.4238
## I^2 (residual heterogeneity / unaccounted variability): 73.09%
## H^2 (unaccounted variability / sampling variability): 3.72
## R^2 (amount of heterogeneity accounted for): 42.67%
##
## Test for Residual Heterogeneity:
## QE(df = 8) = 26.2030, p-val = 0.0010
##
## Test of Moderators (coefficient(s) 2:5):
## QM(df = 4) = 9.5254, p-val = 0.0492
##
## Model Results:
##
## estimate se zval pval ci.lb
## intrcpt -14.4984 38.3943 -0.3776 0.7057 -89.7498
## Latitude -0.0236 0.0132 -1.7816 0.0748 -0.0495
## Year 0.0075 0.0194 0.3849 0.7003 -0.0306
## Allocationrandom -0.3421 0.4180 -0.8183 0.4132 -1.1613
## Allocationsystematic 0.0101 0.4467 0.0226 0.9820 -0.8654
## ci.ub
## intrcpt 60.7531
## Latitude 0.0024 .
## Year 0.0456
## Allocationrandom 0.4772
## Allocationsystematic 0.8856
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(metaReg)
metaReg.lat = rma(yi, vi, mods = ~Latitude, data = dat)
# Print the meta-regression results
metaReg.lat
##
## Mixed-Effects Model (k = 13; tau^2 estimator: REML)
##
## tau^2 (estimated amount of residual heterogeneity): 0.0764 (SE = 0.0591)
## tau (square root of estimated tau^2 value): 0.2763
## I^2 (residual heterogeneity / unaccounted variability): 68.39%
## H^2 (unaccounted variability / sampling variability): 3.16
## R^2 (amount of heterogeneity accounted for): 75.62%
##
## Test for Residual Heterogeneity:
## QE(df = 11) = 30.7331, p-val = 0.0012
##
## Test of Moderators (coefficient(s) 2):
## QM(df = 1) = 16.3571, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt 0.2515 0.2491 1.0095 0.3127 -0.2368 0.7397
## Latitude -0.0291 0.0072 -4.0444 <.0001 -0.0432 -0.0150 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(metaReg.lat)
# Create a new latitude vector
newlat = 0:60
# Using the meta-regression and calculate the predicted values
preds = predict(metaReg.lat, newmods = newlat, transf = exp)
# Use the inverse-variance to create a weighting for the data
wi = 1/sqrt(dat$vi)
size = 1 + 3 * (wi - min(wi))/(max(wi) - min(wi))
# Plot the RR
plot(dat$Latitude, exp(dat$yi), pch = 19, cex = size,
xlab = "Latitude", ylab = "Relative Risk",
las = 1, bty = "l", log = "y")
# Add a thicker line for the meta-regression and the CIs
lines(newlat, preds$pred, lwd=3)
lines(newlat, preds$ci.lb, lty = "dashed")
lines(newlat, preds$ci.ub, lty = "dashed")
# Add a dotted horizontal line for equal-effectiveness
abline(h = 1, lwd=3,lty = "dotted")
