pairmat.RdFunctions to construct a matrix of pairwise contrasts for objects of class "rma".
pairmat(x, btt, btt2, ...)an object of class "rma".
vector of indices to specify for which coefficients pairwise contrasts should be constructed. Can also be a string to grep for. See ‘Details’.
optional argument to specify a second set of coefficients that should also be included in the contrast matrix.
other arguments.
When a meta-regression model includes a categorical moderator variable (i.e., a factor), there is often interest in testing whether the coefficients representing the various levels of the factor differ significantly from each other. The present function constructs the pairwise contrast matrix between all factor levels for a particular factor, which can be used together with the anova function to carry out such tests and the predict function to obtain corresponding confidence intervals.
The x argument is used to specify a meta-regression model and the btt argument the indices of the coefficients for which pairwise contrasts should be constructed. For example, with btt=2:4, contrasts are formed based on the second, third, and fourth coefficient of the model. Instead of specifying the coefficient numbers, one can specify a string for btt. In that case, grep will be used to search for all coefficient names that match the string.
At times, it may be useful to include a second set of coefficients in the contrast matrix (not as pairwise contrasts, but as ‘main effects’). This can be done via the btt2 argument.
When using the present function in a call to the anova or predict functions, argument x does not need to specified, as the function will then automatically construct the contrast matrix based on the model object passed to the anova or predict function. See below for examples.
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
### calculate log risk ratios and corresponding sampling variances
dat <- escalc(measure="RR", ai=tpos, bi=tneg, ci=cpos, di=cneg, data=dat.bcg)
### mixed-effects meta-regression model with the allocation method as a moderator;
### by removing the intercept term, we obtain the estimated average effect for each
### factor level from the model
res <- rma(yi, vi, mods = ~ 0 + alloc, data=dat)
res
#>
#> Mixed-Effects Model (k = 13; tau^2 estimator: REML)
#>
#> tau^2 (estimated amount of residual heterogeneity): 0.3615 (SE = 0.2111)
#> tau (square root of estimated tau^2 value): 0.6013
#> I^2 (residual heterogeneity / unaccounted variability): 88.77%
#> H^2 (unaccounted variability / sampling variability): 8.91
#>
#> Test for Residual Heterogeneity:
#> QE(df = 10) = 132.3676, p-val < .0001
#>
#> Test of Moderators (coefficients 1:3):
#> QM(df = 3) = 15.9842, p-val = 0.0011
#>
#> Model Results:
#>
#> estimate se zval pval ci.lb ci.ub
#> allocalternate -0.5180 0.4412 -1.1740 0.2404 -1.3827 0.3468
#> allocrandom -0.9658 0.2672 -3.6138 0.0003 -1.4896 -0.4420 ***
#> allocsystematic -0.4289 0.3449 -1.2434 0.2137 -1.1050 0.2472
#>
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#>
### construct the contrast matrix for the 'alloc' factor
pairmat(res, btt=1:3)
#> [,1] [,2] [,3]
#> allocalternate-allocrandom -1 1 0
#> allocalternate-allocsystematic -1 0 1
#> allocrandom-allocsystematic 0 -1 1
pairmat(res, btt="alloc")
#> [,1] [,2] [,3]
#> allocalternate-allocrandom -1 1 0
#> allocalternate-allocsystematic -1 0 1
#> allocrandom-allocsystematic 0 -1 1
### test all pairwise contrasts
anova(res, X=pairmat(btt=1:3))
#>
#> Hypotheses:
#> 1: -allocalternate + allocrandom = 0
#> 2: -allocalternate + allocsystematic = 0
#> 3: -allocrandom + allocsystematic = 0
#>
#> Results:
#> estimate se zval pval
#> 1: -0.4478 0.5158 -0.8682 0.3853
#> 2: 0.0890 0.5600 0.1590 0.8737
#> 3: 0.5369 0.4364 1.2303 0.2186
#>
anova(res, X=pairmat(btt="alloc"))
#>
#> Hypotheses:
#> 1: -allocalternate + allocrandom = 0
#> 2: -allocalternate + allocsystematic = 0
#> 3: -allocrandom + allocsystematic = 0
#>
#> Results:
#> estimate se zval pval
#> 1: -0.4478 0.5158 -0.8682 0.3853
#> 2: 0.0890 0.5600 0.1590 0.8737
#> 3: 0.5369 0.4364 1.2303 0.2186
#>
### obtain the corresponding confidence intervals
predict(res, newmods=pairmat(btt="alloc"))
#>
#> pred se ci.lb ci.ub pi.lb pi.ub
#> allocalternate-allocrandom -0.4478 0.5158 -1.4588 0.5632 -2.0005 1.1049
#> allocalternate-allocsystematic 0.0890 0.5600 -1.0086 1.1867 -1.5214 1.6995
#> allocrandom-allocsystematic 0.5369 0.4364 -0.3184 1.3921 -0.9192 1.9929
#>
### test all pairwise contrasts adjusting for multiple testing
anova(res, X=pairmat(btt="alloc"), adjust="bonf")
#>
#> Hypotheses:
#> 1: -allocalternate + allocrandom = 0
#> 2: -allocalternate + allocsystematic = 0
#> 3: -allocrandom + allocsystematic = 0
#>
#> Results:
#> estimate se zval pval
#> 1: -0.4478 0.5158 -0.8682 1.0000
#> 2: 0.0890 0.5600 0.1590 1.0000
#> 3: 0.5369 0.4364 1.2303 0.6557
#>
### fit the same model, but including the intercept term; then 'alternate' is the
### reference level and the coefficients for 'random' and 'systematic' already
### represent pairwise contrasts with this reference level
res <- rma(yi, vi, mods = ~ alloc, data=dat)
res
#>
#> Mixed-Effects Model (k = 13; tau^2 estimator: REML)
#>
#> tau^2 (estimated amount of residual heterogeneity): 0.3615 (SE = 0.2111)
#> tau (square root of estimated tau^2 value): 0.6013
#> I^2 (residual heterogeneity / unaccounted variability): 88.77%
#> H^2 (unaccounted variability / sampling variability): 8.91
#> R^2 (amount of heterogeneity accounted for): 0.00%
#>
#> Test for Residual Heterogeneity:
#> QE(df = 10) = 132.3676, p-val < .0001
#>
#> Test of Moderators (coefficients 2:3):
#> QM(df = 2) = 1.7675, p-val = 0.4132
#>
#> Model Results:
#>
#> estimate se zval pval ci.lb ci.ub
#> intrcpt -0.5180 0.4412 -1.1740 0.2404 -1.3827 0.3468
#> allocrandom -0.4478 0.5158 -0.8682 0.3853 -1.4588 0.5632
#> allocsystematic 0.0890 0.5600 0.1590 0.8737 -1.0086 1.1867
#>
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#>
### in this case, we want to include these coefficients directly in the contrast
### matrix (btt2=2:3) but also include the pairwise contrast between them (btt=2:3)
pairmat(res, btt=2:3, btt2=2:3)
#> [,1] [,2] [,3]
#> allocrandom 0 1 0
#> allocsystematic 0 0 1
#> allocrandom-allocsystematic 0 -1 1
pairmat(res, btt="alloc", btt2="alloc")
#> [,1] [,2] [,3]
#> allocrandom 0 1 0
#> allocsystematic 0 0 1
#> allocrandom-allocsystematic 0 -1 1
### test all pairwise contrasts
anova(res, X=pairmat(btt=2:3, btt2=2:3))
#>
#> Hypotheses:
#> 1: allocrandom = 0
#> 2: allocsystematic = 0
#> 3: -allocrandom + allocsystematic = 0
#>
#> Results:
#> estimate se zval pval
#> 1: -0.4478 0.5158 -0.8682 0.3853
#> 2: 0.0890 0.5600 0.1590 0.8737
#> 3: 0.5369 0.4364 1.2303 0.2186
#>
anova(res, X=pairmat(btt="alloc", btt2="alloc"))
#>
#> Hypotheses:
#> 1: allocrandom = 0
#> 2: allocsystematic = 0
#> 3: -allocrandom + allocsystematic = 0
#>
#> Results:
#> estimate se zval pval
#> 1: -0.4478 0.5158 -0.8682 0.3853
#> 2: 0.0890 0.5600 0.1590 0.8737
#> 3: 0.5369 0.4364 1.2303 0.2186
#>
### obtain the corresponding confidence intervals
predict(res, newmods=pairmat(btt="alloc", btt2="alloc"))
#>
#> pred se ci.lb ci.ub pi.lb pi.ub
#> allocrandom -0.4478 0.5158 -1.4588 0.5632 -2.0005 1.1049
#> allocsystematic 0.0890 0.5600 -1.0086 1.1867 -1.5214 1.6995
#> allocrandom-allocsystematic 0.5369 0.4364 -0.3184 1.3921 -0.9192 1.9929
#>
### meta-regression model with 'ablat' and 'alloc' as moderators
res <- rma(yi, vi, mods = ~ ablat + alloc, data=dat)
res
#>
#> Mixed-Effects Model (k = 13; tau^2 estimator: REML)
#>
#> tau^2 (estimated amount of residual heterogeneity): 0.1446 (SE = 0.1124)
#> tau (square root of estimated tau^2 value): 0.3803
#> I^2 (residual heterogeneity / unaccounted variability): 70.11%
#> H^2 (unaccounted variability / sampling variability): 3.35
#> R^2 (amount of heterogeneity accounted for): 53.84%
#>
#> Test for Residual Heterogeneity:
#> QE(df = 9) = 26.2034, p-val = 0.0019
#>
#> Test of Moderators (coefficients 2:4):
#> QM(df = 3) = 11.0605, p-val = 0.0114
#>
#> Model Results:
#>
#> estimate se zval pval ci.lb ci.ub
#> intrcpt 0.2932 0.4050 0.7239 0.4691 -0.5006 1.0870
#> ablat -0.0273 0.0092 -2.9650 0.0030 -0.0453 -0.0092 **
#> allocrandom -0.2675 0.3504 -0.7633 0.4453 -0.9543 0.4193
#> allocsystematic 0.0585 0.3795 0.1540 0.8776 -0.6854 0.8023
#>
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#>
### test all pairwise contrasts between the 'alloc' levels (while controlling for 'ablat')
anova(res, X=pairmat(btt="alloc", btt2="alloc"))
#>
#> Hypotheses:
#> 1: allocrandom = 0
#> 2: allocsystematic = 0
#> 3: -allocrandom + allocsystematic = 0
#>
#> Results:
#> estimate se zval pval
#> 1: -0.2675 0.3504 -0.7633 0.4453
#> 2: 0.0585 0.3795 0.1540 0.8776
#> 3: 0.3260 0.3104 1.0501 0.2937
#>
anova(res, X=pairmat(btt="alloc", btt2="alloc"))
#>
#> Hypotheses:
#> 1: allocrandom = 0
#> 2: allocsystematic = 0
#> 3: -allocrandom + allocsystematic = 0
#>
#> Results:
#> estimate se zval pval
#> 1: -0.2675 0.3504 -0.7633 0.4453
#> 2: 0.0585 0.3795 0.1540 0.8776
#> 3: 0.3260 0.3104 1.0501 0.2937
#>
### obtain the corresponding confidence intervals
predict(res, newmods=pairmat(btt="alloc", btt2="alloc"))
#>
#> pred se ci.lb ci.ub pi.lb pi.ub
#> allocrandom -0.2675 0.3504 -0.9543 0.4193 -1.2810 0.7460
#> allocsystematic 0.0585 0.3795 -0.6854 0.8023 -0.9945 1.1114
#> allocrandom-allocsystematic 0.3260 0.3104 -0.2824 0.9343 -0.6361 1.2880
#>
### an example of a meta-regression model with more factors levels
dat <- dat.bangertdrowns2004
res <- rma(yi, vi, mods = ~ 0 + factor(grade), data=dat)
res
#>
#> Mixed-Effects Model (k = 48; tau^2 estimator: REML)
#>
#> tau^2 (estimated amount of residual heterogeneity): 0.0539 (SE = 0.0216)
#> tau (square root of estimated tau^2 value): 0.2322
#> I^2 (residual heterogeneity / unaccounted variability): 59.15%
#> H^2 (unaccounted variability / sampling variability): 2.45
#>
#> Test for Residual Heterogeneity:
#> QE(df = 44) = 102.0036, p-val < .0001
#>
#> Test of Moderators (coefficients 1:4):
#> QM(df = 4) = 28.5360, p-val < .0001
#>
#> Model Results:
#>
#> estimate se zval pval ci.lb ci.ub
#> factor(grade)1 0.2639 0.0898 2.9393 0.0033 0.0879 0.4399 **
#> factor(grade)2 -0.1088 0.1450 -0.7504 0.4530 -0.3929 0.1753
#> factor(grade)3 0.2888 0.1027 2.8111 0.0049 0.0874 0.4901 **
#> factor(grade)4 0.2484 0.0735 3.3810 0.0007 0.1044 0.3924 ***
#>
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#>
### test all pairwise contrasts between the 'grade' levels
anova(res, X=pairmat(btt="grade"))
#>
#> Hypotheses:
#> 1: -factor(grade)1 + factor(grade)2 = 0
#> 2: -factor(grade)1 + factor(grade)3 = 0
#> 3: -factor(grade)1 + factor(grade)4 = 0
#> 4: -factor(grade)2 + factor(grade)3 = 0
#> 5: -factor(grade)2 + factor(grade)4 = 0
#> 6: -factor(grade)3 + factor(grade)4 = 0
#>
#> Results:
#> estimate se zval pval
#> 1: -0.3727 0.1705 -2.1856 0.0288 *
#> 2: 0.0248 0.1364 0.1821 0.8555
#> 3: -0.0155 0.1160 -0.1338 0.8935
#> 4: 0.3975 0.1777 2.2375 0.0253 *
#> 5: 0.3572 0.1625 2.1977 0.0280 *
#> 6: -0.0404 0.1263 -0.3197 0.7492
#>
### obtain the corresponding confidence intervals
predict(res, newmods=pairmat(btt="grade"))
#>
#> pred se ci.lb ci.ub pi.lb pi.ub
#> factor(grade)1-factor(grade)2 -0.3727 0.1705 -0.7069 -0.0385 -0.9374 0.1920
#> factor(grade)1-factor(grade)3 0.0248 0.1364 -0.2425 0.2922 -0.5031 0.5528
#> factor(grade)1-factor(grade)4 -0.0155 0.1160 -0.2429 0.2118 -0.5244 0.4933
#> factor(grade)2-factor(grade)3 0.3975 0.1777 0.0493 0.7458 -0.1756 0.9707
#> factor(grade)2-factor(grade)4 0.3572 0.1625 0.0386 0.6757 -0.1984 0.9127
#> factor(grade)3-factor(grade)4 -0.0404 0.1263 -0.2879 0.2071 -0.5585 0.4778
#>
### test all pairwise contrasts adjusting for multiple testing
anova(res, X=pairmat(btt="grade"), adjust="bonf")
#>
#> Hypotheses:
#> 1: -factor(grade)1 + factor(grade)2 = 0
#> 2: -factor(grade)1 + factor(grade)3 = 0
#> 3: -factor(grade)1 + factor(grade)4 = 0
#> 4: -factor(grade)2 + factor(grade)3 = 0
#> 5: -factor(grade)2 + factor(grade)4 = 0
#> 6: -factor(grade)3 + factor(grade)4 = 0
#>
#> Results:
#> estimate se zval pval
#> 1: -0.3727 0.1705 -2.1856 0.1731
#> 2: 0.0248 0.1364 0.1821 1.0000
#> 3: -0.0155 0.1160 -0.1338 1.0000
#> 4: 0.3975 0.1777 2.2375 0.1515
#> 5: 0.3572 0.1625 2.1977 0.1678
#> 6: -0.0404 0.1263 -0.3197 1.0000
#>