lsmeans compatibility with mixed() and type 2 tests

Home forums Mixed Models lsmeans compatibility with mixed() and type 2 tests

Tagged: , ,

This topic contains 6 replies, has 2 voices, and was last updated by  henrik 1 month, 3 weeks ago.

  • Author
    Posts
  • #114

    Elvio Blini
    Participant

    Hi Henrik,

    thanks for your time.
    I’m unsure of lsmeans compatibility with mixed objects following type 2 tests. Here’s an example:

    library(afex)
    
    sessionInfo()
    
    # R version 3.4.0 (2017-04-21)
    # Platform: x86_64-w64-mingw32/x64 (64-bit)
    # Running under: Windows >= 8 x64 (build 9200)
    # 
    # Matrix products: default
    # 
    # locale:
    #   [1] LC_COLLATE=English_United Kingdom.1252  LC_CTYPE=English_United Kingdom.1252    LC_MONETARY=English_United Kingdom.1252
    # [4] LC_NUMERIC=C                            LC_TIME=English_United Kingdom.1252    
    # 
    # attached base packages:
    #   [1] stats     graphics  grDevices utils     datasets  methods   base     
    # 
    # other attached packages:
    #   [1] afex_0.18-0      lsmeans_2.26-3   estimability_1.2 lme4_1.1-13      Matrix_1.2-9    
    # 
    # loaded via a namespace (and not attached):
    #   [1] zoo_1.8-0           modeltools_0.2-21   coin_1.1-2          reshape2_1.4.2      splines_3.4.0       lmerTest_2.0-33     lattice_0.20-35    
    # [8] colorspace_1.3-2    htmltools_0.3.6     stats4_3.4.0        mgcv_1.8-17         base64enc_0.1-3     survival_2.41-3     rlang_0.1.1        
    # [15] nloptr_1.0.4        foreign_0.8-69      RColorBrewer_1.1-2  multcomp_1.4-6      plyr_1.8.4          stringr_1.2.0       MatrixModels_0.4-1 
    # [22] munsell_0.4.3       gtable_0.2.0        htmlwidgets_0.9     mvtnorm_1.0-6       coda_0.19-1         codetools_0.2-15    knitr_1.16         
    # [29] latticeExtra_0.6-28 SparseM_1.77        quantreg_5.33       pbkrtest_0.4-7      parallel_3.4.0      htmlTable_1.9       TH.data_1.0-8      
    # [36] Rcpp_0.12.4         xtable_1.8-2        acepack_1.4.1       scales_0.4.1        backports_1.1.0     checkmate_1.8.3     Hmisc_4.0-3        
    # [43] gridExtra_2.2.1     digest_0.6.12       ggplot2_2.2.1       stringi_1.1.5       grid_3.4.0          tools_3.4.0         magrittr_1.5       
    # [50] sandwich_2.4-0      lazyeval_0.2.0      tibble_1.3.3        Formula_1.2-2       cluster_2.0.6       car_2.1-5           MASS_7.3-47        
    # [57] data.table_1.10.4   minqa_1.2.4         rpart_4.1-11        nnet_7.3-12         nlme_3.1-131        compiler_3.4.0 
    
    #example 1
    
    data(md_16.4)
    str(md_16.4)
    
    #works
    (mixed1 <- mixed(induct ~ cond + (1|room:cond), md_16.4))
    lsmeans(mixed1, specs = c("cond"))
    lsmeans(mixed1, pairwise ~ cond)$contrast
    
    #does not work
    (mixed2 <- mixed(induct ~ cond + (1|room:cond), md_16.4, type=2, method="LRT"))
    lsmeans(mixed2, specs = c("cond"))
    lsmeans(mixed2, ~cond)
    lsmeans(mixed2, pairwise ~ cond)$contrast
    
    # 
    # Error in UseMethod("recover.data") : 
    #   no applicable method for 'recover.data' applied to an object of class "list"
    # Error in ref.grid(object = list(anova_table = list(Df = 3, Chisq = 3.52370343128746,  : 
    #                                                      Perhaps a 'data' or 'params' argument is needed
    
    #providing full model only does not work either
    lsmeans(mixed2$full_model, "cond")
    # Error in ref.grid(object = list(<S4 object of class "merModLmerTest">)) : 
    #   Can't handle an object of class  “list” 
    #  Use help("models", package = "lsmeans") for information on supported models.
    
    #works
    (mixed3 <- mixed(induct ~ cond + (1|room:cond), md_16.4, type=3, method="LRT"))
    lsmeans(mixed3, specs = c("cond"))
    
    #doesn't work for type 2 tests?
    
    

    Hope it will work out…

    Thanks again!

    Elvio

  • #115

    henrik
    Keymaster

    That is a indeed a bug when passing mixed objects with type=2 and method="LRT" (or "PB"). It is now corrected in the development version of afex which can be installed via: devtools::install_github("singmann/afex@master")

    If it is not possible to use the development version, a work-around along the following lines should work:

    
    require(afex)
    data(md_16.4)
    (mixed2 <- mixed(induct ~ cond + (1|room:cond), md_16.4, type=2, method="LRT"))
    lsmeans(mixed2$full_model[[1]], specs = c("cond"))
    lsmeans(mixed2$full_model[[1]], ~cond)
    lsmeans(mixed2$full_model[[1]], pairwise ~ cond)$contrast

    As you can see, we simply directly pass the full model. The problem for type=2 is that in case in which nested model comparisons are used (e.g., when method = "LRT") the full_model slot is a list of lists, so we need to use the first element of the slot.

    Thanks again for the report.

  • #116

    Elvio Blini
    Participant

    Oh, I see! Thank you very much for your quick reply and fix!

  • #119

    Elvio Blini
    Participant

    Hi!

    Just to add up: in the case of multiple independent factors the work-around would be:

    
    lsmeans(tail(svv_models$full_model, 1)[[1]], specs = c("cond"))
    

    Correct?

    This would avoid this error:

    
    # Error in array(seq_len(nrow(RG@linfct)), dims) : 
    #   'dims' cannot be of length 0
    # In addition: Warning message:
    # In model.matrix.default(trms, m, contrasts.arg = contrasts) :
    #   problem with term 1 in model.matrix: no columns are assigned
    

    (I’ve switched to the github version but I still need this work-around when testing multiple IV, 🙂 )

    Best,
    Elvio

    • #120

      henrik
      Keymaster

      Can you please post example code that triggers the error? I never use type II tests so I do not understand what you are up to. If I can recreate the error I will most likely be able to fix it in the github version.

  • #121

    Elvio Blini
    Participant

    Yes, sorry. Actually, my solution is working for my data and models, but then with built-in data another one seem to pop-out…

    library(afex)
    
    sessionInfo()
     
    # R version 3.4.0 (2017-04-21)
    # Platform: x86_64-w64-mingw32/x64 (64-bit)
    # Running under: Windows >= 8 x64 (build 9200)
    # 
    # Matrix products: default
    # 
    # locale:
    #   [1] LC_COLLATE=English_United Kingdom.1252  LC_CTYPE=English_United Kingdom.1252    LC_MONETARY=English_United Kingdom.1252
    # [4] LC_NUMERIC=C                            LC_TIME=English_United Kingdom.1252    
    # 
    # attached base packages:
    #   [1] stats     graphics  grDevices utils     datasets  methods   base     
    # 
    # other attached packages:
    #   [1] afex_0.18-0      lsmeans_2.26-3   estimability_1.2 lme4_1.1-13      Matrix_1.2-9    
    # 
    # loaded via a namespace (and not attached):
    #   [1] zoo_1.8-0           modeltools_0.2-21   coin_1.2-1          reshape2_1.4.2      splines_3.4.0       lmerTest_2.0-33     lattice_0.20-35    
    # [8] colorspace_1.3-2    htmltools_0.3.6     stats4_3.4.0        mgcv_1.8-17         base64enc_0.1-3     survival_2.41-3     rlang_0.1.1        
    # [15] nloptr_1.0.4        foreign_0.8-69      RColorBrewer_1.1-2  multcomp_1.4-6      plyr_1.8.4          stringr_1.2.0       MatrixModels_0.4-1 
    # [22] munsell_0.4.3       gtable_0.2.0        htmlwidgets_0.9     mvtnorm_1.0-6       coda_0.19-1         codetools_0.2-15    knitr_1.16         
    # [29] latticeExtra_0.6-28 SparseM_1.77        quantreg_5.33       pbkrtest_0.4-7      parallel_3.4.0      htmlTable_1.9       TH.data_1.0-8      
    # [36] Rcpp_0.12.12        xtable_1.8-2        acepack_1.4.1       scales_0.4.1        backports_1.1.0     checkmate_1.8.3     Hmisc_4.0-3        
    # [43] gridExtra_2.2.1     digest_0.6.12       ggplot2_2.2.1       stringi_1.1.5       grid_3.4.0          tools_3.4.0         magrittr_1.5       
    # [50] sandwich_2.4-0      lazyeval_0.2.0      tibble_1.3.3        Formula_1.2-2       cluster_2.0.6       car_2.1-5           MASS_7.3-47        
    # [57] data.table_1.10.4   minqa_1.2.4         rpart_4.1-11        nnet_7.3-12         nlme_3.1-131        compiler_3.4.0 
    
    #example 1
    
    data("sk2011.2")
    
    # use only affirmation problems (S&K also splitted the data like this)
    sk2_aff <- droplevels(sk2011.2[sk2011.2$what == "affirmation",])
    
    # set up model with maximal by-participant random slopes 
    sk_m1 <- mixed(response ~ inference*type+(inference*type|id), sk2_aff, 
                   type= 2, method= "LRT")
    
    #fails
    lsmeans(sk_m1, ~type)
    # 
    # Error in array(seq_len(nrow(RG@linfct)), dims) : 
    #   'dims' cannot be of length 0
    # In addition: Warning message:
    #   In model.matrix.default(trms, m, contrasts.arg = contrasts) :
    #   problem with term 1 in model.matrix: no columns are assigned
    
    #same as before
    lsmeans(sk_m1$full_model[[1]], ~type)
    
    #fails again...
    lsmeans(tail(sk_m1$full_model, 1)[[1]], ~type)
     
    # NOTE: Results may be misleading due to involvement in interactions
    # Error in calculation of the Satterthwaite's approximation. The output of lme4 package is returned
    # Error in base::chol2inv(x, ...) : 'a' must be a numeric matrix
    
  • #122

    henrik
    Keymaster

    Okay, my bug fix did not work for models with interactions. I have provided a new bug fix that should work now, simply download the newest github version.

    The last error (Error in calculation of the Satterthwaite's approximation. The output of lme4 package is returned) however has nothing to do with afex. This is a problem of lsmeans. You need to switch how the df are calculated:

    lsm.options(lmer.df = "Kenward-Roger") # set df for lsmeans to KR
    lsm.options(lmer.df = "Satterthwaite") # the default (does not work in some cases)
    lsm.options(lmer.df = "asymptotic") # the fastest, no df

You must be logged in to reply to this topic.