| mono.con {mgcv} | R Documentation |
Finds linear constraints sufficient for monotonicity (and
optionally upper and/or lower boundedness) of a cubic regression
spline. The basis representation assumed is that given by the
gam, "cr" basis: that is the spline has a set of knots,
which have fixed x values, but the y values of which constitute the
parameters of the spline.
mono.con(x,up=TRUE,lower=NA,upper=NA)
x |
The array of knot locations. |
up |
If TRUE then the constraints imply increase, if
FALSE then decrease. |
lower |
This specifies the lower bound on the spline unless it is
NA in which case no lower bound is imposed. |
upper |
This specifies the upper bound on the spline unless it is
NA in which case no upper bound is imposed. |
Consider the natural cubic spline passing through the points:
(x_i,p_i), i=1..n. Then it is possible
to find a relatively small set of linear constraints on p
sufficient to ensure monotonicity (and bounds if required):
Ap>=b. Details are given in Wood (1994).
This function returns a list containing A and b.
The function returns a list containing constraint matrix
A and constraint vector b.
Simon N. Wood simon@stats.gla.ac.uk
Gill, P.E., Murray, W. and Wright, M.H. (1981) Practical Optimization. Academic Press, London.
Wood, S.N. (1994) Monotonic smoothing splines fitted by cross validation SIAM Journal on Scientific Computing 15(5):1126-1133
http://www.stats.gla.ac.uk/~simon/
# Fit a monotonic penalized regression spline ..... # Generate data from a monotonic truth. set.seed(10);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, constriants etc. for monotonic spline.... gam(y~s(x,k=10,bs="cr")-1,data=dat,fit=FALSE)->G F<-mono.con(G$smooth[[1]]$xp); G$Ain<-F$A;G$bin<-F$b;G$C<-matrix(0,0,0);G$sp<-f.ug$sp; G$p<-G$smooth[[1]]$xp;G$y<-y;G$w<-y*0+1 p<-pcls(G); # fit spline (using s.p. from unconstrained fit) # now modify the gam object from unconstrained fit a little, to use it # for predicting and plotting constrained fit. p<-c(0,p);f.ug$coefficients<-p; x<-seq(min(x),max(x),length=200) lines(x,predict.gam(f.ug,newdata=data.frame(x=x)),col=2)