The metafor Package

A Meta-Analysis Package for R

User Tools

Site Tools

tips:computing_adjusted_effects

Differences

This shows you the differences between two versions of the page.

Link to this comparison view

Both sides previous revisionPrevious revision
Next revision		Previous revision
Next revisionBoth sides next revision
tips:computing_adjusted_effects [2020/08/18 19:53] – Wolfgang Viechtbauertips:computing_adjusted_effects [2022/07/09 09:57] – Wolfgang Viechtbauer
Line 10: Line 10:
 dat <- dat.bcg dat <- dat.bcg
 dat <- escalc(measure="RR", ai=tpos, bi=tneg, ci=cpos, di=cneg, data=dat, dat <- escalc(measure="RR", ai=tpos, bi=tneg, ci=cpos, di=cneg, data=dat,
-              slab=paste0(dat$author, ", ", dat$year))+              slab=paste0(author, ", ", year))
 dat dat
 </code> </code>
Line 37: Line 37:
  
 <code rsplus> <code rsplus>
-res <- rma(yi, vi, data=dat) +res1 <- rma(yi, vi, data=dat) 
-res+res1
 </code> </code>
 <code output> <code output>
Line 53: Line 53:
 Model Results: Model Results:
  
-estimate      se     zval    pval    ci.lb    ci.ub +estimate      se     zval    pval    ci.lb    ci.ub       
- -0.7145  0.1798  -3.9744  <.0001  -1.0669  -0.3622  ***+ -0.7145  0.1798  -3.9744  <.0001  -1.0669  -0.3622  *** 
  
 --- ---
-Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1+Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
 </code> </code>
  
