mgcv: pcls() makes everything linear
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.