CTmeta | R Documentation |
Continuous-time meta-analysis (CTmeta) on standardized lagged effects (Phi) taking into account the various time-intervals used in the primary studies. There is also an interactive web application on my website to perform CTmeta: https://www.uu.nl/staff/RMKuiper/Websites%20%2F%20Shiny%20apps.
CTmeta(
N,
DeltaT,
DeltaTStar,
Phi,
SigmaVAR = NULL,
Gamma = NULL,
Moderators = 0,
Mod = NULL,
FEorRE = 1,
BetweenLevel = NULL,
Label = NULL,
alpha = 0.05,
PrintPlot = FALSE
)
N |
Number of persons (panel data) or number of measurement occasions - 1 (time series data) used in the S primary studies. Matrix of size S times 1. |
DeltaT |
The time intervals used in the S primary studies. Matrix of size S times 1. Note that all the time intervals should be on the same scale (e.g., two time-intervals of 60 minutes and 2 hours, should be either 60 and 120 or 1 and 2). |
DeltaTStar |
The time interval (scalar) to which the standardized lagged effects matrix should be transformed to. |
Phi |
Stacked matrix of size S*q times q or array with dimensions q times q times S of (un)standardized lagged effects matrices for all S primary studies in the meta-analysis; with q the number of variables (leading to an q times q lagged effects matrix in a single primary study). Note: In case primary studies report (lagged) correlation matrices, one can use the function 'TransPhi_Corr' to transform those to the corresponding standardized lagged effects matrices (see ?TransPhi_Corr and examples below). |
SigmaVAR |
Stacked matrix of size S*q times q or array with dimensions q times q times S of residual covariance matrices of the first-order discrete-time vector autoregressive (DT-VAR(1)) model. |
Gamma |
Optional (either SigmaVAR or Gamma). Stacked matrix of size S*q times q or array with dimensions q times q times S of stationary covariance matrices, that is, the contemporaneous covariance matrices of the data sets. Note that if Phi and Gamma are known, SigmaVAR can be calculated. Hence, only SigmaVAR or Gamma is needed as input (if only Gamma, then use 'Gamma = Gamma' or set SigmaVAR to NULL, see examples below). |
Moderators |
Optional. Indicator (TRUE/FALSE or 1/0) whether there are moderators to be included (TRUE or 1) or not (FALSE or 0). By default, Moderators = 0. |
Mod |
Optional. An S x m matrix of m moderators to be included in the analysis when Moderators == 1. By default, Mod = NULL. |
FEorRE |
Optional. Indicator (1/2) whether continuous-time meta-analysis should use a fixed-effects model (1) or random-effects model (2). By default, FEorRE = 1. |
BetweenLevel |
Optional. Needed in case of a 2-level multilevel meta-analysis. BetweenLevel should be an S-vector or S x 1 matrix (namely one value for each study). It will only be used, if FEorRE == 2 (i.e., in a random-effects model). Then, one can add a between-level to the random part (e.g., sample number if multiple studies use the sample such that there is dependency between those studies). Note that the within level is Study number. By default, BetweenLevel = NULL. |
Label |
Optional. If one creates, for example, a funnel or forest plot the labeling used in the rma.mv function is used. Label should be an q*q*S-vector (namely one value for each of the elements in a study-specific Phi (of size q x q) and for each study). It will only be used, in case the multivariate approach can be used (in case of the univariate approach, it will always use the labeling Study 1 to Study S). By default, Label = NULL; then it will use the labeling Study 1 Phi11, Study Phi 12, ..., Study 1 Phi qq, ... Study S Phi11, ..., Study S Phiqq. |
alpha |
Optional. The alpha level in determining the (1-alpha)*100% confidence interval (CI). By default, alpha = 0.05; resulting in a 95% CI. |
PrintPlot |
Optional. Indicator (TRUE/FALSE or 1/0) for rendering a Phi-plot (TRUE or 1) or not (FALSE or 0). By default, PrintPlot = FALSE. |
The output comprises, among others, the overall vectorized transformed standardized lagged effects, their covariance matrix, and the corresponding elliptical/multivariate 95% CI.
# library(CTmeta)
##################################################################################################
# Input needed in examples below with q=2 variables and S=3 primary studies
#
N <- matrix(c(643, 651, 473))
DeltaT <- matrix(c(2, 3, 1))
DeltaTStar <- 1
#
# I will use the example matrices stored in the package:
Phi <- myPhi
SigmaVAR <- mySigmaVAR
Gamma <- myGamma # Note: CTmeta does not need both SigmaVAR and Gamma, as denomstrated below.
# These are all three stacked matrices of size S*q times q.
# The CTmeta function will standardize these matrices (to make comparison of effects meaningful).
#
Moderators = 0 # By default set to 0. Hence, not per se needed, as demonstrated below.
##################################################################################################
# Below, you can find example code. Note that, here, only 3 primary studies are used.
# In practice, one would normally have (many) more, but the code stays (more or less) the same.
### Examples without comments ###
## Example without moderators ##
# Fixed effects model #
# Run CTmeta
CTma <- CTmeta(N, DeltaT, DeltaTStar, Phi, SigmaVAR)
CTma
summary(CTma)
# Random effects (RE) model #
CTmeta(N, DeltaT, DeltaTStar, Phi, SigmaVAR, FEorRE = 2)
BetweenLevel <- c(1, 1, 2) # Assuming the first two studies used the same sample/dataset
CTmeta(N, DeltaT, DeltaTStar, Phi, SigmaVAR, FEorRE = 2, BetweenLevel = BetweenLevel) # Two-level RE meta-analysis example
## Example with moderators ##
Mod <- matrix(c(64,65,47)) # 1 moderator
CTma <- CTmeta(N, DeltaT, DeltaTStar, Phi, SigmaVAR, Moderators = 1, Mod = Mod) # fixed effects model
summary(CTma)
q <- dim(Phi)[2]; Mod <- matrix(cbind(c(64,65,47), c(78,89,34)), ncol = q); colnames(Mod) <- c("Mod1", "Mod2") # two moderators, in each column 1
CTma <- CTmeta(N, DeltaT, DeltaTStar, Phi, SigmaVAR, Moderators = 1, Mod = Mod, FEorRE = 2); CTma$tau2 # random effects model
summary(CTma)
Mod <- matrix(c(64,65,47)) # 1 moderator
BetweenLevel <- c(1, 1, 2) # Assuming the first two studies used the same sample/dataset.
CTmeta(N, DeltaT, DeltaTStar, Phi, SigmaVAR, Moderators = 1, Mod = Mod, FEorRE = 2, BetweenLevel = BetweenLevel) # Two-level RE meta-analysis example
## funnel and forest plots ##
CTma <- CTmeta(N, DeltaT, DeltaTStar, Phi, SigmaVAR, FEorRE = 2)
funnel(CTma$summaryMetaAnalysis, label = 'out')
forest(CTma$summaryMetaAnalysis)
# One can do the same for the two-level analysis.
#
# Note, in case a univariate approach had to be taken, leading to multiple analyses, then one should create a plot per analysis:
# lapply(CTma$summaryMetaAnalysis_jk, funnel, label = 'out')
#
# In case you want to create a plot per element of the study-specific Phi's:
#elt <- rep(F, (q*q))
#elt[1] <- T # First element out of q*q true, so referring to element Phi11.
#yi_Phi11 <- CTma$summaryMetaAnalysis$yi[elt]
#vi_Phi11 <- CTma$summaryMetaAnalysis$vi[elt]
#funnel(yi_Phi11, vi_Phi11, label = 'out')
#forest(yi_Phi11, vi_Phi11)
elt <- rep(F, (q*q))
elt[2] <- T # Second element out of q*q true, so referring to element Phi12.
yi_Phi12 <- CTma$summaryMetaAnalysis$yi[elt]
vi_Phi12 <- CTma$summaryMetaAnalysis$vi[elt]
funnel(yi_Phi12, vi_Phi12, label = 'out')
forest(yi_Phi12, vi_Phi12)
## Make customized Phi-plot of resulting overall Phi ##
CTma <- CTmeta(N, DeltaT, DeltaTStar, Phi, SigmaVAR, PrintPlot = TRUE)
CTma$PhiPlot
out_CTmeta <- CTmeta(N, DeltaT, DeltaTStar, Phi, SigmaVAR)
overallPhi <- out_CTmeta$Overall_standPhi
Title <- as.list(expression(paste0(Phi(Delta[t]), " plot:"),
"How do the overall lagged parameters vary as a function of the time-interval"))
PhiPlot(DeltaTStar, overallPhi, Min = 0, Max = 40, Step = 0.5, Title = Title)
# or
ggPhiPlot <- ggPhiPlot(DeltaTStar, overallPhi, Min = 0, Max = 40, Step = 0.5, Title = Title)
ggPhiPlot$PhiPlot
## Evaluate dominance of (absolute values of) cross-lagged effects ##
if (!require("restriktor")) install.packages("restriktor") # Use restriktor package for function goric().
# Authors of goric(): Vanbrabant and Kuiper.
library(restriktor)
out_CTmeta <- CTmeta(N, DeltaT, DeltaTStar, Phi, SigmaVAR)
H1 <- "abs(overallPhi12) < abs(overallPhi21)"
goric(out_CTmeta, hypotheses = list(H1), type = "gorica", comparison = "complement")
out_CTmeta <- CTmeta(N, DeltaT, DeltaTStar, Phi, SigmaVAR)
H1 <- "abs(overallPhi12) < abs(overallPhi21); abs(overallPhi11) < abs(overallPhi22)"
goric(out_CTmeta, hypotheses = list(H1), type = "gorica", comparison = "complement")
out_CTmeta <- CTmeta(N, DeltaT, DeltaTStar, Phi, SigmaVAR)
est <- coef(out_CTmeta) # or: est <- out_CTmeta$Overall_vecStandPhi_DeltaTStar
VCOV <- vcov(out_CTmeta) # or: VCOV <- out_CTmeta$CovMx_OverallPhi_DeltaTStar
#goric(est, VCOV = VCOV, hypotheses = list(H1), type = "gorica", comparison = "complement")
# With this input, one can only obtain the GORICA, so one can also use:
goric(est, VCOV = VCOV, hypotheses = list(H1), comparison = "complement")
## What if primary studies report a (lagged) correlation matrix ##
q <- 2
corr_YXYX <- matrix(c(1.00, 0.40, 0.63, 0.34,
0.40, 1.00, 0.31, 0.63,
0.63, 0.31, 1.00, 0.41,
0.34, 0.63, 0.41, 1.00), byrow = T, ncol = 2*q)
N <- matrix(c(643, 651, 473))
DeltaT <- matrix(c(2, 3, 1))
DeltaTStar <- 1
#
# Create input:
out_1 <- TransPhi_Corr(DeltaTStar = DeltaT[1], DeltaT = 1, N = N[1], corr_YXYX)
Phi_1 <- out_1$standPhi_DeltaTStar
SigmaVAR_1 <- out_1$standSigmaVAR_DeltaTStar
out_2 <- TransPhi_Corr(DeltaTStar = DeltaT[2], DeltaT = 1, N = N[2], corr_YXYX)
Phi_2 <- out_2$standPhi_DeltaTStar
SigmaVAR_2 <- out_2$standSigmaVAR_DeltaTStar
out_3 <- TransPhi_Corr(DeltaTStar = DeltaT[3], DeltaT = 1, N = N[3], corr_YXYX)
Phi_3 <- out_3$standPhi_DeltaTStar
SigmaVAR_3 <- out_3$standSigmaVAR_DeltaTStar
#
Phi <- rbind(Phi_1, Phi_2, Phi_3) # This, returns a stacked matrix of size S q times q.
SigmaVAR <- rbind(SigmaVAR_1, SigmaVAR_2, SigmaVAR_3)
out_CTmeta <- CTmeta(N, DeltaT, DeltaTStar, Phi, SigmaVAR)
out_CTmeta$Overall_standPhi
##############################
### Examples with comments ###
## Example without moderators ##
# Fixed effects model #
# Run CTmeta with, for instance,
CTmeta(N, DeltaT, DeltaTStar, Phi, SigmaVAR)
# There are multiple options; use one of the following:
CTmeta(N, DeltaT, DeltaTStar, Phi, SigmaVAR, Gamma, Moderators, Mod, 1) # The 1, here, says FEorRE = 1
CTmeta(N, DeltaT, DeltaTStar, Phi, SigmaVAR, Gamma, Moderators, Mod) # Notably, if Moderators = 0, 'Mod' is not inspected.
CTmeta(N, DeltaT, DeltaTStar, Phi, SigmaVAR, Gamma, Moderators)
CTmeta(N, DeltaT, DeltaTStar, Phi, SigmaVAR, Gamma)
CTmeta(N, DeltaT, DeltaTStar, Phi, SigmaVAR) # Says, implicitly, Gamma = NULL
CTmeta(N, DeltaT, DeltaTStar, Phi, Gamma = Gamma) # Says, implicitly, SigmaVAR = NULL
CTmeta(N, DeltaT, DeltaTStar, Phi, NULL, Gamma) # Says SigmaVAR = NULL
# Note: Do NOT use
#CTmeta(N, DeltaT, DeltaTStar, Phi, Gamma)
# Then, CTmeta incorrectly uses SigmaVAR = Gamma.
# Different types of output options are possible:
CTma <- CTmeta(N, DeltaT, DeltaTStar, Phi, SigmaVAR)
CTma # gives same as print(CTma)
summary(CTma)
print(CTma, digits = 4)
summary(CTma, digits = 4)
# In Rstudio, use 'CTma$' to see what output there is. For example:
CTma$summaryMetaAnalysis
# Random effects model #
# Add "FEorRE = 2"; e.g.,
CTmeta(N, DeltaT, DeltaTStar, Phi, SigmaVAR, FEorRE = 2)
# Two-level RE meta-analysis example
BetweenLevel <- c(1, 1, 2) # Assuming the first two studies used the same sample/dataset
CTmeta(N, DeltaT, DeltaTStar, Phi, SigmaVAR, FEorRE = 2, BetweenLevel = BetweenLevel) # Two-level RE meta-analysis example
# Note, one can also use this in case there are moderators (as in the example below).
## Example with moderators ##
Mod <- matrix(c(64,65,47)) # 1 moderator
#q <- dim(Phi)[2]; Mod <- matrix(cbind(c(64,65,47), c(78,89,34)), ncol = q); colnames(Mod) <- c("Mod1", "Mod2") # two moderators, in each column 1
CTma <- CTmeta(N, DeltaT, DeltaTStar, Phi, SigmaVAR, Moderators = 1, Mod = Mod) # fixed effects model
#CTma <- CTmeta(N, DeltaT, DeltaTStar, Phi, SigmaVAR, Moderators = 1, Mod = Mod, FEorRE = 2); CTma$tau2 # random effects model
summary(CTma)
# Two-level RE meta-analysis example
Mod <- matrix(c(64,65,47)) # 1 moderator
BetweenLevel <- c(1, 1, 2) # Assuming the first two studies used the same sample/dataset
CTmeta(N, DeltaT, DeltaTStar, Phi, SigmaVAR, Moderators = 1, Mod = Mod, FEorRE = 2, BetweenLevel = BetweenLevel) # Two-level RE meta-analysis example
## funnel and forest plots ##
CTma <- CTmeta(N, DeltaT, DeltaTStar, Phi, SigmaVAR, FEorRE = 2)
funnel(CTma$summaryMetaAnalysis, label = 'out')
forest(CTma$summaryMetaAnalysis)
# One can do the same for the two-level analysis.
#
# Note, in case a univariate approach had to be taken, leading to multiple analyses, then one should create a plot per analysis:
# lapply(CTma$summaryMetaAnalysis_jk, funnel, label = 'out')
#
# Notes on funnel and forest:
# - These plots are now based on the q*q elements in the study-specific Phi's and the S studies.
# See below, how you can create these plots per element of Phi. This should then be done separately for all q*q elements.
# - In case label names are too long or not insightful, one can change them by using the Label argument.
# In case of the funnel plot, one can also decide to not use the "label = 'out'" part.
#
# Notes on forest:
# - For random-effects models of class "rma.mv" (see rma.mv) with multiple values, the addpred argument can be used to specify for which level of the inner factor the prediction interval should be provided (since the intervals differ depending on the value).
# - One can also look into the functionality of addpoly().
#
# In case you want to create a plot per element of the study-specific Phi's:
#elt <- rep(F, (q*q))
#elt[1] <- T # First element out of q*q true, so referring to element Phi11.
#yi_Phi11 <- CTma$summaryMetaAnalysis$yi[elt]
#vi_Phi11 <- CTma$summaryMetaAnalysis$vi[elt]
#funnel(yi_Phi11, vi_Phi11, label = 'out')
#forest(yi_Phi11, vi_Phi11)
elt <- rep(F, (q*q))
elt[2] <- T # Second element out of q*q true, so referring to element Phi12.
yi_Phi12 <- CTma$summaryMetaAnalysis$yi[elt]
vi_Phi12 <- CTma$summaryMetaAnalysis$vi[elt]
funnel(yi_Phi12, vi_Phi12, label = 'out')
forest(yi_Phi12, vi_Phi12)
## Make customized Phi-plot of resulting overall Phi ##
# Option 1: Using the plot option in the function:
CTmeta(N, DeltaT, DeltaTStar, Phi, SigmaVAR, PrintPlot = TRUE)
# The plot can be stored and retrieved as an object as follows:
# CTma <- CTmeta(N, DeltaT, DeltaTStar, Phi, SigmaVAR, PrintPlot = TRUE)
# CTma$PhiPlot
# Option 2: A customized Phi-plot can be made using the function 'PhiPlot' (see below) or by using the interactive web app from my website (\url{https://www.uu.nl/staff/RMKuiper/Websites\%20\%2F\%20Shiny\%20apps}).
# Alternatively, one can use the function 'ggPhiPlot' instead of 'PhiPlot'.
# Extract the (q times q) overall Phi matrix
out_CTmeta <- CTmeta(N, DeltaT, DeltaTStar, Phi, SigmaVAR)
# resulting overall Phi:
overallPhi <- out_CTmeta$Overall_standPhi
# Make Phi-plot:
Title <- as.list(expression(paste0(Phi(Delta[t]), " plot:"),
"How do the overall lagged parameters vary as a function of the time-interval"))
PhiPlot(DeltaTStar, overallPhi, Min = 0, Max = 40, Step = 0.5, Title = Title)
# or
ggPhiPlot(DeltaTStar, overallPhi, Min = 0, Max = 40, Step = 0.5, Title = Title)
# The plot can be stored and retrieved as an object as follows:
# ggPhiPlot <- ggPhiPlot(DeltaTStar, overallPhi, Min = 0, Max = 40, Step = 0.5, Title = Title)
# ggPhiPlot$PhiPlot
## Evaluate dominance of (absolute values of) cross-lagged effects ##
if (!require("restriktor")) install.packages("restriktor") # Use restriktor package for function goric().
# Authors of goric(): Vanbrabant and Kuiper.
#
# If version from github needed:
#if (!require("devtools")) install.packages("devtools")
#library(devtools)
#install_github("LeonardV/restriktor")
#
library(restriktor)
# Option 1
# Use CTmeta object
out_CTmeta <- CTmeta(N, DeltaT, DeltaTStar, Phi, SigmaVAR)
#
# Example 1: dominance of (absolute values of) cross-lagged effects
# Specify hypothesis
H1 <- "abs(overallPhi12) < abs(overallPhi21)"
#H2 <- "abs(overallPhi12) > abs(overallPhi21)" # = complement of H1 and does not need to be specified, see below.
# Evaluate dominance of cross-lagged effects via AIC-type criterion called the GORICA (Altinisik, Nederhof, Hoijtink, Oldehinkel, Kuiper, accepted 2021).
#goric(out_CTmeta, hypotheses = list(H1, H2), type = "gorica", comparison = "none")
# or, since H2 is complement of H1, equivalently:
goric(out_CTmeta, hypotheses = list(H1), type = "gorica", comparison = "complement")
#
# Example 2: dominance of (absolute values of) cross-lagged effects and autoregressive effects
H1 <- "abs(overallPhi12) < abs(overallPhi21); abs(overallPhi11) < abs(overallPhi22)"
# Note that, now, specification of complement 'by hand' harder.
goric(out_CTmeta, hypotheses = list(H1), type = "gorica", comparison = "complement")
#
#
# Option 2
# Extract the vectorized overall standardized overallPhi matrix and its covariance matrix
# using the functions coef() and vcov()
out_CTmeta <- CTmeta(N, DeltaT, DeltaTStar, Phi, SigmaVAR)
est <- coef(out_CTmeta) # or: est <- out_CTmeta$Overall_vecStandPhi_DeltaTStar
VCOV <- vcov(out_CTmeta) # or: VCOV <- out_CTmeta$CovMx_OverallPhi_DeltaTStar
goric(est, VCOV = VCOV, hypotheses = list(H1), type = "gorica", comparison = "complement")
## What if primary studies report a (lagged) correlation matrix ##
# Suppose, for ease, that all S=3 primary studies reported the following lagged correlation matrix:
q <- 2
corr_YXYX <- matrix(c(1.00, 0.40, 0.63, 0.34,
0.40, 1.00, 0.31, 0.63,
0.63, 0.31, 1.00, 0.41,
0.34, 0.63, 0.41, 1.00), byrow = T, ncol = 2*q)
# In the example below, the same N and DeltaT(Star) values are used:
N <- matrix(c(643, 651, 473))
DeltaT <- matrix(c(2, 3, 1))
DeltaTStar <- 1
# Use the function 'TransPhi_Corr' to calculate the corresponding standardized lagged effects matrix per primary study.
# Note that one can already make the time-intervals equal via the arguments DeltaTStar and DeltaT, but CTmeta can as well.
# In this example, I deliberately make the time-intervals unequal, such that the example is in line with the input (i.e., DeltaT <- matrix(c(2, 3, 1))) and such the resulting overall Phi should equal the Phi that underlies this lagged correlation matrix (which I check at the end).
out_1 <- TransPhi_Corr(DeltaTStar = DeltaT[1], DeltaT = 1, N = N[1], corr_YXYX)
Phi_1 <- out_1$standPhi_DeltaTStar
SigmaVAR_1 <- out_1$standSigmaVAR_DeltaTStar
out_2 <- TransPhi_Corr(DeltaTStar = DeltaT[2], DeltaT = 1, N = N[2], corr_YXYX)
Phi_2 <- out_2$standPhi_DeltaTStar
SigmaVAR_2 <- out_2$standSigmaVAR_DeltaTStar
out_3 <- TransPhi_Corr(DeltaTStar = DeltaT[3], DeltaT = 1, N = N[3], corr_YXYX)
Phi_3 <- out_3$standPhi_DeltaTStar
SigmaVAR_3 <- out_3$standSigmaVAR_DeltaTStar
# Make Phi
Phi <- rbind(Phi_1, Phi_2, Phi_3) # This, returns a stacked matrix of size S q times q.
SigmaVAR <- rbind(SigmaVAR_1, SigmaVAR_2, SigmaVAR_3)
# For more details, see ?TransPhi_Corr
# The example CTmeta() code above can be run using this Phi; e.g.,
CTmeta(N, DeltaT, DeltaTStar, Phi, SigmaVAR)
# The overall q-times-q (here, 2x2) lagged effects matrix Phi
out_CTmeta <- CTmeta(N, DeltaT, DeltaTStar, Phi, SigmaVAR)
out_CTmeta$Overall_standPhi
#
# As a check, to see indeed that CTmeta works properly (where the resulting Phi is independent of the choise of N).
TransPhi_Corr(DeltaTStar = 1, DeltaT = 1, N = 100, corr_YXYX)$standPhi_DeltaTStar
# Note that is normally not a check you would do.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.