-To make the results easier to interpret, we can exponentiate the estimated average log risk ratio (i.e., $-0.7145$) to obtain an estimate of the average risk ratio. This can be done with the ''predict()'' function as follows.+To make the results easier to interpret, we can exponentiate the estimated average log risk ratio (i.e., $-0.7145$) to obtain an estimate of the average risk ratio.((Strictly speaking, since exponentiation is a non-linear transformation, the back-transformed value is an estimate of the //median// risk ratio, but we will gloss over this fact.)) This can be done with the ''predict()'' function as follows.
  
 <code rsplus> <code rsplus>
-predict(res, transf=exp, digits=2)+predict(res1, transf=exp, digits=2)
 </code> </code>
 <code output> <code output>
Line 75: Line 75:
  
 <code rsplus> <code rsplus>
-forest(res, xlim=c(-6.8,3.8), header=TRUE, atransf=exp,+forest(res1, xlim=c(-6.8,3.8), header=TRUE, atransf=exp,
        at=log(c(1/16, 1/8, 1/4, 1/2, 1, 2, 4, 8)), digits=c(2L,4L))        at=log(c(1/16, 1/8, 1/4, 1/2, 1, 2, 4, 8)), digits=c(2L,4L))
 </code> </code>
Line 86: Line 86:
  
 <code rsplus> <code rsplus>
-res <- rma(yi, vi, mods = ~ ablat, data=dat) +res2 <- rma(yi, vi, mods = ~ ablat, data=dat) 
-res+res2
 </code> </code>
 <code output> <code output>
Line 106: Line 106:
 Model Results: Model Results:
  
-         estimate      se     zval    pval    ci.lb    ci.ub +         estimate      se     zval    pval    ci.lb    ci.ub       
-intrcpt    0.2515  0.2491   1.0095  0.3127  -0.2368   0.7397 +intrcpt    0.2515  0.2491   1.0095  0.3127  -0.2368   0.7397       
-ablat     -0.0291  0.0072  -4.0444  <.0001  -0.0432  -0.0150  ***+ablat     -0.0291  0.0072  -4.0444  <.0001  -0.0432  -0.0150  *** 
  
 --- ---
-Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1+Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
 </code> </code>
  
-The results indeed suggest that the size of the average log risk ratio depends on absolute latitude.+The results indeed suggest that the size of the average log risk ratio is a function of the absolute latitude of the study locations.
  
-In such a meta-regression model, there is no longer **the** average effect, so if we want to estimate the size of the effect, we have to specify a value for ''ablat''. One possibility is to estimate the average risk ratio based on the average absolute latitude of the 13 trials included in the meta-analysis. We can do this as follows.+In such a meta-regression model, there is no longer a single average effect (since it depends on the value of the moderator), so if we want to estimate the size of the effect, we have to specify a value for ''ablat''. One possibility is to estimate the average risk ratio based on the average absolute latitude of the 13 trials included in the meta-analysis. We can do this as follows.
  
 <code rsplus> <code rsplus>
-predict(res, newmods = mean(dat$ablat), transf=exp, digits=2)+predict(res2, newmods = mean(dat$ablat), transf=exp, digits=2)
 </code> </code>
 <code output> <code output>
Line 131: Line 131:
  
 <code rsplus> <code rsplus>
-forest(res, xlim=c(-6.8,3.8), header=TRUE, atransf=exp,+forest(res2, xlim=c(-6.8,3.8), header=TRUE, atransf=exp,
        at=log(c(1/16, 1/8, 1/4, 1/2, 1, 2, 4, 8)), digits=c(2L,4L),        at=log(c(1/16, 1/8, 1/4, 1/2, 1, 2, 4, 8)), digits=c(2L,4L),
-       ilab=dat$ablat, ilab.xpos=-3.5, order=order(dat$ablat), ylim=c(-1.5,15))+       ilab=ablat, ilab.xpos=-3.5, order=ablat, ylim=c(-2.5,15))
 text(-3.5, 15, "Lattitude", font=2) text(-3.5, 15, "Lattitude", font=2)
 abline(h=0) abline(h=0)
-sav <- predict(res, newmods = mean(dat$ablat)) +sav1 <- predict(res1) 
-addpoly(sav$pred, sei=sav$se, atransf=exp, digits=2, mlab="Adjusted Effect")+addpoly(sav1, row=-1, mlab="Random-Effects Model") 
 +sav2 <- predict(res2, newmods = mean(dat$ablat)) 
 +addpoly(sav2, row=-2, mlab="Meta-Regression Model (Adjusted Effect)")
 </code> </code>
  
-In the forest plot, the values of the moderator are added as an extra column and the studies are ordered based on their absolute latitude value to make this clearer. Some extra space was also added to the forest plot to add the 'adjusted effect' at the bottom.+In the forest plot, the values of the moderator are added as an extra column and the studies are ordered based on their absolute latitude value to make this clearer. Some extra space was also added to the forest plot to add the estimate from the standard random-effects model and the 'adjusted effect' from the meta-regression model at the bottom.
  
 {{ tips:adjusted_effects_forest_me.png?nolink }} {{ tips:adjusted_effects_forest_me.png?nolink }}
Line 149: Line 151:
  
 <code rsplus> <code rsplus>
-res <- rma(yi, vi, mods = ~ alloc, data=dat) +res3 <- rma(yi, vi, mods = ~ alloc, data=dat) 
-res+res3
 </code> </code>
 <code output> <code output>
Line 178: Line 180:
 </code> </code>
  
-The intercept reflects the estimated average log risk ratio for level ''alternate'' (the reference level), while the other two coefficients estimate how much the average log risk ratios for levels ''random'' and ''systematic'' differ from the reference level. Note that the omnibus test of these two coefficients is not significant ($Q_M = 1.77, df = 2, p = .41$), which indicates that there is insufficient evidence that the average log risk ratio actually differs across the three levels. However, for illustration purposes, we'll proceed with our analysis of this moderator.+The intercept reflects the estimated average log risk ratio for level ''alternate'' (the reference level), while the other two coefficients estimate how much the average log risk ratios for levels ''random'' and ''systematic'' differ from the reference level. Note that the omnibus test of these two coefficients is actually not significant ($Q_M = 1.77, \mbox{df} = 2, p = .41$), which indicates that there is insufficient evidence to conclude that the average log risk ratio actually differs across the three levels. However, for illustration purposes, we'll proceed with our analysis of this moderator.
  
 First, we can compute the estimated average risk ratios for the three levels with: First, we can compute the estimated average risk ratios for the three levels with:
  
 <code rsplus> <code rsplus>
-predict(res, newmods = c(0,0), transf=exp, digits=2) +predict(res3, newmods = c(0,0), transf=exp, digits=2) 
-predict(res, newmods = c(1,0), transf=exp, digits=2) +predict(res3, newmods = c(1,0), transf=exp, digits=2) 
-predict(res, newmods = c(0,1), transf=exp, digits=2)+predict(res3, newmods = c(0,1), transf=exp, digits=2)
 </code> </code>
 <code output> <code output>
Line 194: Line 196:
 </code> </code>
  
-Note: The output above was obtained with ''predict(res, newmods = rbind(c(0,0), c(1,0), c(0,1)), transf=exp, digits=2)'', which provides the three estimates in a single line of code. However, the code above is a bit easier to read and shows how we need to set the two dummy variables (that are created for the ''random'' and ''systematic'' levels) to 0 or 1 to obtain the estimates for the three levels.+Note: The output above was obtained with ''predict(res3, newmods = rbind(c(0,0), c(1,0), c(0,1)), transf=exp, digits=2)'', which provides the three estimates in a single line of code. However, the code above is a bit easier to read and shows how we need to set the two dummy variables (that are created for the ''random'' and ''systematic'' levels) to 0 or 1 to obtain the estimates for the three levels.
  
 But what should we do if we want to compute an 'adjusted effect' again? In other words, what values should we plug into our model equation (and hence into ''newmods'') to obtain such an estimate? One common approach is to use the mean of the respective dummy variables. We can obtain these values from the 'model matrix' and taking means across columns (leaving out the intercept). But what should we do if we want to compute an 'adjusted effect' again? In other words, what values should we plug into our model equation (and hence into ''newmods'') to obtain such an estimate? One common approach is to use the mean of the respective dummy variables. We can obtain these values from the 'model matrix' and taking means across columns (leaving out the intercept).
  
 <code rsplus> <code rsplus>
-colMeans(model.matrix(res))[-1]+colMeans(model.matrix(res3))[-1]
 </code> </code>
 <code output> <code output>
Line 209: Line 211:
  
 <code rsplus> <code rsplus>
-predict(res, newmods = colMeans(model.matrix(res))[-1], transf=exp, digits=2)+predict(res3, newmods = colMeans(model.matrix(res3))[-1], transf=exp, digits=2)
 </code> </code>
 <code output> <code output>
Line 223: Line 225:
  
 <code rsplus> <code rsplus>
-predict(res, newmods = c(1/3,1/3), transf=exp, digits=2)+predict(res3, newmods = c(1/3,1/3), transf=exp, digits=2)
 </code> </code>
 <code output> <code output>
Line 235: Line 237:
  
 <code rsplus> <code rsplus>
-res <- rma(yi, vi, mods = ~ ablat + alloc, data=dat) +res4 <- rma(yi, vi, mods = ~ ablat + alloc, data=dat) 
-res+res4
 </code> </code>
 <code output> <code output>
Line 268: Line 270:
  
 <code rsplus> <code rsplus>
-predict(res, newmods = colMeans(model.matrix(res))[-1], transf=exp, digits=2)+predict(res4, newmods = colMeans(model.matrix(res4))[-1], transf=exp, digits=2)
 </code> </code>
 <code output> <code output>
tips/computing_adjusted_effects.txt · Last modified: 2024/06/07 12:32 by Wolfgang Viechtbauer