The function r.vcov computes effect sizes and variance-covariance matrices between effect sizes of the same of different types. Effect sizes include mean difference (MD), standardized mean difference (SMD), log odds ratio (logOR), log risk ratio (logRR), and risk difference (RD). Formulas are in Table I of Wei et al.'s paper (2013).

mix.vcov(d, r, nt, nc, st, sc, n_rt = NA, n_rc = NA, sdt, sdc, type, name = NULL, na.impute = NA)

Arguments

d

A \(N \times p\) matrix or data frame with mean differences (MD) and/or standard mean differences (SMD) from the \(N\) studies. d[i,j] is the value from study \(i\) for outcome \(j\). If outcome \(j\) is not MD or SMD, NA has to be imputed in column \(j\).

r

A \(N\)-dimensional list of \(p \times p\) correlation matrices for the \(p\) outcomes from the \(N\) studies. r[[k]][i,j] is the correlation coefficient between outcome \(i\) and outcome \(j\) from study \(k\).

nt

A \(N \times p\) matrix storing sample sizes in the treatment group reporting the \(p\) outcomes. nt[i,j] is the sample size from study \(i\) reporting outcome \(j\).

nc

A matrix defined in a similar way as nt for the control group.

st

A \(N \times p\) matrix recording number of participants with event for all outcomes (dichotomous) in treatment group. st[i,j] reports number of participants with event for outcome \(j\) in treatment group for study \(i\). If outcome \(j\) is not dichotomous, NA has to be imputed in column \(j\).

sc

Defined in a similar way as st for the control group.

n_rt

A \(N\)-dimensional list of \(p \times p\) matrices storing sample sizes in the treatment group reporting pairwise outcomes in the off-diagonal elements. n_rt[[k]][i,j] is the sample size reporting both outcome \(i\) and outcome \(j\) from study \(k\). Diagonal elements of these matrices are discarded. The default value is NA, which means that the smaller sample size reporting the corresponding two outcomes is imputed: i.e. n_rt[[k]][i,j]=min(nt[k,i],nt[k,j]).

n_rc

A list defined in a similar way as n_rt for the control group.

sdt

A \(N \times p\) matrix storing sample standard deviations for each outcome from treatment group. sdt[i,j] is the sample standard deviation from study \(i\) for outcome \(j\). If outcome \(j\) is not continuous such as MD or SMD, NA has to be imputed in the \(j\)th column.

sdc

A matrix defined in a similar way as sdt for the control group.

type

A \(p\)-dimensional vector indicating types of effect sizes. "MD" stands for mean difference, "SMD" stands for standardized mean difference, "logOR" stands for log odds ratio, "logRR" stands for log risk ratio, and "RD" stands for risk difference.

name

A \(p\)-dimensional vector storing names for the effect sizes.

na.impute

Missing values in d can be imputed by a numeric value, such as zero by setting na.impute = 0. With the default setting na.impute = NA, missing values are not imputed. If specifying na.impute = "average", missing values are imputed by the mean of d that is sample-size weighted from the complete records.

Author

Min Lu

Value

ef

A \(N \times p\) data frame whose columns are computed effect sizes according to the input argument "type". SMD will be converted to Hedges's g (Wei and Higgins, 2013).

list.vcov

A \(N\)-dimensional list of \(p(p+1)/2 \times p(p+1)/2\) matrices of computed variance-covariance matrices.

matrix.vcov

A \(N \times p(p+1)/2\) matrix whose rows are computed variance-covariance vectors.

References

Wei, Y., & Higgins, J. (2013). Estimating within study covariances in multivariate meta-analysis with multiple outcomes. Statistics in Medicine, 32(7), 119-1205.

Examples

