tests/testthat/test-profileLikelihood.R

In Cyclops: Cyclic Coordinate Descent for Logistic, Poisson and Survival Analysis

library("testthat")
library("survival")

suppressWarnings(RNGversion("3.5.0"))

test_that("Check evaluate profile likelihood with two variables", {
    test <- read.table(header=T, sep = ",", text = "
start, length, event, x1, x2
0, 4,  1,0,0
0, 3.5,1,2,0
0, 3,  0,0,1
0, 2.5,1,0,1
0, 2,  1,1,1
0, 1.5,0,1,0
0, 1,  1,1,0
")
    data <- createCyclopsData(Surv(length, event) ~ x1 + x2, data = test,
                              modelType = "cox")
    fit <- fitCyclopsModel(data)

    x <- seq(from = 0.8, to = 0.9, length = 100)
    out <- getCyclopsProfileLogLikelihood(fit, "x1", x)

    argMax <- out$point[which(out$value == max(out$value))]
    expect_equivalent(coef(fit)["x1"], argMax, tolerance = 0.01)
})

test_that("Check evaluate profile likelihood with multiple threads", {
    test <- read.table(header=T, sep = ",", text = "
start, length, event, x1, x2
0, 4,  1,0,0
0, 3.5,1,2,0
0, 3,  0,0,1
0, 2.5,1,0,1
0, 2,  1,1,1
0, 1.5,0,1,0
0, 1,  1,1,0
")
    data <- createCyclopsData(Surv(length, event) ~ x1 + x2, data = test,
                              modelType = "cox")
    fit <- fitCyclopsModel(data, control = createControl(threads = 2))

    x <- seq(from = 0.8, to = 0.9, length = 100)
    out <- getCyclopsProfileLogLikelihood(fit, "x1", x)

    argMax <- out$point[which(out$value == max(out$value))]
    expect_equivalent(coef(fit)["x1"], argMax, tolerance = 0.01)
})

test_that("Check evaluate profile likelihood with one variable; cold-start", {
    test <- read.table(header=T, sep = ",", text = "
start, length, event, x1, x2
0, 4,  1,0,0
0, 3.5,1,2,0
0, 3,  0,0,1
0, 2.5,1,0,1
0, 2,  1,1,1
0, 1.5,0,1,0
0, 1,  1,1,0
")
    data <- createCyclopsData(Surv(length, event) ~ x1, data = test,
                              modelType = "cox")
    fit <- fitCyclopsModel(data, fixedCoefficients = c(TRUE))

    x <- seq(from = 0.3, to = 0.4, length = 100)
    out <- getCyclopsProfileLogLikelihood(fit, "x1", x)
    argMax <- out$point[which(out$value == max(out$value))]

    gold <- fitCyclopsModel(data, forceNewObject = TRUE)
    expect_equivalent(coef(gold)["x1"], argMax, tolerance = 0.01)
    expect_equivalent(coef(fit)["x1"], 0)
})

test_that("Check adapative profiling likelihood", {
    test <- read.table(header=T, sep = ",", text = "
start, length, event, x1, x2
0, 4,  1,0,0
0, 3.5,1,2,0
0, 3,  0,0,1
0, 2.5,1,0,1
0, 2,  1,1,1
0, 1.5,0,1,0
0, 1,  1,1,0
")
    data <- createCyclopsData(Surv(length, event) ~ x1 + x2, data = test,
                              modelType = "cox")
    fit <- fitCyclopsModel(data)

    expect_error(getCyclopsProfileLogLikelihood(fit, "x1"))
    expect_error(getCyclopsProfileLogLikelihood(fit, "x1", x = c(0,1), bounds = c(0.8, 0.9)))
    expect_error(getCyclopsProfileLogLikelihood(fit, "x1", bounds = c(0.8, 0.8)))

    out <- getCyclopsProfileLogLikelihood(fit, "x1", bounds = c(0, 2), initialGridSize = 10)
    expect_gt(nrow(out), 10)

    argMax <- out$point[which(out$value == max(out$value))]
    expect_equal(coef(fit)["x1"], argMax)
})

test_that("Check adapative profiling likelihood, other covariate is perfect predictor", {

    test <- read.table(header=T, sep = ",", text = "
start, length, event, x1, x2
0, 4,  1,0,1
0, 4,1,2,1
0, 4,  0,0,0
0, 4,1,0,1
0, 4,  1,1,1
0, 4,0,1,0
0, 4,  1,1,1
")
    data <- createCyclopsData(event ~ x1 + x2, data = test,
                              modelType = "cpr")

    # Recreating combination of circumstances that lead to infinite log likelihood
    # Note: infinite variance would have been selected by cross-validation if we had enough data:
    fit <- fitCyclopsModel(data,
                           prior = createPrior("laplace", variance = Inf, exclude = "x1"))
    expect_warning(getCyclopsProfileLogLikelihood(fit, "x1", bounds = c(0, 2), initialGridSize = 10),
                   "Failing to compute likelihood at entire initial grid")

    # This just produces NaN likelihood:
    fit <- fitCyclopsModel(data)
    expect_warning(getCyclopsProfileLogLikelihood(fit, "x1", bounds = c(0, 2), initialGridSize = 10),
                   "Failing to compute likelihood at entire initial grid")
})