Support models with monotonic predictors

[aaelony]

Hello, I have brms model with a number of mo(varname) terms.

It seems that projpred does not support these terms or is it that brms is unaware?

https://github.com/paul-buerkner/brms/blob/master/R/projpred.R#L117-L120

vs_11 <- varsel(brm_model_11, method = "forward")

! Projpred only supports standard multilevel and smoothing terms as well as offsets.
Error in `get_refmodel.brmsfit()`:

Thanks in advance

[fweber144]

Hi @aaelony,

It seems that projpred does not support these terms or is it that brms is unaware?

projpred does not support these terms, sorry.

@aloctavodia: Should I add a label for this (e.g., labeling this as a feature request)?

[aaelony]

Hi @fweber144 ,

Thanks for confirming. Can you sketch out the reasons why these are not supported and what would need to be implemented in projpred (even at a high level)?

Thanks in advance

[fweber144]

I think the bottleneck is to find a maximum-likelihood (or weighted maximum-likelihood, in case of the augmented-data projection) fitter that supports monotonic predictors (and the response families already available in projpred) and then probably also to perform all the necessary adjustments in projpred's codebase.

[aaelony]

Alright, please forgive my audacity of allowing Claude a crack at this.

Claude's response below (if it is helpful)... To be honest, it is a bit crazy that it can come up with this, but I would caution that someone more qualified than myself should get a sense if the proposal is reasonable or at all useful. Feel free to disregard if it appears to be nonsensical.

Proposal: ML Fitter for Monotonic Predictors in projpred

Background: what projpred's fitters do

During projection, projpred's divmin() function iterates over posterior draws (or clusters of draws) and for each one calls a submodel fitter that minimises the KL-divergence from that draw's predictive distribution. Concretely, the fitter receives:

• x: the submodel design matrix,
• y: the reference model's fitted values for that draw (the "pseudo-response"),
• obsvar: the reference model's predictive variance (used as IWLS weights),
• family, weights, etc.

and returns a subfit object with coefficients alpha/beta that can be used to generate predictions. The dispatch lives in divmin() and currently recognises four combinations of has_grp × has_add, routing to fit_glm_ridge_callback, fit_glmer_callback, fit_gam_callback, or fit_gamm_callback.

The problem: mo() terms cannot fit into a standard design matrix because the model is:

η_i = X_i β  +  Σ_j  b_j · D_j · (X^{mo}_{ij})ᵀ ζ_j

where b_j ∈ ℝ is a free scalar, D_j = K_j (number of ordinal levels minus 1), X^{mo}_{ij} ∈ {0,1}^{K_j} is the cumulative indicator row (entry k equals 1 iff observation i has category ≥ k), and ζ_j ∈ Δ^{K_j-1} is a simplex (all entries positive, summing to 1). The simplex constraint means this is not a GLM.

---

The key idea: softmax reparameterisation makes the problem unconstrained

Replace each simplex ζ_j with the softmax of a free vector φ_j ∈ ℝ^{K_j}, fixing one element to zero for identifiability:

ζ_{j,k} = exp(φ_{j,k}) / Σ_{l=1}^{K_j} exp(φ_{j,l}),   with φ_{j,K_j} ≡ 0.

The full parameter vector is then unconstrained and has dimension p + Σ_j K_j, where p is the dimension of the ordinary linear part (including intercept) and K_j is the number of ordinal levels for the j-th monotonic term. Standard L-BFGS-B (via stats::optim) can optimise this.

---

Implementation plan

1. Formula utilities (e.g. in misc.R or a new mo_utils.R)

# Check whether formula contains any mo() terms
has_mo_terms <- function(formula) {
  any(grepl("\\bmo\\s*\\(", deparse(formula)))
}

# Extract info about each mo() term
parse_mo_terms <- function(formula, data) {
  tt <- terms(formula, specials = "mo")
  mo_idx <- attr(tt, "specials")$mo      # positions among term labels
  mo_labels <- attr(tt, "term.labels")[mo_idx]

  lapply(mo_labels, function(lbl) {
    # lbl is something like "mo(income)" or "mo(income, id = \"i\")"
    var_name <- sub("^mo\\(([^,)]+).*\\)$", "\\1", lbl)
    x_raw <- data[[trimws(var_name)]]
    if (is.ordered(x_raw)) {
      cats <- levels(x_raw)
      x_int <- as.integer(x_raw) - 1L   # 0-based
    } else {
      x_int <- as.integer(x_raw)
      cats <- as.character(sort(unique(x_int)))
    }
    K <- length(cats) - 1L              # D in Bürkner & Charpentier notation
    # Build cumulative indicator matrix  N x K
    Xmo <- outer(x_int, seq_len(K), function(xi, k) as.integer(xi >= k))
    list(label = lbl, var = var_name, K = K, Xmo = Xmo)
  })
}

