search

leeper / margins

An R Port of Stata's 'margins' Command
https://cloud.r-project.org/package=margins
Other
263 stars 40 forks source link

margins() doesn't work for splines/polynomial predictors in (g)lmer models? #171

Open msonderegger opened 3 years ago

msonderegger commented 3 years ago

Please specify whether your issue is about:

Thank you for the excellent package! I think there may be a bug in (g)lmer models with non-linear functions of predictors.

Here is a minimal reproducible example:

library("margins")
library("lme4")

fm1 <- lmer(Reaction ~ ns(Days,2) + (Days|Subject), sleepstudy)
margins(fm1)

Error in data.table::rbindlist(list(d0, d1)) : 
  Column 1 of item 1 is length 180 inconsistent with column 2 which is length 360. Only length-1 columns are recycled.

The same error occurs if I use splines::ns in the model formula instead of ns. The same issue happens using stats::poly instead of splines::ns, but not using log(Days+1) instead of ns(Days,2). However ns, poly, etc. functionality does work for lm models in margins, which makes me think this may be a bug?

Session info and traceback:

> traceback()
11: data.table::rbindlist(list(d0, d1))
10: prediction.merMod(model = model, data = data.table::rbindlist(list(d0, 
        d1)), type = type, calculate_se = FALSE, ...)
9: prediction(model = model, data = data.table::rbindlist(list(d0, 
       d1)), type = type, calculate_se = FALSE, ...)
8: dydx.default(X[[i]], ...)
7: FUN(X[[i]], ...)
6: lapply(c(varslist$nnames, varslist$lnames), dydx, data = data, 
       model = model, type = type, eps = eps, as.data.frame = as.data.frame, 
       ...)
5: marginal_effects.lmerMod(model = model, data = data, variables = variables, 
       type = type, eps = eps, varslist = varslist, ...)
4: marginal_effects(model = model, data = data, variables = variables, 
       type = type, eps = eps, varslist = varslist, ...)
3: build_margins(model = model, data = data_list[[i]], variables = variables, 
       type = type, vcov = vcov, vce = vce, iterations = iterations, 
       unit_ses = unit_ses, eps = eps, varslist = varslist, ...)
2: margins.lmerMod(fm1)
1: margins(fm1)

> sessionInfo()
R version 4.0.3 (2020-10-10)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Big Sur 10.16

Matrix products: default
LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib

Random number generation:
 RNG:     Mersenne-Twister 
 Normal:  Inversion 
 Sample:  Rounding 

locale:
[1] en_CA.UTF-8/en_CA.UTF-8/en_CA.UTF-8/C/en_CA.UTF-8/en_CA.UTF-8

attached base packages:
[1] splines   stats     graphics  grDevices utils     datasets 
[7] methods   base     

other attached packages:
 [1] margins_0.3.26  forcats_0.5.1   stringr_1.4.0  
 [4] dplyr_1.0.4     purrr_0.3.4     readr_1.4.0    
 [7] tidyr_1.1.2     tibble_3.0.6    ggplot2_3.3.3  
[10] tidyverse_1.3.0 ggeffects_1.0.1 lmerTest_3.1-3 
[13] lme4_1.1-26     Matrix_1.3-2    emmeans_1.5.4  

loaded via a namespace (and not attached):
 [1] httr_1.4.2          jsonlite_1.7.2      modelr_0.1.8       
 [4] assertthat_0.2.1    statmod_1.4.35      cellranger_1.1.0   
 [7] numDeriv_2016.8-1.1 pillar_1.4.7        backports_1.2.1    
[10] lattice_0.20-41     glue_1.4.2          digest_0.6.27      
[13] rvest_0.3.6         snakecase_0.11.0    minqa_1.2.4        
[16] colorspace_2.0-0    sandwich_3.0-0      pkgconfig_2.0.3    
[19] broom_0.7.4         haven_2.3.1         xtable_1.8-4       
[22] mvtnorm_1.1-1       scales_1.1.1        generics_0.1.0     
[25] farver_2.0.3        sjlabelled_1.1.7    ellipsis_0.3.1     
[28] TH.data_1.0-10      withr_2.4.1         cli_2.3.0          
[31] survival_3.2-7      magrittr_2.0.1      crayon_1.4.0       
[34] readxl_1.3.1        estimability_1.3    fs_1.5.0           
[37] fansi_0.4.2         nlme_3.1-149        MASS_7.3-53        
[40] xml2_1.3.2          data.table_1.13.6   tools_4.0.3        
[43] hms_1.0.0           lifecycle_0.2.0     multcomp_1.4-15    
[46] munsell_0.5.0       reprex_1.0.0        prediction_0.3.14  
[49] packrat_0.5.0       compiler_4.0.3      rlang_0.4.10       
[52] grid_4.0.3          nloptr_1.2.2.2      rstudioapi_0.13    
[55] labeling_0.4.2      boot_1.3-25         gtable_0.3.0       
[58] codetools_0.2-16    DBI_1.1.1           R6_2.5.0           
[61] zoo_1.8-8           lubridate_1.7.9.2   knitr_1.31         
[64] utf8_1.1.4          insight_0.12.0      stringi_1.5.3      
[67] pscl_1.5.5          Rcpp_1.0.6          vctrs_0.3.6        
[70] dbplyr_2.1.0        tidyselect_1.1.0    xfun_0.20          
[73] coda_0.19-4
msonderegger commented 3 years ago

In case anyone else has this issue -- I think an equivalent manual workaround is to use marginal_effects (which works), then take column means. For the example above:

marginal_effects(fm1, sleepstudy) %>% lapply(.,mean)
$dydx_Days
[1] 10.46729