lsmeans compatibility with mixed() and type 2 tests

[Elvio Blini]

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

[henrik]

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.

[Elvio Blini]

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

[Elvio Blini]

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

[henrik]

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.

[Elvio Blini]

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

[henrik]

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

[Author]

Posts