# Linear model with constraints

I'm quite new to R and I have a following problem:

I have a simple 2-factor linear model:

# factor1 has 8 categorical values, factor2 has 6 categories
Rate ~ factor1 + factor2
model1 <- lm(Rate~factor1+factor2, data=myData)

And want to put constraints SUM of factor1 coefficients = 0, the same for factor2.

None of the manuals gives any clue how to do this.

I found a link to similar problem here but it is different and I couldnt figure out how to modify it...

-
Please do not crosspost simultaneously here and on SO. –  Dirk Eddelbuettel Sep 28 '10 at 20:47
Certainly, Sir :) –  Vytautas Oct 3 '10 at 11:04

You can do this using contrasts:

options(contrasts=c('contr.sum', 'contr.sum'))

UPDATE: A little googling reveals a page which might be a little clearer:

-
Thanks! It did work, question: it did not force the condition that coefficients sum up to 0, but minimized the sums for both factors. Is there a way to force the ==0 condition, or it simply does not do that since it would affect the error ? –  Vytautas Sep 28 '10 at 19:26
@Vytautas: I think you're forgetting the reference level for your factor which is not shown (since it's the baseline for the others). –  ars Sep 28 '10 at 21:14
@Vytautas or take a look at Interpreting model matrix columns when using contr.sum, j.mp/9pNFQe. –  chl Sep 28 '10 at 21:43

There is a trick to be had here. For simplicity, suppose you are trying to build a model of the form $$y = b_1 x_1 + b_2 x_2 + b_3 x_3,$$ subject to $b_1 + b_2 + b_3 = 0$. Simply re-express $b_3$ as $b_3 = - b_1 - b_2$, which is to say you are trying to build a model of the form $$y = b_1 x_1 + b_2 x_2 - (b_1 + b_2) x_3 = b_1 (x_1 - x_3) + b_2 (x_2 - x_3).$$ So create new variables $\tilde{x}_1 = x_1 - x_3,$ and $\tilde{x}_2 = x_2 - x_3$, and perform your regression using these transformed variables as the independent variables.

You should be able to apply this trick independently to the catgories of factor1 and factor2. (I have assumed the data is given as 0/1 indicators of membership in the individual categories.)

-
Huh, factors are well handled using the factor command, so no need to rely on dummy coding (otherwise, it's accessible through model.matrix(lm(x~a*b)) for my toy example). –  chl Sep 28 '10 at 18:01
@chl thanks. I know very little R, but suspected this simple algebraic trick could be applied in any case, or at least illustrate what is going on 'under the hood' –  shabbychef Sep 28 '10 at 18:05
Yes, definitely! this was just a precision. (and you have my +1) –  chl Sep 28 '10 at 18:14

The easiest is to use the appropriate built-in function:

myContrasts <- list(factor1=contr.sum(length(levels(factor1))),
factor2=contr.sum(length(levels(factor2))))

model1 <- lm(Rate ~ factor1 + factor2, data=myData, contrasts=myContrasts)
-

There are packages for it.

One example is glmc - Generalized Linear Models Subject to Constraints.

See documentation here: http://cran.r-project.org/web/packages/glmc/glmc.pdf

A very good overview on optimization techniques with the use of constaints can be found here: http://zoonek.free.fr/blosxom/R/2012-06-01_Optimization.html

Example code (from the glmc documentaion):

library(glmc)
#Specify the data.
n <- rbind(c(5903,230),c(5157,350))
mat <- matrix(0,nrow=sum(n),ncol=2)
mat <- rbind(matrix(1,nrow=n[1,1],ncol=1)%*%c(0,0),
matrix(1,nrow=n[1,2],ncol=1)%*%c(0,1),
matrix(1,nrow=n[2,1],ncol=1)%*%c(1,0),
matrix(1,nrow=n[2,2],ncol=1)%*%c(1,1))
#Specifying the population constraints.
gfr <- .06179*matrix(1,nrow=nrow(mat),ncol=1)
g <- matrix(1,nrow=nrow(mat),ncol=1)

amat <- matrix(mat[,2]*g-gfr,ncol=1)

hrh <- data.frame(birth=mat[,2], child=mat[,1], constraints=amat)
gfit <- glmc(birth~child, data=hrh, family="binomial",emplik.method="Owen",
control=glmc.control(maxit.glm=10,maxit.weights=200,