data1 <- read.csv("data1.csv")
str(data1)
## 'data.frame': 17 obs. of 8 variables:
## $ author: Factor w/ 17 levels "Aishwarya","Anarkalli",..: 4 6 5 2 3 8 7 13 14 16 ...
## $ year : int 1998 1999 2007 1993 1994 1994 1997 2006 2008 1997 ...
## $ Ne : int 23 21 14 21 18 13 11 14 19 25 ...
## $ Me : num 23.5 17.5 22.3 16.9 18.8 ...
## $ Se : num 11.7 15.3 14.5 12.9 11.5 19.1 16.5 19.2 12.6 14.2 ...
## $ Nc : int 16 25 13 16 18 23 11 14 15 16 ...
## $ Mc : num 21.3 21.4 31.4 23.6 29.9 ...
## $ Sc : num 22.2 17.3 21.3 16.5 18.3 ...
data1
## author year Ne Me Se Nc Mc Sc
## 1 Bill 1998 23 23.54 11.70 16 21.33 22.22
## 2 Dick 1999 21 17.50 15.30 25 21.43 17.32
## 3 Chawla 2007 14 22.30 14.50 13 31.43 21.32
## 4 Anarkalli 1993 21 16.90 12.90 16 23.65 16.54
## 5 Biaggi 1994 18 18.80 11.50 18 29.88 18.32
## 6 Lundwig 1994 13 15.30 19.10 23 35.43 17.65
## 7 Dulquer 1997 11 14.90 16.50 11 29.87 18.43
## 8 Nivin 2006 14 24.56 19.20 14 45.44 16.32
## 9 Pauly 2008 19 18.60 12.60 15 43.67 18.43
## 10 Salmaan 1997 25 13.56 14.20 16 29.65 12.21
## 11 Aishwarya 1994 11 14.44 11.54 16 41.27 20.01
## 12 Rajesh 1993 14 19.32 9.21 21 38.65 12.07
## 13 Meena 1995 12 14.00 11.65 17 31.49 15.32
## 14 Mohanlal 1985 13 13.30 11.54 21 38.32 14.27
## 15 Mammooty 1999 11 13.50 19.43 10 45.67 18.54
## 16 Manju 2000 10 19.50 8.54 12 43.21 11.54
## 17 Warrier 2003 9 13.10 9.67 12 32.10 7.86
library(meta)
## Loading 'meta' package (version 4.8-4).
## Type 'help(meta)' for a brief overview.
m <- metacont(Ne, Me, Se, Nc, Mc, Sc, studlab=paste(author, year), data=data1)
mdata <- as.data.frame(m)
mdata$studlab <- as.character(mdata$studlab)
data12 <- mdata[rank(mdata$seTE)<=7, c(7,1:6)]
names(data12) <- c("study", "Ne", "Me", "Se", "Nc", "Mc", "Sc")
data12
## study Ne Me Se Nc Mc Sc
## 2 Dick 1999 21 17.50 15.30 25 21.43 17.32
## 4 Anarkalli 1993 21 16.90 12.90 16 23.65 16.54
## 10 Salmaan 1997 25 13.56 14.20 16 29.65 12.21
## 12 Rajesh 1993 14 19.32 9.21 21 38.65 12.07
## 14 Mohanlal 1985 13 13.30 11.54 21 38.32 14.27
## 16 Manju 2000 10 19.50 8.54 12 43.21 11.54
## 17 Warrier 2003 9 13.10 9.67 12 32.10 7.86
data12$Nem <- rep(0, length(data12$Ne))
data12$Nem[data12$study=="Anarkalli 1993"] <- 8
data12$Nem[data12$study=="Salmaan 1997"] <- 5
data12$Ncm <- data12$Nem
data12$Neo <- data12$Ne - data12$Nem
data12$Nco <- data12$Nc - data12$Ncm
mm1 <- metacont(Neo, Me, Se, Nco, Mc, Sc, data=data12, studlab=study)
data12$TE <- mm1$TE
data12$seTE <- mm1$seTE
data12[, c("study", "TE", "seTE", "Neo", "Nco")]
## study TE seTE Neo Nco
## 2 Dick 1999 -3.93 4.811075 21 25
## 4 Anarkalli 1993 -6.75 6.855452 13 8
## 10 Salmaan 1997 -16.09 4.861594 20 11
## 12 Rajesh 1993 -19.33 3.605030 14 21
## 14 Mohanlal 1985 -25.02 4.465509 13 21
## 16 Manju 2000 -23.71 4.288449 10 12
## 17 Warrier 2003 -19.00 3.941850 9 12
semiss <- function(se, n.e, p.e, n.c, p.c,
mu.e, nu.e, mu.c, nu.c, rho){
V1 <- function(p.e, p.c, nu.e, nu.c, rho)
nu.e^2*p.e^2 + nu.c^2*p.c^2 - 2*rho*nu.e*nu.c*p.e*p.c
V2 <- function(n.e, p.e, n.c, p.c, mu.e, nu.e, mu.c, nu.c)
(mu.e^2+nu.e^2)*p.e*(1-p.e)/n.e +
(mu.c^2+nu.c^2)*p.c*(1-p.c)/n.c
sqrt(se^2 +
V1(p.e, p.c, nu.e, nu.c, rho) +
V2(n.e, p.e, n.e, p.c,
mu.e, nu.e, mu.c, nu.c))
}
# 1. Define parameters
mu.e <- mu.c <- 8
nu.e <- nu.c <- 0
rho <- 0
# 2. Calculate proportion missing in each study arm
data12$Pe <- with(data12, Pe <- (Ne - Neo) / Ne)
data12$Pc <- with(data12, Pc <- (Nc - Nco) / Nc)
# 3. Calculate the mean effect for each study under strategy 1
data12$TEs1 <- with(data12, TE + mu.e*Pe - mu.c*Pc)
data12$seTEs1 <- with(data12,
semiss(seTE, Neo, Pe, Nco, Pc,
mu.e, nu.e, mu.c, nu.c, rho))
# 4. Create indicator for studies with missings
selmiss <- data12$study %in% c("Anarkalli 1993", "Salmaan 1997")
data12[selmiss, c("study", "TE", "seTE", "TEs1", "seTEs1")]
## study TE seTE TEs1 seTEs1
## 4 Anarkalli 1993 -6.75 6.855452 -7.702381 7.027730
## 10 Salmaan 1997 -16.09 4.861594 -16.990000 4.983433
mm1.s1 <- metagen(TEs1, seTEs1, data=data12, studlab=study, comb.fixed=FALSE)
print(summary(mm1.s1), digits=2)
## Number of studies combined: k = 7
##
## 95%-CI z p-value
## Random effects model -17.24 [-22.63; -11.85] -6.27 < 0.0001
##
## Quantifying heterogeneity:
## tau^2 = 31.2408; H = 1.59 [1.05; 2.41]; I^2 = 60.5% [9.4%; 82.8%]
##
## Test of heterogeneity:
## Q d.f. p-value
## 15.18 6 0.0189
##
## Details on meta-analytical method:
## - Inverse variance method
## - DerSimonian-Laird estimator for tau^2
mu.e <- 8
mu.c <- -mu.e
nu.e <- nu.c <- 0
rho <- 0
data12$TEs2 <- with(data12, TE + mu.e*Pe - mu.c*Pc)
data12$seTEs2 <- with(data12, semiss(seTE, Neo, Pe, Nco, Pc, mu.e, nu.e, mu.c, nu.c, rho))
data12[selmiss, c("study", "TE", "seTE", "TEs2", "seTEs2")]
## study TE seTE TEs2 seTEs2
## 4 Anarkalli 1993 -6.75 6.855452 0.297619 7.027730
## 10 Salmaan 1997 -16.09 4.861594 -11.990000 4.983433
mm1.s2 <- metagen(TEs2, seTEs2, data=data12, studlab=study, comb.fixed=FALSE)
print(summary(mm1.s2), digits=2)
## Number of studies combined: k = 7
##
## 95%-CI z p-value
## Random effects model -15.56 [-21.86; -9.26] -4.84 < 0.0001
##
## Quantifying heterogeneity:
## tau^2 = 50.3247; H = 1.86 [1.26; 2.75]; I^2 = 71.1% [37.2%; 86.7%]
##
## Test of heterogeneity:
## Q d.f. p-value
## 20.79 6 0.0020
##
## Details on meta-analytical method:
## - Inverse variance method
## - DerSimonian-Laird estimator for tau^2
mu.e <- mu.c <- 8
nu.e <- nu.c <- sqrt(5)
rho <- 1
# 2. Calculate the mean effect for each study under strategy 2
data12$TEs3 <- with(data12, TE + mu.e*Pe - mu.c*Pc)
data12$seTEs3 <- with(data12, semiss(seTE, Neo, Pe, Nco, Pc, mu.e, nu.e, mu.c, nu.c, rho))
data12[selmiss, c("study", "TE", "seTE", "TEs3", "seTEs3")]
## study TE seTE TEs3 seTEs3
## 4 Anarkalli 1993 -6.75 6.855452 -7.702381 7.046042
## 10 Salmaan 1997 -16.09 4.861594 -16.990000 4.999159
mm1.s3 <- metagen(TEs3, seTEs3, data=data12, studlab=study, comb.fixed=FALSE)
print(summary(mm1.s3), digits=2)
## Number of studies combined: k = 7
##
## 95%-CI z p-value
## Random effects model -17.24 [-22.63; -11.85] -6.27 < 0.0001
##
## Quantifying heterogeneity:
## tau^2 = 31.2394; H = 1.59 [1.05; 2.41]; I^2 = 60.5% [9.4%; 82.7%]
##
## Test of heterogeneity:
## Q d.f. p-value
## 15.17 6 0.0190
##
## Details on meta-analytical method:
## - Inverse variance method
## - DerSimonian-Laird estimator for tau^2
mu.e <- 8
mu.c <- -mu.e
nu.e <- nu.c <- sqrt(5)
rho <- 0
data12$TEs4 <- with(data12, TE + mu.e*Pe - mu.c*Pc)
data12$seTEs4 <- with(data12, semiss(seTE, Neo, Pe, Nco, Pc, mu.e, nu.e, mu.c, nu.c, rho))
data12[selmiss, c("study", "TE", "seTE", "TEs4", "seTEs4")]
## study TE seTE TEs4 seTEs4
## 4 Anarkalli 1993 -6.75 6.855452 0.297619 7.179935
## 10 Salmaan 1997 -16.09 4.861594 -11.990000 5.061284
mm1.s4 <- metagen(TEs4, seTEs4, data=data12, studlab=study, comb.fixed=FALSE)
print(summary(mm1.s4), digits=2)
## Number of studies combined: k = 7
##
## 95%-CI z p-value
## Random effects model -15.61 [-21.90; -9.32] -4.87 < 0.0001
##
## Quantifying heterogeneity:
## tau^2 = 49.7118; H = 1.85 [1.25; 2.73]; I^2 = 70.7% [36.2%; 86.6%]
##
## Test of heterogeneity:
## Q d.f. p-value
## 20.51 6 0.0022
##
## Details on meta-analytical method:
## - Inverse variance method
## - DerSimonian-Laird estimator for tau^2
data13 <- mdata[rank(mdata$seTE)<=7, c(7, 1:6)]
names(data13) <- c("study", "Ne", "Me", "Se", "Nc", "Mc", "Sc")
data13.noManju <- data13
data13.noManju$Se[data13.noManju$study=="Manju 2000"] <- NA
data13.noManju$Sc[data13.noManju$study=="Manju 2000"] <- NA
data13.noManju
## study Ne Me Se Nc Mc Sc
## 2 Dick 1999 21 17.50 15.30 25 21.43 17.32
## 4 Anarkalli 1993 21 16.90 12.90 16 23.65 16.54
## 10 Salmaan 1997 25 13.56 14.20 16 29.65 12.21
## 12 Rajesh 1993 14 19.32 9.21 21 38.65 12.07
## 14 Mohanlal 1985 13 13.30 11.54 21 38.32 14.27
## 16 Manju 2000 10 19.50 NA 12 43.21 NA
## 17 Warrier 2003 9 13.10 9.67 12 32.10 7.86
f.tests <- with(data13.noManju, c(pf(Se^2/Se[1]^2, Ne-1, Ne[1]-1),
pf(Sc^2/Sc[1]^2, Nc-1, Nc[1]-1)))
f.tests <- matrix(f.tests, ncol=2)[-1,]
rownames(f.tests) <- data13.noManju$study[-1]
colnames(f.tests) <- c("pval.e", "pval.c")
round(f.tests, 2)
## pval.e pval.c
## Anarkalli 1993 0.23 0.44
## Salmaan 1997 0.36 0.08
## Rajesh 1993 0.03 0.05
## Mohanlal 1985 0.16 0.19
## Manju 2000 NA NA
## Warrier 2003 0.09 0.00
miss <- data13.noManju
M <- 100
imp.Manju <- data.frame(seTE=rep(NA, M),
TE.fixed=NA, seTE.fixed=NA,
TE.random=NA, seTE.random=NA,
tau=NA)
set.seed(123)
selimp <- !(miss$study %in% c("Manju 2000", "Warrier 2003"))
selManju <- miss$study=="Manju 2000"
S2.e <- with(miss, sum((Ne[selimp]-1)*Se[selimp]^2))
S2.c <- with(miss, sum((Nc[selimp]-1)*Sc[selimp]^2))
df.e <- sum(miss$Ne[selimp]-1)
df.c <- sum(miss$Nc[selimp]-1)
miss$Se.Manju <- miss$Se
miss$Sc.Manju <- miss$Sc
for (m in 1:M) {
sigma2.e <- S2.e / rchisq(1, df=df.e)
sigma2.c <- S2.c / rchisq(1, df=df.c)
sd.e.Manju <- sigma2.e * rchisq(1, df=miss$Ne[selManju]-1)/(miss$Ne[selManju]-1)
sd.c.Manju <- sigma2.c * rchisq(1, df=miss$Nc[selManju]-1)/(miss$Nc[selManju]-1)
imp.Manju$seTE[m] <- sqrt(sd.e.Manju/miss$Ne[selManju] + sd.c.Manju/miss$Nc[selManju])
miss$Se.Manju[selManju] <- sqrt(sd.e.Manju)
miss$Sc.Manju[selManju] <- sqrt(sd.c.Manju)
m.Manju <- metacont(n.e=Ne, mean.e=Me, sd.e=Se.Manju, n.c=Nc, mean.c=Mc, sd.c=Sc.Manju, data=miss)
imp.Manju$TE.fixed[m] <- m.Manju$TE.fixed
imp.Manju$seTE.fixed[m] <- m.Manju$seTE.fixed
imp.Manju$TE.random[m] <- m.Manju$TE.random
imp.Manju$seTE.random[m] <- m.Manju$seTE.random
imp.Manju$tau[m] <- m.Manju$tau
}
s2.b.fixed <- var(imp.Manju$TE.fixed)
s2.b.random <- var(imp.Manju$TE.random)
s2.w.fixed <- mean(imp.Manju$seTE.fixed^2)
s2.w.random <- mean(imp.Manju$seTE.random^2)
M <- length(imp.Manju$TE.fixed)
TE.fixed.imp <- mean(imp.Manju$TE.fixed)
seTE.fixed.imp <- sqrt(var(imp.Manju$TE.fixed)*(1 + 1/M) + mean(imp.Manju$seTE.fixed^2))
TE.random.imp <- mean(imp.Manju$TE.random)
seTE.random.imp <- sqrt(var(imp.Manju$TE.random)*(1 + 1/M) + mean(imp.Manju$seTE.random^2))
df.fixed <- (M-1)*(1 + s2.w.fixed /((1+1/M)*s2.b.fixed) )^2
df.random <- (M-1)*(1 + s2.w.random/((1+1/M)*s2.b.random))^2
df.fixed
## [1] 322049.9
df.random
## [1] 15881824
round(unlist(ci(TE.fixed.imp, seTE.fixed.imp))[1:5], 2)
## TE seTE lower upper z
## -16.70 1.67 -19.98 -13.42 -9.98
round(unlist(ci(TE.random.imp, seTE.random.imp))[1:5], 2)
## TE seTE lower upper z
## -16.35 2.83 -21.90 -10.80 -5.78
mm2 <- metacont(Ne, Me, Se, Nc, Mc, Sc, studlab=study, data=data13)
TE.Manju <- mm2$TE[mm2$studlab=="Manju 2000"]
seTE.Manju <- mm2$seTE[mm2$studlab=="Manju 2000"]
funnel(mm2, xlim=c(-30, 0), ylim=c(max(imp.Manju$seTE), 0))
text(TE.Manju-0.25, seTE.Manju-0.1, "Manu 2000", cex=0.8, adj=1)
set.seed(456)
points(jitter(rep(TE.Manju, length(imp.Manju$seTE)), 0.5),
imp.Manju$seTE, pch=2, cex=0.2)