# Build the "plain" design matrix with mo() terms dropped
design_matrix_no_mo <- function(formula, data) {
  # Replace mo(x) with a zero column so model.matrix doesn't error
  f2 <- update(formula, reformulate(
    c("1", setdiff(attr(terms(formula, specials="mo"), "term.labels"),
                   grep("^mo\\(", attr(terms(formula, specials="mo"),
                                       "term.labels"), value=TRUE))),
    response = NULL))
  model.matrix(f2, data)
}

2. The new fitter: fit_mo_callback() (in divergence_minimizers.R)

fit_mo_callback <- function(formula, data, family, weights, y, obsvar,
                             regul = 1e-4, ...) {
  ## ---- Parse -------------------------------------------------------
  mo_info  <- parse_mo_terms(formula, data)          # list of mo term info
  Xreg     <- design_matrix_no_mo(formula, data)     # N x p (includes intercept)
  n        <- nrow(Xreg)
  p        <- ncol(Xreg)
  J        <- length(mo_info)

  ## Indices in the parameter vector
  ## [1..p]  -> beta (ordinary coefficients, including intercept)
  ## then for j = 1..J:
  ##   [p + offset_j + 1]               -> b_j   (scalar)
  ##   [p + offset_j + 2 .. + K_j]      -> phi_j (K_j - 1 free log-ratio params)
  offsets <- c(0, cumsum(sapply(mo_info, function(m) 1L + m$K - 1L)))
  n_par   <- p + offsets[J + 1]

  ## ---- Eta from params ---------------------------------------------
  compute_eta <- function(params) {
    beta <- params[seq_len(p)]
    eta  <- Xreg %*% beta
    for (j in seq_len(J)) {
      K_j    <- mo_info[[j]]$K
      Xmo_j  <- mo_info[[j]]$Xmo
      b_j    <- params[p + offsets[j] + 1L]
      phi_j  <- c(params[p + offsets[j] + seq(2L, K_j)], 0)  # fix last to 0
      zeta_j <- exp(phi_j) / sum(exp(phi_j))                  # softmax -> simplex
      eta    <- eta + b_j * K_j * (Xmo_j %*% zeta_j)
    }
    as.vector(eta)
  }

  ## ---- Objective: IWLS-style weighted deviance (matches glm_ridge logic) --
  ## projpred uses a Gaussian pseudo-likelihood with
  ##   target y_s  (reference model fitted values for draw s)
  ##   weight 1/obsvar_s
  ## The KL-minimising submodel coefficients satisfy exactly this WLS problem
  ## for Gaussian pseudo-response; for GLMs the IRLS iteration does it
  ## iteratively. Here we optimise the full objective directly.
  ##
  ## For the traditional projection the objective is the weighted squared loss
  ## (Gaussian working response):
  ##   L(θ) = Σ_i  w_i / (2 * obsvar_i) * (y_i - η_i)^2  + λ/2 * ||β||^2
  ## where y and obsvar are the per-draw reference values.
  ## For a general family we use the family's deviance/variance function in IRLS.

  neg_obj <- function(params) {
    eta <- compute_eta(params)
    mu  <- family$linkinv(eta)
    ## Working response and weights for IRLS step (same as in glm_ridge):
    ## Here obsvar plays the role of 1/prior-weight * V(mu) * (g'(mu))^2
    resid_sq <- sum(weights * (y - mu)^2 / obsvar)
    ## L2 regularisation on beta (not on b_j or phi_j to keep geometry clean)
    pen <- regul * sum(params[seq_len(p)]^2)
    resid_sq + pen
  }

  ## gradient (optional but speeds up BFGS greatly)
  grad_obj <- function(params) {
    eta  <- compute_eta(params)
    mu   <- family$linkinv(eta)
    dmu_deta <- family$mu.eta(eta)
    ## delta_i = -2 * w_i * (y_i - mu_i) * dmu/deta_i / obsvar_i
    delta <- -2 * weights * (y - mu) * dmu_deta / obsvar
    g <- numeric(n_par)
    ## gradient w.r.t. beta
    g[seq_len(p)] <- t(Xreg) %*% delta + 2 * regul * params[seq_len(p)]
    ## gradient w.r.t. b_j and phi_j
    for (j in seq_len(J)) {
      K_j    <- mo_info[[j]]$K
      Xmo_j  <- mo_info[[j]]$Xmo
      b_j    <- params[p + offsets[j] + 1L]
      phi_j  <- c(params[p + offsets[j] + seq(2L, K_j)], 0)
      zeta_j <- exp(phi_j) / sum(exp(phi_j))
      v_j    <- as.vector(Xmo_j %*% zeta_j)      # N-vector

      ## d/d(b_j)
      g[p + offsets[j] + 1L] <- sum(delta * K_j * v_j)

      ## d/d(phi_{j,k}) for k = 1..K_j-1 (last is fixed to 0)
      ## d(zeta_l)/d(phi_k) = zeta_l * (I(l==k) - zeta_k)
      ## d(v_j)/d(phi_k) = b_j * K_j * Xmo_j %*% d(zeta)/d(phi_k)
      for (k in seq_len(K_j - 1L)) {
        dzeta_dphi_k <- zeta_j * (seq_len(K_j) == k) - zeta_j[k] * zeta_j
        dv_dphi_k    <- as.vector(Xmo_j %*% dzeta_dphi_k)
        g[p + offsets[j] + 1L + k] <- sum(delta * b_j * K_j * dv_dphi_k)
      }
    }
    g
  }

  ## ---- Initialise --------------------------------------------------
  ## Start betas from ridge GLM ignoring mo terms (fast warm start)
  init      <- numeric(n_par)
  ## b_j = 0 -> simplex parameters don't matter; phi_j = 0 -> uniform simplex

  ## ---- Optimise ----------------------------------------------------
  opt <- tryCatch(
    optim(par    = init,
          fn     = neg_obj,
          gr     = grad_obj,
          method = "BFGS",
          control = list(maxit = 500, reltol = 1e-10)),
    error = function(e) {
      warning("fit_mo_callback: optim failed (", conditionMessage(e),
              "); returning zero-coefficient subfit.")
      list(par = init, convergence = 99)
    }
  )
  if (opt$convergence != 0)
    warning("fit_mo_callback: optim did not converge (code ", opt$convergence, ")")

  params_hat <- opt$par
  eta_hat    <- compute_eta(params_hat)
  mu_hat     <- family$linkinv(eta_hat)

  ## ---- Reconstruct named coefficient list --------------------------
  ## Pack into a format analogous to the glm_ridge subfit:
  ##   alpha = intercept,  beta = named vector of remaining coefficients
  beta_reg <- params_hat[seq_len(p)]
  names(beta_reg) <- colnames(Xreg)

  mo_coefs <- lapply(seq_len(J), function(j) {
    K_j   <- mo_info[[j]]$K
    b_j   <- params_hat[p + offsets[j] + 1L]
    phi_j <- c(params_hat[p + offsets[j] + seq(2L, K_j)], 0)
    zeta_j <- exp(phi_j) / sum(exp(phi_j))
    list(b = b_j, zeta = zeta_j, K = K_j, label = mo_info[[j]]$label)
  })

  sub <- list(
    alpha    = beta_reg["(Intercept)"],
    beta     = beta_reg[names(beta_reg) != "(Intercept)"],
    mo_coefs = mo_coefs,          # <-- new field carrying simplex params
    w        = mu_hat,
    formula  = formula,
    x        = Xreg,
    y        = y
  )
  class(sub) <- "subfit"
  return(sub)
}

