search

amices / mice

Multivariate Imputation by Chained Equations
https://amices.org/mice/
GNU General Public License v2.0
447 stars 108 forks source link

Error in eval(expr, p): object 'xxx' not found #677

Closed benebeer closed 1 week ago

benebeer commented 2 weeks ago

Describe the bug The variable "intervar" is supposed to work as a dynamic variable in the negative binomial regression. It represents the interaction variables of interest.

intervar can be called within a caption: cat(paste0("\n**---Interaction analysis ", intervar, " ---**\n"))

However, it cannot be found within the environment of the with-function:

      binomial_inter <- with(reprexdf_imputed_weighted, {
        glm.nb(as.formula(
          paste(
            "sos_out_countdeathcvhosphf5yr ~ shf_bbl * ",
            intervar,
            " + offset(log(followup_overall))"
          )
        ))
      })

Error in eval(expr, p): object 'intervar' not found

Creating the formula outside of the with-environment and then pasting it creates the same error. Similar code has worked in previous projects, though. Thus, there might be a problem with the current mice version?

Expected behavior Without the bug, this binomial regression would work just as if I had called the interaction variable directly, e.g., shf_age_cat75:

      binomial_inter <- with(
      reprexdf_imputed_weighted,
      glm.nb(
      sos_out_countdeathcvhosphf5yr ~ shf_bbl * shf_age_cat75 + offset(log(followup_overall))
      )
      )

This piece of code works - but calling each interaction variable separately produces a large amount of repetitive code...

Reprex

#Packages
library(tidyverse)
library(tableone)
library(labelled)
library(openxlsx)
library(mice)
#> 
#> Attaching package: 'mice'
#> The following object is masked from 'package:stats':
#> 
#>     filter
#> The following objects are masked from 'package:base':
#> 
#>     cbind, rbind
library(future)
library(furrr)
library(survival)
#> 
#> Attaching package: 'survival'
#> The following object is masked from 'package:future':
#> 
#>     cluster
library(MatchThem)
#> 
#> Attaching package: 'MatchThem'
#> The following objects are masked from 'package:mice':
#> 
#>     cbind, pool
#> The following object is masked from 'package:base':
#> 
#>     cbind
library(cobalt)
#>  cobalt (Version 4.5.5, Build Date: 2024-04-02)
library(SuperLearner)
#> Loading required package: nnls
#> Loading required package: gam
#> Loading required package: splines
#> Loading required package: foreach
#> 
#> Attaching package: 'foreach'
#> The following objects are masked from 'package:purrr':
#> 
#>     accumulate, when
#> Loaded gam 1.22-5
#> Super Learner
#> Version: 2.0-29
#> Package created on 2024-02-06
library(spatstat)
#> Loading required package: spatstat.data
#> Loading required package: spatstat.univar
#> spatstat.univar 3.0-1
#> Loading required package: spatstat.geom
#> spatstat.geom 3.3-3
#> Loading required package: spatstat.random
#> spatstat.random 3.3-2
#> Loading required package: spatstat.explore
#> Loading required package: nlme
#> 
#> Attaching package: 'nlme'
#> The following object is masked from 'package:dplyr':
#> 
#>     collapse
#> spatstat.explore 3.3-2
#> 
#> Attaching package: 'spatstat.explore'
#> The following object is masked from 'package:MatchThem':
#> 
#>     pool
#> The following object is masked from 'package:mice':
#> 
#>     pool
#> Loading required package: spatstat.model
#> Loading required package: rpart
#> spatstat.model 3.3-2
#> 
#> Attaching package: 'spatstat.model'
#> The following object is masked from 'package:mice':
#> 
#>     ic
#> Loading required package: spatstat.linnet
#> spatstat.linnet 3.2-2
#> 
#> spatstat 3.2-1 
#> For an introduction to spatstat, type 'beginner'
library(WeightIt)
library(arm)
#> Loading required package: MASS
#> 
#> Attaching package: 'MASS'
#> The following object is masked from 'package:spatstat.geom':
#> 
#>     area
#> The following object is masked from 'package:dplyr':
#> 
#>     select
#> Loading required package: Matrix
#> 
#> Attaching package: 'Matrix'
#> The following objects are masked from 'package:tidyr':
#> 
#>     expand, pack, unpack
#> Loading required package: lme4
#> 
#> Attaching package: 'lme4'
#> The following object is masked from 'package:nlme':
#> 
#>     lmList
#> 
#> arm (Version 1.14-4, built: 2024-4-1)
#> Working directory is /private/var/folders/g7/syjhtqd51z7gj81rl96m7bjc0000gn/T/RtmpaVYgfZ/reprex-278e6ba19228-irate-booby
#> 
#> Attaching package: 'arm'
#> The following object is masked from 'package:spatstat.geom':
#> 
#>     rescale
library(bigmemory)
library(ncvreg)
library(biglasso)
library(earth)
#> Loading required package: Formula
#> Loading required package: plotmo
#> Loading required package: plotrix
#> 
#> Attaching package: 'plotrix'
#> The following object is masked from 'package:arm':
#> 
#>     rescale
#> The following objects are masked from 'package:spatstat.geom':
#> 
#>     hexagon, rescale
library(KernelKnn)
library(kernlab)
#> 
#> Attaching package: 'kernlab'
#> The following object is masked from 'package:mice':
#> 
#>     convergence
#> The following object is masked from 'package:purrr':
#> 
#>     cross
#> The following object is masked from 'package:ggplot2':
#> 
#>     alpha
library(randomForest)
#> randomForest 4.7-1.2
#> Type rfNews() to see new features/changes/bug fixes.
#> 
#> Attaching package: 'randomForest'
#> The following object is masked from 'package:dplyr':
#> 
#>     combine
#> The following object is masked from 'package:ggplot2':
#> 
#>     margin
library(speedglm)
#> Loading required package: biglm
#> Loading required package: DBI
library(xgboost)
#> 
#> Attaching package: 'xgboost'
#> The following object is masked from 'package:dplyr':
#> 
#>     slice
library(caret)
#> Loading required package: lattice
#> 
#> Attaching package: 'lattice'
#> The following object is masked from 'package:spatstat.model':
#> 
#>     panel.histogram
#> The following object is masked from 'package:spatstat.explore':
#> 
#>     panel.histogram
#> 
#> Attaching package: 'caret'
#> The following object is masked from 'package:survival':
#> 
#>     cluster
#> The following object is masked from 'package:future':
#> 
#>     cluster
#> The following object is masked from 'package:purrr':
#> 
#>     lift
library(lattice)
library(glue)
#> 
#> Attaching package: 'glue'
#> The following object is masked from 'package:WeightIt':
#> 
#>     trim
#> The following object is masked from 'package:MatchThem':
#> 
#>     trim
library(htmlwidgets)
library(DiagrammeR)
library(survminer)
#> Loading required package: ggpubr
#> 
#> Attaching package: 'ggpubr'
#> The following objects are masked from 'package:spatstat.geom':
#> 
#>     border, rotate
#> 
#> Attaching package: 'survminer'
#> The following object is masked from 'package:survival':
#> 
#>     myeloma
library(adjustedCurves)
library(patchwork)
#> 
#> Attaching package: 'patchwork'
#> The following object is masked from 'package:MASS':
#> 
#>     area
#> The following object is masked from 'package:spatstat.geom':
#> 
#>     area
library(data.table)
#> 
#> Attaching package: 'data.table'
#> The following object is masked from 'package:spatstat.geom':
#> 
#>     shift
#> The following objects are masked from 'package:lubridate':
#> 
#>     hour, isoweek, mday, minute, month, quarter, second, wday, week,
#>     yday, year
#> The following objects are masked from 'package:dplyr':
#> 
#>     between, first, last
#> The following object is masked from 'package:purrr':
#> 
#>     transpose
library(forestplot)
#> Loading required package: grid
#> 
#> Attaching package: 'grid'
#> The following object is masked from 'package:spatstat.geom':
#> 
#>     as.mask
#> Loading required package: checkmate
#> Loading required package: abind
library(gridExtra)
#> 
#> Attaching package: 'gridExtra'
#> The following object is masked from 'package:randomForest':
#> 
#>     combine
#> The following object is masked from 'package:dplyr':
#> 
#>     combine

