tests/testthat/test-ANOVA-&-lmerTest.R

In spaMM: Mixed-Effect Models, with or without Spatial Random Effects

cat(crayon::yellow("\ntest lmerTest interface and other ANOVA tables:")) 
# See test-rank for additional tests of anova()

{ 
  { # LM (more LMs with sleepstudy data below)
    fit <- lm(sr ~ ., data = LifeCycleSavings)
    spfit <- fitme(sr ~ pop15+pop75+dpi+ddpi , data = LifeCycleSavings)
    testthat::test_that("whether lm() and fitme ML yield identical anova tables",
                        testthat::expect_true(diff(range(na.omit(unlist(anova(fit))-unlist(anova(spfit)))))<1e-10))
    spfitRE <- fitme(sr ~ pop15+pop75+dpi+ddpi , data = LifeCycleSavings, method="REML")
    testthat::test_that("whether fitme ML and fitme REML yield identical anova tables for LM",
                        testthat::expect_true(diff(range(na.omit(unlist(anova(spfit))-unlist(anova(spfitRE)))))<1e-10))
  }
  
  { # GLM
    dat <- data.frame(treatment=gl(3,3), outcome= gl(3,1,9), 
                      counts=c(18,17,15,20,10,20,25,13,12)) # showing data
    glm.D93 <- glm(counts ~ outcome + treatment, family = poisson(), data=dat)
    spglm <- fitme(counts ~ outcome + treatment, family = poisson(), data=dat)
    # no F test for poisson
    testthat::test_that("whether fitme and poisson glm yield identical anova Cp tables",
                        testthat::expect_true(diff(range(na.omit(unlist(anova(glm.D93, test = "Cp"))-unlist(anova(spglm, test = "Cp")))))<1e-10))
    testthat::test_that("whether fitme and poisson glm yield identical anova Chisq tables",
                        testthat::expect_true(diff(range(na.omit(unlist(anova(glm.D93, test = "Chisq"))-unlist(anova(spglm, test = "Chisq")))))<1e-10))
    testthat::test_that("whether fitme and poisson glm yield identical anova Rao [score test] tables",
                        testthat::expect_true(diff(range(na.omit(unlist(anova(glm.D93, test = "Rao"))-unlist(anova(spglm, test = "Rao")))))<1e-5)) # less accurate
    
    clotting <- data.frame(
      u = c(5,10,15,20,30,40,60,80,100),
      lot1 = c(118,58,42,35,27,25,21,19,18),
      lot2 = c(69,35,26,21,18,16,13,12,12))
    fit <- glm(lot1 ~ log(u), data = clotting, family = Gamma)
    spglm <- fitme(lot1 ~ log(u), data = clotting, family = Gamma, method="REML")
    # F test returned for Gamma, but large differences are expected from different phi estimation
    #testthat::test_that("whether fitme and poisson glm yield identical anova F tables",
    #                    testthat::expect_true(diff(range(na.omit(unlist(anova(fit, test = "F"))-unlist(anova(spglm, test = "F")))))<1e-3)) 
    testthat::test_that("whether fitme and Gamma glm yield identical anova Cp tables",
                        testthat::expect_true(diff(range(na.omit(unlist(anova(fit, test = "Cp"))-unlist(anova(spglm, test = "Cp")))))<1e-3)) # some difference -- expected from different phi estimation
    # spglm <- fitme(lot1 ~ log(u), data = clotting, family = Gamma, method="ML") # differs more on Cp
    testthat::test_that("whether fitme and Gamma glm yield identical anova Chisq tables",
                        testthat::expect_true(diff(range(na.omit(unlist(anova(fit, test = "Chisq"))-unlist(anova(spglm, test = "Chisq")))))<1e-10))
    testthat::test_that("whether fitme and Gamma glm yield identical anova Rao [score test] tables",
                        testthat::expect_true(diff(range(na.omit(unlist(anova(fit, test = "Rao"))-unlist(anova(spglm, test = "Rao")))))<1e-5)) 
    testthat::test_that("whether fitme and Gamma glm yield identical drop1 F tables (apart from AICs)",
                        testthat::expect_true(diff(range(na.omit(unlist(drop1(fit, test = "F"))-unlist(drop1(spglm, test = "F")))[-(4:5)]))<1e-3)) 
    testthat::test_that("whether fitme and Gamma glm yield identical drop1 Chisq tables (apart from AICs and scaled dev)",
                        testthat::expect_true(diff(range(na.omit(unlist(drop1(fit, test = "Chisq"))-unlist(drop1(spglm, test = "Chisq")))[-(4:6)]))<1e-3)) 
    testthat::test_that("whether fitme and Gamma glm yield identical drop1 Rao [score test] tables (apart from AICs and and scaled dev)",
                        testthat::expect_true(diff(range(na.omit(unlist(drop1(fit, test = "Rao"))-unlist(drop1(spglm, test = "Rao")))[-(4:6)]))<1e-3)) # some fifference -- expected from different phi estimation
  }

  # LM and LMM
  if (spaMM.getOption("example_maxtime")>3) { # not an actual timing
    if(requireNamespace("lme4", quietly = TRUE)) {
      data("sleepstudy", package="lme4")
      
      { # LM
        lmf <- lm(formula = Reaction ~ Days+I(Days^2), data=sleepstudy)
        glmf <- glm(formula = Reaction ~ Days+I(Days^2), data=sleepstudy)
        splm <- fitme(formula = Reaction ~ Days+I(Days^2), data=sleepstudy, method="REML")
        anova(lmf)
        ano1 <- anova(glmf, test="F")
        ano2 <- anova(splm)
        crit <- diff(range(ano1[2:3,"Pr(>F)"]-ano2[1:2,"Pr(>F)"]))
        testthat::test_that("equivalent anova glm() and fitme()",
                            testthat::expect_true(crit<1e-10))
      }
      
      if(requireNamespace("lmerTest", quietly = TRUE)) { # LMM
        
        spfit <- fitme(Reaction ~ Days + (1|Subject), sleepstudy)
        spfit_lmlt <- as_LMLT(spfit) 
        ano3 <- lmerTest::contest(spfit_lmlt, L=c(1,0))
        spfit_lmlt <- as_LMLT(spfit, transf=FALSE) 
        ano2 <- lmerTest::contest(spfit_lmlt, L=c(1,0))
        fm <- lme4::lmer(Reaction ~ Days + (1|Subject), sleepstudy, REML=FALSE)
        ano1 <- lmerTest::contest(fm, L=c(1,0))
        crit <- abs(1-ano1[,"F value"]/ano2[,"F value"]) # p-value too low for informative comparison; compar of relative F values
        testthat::test_that("contest() results unchanged", testthat::expect_true(crit<1e-6)) 
        
        spfit <- fitme(Reaction ~ Days + (1+Days|Subject), sleepstudy)
        spfit_lmlt <- as_LMLT(spfit) 
        ano3 <- lmerTest::contest(spfit_lmlt, L=c(1,0)) 
        spfit_lmlt <- as_LMLT(spfit, transf=FALSE) 
        ano2 <- lmerTest::contest(spfit_lmlt, L=c(1,0)) 
        fm <- lme4::lmer(Reaction ~ Days + (1+Days|Subject), sleepstudy, REML=FALSE)
        ano1 <- lmerTest::contest(fm, L=c(1,0))
        crit <- diff(range(ano3[,"F value"],ano2[,"F value"], 1436.88833634)) # but slight difference with lmer; already observed with brute refit rather than hlcorcall hack; affected by xtol_rel
        testthat::test_that("equivalent contest() merMod and fitme()",
                            testthat::expect_true(crit<1e-3))
        
        fm <- lme4::lmer(Reaction ~ Days + I(Days^2) + (1|Subject) + (0+Days|Subject), sleepstudy)
        spfit <- fitme(Reaction ~ Days + I(Days^2) + (1|Subject) + (0+Days|Subject), sleepstudy, method="REML")
        spfit_lmlt <- as_LMLT(spfit) 
        ano3 <- anova(spfit_lmlt, type="1")
        ano2 <- anova(spfit_lmlt, type="1", transf=FALSE)
        ano1 <- anova(lmerTest::as_lmerModLmerTest(fm), type="1")
        crit <- diff(range(ano1[,"Pr(>F)"]-ano2[,"Pr(>F)"]))
        testthat::test_that("equivalent type-1 anova merMod and fitme()",
                            testthat::expect_true(crit<1e-6))
        ano3 <- anova(spfit_lmlt, type="2")
        ano2 <- anova(spfit_lmlt, type="2", transf=FALSE)
        ano1 <- anova(lmerTest::as_lmerModLmerTest(fm), type="2")
        crit <- diff(range(ano1[,"Pr(>F)"]-ano2[,"Pr(>F)"]))
        testthat::test_that("equivalent type-2 anova merMod and fitme()",
                            testthat::expect_true(crit<1e-6))
        ano3 <- anova(spfit_lmlt, type="3")
        ano2 <- anova(spfit_lmlt, type="3", transf=FALSE)
        ano1 <- anova(lmerTest::as_lmerModLmerTest(fm), type="3")
        crit <- diff(range(ano1[,"Pr(>F)"]-ano2[,"Pr(>F)"]))
        testthat::test_that("equivalent type-3 anova merMod and fitme()",
                            testthat::expect_true(crit<1e-6))
      } else message("Package 'lmerTest' not available for testing as_LMLT().")
    } 
  }
  
  { # control of input
    data("wafers")
    m1 <- HLfit(y ~X1+X2+X1*X3+X2*X3+I(X2^2)+(1|batch),family=Gamma(log),
                resid.model = ~ X3+I(X3^2) ,data=wafers,method="ML")
    # m2 <- HLfit(y ~X1+X2+X1*X3+X2*X3+I(X2^2)+(1|batch),family=Gamma(log),
    #             resid.model = ~ X3 ,data=wafers,method="ML")
    #  # LRT
    # crit <- diff(range(unlist(anova(m1,m2)$basicLRT)-c(66.33115,  1, 3.330669e-16)))
    # testthat::test_that("anova(m1,m2) gives expected result", testthat::expect_true(crit<1e-5))
    
    m2 <- HLfit(y ~X1+X2+X1*X3+X2*X3+I(X2^2)+(1|batch),family=Gamma(log),
                resid.model = ~ 1 + (1|batch) ,data=wafers,method="ML")
    # detection of mixed-effect phimodel:
    testthat::test_that("anova(m1,m2) gives expected result", 
                        testthat::expect_true(is.null(suppressWarnings(anova(m1,m2)))))

  }
}