There's an example at the end of the pcls help file in version 1.5-0 --- just submitted to CRAN.
best, Simon On Wednesday 25 February 2009 18:48, Benjamin Tyner wrote: > Hi, > > Does anyone know how to fit a GAM where one or more smooth terms are > constrained to be monotonic, in the presence of "by" variables or > other terms? I looked at the example in ?pcls but so far have not been > able to adapt it to the case where there is more than one predictor. > For example, > > require(mgcv) > > set.seed(0) > n<-100 > > > # Generate data from a monotonic truth. > x<-runif(100)*4-1 > x<-sort(x) > > m <- structure(rep(1:2,50), .Label=c("one","two"), class="factor") > > f<-as.integer(m) * exp(4*x)/(1+exp(4*x)) > y<-f+rnorm(100)*0.1 > plot(x,y,xlim=c(min(x), max(x)*2)) > dat<-data.frame(x=x,m=m,y=y) > > > > # Show regular spline fit (and save fitted object) > f.ug<-gam(y~m+s(x,k=10,by=m,bs="cr")) > > bool <- m=="one" > yhat <- fitted(f.ug) > > lines(x[bool],yhat[bool]) > lines(x[!bool],yhat[!bool]) > > xx <- seq(max(x), 2*max(x), length=100) > mm <- structure(rep(1:2,50), .Label=c("one","two"), class="factor") > yy <- predict(f.ug, newdata=data.frame(m=mm,x=xx)) > > bool <- mm=="one" > > lines(xx[bool],yy[bool], lty=2) # show effect of extrapolation > lines(xx[!bool],yy[!bool], lty=2) > > > # this is where I start running into trouble > sm<-smoothCon(s(x,k=10,by=m,bs="cr"),dat,knots=NULL)[[1]] > FF<-mono.con(sm$xp); # get constraints > G<-list(y=y, > w=rep(1, n), > X=sm$X, > C=matrix(0,0,0), > S = sm$S, > off = 0, > sp=f.ug$sp, > p=sm$xp, > Ain = FF$A, > bin = FF$b > ) > > > p<-pcls(G) # fit spline (using s.p. from unconstrained fit) > fv<-Predict.matrix(sm,data.frame(x=x))%*%p > > # can we do this without calling smoothCon directly ? > # also having trouble here. > f.nofit<-gam(y~m+s(x,k=10,by=m,bs="cr"),fit=FALSE) > FF2 <- mono.con(f.nofit$smooth[[1]]$xp) > > stopifnot(identical(FF, FF2)) > > G2 <- list(y = f.nofit$y, > w = f.nofit$w, > X = f.nofit$X, > C = f.nofit$C, > S = f.nofit$smooth[[1]]$S, > off = f.nofit$off, > sp = f.ug$sp, > p = f.nofit$smooth[[1]]$xp, > Ain = FF2$A, > bin = FF2$b > ) > > p2 <- pcls(G2) > fv2<-Predict.matrix(f.nofit$smooth[[1]],data.frame(x=x))%*%p2 > > > > Many thanks > -Ben > > ______________________________________________ > R-help@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 Sciences, University of Bath, Bath, BA2 7AY UK > +44 1225 386603 www.maths.bath.ac.uk/~sw283 ______________________________________________ R-help@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.