tests/testthat/test-predict3.R

In merTools: Tools for Analyzing Mixed Effect Regression Models

# Special cases - rank deficiency----

test_that("Prediction intervals are accurate with interaction terms and rank deficiency", {
  skip_on_ci()
  skip_on_cran()
  set.seed(54656)
  n <- 20
  x <- y <- rnorm(n)
  z <- rnorm(n)
  r <- sample(1:5, size=n, replace=TRUE)
  d <- data.frame(x,y,z,r)
  d2 <- expand.grid(a=factor(1:4),b=factor(1:4),rep=1:10)
  n <- nrow(d2)
  d2 <- transform(d2,r=sample(1:5, size=n, replace=TRUE),
                  z=rnorm(n))
  d2 <- subset(d2,!(a=="4" & b=="4"))

  fm <- lmer( z ~ a*b + (1|r), data=d2) |> suppressMessages()
  expect_s3_class(predictInterval(fm, newdata = d2[1:10, ]), "data.frame")

  newPred <- predictInterval(fm, newdata = d2, level = 0.8, n.sims = 1500,
                             stat = 'median', include.resid.var = FALSE,
                             seed = 2342)
  truPred <- predict(fm, newdata = d2)
  expect_equal(mean(newPred$fit - truPred), 0, tolerance = sd(truPred)/15)
  fm2 <- lmer( z ~ a*b + (1+b|r), data=d2) |> suppressMessages()
  #In the call below we are getting warnings because our call to mvtnorm::rmvnorm
  #is shotting a warning when mean and sigma of multivariate distribution are
  #zero using the method="chol
  newPred <- suppressWarnings(predictInterval(fm2, newdata = d2, level = 0.8,
                                              n.sims = 1000,
                                              stat = 'median', include.resid.var = FALSE))
  truPred <- predict(fm2, newdata = d2)
  expect_s3_class(newPred, "data.frame")
  expect_equal(mean(newPred$fit - truPred), 0, tolerance = sd(truPred)/10)
})

# Test the simResults----


test_that("simResults option behaves", {
  skip_on_cran()
  m1 <- lmer(Reaction ~ Days + (1 | Subject), sleepstudy)
  preds1 <- predictInterval(m1, newdata = sleepstudy[1:5, ])
  preds2 <- predictInterval(m1, newdata = sleepstudy[1:5, ],
                            returnSims = TRUE)
  expect_null(attr(preds1, "sim.results"))
  expect_true(is.matrix(attr(preds2, "sim.results")))
  out <- attr(preds2, "sim.results")
  expect_equal(ncol(out), 1000)
  expect_equal(nrow(out), 5)
  m1 <- lmer(Reaction ~ Days + (1 | Subject), sleepstudy)
  preds1 <- predictInterval(m1, newdata = sleepstudy[1:5, ],
                            returnSims = TRUE,
                            which = "random", seed = 23151)
  preds2 <- predictInterval(m1, newdata = sleepstudy[1:5, ], which = "fixed",
                            returnSims = TRUE, seed = 23151)
  preds3 <- predictInterval(m1, newdata = sleepstudy[1:5, ], which = "all",
                            returnSims = TRUE, seed = 23151)
  preds4 <- predictInterval(m1, newdata = sleepstudy[1:5, ],
                            returnSims = TRUE, seed = 23151)
  preds1b <- predictInterval(m1, newdata = sleepstudy[1:5, ],
                             returnSims = TRUE,
                             which = "random", seed = 23151,
                             include.resid.var = FALSE)
  preds2b <- predictInterval(m1, newdata = sleepstudy[1:5, ], which = "fixed",
                             returnSims = TRUE, seed = 23151,
                             include.resid.var = FALSE)
  preds3b <- predictInterval(m1, newdata = sleepstudy[1:5, ], which = "all",
                             returnSims = TRUE, seed = 23151,
                             include.resid.var = FALSE)
  preds4b <- predictInterval(m1, newdata = sleepstudy[1:5, ],
                             returnSims = TRUE, seed = 23151,
                             include.resid.var = FALSE)
  expect_true(is.matrix(attr(preds1, "sim.results")))
  expect_gt(abs(mean(attr(preds1, "sim.results") - attr(preds2, "sim.results"))),
            200)
  expect_true(is.matrix(attr(preds2, "sim.results")))
  expect_gt(abs(mean(attr(preds4, "sim.results") - attr(preds1, "sim.results"))),
            200)
  expect_gt(abs(mean(attr(preds4, "sim.results") - attr(preds1, "sim.results"))),
            abs(mean(attr(preds1, "sim.results") - attr(preds2, "sim.results"))))
  expect_type(attr(preds3, "sim.results"), "list")
  expect_true(is.list(attr(preds3, "sim.results")))
  expect_true(is.matrix(attr(preds1b, "sim.results")))
  expect_true(is.matrix(attr(preds2b, "sim.results")))
  expect_true(is.list(attr(preds3b, "sim.results")))
  expect_true(is.matrix(attr(preds4b, "sim.results")))
  # Check that samples are wider for include.resid.var = TRUE
  expect_gt(quantile(attr(preds1, "sim.results"), probs = 0.9) - quantile(attr(preds1b, "sim.results"), probs = 0.9),
            20)
  expect_lt(quantile(attr(preds1, "sim.results"), probs = 0.1) - quantile(attr(preds1b, "sim.results"), probs = 0.1),
            -20)
  expect_gt(quantile(attr(preds2, "sim.results"), probs = 0.9) - quantile(attr(preds2b, "sim.results"), probs = 0.9),
            20)
  expect_lt(quantile(attr(preds2, "sim.results"), probs = 0.1) - quantile(attr(preds2b, "sim.results"), probs = 0.1),
            -20)
  expect_gt(quantile(attr(preds4, "sim.results"), probs = 0.9) - quantile(attr(preds4b, "sim.results"), probs = 0.9),
            15)
  expect_lt(quantile(attr(preds4, "sim.results"), probs = 0.1) - quantile(attr(preds4b, "sim.results"), probs = 0.1),
            -15)
})

