Title | Fitting a linear regression with interval (inequality) constraints using nl | |

Author | Isabel Canette, StataCorp |

If you need to fit a linear model with linear constraints, you can use the
Stata command
**cnsreg**.
If you need to fit a nonlinear
model with interval constraints, you can use the
**ml**
command, as explained in the FAQ
How do I fit a regression with interval (inequality) constraints in Stata?
However, if you have a linear regression, the simplest way to include these
kinds of constraints is by using the
**nl** command.

- Introduction
- Example 1: Constraints of the form a > 0
- Example 2: Constraints of the form 0 < a < 1
- Example 3: Constraints of the form -1 < a < 1
- Example 4: Constraints of the form 0 < a < b
- Example 5: Constraints of the form k1.a < k2.b < k3.c
- Example 6: Setting a model where parameters are proportions

First, let's review how to fit a linear regression using
**nl**. We will use this command to fit a regression of
**mpg2** on **price** and **turn**:

. sysuse auto(1978 Automobile Data). generate mpg2 = -mpg . nl (mpg2 = {a}*price + {b}*turn + {c})(obs = 74) Iteration 0: residual SS = 1016.186 Iteration 1: residual SS = 1016.186

Source | SS df MS | ||

Number of obs = 74 | |||

Model | 1427.2735 2 713.636766 | R-squared = 0.5841 | |

Residual | 1016.1859 71 14.3124778 | Adj R-squared = 0.5724 | |

Root MSE = 3.783184 | |||

Total | 2443.4595 73 33.4720474 | Res. dev. = 403.8641 |

mpg2 | Coefficient Std. err. t P>|t| [95% conf. interval] | |

/a | .0005335 .0001579 3.38 0.001 .0002187 .0008483 | |

/b | .8350376 .1058498 7.89 0.000 .6239791 1.046096 | |

/c | -57.69477 4.027951 -14.32 0.000 -65.72627 -49.66326 | |

We will set inequality constraints (and interval constraints) via
transformations. For example, let's assume that we want parameter
**a** to be positive. This can be achieved by expressing
**a** as an exponential. We can, therefore, estimate
**lna = ln(a)** and then recover **a = exp(lna)
** after the estimation. The trick is to use a transformation whose
range is the interval over which we want to restrict the parameter.

Type

. help math functions

to see the mathematical functions available in Stata.

There are many ways to set interval constraints. The following examples show some possibilities.

As stated before, we will estimate **ln(a)** instead of
**a**.

. nl (mpg2 = exp({lna})*price + {b}*turn + {c}), nolog(obs = 74)

Source | SS df MS | ||

Number of obs = 74 | |||

Model | 1427.2735 2 713.636766 | R-squared = 0.5841 | |

Residual | 1016.1859 71 14.3124778 | Adj R-squared = 0.5724 | |

Root MSE = 3.783184 | |||

Total | 2443.4595 73 33.4720474 | Res. dev. = 403.8641 |

mpg2 | Coefficient Std. err. t P>|t| [95% conf. interval] | |

/lna | -7.535992 .2959172 -25.47 0.000 -8.126034 -6.94595 | |

/b | .8350376 .1058498 7.89 0.000 .6239791 1.046096 | |

/c | -57.69477 4.027951 -14.32 0.000 -65.72627 -49.66326 | |

The output shows the parameter **lna**. To recover
**a**, we can use the **nlcom** command; we can
always call **nl** (or any estimation command) with the
**coeflegend** option to see how to refer to the parameters in
further expressions.

. nl, coeflegend

Source | SS df MS | ||

Number of obs = 74 | |||

Model | 1427.2735 2 713.636766 | R-squared = 0.5841 | |

Residual | 1016.1859 71 14.3124778 | Adj R-squared = 0.5724 | |

Root MSE = 3.783184 | |||

Total | 2443.4595 73 33.4720474 | Res. dev. = 403.8641 |

mpg2 | Coefficient Legend | |

/lna | -7.535992 _b[lna:_cons] | |

/b | .8350376 _b[b:_cons] | |

/c | -57.69477 _b[c:_cons] | |

mpg2 | Coefficient Std. err. z P>|z| [95% conf. interval] | |

a | .0005335 .0001579 3.38 0.001 .0002241 .000843 | |

Because the range of the inverse logit function is the interval (0,1), we
can use the Stata function **invlogit()** to set this
restriction.

. nl (mpg2 = invlogit({lgta})*price + {b}*turn + {c}), nolog(obs = 74)

