How best to estimate regression parameters subject to constraints?

Question

The setting

Consider a least squares model for $y$ as a function of $x,$ possibly nonlinear in the parameters. Abstractly this can be expressed as

$$y = f(x;\theta) + \varepsilon$$

with the usual assumptions: the $\varepsilon$ are iid zero-mean errors of finite variance $\sigma^2.$

Further suppose this model represents a well-established mathematical or physical law that constrains the "reasonable" values of the parameter $\theta.$

As a concrete example, consider a standard model in radiological measurements such as a (simplified) Lorentzian peak at $x=0$ added to a constant background:

$$f(x,(\alpha,\beta,\gamma)) = \frac{\alpha}{\beta + x^2} + \gamma.$$

All three parameters $\theta = (\alpha,\beta,\gamma)$ are of interest and should be estimated as accurately as possible. However, the only physically realistic values of the parameters are non-negative (and $\beta$ is nonzero). Here is an illustration.

Generally, suppose there is a known nontrivial set $\Theta$ of realistic values of $\theta.$ I will refer to $\Theta$ as the constraints.

A fundamental question is whether the peak is present at all. That can be tested with the null hypothesis $H_0:\alpha=0,$ a situation that often holds in reality (maybe the peak is produced by a contaminant that usually is not present, for instance).

When fitting such a model, it is possible for any of the parameter estimates $\hat\theta$ to lie outside the constraints. (That is, $\hat\theta\notin \Theta.$) For instance, $\hat\alpha$ might be negative when (a) the true value $\alpha$ is much smaller than the error standard deviation $\sigma$ and (b) due to random (bad) luck, it looks like the peak is inverted. This happens frequently (about half the time) when $\alpha$ truly is zero.

The question

One has two options in these circumstances:

1. Use whatever estimates $\hat\theta$ the least squares procedure provides, even if they lie outside the constraints.

2. Force the least squares solution to conform to the constraints.

There are cogent arguments for and against both options: (1) can yield unrealistic estimates, but when applied to a large number of datasets might provide better average estimates on the whole; whereas (2) is forced to yield realistic estimates, but that can introduce a bias. For instance, forcing $\alpha$ to be nonnegative will produce positive nonzero estimates from time to time even when $\alpha$ is truly zero. In general, we're considering the possibility $H_0: \theta\in\partial\Theta$ that the parameter lies on the boundary of the constraints and wish to test the alternative $H_A:\theta\in \Theta \setminus \partial\Theta$ (the interior of the constraint space).

Is one of these options preferable? If so, which one and why? Is one of these options so problematic that it should not be used? If so, which one and why?

I also invite suggestions for additional options (such as regularization, penalized least squares, and so on), provided they can be shown superior to both (1) and (2).

---

I welcome answers from varying perspectives (classical, Bayes, etc.) and have enough bounty points to reward multiple good answers ;-).

Answer 1

Is one of these options preferable? If so, which one and why? Is one of these options so problematic that it should not be used? If so, which one and why?

Depends on your goal. If you trust your priors, model specification and error distribution, you would like to include as much of those prior beliefs in the model as you can.

If you want to perform an analysis on a large number of samples, on the other hand, you might want to perform a different analysis altogether. Of course, if you need to find out whether there is some systematic bias in your measurements, you need to use an estimator that has minimal bias -- but then could choose to model the bias directly.

So in general, I would dismiss your pro towards the method (1): can yield unrealistic estimates, but when applied to a large number of datasets might provide better average estimates on the whole, and only worry about the bias if it is the quantity you want to model.

Frequentist vs. Bayesian interpretation

There was a paper I found comparing frequentist constraints and bayesian priors, the difference of which comes down to interpretation as well as practical usage. Since the frequentist approach treats the parameter $\theta$ as unknown, but fixed, element of $\Theta$, it does not make much sense to me to allow values that are physically impossible -- i.e. there is no sensible interpretation for those. In practice, then, this means that the estimation must be restricted to a set of parameters that satisfy the constraints, if you wish to set them. Alternatively, you simply accept that your model doesn't actually model the phenomena.

In the Bayesian framework, the parameter itself is a random variable, for which we know at least that $\theta \in \Theta$. Of course, there are infinitely many distributions that can be used to assign total probability of 1 to that space. The interpretation of prior distributions to restrict parameter values is straightforward, and you get the (posterior) probability distribution of the parameter values, which, again, has quite simple interpretation. It could include physically impossible values, but in that case it is clear that they were, also, included in your prior beliefs.

I would say, then, that in the Bayesian case it is quite OK to not restrict the parameter to a certain range, and this has a simple explanation: We chose to allow the parameter to take those values, even if we knew that they were impossible. Conversely if you choose to restrict the parameter space by using a prior, the meaning of that is clear and reflects the analyst's prior beliefs: There is zero probability that the parameter could take these values.

