tests/testthat/test-logisregr.R

In lrstat: Power and Sample Size Calculation for Non-Proportional Hazards and Beyond

testthat::test_that("logisregr: freq works", {
  fit1 <- logisregr(ingots, event="NotReady", covariates="Heat*Soak",
                    freq="Freq")

  # expand data into individual subjects
  ingots2 <- ingots[rep(1:nrow(ingots), ingots$Freq),]
  fit2 <- logisregr(ingots2, event="NotReady", covariates="Heat*Soak")

  testthat::expect_equal(fit1$sumstat, fit2$sumstat)
  testthat::expect_equal(fit1$parest, fit2$parest)
})


testthat::test_that("logisregr: firth with plci", {
  fit1 <- logisregr(sexagg, event="case",
                    covariates=c("age", "oc", "vic", "vicl", "vis", "dia"),
                    freq="COUNT", firth=TRUE, plci=TRUE)

  # the following values are almost identical to those from SAS PROC LOGISTIC
  # with numerical differences due to different convergence criteria used
  coef = c(0.1203, -1.1060, -0.0688, 2.2689, -2.1114, -0.7883, 3.0959)
  ci.lower = c(-0.8186, -1.9738, -0.9414, 1.2730, -3.2609, -1.6081, 0.7746)
  ci.upper = c(1.0732, -0.3074, 0.7892, 3.4354, -1.1177, 0.0152, 8.0303)

  testthat::expect_equal(round(fit1$parest$beta, 4), coef)
  testthat::expect_equal(round(fit1$parest$lower, 4), ci.lower)
  testthat::expect_equal(round(fit1$parest$upper, 4), ci.upper)
})


testthat::test_that("logisregr: firth with flic", {
  fit1 <- logisregr(sexagg, event="case",
                    covariates=c("age", "oc", "vic", "vicl", "vis", "dia"),
                    freq="COUNT", firth=TRUE, flic=TRUE)

  # the following value is from the flic function in the logistf package
  # of note, the logistf function with the flic option yields a different
  # value, and the variance of the intercept is off in the logistf package
  intercept = 0.1315

  testthat::expect_equal(round(fit1$parest$beta[1], 4), 0.1315)
})


testthat::test_that("logisregr: robust variance", {
  options(contrasts = c("contr.SAS", "contr.poly"))

  fit1 <- logisregr(six, event="wheeze",
                    covariates=c("city", "age", "smoke"),
                    id="case", robust=TRUE)

  # the following values can be obtained from SAS PROC GENMOD statements:
  #   proc genmod data=six;
  #     class case city;
  #     model wheeze(event='1') = city age smoke / dist=bin;
  #     repeated subject=case / type=ind;
  #   run;
  # here the standard error is based on the robust variance estimate
  coef = c(1.2597, 0.1391, -0.2003, -0.1284)
  stderr = c(3.0645, 0.6859, 0.2820, 0.3926)

  testthat::expect_equal(round(fit1$parest$beta, 4), coef)
  testthat::expect_equal(round(fit1$parest$sebeta, 4), stderr)
})