runehaubo / lmerTestR

Repository for the R-package lmerTest
49 stars 9 forks source link

new lmerTest does not support models fitted with old version #5

Closed singmann closed 5 years ago

singmann commented 6 years ago

Hi Rune,

I noticed one more problem when preparing my CRAN submission of afex. The current lmerTest does not support models estimated in the current CRAN version. For example:

install.packages("lmerTest") # current CRAN version
library("lmerTest")
fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy)
anova(fm1)
# Analysis of Variance Table of type III  with  Satterthwaite 
# approximation for degrees of freedom
#      Sum Sq Mean Sq NumDF DenDF F.value    Pr(>F)    
# Days  30031   30031     1    17  45.853 3.264e-06 ***
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
save(fm1, file = "cranversion.rda")

After restarting and installing the new version, this object is essentially useless:

devtools::install_github("runehaubo/lmerTestR")
library("lmerTest")
load("cranversion.rda")
fm1 ## no print method
anova(fm1)
# Error in UseMethod("anova") : 
#   no applicable method for 'anova' applied to an object of class "merModLmerTest"

This seems like a serious shortcoming to me. Sometimes estimating models is quite expensive and working with fitted models is quite convenient. For example, for many of my published papers the fitted models are available to download to reproduce the analysis. With this update they would become pretty useless.

I am not sure how easy it would be to make the appropriate changes as the internals appear to have changed somewhat. but maybe already reverting the decision to change the clas name might work. At least for my real problem, the interplay with emmeans, this appear to work enough.

library(emmeans)
emmeans(fm1, "Days")
# Error in ref_grid(object, ...) : 
#   Can't handle an object of class  “merModLmerTest” 
#  Use help("models", package = "emmeans") for information on supported models.

class(fm1) <- "lmerModLmerTest"
emmeans(fm1, "Days")
#  Days   emmean       SE df lower.CL upper.CL
#   4.5 298.5079 9.049936 17 279.4142 317.6016
# 
# Degrees-of-freedom method: kenward-roger 
# Confidence level used: 0.95 

However, whereas print then also work, anova still does not:

fm1
# Linear mixed model fit by REML ['lmerModLmerTest']
# Formula: Reaction ~ Days + (Days | Subject)
#    Data: sleepstudy
# REML criterion at convergence: 1743.628
# Random effects:
#  Groups   Name        Std.Dev. Corr
#  Subject  (Intercept) 24.740       
#           Days         5.922   0.07
#  Residual             25.592       
# Number of obs: 180, groups:  Subject, 18
# Fixed Effects:
# (Intercept)         Days  
#      251.41        10.47
anova(fm1)
# Error in qform(L, model@vcov_beta) : 
#   no slot of name "vcov_beta" for this object of class "lmerModLmerTest"

Let me know what you think, Henrik

singmann commented 6 years ago

Sorry to bother you again, but I have now added a test to afex to check for this backwards compatibility. The test loads an object created with the current CRAN version of lmerTest and then calls emmeans. With the current CRAN version this works obviously:

library(afex)
load(url("http://singmann.org/download/r/m_machines_lmerTest-pre3.0.rda"))
emmeans(m_machines, "Machine", lmer.df = "asymptotic")
#  Machine   emmean       SE  df asymp.LCL asymp.UCL
#  A       52.35556 2.397324 Inf  47.65689  57.05422
#  B       60.32222 2.633283 Inf  55.16108  65.48336
#  C       66.27222 2.836796 Inf  60.71220  71.83224
# 
# Degrees-of-freedom method: asymptotic 
# Confidence level used: 0.95 

With the current development version this fails again with:

> emmeans(m_machines, "Machine", lmer.df = "asymptotic")
Error in object[[full_model_name]][[1]] : 
  this S4 class is not subsettable
Error in ref_grid(object, ...) : 
  Perhap a 'data' or 'params' argument is needed

Is there any chance I can convince you to perhaps not change the class name from merModLmerTest to lmerModLmerTest? This would at least solve my problems.

