1
$\begingroup$

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$.

$\endgroup$

1 Answer 1

2
$\begingroup$

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)

enter image description here

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()

enter image description here

$\endgroup$
9
  • $\begingroup$ 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. $\endgroup$
    – Xu Wang
    May 24, 2020 at 14:37
  • $\begingroup$ 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. $\endgroup$
    – usεr11852
    May 24, 2020 at 15:13
  • $\begingroup$ If it is the smoothness constraint that breaks the deal, I will write something with cgam. $\endgroup$
    – usεr11852
    May 24, 2020 at 15:22
  • $\begingroup$ 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. $\endgroup$
    – Xu Wang
    May 24, 2020 at 17:18
  • $\begingroup$ 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. $\endgroup$
    – usεr11852
    May 24, 2020 at 18:36

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.