#Show session info
sessionInfo()
#> R version 4.4.1 (2024-06-14)
#> Platform: aarch64-apple-darwin20
#> Running under: macOS 15.1
#> 
#> Matrix products: default
#> BLAS:   /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRblas.0.dylib 
#> LAPACK: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.0
#> 
#> locale:
#> [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
#> 
#> time zone: Europe/Stockholm
#> tzcode source: internal
#> 
#> attached base packages:
#> [1] grid      splines   stats     graphics  grDevices utils     datasets 
#> [8] methods   base     
#> 
#> other attached packages:
#>  [1] gridExtra_2.3          forestplot_3.1.5       abind_1.4-8           
#>  [4] checkmate_2.3.2        data.table_1.16.0      patchwork_1.3.0       
#>  [7] adjustedCurves_0.11.2  survminer_0.4.9        ggpubr_0.6.0          
#> [10] DiagrammeR_1.0.11      htmlwidgets_1.6.4      glue_1.8.0            
#> [13] caret_6.0-94           lattice_0.22-6         xgboost_1.7.8.1       
#> [16] speedglm_0.3-5         biglm_0.9-3            DBI_1.2.3             
#> [19] randomForest_4.7-1.2   kernlab_0.9-33         KernelKnn_1.1.5       
#> [22] earth_5.3.3            plotmo_3.6.4           plotrix_3.8-4         
#> [25] Formula_1.2-5          biglasso_1.6.0         ncvreg_3.14.3         
#> [28] bigmemory_4.6.4        arm_1.14-4             lme4_1.1-35.5         
#> [31] Matrix_1.7-0           MASS_7.3-61            WeightIt_1.3.0        
#> [34] spatstat_3.2-1         spatstat.linnet_3.2-2  spatstat.model_3.3-2  
#> [37] rpart_4.1.23           spatstat.explore_3.3-2 nlme_3.1-166          
#> [40] spatstat.random_3.3-2  spatstat.geom_3.3-3    spatstat.univar_3.0-1 
#> [43] spatstat.data_3.1-2    SuperLearner_2.0-29    gam_1.22-5            
#> [46] foreach_1.5.2          nnls_1.5               cobalt_4.5.5          
#> [49] MatchThem_1.2.1        survival_3.7-0         furrr_0.3.1           
#> [52] future_1.34.0          mice_3.16.16           openxlsx_4.2.7.1      
#> [55] labelled_2.13.0        tableone_0.13.2        lubridate_1.9.3       
#> [58] forcats_1.0.0          stringr_1.5.1          dplyr_1.1.4           
#> [61] purrr_1.0.2            readr_2.1.5            tidyr_1.3.1           
#> [64] tibble_3.2.1           ggplot2_3.5.1          tidyverse_2.0.0       
#> 
#> loaded via a namespace (and not attached):
#>  [1] polyclip_1.10-7       hardhat_1.4.0         pROC_1.18.5          
#>  [4] reprex_2.1.1          lifecycle_1.0.4       rstatix_0.7.2        
#>  [7] globals_0.16.3        backports_1.5.0       survey_4.4-2         
#> [10] magrittr_2.0.3        rmarkdown_2.28        yaml_2.3.10          
#> [13] doRNG_1.8.6           zip_2.3.1             spatstat.sparse_3.1-0
#> [16] minqa_1.2.8           RColorBrewer_1.1-3    nnet_7.3-19          
#> [19] ipred_0.9-15          lava_1.8.0            KMsurv_0.1-5         
#> [22] listenv_0.9.1         spatstat.utils_3.1-0  goftest_1.2-3        
#> [25] parallelly_1.38.0     codetools_0.2-20      tidyselect_1.2.1     
#> [28] shape_1.4.6.1         farver_2.1.2          stats4_4.4.1         
#> [31] jsonlite_1.8.9        mitml_0.4-5           iterators_1.0.14     
#> [34] tools_4.4.1           Rcpp_1.0.13-1         prodlim_2024.06.25   
#> [37] pan_1.9               xfun_0.47             mgcv_1.9-1           
#> [40] withr_3.0.2           fastmap_1.2.0         mitools_2.4          
#> [43] boot_1.3-31           fansi_1.0.6           digest_0.6.37        
#> [46] timechange_0.3.0      R6_2.5.1              colorspace_2.1-1     
#> [49] tensor_1.5            utf8_1.2.4            generics_0.1.3       
#> [52] recipes_1.1.0         class_7.3-22          ModelMetrics_1.2.2.2 
#> [55] pkgconfig_2.0.3       gtable_0.3.5          timeDate_4041.110    
#> [58] survMisc_0.5.6        htmltools_0.5.8.1     carData_3.0-5        
#> [61] scales_1.3.0          gower_1.0.1           bigmemory.sri_0.1.8  
#> [64] knitr_1.48            km.ci_0.5-6           rstudioapi_0.16.0    
#> [67] tzdb_0.4.0            reshape2_1.4.4        uuid_1.2-1           
#> [70] coda_0.19-4.1         visNetwork_2.1.2      nloptr_2.1.1         
#> [73] zoo_1.8-12            parallel_4.4.1        pillar_1.9.0         
#> [76] vctrs_0.6.5           car_3.1-2             jomo_2.7-6           
#> [79] xtable_1.8-4          evaluate_1.0.0        cli_3.6.3            
#> [82] compiler_4.4.1        rlang_1.1.4           crayon_1.5.3         
#> [85] rngtools_1.5.2        MatchIt_4.5.5         future.apply_1.11.2  
#> [88] ggsignif_0.6.4        plyr_1.8.9            fs_1.6.4             
#> [91] stringi_1.8.4         deldir_2.0-4          munsell_0.5.1        
#> [94] glmnet_4.1-8          hms_1.1.3             haven_2.5.4          
#> [97] broom_1.0.7
system("sysctl -n machdep.cpu.brand_string", intern = TRUE)
#> [1] "Apple M2"

