monotone splines using package mgcv

monotone_spline.R

library(mgcv)
#library(modules) #devtools::install_github(klmr/modules)
#mgcv<-import_package("mgcv")

mspline<-function(x,y,k=10,lower=NA,upper=NA){
  #fits a monotonic spline to data
  #small values of k= more smoothing (flatter curves)
  #large values of k= more flexible (wiggly curves)
  #k is related to effective degrees of freedom and number of knots
  #use unconstrained gam to get rough parameter estimates
  #lower, upper optional bounds on the function
  #basically a slight modification of an example in the mgcv::pcls documentation
  dat<-data.frame(x=x,y=y)
  init_gam <- gam(y~s(x,k=k,bs="cr"))
  # Create Design matrix, constraints etc. for monotonic spline....
  sm <- smoothCon(s(x,k=k,bs="cr"),dat,knots=NULL)[[1]]
  mc <- mono.con(sm$xp,lower=lower,upper=upper) # monotonicity constraints
  M <- list(X=sm$X, y=y, #design matrix, outcome
            C=matrix(0,0,0), #equality constraints (none)
            Ain=mc$A, bin=mc$b, #inequality constraints
            sp=init_gam$sp, p=sm$xp, #initial guesses for param estimates
            S=sm$S, #smoothness penalty matrix
            w=y*0+1, off=0 #weights, offset
            )
  #fit spine using penalized constrained least squares
  p<-pcls(M)
  return(list(sm=sm,p=p))
}

predict.mspline<-function(msp,x){
  #using the monotone spline msp, predict values for the vector x
  Predict.matrix(msp$sm,data.frame(x=x))%*%msp$p
}

# Generate data from a monotonic truth.
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)
fv<-mspline(x,y,5)
lines(x,predict.mspline(fv,x),col="red")
fv<-mspline(x,y,10)
lines(x,predict.mspline(fv,x),col="blue")
legend("bottomright",lty=1,paste0("k=",c(5,10)),col=c("red","blue"))