[R] Use pcls in "mgcv" package to achieve constrained cubic spline

[Simon Wood]

Actually, you might as well use "gam" directly for this (has the 
advantage that the smoothing parameter will be chosen correctly subject 
to the constraint). Here is some code. Key idea is to set the basis and 
penalty for the spline up first, apply the constraint, and then use gam 
to fit it...

best,
Simon

## Example constraining spline to pass through a
## particular point (0,.6)...

## Fake some data...

library(mgcv)
set.seed(0)
n <- 100
x <- runif(n)*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)

## Create a spline basis and penalty, making sure there is a knot
## at the constraint point, (0 here, but could be anywhere)
knots <- data.frame(x=seq(-1,3,length=9)) ## create knots
## set up smoother...
sm <- smoothCon(s(x,k=9,bs="cr"),dat,knots=knots)[[1]]

## 3rd parameter is value of spline at knot location 0,
## set it to 0 by dropping...
X <- sm$X[,-3]        ## spline basis
S <- sm$S[[1]][-3,-3] ## spline penalty
off <- y*0 + .6       ## offset term to force curve through (0, .6)

## fit spline constrained through (0, .6)...
b <- gam(y ~ X - 1 + offset(off),paraPen=list(X=list(S)))
lines(x,predict(b))

## compare to unconstrained fit...

b.u <- gam(y ~ s(x,k=9),data=dat,knots=knots)
lines(x,predict(b.u))


On 11/03/13 18:06, Victor hyk wrote:
> Hello everyone,
>           Dr. wood told me that I can adapting his example to force cubic spline to pass through certain point.
>           I still have no idea how to achieve this. Suppose we want to force the cubic spline to pass (1,1), how can
> I achieve this by adapting the following code?
>
> # Penalized example: monotonic penalized regression spline .....
> # 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)
> 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)
> lines(x,f,col="blue")
>             Thanks a lot!!
>
>              Victor
> 	[[alternative HTML version deleted]]
>
>
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>


-- 
Simon Wood, Mathematical Science, University of Bath BA2 7AY UK
+44 (0)1225 386603               http://people.bath.ac.uk/sw283