#Show mice version
packageVersion("mice")
#> [1] '3.16.16'

#Create dummy dataframe reprexdf
set.seed(42)
n_obs <- 100

lopnr <- 1:n_obs

shf_age_cat75 <- sample(c("≤ 75", "> 75"), n_obs, replace = TRUE)
shf_age_cat75[sample(1:n_obs, size = floor(0.1 * n_obs))] <- NA  # 10% fehlende Werte

shf_sex <- sample(c("Female", "Male"), n_obs, replace = TRUE)
shf_sex[sample(1:n_obs, size = floor(0.1 * n_obs))] <- NA  # 10% fehlende Werte

shf_ef_hfref <- sample(c("< 30", "30-39"), n_obs, replace = TRUE)
shf_ef_hfref[sample(1:n_obs, size = floor(0.1 * n_obs))] <- NA  # 10% fehlende Werte

sos_out_countdeathcvhosphf5yr <- sample(0:10, n_obs, replace = TRUE)

shf_bbl <- sample(c("No", "Yes"), n_obs, replace = TRUE)

followup_overall <- round(runif(n_obs, 1.0, 5.0), 2)

reprexdf <- data.frame(
  lopnr = lopnr,
  shf_age_cat75 = factor(shf_age_cat75),
  shf_sex = factor(shf_sex),
  shf_ef_hfref = factor(shf_ef_hfref),
  sos_out_countdeathcvhosphf5yr = sos_out_countdeathcvhosphf5yr,
  shf_bbl = factor(shf_bbl),
  followup_overall = followup_overall
)

