 »  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

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   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

Note: Parameter c is used as a constant term during estimation.

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   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

Note: Parameter c is used as a constant term during estimation.

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]

Note: Parameter c is used as a constant term during estimation.

. nlcom a: exp(_b[lna:_cons])

a:  exp(_b[lna:_cons])

mpg2   Coefficient  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   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

Note: Parameter c is used as a constant term during estimation.

. nlcom  a: invlogit(_b[lgta:_cons])

a:  invlogit(_b[lgta:_cons])

mpg2   Coefficient  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   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

Note: Parameter c is used as a constant term during estimation.

. nlcom  a: tanh(_b[atanha:_cons])

a:  tanh(_b[atanha:_cons])

mpg2   Coefficient  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   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

Note: Parameter c is used as a constant term during estimation.

. 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   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



### 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   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

Note: Parameter lnd2 is used as a constant term during estimation.

. 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   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



### 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
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

. 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   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

Note: Parameter a4 is used as a constant term during estimation.

. 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   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.