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.

The Stata Blog: Not Elsewhere Classified Find us on Facebook Follow us on Twitter LinkedIn Google+ Watch us on YouTube