Search
>> Home >> Resources & support >> FAQs >> Fitting a linear regression with interval (inequality) constraints using nl

### How do I fit a linear regression with interval (inequality) constraints in Stata?

 Title Fitting a linear regression with interval (inequality) constraints using nl Author Isabel Canette, StataCorp Date October 2011; minor revisions April 2015

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

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 Coef. 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
Parameter c taken as constant term in model & ANOVA table

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.

### Example 1: Constraints of the form a > 0

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 Coef. 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
Parameter c taken as constant term in model & ANOVA table

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 Coef. Legend /lna -7.535992 _b[lna:_cons] /b .8350376 _b[b:_cons] /c -57.69477 _b[c:_cons]
Parameter c taken as constant term in model & ANOVA table . nlcom a: exp(_b[lna:_cons]) a: exp(_b[lna:_cons])
 mpg2 Coef. Std. Err. z P>|z| [95% Conf. Interval] a .0005335 .0001579 3.38 0.001 .0002241 .000843

### Example 2: Constraints of the form 0 < a < 1

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 Coef. 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
Parameter c taken as constant term in model & ANOVA table . nlcom a: invlogit(_b[lgta:_cons]) a: invlogit(_b[lgta:_cons])
 mpg2 Coef. Std. Err. z P>|z| [95% Conf. Interval] a .0005335 .0001579 3.38 0.001 .0002241 .000843

### Example 3: Constraints of the form -1 < a < 1

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 Coef. 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
Parameter c taken as constant term in model & ANOVA table . nlcom a: tanh(_b[atanha:_cons]) a: tanh(_b[atanha:_cons])
 mpg2 Coef. Std. Err. z P>|z| [95% Conf. Interval] a .0005335 .0001579 3.38 0.001 .0002241 .000843

### Example 4: Constraints of the form 0 < a < b

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 Coef. 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
Parameter c taken as constant term in model & ANOVA table . nlcom (a: exp(_b[lna:_cons])) (b: exp(_b[lna:_cons])+exp(_b[lndiff:_cons])) a: exp(_b[lna:_cons]) b: exp(_b[lna:_cons])+exp(_b[lndiff:_cons])
 mpg2 Coef. 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

### Example 5: Constraints of the form k1.a < k2.b < k2.c

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 Coef. 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
Parameter lnd2 taken as constant term in model & ANOVA table . nlcom (a: _b[a:_cons]) /// > (b: 0.1*(exp(_b[lnd1:_cons])+0.5*_b[a:_cons])) /// > (c: 0.5*(exp(_b[lnd2:_cons]) + /// > exp(_b[lnd1:_cons])+0.5*_b[a:_cons])) a: _b[a:_cons] b: 0.1*(exp(_b[lnd1:_cons])+0.5*_b[a:_cons] ) c: 0.5*(exp(_b[lnd2:_cons]) + exp(_b[lnd1:_cons])+0.5*_b[a:_cons])
 turn Coef. 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

### Example 6: Setting a model where parameters are proportions

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 Coef. 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
. display "sum of coefficients = " _b[x1] + _b[x2] + _b[x3] sum of coefficients = 1.0825702

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 Coef. 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
Parameter a4 taken as constant term in model & ANOVA table . local na2 exp(_b[t2:_cons])/(1+exp(_b[t2:_cons])+exp(_b[t3:_cons])) . local na3 exp(_b[t3:_cons])/(1+exp(_b[t2:_cons])+exp(_b[t3:_cons])) . local na1 1/(1+exp(_b[t2:_cons])+exp(_b[t3:_cons])) . nlcom (a1: na1') (a2: na2') (a3: na3') a1: 1/(1+exp(_b[t2:_cons])+exp(_b[t3:_cons])) a2: exp(_b[t2:_cons])/(1+exp(_b[t2:_cons])+exp(_b[t3:_cons])) a3: exp(_b[t3:_cons])/(1+exp(_b[t2:_cons])+exp(_b[t3:_cons]))
 y Coef. 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.