Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Force certain parameters to have positive coefficients in lm()

I would like to know how to constrain certain parameters in lm() to have positive coefficients. There are a few packages or functions (e.g. display) that can make all coefficients, and the intercept, positive.

For instance, in this example, I would like to force only x1 and x2 to have positive coefficients.

    x1=c(NA,rnorm(99)*10)
    x2=c(NA,NA,rnorm(98)*10)
    x3=rnorm(100)*10
    y=sin(x1)+cos(x2)-x3+rnorm(100)

    lm(y~x1+x2+x3)

    Call:
      lm(formula = y ~ x1 + x2 + x3)       
    Coefficients:
      (Intercept)           x1           x2           x3  
    -0.06278      0.02261     -0.02233     -0.99626

I have tried function nnnpls() in package nnls, it can control the coefficient sign easily. Unfortunately I can't use it due to issues with NAs in the data as this function doesn't allow NA.

I saw function glmc() can be used to apply constraints but I couldn't get it working.

Could someone let me know what should I do?

like image 725
Yu Deng Avatar asked Dec 02 '14 08:12

Yu Deng


1 Answers

You can use package penalized:

set.seed(1)

x1=c(NA,rnorm(99)*10)
x2=c(NA,NA,rnorm(98)*10)
x3=rnorm(100)*10
y=sin(x1)+cos(x2)-x3+rnorm(100)
DF <- data.frame(x1,x2,x3,y)

lm(y~x1+x2+x3, data=DF)
#Call:
#lm(formula = y ~ x1 + x2 + x3, data = DF)
#
#Coefficients:
#(Intercept)           x1           x2           x3  
#   -0.02438     -0.01735     -0.02030     -0.98203  

This gives the same:

library(penalized)

mod1 <- penalized(y, ~ x1 + x2 + x3, ~1, 
                  lambda1=0, lambda2=0, positive = FALSE, data=na.omit(DF))
coef(mod1)
#(Intercept)          x1          x2          x3 
#-0.02438357 -0.01734856 -0.02030120 -0.98202831 

If you constraint the coefficients of x1 and x2 to be positive, they become zero (as expected):

mod2 <- penalized(y, ~ x1 + x2 + x3, ~1, 
                  lambda1=0, lambda2=0, positive = c(T, T, F), data=na.omit(DF))
coef(mod2)
#(Intercept)          x3 
#-0.03922266 -0.98011223 
like image 187
Roland Avatar answered Oct 20 '22 13:10

Roland