compare_analytical_and_numerical_grad <- function(
  model_name, n_obs = 32, n_pred = 4, n_test = 10, data_seed = 1918, loc_seed = 615
) {
  n_obs <- 32; n_pred <- 4
  data <- simulate_data(n_obs, n_pred, model_name, seed = data_seed)
  design <- data$design; outcome <- data$outcome
  model <- new_regression_model(design, outcome, model_name) 
  loglik_func <- function (coef) { 
    calc_loglik(model, coef)
  }
  grad_func <- function (coef) {
    calc_grad(model, coef)
  }
  set.seed(loc_seed)
  grads_are_close <- TRUE
  for (i in 1:n_test) {
    if (!grads_are_close) break
    regcoef <- rnorm(n_pred)
    analytical_grad <- grad_func(regcoef)
    numerical_grad <- approx_grad_via_finite_diff(loglik_func, regcoef)
    grads_are_close <- are_all_close(
      analytical_grad, numerical_grad, abs_tol = Inf, rel_tol = 1e-3
    )
  }
  return(grads_are_close)
}

test_that("model's analytical gradient is close to numerical one", {
  models <- c("linear", "logit", "poisson")
  for (model_name in models) {
    expect_true(compare_analytical_and_numerical_grad(model_name),
                label = paste0("Testing model: ", model_name))
  }
})