mgcv (version 0.9-6)

pcls: Penalized Constrained Least Squares Fitting

Description

Solves least squares problems with quadratic penalties subject to linear equality and inequality constraints using quadratic programming.

Usage

pcls(M)

Arguments

M

is the single list argument to pcls. It should have the following elements (last 3 are not in argument for mgcv) :

• y

{The response data vector.} w{A vector of weights for the

Value

• The function returns an array containing the estimated parameter vector.

Details

This solves the problem: $$minimise~ \| { \bf W}^{1/2} ({ \bf Xp - y} ) \|^2 + \sum_{i=1}^m \lambda_i {\bf p^\prime S}_i{\bf p}$$ subject to constraints ${\bf Cp}={\bf 0}$ and ${\bf A}_{in}{\bf p}>{\bf b}_{in}$, w.r.t. $\bf p$ given the smoothing parameters $\lambda_i$. ${\bf X}$ is a design matrix, $\bf p$ a parameter vector, $\bf y$ a data vector, $\bf W$ a diagonal weight matrix, ${\bf S}_i$ a positive semi-definite matrix of coefficients defining the ith penalty and $\bf C$ a matrix of coefficients defining the linear equality constraints on the problem. The smoothing parameters are the $\lambda_i$. Note that ${\bf X}$ must be of full column rank, at least when projected into the null space of any equality constraints. ${\bf A}_{in}$ is a matrix of coefficients defining the inequality constraints, while ${\bf b}_{in}$ is a vector involved in defining the inequality constraints.

Quadratic programming is used to perform the solution. The method used is designed for maximum stability with least squares problems: i.e. ${\bf X}^\prime {\bf X}$ is not formed explicitly. See Gill et al. 1981.

References

Gill, P.E., Murray, W. and Wright, M.H. (1981) Practical Optimization. Academic Press, London.

Wood, S.N. (1994) Monotonic smoothing splines fitted by cross validation SIAM Journal on Scientific Computing 15(5):1126-1133

http://www.stats.gla.ac.uk/~simon/

See Also

mgcv mono.con

Examples

# first an un-penalized example - fit E(y)=a+bx subject to a>0
set.seed(0)
n<-100
x<-runif(n);y<-x-0.2+rnorm(n)*0.1
M<-list(X=matrix(0,n,2),p=c(0.1,0.5),off=array(0,0),df=array(0,0),S=0,
Ain=matrix(0,1,2),bin=0,C=matrix(0,0,0),sp<-0,y=y,w=y*0+1)
M$X[,1]<-1;M$X[,2]<-x;M$Ain[1,]<-c(1,0)
pcls(M)->M$p
plot(x,y);abline(M$p,col=2);abline(coef(lm(y~x)),col=3)

# and now a penalized example: a monotonic penalized regression spline .....

# Generate data from a monotonic truth.
x<-runif(100)*4-1;x<-sort(x);
f<-exp(4*x)/(1+exp(4*x));y<-f+rnorm(100)*0.1;plot(x,y)
dat<-data.frame(x=x,y=y)
# Show regular spline fit (and save fitted object)
f.ug<-gam(y~s(x,k=10,bs="cr"));lines(x,fitted(f.ug))
# Create Design matrix, constriants etc. for monotonic spline....
gam.setup(y~s(x,k=10,bs="cr")-1,dat,fit.method="mgcv")->G;
GAMsetup(G)->G;F<-mono.con(G$xp);
G$Ain<-F$A;G$bin<-F$b;G$C<-matrix(0,0,0);G$sp<-f.ug$sp;
G$p<-G$xp;G$y<-y;G$w<-y*0+1;

p<-pcls(G);  # fit spline (using s.p. from unconstrained fit)

# now modify the gam object from unconstrained fit a little, to use it
# for predicting and plotting constrained fit. 
p<-c(0,p);f.ug$coefficients<-p; 
lines(x,predict.gam(f.ug,newdata=data.frame(x=x)),col=2)

# now a tprs example of the same thing....

f.ug<-gam(y~s(x,k=10));lines(x,fitted(f.ug))
# Create Design matrix, constriants etc. for monotonic spline....
gam.setup(y~s(x,k=10),dat,fit.method="mgcv")->G;
GAMsetup(G)->G;
nc<-40         # number of constraints
xc<-0:nc/nc # points on [0,1]  
xc<-xc*4-1  # points at which to impose constraints
A0<-predict.gam(f.ug,data.frame(x=xc),type="lpmatrix") 
# ... A0%*%p evaluates spline at xc points
A1<-predict.gam(f.ug,data.frame(x=xc+1e-6),type="lpmatrix") 
A<-(A1-A0)/1e-6    
# ... approx. constraint matrix (A%*%p is -ve spline gradient at points xc)
G$Ain<-A;    # constraint matrix
G$bin<-rep(0,nc);  # constraint vector
G$sp<-f.ug$sp; # use smoothing parameters from un-constrained fit
G$p<-rep(0,11);G$p[11]<-0.1  
# ... monotonic start params, got by setting coefs of polynomial part
G$p[10]<- -mean(0.1*x) 
# ... must meet gam side conditions: sum of smooth over x's is zero 
G$y<-y;G$w<-y*0+1
p<-pcls(G);  # fit spline (using s.p. from unconstrained fit)

# now modify the gam object from unconstrained fit a little, to use it
# for predicting and plotting constrained fit. 
f.ug$coefficients<-p; 
lines(x,predict.gam(f.ug,newdata=data.frame(x=x)),col=3)