3. Dispatch in divmin() (in divergence_minimizers.R)

Add a has_mo flag alongside the existing has_grp / has_add checks:

has_mo <- has_mo_terms(formula)

if (!has_grp && !has_add && !has_mo) {
  sdivmin <- get(getOption("projpred.glm_fitter", "fit_glm_ridge_callback"),
                 mode = "function")
} else if (!has_grp && !has_add && has_mo) {
  sdivmin <- fit_mo_callback
} else if (has_grp && !has_add && !has_mo) {
  sdivmin <- fit_glmer_callback
  # ... (existing PQL branch)
} else if (!has_grp && has_add && !has_mo) {
  sdivmin <- fit_gam_callback
} else if (has_grp && has_add && !has_mo) {
  sdivmin <- fit_gamm_callback
} else {
  ## Combinations with mo() + grp or add: could fall back to fit_mo_callback
  ## extended to handle random effects via lme4 (future work), or error out.
  stop("Combinations of mo() with group or additive terms are not yet supported.")
}

4. Prediction from a mo subfit (in predict.subfit / proj_linpred)

The existing predict.subfit uses x %*% beta + alpha. For mo subfits, a small extension is needed:

predict.subfit <- function(object, newdata = NULL, ...) {
  if (!is.null(newdata)) {
    Xreg <- design_matrix_no_mo(object$formula, newdata)
    eta  <- as.vector(Xreg %*% c(object$alpha, object$beta))
    if (!is.null(object$mo_coefs)) {
      mo_info_new <- parse_mo_terms(object$formula, newdata)
      for (j in seq_along(object$mo_coefs)) {
        mc  <- object$mo_coefs[[j]]
        eta <- eta + mc$b * mc$K * as.vector(mo_info_new[[j]]$Xmo %*% mc$zeta)
      }
    }
  } else {
    eta <- as.vector(object$x %*% c(object$alpha, object$beta))
    # mo contribution already baked into object$w via mu_hat
    return(object$w)
  }
  object$family$linkinv(eta)
}