class(reprexdf)
#> [1] "data.frame"
summary(reprexdf)
#>      lopnr        shf_age_cat75   shf_sex   shf_ef_hfref
#>  Min.   :  1.00   > 75:51       Female:37   < 30 :53    
#>  1st Qu.: 25.75   ≤ 75:39       Male  :53   30-39:37    
#>  Median : 50.50   NA's:10       NA's  :10   NA's :10    
#>  Mean   : 50.50                                         
#>  3rd Qu.: 75.25                                         
#>  Max.   :100.00                                         
#>  sos_out_countdeathcvhosphf5yr shf_bbl  followup_overall
#>  Min.   : 0.00                 No :51   Min.   :1.210   
#>  1st Qu.: 2.00                 Yes:49   1st Qu.:2.215   
#>  Median : 5.00                          Median :2.920   
#>  Mean   : 5.13                          Mean   :3.015   
#>  3rd Qu.: 8.00                          3rd Qu.:3.860   
#>  Max.   :10.00                          Max.   :4.980

#Impute data
reprexdf_imputed <- futuremice(
  reprexdf,
  m = 5,
  maxit = 5,
  defaultMethod = c("pmm", "logreg", "polyreg", "polr"),
  parallelseed = 234,
  n.core = (parallel::detectCores() - 1)
)
#> Warning in check.cores(n.core, available, m): 'n.core' exceeds the maximum
#> number of available cores on your machine or the number of imputations, and is
#> set to 5

class(reprexdf_imputed)
#> [1] "mids"

#Apply IPTW (SuperLearner)
propensity_covariates <- c(
  "shf_age_cat75",
  "shf_sex",
  "shf_ef_hfref",
  "sos_out_countdeathcvhosphf5yr",
  "followup_overall"
)

plan(multisession, workers = parallel::detectCores() - 1)

reprexdf_imputed_weighted <-
  weightthem(
    as.formula(paste(
      "shf_bbl ~ ", paste(propensity_covariates, collapse = " + ")
    )),
    data = reprexdf_imputed,
    approach = "across",
    method = "super",
    stabilize = FALSE,
    estimand = "ATO",
    SL.library = c(
      "SL.glmnet",
      "SL.randomForest",
      "SL.xgboost",
      "SL.speedglm",
      "SL.mean",
      "SL.bayesglm",
      "SL.lm",
      "SL.glm"
    ),
    SL.method = "method.balance",
    criterion = "smd.mean"
  )
#> Estimating distances | dataset: #1
#>  #2 #3 #4 #5
#> Estimating weights     | dataset: #1 #2 #3 #4 #5

class(reprexdf_imputed_weighted)
#> [1] "wimids"

#Create main analysis function
test_number_function <- function(reprexdf_imputed_weighted) {
  interaction_variables <- c("shf_age_cat75", "shf_sex", "shf_ef_hfref")

  for (i in 1:length(interaction_variables)) {
    intervar <- interaction_variables[[i]]

    cat(paste0("\n**---Interaction analysis ", intervar, " ---**\n"))

    #THIS PIECE OF CODE WORKS (no dynamic variable)
    # binomial_inter <- with(
    #   reprexdf_imputed_weighted,
    #   glm.nb(
    #     sos_out_countdeathcvhosphf5yr ~ shf_bbl * shf_age_cat75 + offset(log(followup_overall))
    #   )
    # )

    #THIS PIECE OF CODE DOES NOT WORK (dynamic variable)
    binomial_inter <- with(reprexdf_imputed_weighted, {
      glm.nb(as.formula(
        paste(
          "sos_out_countdeathcvhosphf5yr ~ shf_bbl * ",
          intervar,
          " + offset(log(followup_overall))"
        )
      ))
    })

    summary(mice::pool(binomial_inter),
            conf.int = TRUE,
            exp = TRUE) %>% print()

    i <- i + 1
  }
}

#Call the function
test_number_function(reprexdf_imputed_weighted)
#> 
#> **---Interaction analysis shf_age_cat75 ---**
#> Error in eval(expr, p): object 'intervar' not found
stefvanbuuren commented 2 weeks ago

I did not foresee this use case. Did you already try to define intervar in the .GlobalEnv?

benebeer commented 2 weeks ago

I created an empty character variable in the GlobalEnv now, before starting the for-loop: intervar <- character()

Works perfectly now - thanks so much! Didn't think about this workaround...

benebeer commented 2 weeks ago

Unfortunately, I need to reopen this issue. Creating intervar in the GlobalEnv has "worked" in the sense that I do not receive the error "object 'intervar' not found" anymore. The code runs. BUT: When I look at the results (main analysis without interaction term and interaction analyses using the for-loop), they are all the same. And the interaction estimate is missing. It seems that the formula is created correctly, but the glm.nb-regression within the with-function might just use the intervar in the GlobalEnv (= empty)?

(I am not sure if this still qualifies as a bug or rather as a feature request. Please feel free to change the tag if necessary.)

Reprex

#Packages
  library(tidyverse)
  library(tableone)
  library(labelled)
  library(openxlsx)
  library(mice)
