R help - Mar 2013 - linear model with equality and inequality (redundant) constraints

Dear R-users,

in the last days I have been trying to estimate a normal linear model with
equality and inequality constraints.

Please find below a simple example of my problem. 

Of course, one could easily see that, though the constraints are consistent,
there is some redundancy in the specific constraints. Nevertheless my actual
applications can get much larger and I would not like to manually check for
presence of redundant constraints in the setting.

First question: would it be possible to estimate a linear model with equality
and inequality redundant constraints?

I also understand that quadratic programming is the key for that, but I was not
able to ?translate? my simple linear problem with constraints for using
solve.QP(quadprog) and ic.est(ic.infer). I have found pcls(mgcv) more
user-friendly. And below you could find my attempt to fit the toy-data with this
function with and without redundant constraints.

Thanks in advance for any help you could provide,
Carlo Giovanni Camarda


rm(list = ls())
library(mgcv)

## estimate p: E(r) = X p, s.t.
##          C p == c
##          A p >= b (though pcls allows only for A p > b)

## response
r <- c(4277.1265,57004.1177,7135.9204,541.6218,1455.2135)
## design matrix
x <- c(6029,0,0,0,0,0,6029,0,
       0,0,0,44453,0,0,0,0,
       1284,0,0,0,0,3042,0,
       0,0,0,10132,0,0,0,0,0,
       10132,0,0,0,0,1955,0,
       0,0,0,3519,0,0,0,0,0,
       10132,0,0,0,0,0,44453,
       0,0,0,0,1955)
X <- matrix(x, 5, 12)
## equality constraints matrix
C <- matrix(0,7,12)
C[1,1:2] <- 1
C[2,c(3,11)] <- 1
C[3,4] <- 1
C[4,5] <- 1
C[5,c(6,7,10)] <- 1
C[6,c(8,12)] <- 1
C[7,9] <- 1
## equality constraints vector, which is not used in pcls
c <- rep(1,7)
## inequality constraints matrix: p \in [0,1] 
A <- rbind(diag(12),
           -1*diag(12))
## inequality constraints vector
b <- c(rep(0,12), rep(-1,12))
## starting values
p.st <- c(0.5,0.5,0.5,1,1,0.2,0.3,0.5,1,0.5,0.5,0.5)
## test starting values
table(C%*%p.st == c)
table(A%*%p.st > b) ## actually I would like to have A p >= b
## list for pcls
M <- list(X=X,
          p=p.st,
          off=array(0,0),
          S=list(),
          Ain=A,
          bin=b,
          C=C,
          sp=array(0,0),
          y=r,
          w=r*0+1)
## estimation????
(p.hat <- pcls(M))

######## avoid redundant constraints
## response
r1 <- c(4277.1265,46649.1177,3616.9204,-9590.3782,-44952.7865)
## design matrix
X1 <- c(6029,-6029,0,0,0,0,44453,
        0,0,-44453,0,10132,0,-10132,
        0,0,0,10132,-10132,0,0,0,
        1955,0,-1955)
dim(X1) <- c(5,5)
## inequality constraints matrix
A1 <- rbind(diag(5),
            -1*diag(5),
            c(0,0,-1,-1,0))
## inequality constraints vector
b1 <- c(rep(0,5), rep(-1,5), -1)
p.st1 <- runif(5, 0.1, 0.2)
## check constraints
A1%*%p.st1 > b1
M1 <- list(X=X1,
           p=p.st1,
           off=array(0,0),
           S=list(),
           Ain=A1,
           bin=b1,
           C=matrix(0,0,0),
           sp=array(0,0),
           y=r1,
           w=r1*0+1)
(p.hat1 <- pcls(M1))

A1%*%p.hat1 > b1

----------
This mail has been sent through the MPI for Demographic Research.  Should you
receive a mail that is apparently from a MPI user without this text displayed,
then the address has most likely been faked. If you are uncertain about the
validity of this message, please check the mail header or ask your system
administrator for assistance.