[R] mgcv: pcls() makes everything linear

[Simon Wood]

I can't see anything immediately wrong except:

1. presumably there are repeated values in 'road_quiet' aren't there? If 
so then your inequality constraint matrix
will contain constraints that are *exact* copies of each other. I'm not 
sure, and don't immediately have time to try it
out, but it could be that this will cause a problem (what might happen 
is that you hit a constraint, project into its null space, but 
numerically the search direction then appears to be about to violate a 
copy of the constraint, so the copy gets added to the set of active 
constraints, and then another copy, and so on until it appears that the 
problem is fully constrained). Could you try replacing your constraints 
with some based on a grid  of e.g. 50 points spread evenly over the 
range of 'road_quiet' values. (If that does solve the problem please let 
me know, and I can add a check for constraint uniqueness to pcls.)

2. A long shot, but what is the range of values of road_quiet? eps 
should probably be around 1e-7 times the typical size of road_quiet (by 
'around' I mean something in the range 1e-4 to 1e-9, say). If eps is 
much too small or much too large
then the constraints become poor approximations to the gradient of the 
spline. I doubt this is the problem, however.

If those two fail you could send me some data (offline) with which I can 
reproduce the problem (smaller the better!), and I can dig further (on 
usual - only for this purpose - basis).

best,

Simon

On 23/07/13 10:19, Kathrine Veie wrote:
> 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
> ______________________________________________
> 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.