Estimating error for parameters from multiple regression with linear constraints

Question

I am working on a multiple linear regression problem where I would like to constrain only some of the parameters to non-negative values. There have been discussions of how to solve for the parameters posted on SE here and on another site. Both solutions used optim() or constrOptim() to minimize the residual sum of squares, which worked very well for me and gives the same values as lm() for the unconstrained version of the problem.

But I would like now to take this a step further and estimate standard errors of those parameters, similar to what I'd get if I used lm(). However, there is nothing in the optim() object that would suggest a route to estimating errors, and I haven't found anything in the control settings that would suggest a route either. So I'm a little bit at a loss as to how to proceed. My hunch is that this approach - an optimization problem that assumes the data values are fixed and the parameter solutions represent a global minimum - does not allow for error. Is there any validity to that, or am I just missing something basic?

This is a small reproducible example, adapted from the example provided in the first link to work with the optim() function, rather than constrOptim().

min.RSS <- function(data, par){
      with(data, sum((par[1]*x1 + par[2]*x2 - y)^2))
    }
dat = data.frame(x1=c(1,2), x2=c(2,3), y=c(5,6))
result = optim(par=c(0,1), min.RSS, data=dat, method="L-BFGS-B", 
                lower=c(-Inf,0), upper=c(0,Inf))
result$convergence
result$par

Answer 1

I am assuming that you have really good reasons for constraining the coefficients and have thought hard about the potential problems with doing so. Note that a linear model is almost always an approximation to reality, and as a result coefficients that "ought to be non-negative" may end up negative in a linear model based on data with correlated predictors. Forcing coefficients to be non-negative poses a serious risk of generating a suboptimal model.

One way to approach your problem would be to repeat your analysis on multiple bootstrap samples from your data set, and use the distributions of the coefficients to provide the desired estimates. The boot package in R is one tool for this approach, although it might take some effort to write the needed function properly.

If you look carefully at the results of your bootstrap analyses, you will almost certainly find that the mean values, over those bootstrapped samples, of coefficients you constrained will differ from their values in your original analysis. That's because constraining non-negativity on a coefficient that might be negative will probably introduce bias in the technical statistical sense. You will need to consider how you will deal with that eventuality, based on your understanding of the subject matter.