CODEWITHSUNDEEP
rregression
tidyx
tidyx

Reputation: 45

how to solve linear regression with conditions in R

I would like to do a linear regression with conditions, like: Y = ax1 + bx2 +c, a>0 , b <1 and c>=0

Y <- c(167, 136, 195, 174, 144, 135, 89, 81, 114, 111)
x1 <- c(2.9, 3.4, 0.7, 1.1, 3.5, 5.0, 6.7, 4.7, 3.7, 8.8)
X2 <- c(60, 47, 63, 62, 40, 60, 50, 35, 40, 40)
mydata <- data.frame(Y,X1,X2)

Is there an easy way to do this in R? .

Upvotes: 1

Views: 254

Answers (1)

Erwin Kalvelagen
Erwin Kalvelagen

Reputation: 16714

Here is how you can do this with CVXR:

> library(CVXR)
> 
> #
> # data
> #
> Y <- c(167, 136, 195, 174, 144, 135, 89, 81, 114, 111)
> X1 <- c(2.9, 3.4, 0.7, 1.1, 3.5, 5.0, 6.7, 4.7, 3.7, 8.8)
> X2 <- c(60, 47, 63, 62, 40, 60, 50, 35, 40, 40)
> 
> #
> # organize as matrix
> #
> X <- cbind(1,X1,X2)
> X
         X1 X2
 [1,] 1 2.9 60
 [2,] 1 3.4 47
 [3,] 1 0.7 63
 [4,] 1 1.1 62
 [5,] 1 3.5 40
 [6,] 1 5.0 60
 [7,] 1 6.7 50
 [8,] 1 4.7 35
 [9,] 1 3.7 40
[10,] 1 8.8 40
> 
> 
> #
> # standard regression
> #
> lm(Y~X1+X2)

Call:
lm(formula = Y ~ X1 + X2)

Coefficients:
(Intercept)           X1           X2  
     82.616       -7.716        1.675  

> 
> #
> # standard regression as QP
> #
> beta <- Variable(3)
> 
> problem <- Problem(Minimize(sum_squares(Y - X %*% beta)))
> result <- solve(problem, verbose=T)
-----------------------------------------------------------------
           OSQP v0.6.0  -  Operator Splitting QP Solver
              (c) Bartolomeo Stellato,  Goran Banjac
        University of Oxford  -  Stanford University 2019
-----------------------------------------------------------------
problem:  variables n = 13, constraints m = 10
          nnz(P) + nnz(A) = 50
settings: linear system solver = qdldl,
          eps_abs = 1.0e-05, eps_rel = 1.0e-05,
          eps_prim_inf = 1.0e-04, eps_dual_inf = 1.0e-04,
          rho = 1.00e-01 (adaptive),
          sigma = 1.00e-06, alpha = 1.60, max_iter = 10000
          check_termination: on (interval 25),
          scaling: on, scaled_termination: off
          warm start: on, polish: on, time_limit: off

iter  objective    pri res    dua res    rho        time
   1   0.0000e+00   1.95e+02   1.70e+05   1.00e-01   3.47e-04s
  50   2.8959e+03   3.20e-07   2.29e-06   1.00e-01   8.43e-04s
plsh   2.8959e+03   4.21e-14   4.51e-13  ---------   1.36e-03s

status:               solved
solution polish:      successful
number of iterations: 50
optimal objective:    2895.8618
run time:             1.36e-03s
optimal rho estimate: 4.57e-02

> result$getValue(beta)
          [,1]
[1,] 82.616486
[2,] -7.716095
[3,]  1.674722
> 
> #
> # add constraints
> #      beta[1] >= 0       (constant term)
> #      beta[2] >= 0.0001  (instead of > we do >= with a small tolerance)
> #      beta[3] <= 0.9999  (instead of < we do <= with a small tolerance)
> #
> problem <- Problem(Minimize(sum_squares(Y - X %*% beta)),
+                    list(beta[1]>=0,
+                         beta[2]>=0.0001,
+                         beta[3]<=0.9999))
> result <- solve(problem, verbose=T)
-----------------------------------------------------------------
           OSQP v0.6.0  -  Operator Splitting QP Solver
              (c) Bartolomeo Stellato,  Goran Banjac
        University of Oxford  -  Stanford University 2019
-----------------------------------------------------------------
problem:  variables n = 13, constraints m = 13
          nnz(P) + nnz(A) = 53
settings: linear system solver = qdldl,
          eps_abs = 1.0e-05, eps_rel = 1.0e-05,
          eps_prim_inf = 1.0e-04, eps_dual_inf = 1.0e-04,
          rho = 1.00e-01 (adaptive),
          sigma = 1.00e-06, alpha = 1.60, max_iter = 10000
          check_termination: on (interval 25),
          scaling: on, scaled_termination: off
          warm start: on, polish: on, time_limit: off

iter  objective    pri res    dua res    rho        time
   1   0.0000e+00   1.95e+02   1.71e+05   1.00e-01   5.30e-04s
 150   7.8534e+03   9.37e-05   5.05e-04   4.62e+00   1.51e-03s
plsh   7.8533e+03   5.55e-14   2.84e-14  ---------   1.94e-03s

status:               solved
solution polish:      successful
number of iterations: 150
optimal objective:    7853.3365
run time:             1.94e-03s
optimal rho estimate: 2.98e+00

> result$getValue(beta)
         [,1]
[1,] 84.90456
[2,]  0.00010
[3,]  0.99990
>

Upvotes: 2

Related Questions