Functions to construct a matrix of pairwise contrasts for objects of class "rma".

pairwise(x, btt, btt2, ...)

Arguments

x

an object of class "rma".

btt

vector of indices to specify for which coefficients pairwise contrasts should be constructed. Can also be a string to grep for. See ‘Details’.

btt2

optional argument to specify a second set of coefficients that should also be included in the contrast matrix.

...

other arguments.

Value

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.

References

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

See also

rma.uni, rma.glmm, and rma.mv for functions to fit meta-regression models for which pairwise contrasts may be useful.

anova for a function to carry out tests of the pairwise contrasts and predict to obtain corresponding confidence/prediction intervals.

Examples

### 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
pairwise(res, btt=1:3)
#>                                [,1] [,2] [,3]
#> allocalternate-allocrandom       -1    1    0
#> allocalternate-allocsystematic   -1    0    1
#> allocrandom-allocsystematic       0   -1    1
pairwise(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=pairwise(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=pairwise(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=pairwise(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=pairwise(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)
pairwise(res, btt=2:3, btt2=2:3)
#>                             [,1] [,2] [,3]
#> allocrandom                    0    1    0
#> allocsystematic                0    0    1
#> allocrandom-allocsystematic    0   -1    1
pairwise(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=pairwise(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=pairwise(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=pairwise(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=pairwise(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=pairwise(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=pairwise(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=pairwise(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=pairwise(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=pairwise(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   
#>