1 item has been added to your cart.
| 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.
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 1000 number 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.
