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.

[Brian]

Hi Finn,

Following up on the above thread where you confirmed the increments-via-bru_mapper_marginal construction for a monotone non-positive function of cumulative treatment dose. That approach worked well. But because f(dose) lives on the logit scale, the interpretability is a challenge since the same increment represents a different absolute reduction depending on the baseline. For the application we need the output the prevalence scale directly. I have therefore moved to a prevalence-scale multiplicative structure where the observed prevalence is the counterfactual prevalence times a retained fraction:

P_i = expit(eta*_i) * M(d_i)

logit(P*_i) = eta*_i = intercept + covariates + SPDE field

with d_i the coverage-weighted effective cumulative dose. For M(d) I am using a Weibull retained fraction: M(d) = exp(-(kappa * d)^alpha),   kappa > 0,  alpha > 0 with kappa = exp(lambda) and alpha = exp(nu), both implemented as unconstrained Gaussian scalars. This gives M(0) = 1 and M strictly decreasing. The full predictor on the logit (link) scale is:

logit(P_i) = logit( expit(eta*_i) * exp(-(exp(lambda)*d_i)^exp(nu)) ) = qlogis( plogis(eta*_i) * exp(-(exp(lambda)*d_i)^exp(nu)) )

The implementation:

```

n <- 800
M <- 6
kappa_true <- 0.55
alpha_true <- 0.65
beta0_true <- -0.4
beta1_true <- 0.8

# We simulate dose with high doses sparse
dose_probs <- c(0.25, 0.20, 0.18, 0.15, 0.10, 0.07, 0.05)
dose <- sample(0:M, n, replace = TRUE, prob = dose_probs)

x1 <- rnorm(n)
dt <- data.table(x1 = x1, dose = dose)
for (j in 1:M) dt[, paste0("c_", j) := runif(n, 0.60, 0.90)] # These are just the coverages of the treatment rounds

dt[, d_eff := Reduce(`+`, map(1:M, \(j) {
  as.integer(dose >= j) * dt[[paste0("c_", j)]]
}))]

p_star <- plogis(beta0_true + beta1_true * x1)
M_true <- exp(-(kappa_true * dt$d_eff)^alpha_true)
p_obs <- p_star * M_true
dt[, N := sample(40:60, n, replace = TRUE)]
dt[, y := rbinom(n, N, p_obs)]

# We declare two Gaussian scalars: lambda (log kappa) and nu (log alpha)
cmp <- ~
  x1(x1, model = "linear") +
  lambda(1, model = "linear", mean.linear = -0.5, prec.linear = 1) +
  nu(1, model = "linear", mean.linear = 0,   prec.linear = 1) +
  Intercept(1)

# The nonlinear predictor on the logit scale
form <- y ~ qlogis(plogis(Intercept + x1) * exp(-(exp(lambda) * d_eff)^exp(nu)))

fit  <- bru(
  cmp,
  formula = form,
  family = "betabinomial",

  data = dt,
  Ntrials = dt$N
)

# Getting the estimated kappa and alpha
kappa_marg <- inla.tmarginal(exp, fit$marginals.fixed[["lambda"]])
alpha_marg <- inla.tmarginal(exp, fit$marginals.fixed[["nu"]])

kappa_s <- inla.zmarginal(kappa_marg)
alpha_s <- inla.zmarginal(alpha_marg)

data.frame(
  parameter = c("kappa", "alpha"),
  truth = c(kappa_true, alpha_true),
  post_mean = c(kappa_s$mean, alpha_s$mean),
  lo95 = c(kappa_s$quant0.025, alpha_s$quant0.025),
  hi95 = c(kappa_s$quant0.975, alpha_s$quant0.975)
)

  parameter truth post_mean      lo95      hi95
1     kappa  0.55 0.5539854 0.5246743 0.5844259
2     alpha  0.65 0.6869553 0.6358073 0.7409575

```

Is expressing the multiplicative prevalence predictor as qlogis(plogis(eta*) * exp(-(exp(lambda)*d)^exp(nu))) the correct approach in inlabru? That is, returning the result on the logit (link) scale so the Beta-Binomial logit link inverts it, or is there a cleaner route such as a custom link or a mapper-based construction that would keep more of the predictor linear and reduce the linearisation cost and generally the implementation.

[Finn Lindgren]

Hi,

Yes, the
 ~ qlogis(plogis(eta) * M(d))

approach is the most straightforward approach.

The next inlabru release (2.15.0, planned for July) will be much faster in the linearisation step for nonlinear models like this (potentially the user needs to let bru() know what assumptions it can make, like the current "is_rowwise" argument to bru_obs), so long-term the speed shouldn't be an issue.  Even the current (2.14.1) version may be faster if you explicitly set is_rowwise=TRUE in the bru_obs() call, as otherwise the automatic structure detection might make the safe but slower assumption that is_rowwise=FALSE.

Finn

To view this discussion, visit https://groups.google.com/d/msgid/r-inla-discussion-group/ea6e8284-4338-4d14-b8c4-c7e8697a5236n%40googlegroups.com.

--

Finn Lindgren
email: finn.l...@gmail.com

[Brian]

Thank you,  for the confirmation and explanation. I have added is_rowwise = TRUE to the bru_obs() call. From the documentation I understand this commits to the predictor at row i depending only on row i of the data which holds for my setup since the predictor qlogis(plogis(Intercept + x1_i + field_i) * exp(-(exp(lambda) * d_i)^exp(nu))) involves only the SPDE field evaluated at s_i and the scalar covariate d_i for that observation, with no cross-row dependencies. Looking forward for the 2.15.0 release.

Thank you again.