## By following these examples, you can appropriately specify hypotheses based on
## your research questions and analytical framework.
# The hypotheses (i.e., constraints) have to be in a list. It is recommended to name
# each hypothesis in the list. Otherwise, the hypotheses are named accordingly 'H1', 'H2', \ldots.
# Example using text-based syntax (the labels x1, x2, and x3 are the names of
# coef(model) or names(vector))
H1 <- '(x1, x2, x3) > 0'
H2 <- '(x1, x3) > 0; x2 < 0'
H3 <- 'x1 > 0; x2 < 0; x3 = 0'
hypotheses = list(hypo1 = H1, hypo2 = H2, hypo3 = H3)
# Here, it is assumed that each study uses the same hypothesis specification
# Example using mixed syntax:
H1 <- 'x1 > x2 > x3 > 0'
hypotheses = list(Ha = H1, Hb = 'x1 = x2 > x3')
# Here, it is assumed that each study uses the same hypothesis specification
# Example using a different set of hypotheses for each study:
# Set for Study 1
H11 <- 'group1 = group2 > group3'
H12 <- 'group2 > group1 > group3'
# Set for Study 2
H21 <- 'gr1 = gr2 > gr3'
H22 <- 'gr2 > gr1 > gr3'
#
# correct
hypotheses = list(set1 = list(H11, H12), set2 = list(H21, H22))
# note that the list names set1 and set2 are redundant and can be left out.
# It is crucial to ensure that the hypotheses across each set are ordered in a similar manner.
#
# NOT correct
#hypotheses = list(set1 = list(H12, H11), set2 = list(H21, H22))
## Example 1: based on one estimate and its variance
# 4 studies
# estimates (beta) (from the 4 primary studies)
est_1 <- c(beta1 = 0.09)
est_2 <- c(beta1 = 0.14)
est_3 <- c(beta1 = 1.09)
est_4 <- c(beta1 = 1.781)
Param_studies <- list(est_1, est_2, est_3, est_4)
# standard error of the beta's (from the 4 primary studies)
vcov_est_1 <- matrix(c(0.029^2), nrow = 1)
vcov_est_2 <- matrix(c(0.054^2), nrow = 1)
vcov_est_3 <- matrix(c(0.093^2), nrow = 1)
vcov_est_4 <- matrix(c(0.179^2), nrow = 1)
CovMx_studies <- list(vcov_est_1, vcov_est_2, vcov_est_3, vcov_est_4)
# Set of hypotheses for each study
# Note: in this case the same for each study
H0 <- "beta1 = 0"
Hpos <- "beta1 > 0"
Hneg <- "beta1 < 0"
hypotheses <- list(H0 = H0, Hpos = Hpos, Hneg = Hneg)
# Since this covers the whole space / covers all theories, we do not need a safeguard-hypothesis:
comparison <- "none"
evS4_added <- evSyn(object = Param_studies, VCOV = CovMx_studies,
hypotheses = hypotheses,
type_ev = "added",
comparison = "none")
evS4_added
summary(evS4_added)
plot(evS4_added)
evS4_equal <- evSyn(object = Param_studies, VCOV = CovMx_studies,
hypotheses = hypotheses,
type_ev = "equal",
comparison = "none")
evS4_equal
summary(evS4_equal)
plot(evS4_equal)
## Example 2: based on (multiple) estimates and their covariance matrices
# 2 studies
# estimates and their covariance matrices (from the 2 primary studies)
est_1 <- c(1.88, 2.54, 0.02)
names(est_1) <- c("group1", "group2", "group3")
vcov_est_1 <- diag(c(0.2149074, 0.2149074, 0.1408014))
est_2 <- c(0.98, 0.02, 0.27)
names(est_2) <- c("gr1", "gr2", "gr3")
vcov_est_2 <- diag(c(0.1382856, 0.1024337, 0.0987754))
# list of estimates
object <- list(est_1, est_2)
# list of covariance matrices of the estimates
VCOV <- CovMx_studies <- list(vcov_est_1, vcov_est_2)
# Hypotheses (study-specific)
# names(est_1) # Specify restrictions using those names
H11 <- 'group1 = group2 > group3'
H12 <- 'group2 > group1 > group3'
# names(est_2) # Specify restrictions using those names
H21 <- 'gr1 = gr2 > gr3'
H22 <- 'gr2 > gr1 > gr3'
# list of (competing) hypotheses
hypotheses <- list(H1 = list(H11, H12), H2 = list(H21, H22))
# Evidence synthesis (added-evidence)
evS2_added <- evSyn(object, VCOV = VCOV, hypotheses = hypotheses,
type_ev = "added", comparison = "unconstrained")
evS2_added
plot(evS2_added)
## Example 3: based on (multiple) estimates and their covariance matrices
# 3 studies
# generate data
ratio <- c(1,1.1,1.2)
n <- c(30, 50, 100)
# Generate data1
n1 <- n[1]
x11 <- rnorm(n1)
x12 <- rnorm(n1)
x13 <- rnorm(n1)
data <- cbind(x11, x12, x13)
# Standardize data - since parameters for continuous variables will be compared
data1 <- as.data.frame(scale(data))
y1 <- ratio[1]*data1$x11 + ratio[2]*data1$x12 + ratio[3]*data1$x13 + rnorm(n1)
# Note: since there is one outcome, the outcome does not need to be standardized.
fit.lm1 <- lm(y1 ~ 1 + x11 + x12 + x13, data = data1)
# Generate data2
n2 <- n[2]
x21 <- rnorm(n2)
x22 <- rnorm(n2)
x23 <- rnorm(n2)
data <- cbind(x21, x22, x23)
data2 <- as.data.frame(scale(data))
y2 <- ratio[1]*data2$x21 + ratio[2]*data2$x22 + ratio[3]*data2$x23 + rnorm(n2)
fit.lm2 <- lm(y2 ~ 1 + x21 + x22 + x23, data = data2)
# Generate data3
n3 <- n[3]
x31 <- rnorm(n3)
x32 <- rnorm(n3)
x33 <- rnorm(n3)
data <- cbind(x31, x32, x33)
data3 <- as.data.frame(scale(data))
y3 <- ratio[1]*data3$x31 + ratio[2]*data3$x32 + ratio[3]*data3$x33 + rnorm(n3)
fit.lm3 <- lm(y3 ~ 1 + x31 + x32 + x33, data = data3)
# Extract estimates and their covariance matrix (per study)
est_1 <- coef(fit.lm1)
est_2 <- coef(fit.lm2)
est_3 <- coef(fit.lm3)
vcov_est_1 <- vcov(fit.lm1)
vcov_est_2 <- vcov(fit.lm2)
vcov_est_3 <- vcov(fit.lm3)
names(est_1) <- names(est_2) <- names(est_3) <- c("intercept", "x1", "x2", "x3")
# Parameter estimate values from the primary studies
Param_studies <- list(est_1, est_2, est_3)
# standard error of the beta's
CovMx_studies <- list(vcov_est_1, vcov_est_2, vcov_est_3)
# Set of hypotheses for each study. Note: in this case the same for each study
hypothesis <- 'x1 < x2 < x3'
# In our theory, we compare estimates of continuous variables, so we standardized
# the data beforehand to ensure comparability. In 'Param_studies' and 'CovMx_studies',
# the intercept can be omitted without affecting the GORIC(A) weights, as there are
# no restrictions on it. Since we have only one theory-based hypothesis, we will
# utilize the more powerful complement of the hypothesis (Vanbrabant, Van Loey, Kuiper, 2019).
# The complement represents the remaining 11 theories, while the unconstrained
# scenario includes all 12 possible theories, including H1.
# Evidence synthesis (added-evidence)
evS3 <- evSyn(object = Param_studies, VCOV = CovMx_studies,
hypotheses = list(H1 = hypothesis),
type_ev = "added",
comparison = "complement")
evS3
plot(evS3)
## Example 4: based on loglikelihood values and penalty values
# Make LL and PT values a list
LL <- as.list(data.frame(t(myLLs)))
penalty.values <- as.list(data.frame(t(myPTs)))
# Evidence synthesis
#
# Added-evidence
evS_LL_added <- evSyn(object = LL, PT = penalty.values,
hypo_names = colnames(myLLs),
type_ev = "added")
evS_LL_added
#
# Equal-evidence
evS_LL_equal <- evSyn(object = LL, PT = penalty.values,
hypo_names = colnames(myLLs),
type_ev = "equal")
evS_LL_equal
## Example 5: based on AIC, ORIC, GORIC(A) values
# Make GORIC values a list
goric.values <- as.list(data.frame(t(myGORICs)))
# Evidence synthesis (added-evidence)
evS_Gv <- evSyn(goric.values,
hypo_names = colnames(myGORICs))
evS_Gv
## Example 6: based on an escalc object
## Example 6a: one outcome, default labeling
if (requireNamespace("metadat", quietly = TRUE)) {
# access data within metafor
data_a <- metadat::dat.berkey1998[metadat::dat.berkey1998$outcome == "PD", ][, c(1:4, 6:7)]
# yi gives the estimates and vi their variances.
# So, here:
# cluster_col = "trial" using the first variable with a cluster variable
# name (year and author are other possibilities)
# outcome_col = NULL # No 'outcome' variable, so 'theta' is used as labeling
# yi_col = "yi"
# vi_cols = "vi"
#
# The hypothesis should thus now be based on the default labelling theta:
H1 <- "abs(theta) > 0.5"
evS_escalc_a <- evSyn(
data_a,
hypotheses = list(H1 = H1),
type_ev = "equal",
comparison = "complement"
)
evS_escalc_a
}
## Example 6b: one outcome, default labeling coming from 'outcome' variable
if (requireNamespace("metadat", quietly = TRUE)) {
data_b <- metadat::dat.berkey1998[metadat::dat.berkey1998$outcome == "PD", ][, c(1:7)]
# yi gives the estimates and vi their variances.
# So, here:
# cluster_col = "trial" # using the first variable with a cluster variable
# name (year and author are other possibilities)
# outcome_col = "outcome" # was NULL, but found 'outcome' variable, so use that
# labeling, that is, "PD"
# yi_col = "yi"
# vi_cols = "vi"
#
# The hypothesis should thus now be based on the labeling based on the 'outcome' variable:
H1 <- "abs(PD) > 0.5"
evS_escalc_b <- evSyn(data_b,
hypotheses = list(H1 = H1),
type_ev = "equal",
comparison = "complement")
evS_escalc_b
}
## Example 6c: 2 outcomes, default labeling coming from 'outcome' variable
if (requireNamespace("metadat", quietly = TRUE)) {
data_c <- metadat::dat.berkey1998
# yi gives the estimates (for both outcomes), vi their variances,
# and v1i & v2i combined their covariance matrices.
# So, here:
# cluster_col = "trial" # using the first variable with a cluster variable
# name (year and author are other possibilities)
# outcome_col = "outcome" # was NULL, but found 'outcome' variable, so use
# that labeling, that is, "PD"
# yi_col = "yi"
# vi_cols = c("v1i", "v2i")
#
# The hypothesis should thus now be based on the labeling based on the 'outcome' variable:
H1 <- "abs(AL) > 0.5, abs(PD) > 0.5"
evS_escalc_c <- evSyn(data_c,
hypotheses = list(H1 = H1),
type_ev = "equal",
comparison = "complement")
evS_escalc_c
}
Run the code above in your browser using DataLab