rconstraintsgammgcv

How to correctly specify lower and upper bounds in mgcv mono.con function

In the documentation for mgcv::pcls(), an example of constructing a monotonically increasing gam model is provided with the following code:

set.seed(1234) # added for reproducibility in this post only
x <- runif(100)*4-1;x <- sort(x);
f <- exp(4*x)/(1+exp(4*x)); y <- f+rnorm(100)*0.1; plot(x,y)
dat <- data.frame(x=x,y=y)
# Show regular spline fit (and save fitted object)
f.ug <- gam(y~s(x,k=10,bs="cr")); lines(x,fitted(f.ug))
# Create Design matrix, constraints etc. for monotonic spline....
sm <- smoothCon(s(x,k=10,bs="cr"),dat,knots=NULL)[[1]]
F <- mono.con(sm$xp);   # get constraints
G <- list(X=sm$X,C=matrix(0,0,0),sp=f.ug$sp,p=sm$xp,y=y,w=y*0+1)
G$Ain <- F$A;G$bin <- F$b;G$S <- sm$S;G$off <- 0

p <- pcls(G);  # fit spline (using s.p. from unconstrained fit)

fv<-Predict.matrix(sm,data.frame(x=x))%*%p
lines(x,fv,col=2)

This produces the follow plot (if set.seed(1234) is used):

enter image description here

Question:

How does one correctly constrain the monotonically increasing model so that the lower bound is 0.1 and the upper bound is 1.0? I naively initially thought one could set these values as the lower and upper bound values in mono.con. For example F<-mono.com(sm$xp, lower=0.1, upper=1.0), but this will return the following value for F$b:

 [1]  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
[25]  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.1 -1.0

(i.e. lower and -1*upper have been appended to F$b (and dim(F$A)[1] has similarly increased by 2), but the call to pcls(G) returns the following: Error in pcls(G) : initial parameters not feasible

Is there a modification to example code in ?pcls that can constrain the monototonically increasing gam to between the blue lines in this image below:

enter image description here

Solution

  • Here is the relevant code:

    F<-mono.con(sm$xp, lower=0.1,upper=1)
F$A[4*(10-1)+1,1] <- 1
G <- list(X=sm$X,C=matrix(0,0,0),sp=f.ug$sp,p=scales::rescale(sm$xp,to=c(0.4,0.7)),y=y,w=y*0+1)
G$Ain <- F$A;G$bin <- F$b;G$S <- sm$S;G$off <- 0

p <- pcls(G)

> p
 [1] 0.1000000 0.1346833 0.2752035 0.5562702 0.7425019 0.9762224 0.9891460 1.0000000 1.0000000 1.0000000

    Few notes:

    1. The error message mentions initial parameters not feasible. Therefore, we can rescale the initial parameter (G$p) to satisfy the constraint.

    2. There is a bug when constructing inequality matrix in mono.con() as the lower bound constraint becomes [0 0 0 0 0...]x >= 0.1 instead of [1 0 0 0 0]x >= 0.1. I think the issues is this where 4*n should be 4*n+lo. This is why I temporarily add F$A[4*(10-1)+1,1] <- 1 .

    For decreasing constraint, we would do the following:

    F<-mono.con(sm$xp,up=FALSE, lower=0.1,upper=1)
F$A[4*(10-1)+1,10] <- 1
F<-mono.con(sm$xp,up=FALSE)
G <- list(X=sm$X,C=matrix(0,0,0),sp=f.ug$sp,p=scales::rescale(sm$xp,to=c(0.7,0.4)),y=y,w=y*0+1)
G$Ain <- F$A;G$bin <- F$b;G$S <- sm$S;G$off <- 0

p <- pcls(G);  # fit spline (using s.p. from unconstrained fit)
p
[1] 1.0000000 1.0000000 1.0000000 0.9891460 0.9762224 0.7425019 0.5562702 0.2752035 0.1346833 0.1000000

    Namely, sm$xp is rescaled to c(0.7,0.4) and the lower bound constraint becomes [0 0 ... 0 0 1]x>=0.1 instead of [0 0 ... 0 0 0]x>=1.