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 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 Sciences, University of Bath, Bath, BA2 7AY UK
> +44 1225 386603 www.maths.bath.ac.uk/~sw283