data14 <- read.csv("dataset14.csv")
data14
## author year mean.das se.das mean.haq se.haq
## 1 Anette 2004 -1.9 0.24 -0.33 0.16
## 2 Jane 2010 -1.3 0.11 -0.36 0.04
## 3 Bank 2002 -1.6 0.03 -0.42 0.01
## 4 Jiten 2015 -1.2 0.17 -0.25 0.06
## 5 Velo 2017 -1.1 0.28 -0.23 0.07
library(meta)
## Loading 'meta' package (version 4.8-4).
## Type 'help(meta)' for a brief overview.
m.das <- metagen(mean.das, se.das, data=data14, sm="MD", studlab=paste(author, year), comb.random=FALSE)
m.haq <- metagen(mean.haq, se.haq, data=data14, sm="MD", studlab=paste(author, year), comb.random=FALSE)
forest(m.das, hetstat=FALSE)
forest(m.haq, hetstat=FALSE)
library(mvmeta)
## This is mvmeta 0.4.7. For an overview type: help('mvmeta-package').
args(mvmeta)
## function (formula, S, data, subset, method = "reml", bscov = "unstr",
## model = TRUE, contrasts = NULL, offset, na.action, control = list())
## NULL
cor2cov <- function(sd1, sd2, cor) sd1*sd2*cor
theta <- cbind(data14$mean.das, data14$mean.haq)
dimnames(theta) <- list(data14$author, c("mean.das", "mean.haq"))
rho <- 0
S.arth <- cbind(data14$se.das^2, cor2cov(data14$se.das, data14$se.haq, rho), data14$se.haq^2)
dimnames(S.arth) <- list(data14$author, c("var.das", "cov", "var.haq"))
S.arth
## var.das cov var.haq
## Anette 0.0576 0 0.0256
## Jane 0.0121 0 0.0016
## Bank 0.0009 0 0.0001
## Jiten 0.0289 0 0.0036
## Velo 0.0784 0 0.0049
m.arth <- mvmeta(theta, S.arth, method="fixed")
print(summary(m.arth), digits=2)
## Call: mvmeta(formula = theta ~ 1, S = S.arth, method = "fixed")
##
## Multivariate fixed-effects meta-analysis
## Dimension: 2
##
## Fixed-effects coefficients
## Estimate Std. Error z Pr(>|z|) 95%ci.lb 95%ci.ub
## mean.das -1.57 0.03 -55.64 0.00 -1.62 -1.51 ***
## mean.haq -0.41 0.01 -43.14 0.00 -0.43 -0.39 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Multivariate Cochran Q-test for heterogeneity:
## Q = 32.98 (df = 8), p-value = 0.00
## I-square statistic = 75.7%
##
## 5 studies, 10 observations, 2 fixed and 0 random-effects parameters
## logLik AIC BIC
## -0.36 4.73 5.33
rho <- 0.9
S.arth2 <- cbind(data14$se.das^2, cor2cov(data14$se.das, data14$se.haq, rho), data14$se.haq^2)
dimnames(S.arth2) <- dimnames(S.arth)
rho <- -0.9
S.arth3 <- cbind(data14$se.das^2, cor2cov(data14$se.das, data14$se.haq, rho), data14$se.haq^2)
dimnames(S.arth3) <- dimnames(S.arth)
rho <- c(0.5, 0.2, 0.1, 0.6, 0.4)
S.arth4 <- cbind(data14$se.das^2, cor2cov(data14$se.das, data14$se.haq, rho), data14$se.haq^2)
dimnames(S.arth4) <- dimnames(S.arth)
m.arth2 <- mvmeta(theta, S.arth2, method="fixed")
m.arth3 <- mvmeta(theta, S.arth3, method="fixed")
m.arth4 <- mvmeta(theta, S.arth4, method="fixed")
round(coef(m.arth2), 3)
## mean.das.(Intercept) mean.haq.(Intercept)
## -1.575 -0.411
vcov(m.arth2)
## mean.das.(Intercept) mean.haq.(Intercept)
## mean.das.(Intercept) 0.0007777489 2.346058e-04
## mean.haq.(Intercept) 0.0002346058 8.781479e-05
round(sqrt(diag(vcov(m.arth2))), 4)
## mean.das.(Intercept) mean.haq.(Intercept)
## 0.0279 0.0094
rho <- 0.9
sample.sizes <- c(45,169,820,102,55)
S.arth.common <- matrix(1/sample.sizes, ncol=1) %*% matrix(c(2.146, rho*sqrt(2.146*0.352), 0.352), nrow=1)
dimnames(S.arth.common) <- list(data14$author, c("var.das", "cov", "var.haq"))
round(S.arth.common, 4)
## var.das cov var.haq
## Anette 0.0477 0.0174 0.0078
## Jane 0.0127 0.0046 0.0021
## Bank 0.0026 0.0010 0.0004
## Jiten 0.0210 0.0077 0.0035
## Velo 0.0390 0.0142 0.0064
m.arth.common <- mvmeta(theta, S.arth.common, method="fixed")
print(summary(m.arth.common), digits=2)
## Call: mvmeta(formula = theta ~ 1, S = S.arth.common, method = "fixed")
##
## Multivariate fixed-effects meta-analysis
## Dimension: 2
##
## Fixed-effects coefficients
## Estimate Std. Error z Pr(>|z|) 95%ci.lb 95%ci.ub
## mean.das -1.51 0.04 -35.61 0.00 -1.59 -1.43 ***
## mean.haq -0.38 0.02 -22.38 0.00 -0.42 -0.35 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Multivariate Cochran Q-test for heterogeneity:
## Q = 52.31 (df = 8), p-value = 0.00
## I-square statistic = 84.7%
##
## 5 studies, 10 observations, 2 fixed and 0 random-effects parameters
## logLik AIC BIC
## -6.21 16.43 17.03
funnel(m.das, xlab="Mean difference from baseline")
sel <- m.das$studlab %in% c("Anette 2004", "Bank 2002", "Jiten 2015")
text(m.das$TE[sel]-0.015, m.das$seTE[sel]-0.002, m.das$studlab[sel], adj=1, cex=0.8)
text(m.das$TE[!sel]+0.015, m.das$seTE[!sel]-0.002, m.das$studlab[!sel], adj=0, cex=0.8)
funnel(m.haq, xlab="Mean difference from baseline")
text(m.haq$TE-0.005, m.haq$seTE-0.002, m.haq$studlab, adj=1, cex=0.8)
library(ellipse)
plot(data14$mean.das, data14$mean.haq, xlim=c(-2.5, -0.5), ylim=c(-0.7, 0), pch="+",
xlab="Mean difference in DAS-28 score from baseline",
ylab="Mean difference in HAQ score from baseline")
rho <- 0.25
with(data14,
for (i in seq(along=mean.das)){
S <- matrix(c(se.das[i]^2,
se.das[i]*se.haq[i]*rho,
se.das[i]*se.haq[i]*rho,
se.haq[i]^2), byrow=TRUE, ncol=2)
lines(ellipse(S, centre=c(mean.das[i], mean.haq[i]), level=0.95 ), col="grey")
})
studlab <- paste(data14$author, data14$year)
velo <- data14$author=="Velo"
studlab[velo] <- paste(data14$author[velo], "\n", data14$year[velo])
text(data14$mean.das+0.015, data14$mean.haq, studlab, adj=0)
theta.miss <- theta
selnava <- rownames(theta.miss)=="Jiten"
theta.miss[selnava, "mean.das"] <- NA
theta.miss
## mean.das mean.haq
## Anette -1.9 -0.33
## Jane -1.3 -0.36
## Bank -1.6 -0.42
## Jiten NA -0.25
## Velo -1.1 -0.23
m.arth.miss <- mvmeta(theta.miss, S.arth, method="fixed")
print(summary(m.arth.miss), digits=2)
## Call: mvmeta(formula = theta.miss ~ 1, S = S.arth, method = "fixed")
##
## Multivariate fixed-effects meta-analysis
## Dimension: 2
##
## Fixed-effects coefficients
## Estimate Std. Error z Pr(>|z|) 95%ci.lb 95%ci.ub
## mean.das -1.58 0.03 -55.23 0.00 -1.63 -1.52 ***
## mean.haq -0.41 0.01 -43.14 0.00 -0.43 -0.39 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Multivariate Cochran Q-test for heterogeneity:
## Q = 28.15 (df = 7), p-value = 0.00
## I-square statistic = 75.1%
##
## 5 studies, 9 observations, 2 fixed and 0 random-effects parameters
## logLik AIC BIC
## 1.20 1.61 2.00
rho <- 0.9
S.arth.r <- cbind(data14$se.das^2, cor2cov(data14$se.das, data14$se.haq, rho), data14$se.haq^2)
m.arth.reml <- mvmeta(theta, S.arth.r, method="reml")
m.arth.ml <- mvmeta(theta, S.arth.r, method="ml")
m.arth.mm <- mvmeta(theta, S.arth.r, method="mm")
print(summary(m.arth.reml), digits=2)
## Call: mvmeta(formula = theta ~ 1, S = S.arth.r, method = "reml")
##
## Multivariate random-effects meta-analysis
## Dimension: 2
## Estimation method: REML
##
## Fixed-effects coefficients
## Estimate Std. Error z Pr(>|z|) 95%ci.lb 95%ci.ub
## mean.das -1.48 0.16 -9.10 0.00 -1.80 -1.16 ***
## mean.haq -0.35 0.04 -7.87 0.00 -0.44 -0.26 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Between-study random-effects (co)variance components
## Structure: General positive-definite
## Std. Dev Corr
## mean.das 0.33 mean.das
## mean.haq 0.08 0.83
##
## Multivariate Cochran Q-test for heterogeneity:
## Q = 45.77 (df = 8), p-value = 0.00
## I-square statistic = 82.5%
##
## 5 studies, 10 observations, 2 fixed and 3 random-effects parameters
## logLik AIC BIC
## 2.22 5.56 5.96
blup(m.arth.mm)
## mean.das mean.haq
## Anette -1.938147 -0.4474356
## Jane -1.297128 -0.3493406
## Bank -1.596638 -0.4183855
## Jiten -1.324239 -0.3040035
## Velo -1.304442 -0.2975754
blup(m.arth.mm) - fitted(m.arth.mm)
## mean.das mean.haq
## Anette -0.4460285 -0.08408749
## Jane 0.1949909 0.01400755
## Bank -0.1045196 -0.05503737
## Jiten 0.1678800 0.05934457
## Velo 0.1876772 0.06577274