-
Notifications
You must be signed in to change notification settings - Fork 21
/
Copy pathlipsey2001.rmd
212 lines (152 loc) · 8.37 KB
/
lipsey2001.rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
---
title: "R Code Corresponding to the Book *Practical Meta-Analysis* by Lipsey and Wilson (2001)"
author: |
| Wolfgang Viechtbauer
| Maastricht University
date: "`r format(Sys.time(), '%d %B, %Y')`"
output:
html_document:
# code_download: true
df_print: default
toc: true
number_sections: false
toc_depth: 3
toc_float:
collapsed: true
theme: default
# lots of nice themes can be used: https://bootswatch.com/
highlight: haddockadj.theme
# rmarkdown::github_document
# pdf_document:
# toc: true
# number_sections: false
# toc_depth: 3
# word_document
fig_caption: no
# bibliography: references.bib
---
```{r klippy, echo=FALSE, include=TRUE}
# remotes::install_github("rlesur/klippy")
klippy::klippy(position = c("top", "right"), color="gray20")
```
```{r crayon, echo=FALSE, message=FALSE, include=TRUE}
library(crayon)
options(crayon.enabled = TRUE)
knitr::knit_hooks$set(output = function(x, options){
paste0(
'<pre class="r-output"><code>',
fansi::sgr_to_html(x = htmltools::htmlEscape(x), warn = FALSE),
'</code></pre>'
)
})
```
## General Notes / Setup
The book *Practical Meta-Analysis* by Lipsey and Wilson (2001) is an excellent introduction to meta-analysis and covers (in less than 250 pages) the most important aspects of conducting a meta-analysis. This document provides the R code to reproduce the analyses conducted in chapter 7 (which is focused on the computational aspects of a meta-analysis, including the equal-, random-, and the mixed-effects model) using the `metafor` package. To read more about the package, see the [package website](https://www.metafor-project.org/) and the [package documentation](https://wviechtb.github.io/metafor/).
The package can be installed with:
```{r, eval=FALSE}
install.packages("metafor")
```
Once the package is installed, we can load it with:
```{r, message=FALSE}
library(metafor)
```
```{r, echo=FALSE, message=FALSE, warning=FALSE}
setmfopt(space=FALSE)
setmfopt(style=list(legend=make_style("gray90"), warning=strip_style))
pointsize <- 14
options(width=94)
```
***
## 1) Data
The data used for the analyses are provided in the book on page 130 in Table 7.1. We can create the same dataset with:
```{r}
dat <- data.frame(
id = c(100, 308, 1596, 2479, 9021, 9028, 161, 172, 537, 7049),
yi = c(-0.33, 0.32, 0.39, 0.31, 0.17, 0.64, -0.33, 0.15, -0.02, 0.00),
vi = c(0.084, 0.035, 0.017, 0.034, 0.072, 0.117, 0.102, 0.093, 0.012, 0.067),
random = c(0, 0, 0, 0, 0, 0, 1, 1, 1, 1),
intensity = c(7, 3, 7, 5, 7, 7, 4, 4, 5, 6))
dat
```
The `yi` values are standardized mean differences based on studies examining the effectiveness of 'challenge programs' to treat juvenile delinquency (presumably with positive values indicating that delinquent behaviors were reduced in the group of juveniles participating in such a challenge program compared to those who did not). The `vi` values are the corresponding sampling variances. The dummy variable `random` indicates whether juveniles were randomly assigned to the two conditions with a study (1 = yes, 0 = no). Finally, `intensity` is a variable coded for the purposes of the meta-analysis which indicates (based on assessment of the coder) the intensity of the challenge program (coded from 1 = 'very low' to 7 = 'very high'). These last two variables are potential moderators of the treatment effectiveness.
***
## 2) Equal-Effects Model
We can fit an equal-effects model to these data with:
```{r}
res.ee <- rma(yi, vi, data=dat, method="EE")
res.ee
```
These results match what is reported on page 132-133 (see Exhibit 7.3).
***
## 3) Random-Effects Model
A random-effects model (using the DerSimonian-Laird estimator for $\tau^2$) can be fitted with:
```{r}
res.re <- rma(yi, vi, data=dat, method="DL")
res.re
```
Again, these results match what is shown in Exhibit 7.3 and on pages 134-135.
A forest plot showing the results of the individual studies and the results from the random-effects model can be created as follows:
```{r, forestplot, fig.width=8, fig.height=6, dev.args=list(pointsize=pointsize), fig.align='center'}
par(mar=c(4,4,2,2))
forest(res.re, xlim=c(-2,3), header=c("Study", "SMD [95% CI]"), top=2,
xlab="Standardized Mean Difference")
```
***
## 4) Analog to the ANOVA
Next, we can carry out the ANOVA-type analysis described on pages 135-138. In particular, a fixed-effects meta-regression model with the `random` variable as moderator can be fitted with:
```{r}
res.anova <- rma(yi, vi, mods = ~ random, data=dat, method="FE")
res.anova
```
Note that the $Q_E$ and $Q_M$ statistics add up to the $Q$ statistic given earlier for the equal-effects models (i.e., $`r fmtx(res.anova$QE)` + `r fmtx(res.anova$QM)` = `r fmtx(res.ee$QE)`$).
The tests for heterogeneity within the two subgroups formed by the moderator can be obtained with:
```{r}
res.r0 <- rma(yi, vi, data=dat, method="EE", subset=random==0)
res.r0
res.r1 <- rma(yi, vi, data=dat, method="EE", subset=random==1)
res.r1
```
Here, we note that the two separate $Q$ statistics add up to the $Q_E$ statistic given earlier (i.e., $`r fmtx(res.r0$QE)` + `r fmtx(res.r1$QE)` = `r fmtx(res.anova$QE)`$).
Finally, the predicted/estimated effect for 'random = 0' and 'random = 1' can be obtained with:
```{r}
predict(res.anova, newmods=c(0,1))
```
All of the values shown in Exhibit 7.4 have now been computed.
***
## 5) Weighted Regression Analysis
The 'weighted regression analysis' described on pages 138-140 is again nothing else than a fixed-effects meta-regression model with moderators. Here, the authors are including both the `random` and the `intensity` moderators in the model. This can be done with:
```{r}
res.wr <- rma(yi, vi, mods = ~ random + intensity, data=dat, method="FE")
res.wr
```
Exhibit 7.6 (page 141) shows the same results.
***
## 6) Mixed-Effects Model
Finally, the mixed-effects meta-regression model is described on pages 140-142. These results can be obtained with:
```{r}
res.me <- rma(yi, vi, mods = ~ random + intensity, data=dat, method="DL")
res.me
```
These are the same results as shown in Exhibit 7.7 (page 141).
In essence, the model describes the (linear) relationship between the treatment intensity and the standardized mean differences and allows the intercept of the regression line to differ depending on whether random assignment was used within the studies or not. We can visualize the results as follows:
```{r, scatterplot, fig.width=8, fig.height=6.5, dev.args=list(pointsize=pointsize), fig.align='center'}
par(mar=c(5,4,2,2))
regplot(res.me, mod="intensity", xlim=c(1,7), pred=FALSE, ci=FALSE,
xlab="Intensity", ylab="Standardized Mean Difference",
pch=ifelse(dat$random == 0, 1, 19))
abline(a=coef(res.me)[1], b=coef(res.me)[3], lwd=2, lty="dotted")
abline(a=coef(res.me)[1] + coef(res.me)[2], b=coef(res.me)[3], lwd=2)
legend("topleft", inset=.02, lty=c("dotted", "solid"), lwd=2,
pch=c(21,19) ,legend=c("No Random Assignment", "Random Assignment"))
```
The points are drawn in size proportional to the weight the various studies received according to the model.
***
## 7) Heterogeneity Accounted For
The amount of heterogeneity accounted for by the moderator(s) included in the model is typically computed by comparing the amount of heterogeneity from the random-effects model with the amount of residual heterogeneity from the mixed-effects model (for more details, see the [here](https://www.metafor-project.org/doku.php/faq#for_mixed-effects_models_how_i)). The value is given in the output above. We can also carry out this computation explicitly by doing a full versus reduced model comparison with (note that a warning will be issued with respect to the likelihood ratio test that is also part of this output since LRTs should be based on ML/REML estimation):
```{r}
anova(res.re, res.me)
```
The value under `R^2` is the amount of variance (heterogeneity) accounted for. In this particular example, we estimate that `r fmtx(res.me$R2, digits=1)`% of the total amount of heterogeneity is accounted for by the two moderators. Note that Lipsey and Wilson calculate $R^2$ in a different way than is usually done in meta-analyses, so the value for $R^2$ shown in Exhibit 7.7 (page 141) is different from the value computed above.
***
## License
This documented is licensed under the following license: [CC Attribution-Noncommercial-Share Alike 4.0 International](http://creativecommons.org/licenses/by-nc-sa/4.0/).