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.