tests/testthat/test-tidy-classes.R

In latents: Multilevel Latent Profile, Class and Transition Models

## Every result class this package returns must answer print(), summary(),
## plot() and as.data.frame(). These tests cover the classes produced by
## multilpa() and enumerate_classes(), plus the bootstrap comparison, and the
## summary classes they return in turn.
##
## Methods added in this sweep reach the generic only after the NAMESPACE has
## been regenerated by document(). Until then they are exercised by explicit
## call; the last block asserts dispatch itself and skips while registration is
## outstanding.

.tidy_covariate_fit <- function() {
  set.seed(5)
  school <- rep(seq_len(16), each = 8)
  high_class <- rep(rep(c(FALSE, TRUE), length.out = 16), each = 8)
  x <- rnorm(128)
  profile <- ifelse(runif(128) < plogis(-1 + 2 * high_class + 0.8 * x), 2L, 1L)
  data <- data.frame(
    school = school, x = x,
    y1 = rnorm(128, ifelse(profile == 2L, 2, -2), 0.7),
    y2 = rnorm(128, ifelse(profile == 2L, 1.5, -1.5), 0.7))
  list(data = data,
       fit = multilpa(data, c("y1", "y2"), "school", n_profiles = 2,
                            n_group_classes = 2, profile_covariates = "x",
                            n_starts = 2, seed = 1))
}

.tidy_enumeration <- function() {
  set.seed(7)
  data <- data.frame(school = rep(seq_len(12), each = 10),
                     score_a = rnorm(120), score_b = rnorm(120))
  enumerate_classes(data, c("score_a", "score_b"), "school", n_profiles = 1:3,
                    n_group_classes = 1, model = NULL, n_starts = 2,
                    seed = 1)
}

test_that("every result class of this sweep defines all four verbs", {
  classes <- c("multilpa_covariates", "multilpa_enumeration",
               "multilpa_bootstrap_lrt")
  verbs <- c("print", "summary", "plot", "as.data.frame")
  defined <- expand.grid(verb = verbs, class = classes,
                         stringsAsFactors = FALSE)
  defined$method <- paste(defined$verb, defined$class, sep = ".")
  found <- vapply(defined$method, exists, logical(1),
                  envir = asNamespace("latents"), inherits = FALSE)
  expect_true(all(found), info = paste("missing:",
    paste(subset(defined, !found)$method, collapse = ", ")))
})

test_that("a covariate summary is a classed, tidy, fully named object", {
  fixture <- .tidy_covariate_fit()
  summary_object <- summary.multilpa_covariates(fixture$fit)
  expect_s3_class(summary_object, "summary_multilpa_covariates")
  # summary.multilpa() builds its result by subsetting the fit with a vector of
  # field names; a name the fit does not carry silently becomes a NULL element
  # keyed NA. Building field by field cannot do that, and this asserts it.
  expect_false(anyNA(names(summary_object)))
  expect_output(print.summary_multilpa_covariates(summary_object),
                "Multilevel LPA with covariates")
  model <- get_results(summary_object, "model")
  expect_s3_class(model, "data.frame")
  expect_identical(nrow(model), 1L)
  expect_identical(model$n_observations, fixture$fit$n_observations)
  expect_identical(model$n_profiles, fixture$fit$n_profiles)
  expect_equal(model$log_likelihood, fixture$fit$log_likelihood)
  expect_equal(model$bic_groups, fixture$fit$bic)
  coefficients <- get_results(summary_object, "coefficients")
  expect_identical(names(coefficients),
                   c("level", "outcome", "term", "parameter", "estimate"))
  expect_identical(nrow(get_results(summary_object, "starts")), 2L)
  expect_identical(nrow(get_results(summary_object, "profiles")), 4L)
  expect_error(get_results(summary_object, "nonsense"),
               class = "latents_bad_argument")
})

test_that("a covariate fit carries the data its inference verbs ask for", {
  fixture <- .tidy_covariate_fit()
  rebuilt <- .multilpa_cov_stored_data(fixture$fit)
  expect_s3_class(rebuilt, "data.frame")
  expect_identical(nrow(rebuilt), fixture$fit$n_observations)
  expect_true(all(c("y1", "y2", "school", "x") %in% names(rebuilt)))
  # Bit-for-bit, not merely close: the object stores everything the observed
  # information is built from, so `data` is redundant for this fit.
  expect_identical(vcov(fixture$fit, rebuilt), vcov(fixture$fit, fixture$data))
  expect_identical(parameter_inference(fixture$fit, rebuilt),
                   parameter_inference(fixture$fit, fixture$data))
})

