deltamethod.RdFunction to apply the (multivariate) delta method to a set of estimates.
deltamethod(x, vcov, fun, level, H0=0, digits)either a vector of estimates or a model object from which model coefficients can be extracted via coef(x).
when x is a vector of estimates, the corresponding variance-covariance matrix (ignored when x is a model object, in which case vcov(x) is used to extract the variance-covariance matrix).
a function to apply to the estimates.
numeric value between 0 and 100 to specify the confidence interval level (see here for details). If unspecified, this either defaults to 95 or, if possible, the corresponding value from the model object.
numeric value to specify the value under the null hypothesis for the Wald-type test(s) (the default is 0). Can also be a vector.
optional integer to specify the number of decimal places to which the printed results should be rounded.
Let \(\hat{\theta}\) denote a vector of \(p\) estimates which can be specified via the x argument and let \(\Sigma\) denote the corresponding \(p \times p\) variance-covariance matrix, which can be specified via the vcov argument. If x is not an vector with estimates, then the function assumes that x is a model object and will try to use coef(x) and vcov(x) to extract the model coefficients and the corresponding variance-covariance matrix (in this case, the vcov argument is ignored).
Let \(f(\cdot)\) be a function, specified via the fun argument, with \(p\) inputs/arguments (or with a single argument that is assumed to be a vector of length \(p\)), which returns a numeric (and atomic) vector of \(q\) transformed estimates. Then the function computes \(f(\hat{\theta})\) and the corresponding variance-covariance matrix of the transformed estimates using the multivariate delta method (e.g., van der Vaart, 1998) with \[\textrm{Var}[f(\hat{\theta})] = \nabla f(\hat{\theta})' \cdot \Sigma \cdot \nabla f(\hat{\theta})\] where \(\nabla f(\hat{\theta})\) denotes the gradient of \(f(\cdot)\) evaluated at \(\hat{\theta}\). The function computes the gradient numerically using the derivative function from the calculus package.
The function also computes Wald-type tests and confidence intervals for the \(q\) transformed estimates. The level argument can be used to control the confidence interval level.
An object of class "deltamethod". The object is a list containing the following components:
a data frame with the transformed estimates, standard errors, test statistics, p-values, and lower/upper confidence interval bounds.
the variance-covariance matrix of the transformed estimates.
some additional elements/values.
The results are formatted and printed with the print function. Extractor functions include coef and vcov.
van der Vaart, A. W. (1998). Asymptotic statistics. Cambridge, UK: Cambridge University Press.
Viechtbauer, W. (2010). Conducting meta-analyses in R with the metafor package. Journal of Statistical Software, 36(3), 1–48. https://doi.org/10.18637/jss.v036.i03
conv.delta for a function to apply the (univariate) delta method to observed effect sizes or outcomes and their sampling variances.
############################################################################
### copy data into 'dat'
dat <- dat.craft2003
### construct dataset and var-cov matrix of the correlations
tmp <- rcalc(ri ~ var1 + var2 | study, ni=ni, data=dat)
V <- tmp$V
dat <- tmp$dat
### turn var1.var2 into a factor with the desired order of levels
dat$var1.var2 <- factor(dat$var1.var2,
levels=c("acog.perf", "asom.perf", "conf.perf", "acog.asom", "acog.conf", "asom.conf"))
### multivariate random-effects model
res <- rma.mv(yi, V, mods = ~ 0 + var1.var2, random = ~ var1.var2 | study, struct="UN", data=dat)
#> Warning: 9 rows with NAs omitted from model fitting.
res
#>
#> Multivariate Meta-Analysis Model (k = 51; method: REML)
#>
#> Variance Components:
#>
#> outer factor: study (nlvls = 9)
#> inner factor: var1.var2 (nlvls = 6)
#>
#> estim sqrt k.lvl fixed level
#> tau^2.1 0.1611 0.4014 9 no acog.perf
#> tau^2.2 0.0604 0.2459 9 no asom.perf
#> tau^2.3 0.0468 0.2163 8 no conf.perf
#> tau^2.4 0.0047 0.0683 9 no acog.asom
#> tau^2.5 0.0125 0.1119 8 no acog.conf
#> tau^2.6 0.0111 0.1052 8 no asom.conf
#>
#> rho.acg.p rho.asm.p rho.cnf. rho.acg.s rho.acg.c rho.asm.c acg.p asm.p cnf.
#> acog.perf 1 - 9 8
#> asom.perf 0.9497 1 no - 8
#> conf.perf -0.6178 -0.5969 1 no no -
#> acog.asom 0.5491 0.4604 -0.9345 1 no no no
#> acog.conf 0.0432 -0.0495 0.7023 -0.6961 1 no no no
#> asom.conf 0.3532 0.2688 -0.1311 -0.0891 0.4193 1 no no no
#> acg.s acg.c asm.c
#> acog.perf 9 8 8
#> asom.perf 9 8 8
#> conf.perf 8 8 8
#> acog.asom - 8 8
#> acog.conf no - 8
#> asom.conf no no -
#>
#> Test for Residual Heterogeneity:
#> QE(df = 45) = 334.8358, p-val < .0001
#>
#> Test of Moderators (coefficients 1:6):
#> QM(df = 6) = 596.7705, p-val < .0001
#>
#> Model Results:
#>
#> estimate se zval pval ci.lb ci.ub
#> var1.var2acog.perf -0.0600 0.1408 -0.4264 0.6698 -0.3359 0.2159
#> var1.var2asom.perf -0.1423 0.0917 -1.5527 0.1205 -0.3220 0.0373
#> var1.var2conf.perf 0.3167 0.0847 3.7393 0.0002 0.1507 0.4827 ***
#> var1.var2acog.asom 0.5671 0.0367 15.4640 <.0001 0.4953 0.6390 ***
#> var1.var2acog.conf -0.4888 0.0509 -9.6048 <.0001 -0.5886 -0.3891 ***
#> var1.var2asom.conf -0.4750 0.0506 -9.3901 <.0001 -0.5741 -0.3758 ***
#>
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#>
### restructure estimated mean correlations into a 4x4 matrix
R <- vec2mat(coef(res))
rownames(R) <- colnames(R) <- c("perf", "acog", "asom", "conf")
round(R, digits=3)
#> perf acog asom conf
#> perf 1.000 -0.060 -0.142 0.317
#> acog -0.060 1.000 0.567 -0.489
#> asom -0.142 0.567 1.000 -0.475
#> conf 0.317 -0.489 -0.475 1.000
### check that order in vcov(res) corresponds to order in R
round(vcov(res), digits=4)
#> var1.var2acog.perf var1.var2asom.perf var1.var2conf.perf var1.var2acog.asom
#> var1.var2acog.perf 0.0198 0.0115 -0.0069 0.0017
#> var1.var2asom.perf 0.0115 0.0084 -0.0043 0.0009
#> var1.var2conf.perf -0.0069 -0.0043 0.0072 -0.0017
#> var1.var2acog.asom 0.0017 0.0009 -0.0017 0.0013
#> var1.var2acog.conf 0.0004 -0.0002 0.0023 -0.0009
#> var1.var2asom.conf 0.0018 0.0010 -0.0004 -0.0004
#> var1.var2acog.conf var1.var2asom.conf
#> var1.var2acog.perf 0.0004 0.0018
#> var1.var2asom.perf -0.0002 0.0010
#> var1.var2conf.perf 0.0023 -0.0004
#> var1.var2acog.asom -0.0009 -0.0004
#> var1.var2acog.conf 0.0026 0.0011
#> var1.var2asom.conf 0.0011 0.0026
### fit regression model with 'perf' as outcome and 'acog', 'asom', and 'conf' as predictors
matreg(1, 2:4, R=R, V=vcov(res))
#>
#> estimate se zval pval ci.lb ci.ub
#> acog 0.1482 0.1566 0.9465 0.3439 -0.1587 0.4550
#> asom -0.0536 0.0768 -0.6979 0.4852 -0.2043 0.0970
#> conf 0.3637 0.0910 3.9985 <.0001 0.1854 0.5419 ***
#>
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#>
### same analysis but using the deltamethod() function
deltamethod(coef(res), vcov(res), fun=function(r1,r2,r3,r4,r5,r6) {
R <- vec2mat(c(r1,r2,r3,r4,r5,r6))
setNames(c(solve(R[-1,-1]) %*% R[2:4,1]), c("acog","asom","conf"))
})
#>
#> estimate se zval pval ci.lb ci.ub
#> acog 0.1482 0.1566 0.9465 0.3439 -0.1587 0.4550
#> asom -0.0536 0.0768 -0.6979 0.4852 -0.2043 0.0970
#> conf 0.3637 0.0910 3.9985 <.0001 0.1854 0.5419 ***
#>
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#>
### using a function that takes a vector as input
deltamethod(coef(res), vcov(res), fun=function(r) {
R <- vec2mat(r)
setNames(c(solve(R[-1,-1]) %*% R[2:4,1]), c("acog","asom","conf"))
})
#>
#> estimate se zval pval ci.lb ci.ub
#> acog 0.1482 0.1566 0.9465 0.3439 -0.1587 0.4550
#> asom -0.0536 0.0768 -0.6979 0.4852 -0.2043 0.0970
#> conf 0.3637 0.0910 3.9985 <.0001 0.1854 0.5419 ***
#>
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#>
############################################################################
### construct dataset and var-cov matrix of the r-to-z transformed correlations
dat <- dat.craft2003
tmp <- rcalc(ri ~ var1 + var2 | study, ni=ni, data=dat, rtoz=TRUE)
V <- tmp$V
dat <- tmp$dat
### turn var1.var2 into a factor with the desired order of levels
dat$var1.var2 <- factor(dat$var1.var2,
levels=c("acog.perf", "asom.perf", "conf.perf", "acog.asom", "acog.conf", "asom.conf"))
### multivariate random-effects model
res <- rma.mv(yi, V, mods = ~ 0 + var1.var2, random = ~ var1.var2 | study, struct="UN", data=dat)
#> Warning: 9 rows with NAs omitted from model fitting.
res
#>
#> Multivariate Meta-Analysis Model (k = 51; method: REML)
#>
#> Variance Components:
#>
#> outer factor: study (nlvls = 9)
#> inner factor: var1.var2 (nlvls = 6)
#>
#> estim sqrt k.lvl fixed level
#> tau^2.1 0.1747 0.4180 9 no acog.perf
#> tau^2.2 0.0581 0.2410 9 no asom.perf
#> tau^2.3 0.0594 0.2438 8 no conf.perf
#> tau^2.4 0.0095 0.0974 9 no acog.asom
#> tau^2.5 0.0166 0.1290 8 no acog.conf
#> tau^2.6 0.0110 0.1048 8 no asom.conf
#>
#> rho.acg.p rho.asm.p rho.cnf. rho.acg.s rho.acg.c rho.asm.c acg.p asm.p cnf.
#> acog.perf 1 - 9 8
#> asom.perf 0.9426 1 no - 8
#> conf.perf -0.5875 -0.5611 1 no no -
#> acog.asom 0.6497 0.4954 -0.9122 1 no no no
#> acog.conf -0.0622 -0.1297 0.8080 -0.6220 1 no no no
#> asom.conf 0.2732 0.2499 -0.1899 0.1025 0.1724 1 no no no
#> acg.s acg.c asm.c
#> acog.perf 9 8 8
#> asom.perf 9 8 8
#> conf.perf 8 8 8
#> acog.asom - 8 8
#> acog.conf no - 8
#> asom.conf no no -
#>
#> Test for Residual Heterogeneity:
#> QE(df = 45) = 272.3337, p-val < .0001
#>
#> Test of Moderators (coefficients 1:6):
#> QM(df = 6) = 400.3522, p-val < .0001
#>
#> Model Results:
#>
#> estimate se zval pval ci.lb ci.ub
#> var1.var2acog.perf -0.0746 0.1493 -0.4995 0.6174 -0.3672 0.2180
#> var1.var2asom.perf -0.1597 0.0927 -1.7226 0.0850 -0.3414 0.0220 .
#> var1.var2conf.perf 0.3460 0.0965 3.5839 0.0003 0.1568 0.5352 ***
#> var1.var2acog.asom 0.6203 0.0526 11.7875 <.0001 0.5172 0.7234 ***
#> var1.var2acog.conf -0.5135 0.0623 -8.2369 <.0001 -0.6356 -0.3913 ***
#> var1.var2asom.conf -0.4872 0.0580 -8.3964 <.0001 -0.6010 -0.3735 ***
#>
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#>
### estimate the difference between r(acog,perf) and r(asom,perf)
deltamethod(res, fun=function(z1,z2,z3,z4,z5,z6) {
transf.ztor(z1) - transf.ztor(z2)
})
#>
#> estimate se zval pval ci.lb ci.ub
#> 0.0839 0.0828 1.0137 0.3107 -0.0784 0.2462
#>
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#>
### using a function that takes a vector as input
deltamethod(res, fun=function(z) {
transf.ztor(z[1]) - transf.ztor(z[2])
})
#>
#> estimate se zval pval ci.lb ci.ub
#> 0.0839 0.0828 1.0137 0.3107 -0.0784 0.2462
#>
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#>
############################################################################