#> 
#> Attaching package: 'mice'
#> The following object is masked from 'package:stats':
#> 
#>     filter
#> The following objects are masked from 'package:base':
#> 
#>     cbind, rbind
  library(future)
  library(furrr)
  library(survival)
#> 
#> Attaching package: 'survival'
#> The following object is masked from 'package:future':
#> 
#>     cluster
  library(MatchThem)
#> 
#> Attaching package: 'MatchThem'
#> The following objects are masked from 'package:mice':
#> 
#>     cbind, pool
#> The following object is masked from 'package:base':
#> 
#>     cbind
  library(cobalt)
#>  cobalt (Version 4.5.5, Build Date: 2024-04-02)
  library(SuperLearner)
#> Loading required package: nnls
#> Loading required package: gam
#> Loading required package: splines
#> Loading required package: foreach
#> 
#> Attaching package: 'foreach'
#> The following objects are masked from 'package:purrr':
#> 
#>     accumulate, when
#> Loaded gam 1.22-5
#> Super Learner
#> Version: 2.0-29
#> Package created on 2024-02-06
  library(spatstat)
#> Loading required package: spatstat.data
#> Loading required package: spatstat.univar
#> spatstat.univar 3.0-1
#> Loading required package: spatstat.geom
#> spatstat.geom 3.3-3
#> Loading required package: spatstat.random
#> spatstat.random 3.3-2
#> Loading required package: spatstat.explore
#> Loading required package: nlme
#> 
#> Attaching package: 'nlme'
#> The following object is masked from 'package:dplyr':
#> 
#>     collapse
#> spatstat.explore 3.3-2
#> 
#> Attaching package: 'spatstat.explore'
#> The following object is masked from 'package:MatchThem':
#> 
#>     pool
#> The following object is masked from 'package:mice':
#> 
#>     pool
#> Loading required package: spatstat.model
#> Loading required package: rpart
#> spatstat.model 3.3-2
#> 
#> Attaching package: 'spatstat.model'
#> The following object is masked from 'package:mice':
#> 
#>     ic
#> Loading required package: spatstat.linnet
#> spatstat.linnet 3.2-2
#> 
#> spatstat 3.2-1 
#> For an introduction to spatstat, type 'beginner'
  library(WeightIt)
  library(arm)
#> Loading required package: MASS
#> 
#> Attaching package: 'MASS'
#> The following object is masked from 'package:spatstat.geom':
#> 
#>     area
#> The following object is masked from 'package:dplyr':
#> 
#>     select
#> Loading required package: Matrix
#> 
#> Attaching package: 'Matrix'
#> The following objects are masked from 'package:tidyr':
#> 
#>     expand, pack, unpack
#> Loading required package: lme4
#> 
#> Attaching package: 'lme4'
#> The following object is masked from 'package:nlme':
#> 
#>     lmList
#> 
#> arm (Version 1.14-4, built: 2024-4-1)
#> Working directory is /private/var/folders/g7/syjhtqd51z7gj81rl96m7bjc0000gn/T/Rtmpl341Ic/reprex-c8ee8157e0b-daft-loon
#> 
#> Attaching package: 'arm'
#> The following object is masked from 'package:spatstat.geom':
#> 
#>     rescale
  library(bigmemory)
  library(ncvreg)
  library(biglasso)
  library(earth)
#> Loading required package: Formula
#> Loading required package: plotmo
#> Loading required package: plotrix
#> 
#> Attaching package: 'plotrix'
#> The following object is masked from 'package:arm':
#> 
#>     rescale
#> The following objects are masked from 'package:spatstat.geom':
#> 
#>     hexagon, rescale
  library(KernelKnn)
  library(kernlab)
#> 
#> Attaching package: 'kernlab'
#> The following object is masked from 'package:mice':
#> 
#>     convergence
#> The following object is masked from 'package:purrr':
#> 
#>     cross
#> The following object is masked from 'package:ggplot2':
#> 
#>     alpha
  library(randomForest)
#> randomForest 4.7-1.2
#> Type rfNews() to see new features/changes/bug fixes.
#> 
#> Attaching package: 'randomForest'
#> The following object is masked from 'package:dplyr':
#> 
#>     combine
#> The following object is masked from 'package:ggplot2':
#> 
#>     margin
  library(speedglm)
#> Loading required package: biglm
#> Loading required package: DBI
  library(xgboost)
#> 
#> Attaching package: 'xgboost'
#> The following object is masked from 'package:dplyr':
#> 
#>     slice
  library(caret)
#> Loading required package: lattice
#> 
#> Attaching package: 'lattice'
#> The following object is masked from 'package:spatstat.model':
#> 
#>     panel.histogram
#> The following object is masked from 'package:spatstat.explore':
#> 
#>     panel.histogram
#> 
#> Attaching package: 'caret'
#> The following object is masked from 'package:survival':
#> 
#>     cluster
#> The following object is masked from 'package:future':
#> 
#>     cluster
#> The following object is masked from 'package:purrr':
#> 
#>     lift
  library(lattice)
  library(glue)
#> 
#> Attaching package: 'glue'
#> The following object is masked from 'package:WeightIt':
#> 
#>     trim
#> The following object is masked from 'package:MatchThem':
#> 
#>     trim
  library(htmlwidgets)
  library(DiagrammeR)
  library(survminer)
