[R] monotonic GAM with more than one term

Simon Wood s.wood at bath.ac.uk
Tue Mar 3 11:53:37 CET 2009


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




More information about the R-help mailing list