If you are concerned about systematic biases introduced by the prior, one option would be to use cross-validation to check the distribution of the errors.

Bayesian approach in practice

Most common (and trivial) case where you restrict the parameter space by a prior functions is when setting a prior for variance of normal distribution: Inverse gamma distribution is often used. For most parameters there is often enough data that accurate prior distributions matter little, and noninformative, computationally easy (e.g. conjugate) priors are often chosen.

More relevant to your question, perhaps, is this is done with physical models in practice. One paper that might be of interest is A Bayesian Search for the Higgs Particle, where one distribution of the prior is based on the theoretical, expected mass of the Higgs Boson -- and another on what we might see if it does not exist.

The same author also had another paper that is more general and thus maybe more relevant to your question.

Answer 2

Option 1 is better. Here's why:

My favorite professor (Herman Friedman) used to say ``if you're not surprised, you haven't learned anything.''

Option 2 is a way of telling the data not to surprise you. It's the antithesis of science.

But what if option 1 gives you surprising values? Then you figure out why. You learn something. Maybe it's something mundane (often, data entry errors or coding errors) and you just fix the problem. But once in a wonderful while, it isn't mundane and you discover something. The model is wrong. Or, at least, wrong in your application.

Suppose you are interested in the heights of men and women. Maybe you have some theory about covariates. But you know men are taller than women, so you constrain the estimates. But what if you discovered some subgroup of people where men are not taller? That's really, really interesting. And someone else will find it, not you, and there goes your great discovery.

(I don't understand your example well enough to use it, but if it is really an absolute physical constraint, then you are in the first paragraph, and finding the weird estimates lets you discover that.)

Answer 3

During reading What are some motivations for using nonnegative least squares? I came across nonnegative least squares which imposes the constraint of all coefficients being nonnegative. The main argument seems to be that prior knowledge is encoded in the constraints and solutions which violate the constraints are known to be unreasonable.

Answer 4

One dimensional simple model

A simple model that illustrates a bit what is going one is a single parameter model

$$X_i \sim N(\theta,1) \quad \text{with the constraint $\theta \geq 0$}$$

For a sample of size $n$

• the sufficient statistic is $ T=\bar{X}$ (and we could use this as an unconstrained estimate $\hat{\theta}_{unconstrained} = T$)
• the constrained maximum likelihood estimate is $\hat{\theta}_{constrained} = \max(0,\bar{X})$.

Histograms for the distributions for these, when $n=5$ and $\theta = 0.5$, look like below. The distribution of the constrained estimate is a rectified Gaussian distribution.

That peak at $mle=0$ is a point mass and represents the tail from $T<0$ that is now closer to the true value $\theta$. In this case the boundaries are making the mean squared error smaller.

Variance, bias, and combining multiple experiments

The sufficient statistic $T$ has no bias but the variance is large and the mean squared error will be larger.

$$\text{mean}(T) = \theta\\ \text{var}(T) = \frac{1}{n}\\ MSE(T) = \frac{1}{n}$$

The constrained maximum likelihood estimate $\hat{\theta}_{constrained}$ has a bias, but due to a smaller variance the mean squared error will be smaller.

(for derivations see: Expectation and Variance of Gaussian going through Rectified Linear or Sigmoid function)

$$\text{mean}(\hat{\theta}_{constrained}) = \theta (1-p) + \frac{\phi(z)}{\sqrt{n}}\\ \text{var}(\hat{\theta}_{constrained}) = \theta^2 p(1-p) + \frac{\theta}{\sqrt{n}} \phi(z)(2p-1) + \frac{1}{n} (1-p-\phi(z)^2) \\ MSE(\hat{\theta}_{constrained}) = \frac{1}{n} (1-p) - \frac{\theta}{\sqrt{n}} \phi(z) + \theta^2 p$$

Where $z = -\theta \sqrt{n}$ and $p = \Phi(z)$ and $\phi$ and $\Phi$ are the standard Gaussian pdf and cdf.

The issue/consideration for the one dimensional simple model is

(1) can yield unrealistic estimates, but when applied to a large number of datasets might provide better average estimates on the whole;

When we have multiple estimates then the variance will reduce and the error of the average of estimates will approach the bias.

So when we use averages of estimates then using biased estimates can be a disadvantage, albeit these estimates having a smaller error when compared individually.

Depending on how the model is used we may prefer methods 1 or 2.

• If we are not using averages then method 2 is better because the mean squared error is larger.
• If we will use some averaging of several estimates (and we are not able to use all the data at once making a single estimate instead of averaging multiple estimates) then model 1 might be better. It will depend on the ratio between variance and bias.

Method 2 might give a better estimate, but that estimate might not be the best statistic for further analysis. In the example above this has been made more salient by method 1 being equivalent to the sufficient statistic, and method 2 is effectively loosing information. The same principle might be true for other methods (like regularised regression) where some estimate has 'improved' properties, but results in a loss of information.