-
Notifications
You must be signed in to change notification settings - Fork 1
/
multipleImputation.Rmd
461 lines (343 loc) · 11.6 KB
/
multipleImputation.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
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
---
title: "Multiple Imputation"
output:
html_document:
code_folding: show
---
```{r setup, include = FALSE}
knitr::opts_chunk$set(
echo = TRUE,
error = TRUE,
comment = "",
class.source = "fold-show")
```
# Preamble
## Install Libraries
```{r}
#install.packages("remotes")
#remotes::install_github("DevPsyLab/petersenlab")
```
## Load Libraries
```{r, message = FALSE, warning = FALSE}
library("tidyverse")
library("psych")
library("mice")
library("micemd")
library("miceadds")
library("mitml")
library("Amelia")
library("jomo")
library("parallel")
library("future")
library("furrr")
library("nlme")
library("MplusAutomation")
```
## Data
```{r}
data(Oxboys, package = "nlme")
Oxboys_addNA <- data.frame(Oxboys)
Oxboys_addNA <- Oxboys_addNA %>%
rename(id = Subject)
Oxboys_addNA$id <- as.integer(Oxboys_addNA$id)
Oxboys_addNA$Occasion <- as.integer(Oxboys_addNA$Occasion)
set.seed(52242)
Oxboys_addNA[sample(1:nrow(Oxboys_addNA), 25), "height"] <- NA
dataToImpute <- Oxboys_addNA
varsToImpute <- c("height")
Y <- varsToImpute
```
# Types of Missingness
1. Missing Completely at Random (MCAR)
- the probability of being missing is the same for all cases
- missingness is not related to variables in the model
1. (Conditionally) Missing at Random (MAR)
- missingness is related to variables in the model, but once we condition on the variables in the model, the missingness is haphazard
- the unobserved values do not play a role
1. Missing Not at Random (MNAR)
- missingness is related to variables that are not in the model (i.e., for reasons that are unknown)
- unobserved variables carry information about whether a case will have missing data
# Approaches to Handle Missing Data
## For MCAR/MAR Missingness
1. Full Information Maximum Likelihood (FIML)
1. Multiple Imputation
## For MNAR Missingness
https://stefvanbuuren.name/fimd/sec-nonignorable.html (archived at https://perma.cc/N7WW-HDZF)
https://cran.r-project.org/web/packages/missingHE/vignettes/Fitting_MNAR_models_in_missingHE.html (archived at https://perma.cc/9X25-5D8G)
1. find more data about the causes for the missingness
1. sensitivity analyses (what-if analyses) to see how sensitive the results are under various scenarios
- https://stefvanbuuren.name/fimd/sec-MCAR.html (archived at https://perma.cc/9KM9-3NX3)
1. selection models
- simultaneously estimate the focal model and a missingness model, where the missingness model has the missing data indicator as a dependent variable, as predicted by the original dependent variable, the original predictors, and any additional covariates etc.
- if the missingness model is approximately correct, the focal model adjusts in way that removes nonresponse bias
- similar to a mediation process
- X → Y
- Y → missingness
- X → missingness
- https://quantitudepod.org/s4e08-craig-enders/ (archived at https://perma.cc/FY9L-L3F7)
1. pattern-mixture models
- missing data indicator is a predictor variable
- similar to a moderation process; subgroups of cases have different parameter estimates
# Approaches to Multiple Imputation
- multiple imputation by joint modeling
- assumes that the variables in follow a joint distribution, e.g.:
- multivariate normal (i.e., multivariate normal imputation)
- `R` packages:
- `jomo`
- `pan`
- `Amelia`
- `Mplus`
- multiple imputation by chained equations
- aka:
- fully conditional specification multiple imputation
- sequential regression multiple imputation
- `R` packages:
- `mice`
- predictive mean matching (`?mice::mice.impute.pmm`) can be useful to obtain bounded imputations for non-normally distributed variables
# Describe Missing Data
```{r}
describe(dataToImpute)
md.pattern(dataToImpute, rotate.names = TRUE)
```
# Pre-Imputation Setup
## Specify Number of Imputations
```{r}
numImputations <- 5 # generally use 100 imputations; this example uses 5 for speed
```
## Detect Cores
```{r}
numCores <- parallel::detectCores() - 1
```
# Multilevel Multiple Imputation
## Three-Level and Cross-Classified Data
https://simongrund1.github.io/posts/multiple-imputation-for-three-level-and-cross-classified-data/ (archived at https://perma.cc/N4PP-A3V6)
## Methods
### `mice` {#mice}
#### Types
##### Continuous Outcomes
https://stefvanbuuren.name/fimd/sec-multioutcome.html#methods (archived at https://perma.cc/8CDA-TS3K)
```{r, eval = FALSE}
?mice::mice.impute.2l.norm
?mice::mice.impute.2l.pan
?mice::mice.impute.2l.lmer
?miceadds::mice.impute.2l.pmm
?miceadds::mice.impute.2l.contextual.pmm
?miceadds::mice.impute.2l.continuous
?micemd::mice.impute.2l.2stage.norm
?micemd::mice.impute.2l.2stage.pmm
?micemd::mice.impute.2l.glm.norm
?micemd::mice.impute.2l.jomo
```
##### Binary Outcomes
https://stefvanbuuren.name/fimd/sec-catoutcome.html#methods-1 (archived at https://perma.cc/5QHF-YRP6)
```{r, eval = FALSE}
?mice::mice.impute.2l.bin
?miceadds::mice.impute.2l.binary
?miceadds::mice.impute.2l.pmm
?miceadds::mice.impute.2l.contextual.pmm
?micemd::mice.impute.2l.2stage.bin
?micemd::mice.impute.2l.glm.bin
```
##### Ordinal Outcomes
https://stefvanbuuren.name/fimd/sec-multioutcome.html#methods (archived at https://perma.cc/8CDA-TS3K)
```{r, eval = FALSE}
?miceadds::mice.impute.2l.pmm
?miceadds::mice.impute.2l.contextual.pmm
?micemd::mice.impute.2l.2stage.pmm
```
##### Count Outcomes
https://stefvanbuuren.name/fimd/sec-catoutcome.html#methods-1 (archived at https://perma.cc/5QHF-YRP6)
```{r, eval = FALSE}
?micemd::mice.impute.2l.2stage.pois
?micemd::mice.impute.2l.glm.pois
?countimp::mice.impute.2l.poisson
?countimp::mice.impute.2l.nb2
?countimp::mice.impute.2l.zihnb
```
#### Specifying the Imputation Method
```{r}
meth <- make.method(dataToImpute)
meth[1:length(meth)] <- ""
meth[Y] <- "2l.pmm" # specify the imputation method here; this can differ by outcome variable
```
#### Specifying the Predictor Matrix
A predictor matrix is a matrix of values, where:
- columns with non-zero values are predictors of the variable specified in the given row
- the diagonal of the predictor matrix should be zero because a variable cannot predict itself
The values are:
- NOT a predictor of the outcome: `0`
- cluster variable: `-2`
- fixed effect of predictor: `1`
- fixed effect and random effect of predictor: `2`
- include cluster mean of predictor in addition to fixed effect of predictor: `3`
- include cluster mean of predictor in addition to fixed effect and random effect of predictor: `4`
```{r}
pred <- make.predictorMatrix(dataToImpute)
pred[1:nrow(pred), 1:ncol(pred)] <- 0
pred[Y, "id"] <- (-2) # cluster variable
pred[Y, "Occasion"] <- 1 # fixed effect predictor
pred[Y, "age"] <- 2 # random effect predictor
pred[Y, Y] <- 1 # fixed effect predictor
diag(pred) <- 0
pred
```
#### Syntax
```{r}
mi_mice <- mice(
as.data.frame(dataToImpute),
method = meth,
predictorMatrix = pred,
m = numImputations,
maxit = 5, # generally use 100 maximum iterations; this example uses 5 for speed
seed = 52242)
```
### `jomo` {#jomo}
```{r}
level1Vars <- c("height")
level2Vars <- c("v3","v4")
clusterVars <- c("id")
fullyObservedCovariates <- c("age","Occasion")
set.seed(52242)
mi_jomo <- jomo(
Y = data.frame(dataToImpute[, level1Vars]),
#Y2 = data.frame(dataToImpute[, level2Vars]),
X = data.frame(dataToImpute[, fullyObservedCovariates]),
clus = data.frame(dataToImpute[, clusterVars]),
nimp = numImputations,
meth = "random"
)
```
### `Amelia` {#amelia}
- `!` in the console output indicates that the current estimated complete data covariance matrix is not invertible
- `*` in the console output indicates that the likelihood has not monotonically increased in that step
```{r}
boundVars <- c("height")
boundCols <- which(names(dataToImpute) %in% boundVars)
boundLower <- 100
boundUpper <- 200
varBounds <- cbind(boundCols, boundLower, boundUpper)
set.seed(52242)
mi_amelia <- amelia(
dataToImpute,
m = numImputations,
ts = "age",
cs = "id",
polytime = 1,
intercs = TRUE,
#ords = ordinalVars,
#bounds = varBounds,
parallel = "snow",
#ncpus = numCores,
empri = .01*nrow(dataToImpute)) # ridge prior for numerical stability
```
### `Mplus` {#mplus}
#### Save `Mplus` Data File
Save `R` object as `Mplus` data file:
```{r}
prepareMplusData(dataToImpute, file.path("dataToImpute.dat"))
```
#### `Mplus` Syntax for Multilevel Imputation
`Mplus` syntax for multilevel imputation:
```Mplus
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!! MPLUS SYNTAX LINES CANNOT EXCEED 90 CHARACTERS;
!!!!! VARIABLE NAMES AND PARAMETER LABELS CANNOT EXCEED 8 CHARACTERS EACH;
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
TITLE: Model Title
DATA: FILE = "dataToImpute.dat";
VARIABLE:
NAMES = id age height Occasion;
MISSING = .;
USEVARIABLES ARE age height Occasion;
!CATEGORICAL ARE INSERT_NAMES_OF_CATEGORICAL_VARIABLES_HERE;
CLUSTER = id;
ANALYSIS:
TYPE = twolevel basic;
bseed = 52242;
PROCESSORS = 2;
DATA IMPUTATION:
IMPUTE = age(0-100) height Occasion; !put ' (c)' after categorical vars
NDATASETS = 100;
SAVE = imp*.dat
```
Putting a range of values after a variable (e.g., `0-100`) sets the lower and upper bounds of the imputed values.
This would save a `implist.dat` file that can be used to run the model on the multiply imputed data, as shown [here](mplus.html#multipleImputation).
# Parallel Processing
https://www.gerkovink.com/miceVignettes/futuremice/Vignette_futuremice.html (archived at https://perma.cc/4SNE-RCSR)
```{r}
mi_parallel_mice <- futuremice(
dataToImpute,
method = meth,
predictorMatrix = pred,
m = numImputations,
maxit = 5, # generally use 100 maximum iterations; this example uses 5 for speed
parallelseed = 52242,
n.core = numCores,
packages = "miceadds")
```
# Imputation Diagnostics
## Logged Events
```{r}
mi_mice$loggedEvents
```
## Trace Plots
On convergence, the streams of the trace plots should intermingle and be free of any trend (at the later iterations).
Convergence is diagnosed when the variance between different sequences is no larger than the variance within each individual sequence.
- https://stefvanbuuren.name/fimd/sec-algoptions.html (archived at https://perma.cc/4S54-465R)
```{r}
plot(mi_mice, c("height"))
```
## Density Plots
```{r}
densityplot(mi_mice)
densityplot(mi_mice, ~ height)
```
## Strip Plots
```{r}
stripplot(mi_mice)
stripplot(mi_mice, height ~ .imp)
```
# Post-Processing
## Modify/Create New Variables
### `mice`
```{r}
mi_mice_long <- complete(
mi_mice,
action = "long",
include = TRUE)
mi_mice_long$newVar <- mi_mice_long$age * mi_mice_long$height
```
### `jomo`
```{r}
mi_jomo <- mi_jomo %>%
rename(height = dataToImpute...level1Vars.)
mi_jomo$newVar <- mi_jomo$age * mi_jomo$height
```
### `Amelia`
```{r}
for(i in 1:length(mi_amelia$imputations)){
mi_amelia$imputations[[i]]$newVar <- mi_amelia$imputations[[i]]$age * mi_amelia$imputations[[i]]$height
}
```
## Convert to `mids` object
### `mice`
```{r}
mi_mice_mids <- as.mids(mi_mice_long)
```
### `jomo`
```{r}
mi_jomo_mids <- miceadds::jomo2mids(mi_jomo)
```
### `Amelia`
```{r}
mi_amelia_mids <- miceadds::datlist2mids(mi_amelia$imputations)
```
## Export for `Mplus`
```{r, eval = FALSE}
mids2mplus(mi_mice_mids, path = file.path("InsertFilePathHere", fsep = ""))
mids2mplus(mi_jomo_mids, path = file.path("InsertFilePathHere", fsep = ""))
mids2mplus(mi_amelia_mids, path = file.path("InsertFilePathHere", fsep = ""))
```
# Resources
https://stefvanbuuren.name/fimd (archived at https://perma.cc/46U2-QTM6)
https://www.gerkovink.com/miceVignettes/Multi_level/Multi_level_data.html (archived at https://perma.cc/SF32-D7ZF)