Repository navigation
Different results within mice compared to directly using method functions #645
Description
Describe the bug
I wanted to conduct a simulation study with one missing variable using mice.impute.lasso.select.norm. However, this method consistently failed to select the correct variables when used within the mice, despite working correctly when used directly. To address this issue, I created my own simplified wrapper function (called it short_mice), which should be sufficient for this situation. I then compared its performance to the mice wrapper through a short simulation using several methods (not just lasso.select.norm) and observed quite different results.
To Reproduce
Nsim = 10 and m = 5 should be increased!
rm(list = ls())
require(tidyverse)
require(mice)
require(mvtnorm)
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# Data generation ----
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
dgp <- function(n=50, p = 100, b_trt = 5, b_1 = 2, sig = 1.5, perc_missing = 0.3, miss_form = "~VOI*X1"){
X <- mvtnorm::rmvnorm(n= n, sigma = diag(p))
x1 <- X[,1]
colnames(X) <- paste0("X", 1:ncol(X))
noise <- rnorm(n=n)
VOI <- rep(c(0,1),ceiling(length(x1)/2) )[1:length(x1)]
y <- 0 + b_1*x1 + b_trt*VOI + noise*sig
data <- data.frame(y,VOI,X)
# Missing generation (MAR) only in the outcome y
m_mis <- glm(family = "binomial",
data = data,
formula = formula(paste0("I(y < quantile(data$y,perc_missing))",miss_form))
)
mi <- predict(m_mis, newdata = data, type = "response")# 1 = Missing
mi_bin <- rbinom(n = nrow(data), size = 1, prob = mi)
amp <- data
amp$y[mi_bin==1] <- NA
return(list("data" = data, "amp" = amp))
}
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# Shortcut mice ----
# can only impute one missing variable
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# data = data$amp
# m = 3
# method = "norm"
short_mice <- function(data, m, method,printFlag=FALSE, ...){
state <- list(it = 0, im = 0, dep = "", meth = "", log = FALSE)
loggedEvents <- data.frame(it = 0, im = 0, dep = "", meth = "", out = "")
y <- data$y
x <- as.matrix(data %>% dplyr::select(-y))
ry <- !is.na(y)
wy = NULL
mi <- lapply(0:m,function(i){
runif(1)
if(i==0){
mi <- data.frame(.imp = i, .id = 1:nrow(x) , y=y, x)
}else{
eval(parse(
text = paste0("ymis <- mice.impute.",method, "(y,ry,x,wy, ...)")
))
mi <- data.frame(.imp = i, .id = 1:nrow(x) , y=y, x)
mi$y[!ry] <- ymis
}
return(mi)
})
mi <- do.call("rbind", mi) %>% as.mids()
return(mi)
}
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# Run small simulation ----
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# run with higher numbers!
Nsim = 10 #100 ## INCREASE!!
m = 5 # 30 ## INCREASE!!
maxit = 1
analysis_form = "y~VOI"
method = c("lasso.select.norm", "norm")
wrapper = c("mice", "short_mice")
method_fact <- expand.grid(method = c("norm","lasso.select.norm","lasso.norm","rf")
,wrapper = c("mice","short_mice")
,stringsAsFactors = FALSE
)
set.seed(32325329)
res <- lapply(1:Nsim,function(s){
# print(s)
data <- dgp()
res <- lapply(1:nrow(method_fact), function(i){
eval(parse(
text = paste0("mice_m <- ",method_fact$wrapper[i]
,"(data = data$amp, m = m, maxit = maxit, method = '"
,method_fact$method[i]
,"',ls.meth = 'ridge', printFlag = FALSE)"
)
))
imp_info <- mice_m %>% with(data = ., lm(formula = formula(analysis_form) )) %>% pool()
mice_res <- imp_info %>% summary() %>%
subset(., term=="VOI", c("estimate","std.error")) %>%
data.frame()
return(mice_res)
})
res <- do.call("rbind", res)
# complete case
mod <- lm(data = data$amp, formula(analysis_form))
cc_res <- summary(mod)$coef["VOI",1:2]
cc_res <- data.frame(t(cc_res), "cc", NA)
names(cc_res) <- c("estimate","std.error", "method", "wrapper")
# imputing true model
mod <- lm(data = data$data, formula(analysis_form))
tc_res <- summary(mod)$coef["VOI",1:2]
tc_res <- data.frame(t(tc_res) )
names(tc_res) <- c("true_estimate","true_std.error")
# combine
res <- cbind(rbind(cbind(res,method_fact),cc_res), tc_res)
res$id <- s
# return
return(res)
})
res <- do.call("rbind", res)
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# Visualization ----
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
res %>%
ggplot(aes( y =paste0(method,"+",wrapper )
,x = (estimate -true_estimate)/true_estimate
,fill = method
) )+
geom_boxplot()+
geom_point()+
theme_bw()+
geom_vline(xintercept = 0)
Created on 2024-05-28 with reprex v2.1.0
Session info
sessioninfo::session_info()
#> ─ Session info ───────────────────────────────────────────────────────────────
#> setting value
#> version R version 4.3.2 (2023-10-31 ucrt)
#> os Windows 10 x64 (build 19045)
#> system x86_64, mingw32
#> ui RTerm
#> language (EN)
#> collate German_Switzerland.utf8
#> ctype German_Switzerland.utf8
#> tz Europe/Zurich
#> date 2024-05-28
#> pandoc 3.1.1 @ C:/Program Files/RStudio/resources/app/bin/quarto/bin/tools/ (via rmarkdown)
#>
#> ─ Packages ───────────────────────────────────────────────────────────────────
#> package * version date (UTC) lib source
#> backports 1.4.1 2021-12-13 [1] CRAN (R 4.3.1)
#> boot 1.3-28.1 2022-11-22 [2] CRAN (R 4.3.2)
#> broom 1.0.5 2023-06-09 [1] CRAN (R 4.3.2)
#> cli 3.6.1 2023-03-23 [1] CRAN (R 4.3.2)
#> codetools 0.2-19 2023-02-01 [2] CRAN (R 4.3.2)
#> colorspace 2.1-0 2023-01-23 [1] CRAN (R 4.3.2)
#> curl 5.1.0 2023-10-02 [1] CRAN (R 4.3.2)
#> digest 0.6.35 2024-03-11 [1] CRAN (R 4.3.3)
#> dplyr * 1.1.3 2023-09-03 [1] CRAN (R 4.3.2)
#> evaluate 0.23 2023-11-01 [1] CRAN (R 4.3.2)
#> fansi 1.0.6 2023-12-08 [1] CRAN (R 4.3.3)
#> farver 2.1.1 2022-07-06 [1] CRAN (R 4.3.2)
#> fastmap 1.1.1 2023-02-24 [1] CRAN (R 4.3.2)
#> forcats * 1.0.0 2023-01-29 [1] CRAN (R 4.3.2)
#> foreach 1.5.2 2022-02-02 [1] CRAN (R 4.3.2)
#> fs 1.6.3 2023-07-20 [1] CRAN (R 4.3.2)
#> generics 0.1.3 2022-07-05 [1] CRAN (R 4.3.2)
#> ggplot2 * 3.5.0 2024-02-23 [1] CRAN (R 4.3.3)
#> glmnet 4.1-8 2023-08-22 [1] CRAN (R 4.3.2)
#> glue 1.7.0 2024-01-09 [1] CRAN (R 4.3.3)
#> gtable 0.3.4 2023-08-21 [1] CRAN (R 4.3.2)
#> highr 0.10 2022-12-22 [1] CRAN (R 4.3.2)
#> hms 1.1.3 2023-03-21 [1] CRAN (R 4.3.2)
#> htmltools 0.5.7 2023-11-03 [1] CRAN (R 4.3.2)
#> iterators 1.0.14 2022-02-05 [1] CRAN (R 4.3.2)
#> jomo 2.7-6 2023-04-15 [1] CRAN (R 4.3.2)
#> knitr 1.45 2023-10-30 [1] CRAN (R 4.3.2)
#> labeling 0.4.3 2023-08-29 [1] CRAN (R 4.3.1)
#> lattice 0.21-9 2023-10-01 [2] CRAN (R 4.3.2)
#> lifecycle 1.0.4 2023-11-07 [1] CRAN (R 4.3.3)
#> lme4 1.1-35.1 2023-11-05 [1] CRAN (R 4.3.2)
#> lubridate * 1.9.3 2023-09-27 [1] CRAN (R 4.3.2)
#> magrittr 2.0.3 2022-03-30 [1] CRAN (R 4.3.2)
#> MASS 7.3-60 2023-05-04 [2] CRAN (R 4.3.2)
#> Matrix 1.6-1.1 2023-09-18 [2] CRAN (R 4.3.2)
#> mice * 3.16.0 2023-06-05 [1] CRAN (R 4.3.2)
#> minqa 1.2.6 2023-09-11 [1] CRAN (R 4.3.2)
#> mitml 0.4-5 2023-03-08 [1] CRAN (R 4.3.2)
#> munsell 0.5.0 2018-06-12 [1] CRAN (R 4.3.2)
#> mvtnorm * 1.2-3 2023-08-25 [1] CRAN (R 4.3.2)
#> nlme 3.1-163 2023-08-09 [2] CRAN (R 4.3.2)
#> nloptr 2.0.3 2022-05-26 [1] CRAN (R 4.3.2)
#> nnet 7.3-19 2023-05-03 [2] CRAN (R 4.3.2)
#> pan 1.9 2023-08-21 [1] CRAN (R 4.3.2)
#> pillar 1.9.0 2023-03-22 [1] CRAN (R 4.3.2)
#> pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 4.3.2)
#> purrr * 1.0.2 2023-08-10 [1] CRAN (R 4.3.2)
#> R.cache 0.16.0 2022-07-21 [1] CRAN (R 4.3.3)
#> R.methodsS3 1.8.2 2022-06-13 [1] CRAN (R 4.3.3)
#> R.oo 1.26.0 2024-01-24 [1] CRAN (R 4.3.3)
#> R.utils 2.12.3 2023-11-18 [1] CRAN (R 4.3.3)
#> R6 2.5.1 2021-08-19 [1] CRAN (R 4.3.2)
#> ranger 0.16.0 2023-11-12 [1] CRAN (R 4.3.2)
#> Rcpp 1.0.12 2024-01-09 [1] CRAN (R 4.3.3)
#> readr * 2.1.4 2023-02-10 [1] CRAN (R 4.3.2)
#> reprex 2.1.0 2024-01-11 [1] CRAN (R 4.3.3)
#> rlang 1.1.2 2023-11-04 [1] CRAN (R 4.3.2)
#> rmarkdown 2.25 2023-09-18 [1] CRAN (R 4.3.2)
#> rpart 4.1.21 2023-10-09 [2] CRAN (R 4.3.2)
#> rstudioapi 0.15.0 2023-07-07 [1] CRAN (R 4.3.2)
#> scales 1.3.0 2023-11-28 [1] CRAN (R 4.3.3)
#> sessioninfo 1.2.2 2021-12-06 [1] CRAN (R 4.3.2)
#> shape 1.4.6.1 2024-02-23 [1] CRAN (R 4.3.2)
#> stringi 1.8.3 2023-12-11 [1] CRAN (R 4.3.2)
#> stringr * 1.5.0 2022-12-02 [1] CRAN (R 4.3.2)
#> styler 1.10.3 2024-04-07 [1] CRAN (R 4.3.3)
#> survival 3.5-7 2023-08-14 [2] CRAN (R 4.3.2)
#> tibble * 3.2.1 2023-03-20 [1] CRAN (R 4.3.2)
#> tidyr * 1.3.0 2023-01-24 [1] CRAN (R 4.3.2)
#> tidyselect 1.2.1 2024-03-11 [1] CRAN (R 4.3.3)
#> tidyverse * 2.0.0 2023-02-22 [1] CRAN (R 4.3.2)
#> timechange 0.2.0 2023-01-11 [1] CRAN (R 4.3.2)
#> tzdb 0.4.0 2023-05-12 [1] CRAN (R 4.3.2)
#> utf8 1.2.4 2023-10-22 [1] CRAN (R 4.3.2)
#> vctrs 0.6.5 2023-12-01 [1] CRAN (R 4.3.3)
#> withr 3.0.0 2024-01-16 [1] CRAN (R 4.3.3)
#> xfun 0.41 2023-11-01 [1] CRAN (R 4.3.2)
#> xml2 1.3.5 2023-07-06 [1] CRAN (R 4.3.2)
#> yaml 2.3.7 2023-01-23 [1] CRAN (R 4.3.2)
#>
#> [1] C:/Users/SepinJ/AppData/Local/R/win-library/4.3
#> [2] C:/Program Files/R/R-4.3.2/library
#>
#> ──────────────────────────────────────────────────────────────────────────────Expected behavior
I would expect the outputs to be more similar.
Thank you very much for any clarification!
Activity
EdoardoCostantini commentedon May 31, 2024
The difference comes from the pre-processing done by the
mice()function. In particular, theremove.lindep()function is run by default in mice to try to stabilize the imputation models. When trying to impute the data you are generating, many variables are dropped out of the imputation models before themice.impute.univariate imputation methods are called. You can check this by running the following:In the logged events, you should see that many variables have been "kicked out," meaning they never made it to the
mice.impute.lasso.select.norm()function. With the wrapper you wrote, the linear dependency checks are not performed, and the lasso is the only way you handle variable selection. To avoid this pre-processing inmice()you can set theepsargument in themice()call to 0. In your code, change the eval statement to:eval(parse( text = paste0( "mice_m <- ", method_fact$wrapper[i], "(data = data$amp, m = m, maxit = maxit, method = '", method_fact$method[i], "',ls.meth = 'ridge', printFlag = FALSE, eps = 0)" ) ))I ran your code with this change, and the resulting plot is the following:
There are no more differences between your wrapper and the standard
mice().I think I should clarify that
eps = 0is recommended when running lasso-based imputation methods. I will make a pull request to update the documentation soon.Let me know if this solves your issue entirely or if there is anything else I can help with.
stefvanbuuren commentedon Aug 16, 2024
@EdoardoCostantini thanks. No really a bug or problem. I will convert the issue to a Discussion.