Monotone non-positive function of an integer covariate in inlabru

[Brian]

Hi,

I want to fit a smooth, monotone-decreasing, non-positive function of an integer covariate in inlabru. The simplified linear predictor is:

logit(p_it) = intercept + covariates_t + psi(s_i) + f(dose_it)

where psi(s_i) is an SPDE spatial random field and dose_it is the cumulative number of integer treatment rounds delivered to location i prior to year t. f() must satisfy:

f(0) = 0,   f(m) <= f(m-1) <= 0   for all m >= 1 

f() is not assumed linear because we expect diminishing returns such that the drop
from 0 to 1 round is large, but additional rounds yield progressively smaller reductions. I am thinking of:

(A) "clinear" on neg_dose = -dose_it with range = c(0, Inf): 

 f_dose(neg_dose, model = "clinear", range = c(0, Inf))

This gives f(m) = alpha * (-m), alpha >= 0, so f(m) <= 0 and f(0) = 0 are satisfied. But equal gaps are forced and no diminishing returns.

(B) RW1 on integer dose with constr = TRUE: 

f_dose(dose, model = "rw1", scale.model = TRUE, constr = TRUE, ...)

This gives a smooth flexible shape but does NOT enforce monotonicity (can wiggle in sparse-data settings) and does NOT anchor at f(0) = 0 (sum-to-zero constraint floats the level). 

Any guidance would be greatly appreciated.

[Olivier Supplisson]

Hi, 

I think you will find some answers in this related thread: https://groups.google.com/g/r-inla-discussion-group/c/BYu22Wol_gw/m/MBMesNT0AwAJ

Best,

Olivier

[Brian]

Thanks very much for the pointer. That thread was helpful, and I got that strict ordering on the latent field is not feasible directly but that reparameterising via non-negative increments is a workable route. 

Building on that, I have implemented a version using bm_marginal that appears to work, and I would appreciate confirmation if it is a sensible use of the feature. The idea is, I want to fit a monotone-decreasing, non-positive function of an integer covariate. Instead of trying to estimate the function values f(0), f(1), f(2),... directly and constraining them to be ordered, we model the increments: delta_m = f(m) - f(m-1) directly, each with a negative-Exponential prior so each step is non-positive by construction: delta_m ~ -Exponential(lambda) for m = 1, ..., M. We then recover f via cumulative summation: f(m) = sum_{j=1}^m delta_j, with f(0) = 0 trivially.

Each delta_m is implemented as a single scalar component with N(0,1) latent and a marginal mapper applying the negative-Exponential PIT:

```

delta_m(
  1,
  model = "iid",
  hyper = list(prec = list(initial = 0, fixed = TRUE)),
  marginal = bru_mapper_marginal(qfun = qneg_exp, pfun = pneg_exp, rate = 1)
)

```

The cumulative sum is built via indicator columns ind_j = I(dose >= j) so the predictor is:

formula = y ~ Intercept + other_covariates + delta_1*ind_1 + delta_2*ind_2 + ... + delta_M*ind_M

This means each delta_j contributes to the sum only when the observation's dose is at least j, naturally producing f(dose_i).

The code:

```

n <- 500
dose_levels <- 0:6
M <- max(dose_levels)
f_true <- c(0, -0.8, -1.4, -1.8, -2.0, -2.1, -2.15)
dose <- sample(
  dose_levels,
  n,
  replace = TRUE,
  prob = c(0.25, 0.20, 0.18, 0.15, 0.10, 0.07, 0.05)
)
N <- sample(40:60, n, replace = TRUE)
eta_true <- -1.2 + f_true[dose + 1]
y <- rbinom(n, size = N, prob = plogis(eta_true))
dt <- data.table(y = y, N = N, dose = dose)

for (j in 1:M) dt[, paste0("ind_", j) := as.integer(dose >= j)]

qneg_exp <- \(p, rate = 1, lower.tail = TRUE, log.p = FALSE) {
  -qexp(p, rate = rate, lower.tail = !lower.tail, log.p = log.p)
}
pneg_exp <- \(q, rate = 1, lower.tail = TRUE, log.p = FALSE) {
  pexp(-q, rate = rate, lower.tail = !lower.tail, log.p = log.p)
}

u_prior <- list(prec = list(initial = 0, fixed = TRUE))
cmp_terms <- c(
  "Intercept(1)",
  map_chr(1:M, \(j) {
    sprintf(
      "delta_%d(1, model = 'iid', hyper = u_prior, marginal = bru_mapper_marginal(qfun = qneg_exp, pfun = pneg_exp, rate = 0.5))",
      j)
    })
)
   
cmp <- as.formula(paste("~ -1 +", paste(cmp_terms, collapse = " + ")))

form_terms <- c("Intercept", map_chr(1:M, \(j) sprintf("delta_%d * ind_%d", j, j)))
form <- as.formula(paste("y ~", paste(form_terms, collapse = " + ")))

lik <- bru_obs(family = "betabinomial", formula = form, data = dt, Ntrials = dt$N)
fit <- bru(cmp, lik)

# Predicting to compare with the truth
pred_grid <- data.table(dose = dose_levels)
for (j in 1:M) pred_grid[, paste0("ind_", j) := as.integer(dose >= j)]
pred <- predict(
  fit,
  pred_grid,
  as.formula(paste("~", paste( map_chr(1:M, \(j) sprintf("delta_%d * ind_%d", j, j)), collapse = " + ")))
  )
data.frame(dose = dose_levels, truth = f_true, fitted = pred$mean)

  dose truth     fitted
1    0  0.00  0.0000000
2    1 -0.80 -0.7602972
3    2 -1.40 -1.4090189
4    3 -1.80 -1.7611553
5    4 -2.00 -2.1433767
6    5 -2.10 -2.2328596
7    6 -2.15 -2.5400197

```

The fitted curve is monotone non-positive by construction, anchored at f(0) = 0 exactly, and recovers the simulated truth closely. I wanted to enquire if this a sensible use of bm_marginal and whether the cumulative-sum-via-indicators construction has no hidden issue with M separate scalar components that share no information across increments.

Thank you

[Finn Lindgren]

Hi Brian,

Yes, that approach makes sense!

Note: the accumulation of the delta*ind values in your construction is a very similar problem to one recently discussed on the inlabru github discussion page for adding multiple “grouped” fields. I’ll take your use case in mind when working on improving that other use case, as they are similar but not identical.

Finn

On 5 May 2026, at 15:53, Brian <brian...@gmail.com> wrote:

Thanks very much for the pointer. That thread was helpful, and I got that strict ordering on the latent field is not feasible directly but that reparameterising via non-negative increments is a workable route. 

--
You received this message because you are subscribed to the Google Groups "R-inla discussion group" group.
To unsubscribe from this group and stop receiving emails from it, send an email to r-inla-discussion...@googlegroups.com.
To view this discussion, visit https://groups.google.com/d/msgid/r-inla-discussion-group/ddb0314c-eb70-42e2-85c4-24d2ac156c5en%40googlegroups.com.

[Brian]

Hi Finn,

Thank you very much for confirming. Glad to hear that the use case may inform the "grouped" fields work; I will follow the github discussions with interest.