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@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.
______________________________________________
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.
