[R] linear model with equality and inequality (redundant) constraints

[Camarda, Carlo Giovanni]

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.