# Test out of sample predictions----

test_that("predictInterval makes predictions without observed outcome", {
  skip_on_ci()
  skip_on_cran()
  possNames <- expand.grid(letters,LETTERS)
  possNames <- paste(possNames[, 1], possNames[, 2])
  newFac <- sample(possNames, 32)
  modData <- data.frame(
    y = rnorm(500),
    x = rnorm(500),
    team_name = sample(newFac, 500, replace = TRUE)
  )
  modData$y[251:500] <- rep(NA, 250)
  m0 <- lmer(y ~ x + (1|team_name), data = modData[1:250,]) |> suppressMessages()
  #In the calls below we are getting warnings because our call to mvtnorm::rmvnorm
  #is shotting a warning when mean and sigma of multivariate distribution are
  #zero using the method="chol
  testPreds1 <- suppressWarnings(predictInterval(m0, newdata = modData[, c(3, 2, 1)]))
  testPreds2 <- suppressWarnings(predictInterval(m0, newdata = modData[1:250, c(2, 3, 1)]))
  testPreds3 <- suppressWarnings(predictInterval(m0, newdata = modData[251:500,]))
  expect_s3_class(testPreds1, "data.frame")
  expect_s3_class(testPreds2, "data.frame")
  expect_s3_class(testPreds3, "data.frame")
})

# Input validation checks----



test_that("dplyr objects are successfully coerced", {
  skip_on_cran()
  set.seed(101)
  data(sleepstudy)
  m1 <- lmer(Reaction ~ Days + (1 | Subject), sleepstudy)
  predData <- sleepstudy |>
    dplyr::group_by(Subject) |>
    dplyr::summarise(Days = mean(Days))
  expect_warning(predictInterval(m1, newdata = predData),
                 regexp = "newdata is tbl_df or tbl object from dplyr package")
  #Suppress the warning that we tested for above
  preds2 <- suppressWarnings(predictInterval(m1, newdata = predData, n.sims=2000))
  expect_s3_class(preds2, "data.frame")
  predData2 <- as.data.frame(predData)
  preds1 <- predictInterval(m1, newdata = predData2, n.sims=2000)
  expect_true(sum(preds1$fit - preds2$fit) > -50 & sum(preds1$fit - preds2$fit) < 50)
})

# Model type warnings for non-binomial GLMM----


test_that("Warnings issued", {
  skip_on_cran()
  d <- expand.grid(fac1=LETTERS[1:5], grp=factor(1:10),
                   obs=1:50)
  d$y <- simulate(~fac1+(1|grp),family = poisson,
                  newdata=d,
                  newparams=list(beta=c(2,-1,3,-2,1.2), theta=c(.33)),
                  seed = 5636)[[1]] |> suppressMessages()
  g1 <- glmer(y~fac1+(1|grp), data=d, family = 'poisson')
  expect_warning(predictInterval(g1, newdata = d[1:100,]),
                 regexp = "Prediction for NLMMs or GLMMs that are not mixed")
  rm(list = ls())
})