test_that("unregularized logistic regression matches output from glm()", {
  set.seed(1)
  X <- scale(matrix(rnorm(3000), ncol = 3))
  x1 <- X[, 1]
  x2 <- X[, 2]
  x3 <- X[, 3]
  z <- 1 + 2 * x1 + 3 * x2 + x3
  pr <- 1 / (1 + exp(-z))
  y <- rbinom(1000, 1, pr)

  df <- data.frame(y = y, x1 = x1, x2 = x2)
  glm_fit <- glm(y ~ x1 + x2 + x3, data = df, family = "binomial")

  g_model <- SLOPE(cbind(x1, x2, x3), y,
    family = "binomial",
    diagnostics = TRUE,
    alpha = 1e-7
  )

  expect_equivalent(coef(glm_fit),
    coef(g_model),
    tol = 1e-3
  )
})

test_that("regularized slope logistic regression picks out correct features", {
  set.seed(2)
  p <- 10
  n <- 200
  k <- 3

  x <- matrix(rnorm(p * n), n, p)
  # x <- scale(x)
  beta <- double(p)
  nz <- sample(p, k)

  beta[nz] <- 10
  z <- x %*% beta + 1
  prob <- 1 / (1 + exp(-z))

  y <- rbinom(n, 1, prob)

  SLOPE_fit <- SLOPE(x, y, family = "binomial", alpha = 1 / sqrt(n))

  expect_setequal(nz, which(SLOPE_fit$nonzeros))
})