Unfortunately, John Fox of the car package has now also contacted me and told me that CRAN will not yet accept the new version of car because it breaks the current afex. For me to update afex I would however really like to now if you change the class name or not. As I said, not changing and preserving some sort of backwards compatibility sounds like the better options in my ears.

runehaubo commented 6 years ago

Hi Henrik,

Hmm, where to start. Basically I think software evolution is a natural thing; you try to do something because you think it is right, but later realise that it wasn't and you have to fix or do it some other way. Preserving backward compatibility is important but sometimes it is less important than fixing a bad construction. R also evolves and so do other packages, so if you want previous analyses to work 'smoothly' then you basically need an 'image' of the software used at the time - there has been several recent JSS articles about such concepts - also providing tools to actually handle it.

That doesn't mean that we shouldn't strive to provide backward compatibility where possible and I am happy to work with you to find a reasonable solution. I have a suggestion for a solution - see below, but at this point I am not convinced that we should change back to merModLmerTest.

Using your example I get the same thing this far:

install.packages("lmerTest") # current CRAN version
library("lmerTest")
packageVersion("lmerTest")
fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy)
anova(fm1)
# Analysis of Variance Table of type III  with  Satterthwaite 
# approximation for degrees of freedom
#      Sum Sq Mean Sq NumDF DenDF F.value    Pr(>F)    
# Days  30031   30031     1    17  45.853 3.264e-06 ***
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
save(fm1, file = "../misc/cranversion.rda")
# restart R
devtools::install_github("runehaubo/lmerTestR")
library("lmerTest")
load("../misc/cranversion.rda")
fm0 <- fm1
# fm1 ## no print method
anova(fm1) # Error

The 'old' fm1 object is essentially an lmerMod object with a different class, that is, there is no actual change to the internal object. We can therefore safely coerce to lmerMod which provides access to all the lme4::merMod methods:

fm1 ## print.merMod
anova(fm1) # lme4::anova.merMod
library(emmeans)
emmeans(fm1, "Days") # works

We can also coerce to lmerModLmerTest to use lmerTest methods:

fm2 <- as_lmerModLmerTest(fm1)
anova(fm2)   # works
summary(fm2)  # works
emmeans(fm2, "Days")  # works

We can also consider adding anova.merModLmerTest and summary.merModLmerTest to directly provide the backward compatibility at a user-level. Something along the following lines might do the trick:

anova.merModLmerTest <- function(object, ..., type = c("III", "II", "I", "3", "2", "1"), 
                                 ddf = c("Satterthwaite", "Kenward-Roger", "lme4")) {
  class(object) <- "lmerMod"
  dots <- list(...)
  models <- if (length(dots)) 
    sapply(dots, is, "lmerModLmerTest") | sapply(dots, is, 
                                                 "merMod") | sapply(dots, is, "lm")
  else logical(0)
  if(any(models)) return(NextMethod())
  df <- match.arg(ddf)
  if (df == "lme4") 
    return(anova(object, ...))

  object <- as_lmerModLmerTest(object)
  anova(object, type=type, ddf=ddf)
}
> anova(fm0)
Type III Analysis of Variance Table with Satterthwaite's method
     Sum Sq Mean Sq NumDF DenDF F value    Pr(>F)    
Days  30031   30031     1    17  45.853 3.264e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Notice that fm0 is a copy of the object made with the old version of lmerTest.

Henrik, let me know what you think about the proposed solutions. Concerning anova and summary methods for the old merModLmerTest objects, I see pros and cons: On the downside it is an extra piece of code that need to be maintained with no future utility and with a fairly straight forward workaround. On the other hand they might be convenient... Feel free to persuade me either way ;-)

The last error you mention,

> library(afex)
> load(url("http://singmann.org/download/r/m_machines_lmerTest-pre3.0.rda"))
> class(m_machines)
[1] "mixed"
> emmeans(m_machines, "Machine", lmer.df = "asymptotic")
Error in object[[full_model_name]][[1]] : 
  this S4 class is not subsettable
Error in ref_grid(object, ...) : 
  Perhaps a 'data' or 'params' argument is needed

