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.