library(tidyverse)
library(knitr)
library(broom)
library(skimr)
nhanes <- mice::nhanes2
nhanes
##      age  bmi  hyp chl
## 1  20-39   NA <NA>  NA
## 2  40-59 22.7   no 187
## 3  20-39   NA   no 187
## 4  60-99   NA <NA>  NA
## 5  20-39 20.4   no 113
## 6  60-99   NA <NA> 184
## 7  20-39 22.5   no 118
## 8  20-39 30.1   no 187
## 9  40-59 22.0   no 238
## 10 40-59   NA <NA>  NA
## 11 20-39   NA <NA>  NA
## 12 40-59   NA <NA>  NA
## 13 60-99 21.7   no 206
## 14 40-59 28.7  yes 204
## 15 20-39 29.6   no  NA
## 16 20-39   NA <NA>  NA
## 17 60-99 27.2  yes 284
## 18 40-59 26.3  yes 199
## 19 20-39 35.3   no 218
## 20 60-99 25.5  yes  NA
## 21 20-39   NA <NA>  NA
## 22 20-39 33.2   no 229
## 23 20-39 27.5   no 131
## 24 60-99 24.9   no  NA
## 25 40-59 27.4   no 186
  1. Explore missingness
library(nnet)
m <- multinom(age ~ bmi + hyp + chl, data = mice::nhanes2)
## # weights:  15 (8 variable)
## initial  value 14.281960 
## iter  10 value 4.305178
## iter  20 value 3.803287
## iter  30 value 3.487554
## iter  40 value 3.241432
## iter  50 value 2.782424
## iter  60 value 2.775229
## iter  70 value 2.653552
## iter  80 value 2.646572
## iter  90 value 2.574802
## iter 100 value 2.572907
## final  value 2.572907 
## stopped after 100 iterations
knitr::kable(tidy(m, conf.int = T),format = "markdown")
y.level term estimate std.error statistic p.value conf.low conf.high
40-59 (Intercept) 2.659563e+23 1.4471141 37.272541 0.0000000 1.559641e+22 4.535196e+24
40-59 bmi 3.640000e-04 5.6516387 -1.401096 0.1611854 0.000000e+00 2.353573e+01
40-59 hypyes 2.213278e+03 1.2615471 6.105384 0.0000000 1.867261e+02 2.623414e+04
40-59 chl 2.449175e+00 0.8347463 1.073082 0.2832344 4.769605e-01 1.257642e+01
60-99 (Intercept) 6.574010e+31 0.4127262 177.510570 0.0000000 2.927637e+31 1.476194e+32
60-99 bmi 9.320000e-05 5.7171151 -1.623397 0.1045046 0.000000e+00 6.849616e+00
60-99 hypyes 7.679312e+05 1.2615489 10.741918 0.0000000 6.478731e+04 9.102375e+06
60-99 chl 2.565845e+00 0.8357393 1.127490 0.2595353 4.987097e-01 1.320119e+01
  1. Complete-case analyis
complete <- nhanes %>% drop_na()
m <- multinom(age ~ bmi + hyp + chl, data = complete)
## # weights:  15 (8 variable)
## initial  value 14.281960 
## iter  10 value 4.305178
## iter  20 value 3.803287
## iter  30 value 3.487554
## iter  40 value 3.241432
## iter  50 value 2.782424
## iter  60 value 2.775229
## iter  70 value 2.653552
## iter  80 value 2.646572
## iter  90 value 2.574802
## iter 100 value 2.572907
## final  value 2.572907 
## stopped after 100 iterations
knitr::kable(tidy(m, conf.int = T),format = "markdown")
y.level term estimate std.error statistic p.value conf.low conf.high
40-59 (Intercept) 2.659563e+23 1.4471141 37.272541 0.0000000 1.559641e+22 4.535196e+24
40-59 bmi 3.640000e-04 5.6516387 -1.401096 0.1611854 0.000000e+00 2.353573e+01
40-59 hypyes 2.213278e+03 1.2615471 6.105384 0.0000000 1.867261e+02 2.623414e+04
40-59 chl 2.449175e+00 0.8347463 1.073082 0.2832344 4.769605e-01 1.257642e+01
60-99 (Intercept) 6.574010e+31 0.4127262 177.510570 0.0000000 2.927637e+31 1.476194e+32
60-99 bmi 9.320000e-05 5.7171151 -1.623397 0.1045046 0.000000e+00 6.849616e+00
60-99 hypyes 7.679312e+05 1.2615489 10.741918 0.0000000 6.478731e+04 9.102375e+06
60-99 chl 2.565845e+00 0.8357393 1.127490 0.2595353 4.987097e-01 1.320119e+01

  1. Single imputation

  2. Indicator variables