if(require("lme4")){
  library(distfreereg)
  all.equal.distfreereg <- distfreereg:::all.equal.distfreereg
  test_dfr_functions <- distfreereg:::test_dfr_functions
  
  n <- 1e2
  set.seed(20250516)
  df_lmerMod <- data.frame(x = rnorm(n),
                           z = rnorm(n),
                           state = factor(sample(c("CA", "OR", "WA"), size = n,
                                                 replace = TRUE, prob = c(2,1,1))),
                           party = factor(sample(c("D", "R"), size = n, replace = TRUE)),
                           g = rep(1:10, 10))
  state_num <- ifelse(df_lmerMod$state == "CA", 7, ifelse(df_lmerMod$state == "OR", 2, 3))
  party_num <- ifelse(df_lmerMod$party == "D", 1, -1)
  theta <- c(2,1)
  df_lmerMod$y <- theta[1]*df_lmerMod$x + theta[2]*df_lmerMod$z + state_num + rnorm(n)
  
  m <- lmer(y ~ x + (1|state), data = df_lmerMod)
  
  set.seed(20250516)
  dfr_lmerMod <- distfreereg(m, verbose = FALSE, control = list(return_on_error = FALSE))
  print(dfr_lmerMod)
  
  set.seed(20250516)
  dfr_lmerMod_J <- distfreereg(test_mean = m,
                               override = list(J = dfr_lmerMod[["J"]]),
                               control = list(return_on_error = FALSE))
  
  set.seed(20250516)
  dfr_lmerMod_fitted <- distfreereg(test_mean = m,
                                    override = list(fitted_values = dfr_lmerMod[["fitted_values"]]),
                                    control = list(return_on_error = FALSE))
  
  stopifnot(all.equal(dfr_lmerMod, dfr_lmerMod_J))
  stopifnot(all.equal(dfr_lmerMod, dfr_lmerMod_fitted))
  
  set.seed(20250516)
  newdata_lmerMod <- data.frame(y = rnorm(10), x = rnorm(10), z = rnorm(10))
  test_dfr_functions(dfr_lmerMod, newdata = newdata_lmerMod)
  print(signif(predict(dfr_lmerMod, newdata = newdata_lmerMod, re.form = NA), 3))
  
  set.seed(20250516)
  dfr_form_lmerMod <- distfreereg(y ~ x + (1|state), data = df_lmerMod,
                                  method = "lmer", verbose = FALSE)
  
  print(dfr_form_lmerMod)
  
  set.seed(20250516)
  dfr_form_lmerMod_verbose <- distfreereg(y ~ x + (1|state), data = df_lmerMod,
                                          method = "lmer")
  
  stopifnot(all.equal(dfr_lmerMod, dfr_form_lmerMod))
  
  cdfr_form_lmerMod <- asymptotics(dfr_form_lmerMod, reps = 5)
  cdfr_lmerMod <- asymptotics(dfr_lmerMod, reps = 5)
  
  print(signif(rejection(cdfr_form_lmerMod, alpha = c(0.1, 0.5))[,2:3], digits = 3))
  print(signif(rejection(cdfr_lmerMod, alpha = c(0.1, 0.5))[,2:3], digits = 3))
  
  # Orderings
  
  set.seed(20250516)
  dfr_lmerMod_asis <- update(dfr_lmerMod, ordering = "asis")
  set.seed(20250516)
  dfr_form_lmerMod_asis <- update(dfr_form_lmerMod, ordering = "asis")
  stopifnot(all.equal(dfr_lmerMod_asis, dfr_form_lmerMod_asis))
  
  set.seed(20250516)
  dfr_lmerMod_optimal <- update(dfr_lmerMod, ordering = "optimal")
  set.seed(20250516)
  dfr_form_lmerMod_optimal <- update(dfr_form_lmerMod, ordering = "optimal")
  stopifnot(all.equal(dfr_lmerMod_optimal, dfr_form_lmerMod_optimal))
  
  set.seed(20250516)
  dfr_lmerMod_natural <- update(dfr_lmerMod, ordering = "natural")
  set.seed(20250516)
  dfr_form_lmerMod_natural <- update(dfr_form_lmerMod, ordering = "natural")
  stopifnot(all.equal(dfr_lmerMod_natural, dfr_form_lmerMod_natural))
  
  set.seed(20250516)
  dfr_lmerMod_g <- update(dfr_lmerMod, ordering = list("g"))
  set.seed(20250516)
  dfr_form_lmerMod_g <- update(dfr_form_lmerMod, ordering = list("g"))
  stopifnot(all.equal(dfr_lmerMod_g, dfr_form_lmerMod_g))
  
  print(df_lmerMod[dfr_lmerMod_g[["res_order"]],][["g"]])
  
  set.seed(20250516)
  dfr_lmerMod_g_grouped <- update(dfr_lmerMod_g, group = TRUE)
  set.seed(20250516)
  dfr_form_lmerMod_g_grouped <- update(dfr_form_lmerMod_g, group = TRUE)
  stopifnot(all.equal(dfr_lmerMod_g_grouped, dfr_form_lmerMod_g_grouped))
  
  
  
  
  ### Partial output
  
  dfr_lmerMod_partial <- distfreereg(test_mean = m, verbose = FALSE,
                                     control = list(orth_tol = 1e-100))
  names(dfr_lmerMod_partial)
}