6
$\begingroup$

I am using the nnls() function from the nnls package in R to do a linear regression for regressors $x_i$ and observations $y$. The function delivers beta coefficients $\beta_i\geq{0}, \forall i$. However, is it possible to apply the constraints only to some regressors so that

$$\beta_k \geq 0 \quad k \in \{1...,10\}, k\neq i \\ \beta_i \in \mathbb{R} \quad i \in \{1...,10\}$$

given that I have 10 regressor variables?

nnls() offers the possibility to enforce some coefficients to be negative and others to be positive. I only want the positive constraint for some of them, the other ones can be either positive or negative.

$\endgroup$
  • 2
    $\begingroup$ Investigate the optim() function in R. This function will allow you to maximize a likelihood function given some constraints. $\endgroup$ – TrynnaDoStat Feb 6 '15 at 12:36
  • $\begingroup$ Added other option via reformulating the objective as a quadratic programming problem - this should be faster than the other options mentioned here... $\endgroup$ – Tom Wenseleers Jun 28 '19 at 10:00
3
$\begingroup$

Regardless of what your computing platform may be, you can trick any linearly constrained linear model solver into doing what you want without having to modify any code at all. Observe that the model

$$\mathbb{E}(Y) = X \beta$$

is equivalent to the model

$$\mathbb{E}(Y + X\gamma) = X(\beta + \gamma)$$

for any fixed coefficients $\gamma$, because the constant vector $X\gamma$ has been added to both sides. Moreover, fitting the second model will produce the same covariance matrix of estimates, etc., as the first because no change has been made in the errors whatsoever. Consequently, after estimating $\beta+\gamma$ as $\widehat{\beta+\gamma}$ using any linear procedure (not just nnls), the estimate of $\beta$ is obtained as

$$\hat\beta = \widehat{\beta+\gamma} -\gamma.$$

Exploit this flexibility to make sure that some of the coefficients of $\widehat{\beta+\gamma}$ will naturally be positive (and therefore not constrained by nnls). This can be done by estimating what those coefficients might be, performing the adjusted procedure, and checking that it has not constrained the corresponding estimates. If any have been constrained, increase the amount of adjustment and repeat until success is achieved.

I propose using an initial ordinary least squares fit to start the procedure. To be conservative, change the estimates by some small multiple $\rho$ of their standard errors. If iteration is needed, keep increasing $\rho$. The following code doubles $\rho$ at each iteration. The entire algorithm is contained within the short repeat block in the middle of the code example below.

As an example, I generated a problem with $200$ observations of seven variables (and included a constant term). The following tableau summarizes the results:

            AIntercept  AX1  AX2   AX3  AX4  AX5 AX6 AX7
True              -3.5 -2.5 -1.5 -0.50 0.50 1.50 2.5 3.5
OLS               -3.4 -2.5 -1.6 -0.52 0.44 1.49 2.4 3.5
NNLS               0.0  0.0  0.0  0.00 0.63 0.87 2.3 4.0
NNLS.0            -3.5 -2.7  0.0  0.00 0.75 1.56 2.3 3.6
Constraints        0.0  0.0  1.0  1.00 1.00 1.00 1.0 1.0
Bound              0.0  0.0  1.0  1.00 0.00 0.00 0.0 0.0

The true values of $\beta$ are listed first. The first half are negative, the second half positive. These are followed by their OLS estimates, their NNLS estimates, and the modified NNLS estimates wherein the first two coefficients (AIntercept and AX1) were not constrained to be positive. The last two lines summarize the constraints, printing "1.0" where positivity constraints were applied. The constraints actually imposed in the solution, using the same indicator method, appear last. In this case the procedure worked well.

#
# Describe a problem.
#
p <- 7                              # Number of variables
n <- 200                            # Number of observations
beta <- 0:p - p/2                   # True coefficients
constraints <- c(0, 0, rep(1, p-1)) # Positivity constraint indicator
#
# Generate data.
#
set.seed(17)
A <- cbind(Intercept=rep(1, n), 
           matrix(rnorm(p*n), n, dimnames=list(c(), paste0("X", 1:p))))
b <- A %*% beta + rnorm(n)
#
# OLS for reference.
#
fit.lm <- lm(b ~ A - 1)
#
# NNLS.
#
library(nnls)
fit.nnls <- nnls(A, b)
#
# NNLS with selective constraints.
#
rho <- 2 
coefficients <- coef(fit.lm)
se <- coef(summary(fit.lm))[, "Std. Error"]
repeat {
  beta.0 <- (coefficients - rho*se) * (1 - constraints)
  b.0 <- b - A %*% beta.0
  fit.nnls.0 <- nnls(A, b.0)
  if (all(constraints[fit.nnls.0$bound] == 1)) break #$
  if (rho > 1000) stop("Unable to find a solution.")
  rho <- rho*2
}
fit.nnls.0$x <- fit.nnls.0$x + beta.0
#
# Compare.
#
bound <- rep(0, p+1); bound[fit.nnls.0$bound] <- 1
print(rbind(True=beta, OLS=coef(fit.lm), NNLS=coef(fit.nnls),
      NNLS.0=coef(fit.nnls.0), 
      Constraints=constraints, Bound = bound), digits=2)