#> Loading required package: ggpubr
#> 
#> Attaching package: 'ggpubr'
#> The following objects are masked from 'package:spatstat.geom':
#> 
#>     border, rotate
#> 
#> Attaching package: 'survminer'
#> The following object is masked from 'package:survival':
#> 
#>     myeloma
  library(adjustedCurves)
  library(patchwork)
#> 
#> Attaching package: 'patchwork'
#> The following object is masked from 'package:MASS':
#> 
#>     area
#> The following object is masked from 'package:spatstat.geom':
#> 
#>     area
  library(data.table)
#> 
#> Attaching package: 'data.table'
#> The following object is masked from 'package:spatstat.geom':
#> 
#>     shift
#> The following objects are masked from 'package:lubridate':
#> 
#>     hour, isoweek, mday, minute, month, quarter, second, wday, week,
#>     yday, year
#> The following objects are masked from 'package:dplyr':
#> 
#>     between, first, last
#> The following object is masked from 'package:purrr':
#> 
#>     transpose
  library(forestplot)
#> Loading required package: grid
#> 
#> Attaching package: 'grid'
#> The following object is masked from 'package:spatstat.geom':
#> 
#>     as.mask
#> Loading required package: checkmate
#> Loading required package: abind
  library(gridExtra)
#> 
#> Attaching package: 'gridExtra'
#> The following object is masked from 'package:randomForest':
#> 
#>     combine
#> The following object is masked from 'package:dplyr':
#> 
#>     combine

  #Show session info
  sessionInfo()
#> R version 4.4.1 (2024-06-14)
#> Platform: aarch64-apple-darwin20
#> Running under: macOS 15.1
#> 
#> Matrix products: default
#> BLAS:   /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRblas.0.dylib 
#> LAPACK: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.0
#> 
#> locale:
#> [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
#> 
#> time zone: Europe/Stockholm
#> tzcode source: internal
#> 
#> attached base packages:
#> [1] grid      splines   stats     graphics  grDevices utils     datasets 
#> [8] methods   base     
#> 
#> other attached packages:
#>  [1] gridExtra_2.3          forestplot_3.1.5       abind_1.4-8           
#>  [4] checkmate_2.3.2        data.table_1.16.0      patchwork_1.3.0       
#>  [7] adjustedCurves_0.11.2  survminer_0.4.9        ggpubr_0.6.0          
#> [10] DiagrammeR_1.0.11      htmlwidgets_1.6.4      glue_1.8.0            
#> [13] caret_6.0-94           lattice_0.22-6         xgboost_1.7.8.1       
#> [16] speedglm_0.3-5         biglm_0.9-3            DBI_1.2.3             
#> [19] randomForest_4.7-1.2   kernlab_0.9-33         KernelKnn_1.1.5       
#> [22] earth_5.3.3            plotmo_3.6.4           plotrix_3.8-4         
#> [25] Formula_1.2-5          biglasso_1.6.0         ncvreg_3.14.3         
#> [28] bigmemory_4.6.4        arm_1.14-4             lme4_1.1-35.5         
#> [31] Matrix_1.7-0           MASS_7.3-61            WeightIt_1.3.0        
#> [34] spatstat_3.2-1         spatstat.linnet_3.2-2  spatstat.model_3.3-2  
#> [37] rpart_4.1.23           spatstat.explore_3.3-2 nlme_3.1-166          
#> [40] spatstat.random_3.3-2  spatstat.geom_3.3-3    spatstat.univar_3.0-1 
#> [43] spatstat.data_3.1-2    SuperLearner_2.0-29    gam_1.22-5            
#> [46] foreach_1.5.2          nnls_1.5               cobalt_4.5.5          
#> [49] MatchThem_1.2.1        survival_3.7-0         furrr_0.3.1           
#> [52] future_1.34.0          mice_3.16.16           openxlsx_4.2.7.1      
#> [55] labelled_2.13.0        tableone_0.13.2        lubridate_1.9.3       
#> [58] forcats_1.0.0          stringr_1.5.1          dplyr_1.1.4           
#> [61] purrr_1.0.2            readr_2.1.5            tidyr_1.3.1           
#> [64] tibble_3.2.1           ggplot2_3.5.1          tidyverse_2.0.0       
#> 
#> loaded via a namespace (and not attached):
#>  [1] polyclip_1.10-7       hardhat_1.4.0         pROC_1.18.5          
#>  [4] reprex_2.1.1          lifecycle_1.0.4       rstatix_0.7.2        
#>  [7] globals_0.16.3        backports_1.5.0       survey_4.4-2         
#> [10] magrittr_2.0.3        rmarkdown_2.28        yaml_2.3.10          
#> [13] doRNG_1.8.6           zip_2.3.1             spatstat.sparse_3.1-0
#> [16] minqa_1.2.8           RColorBrewer_1.1-3    nnet_7.3-19          
#> [19] ipred_0.9-15          lava_1.8.0            KMsurv_0.1-5         
#> [22] listenv_0.9.1         spatstat.utils_3.1-0  goftest_1.2-3        
#> [25] parallelly_1.38.0     codetools_0.2-20      tidyselect_1.2.1     
#> [28] shape_1.4.6.1         farver_2.1.2          stats4_4.4.1         
#> [31] jsonlite_1.8.9        mitml_0.4-5           iterators_1.0.14     
#> [34] tools_4.4.1           Rcpp_1.0.13-1         prodlim_2024.06.25   
#> [37] pan_1.9               xfun_0.47             mgcv_1.9-1           
#> [40] withr_3.0.2           fastmap_1.2.0         mitools_2.4          
#> [43] boot_1.3-31           fansi_1.0.6           digest_0.6.37        
#> [46] timechange_0.3.0      R6_2.5.1              colorspace_2.1-1     
#> [49] tensor_1.5            utf8_1.2.4            generics_0.1.3       
#> [52] recipes_1.1.0         class_7.3-22          ModelMetrics_1.2.2.2 
#> [55] pkgconfig_2.0.3       gtable_0.3.5          timeDate_4041.110    
#> [58] survMisc_0.5.6        htmltools_0.5.8.1     carData_3.0-5        
#> [61] scales_1.3.0          gower_1.0.1           bigmemory.sri_0.1.8  
#> [64] knitr_1.48            km.ci_0.5-6           rstudioapi_0.16.0    
#> [67] tzdb_0.4.0            reshape2_1.4.4        uuid_1.2-1           
#> [70] coda_0.19-4.1         visNetwork_2.1.2      nloptr_2.1.1         
#> [73] zoo_1.8-12            parallel_4.4.1        pillar_1.9.0         
#> [76] vctrs_0.6.5           car_3.1-2             jomo_2.7-6           
#> [79] xtable_1.8-4          evaluate_1.0.0        cli_3.6.3            
#> [82] compiler_4.4.1        rlang_1.1.4           crayon_1.5.3         
#> [85] rngtools_1.5.2        MatchIt_4.5.5         future.apply_1.11.2  
#> [88] ggsignif_0.6.4        plyr_1.8.9            fs_1.6.4             
#> [91] stringi_1.8.4         deldir_2.0-4          munsell_0.5.1        
#> [94] glmnet_4.1-8          hms_1.1.3             haven_2.5.4          
#> [97] broom_1.0.7
  system("sysctl -n machdep.cpu.brand_string", intern = TRUE)
