# Bivariate isotonic regression on unbalanced grid

I would like to do classic isotonic regression on bivariate data (that is to say, x with two columns). The function biviso of package Iso is very fast but requires a balanced grid. I would like to do isotonic regression on unbalanced grid. For example:

x1 <- runif(100)
x2 <- runif(100)
y <- x1*x2 + rnorm(100)


Is there an R package or custom R code that can do this? By "classic" isotonic regression, I attempt to make reference to the following definition, which is in the Wikipedia article:

$${\displaystyle \min \sum _{i=1}^{n}w_{i}(x_{i}-a_{i})^{2}}$$

$${{\text{subject to }}x_{i}\leq x_{j}{\text{ for all }}(i,j).}$$

The weights $$w_i$$ are not important for my usage so we can define $$w_i=1$$.

The easiest way I think of is using monotonic splines through scam. That way we will constraint the fit of each term to be a monotonic function giving us a monotonic response. I prefer it to other functions that solve the constrained optimisation problems directly because it gives us access to (most of) the diagnostic goodies we get from mgcv::gam.

set.seed(5)
x1 <- runif(100)
x2 <- runif(100)
y <- x1*x2 + rnorm(100)
my_data = data.frame(x1=x1, x2=x2, y=y)
library(scam)
my_model_sc <- scam(y ~ s(x1, bs="mpi")+ s(x2, bs="mpi"), data=my_data)


A quick visual inspection confirm that our fit is monotonic along x1 and x2 too.

library(pdp)
pd <- partial(my_model_sc, pred.var = c("x1", "x2"), smooth = FALSE,
grid.resolution = 100, progress = "text" )
plotPartial(pd)


If we are looking to replicate something closer to the step changes seen in QP solutions, we can look into cgam.

library(cgam)
my_model_cg <- cgam(y ~ incr(x1)+ incr(x2), data=my_data)


Here we can specify an increasing but smoothed fit through incr (they are other options to have smooth convex or concave fits, etc).

I append a bit of code below showing how both functions allow us to have a monotonic fit across values of x2 for a fixed x1.

my_data_2=my_data[order(my_data$$x2),] my_data_2$$x1=0.5;

par(mfrow=c(1,2))
plot(x=my_data_2$$x2, y = predict(my_model_cg, newData = my_data_2)$$fit, lwd=2,
t='l', xlab="x2", ylab="f(x1=0.5, x2)", ylim=c(-0.4,0.5), main="CGAM fit")
grid()
plot(x=my_data_2\$x2, y = predict(my_model_sc,newdata =my_data_2), lwd=2,
t='l', xlab="x2", ylab="f(x1=0.5, x2)", ylim=c(-0.4,0.5), main="SCAM fit")
grid()


• Thank you for your kind time! I am hoping to use classical isotonic regression as defined on Wikipedia (en.wikipedia.org/wiki/Isotonic_regression). I should have added link as reference. I apologize for my carelessness. Thank you for this new technique. I also am interested in studying it one day. May 24, 2020 at 14:37
• Nothing to apologize for. That said, this is: "fitting a free-form line to a sequence of observations such that the fitted line is non-decreasing (or non-increasing) everywhere, and lies as close to the observations as possible." (as Wikipedia defines it) We use monotonic regression like Dette et al. (2006) in monreg were we circumvent the quadratic programming formulation. Just instead of smoothing via Nadaraya-Watson (as in that paper) we smooth with splines. May 24, 2020 at 15:13
• If it is the smoothness constraint that breaks the deal, I will write something with cgam. May 24, 2020 at 15:22
• You are right. I did not know "isotonic regression" was not well-defined (?). I have updated the question to refer to "classic" isotonic regression and given definition. I would be interested in how you remove smoothness constraint with cgam. It would be interesting to compare performance with the classical isotonic regression. However, I am still interested in classical isotonic regression. If there is an answer that allows classical isotonic regression I will choose it. Thank you for your time and flexibility. +1. May 24, 2020 at 17:18
• Cool, I am glad I could help. I added a relevant example. Please note that both examples provided, do minimise exactly what is described in the Wikipedia page subject to smoothness constraints (in the case of scam) on the isotonic estimator. May 24, 2020 at 18:36