comes from emmeans and I don't directly see the connection to lmerTest? It may still root in lmerTest, but at the surface it appears to be related to how emmeans handle mixed objects...

In terms of CRAN submission, my preference would be to have a version of afex which is compatible with the new lmerTest on CRAN rather sooner than later.

Cheers, Rune

singmann commented 6 years ago

Hi Rune,

Thanks for your response. And whereas I totally agree that software needs to move forward at some point, the whole problem here seems to come from the fact that the class name was renamed from merModLmerTest to lmerModLmerTest. Without the renaming there would be no real problem. My issue really is that this renaming breaks the emmeans methods for existing merModLmerTest objects. This is what failed in my last example.

I have now changed my emmeans methods such that it automatically changes the class of any merModLmerTest or lmerModLmerTest to lmerMod (e.g., https://github.com/singmann/afex/blob/master/R/mixed.R#L888). This object is then simply passed on which works nicely. Because I have already received a warning mail from CRAN, I will push this version to CRAN asap. Just running some last checks.

Thanks, Henrik

runehaubo commented 6 years ago

Not sure if you just want to leave this as it is given the resent update to CRAN, but a few comments just in case:

Regarding the merModLmerTest versus lmerModLmerTest issue, why is it not possible to have 'similar' emmeans methods for each class? I would think that an emmeans method for merModLmerTest objects would just have to coerce merModLmerTest -> lmerMod -> lmerModLmerTest and then call the emmeans method for lmerModLmerTest objects??

Also, I wonder why you want to coerce lmerModLmerTest objects to lmerMod - it appears inefficient since the computationally 'heavy' parts of computing the Satterthwaite df will be performed twice. Maybe I'm missing an important point here, but why not coerce the other way from lmerMod to lmerModLmerTest or just leave those two at their respective classes? After all lmerModLmerTest inherits from lmerMod so having just a method for lmerMod should be sufficient and avoids removing the extra information in the 'lmerModLmerTest' objects.

Cheers, Rune

aminadibi commented 6 years ago

I had the same issue, and it was not clear from the error Warning: Error in UseMethod: no applicable method for 'fixef' applied to an object of class "merModLmerTest" that the problem was backward compatibility with previously fitted model objects. I'd say at least make it clear to users that they cannot reuse previously-fitted models.

singmann commented 6 years ago

I still believe it would be a good idea to provide at least some sort of backward compatibility. Either by keeping the support for the old models or at least some relevant methods (print, anova, and the like). Or maybe just an informative error message.

runehaubo commented 6 years ago

Thanks for bringing this up.

I'll work on a solution along the lines of adding anova.merModLmerTest and summary.merModLmerTest as outlined above - what do you think of this? Would this be sufficient?

lmerTest has never defined basic methods such as print, fixef, coef etc. but in this respect only overloaded anova and summary

runehaubo commented 6 years ago

Legacy lmerTest support has been added in the legacy-support branch.

Comments and suggestions are highly appreciate before I merge into master.

Vadishares commented 6 years ago

Hi Rune,

Our team also faced this issue of compatibility, we were not able to get a summary of the model and also not to predict.

Error in UseMethod("predict") : no applicable method for 'predict' applied to an object of class "merModLmerTest".

The simplest workaround solution for this is to change the class of the model like the following, class(model) <- "lmerMod" ( Originally posted in StackOverflow :https://stackoverflow.com/questions/31319030/r-stargazer-lme4-and-lmertest-incompatibility)

Hope that helps.

Thanks Vadiraj

runehaubo commented 6 years ago

Thanks for your comments. As a previous comment mentions, support for legacy fits has been added to the legacy-support branch so I would expect the predict method for merModLmerTest objects to work if you install that branch (I think devtools::install_github("runehaubo/lmerTestR", ref="legacy-support") should do that). Alternatively I suppose you could also just run merModLmerTest <- setClass("merModLmerTest", contains = c("lmerMod")) to define the class to contain the lmerMod class and continue to use the master branch or CRAN version of lmerTest.

Cheers Rune

runehaubo commented 5 years ago

Closed by pull request merge #22