Merge pull request #3 from kylehamilton/tests - Import metafor tests

For data that was in metafor I moved the tests over to metadat and changed the tests I wrote for the other datasets to match the correct format

tests/testthat/test-dat.begg1989.R → tests/testthat/test_analysis_example_begg1989.r

@@ -1,4 +1,4 @@
-context("test-dat.begg1989.R")
+context("Checking analysis example: dat.begg1989")
 library(metadat)
 library(metafor)
 

tests/testthat/test_analysis_example_berkey1998.r

@@ -0,0 +1,74 @@
+### library(metafor); library(testthat); Sys.setenv(NOT_CRAN="true")
+
+### see also: http://www.metafor-project.org/doku.php/analyses:berkey1998
+
+library(metafor)
+source("tolerances.r") # read in tolerances
+
+context("Checking analysis example: berkey1998")
+
+### load data
+dat <- metadat::dat.berkey1998
+
+### construct variance-covariance matrix of the observed outcomes
+V <- bldiag(lapply(split(dat[,c("v1i", "v2i")], dat$trial), as.matrix))
+
+test_that("results are correct for the multiple outcomes random-effects model.", {
+
+   ### multiple outcomes random-effects model (with ML estimation)
+   res <- rma.mv(yi, V, mods = ~ outcome - 1, random = ~ outcome | trial, struct="UN", data=dat, method="ML")
+   out <- capture.output(print(res)) ### so that print.rma.mv() is run (at least once)
+
+   ### (results for this model not given in paper)
+   expect_equivalent(coef(res), c(-0.3379, 0.3448), tolerance=.tol[["coef"]])
+   expect_equivalent(res$se, c(0.0798, 0.0495), tolerance=.tol[["se"]])
+   expect_equivalent(res$tau2, c(0.0261, 0.0070), tolerance=.tol[["var"]])
+   expect_equivalent(res$rho, 0.6992, tolerance=.tol[["cor"]])
+
+})
+
+test_that("results are correct for the multiple outcomes mixed-effects (meta-regression) model.", {
+
+   ### multiple outcomes mixed-effects (meta-regression) model (with ML estimation)
+   res <- rma.mv(yi, V, mods = ~ outcome + outcome:I(year - 1983) - 1, random = ~ outcome | trial, struct="UN", data=dat, method="ML")
+
+   ### compare with results on page 2545 (Table II)
+   expect_equivalent(coef(res), c(-0.3351, 0.3479, -0.0108, 0.0010), tolerance=.tol[["coef"]])
+   expect_equivalent(res$se, c(0.0787, 0.0520, 0.0243, 0.0154), tolerance=.tol[["se"]])
+   expect_equivalent(res$tau2, c(0.0250, 0.0080), tolerance=.tol[["var"]])
+   expect_equivalent(res$rho, 0.6587, tolerance=.tol[["cor"]])
+
+   ### compute the covariance
+   tmp <- res$rho*sqrt(res$tau2[1]*res$tau2[2])
+   expect_equivalent(tmp, 0.0093, tolerance=.tol[["cov"]])
+
+   ### test the difference in slopes
+   res <- rma.mv(yi, V, mods = ~ outcome*I(year - 1983) - 1, random = ~ outcome | trial, struct="UN", data=dat, method="ML")
+
+   ### (results for this model not given in paper)
+   expect_equivalent(coef(res), c(-0.3351, 0.3479, -0.0108, 0.0118), tolerance=.tol[["coef"]])
+   expect_equivalent(res$se, c(0.0787, 0.0520, 0.0243, 0.0199), tolerance=.tol[["se"]])
+   expect_equivalent(res$pval, c(0.0000, 0.0000, 0.6563, 0.5534), tolerance=.tol[["pval"]])
+
+})
+
+test_that("results are correct when testing var-cov structures against each other with LRTs.", {
+
+   ### test whether the amount of heterogeneity is the same in the two outcomes
+   res1 <- rma.mv(yi, V, mods = ~ outcome - 1, random = ~ outcome | trial, struct="UN", data=dat, method="ML")
+   res0 <- rma.mv(yi, V, mods = ~ outcome - 1, random = ~ outcome | trial, struct="CS", data=dat, method="ML")
+   tmp <- anova(res0, res1)
+   out <- capture.output(print(tmp)) ### so that print.anova.rma() is run (at least once)
+
+   ### (results for this not given in paper)
+   expect_equivalent(tmp$pval, 0.2597, tolerance=.tol[["pval"]])
+
+   ### test the correlation among the true effects
+   res1 <- rma.mv(yi, V, mods = ~ outcome - 1, random = ~ outcome | trial, struct="UN", data=dat, method="ML")
+   res0 <- rma.mv(yi, V, mods = ~ outcome - 1, random = ~ outcome | trial, struct="UN", data=dat, method="ML", rho=0)
+   tmp <- anova(res0, res1)
+
+   ### (results for this not given in paper)
+   expect_equivalent(tmp$pval, 0.2452, tolerance=.tol[["pval"]])
+
+})