#> [1] "Apple M2"

  #Show mice version
  packageVersion("mice")
#> [1] '3.16.16'

  #Create dummy dataframe reprexdf
  set.seed(42)
  n_obs <- 100

  lopnr <- 1:n_obs

  shf_age_cat75 <- sample(c("≤ 75", "> 75"), n_obs, replace = TRUE)
  shf_age_cat75[sample(1:n_obs, size = floor(0.1 * n_obs))] <- NA  # 10% fehlende Werte

  shf_sex <- sample(c("Female", "Male"), n_obs, replace = TRUE)
  shf_sex[sample(1:n_obs, size = floor(0.1 * n_obs))] <- NA  # 10% fehlende Werte

  shf_ef_hfref <- sample(c("< 30", "30-39"), n_obs, replace = TRUE)
  shf_ef_hfref[sample(1:n_obs, size = floor(0.1 * n_obs))] <- NA  # 10% fehlende Werte

  sos_out_countdeathcvhosphf5yr <- sample(0:10, n_obs, replace = TRUE)

  shf_bbl <- sample(c("No", "Yes"), n_obs, replace = TRUE)

  followup_overall <- round(runif(n_obs, 1.0, 5.0), 2)

  reprexdf <- data.frame(
    lopnr = lopnr,
    shf_age_cat75 = factor(shf_age_cat75),
    shf_sex = factor(shf_sex),
    shf_ef_hfref = factor(shf_ef_hfref),
    sos_out_countdeathcvhosphf5yr = sos_out_countdeathcvhosphf5yr,
    shf_bbl = factor(shf_bbl),
    followup_overall = followup_overall
  )

  class(reprexdf)
#> [1] "data.frame"
  summary(reprexdf)
#>      lopnr        shf_age_cat75   shf_sex   shf_ef_hfref
#>  Min.   :  1.00   > 75:51       Female:37   < 30 :53    
#>  1st Qu.: 25.75   ≤ 75:39       Male  :53   30-39:37    
#>  Median : 50.50   NA's:10       NA's  :10   NA's :10    
#>  Mean   : 50.50                                         
#>  3rd Qu.: 75.25                                         
#>  Max.   :100.00                                         
#>  sos_out_countdeathcvhosphf5yr shf_bbl  followup_overall
#>  Min.   : 0.00                 No :51   Min.   :1.210   
#>  1st Qu.: 2.00                 Yes:49   1st Qu.:2.215   
#>  Median : 5.00                          Median :2.920   
#>  Mean   : 5.13                          Mean   :3.015   
#>  3rd Qu.: 8.00                          3rd Qu.:3.860   
#>  Max.   :10.00                          Max.   :4.980

  #Impute data
  reprexdf_imputed <- futuremice(
    reprexdf,
    m = 5,
    maxit = 5,
    defaultMethod = c("pmm", "logreg", "polyreg", "polr"),
    parallelseed = 234,
    n.core = (parallel::detectCores() - 1)
  )