Source | SS df MS | ||

Number of obs = 74 | |||

Model | 1427.2735 2 713.636766 | R-squared = 0.5841 | |

Residual | 1016.1859 71 14.3124778 | Adj R-squared = 0.5724 | |

Root MSE = 3.783184 | |||

Total | 2443.4595 73 33.4720474 | Res. dev. = 403.8641 |

mpg2 | Coefficient Std. err. t P>|t| [95% conf. interval] | |

/lgta | -7.535459 .2960752 -25.45 0.000 -8.125816 -6.945101 | |

/b | .8350376 .1058498 7.89 0.000 .6239791 1.046096 | |

/c | -57.69477 4.027951 -14.32 0.000 -65.72627 -49.66326 | |

mpg2 | Coefficient Std. err. z P>|z| [95% conf. interval] | |

a | .0005335 .0001579 3.38 0.001 .0002241 .000843 | |

We can use the hyperbolic tangent function for constraints like this.

. nl (mpg2 = tanh({atanha})*price + {b}*turn + {c}), nolog(obs = 74)

Source | SS df MS | ||

Number of obs = 74 | |||

Model | 1427.2735 2 713.636766 | R-squared = 0.5841 | |

Residual | 1016.1859 71 14.3124778 | Adj R-squared = 0.5724 | |

Root MSE = 3.783184 | |||

Total | 2443.4595 73 33.4720474 | Res. dev. = 403.8641 |

mpg2 | Coefficient Std. err. t P>|t| [95% conf. interval] | |

/atanha | .0005335 .0001579 3.38 0.001 .0002187 .0008483 | |

/b | .8350376 .1058498 7.89 0.000 .6239791 1.046096 | |

/c | -57.69477 4.027951 -14.32 0.000 -65.72627 -49.66326 | |

mpg2 | Coefficient Std. err. z P>|z| [95% conf. interval] | |

a | .0005335 .0001579 3.38 0.001 .0002241 .000843 | |

We can express **a** as an exponential, as in Example 1, to
ensure that it will be positive. In addition, we want to set the restriction
**b>a**; therefore, we can also express the difference
**b−a** as an exponential.

. nl (mpg2 = exp({lna})*price + (exp({lndiff})+exp({lna}))*turn + {c}), nolog(obs = 74)

Source | SS df MS | ||

Number of obs = 74 | |||

Model | 1427.2735 2 713.636766 | R-squared = 0.5841 | |

Residual | 1016.1859 71 14.3124778 | Adj R-squared = 0.5724 | |

Root MSE = 3.783184 | |||

Total | 2443.4595 73 33.4720474 | Res. dev. = 403.8641 |

mpg2 | Coefficient Std. err. t P>|t| [95% conf. interval] | |

/lna | -7.535992 .2959172 -25.47 0.000 -8.126035 -6.94595 | |

/lndiff | -.1809177 .1269002 -1.43 0.158 -.4339496 .0721142 | |

/c | -57.69477 4.027951 -14.32 0.000 -65.72627 -49.66326 | |

mpg2 | Coefficient Std. err. z P>|z| [95% conf. interval] | |

a | .0005335 .0001579 3.38 0.001 .0002241 .000843 | |

b | .8350376 .1058498 7.89 0.000 .6275758 1.042499 | |

The concepts in Example 4 can be extended to similar cases, like
**a<b<c** or **a<b<2c**.

For example, let's assume that we want to fit the regression

nl (turn = {a}*headroom + {b}*displacement + {c})

with the constraints **0.5a<10b<2c**.

These two inequalities can be presented as

10b - 0.5a > 0 2c - 10b > 0

and we can express the left-hand sides of these as exponentials to ensure that they will turn out positive. We can then estimate the two following parameters

d1 = exp(lnd1) = 10b - 0.5a d2 = exp(lnd2) = 2c - 10b

which imply that we can substitute **b** and **c**
as follows:

b = 0.1(d1 + 0.5a) c = 0.5(d2 + 10b)

Hence, our command line would look like this:

. nl (turn = {a}*headroom + 0.1*(exp({lnd1})+0.5*{a})*displacement>+ 0.5*(exp({lnd2}) + 10*0.1*(exp({lnd1})+0.5*{a}))), nolog(obs = 74)

Source | SS df MS | ||

Number of obs = 74 | |||

Model | 858.16801 2 429.084007 | R-squared = 0.6074 | |

Residual | 554.69685 71 7.8126317 | Adj R-squared = 0.5963 | |