tests/testthat/test_analysis_example_ishak2007.r

@@ -0,0 +1,134 @@
+### library(metafor); library(testthat); Sys.setenv(NOT_CRAN="true")
+
+context("Checking analysis example: ishak2007")
+
+library(metafor)
+source("tolerances.r") # read in tolerances
+
+### load dataset
+dat <- metadat::dat.ishak2007
+
+### create long format dataset
+dat.long <- reshape(dat, direction="long", idvar="study", v.names=c("yi","vi"),
+                    varying=list(c(2,4,6,8), c(3,5,7,9)))
+dat.long <- dat.long[order(dat.long$study, dat.long$time),]
+rownames(dat.long) <- 1:nrow(dat.long)
+
+### remove missing measurement occasions from dat.long
+is.miss  <- is.na(dat.long$yi)
+dat.long <- dat.long[!is.miss,]
+
+### construct the full (block diagonal) V matrix with an AR(1) structure
+rho.within <- .97 ### value as estimated by Ishak et al. (2007)
+V <- lapply(split(with(dat, cbind(v1i, v2i, v3i, v4i)), dat$study), diag)
+V <- lapply(V, function(v) sqrt(v) %*% toeplitz(ARMAacf(ar=rho.within, lag.max=3)) %*% sqrt(v))
+V <- bldiag(V)
+V <- V[!is.miss,!is.miss] ### remove missing measurement occasions from V
+
+test_that("results are correct for diag(V) and struct='DIAG'.", {
+
+   res <- rma.mv(yi, diag(V), mods = ~ factor(time) - 1, random = ~ factor(time) | study,
+                 struct = "DIAG", data = dat.long)
+
+   ### Table 1, column "Time-specific (Independence)"
+   expect_equivalent(coef(res), c(-24.8686, -27.4728, -28.5239, -24.1415), tolerance=.tol[["coef"]])
+   expect_equivalent(res$tau2, c(23.0537, 27.8113, 27.6767, 29.9405), tolerance=.tol[["var"]])
+
+})
+
+test_that("results are correct for diag(V) and random study effects.", {
+
+   res <- rma.mv(yi, diag(V), mods = ~ factor(time) - 1, random = ~ 1 | study, data = dat.long)
+
+   ### Table 1, column "Random study effects"
+   expect_equivalent(coef(res), c(-26.2127, -27.1916, -28.5464, -25.6339), tolerance=.tol[["coef"]])
+   expect_equivalent(res$sigma2, 26.6829, tolerance=.tol[["var"]])
+
+})
+
+test_that("results are correct for diag(V) and struct='ID'.", {
+
+   res <- rma.mv(yi, diag(V), mods = ~ factor(time) - 1, random = ~ factor(time) | study,
+                 struct = "ID", data = dat.long)
+
+   ### not in paper
+   expect_equivalent(coef(res), c(-24.8792, -27.4670, -28.5185, -24.1502), tolerance=.tol[["coef"]])
+   expect_equivalent(res$tau2, 26.6847, tolerance=.tol[["var"]])
+
+})
+
+test_that("results are correct for diag(V) and struct='HAR'.", {
+
+   res <- rma.mv(yi, diag(V), mods = ~ factor(time) - 1, random = ~ time | study,
+                 struct = "HAR", data = dat.long)
+
+   ### Table 1, column "Correlated random time effects"
+   expect_equivalent(coef(res), c(-25.9578, -27.3100, -28.5543, -25.7923), tolerance=.tol[["coef"]]) # -27.5 in Table vs -27.3
+   expect_equivalent(res$tau2, c(20.3185, 35.9720, 26.4233, 30.1298), tolerance=.tol[["var"]]) # 20.4 in Table vs 20.3
+   expect_equivalent(res$rho, 1.0000, tolerance=.tol[["cor"]])
+
+})
+
+test_that("results are correct for struct='HAR'.", {
+
+   res <- rma.mv(yi, V, mods = ~ factor(time) - 1, random = ~ time | study,
+                 struct = "HAR", data = dat.long)
+
+   ### Table 1, column "Multivariate model"
+   expect_equivalent(coef(res), c(-25.9047, -27.4608, -28.6559, -26.4934), tolerance=.tol[["coef"]])
+   expect_equivalent(res$tau2, c(22.7258, 33.7295, 26.1426, 31.1803), tolerance=.tol[["var"]]) # 22.6 in Table vs 22.7; 31.1 in Table vs 31.2
+   expect_equivalent(res$rho, 0.8832, tolerance=.tol[["cor"]])
+
+})
+
+test_that("results are correct for struct='AR'.", {
+
+   res <- rma.mv(yi, V, mods = ~ factor(time) - 1, random = ~ time | study,
+                 struct = "AR", data = dat.long)
+
+   ### not in paper
+   expect_equivalent(coef(res), c(-25.9418, -27.3937, -28.7054, -26.3970), tolerance=.tol[["coef"]])
+   expect_equivalent(res$tau2, 26.6874, tolerance=.tol[["var"]])
+   expect_equivalent(res$rho, 0.8656, tolerance=.tol[["cor"]])
+
+})
+
+test_that("results are correct for struct='HCS'.", {
+
+   res <- rma.mv(yi, V, mods = ~ factor(time) - 1, random = ~ factor(time) | study,
+                 struct = "HCS", data = dat.long)
+
+   ### not in paper
+   expect_equivalent(coef(res), c(-25.8814, -27.3293, -28.6510, -26.6631), tolerance=.tol[["coef"]])
+   expect_equivalent(res$tau2, c(20.8629, 32.7429, 27.6593, 32.1908), tolerance=.tol[["var"]])
+
+})
+
+test_that("results are correct for struct='CAR'.", {
+
+   res <- rma.mv(yi, V, mods = ~ factor(time) - 1, random = ~ time | study,
+                 struct = "CAR", data = dat.long)
+
+   ### not in paper
+   expect_equivalent(coef(res), c(-25.9418, -27.3937, -28.7054, -26.3970), tolerance=.tol[["coef"]])
+   expect_equivalent(res$tau2, 26.6875, tolerance=.tol[["var"]])
+   expect_equivalent(res$rho, 0.8656, tolerance=.tol[["cor"]])
+
+})
+
+test_that("results are correct for struct='CAR' with unequally spaced time points.", {
+
+   dat.long$time[dat.long$time == 4] <- 24/3
+   dat.long$time[dat.long$time == 3] <- 12/3
+   dat.long$time[dat.long$time == 2] <-  6/3
+   dat.long$time[dat.long$time == 1] <-  3/3
+
+   res <- rma.mv(yi, V, mods = ~ factor(time) - 1, random = ~ time | study,
+                 struct = "CAR", data = dat.long)
+
+   ### not in paper
+   expect_equivalent(coef(res), c(-26.0293, -27.3838, -28.7339, -26.0515), tolerance=.tol[["coef"]])
+   expect_equivalent(res$tau2, 26.9825, tolerance=.tol[["var"]])
+   expect_equivalent(res$rho, 0.9171, tolerance=.tol[["cor"]])
+
+})