Multiple imputation MICE categorical variable: problem with pooling after running ANOVA

Question

I am trying to answer a question about satisfaction and its relation with a certain variable (numeric, 1-10). However, my data contains a lot of missing values in the satisfaction outcome, therefore I applied multiple imputation (all other variables are complete, no missings).

So far, I managed to create a MI with 20 datasets. I ran the ANOVA on all imputed datasets which gives me a Mira object (list of 4) that contains the analyses of all 20 sets. However, when I try to pool them I get the following error:

Error in summarize(): ! Problem while computing qbar = mean(.data$estimate). ℹ The error occurred in group 1: term = satisfaction Caused by error in .data$estimate: ! Column estimate not found in .data.

It looks like it needs some sort of mean. However, I am new to R and do not quite understand what I am doing wrong. There indeed is no mean or anything alike from the satisfaction in my Mira object, see screenshot below. Does anybody know what I am doing wrong? Thank you for your response!

Answer 1

The problem is that you need to be pooling the underlying lm models, not the ANOVA summaries of the models. The statistical theory and programs for multiple imputation are built around regression modeling, so we need to use that structure. ANOVA is equivalent to a linear regression using a single categorical variable, so the resulting analysis is equivalent to what you want to do.

The output won't be the same as what you're probably used to; the usual sum of squares decomposition doesn't really make sense after multiple imputation. Instead you have a couple options for doing a hypothesis test/summary. The first option is to go all-in on the regression modeling approach and look at each group/category separately. This method will define one of the groups as the "reference group" and the p-values will correspond to hypothesis tests that the average in that group is the same as the average in the reference group. In this example data reg has 5 levels (north, east, west, south, city) and north is the reference level. The p-values below show whether each of the other 4 regions has an average height that is statistically different than the north region. You can change the reference level to whatever you want.

library(mice)
# Set seed for reproducibility
set.seed(1234)
# Use built-in dataset
dat <- mice::boys
# Make imputations
imps <- mice(dat)
# Fit lm models
fits <- with(imps, lm(hgt ~ reg))
# Pool models
pooled <- pool(fits)
# Get results
summary(pooled)
#>          term  estimate std.error statistic       df      p.value
#> 1 (Intercept) 149.62394  5.158947 29.002809 689.4577 0.0000000000
#> 2     regeast -16.26098  6.294533 -2.583350 728.6964 0.0099782919
#> 3     regwest -20.68640  5.954958 -3.473811 710.4351 0.0005442446
#> 4    regsouth -22.92516  6.171241 -3.714838 667.8366 0.0002203387
#> 5     regcity -25.36682  7.461717 -3.399596 721.9398 0.0007119482

The above procedure outputs many different p-values, each corresponding to a pair-wise comparison. If you want to perform an "omnibus" test of whether the entire set of reg variables is significant overall, you can use a "D1", "D2", "D3" test. If you want details see this textbook by the mice creator, https://stefvanbuuren.name/fimd/sec-multiparameter.html

The author recommends the D1 test in general. In the case when you only have a single independent variable in the lm() models it's quite easy to perform this test.

D1(fits)
#>    test statistic df1      df2 dfcom    p.value         riv
#>  1 ~~ 2  4.311034   4 736.3446   743 0.00187468 0.008875879

This p-value corresponds to the test that the intercept-only model is equally good as the model with reg. The p-value is small, so we conclude that reg significantly affects the hgt variable.

If you're performing an ANCOVA with other variables in the model it's a little more complicated to perform the D1 test, but still quite manageable.