#> Warning in check.cores(n.core, available, m): 'n.core' exceeds the maximum
#> number of available cores on your machine or the number of imputations, and is
#> set to 5

  class(reprexdf_imputed)
#> [1] "mids"

  #Apply IPTW (SuperLearner)
  propensity_covariates <- c(
    "shf_age_cat75",
    "shf_sex",
    "shf_ef_hfref",
    "sos_out_countdeathcvhosphf5yr",
    "followup_overall"
  )

  plan(multisession, workers = parallel::detectCores() - 1)

  reprexdf_imputed_weighted <-
    weightthem(
      as.formula(paste(
        "shf_bbl ~ ", paste(propensity_covariates, collapse = " + ")
      )),
      data = reprexdf_imputed,
      approach = "across",
      method = "super",
      stabilize = FALSE,
      estimand = "ATO",
      SL.library = c(
        "SL.glmnet",
        "SL.randomForest",
        "SL.xgboost",
        "SL.speedglm",
        "SL.mean",
        "SL.bayesglm",
        "SL.lm",
        "SL.glm"
      ),
      SL.method = "method.balance",
      criterion = "smd.mean"
    )
#> Estimating distances | dataset: #1
#>  #2 #3 #4 #5
#> Estimating weights     | dataset: #1 #2 #3 #4 #5

  class(reprexdf_imputed_weighted)
#> [1] "wimids"

  #Create main analysis function
  intervar <- character() ##Create empty intervar in GlobalEnv, will be filled in for-loop

  test_number_function <- function(df_imputed_weighted) {
    interaction_variables <- c("shf_age_cat75", "shf_sex", "shf_ef_hfref")
    results_list <- list()

    cat("\n**---Main interaction analysis---**\n")
    summary(mice::pool(with(
      df_imputed_weighted,
      glm.nb(sos_out_countdeathcvhosphf5yr ~ shf_bbl + offset(log(followup_overall)))
    )), conf.int = TRUE, exp = TRUE) %>% print()

    cat("\n**---Interaction interaction analysis---**\n")
    for (i in 1:length(interaction_variables)) {
      intervar <- interaction_variables[i]

      ##Test formula composition
      formula <- as.formula(paste(
        "sos_out_countdeathcvhosphf5yr ~ shf_bbl * ", intervar, " + offset(log(followup_overall))"
      ))
      print(formula)

      interterm_results <- summary(mice::pool(with(
        df_imputed_weighted, glm.nb(as.formula(
          paste(
            "sos_out_countdeathcvhosphf5yr ~ shf_bbl * ",
            intervar,
            " + offset(log(followup_overall))"
          )
        ))
      )), conf.int = TRUE, exp = TRUE)

      results_list[[paste0("interterm_", intervar)]] <- interterm_results
    }
    return(results_list)
  }
test_number_function(reprexdf_imputed_weighted)

Output (results)

> ---Main interaction analysis---

> term estimate std.error statistic df p.value conf.low

> 1 (Intercept) 1.8631883 0.1470628 4.2314524 96.04978 5.309108e-05 1.3914881

> 2 shf_bblYes 0.9758789 0.2087355 -0.1169747 96.04978 9.071243e-01 0.6448415

> conf.high

> 1 2.494790

> 2 1.476859

>

> ---Interaction interaction analysis---

> sos_out_countdeathcvhosphf5yr ~ shf_bbl * shf_age_cat75 + offset(log(followup_overall))

> <environment: 0x12dcdd190>

> sos_out_countdeathcvhosphf5yr ~ shf_bbl * shf_sex + offset(log(followup_overall))

> <environment: 0x12dcdd190>

> sos_out_countdeathcvhosphf5yr ~ shf_bbl * shf_ef_hfref + offset(log(followup_overall))

> <environment: 0x12dcdd190>

> $interterm_shf_age_cat75

> term estimate std.error statistic df p.value conf.low

> 1 (Intercept) 1.8631883 0.1470628 4.2314524 96.04978 5.309108e-05 1.3914881

> 2 shf_bblYes 0.9758789 0.2087355 -0.1169747 96.04978 9.071243e-01 0.6448415

> conf.high

> 1 2.494790

> 2 1.476859

>

> $interterm_shf_sex

> term estimate std.error statistic df p.value conf.low

> 1 (Intercept) 1.8631883 0.1470628 4.2314524 96.04978 5.309108e-05 1.3914881

> 2 shf_bblYes 0.9758789 0.2087355 -0.1169747 96.04978 9.071243e-01 0.6448415

> conf.high

> 1 2.494790

> 2 1.476859

>

> $interterm_shf_ef_hfref

> term estimate std.error statistic df p.value conf.low

> 1 (Intercept) 1.8631883 0.1470628 4.2314524 96.04978 5.309108e-05 1.3914881

> 2 shf_bblYes 0.9758789 0.2087355 -0.1169747 96.04978 9.071243e-01 0.6448415

> conf.high

> 1 2.494790

> 2 1.476859

stefvanbuuren commented 1 week ago

There is no need to create intervar explicitly.

Try changing the local assigment <- into the global assignment <<-. I guess intervar <<- interaction_variables[i] should do the job.

benebeer commented 1 week ago

intervar <<- interaction_variables[i] solved the issue, indeed. Thanks a lot!!