######################################################################### # Example: Geeganage2010 data # Preparing covariances for a multivariate meta-analysis # Choose variable SBP, DBP, DD, D with effect sizes "MD","MD","RD","lgOR" ######################################################################### data(Geeganage2010) ## set the correlation coefficients list r r12 <- 0.71 r13 <- 0.5 r14 <- 0.25 r23 <- 0.6 r24 <- 0.16 r34 <- 0.16 r <- vecTosm(c(r12, r13, r14, r23, r24, r34)) diag(r) <- 1 mix.r <- lapply(1:nrow(Geeganage2010), function(i){r}) ## compute variance co-variances computvcov <- mix.vcov(type = c("MD", "MD", "RD", "lgOR"), d = cbind(Geeganage2010$MD_SBP, Geeganage2010$MD_DBP, NA, NA), sdt = cbind(Geeganage2010$sdt_SBP, Geeganage2010$sdt_DBP, NA, NA), sdc = cbind(Geeganage2010$sdc_SBP, Geeganage2010$sdc_DBP, NA, NA), nt = as.matrix(subset(Geeganage2010, select = c(nt_SBP, nt_DBP, nt_DD, nt_D))), nc = as.matrix(subset(Geeganage2010, select = c(nc_SBP, nc_DBP, nc_DD, nc_D))), st = cbind(NA, NA, Geeganage2010$st_DD, Geeganage2010$st_D), sc = cbind(NA, NA, Geeganage2010$sc_DD, Geeganage2010$sc_D), r = mix.r, name = c("MD.SBP", "MD.DBP", "RD.DD", "lgOR.D")) # save different effect sizes in y y <- computvcov$ef head(y)
#> MD.SBP MD.DBP RD.DD lgOR.D #> 1 -2.47 -3.44 0.00000000 -1.0986123 #> 2 1.61 -0.34 0.18750000 0.5959834 #> 3 -8.16 -6.44 0.02554455 0.5892102 #> 4 -3.17 -3.41 0.04000000 0.4444945 #> 5 -0.15 -2.39 0.01920750 0.1000835 #> 6 -9.83 1.93 -0.25000000 -0.5108256
# save variances and covariances of all the effect sizes in a matrix S S <- computvcov$matrix.vcov S[1, ]
#> var_MD.SBP cov_MD.SBP_MD.DBP cov_MD.SBP_RD.DD cov_MD.SBP_lgOR.D var_MD.DBP #> 1 87.98831 34.81409 0.9245278 2.278204 27.85141 #> cov_MD.DBP_RD.DD cov_MD.DBP_lgOR.D var_RD.DD cov_RD.DD_lgOR.D var_lgOR.D #> 1 0.6207 0.7907191 0.040625 0.02741618 1.020833
## fixed-effect model MMA_FE <- summary(metafixed(y = y, Slist = computvcov$list.vcov)) ####################################################################### # Running random-effects model using package "mvmeta" or "metaSEM" ####################################################################### # Restricted maximum likelihood (REML) estimator from the mvmeta package library(mvmeta) mvmeta_RE <- summary(mvmeta(cbind(MD.SBP, MD.DBP, RD.DD, lgOR.D) ~., S = S, data = y, method = "reml")) mvmeta_RE
#> Call: mvmeta(formula = cbind(MD.SBP, MD.DBP, RD.DD, lgOR.D) ~ ., S = S, #> data = y, method = "reml") #> #> Multivariate random-effects meta-analysis #> Dimension: 4 #> Estimation method: REML #> #> Fixed-effects coefficients #> Estimate Std. Error z Pr(>|z|) 95%ci.lb 95%ci.ub #> MD.SBP -2.2985 1.6172 -1.4212 0.1552 -5.4682 0.8712 #> MD.DBP -2.5155 1.1940 -2.1067 0.0351 -4.8558 -0.1752 * #> RD.DD 0.0135 0.0195 0.6928 0.4885 -0.0247 0.0517 #> lgOR.D 0.0836 0.1004 0.8327 0.4050 -0.1132 0.2804 #> --- #> 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 #> MD.SBP 4.6739 MD.SBP MD.DBP RD.DD #> MD.DBP 3.9292 0.9776 #> RD.DD 0.0218 0.3666 0.1627 #> lgOR.D 0.0454 0.7369 0.5782 0.8991 #> #> Multivariate Cochran Q-test for heterogeneity: #> Q = 137.4078 (df = 64), p-value = 0.0000 #> I-square statistic = 53.4% #> #> 17 studies, 68 observations, 4 fixed and 10 random-effects parameters #> logLik AIC BIC #> -92.2480 212.4959 242.7203 #>
# maximum likelihood estimators from the metaSEM package library(metaSEM) metaSEM_RE <- summary(meta(y = y, v = S)) metaSEM_RE
#> #> Call: #> meta(y = y, v = S) #> #> 95% confidence intervals: z statistic approximation (robust=FALSE) #> Coefficients: #> Estimate Std.Error lbound ubound z value Pr(>|z|) #> Intercept1 -1.2594e+00 1.1408e+01 -2.3619e+01 2.1100e+01 -0.1104 0.9121 #> Intercept2 -1.6238e+00 1.5133e+01 -3.1283e+01 2.8036e+01 -0.1073 0.9145 #> Intercept3 1.3275e-02 5.8909e+02 -1.1546e+03 1.1546e+03 0.0000 1.0000 #> Intercept4 5.9969e-02 1.2581e+02 -2.4653e+02 2.4665e+02 0.0005 0.9996 #> Tau2_1_1 6.0518e-01 1.0287e+02 -2.0102e+02 2.0223e+02 0.0059 0.9953 #> Tau2_2_1 1.3066e+00 2.5264e+02 -4.9385e+02 4.9647e+02 0.0052 0.9959 #> Tau2_2_2 1.2662e+00 1.5463e+02 -3.0181e+02 3.0434e+02 0.0082 0.9935 #> Tau2_3_1 -6.0376e-03 1.1486e+04 -2.2513e+04 2.2513e+04 0.0000 1.0000 #> Tau2_3_2 -1.9688e-02 1.4137e+04 -2.7709e+04 2.7709e+04 0.0000 1.0000 #> Tau2_3_3 1.5113e-03 NA NA NA NA NA #> Tau2_4_1 -4.5208e-02 2.4796e+03 -4.8599e+03 4.8599e+03 0.0000 1.0000 #> Tau2_4_2 6.2169e-03 3.0646e+03 -6.0065e+03 6.0065e+03 0.0000 1.0000 #> Tau2_4_3 1.0743e-02 NA NA NA NA NA #> Tau2_4_4 2.2724e-03 1.4893e+04 -2.9189e+04 2.9189e+04 0.0000 1.0000 #> #> Q statistic on the homogeneity of effect sizes: 137.4078 #> Degrees of freedom of the Q statistic: 64 #> P value of the Q statistic: 2.790395e-07 #> #> Heterogeneity indices (based on the estimated Tau2): #> Estimate #> Intercept1: I2 (Q statistic) 0.0340 #> Intercept2: I2 (Q statistic) 0.2017 #> Intercept3: I2 (Q statistic) 0.1992 #> Intercept4: I2 (Q statistic) 0.0127 #> #> Number of studies (or clusters): 17 #> Number of observed statistics: 68 #> Number of estimated parameters: 14 #> Degrees of freedom: 54 #> -2 log likelihood: 194.1707 #> OpenMx status1: 5 ("0" or "1": The optimization is considered fine. #> Other values may indicate problems.)
#> Warning: OpenMx status1 is neither 0 or 1. You are advised to 'rerun' it again.
############################################################## # Plotting the result: ############################################################## obj <- MMA_FE # obj <- mvmeta_RE # obj <- metaSEM_RE # pdf("CI.pdf", width = 4, height = 7) plotCI(y = computvcov$ef, v = computvcov$list.vcov, name.y = NULL, name.study = Geeganage2010$studyID, y.all = obj$coefficients[,1], y.all.se = obj$coefficients[,2], hline = c(0, 0, 0, 1))
#> $`Plotting MD.SBP`
#> #> $`Plotting MD.DBP`
#> #> $`Plotting RD.DD`
#> #> $`Plotting lgOR.D`
#>
# dev.off()