[R] mgcv: pcls() makes everything linear

[Kathrine Veie]

Dear R helplisters,

I am trying to implement a mononicity constraint on a smooth term in my generalized additive model with the mgcv package (v. 1.7-18). I adapted the code from an example given in the help file for pcls(). The example runs just as one would expect, but when I adapt it and use the code on my data, the code results in a linear fit on ALL smooth terms in the model even though I only place restrictions on a single one of them. 

Can anyone give me any idea why this is the case?  The basis spline dimension is taken from the unconstrained object, but the coefficients determined by pcls() are such that the fit is linear (basically very, very small coefficients and large coefficients on a single spline component for each covariate as far as I can tell). (see figure which I hope comes through to the help list posting!)

The example is simpler than the model I actually want to estimate, but the problem is the same in the more expanded version. In the simple version, the curve I want to restrict to be monotonically increasing already satisfies the constraint, but in the full model it is more "wiggly" and has some downward sloping parts. 

In the code below I impose the constraint using a subset of the data, but I tried the same thing using the full data set for predictions and except for that being a bit slower the results are all the same. (The full data set has 14,000 observations)

## Preliminary unconstrained gam fit...
G <- gam(ksumphk~s(arlsaml)+s(road_quiet)+s(centrum), data=per1.200m,family=gaussian(link=log), fit=FALSE)
b <- gam(G=G)
c <- b ##Save unconstrained estimates in separate model 

## generate constraints, by finite differencing
## using predict.gam ....
eps <- 1e-7

#Sample from the data set to avoid too many obs in prediction

pd0 <- data.frame(per1.200m[sample(nrow(per1.200m),200),])
pd1 <- pd0
pd1$road_quiet<- pd0$road_quiet +eps

X0 <- predict(b,newdata=pd0,type="lpmatrix")
X1 <- predict(b,newdata=pd1,type="lpmatrix")
Xz <- (X1-X0)/eps


G$Ain <- rbind(Xz) ## inequality constraint matrix
G$bin <- rep(0,nrow(G$Ain))
G$sp <- b$sp
G$C<-matrix(0,0,0)
G$p <- coef(b)
G$off <- G$off-1 ## to match what pcls is expecting
## force inital parameters to meet constraint
G$p[11:18] <- 0.0
p <- pcls(G) ## constrained fit
par(mfrow=c(2,3))
plot(c) ## original fit
b$coefficients <- p
plot(b) ## constrained fit
## note that standard errors


Any help or suggestions on where to read more about pcls() would be very much appreciated!

Kind regards,
Kathrine