tests/testthat/test-coxph.R

In JointAI: Joint Analysis and Imputation of Incomplete Data

library("JointAI")

# Sys.setenv(IS_CHECK = "true")

skip_on_cran()

if (identical(Sys.getenv("NOT_CRAN"), "true")) {

  PBC2 <- PBC[match(unique(PBC$id), PBC$id), ]
  PBC2$center <- cut(as.numeric(PBC2$id), c(-Inf, seq(30, 270, 30), Inf))
  PBC$center <- cut(as.numeric(PBC$id), c(-Inf, seq(30, 270, 30), Inf))

  PBC2$futime2 <- PBC2$futime
  PBC2$status2 <- PBC2$status
  PBC2$futime2[1:10] <- NA
  PBC2$status2[11:20] <- NA


  run_coxph_models <- function() {

    sink(tempfile())
    on.exit(sink())
    invisible(force(suppressWarnings({

      models <- list(
        # no covariates
        m0a = coxph_imp(Surv(futime, status != "censored") ~ 1, data = PBC2,
                        n.adapt = 1, n.iter = 4, seed = 2020,
                        warn = FALSE, mess = FALSE),

        # only complete
        m1a = coxph_imp(Surv(futime, status != "censored") ~ age + sex,
                        data = PBC2, n.adapt = 2, n.iter = 4, seed = 2020,
                        warn = FALSE, mess = FALSE),
        m1b = coxph_imp(Surv(futime, I(status != "censored")) ~ age + sex,
                        data = PBC2, n.adapt = 2, n.iter = 4, seed = 2020,
                        warn = FALSE, mess = FALSE),

        # only incomplete
        m2a = coxph_imp(Surv(futime, status != "censored") ~ copper,
                        data = PBC2, n.adapt = 2, n.iter = 4, seed = 2020,
                        warn = FALSE, mess = FALSE),


        # complex structures
        m3a = coxph_imp(Surv(futime, status != "censored") ~ copper + sex +
                          age + abs(age - copper) + log(trig),
                        data = PBC2, trunc = list(trig = c(0.0001, NA)),
                        n.adapt = 2, n.iter = 10, seed = 2020,
                        warn = FALSE, mess = FALSE),

        m3b = coxph_imp(Surv(futime, status != "censored") ~ copper + sex +
                          age + abs(age - copper) + log(trig) + (1 | center),
                        data = PBC2, trunc = list(trig = c(0.0001, NA)),
                        n.adapt = 2, n.iter = 10, seed = 2020,
                        warn = FALSE, mess = FALSE),


        # time-varying covariates ---------------------------------------------
        # time-dep variables of different types
        m4a = coxph_imp(Surv(futime, status != "censored") ~ age + sex + trt +
                          albumin + platelet + stage + (1 | id), data = PBC,
                        timevar = "day", n.adapt = 2, n.iter = 10, seed = 2020,
                        warn = FALSE, mess = FALSE),

        # interaction between complete baseline and time-dep variable.
        # Note: interaction with incomplete baseline does not work!
        m4b = coxph_imp(Surv(futime, status != "censored") ~ age + sex * trt +
                          albumin + log(platelet) + (1 | id), data = PBC,
                        timevar = "day",
                        n.adapt = 2, n.iter = 10, seed = 2020,
                        warn = FALSE, mess = FALSE),

        # transformation of incomplete time-dep variable
        m4c = coxph_imp(Surv(futime, status != "censored") ~ age + sex +
                          albumin + log(platelet) + (1 | id) + (1 | center),
                        data = PBC, timevar = "day",
                        n.adapt = 2, n.iter = 10, seed = 2020,
                        warn = FALSE, mess = FALSE),

        m4d = coxph_imp(Surv(futime, status != "censored") ~ age + sex +
                          albumin + ns(platelet, df = 2) + (1 | id) +
                          (1 | center),
                        data = PBC, timevar = "day",
                        n.adapt = 2, n.iter = 10, seed = 2020,
                        warn = FALSE, mess = FALSE)
      )
    }
    )))
    models
  }

  models <- run_coxph_models()
  models0 <- set0_list(models)


  test_that("models run", {
    for (k in seq_along(models)) {
      expect_s3_class(models[[k]], "JointAI")
    }
  })


  test_that("there are no duplicate betas/alphas in the jagsmodel", {
    expect_null(unlist(lapply(models[!grepl("m4[[:alpha:]]+", names(models))],
                              find_dupl_parms)))

    expect_equal(
      find_dupl_parms(models[[7]]),
      paste0("beta[", models[[7]]$par_index_main[[1]]["M_lvlone", 'start'] :
               models[[7]]$par_index_main[[1]]["M_lvlone", 'end'], "]")
    )

    expect_equal(
      find_dupl_parms(models[[8]]),
      paste0("beta[", models[[8]]$par_index_main[[1]]["M_lvlone", 'start'] :
               models[[8]]$par_index_main[[1]]["M_lvlone", 'end'], "]")
    )

    expect_equal(
      find_dupl_parms(models[[9]]),
      paste0("beta[", models[[9]]$par_index_main[[1]]["M_lvlone", 'start'] :
               models[[9]]$par_index_main[[1]]["M_lvlone", 'end'], "]")
    )
  })


  test_that("MCMC is mcmc.list", {
    for (i in seq_along(models)) {
      expect_s3_class(models[[i]]$MCMC, "mcmc.list")
    }
  })


  test_that("MCMC samples can be plottet", {
    for (k in seq_along(models)) {
      expect_silent(traceplot(models[[k]]))
      expect_silent(densplot(models[[k]]))
      expect_silent(plot(MC_error(models[[k]])))
    }
  })

  # test_that("data_list remains the same", {
  #   skip_on_cran()
  #   expect_snapshot(lapply(models, "[[", "data_list"))
  # })

  test_that("jagsmodel remains the same", {
    expect_snapshot(lapply(models, "[[", "jagsmodel"))
  })


  test_that("GRcrit and MCerror give same result", {
    expect_snapshot(lapply(models0, GR_crit, multivariate = FALSE))
    expect_snapshot(lapply(models0, MC_error))
  })



  test_that("summary output remained the same", {
    skip_on_cran()
    expect_snapshot(lapply(models0, print))
    expect_snapshot(lapply(models0, coef))
    expect_snapshot(lapply(models0, confint))
    expect_snapshot(lapply(models0, summary))
    expect_snapshot(lapply(models0, function(x) coef(summary(x))))
  })


  test_that("prediction works", {
    expect_s3_class(predict(models$m3a, type = "lp", warn = FALSE)$fitted,
                    "data.frame")
    expect_s3_class(predict(models$m3b, type = "lp", warn = FALSE)$fitted,
                    "data.frame")
    expect_s3_class(predict(models$m4a, type = "lp", warn = FALSE)$fitted,
                    "data.frame")
    expect_s3_class(predict(models$m4b, type = "lp", warn = FALSE)$fitted,
                    "data.frame")
    expect_s3_class(predict(models$m4c, type = "lp", warn = FALSE)$fitted,
                    "data.frame")

    expect_s3_class(predict(models$m3a, type = "survival", warn = FALSE)$fitted,
                    "data.frame")
    expect_s3_class(predict(models$m3b, type = "survival", warn = FALSE)$fitted,
                    "data.frame")
    expect_s3_class(predict(models$m4a, type = "survival", warn = FALSE)$fitted,
                    "data.frame")
    expect_s3_class(predict(models$m4b, type = "survival", warn = FALSE)$fitted,
                    "data.frame")
    expect_s3_class(predict(models$m4c, type = "survival", warn = FALSE)$fitted,
                    "data.frame")

    expect_s3_class(predict(models$m3a, type = "risk", warn = FALSE)$fitted,
                    "data.frame")
    expect_s3_class(predict(models$m3b, type = "risk", warn = FALSE)$fitted,
                    "data.frame")
    expect_s3_class(predict(models$m4a, type = "risk", warn = FALSE)$fitted,
                    "data.frame")
    expect_s3_class(predict(models$m4b, type = "risk", warn = FALSE)$fitted
                    , "data.frame")
    expect_s3_class(predict(models$m4c, type = "risk", warn = FALSE)$fitted,
                    "data.frame")

    expect_s3_class(predict(models$m3a, type = "expected", warn = FALSE)$fitted,
                    "data.frame")
    expect_s3_class(predict(models$m3b, type = "expected", warn = FALSE)$fitted,
                    "data.frame")
    expect_s3_class(predict(models$m4a, type = "expected", warn = FALSE)$fitted,
                    "data.frame")
    expect_s3_class(predict(models$m4b, type = "expected", warn = FALSE)$fitted,
                    "data.frame")
    expect_s3_class(predict(models$m4c, type = "expected", warn = FALSE)$fitted,
                    "data.frame")
  })



  test_that("residuals", {
    # residuals are not yet implemented
    expect_error(residuals(models$m3b, type = "working", warn = FALSE))
    expect_error(residuals(models$m3b, type = "response", warn = FALSE))
  })


  test_that("model can (not) be plottet", {
    for (i in seq_along(models)) {
      expect_error(plot(models[[i]]))
    }
  })


  test_that("wrong models give errors", {
    # more than two event types
    expect_error(coxph_imp(Surv(futime, status) ~ copper + sex,
                           data = PBC2, n.adapt = 2, n.iter = 4, seed = 2020))

    # missing values in event time
    expect_error(coxph_imp(Surv(futime2, status != "censored") ~ copper + sex,
                           data = PBC2, n.adapt = 2, n.iter = 4, seed = 2020))

    # missing values in event status
    expect_error(coxph_imp(Surv(futime, status2 != "censored") ~ copper + sex,
                           data = PBC2, n.adapt = 2, n.iter = 4, seed = 2020))

    # wrong outcome
    expect_error(coxph_imp(futime ~ copper + sex,
                           data = PBC2, n.adapt = 2, n.iter = 4, seed = 2020))

    # no argument formula
    expect_error(coxph_imp(fixed = futime ~ copper + sex,
                           data = PBC2, n.adapt = 2, n.iter = 4, seed = 2020))

    # time-varying covariates without observations for some subjects
    expect_error(coxph_imp(Surv(futime, status != "censored") ~ age + sex + trt +
                             albumin + chol + (1 | id), data = PBC,
                           timevar = "day",
                           n.adapt = 2, n.iter = 4, seed = 2020
    ))

    # interaction between time-varying covariates and incomplete baseline
    # variable
    expect_error(coxph_imp(Surv(futime, status != "censored") ~ age +
                             trig * platelet + (1 | id), data = PBC,
                           timevar = "day",
                           n.adapt = 2, n.iter = 4, seed = 2020)
    )

  })
}
# Sys.setenv(IS_CHECK = "")