Fix novel level value in step_lencode_glm()

NEWS.md

@@ -1,5 +1,7 @@
 # embed (development version)
 
+* Fixed bug where `step_lencode_glm()` calculated the value for novel levels from the mean of the per-level coefficients instead of a pooled model fit to the whole outcome. (#243)
+
 # embed 1.2.2
 
 * Fixed bug on step_umap() where the number of calculated components would be zero. (#271)

R/lencode_glm.R

@@ -40,7 +40,8 @@
 #' factor outcomes are used, the log-odds reflect the event of interest being
 #' the _first_ level of the factor.
 #'
-#' For novel levels, a slightly timmed average of the coefficients is returned.
+#' For novel levels, the coefficient of an intercept-only model fit to the
+#' whole outcome is returned.
 #'
 #' # Tidying
 #'
@@ -168,10 +169,12 @@ glm_coefs <- function(x, y, wts = NULL, ...) {
     x <- as_tibble(x)
   }
 
+  dat <- vec_cbind(x, y)
+
   mod <-
     glm(
       form,
-      data = vec_cbind(x, y),
+      data = dat,
       family = fam,
       weights = wts,
       na.action = na.omit,
@@ -182,7 +185,20 @@ glm_coefs <- function(x, y, wts = NULL, ...) {
   names(coefs) <- gsub("^value", "", names(coefs))
   mean_coef <- mean(coefs, na.rm = TRUE, trim = .1)
   coefs[is.na(coefs)] <- mean_coef
-  coefs <- c(coefs, ..new = mean_coef)
+
+  # For unseen levels, use the coefficient of a pooled model (i.e. the
+  # overall outcome mean on the link scale) rather than a trimmed mean of
+  # the per-level coefficients. See #243.
+  new_form <- as.formula(paste0(names(y), "~ 1"))
+  new_mod <- glm(
+    new_form,
+    data = dat,
+    family = fam,
+    weights = wts,
+    na.action = na.omit,
+    ...
+  )
+  coefs <- c(coefs, ..new = unname(coef(new_mod)))
   if (is.factor(y[[1]])) {
     coefs <- -coefs
   }

tests/testthat/test-lencode_glm.R

@@ -238,6 +238,32 @@ test_that("bad args", {
   )
 })
 
+test_that("new level value is based on a pooled model, not per-level coefs (#243)", {
+  reg_test <- recipe(x1 ~ ., data = ex_dat) |>
+    step_lencode_glm(x3, outcome = vars(x1)) |>
+    prep(training = ex_dat, retain = TRUE)
+
+  ref_mod <- glm(x1 ~ 1, data = ex_dat, family = gaussian)
+
+  key <- reg_test$steps[[1]]$mapping$x3
+  expect_equal(
+    key$..value[key$..level == "..new"],
+    unname(coef(ref_mod))
+  )
+
+  class_test <- recipe(x2 ~ ., data = ex_dat) |>
+    step_lencode_glm(x3, outcome = vars(x2)) |>
+    prep(training = ex_dat, retain = TRUE)
+
+  ref_mod_cls <- glm(x2 ~ 1, data = ex_dat, family = binomial)
+
+  key_cls <- class_test$steps[[1]]$mapping$x3
+  expect_equal(
+    key_cls$..value[key_cls$..level == "..new"],
+    -unname(coef(ref_mod_cls))
+  )
+})
+
 test_that("case weights", {
   wts_int <- rep(c(0, 1), times = c(100, 400))