library(meta)
## Loading 'meta' package (version 4.8-4).
## Type 'help(meta)' for a brief overview.
library(rmeta)
## Loading required package: grid
library(metafor)
## Loading required package: Matrix
## Loading 'metafor' package (version 2.0-0). For an overview
## and introduction to the package please type: help(metafor).
##
## Attaching package: 'metafor'
## The following objects are masked from 'package:meta':
##
## baujat, forest, funnel, funnel.default, labbe, radial,
## trimfill
data11 <- read.csv("data11.csv")
str(data11)
## 'data.frame': 37 obs. of 6 variables:
## $ study: int 1 2 3 4 5 6 7 8 9 10 ...
## $ Ee : int 1 2 3 4 5 6 7 8 9 10 ...
## $ Ne : int 28 73 70 37 12 28 22 32 52 25 ...
## $ Ec : int 11 51 42 15 1 14 8 44 38 17 ...
## $ Nc : int 25 132 111 42 16 48 36 51 88 45 ...
## $ X : int 32 120 110 48 12 55 28 56 81 36 ...
summary(data11$Ee/data11$Ne)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0274 0.2857 0.6757 0.5802 0.8148 0.9487
summary(data11$Ec/data11$Nc)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.2222 0.3571 0.3529 0.4400 0.8627
ms1 <- metabin(Ee, Ne, Ec, Nc, data=data11, sm="OR")
forest(ms1, label.left="Rx worse", label.right="Rx better", ff.lr="bold")
funnel(ms1)
args(funnel.meta)
## function (x, xlim = NULL, ylim = NULL, xlab = NULL, ylab = NULL,
## comb.fixed = x$comb.fixed, comb.random = x$comb.random, axes = TRUE,
## pch = if (!inherits(x, "trimfill")) 21 else ifelse(x$trimfill,
## 1, 21), text = NULL, cex = 1, lty.fixed = 2, lty.random = 9,
## lwd = 1, lwd.fixed = lwd, lwd.random = lwd, col = "black",
## bg = "darkgray", col.fixed = "black", col.random = "black",
## log = "", yaxis = "se", contour.levels = NULL, col.contour,
## ref = ifelse(backtransf & is.relative.effect(x$sm), 1, 0),
## level = x$level, studlab = FALSE, cex.studlab = 0.8, pos.studlab = 2,
## backtransf = x$backtransf, ...)
## NULL
funnel(ms1, comb.random=FALSE, pch=16,
contour=c(0.9, 0.95, 0.99),
col.contour=c("darkgray", "gray","lightgray"))
legend(0.25, 1.25,
c("0.1 > p > 0.05", "0.05 > p > 0.01", "< 0.01"),
fill=c("darkgray", "gray","lightgray"), bty="n")
plot(1/ms1$seTE, ms1$TE/ms1$seTE)
radial(ms1)
metabias(ms1, method="rank")
##
## Rank correlation test of funnel plot asymmetry
##
## data: ms1
## z = 2.8512, p-value = 0.004356
## alternative hypothesis: asymmetry in funnel plot
## sample estimates:
## ks se.ks
## 218.00000 76.45914
metabias(ms1, method="linreg")
##
## Linear regression test of funnel plot asymmetry
##
## data: ms1
## t = 2.9663, df = 35, p-value = 0.005401
## alternative hypothesis: asymmetry in funnel plot
## sample estimates:
## bias se.bias slope
## 4.027707 1.357805 -1.659405
reg <- lm(I(ms1$TE/ms1$seTE) ~ I(1/ms1$seTE))
summary(reg)
##
## Call:
## lm(formula = I(ms1$TE/ms1$seTE) ~ I(1/ms1$seTE))
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.2440 -0.6324 0.9702 1.7010 4.4424
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.0277 1.3578 2.966 0.0054 **
## I(1/ms1$seTE) -1.6594 0.7608 -2.181 0.0360 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.866 on 35 degrees of freedom
## Multiple R-squared: 0.1197, Adjusted R-squared: 0.09451
## F-statistic: 4.757 on 1 and 35 DF, p-value: 0.03598
radial(ms1)
abline(reg)
metabias(ms1, method = "linreg", plotit=TRUE)
##
## Linear regression test of funnel plot asymmetry
##
## data: ms1
## t = 2.9663, df = 35, p-value = 0.005401
## alternative hypothesis: asymmetry in funnel plot
## sample estimates:
## bias se.bias slope
## 4.027707 1.357805 -1.659405
metabias(ms1, method = "mm")
##
## Linear regression test of funnel plot asymmetry (methods of
## moment)
##
## data: ms1
## t = 4.0819, df = 35, p-value = 0.0002462
## alternative hypothesis: asymmetry in funnel plot
## sample estimates:
## bias se.bias slope
## 4.196408 1.028060 -1.765899
metabias(ms1, method = "score")
##
## Linear regression test of funnel plot asymmetry (efficient score)
##
## data: ms1
## t = 4.567, df = 35, p-value = 5.888e-05
## alternative hypothesis: asymmetry in funnel plot
## sample estimates:
## bias se.bias slope
## 8.684086 1.901467 -3.549595
metabias(ms1, method = "peters")
##
## Linear regression test of funnel plot asymmetry (based on sample
## size)
##
## data: ms1
## t = 5.5303, df = 35, p-value = 3.214e-06
## alternative hypothesis: asymmetry in funnel plot
## sample estimates:
## bias se.bias slope
## 258.47072 46.73699 -2.56456
metabias(ms1, method = "count")
##
## Rank correlation test of funnel plot asymmetry (based on counts)
##
## data: ms1
## z = 3.3482, p-value = 0.0008134
## alternative hypothesis: asymmetry in funnel plot
## sample estimates:
## ks se.ks
## 256.00000 76.45914
ms1.asd <- update(ms1, sm="ASD")
summary(ms1.asd)
## Number of studies combined: k = 37
##
## ASD 95%-CI z p-value
## Fixed effect model 0.1027 [0.0666; 0.1388] 5.57 < 0.0001
## Random effects model 0.2535 [0.1015; 0.4056] 3.27 0.0011
##
## Quantifying heterogeneity:
## tau^2 = 0.2072; H = 4.17 [3.76; 4.62]; I^2 = 94.2% [92.9%; 95.3%]
##
## Test of heterogeneity:
## Q d.f. p-value
## 624.89 36 < 0.0001
##
## Details on meta-analytical method:
## - Inverse variance method
## - DerSimonian-Laird estimator for tau^2
funnel(ms1.asd)
metabias(ms1.asd, method = "mm", plotit=TRUE)
##
## Linear regression test of funnel plot asymmetry (methods of
## moment)
##
## data: ms1.asd
## t = 5.6866, df = 35, p-value = 1.999e-06
## alternative hypothesis: asymmetry in funnel plot
## sample estimates:
## bias se.bias slope
## 11.442814 2.012229 -1.137999
tf1 <- trimfill(ms1)
class(tf1)
## [1] "metagen" "meta" "trimfill"
funnel(tf1)
print(tf1, digits=2, comb.fixed=TRUE)
## OR 95%-CI %W(fixed) %W(random)
## 1 0.05 [ 0.01; 0.40] 0.6 1.9
## 2 0.04 [ 0.01; 0.19] 1.4 2.3
## 3 0.07 [ 0.02; 0.25] 2.0 2.3
## 4 0.22 [ 0.06; 0.74] 2.0 2.3
## 5 10.71 [ 1.05; 109.78] 0.6 1.9
## 6 0.66 [ 0.22; 1.98] 2.5 2.4
## 7 1.63 [ 0.50; 5.38] 2.1 2.4
## 8 0.05 [ 0.02; 0.16] 2.3 2.4
## 9 0.28 [ 0.12; 0.63] 4.3 2.5
## 10 1.10 [ 0.40; 2.99] 3.0 2.4
## 11 30.56 [ 5.31; 175.90] 1.0 2.1
## 12 0.51 [ 0.23; 1.12] 4.9 2.5
## 13 3.61 [ 1.23; 10.64] 2.6 2.4
## 14 23.33 [ 4.50; 120.94] 1.1 2.2
## 15 3.28 [ 1.21; 8.88] 3.0 2.4
## 16 6.40 [ 1.89; 21.68] 2.0 2.3
## 17 1.38 [ 0.58; 3.28] 4.0 2.5
## 18 0.71 [ 0.37; 1.38] 6.9 2.5
## 19 0.11 [ 0.05; 0.23] 4.9 2.5
## 20 16.67 [ 4.07; 68.18] 1.5 2.3
## 21 15.60 [ 4.32; 56.31] 1.8 2.3
## 22 94.09 [ 4.77; 1854.23] 0.3 1.6
## 23 11.50 [ 3.42; 38.71] 2.0 2.3
## 24 1.83 [ 0.89; 3.77] 5.7 2.5
## 25 5.06 [ 2.00; 12.80] 3.5 2.4
## 26 2.04 [ 0.67; 6.21] 2.4 2.4
## 27 8.59 [ 1.54; 48.01] 1.0 2.1
## 28 12.83 [ 2.46; 66.92] 1.1 2.2
## 29 12.30 [ 2.95; 51.26] 1.5 2.3
## 30 20.62 [ 6.97; 61.00] 2.5 2.4
## 31 11.19 [ 2.26; 55.33] 1.2 2.2
## 32 429.00 [19.45; 9462.41] 0.3 1.5
## 33 2.78 [ 1.29; 5.98] 5.1 2.5
## 34 1.32 [ 0.70; 2.51] 7.3 2.5
## 35 19.44 [ 4.90; 77.23] 1.6 2.3
## 36 234.00 [19.54; 2802.15] 0.5 1.8
## 37 20.35 [ 3.87; 107.10] 1.1 2.2
## Filled: 35 0.07 [ 0.02; 0.28] 1.6 2.3
## Filled: 37 0.07 [ 0.01; 0.36] 1.1 2.2
## Filled: 30 0.07 [ 0.02; 0.20] 2.5 2.4
## Filled: 14 0.06 [ 0.01; 0.31] 1.1 2.2
## Filled: 11 0.05 [ 0.01; 0.26] 1.0 2.1
## Filled: 22 0.01 [ 0.00; 0.29] 0.3 1.6
## Filled: 36 0.01 [ 0.00; 0.07] 0.5 1.8
## Filled: 32 0.00 [ 0.00; 0.07] 0.3 1.5
##
## Number of studies combined: k = 45 (with 8 added studies)
##
## OR 95%-CI z p-value
## Fixed effect model 1.17 [0.99; 1.40] 1.83 0.0675
## Random effects model 1.35 [0.74; 2.46] 0.99 0.3221
##
## Quantifying heterogeneity:
## tau^2 = 3.5926; H = 3.34 [3.00; 3.71]; I^2 = 91.0% [88.9%; 92.8%]
##
## Test of heterogeneity:
## Q d.f. p-value
## 489.76 44 < 0.0001
##
## Details on meta-analytical method:
## - Inverse variance method
## - DerSimonian-Laird estimator for tau^2
## - Trim-and-fill method to adjust for funnel plot asymmetry
library(metasens)
## Loading 'metasens' package (version 0.3-1).
c1 <- copas(ms1)
c1$TE.slope
## [1] 1.0346685 0.9998886 0.8001915 0.6007937 0.4012605 0.2003563
max(c1$seTE)
## [1] 1.578409
gamma0 <- min(c1$gamma0.range) + c1$x.slope*diff(c1$gamma0.range)
gamma1 <- min(c1$gamma1.range) + c1$y.slope*diff(c1$gamma1.range)
print(pnorm(gamma0 + gamma1/max(c1$seTE)), digits=2)
## [1] 0.97 0.77 0.62 0.50 0.40
print(summary(c1), digits=2)
## Summary of Copas selection model analysis:
##
## publprob OR 95%-CI pval.treat pval.rsb N.unpubl
## 1.00 2.81 [1.41; 5.60] 0.0032 0.0004 0
## 0.97 2.72 [1.43; 5.17] 0.0023 0.0005 0
## 0.77 2.23 [0.68; 7.24] 0.1835 0.0016 5
## 0.62 1.82 [0.70; 4.77] 0.2203 0.0045 12
## 0.50 1.49 [0.59; 3.81] 0.4011 0.0108 20
## 0.40 1.22 [0.48; 3.12] 0.6751 0.0231 31
##
## Copas model (adj)
## Random effects model 2.72 [1.50; 4.92] 0.0009
##
## Significance level for test of residual selection bias: 0.1
##
## Legend:
## publprob - Probability of publishing study with largest standard error
## pval.treat - P-value for hypothesis of overall treatment effect
## pval.rsb - P-value for hypothesis that no selection remains unexplained
## N.unpubl - Approximate number of unpublished studies suggested by model
c1
## Copas selection model analysis
##
## min max
## range of gamma0: -0.59 2.00
## range of gamma1: 0.00 0.53
##
## Largest standard error (SE): 1.5784
##
## Range of probability publishing trial with largest SE:
## min max
## 0.279 0.990
##
##
## Calculation of orthogonal line:
##
## level nobs adj.r.square slope se.slope
## -0.4 4 0.9998755 -3.196301 0.020595548
## -0.2 10 0.9999677 -3.241884 0.006145747
## 0.0 16 0.9999513 -3.279985 0.005910886
## 0.2 22 0.9998984 -3.338760 0.007345775
## 0.4 25 0.9998446 -3.391452 0.008629817
## 0.6 25 0.9998059 -3.457878 0.009834313
## 0.8 26 0.9997080 -3.577249 0.012227355
## 1.0 19 0.9993233 -4.113930 0.025232553
##
##
## Summary of Copas selection model analysis:
##
## publprob OR 95%-CI pval.treat pval.rsb N.unpubl
## 1.00 2.8142 [1.4135; 5.6027] 0.0032 0.0004 0
## 0.97 2.7180 [1.4292; 5.1689] 0.0023 0.0005 0
## 0.77 2.2260 [0.6846; 7.2376] 0.1835 0.0016 5
## 0.62 1.8236 [0.6978; 4.7655] 0.2203 0.0045 12
## 0.50 1.4937 [0.5854; 3.8112] 0.4011 0.0108 20
## 0.40 1.2218 [0.4787; 3.1185] 0.6751 0.0231 31
##
## Copas model (adj)
## Random effects model 2.7207 [1.5047; 4.9193] 0.0009
##
## Significance level for test of residual selection bias: 0.1
##
## Legend:
## publprob - Probability of publishing study with largest standard error
## pval.treat - P-value for hypothesis of overall treatment effect
## pval.rsb - P-value for hypothesis that no selection remains unexplained
## N.unpubl - Approximate number of unpublished studies suggested by model
c2 <- copas(ms1, gamma0.range=c(-1,2), gamma1.range=c(0,1))
print(summary(c2), digits=2)
## Summary of Copas selection model analysis:
##
## publprob OR 95%-CI pval.treat pval.rsb N.unpubl
## 1.00 2.81 [1.41; 5.60] 0.0032 0.0004 0
## 0.96 2.72 [1.42; 5.22] 0.0027 0.0005 0
## 0.72 2.23 [0.82; 6.06] 0.1155 0.0024 5
## 0.56 1.83 [0.78; 4.30] 0.1651 0.0078 12
## 0.44 1.50 [0.65; 3.49] 0.3444 0.0200 20
## 0.35 1.22 [0.52; 2.86] 0.6401 0.0435 31
## 0.28 1.00 [0.42; 2.35] 0.9955 0.0799 43
##
## Copas model (adj)
## Random effects model 2.72 [1.50; 4.92] 0.0009
##
## Significance level for test of residual selection bias: 0.1
##
## Legend:
## publprob - Probability of publishing study with largest standard error
## pval.treat - P-value for hypothesis of overall treatment effect
## pval.rsb - P-value for hypothesis that no selection remains unexplained
## N.unpubl - Approximate number of unpublished studies suggested by model
l1 <- limitmeta(ms1)
print(l1, digits=2)
## Results for individual studies (left: original data; right: shrunken estimates)
##
## OR 95%-CI OR 95%-CI
## 1 0.05 [ 0.01; 0.40] 0.01 [0.00; 0.05]
## 2 0.04 [ 0.01; 0.19] 0.02 [0.00; 0.06]
## 3 0.07 [ 0.02; 0.25] 0.03 [0.01; 0.11]
## 4 0.22 [ 0.06; 0.74] 0.09 [0.03; 0.30]
## 5 10.71 [ 1.05; 109.78] 0.38 [0.04; 3.86]
## 6 0.66 [ 0.22; 1.98] 0.30 [0.10; 0.90]
## 7 1.63 [ 0.50; 5.38] 0.61 [0.19; 2.02]
## 8 0.05 [ 0.02; 0.16] 0.03 [0.01; 0.08]
## 9 0.28 [ 0.12; 0.63] 0.18 [0.08; 0.41]
## 10 1.10 [ 0.40; 2.99] 0.55 [0.20; 1.50]
## 11 30.56 [ 5.31; 175.90] 3.21 [0.56; 18.46]
## 12 0.51 [ 0.23; 1.12] 0.34 [0.16; 0.74]
## 13 3.61 [ 1.23; 10.64] 1.54 [0.52; 4.53]
## 14 23.33 [ 4.50; 120.94] 3.15 [0.61; 16.31]
## 15 3.28 [ 1.21; 8.88] 1.58 [0.59; 4.28]
## 16 6.40 [ 1.89; 21.68] 2.12 [0.63; 7.20]
## 17 1.38 [ 0.58; 3.28] 0.81 [0.34; 1.93]
## 18 0.71 [ 0.37; 1.38] 0.53 [0.27; 1.02]
## 19 0.11 [ 0.05; 0.23] 0.07 [0.03; 0.16]
## 20 16.67 [ 4.07; 68.18] 3.70 [0.91; 15.16]
## 21 15.60 [ 4.32; 56.31] 4.38 [1.21; 15.83]
## 22 94.09 [ 4.77; 1854.23] 0.48 [0.02; 9.37]
## 23 11.50 [ 3.42; 38.71] 3.72 [1.11; 12.53]
## 24 1.83 [ 0.89; 3.77] 1.24 [0.60; 2.57]
## 25 5.06 [ 2.00; 12.80] 2.63 [1.04; 6.65]
## 26 2.04 [ 0.67; 6.21] 0.85 [0.28; 2.60]
## 27 8.59 [ 1.54; 48.01] 1.11 [0.20; 6.19]
## 28 12.83 [ 2.46; 66.92] 1.82 [0.35; 9.49]
## 29 12.30 [ 2.95; 51.26] 2.71 [0.65; 11.28]
## 30 20.62 [ 6.97; 61.00] 8.01 [2.71; 23.69]
## 31 11.19 [ 2.26; 55.33] 1.79 [0.36; 8.85]
## 32 429.00 [19.45; 9462.41] 1.12 [0.05; 24.61]
## 33 2.78 [ 1.29; 5.98] 1.79 [0.83; 3.86]
## 34 1.32 [ 0.70; 2.51] 0.99 [0.52; 1.87]
## 35 19.44 [ 4.90; 77.23] 4.51 [1.14; 17.93]
## 36 234.00 [19.54; 2802.15] 3.19 [0.27; 38.17]
## 37 20.35 [ 3.87; 107.10] 2.70 [0.51; 14.20]
##
## Result of limit meta-analysis:
##
## Random effects model OR 95%-CI z pval
## Adjusted estimate 0.74 [0.44; 1.24] -1.15 0.2496
## Unadjusted estimate 2.72 [1.50; 4.92] 3.31 0.0009
##
## Quantifying heterogeneity:
## tau^2 = 2.8702; I^2 = 90% [87.2%; 92.2%]; G^2 = 97.8%
##
## Test of heterogeneity:
## Q d.f. p-value
## 360.41 36 < 0.0001
##
## Test of small-study effects:
## Q-Q' d.f. p-value
## 73.02 1 < 0.0001
##
## Test of residual heterogeneity beyond small-study effects:
## Q' d.f. p-value
## 287.39 35 < 0.0001
##
## Details on adjustment method:
## - expectation (beta0)
funnel(l1)
