Linear Regression with individual constraints in R 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?
The nnls package also offers the possibility to enforce some coefficients to be negative and others to be positive via the function nnnpls(). However, I only want the positive constraint for some parameters, the other ones can be either positive or negative (i.e., unconstrained).
 A: 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)

A: 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 or the osqp package. As the osqp package is the fastest, and allows the covariate matrix to be either dense or sparse, I'll use that one here.
Using that package, here is a function that does a constrained least squares fit using vectors of lower and upper bounds lower and upper on the fitted coefficients and allowing X to be either dense or sparse:
constrainedLS_osqp <- function(y, X,
                               lower=rep(0, ncol(X)), upper=rep(Inf, ncol(X)),
                               x.start = NULL, y.start = NULL) {
  require(osqp)
  require(Matrix)
  XtX = crossprod(X,X)
  Xty = crossprod(X,y)
  
  settings = osqpSettings(verbose = FALSE, eps_abs = 1e-8, eps_rel = 1e-8, linsys_solver = 0L,
                          warm_start = FALSE)
  pff = .sparseDiagonal(ncol(X))
  model <- osqp(XtX, -Xty, pff, l=lower, u=upper, pars=settings)
  
  if (!is.null(x.start)) model$WarmStart(x=x.start, y=y.start)

  coefs = model$Solve()$x # fitted coefficients

  coefs = pmax(lower, pmin(coefs, upper) ) # fitted coefficients sometimes go very slightly outside constraint zone due to numerical inaccuracies in solver - this is fixed here via clipping

  return(coefs)
}

A: 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.
