forked from psych251/class_examples
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathmeta-analysis_exercise.Rmd
131 lines (94 loc) · 3.12 KB
/
meta-analysis_exercise.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
---
title: "Meta-Analysis Demo"
author: "Psych 251"
date: "10/27/2020"
output: html_document
---
# Preliminaries
Packages needed for the class example.
```{r packages}
library(tidyverse)
library(metafor)
library(here)
```
You'd need to run this next chunk if you wanted to re-download the data from MetaLab, which would also mean installing the `metalabr` package (commented out).
```{r data, eval=FALSE}
#devtools::install_github("langcog/metalabr")
library(metalabr)
ml_dataset_info <- get_metalab_dataset_info()
d <- get_metalab_data(filter(ml_dataset_info, short_name == "mutex"))
d %>%
select(long_cite, short_cite, expt_num, n, d_calc, d_var_calc,
mean_age_months, year) %>%
filter(mean_age_months <= 24) %>% # subset for ease
write_csv(here("data/mutex.csv"))
```
Load in the pre-cached data.
```{r dataloading}
d <- read_csv(here("data/mutex.csv"))
```
# Basic descriptives
Always take a look at the data first.
```{r examine}
d
```
The effect sizes are in `d_calc` (for Cohen's $d$, calculated from the data).
```{r hist}
ggplot(d, aes(x = d_calc)) +
geom_histogram(binwidth = .1) +
xlab("Effect Size (d)")
```
Since these are developmental data, we can plot them against age.
```{r ageplot}
ggplot(d, aes(x = mean_age_months, y = d_calc)) +
geom_point(aes(size = n), alpha = .5) +
geom_smooth(method = "lm") +
xlab("Mean age (months)") +
ylab("Effect size (d)") +
geom_hline(lty = 2, yintercept = 0)
```
# Meta-analysis
Random effects meta-analysis - this is the default for `metafor`. The main command for metafor is `rma` - that's like the `lm` (linear model) or `lmer` (linear mixed effects model of the `lme4` package).
```{r ranef}
random_effects_mod <- rma(yi = d_calc, vi = d_var_calc,
slab = short_cite, data = d)
summary(random_effects_mod)
```
Investigate heterogeneity if we remove the outlier.
```{r}
d_post_hoc_pruned <- filter(d, d_calc < 1 & d_calc > 0)
random_effects_posthoc_mod <- rma(yi = d_calc, vi = d_var_calc,
slab = short_cite, data = d_post_hoc_pruned)
summary(random_effects_posthoc_mod)
```
For kicks, try fixed effects.
```{r fixef}
fixed_effects_mod <- rma(yi = d_calc, vi = d_var_calc,
slab = short_cite, data = d, method = "FE")
summary(fixed_effects_mod)
```
# Forest plot
`metafor` also lets you create forest plots.
```{r forest-ranef}
forest(random_effects_mod)
```
Compare to the forest plot for the fixed effects model.
```{r forest-fixef}
forest(fixed_effects_mod)
```
# Funnel plot
A funnel plot can be used to diagnose publication bias (though it's not the most sensitive way to do so).
```{r funnel}
funnel(random_effects_mod)
```
# Meta-regression
Meta-regression asks whether study-level covariates (like say year of publication or average age of kids) are related to the effect size.
```{r metareg-age}
meta_reg_model <- rma(yi = d_calc, vi = d_var_calc,
mods = ~ mean_age_months,
slab = short_cite, data = d)
summary(meta_reg_model)
```
Try asking if publication `year` is a significant meta-regressor.
```{r metareg-year}
```