test_that("an enumeration summary names the minimising candidate per criterion", {
  candidates <- .tidy_enumeration()
  summary_object <- summary.multilpa_enumeration(candidates)
  expect_s3_class(summary_object, "summary_multilpa_enumeration")
  expect_false(anyNA(names(summary_object)))
  expect_output(print.summary_multilpa_enumeration(summary_object),
                "-- criteria", fixed = TRUE)
  criteria <- get_results(summary_object, "criteria")
  expect_identical(names(criteria), c("criterion", "convention", "n_profiles",
                                      "n_group_classes", "model", "value"))
  expect_identical(nrow(criteria), 14L)
  expect_setequal(unique(criteria$criterion),
                  c("aic", "kic", "bic", "sabic", "caic", "awe", "icl", "clc"))
  grid <- as.data.frame(candidates)
  # The reported value must be the smallest value of that column among the
  # converged candidates, and must sit on the candidate it names.
  column <- ifelse(is.na(criteria$convention), criteria$criterion,
                   paste(criteria$criterion,
                         sub("individuals", "individual", criteria$convention),
                         sep = "_"))
  smallest <- vapply(column, function(name) {
    min(subset(grid, converged)[[name]])
  }, numeric(1))
  expect_equal(criteria$value, unname(smallest))
  named <- vapply(seq_len(nrow(criteria)), function(row) {
    candidate <- subset(grid, n_profiles == criteria$n_profiles[row] &
                          n_group_classes == criteria$n_group_classes[row] &
                          model == criteria$model[row])
    candidate[[column[row]]]
  }, numeric(1))
  expect_equal(criteria$value, named)
  expect_identical(nrow(get_results(summary_object, "candidates")), nrow(grid))
  expect_identical(as.data.frame(summary_object), grid)

  # The same class counts under different covariance models identify two
  # candidates. Criterion minima and the printed disagreement count retain
  # that distinction.
  competing <- candidates
  competing$table <- rbind(grid[1L, , drop = FALSE], grid[1L, , drop = FALSE])
  competing$table$model <- c("EII", "VII")
  competing$table$converged <- TRUE
  competing$table$aic <- c(0, 1)
  competing$table$kic <- c(1, 0)
  competing$fits <- rep(candidates$fits[1L], 2L)
  compared <- summary.multilpa_enumeration(competing)
  winners <- get_results(compared, "criteria")
  expect_identical(winners$model[winners$criterion == "aic"], "EII")
  expect_identical(winners$model[winners$criterion == "kic"], "VII")
  expect_output(print.summary_multilpa_enumeration(compared),
                "2 distinct candidate(s)", fixed = TRUE)
})

test_that("a candidate fit is reached by its class counts, never by position", {
  candidates <- .tidy_enumeration()
  chosen <- candidate_fit(candidates, n_profiles = 2, n_group_classes = 1)
  expect_s3_class(chosen, "multilpa")
  expect_identical(chosen$n_profiles, 2L)
  expect_identical(chosen$n_group_classes, 1L)
  expect_equal(chosen$bic, subset(as.data.frame(candidates),
                                  n_profiles == 2L)$bic_groups)
  expect_error(candidate_fit(candidates, n_profiles = 9),
               class = "latents_unknown_candidate")
  expect_error(candidate_fit(candidates, n_profiles = "two"),
               "single positive integer")
  expect_error(candidate_fit(candidates, n_profiles = 1.5),
               "single positive integer")
  expect_error(candidate_fit(candidates, n_profiles = 2,
                             n_group_classes = 1.5),
               "single positive integer")
})

test_that("a bootstrap comparison is classed and plots its simulated null", {
  skip_on_cran()
  set.seed(53)
  data <- data.frame(g = rep(1:30, each = 8), y = rnorm(240))
  small <- multilpa(data, "y", "g", 1, 1, variance_model = "equal",
                    n_starts = 1, seed = 8)
  large <- multilpa(data, "y", "g", 2, 1, variance_model = "equal",
                    n_starts = 1, seed = 8, max_iter = 3000, tol = 1e-7)
  result <- bootstrap_lrt(small, large, data, iter = 3, n_starts = 1,
                          max_iter = 3000, tol = 1e-7, seed = 42)
  expect_s3_class(result, "multilpa_bootstrap_lrt")
  summary_object <- summary.multilpa_bootstrap_lrt(result)
  expect_s3_class(summary_object, "summary_multilpa_bootstrap_lrt")
  expect_false(anyNA(names(summary_object)))
  test <- get_results(summary_object, "test")
  expect_identical(nrow(test), 1L)
  expect_identical(test$null_profiles, 1L)
  expect_identical(test$alternative_profiles, 2L)
  expect_identical(nrow(get_results(summary_object, "replicates")), 3L)
  file <- tempfile(fileext = ".pdf")
  on.exit(unlink(file), add = TRUE)
  grDevices::pdf(file)
  on.exit(grDevices::dev.off(), add = TRUE, after = FALSE)
  expect_s3_class(plot.multilpa_bootstrap_lrt(result), "multilpa_bootstrap_lrt")
})

test_that("the new methods dispatch through their generics", {
  fixture <- .tidy_covariate_fit()
  # Registration happens when document() rewrites NAMESPACE. Until the
  # integrator has run it, dispatch cannot reach a method added in this sweep,
  # and the explicit calls above are what verify them.
  skip_if_not(inherits(summary(fixture$fit), "summary_multilpa_covariates"),
              "S3 methods are registered by devtools::document(); run it first.")
  candidates <- .tidy_enumeration()
  expect_s3_class(summary(candidates), "summary_multilpa_enumeration")
  expect_s3_class(as.data.frame(summary(fixture$fit)), "data.frame")
  expect_s3_class(as.data.frame(summary(candidates)), "data.frame")
  skip_if_not_installed("ggplot2")
  expect_plot(plot(fixture$fit))
})