dat <- readxl::read_xlsx("DN_2020_Shoulder.xlsx")

require(meta)
settings.meta(digits = 2, digits.tau2 = 2, digits.tau = 2,
              print.I2 = FALSE)


dat1 <- dat[-c(6, 7, 8, 9), ] # choose subset


m0 <- metacont(n.e=Ne, mean.e=Me, sd.e=Se,
               n.c=Nc, mean.c=Mc, sd.c=Sc,
               data = dat1,
               sm = "SMD",
               comb.fixed = FALSE,
               method.tau = "REML")
##
mb0 <- metabind(m0,
                update(m0, hakn = TRUE, adhoc.hakn = ""),
                update(m0, hakn = TRUE, adhoc.hakn = "ci"),
                update(m0, hakn = TRUE, adhoc.hakn = "se"),
                update(m0, hakn = TRUE, adhoc.hakn = "iqwig6"),
                name = c("Classic RE",
                         "Hartung-Knapp",
                         'Hartung-Knapp, adhoc.hakn = "ci"',
                         'Hartung-Knapp, adhoc.hakn = "se"',
                         'Hartung-Knapp, adhoc.hakn = "iqwig6"'))
##
summary(mb0, overall = FALSE, overall.hetstat = FALSE)


## van Aert and Jackson (2019), RSM, equation (5)
##
Hstar <- with(m0,
          1 / (k - 1) * sum(w.random * (TE - TE.random)^2, na.rm = TRUE))
Hstar

## van Aert and Jackson (2019), RSM, p. 525
##
z <- qnorm(0.975)
t <- qt(0.975, df = m0$k - 1)
##
Hstar
Hstar < z / t

set.seed(1234)
m <- m0
m$TE <- rep(1, length(m$TE)) + rnorm(length(m$TE), 0, 0.01)


m1 <- metagen(TE, seTE, data = m, method.tau = "REML")
##
mb1 <- metabind(m1,
                update(m1, hakn = TRUE, adhoc.hakn = ""),
                update(m1, hakn = TRUE, adhoc.hakn = "ci"),
                update(m1, hakn = TRUE, adhoc.hakn = "se"),
                update(m1, hakn = TRUE, adhoc.hakn = "iqwig6"),
                name = c("Classic RE",
                         "Hartung-Knapp",
                         'Hartung-Knapp, adhoc.hakn = "ci"',
                         'Hartung-Knapp, adhoc.hakn = "se"',
                         'Hartung-Knapp, adhoc.hakn = "iqwig6"'))
##
summary(mb1, overall = FALSE, overall.hetstat = FALSE)


m2 <- update(m1, subset = 1:2)
##
mb2 <- metabind(m2,
                update(m2, hakn = TRUE, adhoc.hakn = ""),
                update(m2, hakn = TRUE, adhoc.hakn = "ci"),
                update(m2, hakn = TRUE, adhoc.hakn = "se"),
                update(m2, hakn = TRUE, adhoc.hakn = "iqwig6"),
                name = c("Classic RE",
                         "Hartung-Knapp",
                         'Hartung-Knapp, adhoc.hakn = "ci"',
                         'Hartung-Knapp, adhoc.hakn = "se"',
                         'Hartung-Knapp, adhoc.hakn = "iqwig6"'))
##
summary(mb2, overall = FALSE, overall.hetstat = FALSE)