Root MSE = 2.795109 | |||

Total | 1412.8649 73 19.3543132 | Res. dev. = 359.0653 |

turn | Coefficient Std. err. t P>|t| [95% conf. interval] | |

/a | .3751308 .4392973 0.85 0.396 -.5008032 1.251065 | |

/lnd1 | -1.782972 1.436274 -1.24 0.219 -4.64682 1.080877 | |

/lnd2 | 4.137724 .0388728 106.44 0.000 4.060214 4.215234 | |

turn | Coefficient Std. err. z P>|z| [95% conf. interval] | |

a | .3751308 .4392973 0.85 0.393 -.4858761 1.236138 | |

b | .0355703 .0040468 8.79 0.000 .0276388 .0435018 | |

c | 31.50786 1.214813 25.94 0.000 29.12687 33.88885 | |

Finally, let’s see how to fit a model where the coefficients are proportions; that is, they are all positive and add up to one.

We will fit the linear model

y = a1*x1 + a2*x2 + a3*x3 + a4 + ε

where **a1**, **a2**, and **a3** are
positive, and **a1 + a2 + a3 = 1**.

We will use the transformation implemented in the Stata command
**mlogit**:

a2 = exp(t2)/(1+exp(t2)+exp(t3)) a3 = exp(t3)/(1+exp(t2)+exp(t3)) a1 = 1/(1+exp(t2)+exp(t3))

Here we illustrate the concept with simulated data:

. clear . set seed 12345 . set obs 1000number of observations (_N) was 0, now 1,000. generate x1 = invnormal(runiform()) . generate x2 = invnormal(runiform()) . generate x3 = invnormal(runiform()) . generate ep = invnormal(runiform()) . generate y = .2*x1 + .5*x2 + .3*x3 + 1 + ep

Although the actual coefficients used for the simulation add up to one, the
estimates obtained with **regress** most likely will not.

. regress y x1 x2 x3

Source | SS df MS | Number of obs = 1,000 | |

F(3, 996) = 145.31 | |||

Model | 432.789199 3 144.263066 | Prob > F = 0.0000 | |

Residual | 988.844775 996 .992816039 | R-squared = 0.3044 | |

Adj R-squared = 0.3023 | |||

Total | 1421.63397 999 1.42305703 | Root MSE = .9964 |

y | Coefficient Std. err. t P>|t| [95% conf. interval] | |

x1 | .2536189 .0315703 8.03 0.000 .1916669 .3155708 | |

x2 | .507397 .0317749 15.97 0.000 .4450436 .5697504 | |

x3 | .3215543 .0314128 10.24 0.000 .2599115 .3831972 | |

_cons | .99246 .0315226 31.48 0.000 .9306016 1.054318 | |

Let’s fit the model by setting the restrictions using
**nl**:

. local ma2 (exp({t2})/(1+exp({t2})+exp({t3}))) . local ma3 (exp({t3})/(1+exp({t2})+exp({t3}))) . local ma1 (1/(1+exp({t2})+exp({t3}))) . nl (y = `ma1'*x1 + `ma2'*x2 + `ma3'*x3 + {a4}), delta(1e-7) nolog(obs = 1000)

Source | SS df MS | ||

Number of obs = 1,000 | |||

Model | 430.4658 2 215.232898 | R-squared = 0.3028 | |

Residual | 991.16818 997 .99415063 | Adj R-squared = 0.3014 | |

Root MSE = .997071 | |||

Total | 1421.634 999 1.42305703 | Res. dev. = 2829.006 |

y | Coefficient Std. err. t P>|t| [95% conf. interval] | |

/t2 | .7524281 .1506533 4.99 0.000 .4567942 1.048062 | |

/t3 | .2617823 .1748548 1.50 0.135 -.0813434 .604908 | |

/a4 | .9913625 .0315356 31.44 0.000 .9294787 1.053246 | |

y | Coefficient Std. err. z P>|z| [95% conf. interval] | |

a1 | .2261732 .0259945 8.70 0.000 .175225 .2771214 | |

a2 | .4799727 .0262525 18.28 0.000 .4285188 .5314266 | |

a3 | .2938541 .025686 11.44 0.000 .2435104 .3441978 | |

If you find convergence issues when using **nl** to solve these
problems, you might try using **ml** instead, as explained in
stata.com/support/faqs/statistics/regression-with-interval-constraints.

Using **ml** would allow you to specify analytical derivatives
and to have better control of your optimization process.