5. Variable selection integration

For projpred's search to work, each mo(x) term must be treated as a single atomic unit — selected or not as a whole — since b and ζ are jointly meaningless when separated. The term label "mo(income)" should appear as-is in the solution path, analogous to how s(x) is treated as one smooth term. This requires verifying that solution_terms() and the term-enumeration logic in varsel() already treats special-function term labels atomically (they do for s(), so the same should hold for mo() without further changes, provided formula parsing uses attr(terms(..., specials="mo"), "term.labels")).

6. Changes to brms/R/projpred.R

Remove (or condition) the guard at lines 117–120 that currently rejects mo() terms:

# BEFORE (brms side):
if (any(sapply(bterms$dpars$mu$fe$sp, function(x) x$type == "mo"))) {
  stop("Projpred only supports standard multilevel and smoothing terms as well as offsets.")
}

# AFTER: remove the mo() branch from this condition,
# since projpred now handles it.

---

Design choices and caveats

Why optim/BFGS rather than IRLS? The simplex constraint makes IRLS inapplicable directly. An alternating approach (fix ζ, update β by IRLS; fix β, project ζ onto simplex) would work but converges slowly near the solution. The unconstrained softmax reparameterisation with BFGS is clean, needs no extra dependencies, and converges quadratically near the optimum. The analytic gradient above makes it fast in practice.

Regularlisation: The L2 penalty on β (ordinary coefficients) mirrors what glm_ridge does. It is reasonable not to penalise b_j or φ_j since these are already implicitly regularised by the projpred objective (the pseudo-response pulls them towards zero when the monotonic term is not needed). Optionally, a small ridge on b_j could be added.

Interactions: mo(x)*z creates an interaction term in brms where a separate simplex is estimated per level of z. This can be handled by treating each combination as a separate mo term in parse_mo_terms(), or by extending the Xmo matrix approach. This is a natural extension for a second PR.

Mixed-effects + mo(): The case has_grp && has_mo would require joint optimisation over random-effect parameters too. A pragmatic first step is to use the estimated BLUP offset from a preliminary lmer fit on the non-mo part, then run fit_mo_callback with that offset fixed — a two-step approximation that is consistent with how projpred already handles GLMMs approximately.

Warm starts: Because projpred calls the fitter once per posterior draw (or cluster), caching the previous draw's opt$par as the next starting point dramatically accelerates convergence and is straightforward to add.

Testing: A minimal test would fit brm(y ~ mo(income), ...), call varsel(), and check that (a) no error is thrown, (b) the recovered mo() coefficient is close to the full-model MAP, and (c) a model with mo(income) scores better than one without it on simulated data where the relationship is genuinely monotonic but non-linear.

---

Summary of files to change

File	Change
R/misc.R (or new R/mo_utils.R)	has_mo_terms(), parse_mo_terms(), design_matrix_no_mo()
R/divergence_minimizers.R	fit_mo_callback() + dispatch in divmin()
R/predict_subfit.R (or wherever predict.subfit lives)	Handle mo_coefs field
R/varsel.R (or solution_terms)	Verify atomic treatment of mo() labels
brms/R/projpred.R	Remove the guard rejecting mo() terms

This is a self-contained, pure-R implementation requiring no new package dependencies — optim is base R. The only truly new mathematical idea is the softmax reparameterisation of the simplex, which converts a constrained quadratic-over-simplex problem into an unconstrained smooth one that BFGS handles well.

[fweber144]

Sorry, I'm currently lacking the time to digest this (I'm spending only a small fraction of my free time on projpred; I'm not working as much on projpred as I used to). But I'm tagging @aloctavodia to see if he has an idea how to deal with this.

[aloctavodia]

I will try to take a look at this next week.