$\endgroup$
1
$\begingroup$

The solution is constrOptim(), which allows for individual constraining of coefficients.

For example:

If I want to regress variable $y$ with the two regressors $x_1, x_2$, I can define the linear regression minimization problem (minimizing squared residuals) myself via

min.RSS <- function(data, par){
  with(data, sum((par[1]*x1 + par[2]*x2 - y)^2))
}

All data required is stored in a data frame, column one and two contains regressors $x_1, x_2$, and column three contains the observations $y$:

dat = data.frame(x1=c(1,2), x2=c(2,3), y=c(5,6))

Then the constrained optimization can be conducted via

myObj        = constrOptim(theta, f, grad, data, ui, ci)
coefficients = myObj$par

where theta contains the starting values for par[1] and par[2], f is the function which is to be minimized, grad the gradient of f or NULL, data contains the data frame dat, ui is the constraint matrix, and ci the constraint vector.

$\endgroup$
0
$\begingroup$

The fastest option (faster than the other solutions posted here) is to observe that nnls can be reformulated as a quadratic programming problem, and that when written as a quadratic programming problem, one can apply individual constraints. Such problems in R can be fit using the quadprog package.

Here is a minimal example with a model where in this case the simulated intercept is nonzero:

# simulate some data
n = 100
p = 20
X = cbind(1,matrix(rnorm(n*p, mean=10, sd=1), nrow=n, ncol=p))
p_nonzero = 10 # nr of nonzero coefficients
SN = 10 # signal to noise ratio of nonzero elements
beta.intercept = -1 # coefficient for intercept 
beta.true = c(beta.intercept, rep(SN, p_nonzero), rep(0, p-p_nonzero))
y.true = X %*% beta.true
y = matrix(rnorm(n, mean=as.vector(y.true), sd=1),ncol=1)

First let's do a regular nnls fit with nonnegativity constraints on all coefficients - this then gives the wrong estimate for the intercept:

library(nnls)
nnls(A=X, b=y)$x
# [1]  0.00000000  9.96888742  9.95706223 10.02573854 10.05833035  9.95739332  9.97245103 10.06452289 10.05171408  9.89976862
# [11]  9.82506371  0.00000000  0.00000000  0.00000000  0.00000000  0.03724875  0.09294024  0.00000000  0.00000000  0.00000000
# [21]  0.00000000

Now using the equivalent quadratic programming reformulation with nonnegativity constraints applied only on some coefficients (in the function below specified by vector nonnegative to indicate for which coefficients the nonnegativity constraints apply) :

constrainedLS <- function(y, X, nonnegative=rep(TRUE,ncol(X))) {
  require(quadprog)
  beta_min=nonnegative*1-1
  beta_min[!nonnegative]=-.Machine$double.xmax
  solve.QP(dvec = crossprod(X,y), bvec = beta_min, Dmat=crossprod(X,X), Amat=diag(ncol(X)))$solution
}

If we apply nonnegativity constraints to all coefficients then we get approx. the same solution as nnls:

constrainedLS (y, X) # approximately same solution as nnls
# [1]  1.688868e-15  9.968887e+00  9.957062e+00  1.002574e+01  1.005833e+01  9.957393e+00  9.972451e+00  1.006452e+01  1.005171e+01
# [10]  9.899769e+00  9.825064e+00 -7.628210e-18 -2.152059e-18 -1.118194e-17  1.833571e-17  3.724875e-02  9.294024e-02  0.000000e+00
# [19]  2.234559e-17 -2.230120e-19  4.129302e-17

If we lift the nonnegativity constraint on the intercept we get the intercept in this example correctly estimated:

constrainedLS (y, X, nonnegative=c(FALSE,rep(TRUE,ncol(X)-1)))
# [1] -6.593426e+00  1.004423e+01  1.001232e+01  1.005849e+01  1.009688e+01  1.002808e+01  1.002120e+01  1.010366e+01  1.009439e+01
# [10]  9.939715e+00  9.867267e+00  2.456823e-18  0.000000e+00 -1.433482e-17  1.430082e-17  1.146234e-01  1.493929e-01  3.987769e-02
# [19] -5.454521e-18  3.872165e-18 -1.005829e-17
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.