Pooling Poisson GLM with Robust Variance after MICE imputation

[felippelazar]

Hi Stef,

I can't find a solution for running the poisson GLM with robust variance in mice imputace data-sets and pooling the results.

In a pipeline without imputed data-sets, the workaround in R for running robust variance is using the sandwich package and 'coeftest' from lmtest. Here is how I estimated std. errors. with robust variance in R (using a fictious example):

library(sandwich)
library(lmtest)

model.diabetes <- glm(diabetes ~ age + sex, data=df, family='poisson')
coeftest(model.diabetes, vcov. = sandwich)

The above code would then display estimates, standard errors, z values and p-values for poisson with robust variance.

I know the problem would be mostly related to the glm itself and not mice, but I was wondering if you have any workaround for solving this problem.

Thanks in advance,

Felippe

[stefvanbuuren]

Felippe,

A reprex would help.

• Anything that you can do with complete data, you can also do with imputed data (after extraction with complete()), so "running" should not present problems
• For "pooling", mice() relies on the broom package. For glm there are appropriate tidy() and glance() functions. Find out if such functions also exist for lmtest::coeftest(), and if not, you may need to create these.

Hope this helps, Stef.

[felippelazar]

Hi Stef,

Thanks for the quick reply. Here goes a reproducible code:

library(lmtest)
#> Loading required package: zoo
#> 
#> Attaching package: 'zoo'
#> The following objects are masked from 'package:base':
#> 
#>     as.Date, as.Date.numeric
library(mice)
#> 
#> Attaching package: 'mice'
#> The following objects are masked from 'package:base':
#> 
#>     cbind, rbind
library(sandwich)
library(broom)

# Importing Data-Frame and Changing to 0-1 Hypertension
nhanes2$hyp <- ifelse(nhanes2$hyp == 'no', 0, 1)

# Imputing MICE
imp <- mice(nhanes2, seed = 23109)
#> 
#>  iter imp variable
#>   1   1  bmi  hyp  chl
#>   1   2  bmi  hyp  chl
#>   1   3  bmi  hyp  chl
#>   1   4  bmi  hyp  chl
#>   1   5  bmi  hyp  chl
#>   2   1  bmi  hyp  chl
#>   2   2  bmi  hyp  chl
#>   2   3  bmi  hyp  chl
#>   2   4  bmi  hyp  chl
#>   2   5  bmi  hyp  chl
#>   3   1  bmi  hyp  chl
#>   3   2  bmi  hyp  chl
#>   3   3  bmi  hyp  chl
#>   3   4  bmi  hyp  chl
#>   3   5  bmi  hyp  chl
#>   4   1  bmi  hyp  chl
#>   4   2  bmi  hyp  chl
#>   4   3  bmi  hyp  chl
#>   4   4  bmi  hyp  chl
#>   4   5  bmi  hyp  chl
#>   5   1  bmi  hyp  chl
#>   5   2  bmi  hyp  chl
#>   5   3  bmi  hyp  chl
#>   5   4  bmi  hyp  chl
#>   5   5  bmi  hyp  chl

# Running Multiple Poisson Robust Regressions
fit <- with(imp, coeftest(glm(hyp ~ age + bmi, family='poisson'), vcov. = sandwich))
#> Warning in sqrt(diag(se)): NaNs produced

# Pooling the Results
print(pool(fit))
#> Error: No glance method for objects of class coeftest

# However, there's a tidy method for this object
coef.object <- coeftest(glm(hyp ~ age + bmi, family='poisson', data=nhanes2), vcov. = sandwich)
broom::tidy(coef.object)
#> # A tibble: 4 x 5
#>   term        estimate std.error statistic  p.value
#>   <chr>          <dbl>     <dbl>     <dbl>    <dbl>
#> 1 (Intercept)  -36.6       5.45      -6.72 1.76e-11
#> 2 age40-59      22.6       1.12      20.1  6.04e-90
#> 3 age60-99      23.4       1.39      16.8  2.08e-63
#> 4 bmi            0.490     0.156      3.15 1.64e- 3

As you mentioned, there is already a tidy() function for coeftest objects, however, I couldn't find a glance function.

My doubt is: how do I add this created functions to the mice pipeline for calculating pooled results?

Thanks again,

Felippe

[felippelazar]

Just to follow-up, if the only information the glance function will take is the df.residual (as mentioned), we could probably use the glance function on the glm object and the tidy function on the coeftest.

Heres what I imagined:

library(lmtest)
#> Loading required package: zoo
#> 
#> Attaching package: 'zoo'
#> The following objects are masked from 'package:base':
#> 
#>     as.Date, as.Date.numeric
library(mice)
#> 
#> Attaching package: 'mice'
#> The following objects are masked from 'package:base':
#> 
#>     cbind, rbind
library(sandwich)
library(broom)

# Importing Data-Frame and Changing to 0-1 Hypertension
nhanes2$hyp <- ifelse(nhanes2$hyp == 'no', 0, 1)

glm.poisson <- glm(hyp ~ age + bmi, family='poisson', data=nhanes2)

tidy(coeftest(glm.poisson, vcov. = sandwich))
#> # A tibble: 4 x 5
#>   term        estimate std.error statistic  p.value
#>   <chr>          <dbl>     <dbl>     <dbl>    <dbl>
#> 1 (Intercept)  -36.6       5.45      -6.72 1.76e-11
#> 2 age40-59      22.6       1.12      20.1  6.04e-90
#> 3 age60-99      23.4       1.39      16.8  2.08e-63
#> 4 bmi            0.490     0.156      3.15 1.64e- 3

glance(glm.poisson)
#> # A tibble: 1 x 7
#>   null.deviance df.null logLik   AIC   BIC deviance df.residual
#>           <dbl>   <int>  <dbl> <dbl> <dbl>    <dbl>       <int>
#> 1          11.1      15  -5.75  19.5  22.6     3.51          12

Does it make sense? My problem would still be how to incorporate this alternate method in the pool function pipeline.

Thanks again,

Felippe

[stefvanbuuren]

Fixed in mice 3.8.5.

I have adapted the pool() function. If dfcom is specified, then glance() is not tried, otherwise glance() is tried, followed by two other methods when there is no glance().

The following runs in mice 3.8.5

library(lmtest)
library(mice)
library(sandwich)
library(broom)

# Importing Data-Frame and Changing to 0-1 Hypertension
nhanes2$hyp <- ifelse(nhanes2$hyp == 'no', 0, 1)

# Imputing MICE
imp <- mice(nhanes2, seed = 23109, print = FALSE)

# Running Multiple Poisson Robust Regressions
fit <- with(imp, coeftest(glm(hyp ~ age + bmi, family='poisson'), vcov. = sandwich))

# investigate availability of tidy() and glance()
summary(fit, "tidy")
summary(fit, "glance")

# Pooling the Results ()
print(pool(fit))

# Check: The inf is incorrect here, but gotten from
df.residual(getfit(fit, 1))

# Set correct df
print(pool(fit, dfcom = 21))

Note that detection may fail occasionally, as above. In that case, specify dfcom manually.

This tweak should also fix #142, #148, #190 and #199.

[felippelazar]

Thanks Stef,

Cheers,

Felippe