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