Search
   >> Home >> Resources & support >> FAQs >> Fitting a regression with interval constraints

How do I fit a regression with interval constraints in Stata?

Title   Fitting a regression with interval constraints
Author Isabel Canette, StataCorp
Date October 2006; updated June 2013

Note:   The examples in this tutorial are for illustration purposes only. If you need to fit a linear regression with interval constraints, use the Stata command nl as explained in How do I fit a linear regression with interval (inequality) constraints in Stata?

If you are interested in fitting a linear regression with a linear constraint, you can use the command cnsreg. However, sometimes you may need to fit a model for which one of the coefficients should be positive or between zero and one. You can fit such a model by reparameterizing the coefficients. This tutorial explains the steps to fit this model by using the Stata command ml. For simplicity, I explain how to set interval constraints on a linear regression. However, this method could be applied to other kinds of regression (e.g., probit, logit, Poisson, tobit), provided that they can be fitted by maximum likelihood. As stated above, if you want to set nonlinear constraints on a linear regression, you can use nl, which would not require any programming. Starting with Stata 13, the mlexp command allows you to fit maximum likelihood models by evaluating an expression. At the end of each section, I include the syntax to fit the models using mlexp.

This FAQ is organized as follows:

  1. Introduction
  2. Prerequisites
  3. Our introductory example
  4. Application: Setting a model where parameters are proportions
  5. Final remarks

1. Introduction

I will start by presenting an example on how to use ml to fit a linear regression. The code in this example can be modified to set interval constraints. I will explain how to perform this modification by showing the steps to restrict one coefficient to be positive.

Finally, I will present an example where the coefficients are restricted to be proportions; that is, they are constrained to be positive numbers that add up to one.

If you know (from the nature of the problem) that the constraints you want to set are true, the following method will allow you to set these restrictions in the model. However, if the data do not support the validity of these constraints, the algorithm for the restricted model probably won’t achieve convergence.

There are essentially two situations where we may want to set interval constraints. One case is when the numerical global maximum of the likelihood actually is where we want it to be, but we want to help the algorithm avoid difficult spots somewhere else. This case is similar to the situation of estimating the standard deviation in a linear regression, where the maximum is positive, but we set the restriction to help the algorithm achieve convergence. The other situation is when we know that the true parameter is in a certain region, but for numerical reasons the global maximum of the likelihood may be somewhere else. We are interested not in finding a global maximum but in finding a maximum in a certain subset. This is the case of the example of proportions, considered in section 5.

If you also need to set linear constraints, you can use the option constraint of the command ml.

2. Prerequisites

This tutorial assumes that you are familiar with the command ml. I will review some characteristics of the estimation by using ml here. However, if you have never used the Stata command ml, see [R] ml.

Technical notes explain some details on the output; you can skip them in a first reading.

3. Our introductory example

Let’s assume that we want to reproduce the following estimation,

. sysuse auto, clear
. regress price mpg displacement

but we want to make sure that the coefficient of variable mpg will be negative.

Because we are going to use method d0, we will borrow the program mynormal_d0.ado from the book Maximum Likelihood Estimation with Stata, Fourth Edition, by Gould, Pitblado, and Poi (you can download this file from the book’s website). Although the syntax for method lf is simpler, if we adapt an lf program to fit a model with interval constraints, the good properties of lf will be lost (see technical note below); some preliminary simulations showed that, when adapted to our problem, method d0 is better for convergence.

Technical note:   Whenever the likelihood meets the linear form, method lf is efficient. With lf, the derivatives are computed at the equation level, and the chain rule is applied to obtain the derivatives with respect to the parameters. This method produces high-quality derivatives and therefore accurate estimators. When we modify a maximum likelihood program to meet interval restrictions according to the procedure explained below, the equations are disassembled; then, even when using an lf program, we will force ml to compute numerical derivatives directly with respect to the coefficients, losing the numerical benefits of lf.

3.1. Using method d0 to fit a linear regression

Let’s start by fitting the model without restrictions, using regress.

. sysuse auto, clear
(1978 Automobile Data)

. regress price mpg displacement

Source SS df MS Number of obs = 74
F( 2, 71) = 13.35
Model 173587098 2 86793549.2 Prob > F = 0.0000
Residual 461478298 71 6499694.33 R-squared = 0.2733
Adj R-squared = 0.2529
Total 635065396 73 8699525.97 Root MSE = 2549.4
price Coef. Std. Err. t P>|t| [95% Conf. Interval]
mpg -121.1833 72.78844 -1.66 0.100 -266.3193 23.95276
displacement 10.50885 4.58548 2.29 0.025 1.365658 19.65203
_cons 6672.766 2299.72 2.90 0.005 2087.254 11258.28

Now let’s fit the same model by using ml:

/* begin do file */
cscript
program mynormal_d0
        version 13.0
        args todo b lnf
        tempname lnsigma sigma
        tempvar mu 
        mleval `mu' = `b', eq(1)
        mleval `lnsigma' = `b', eq(2) scalar
        quietly {
                scalar `sigma' = exp(`lnsigma')
                mlsum `lnf' = ln( normalden($ML_y1,`mu',`sigma') )
        }
end

sysuse auto 
ml model d0 mynormal_d0 (xb: price = mpg displacement) (lnsigma:), ///
diparm(lnsigma, exp label(sigma))  
ml maximize
/* end do file */

Here are the results:


                                                  Number of obs   =         74
                                                  Wald chi2(2)    =      27.84
Log likelihood = -683.89903                       Prob > chi2     =     0.0000

price Coef. Std. Err. z P>|z| [95% Conf. Interval]
xb
mpg -121.1831 71.29769 -1.70 0.089 -260.924 18.55784
displacement 10.50884 4.491567 2.34 0.019 1.705527 19.31215
_cons 6672.764 2252.621 2.96 0.003 2257.709 11087.82
lnsigma
_cons 7.82294 .0821995 95.17 0.000 7.661832 7.984048
sigma 2497.236 205.2714 2125.648 2933.782

Technical note: You will see some differences between the output of ml and the output of regress. First, the standard errors are different because a different number of degrees of freedom was used to compute the variance of the errors. The variance for the errors in regress and ml is estimated, respectively, as v1 and v2, where

v1 = e'e/(N-k)
v2 = e'e/(N)

Here e is the vector of residuals, N is the number of observations, and k is the number of regressors, including the constant.

Second, the confidence intervals reported by regress are based on the t distribution, and ml uses the Gaussian distribution to compute them.

Third, regress performs an F test, and ml performs a Wald test. See the following FAQ to learn how these tests relate to each other: Why does test sometimes produce chi-squared and other times F statistics?

Let’s have a closer look at the code for ml; the line

ml model d0 mynormal_d0 (xb: price = mpg displacement) (lnsigma:)

introduces two equations. The first equation,

(xb: price = mpg displacement)

indicates the following:

The second equation,

(lnsigma:)

indicates that

Because this is a constant-only equation, there is no need to generate a whole variable for it in the evaluation. Evaluating it as scalar is more efficient, as in

mleval `lnsigma' = `b', eq(2) scalar

This is essentially the same as

scalar `lnsigma' = `b'[1,4]

This parameter is used in the line

scalar double `sigma' = exp(`lnsigma')

We are not estimating the standard deviation (sigma), but its logarithm (lnsigma). This trick allows us to express the standard deviation as an exponent (exp(`lnsigma')) and guarantees that sigma will be positive. We make sure that the algorithm will not search for sigma in the negative subspace. We want to use this trick to apply an inequality constraint to a coefficient. The model above could also be fit by using the mlexp syntax:

mlexp (ln(normalden(price, {xb:mpg displacement}+{b0}, exp({lnsigma}))))

3.2. Setting separate equations for the coefficients

In this step, we set up a separate equation for the coefficient that we will restrict.

In the program mynormal_d0, the coefficient of mpg is included in equation 1. To set a restriction on it, we need to extract this coefficient from the equation. More precisely, we need to include this coefficient as a scalar parameter in the model, that is

ml model d0 mynormal_b (a: ) (xb: price = displacement ) (lnsigma:)

where the constant-only equation (a: ) corresponds to the coefficient of mpg. The question of how to include the variable itself remains; as we cannot introduce this variable in any of the equations, we will need to use a global macro. Look at the following code:

/* begin do file */
cscript
program mynormal_b
        version 13.0
        args todo b lnf
        tempname a lnsigma sigma
        tempvar xb mu 
        mleval `a' = `b', eq(1) scalar
        mleval `xb' = `b', eq(2)
        mleval `lnsigma' = `b', eq(3) scalar
        quietly {
                generate double `mu' = `xb' +`a'*$x2
                scalar `sigma' = exp(`lnsigma')
                mlsum `lnf' = ln( normalden($ML_y1,`mu',`sigma') )
        }
end

sysuse auto
global x2  mpg

ml model d0 mynormal_b (a: ) (xb: price = displacement ) (lnsigma:), ///
diparm(lnsigma, exp label(sigma))
ml maximize
/* end do file */

Here are the estimation results:

                                                  Number of obs   =         74
                                                  Wald chi2(0)    =          .
Log likelihood = -683.89903                       Prob > chi2     =          .

price Coef. Std. Err. z P>|z| [95% Conf. Interval]
a
_cons -121.1833 71.29771 -1.70 0.089 -260.9242 18.55769
xb
displacement 10.50885 4.491569 2.34 0.019 1.705532 19.31216
_cons 6672.766 2252.621 2.96 0.003 2257.709 11087.82
lnsigma
_cons 7.82294 .0821995 95.17 0.000 7.661832 7.984048
sigma 2497.237 205.2716 2125.649 2933.784

Technical note:   When we use mynormal_b, missing values are reported for the Wald chi2 test. ml, by default, reports this statistic for all the coefficients except the intercept in the first equation; the number of degrees of freedom is reported as zero because the first equation has only the intercept.

To perform a Wald test on the coefficients, we can use the command test after the estimation:

. matrix list e(b) 

e(b)[1,4]
               a:           xb:           xb:      lnsigma:
           _cons  displacement         _cons         _cons
y1    -121.18326     10.508845     6672.7659     7.8229403

. test ([a]_cons) ([xb]displacement) 

 ( 1)  [a]_cons = 0
 ( 2)  [xb]displacement = 0

           chi2(  2) =   27.84
         Prob > chi2 =    0.0000

These are the values reported by ml when running mynormal_d0.

The model above could also be fit by using the following mlexp syntax:

mlexp (ln(normalden(price, {a}*mpg+{b1}*displacement+{b0}, exp({lnsigma}))))

3.3. Including the restrictions in the model

Now we only have to change one line in the program: we substitute `a' by -exp(`a') in the evaluation and display the transformed parameter after the estimation:

/* begin do file */ 
cscript
program mynormal_c
        version 13.0
        args todo b lnf
        tempname a  lnsigma sigma
        tempvar xb mu
        mleval `a' = `b', eq(1) scalar
        mleval `xb' = `b', eq(2)
        mleval `lnsigma' = `b', eq(3)
          quietly {
                generate double `mu' = `xb' -exp(`a')*$x2
                scalar `sigma' = exp(`lnsigma')
                mlsum `lnf' = ln( normalden($ML_y1,`mu',`sigma') )
        }
end


sysuse auto
global x2  mpg
ml model d0 mynormal_c (a: ) (xb: price = displacement ) (lnsigma:), ///
diparm(lnsigma, exp label(sigma))
ml init lnsigma:_cons=7.8
ml maximize
/* end do file */

Here we used the global macro x2 to incorporate the variable mpg; when including a variable manually (i.e., not through the standard syntax of ml), you must check for missing values.

Here are the results:

                                                  Number of obs   =         74
                                                  Wald chi2(0)    =          .
Log likelihood = -683.89903                       Prob > chi2     =          .

price Coef. Std. Err. z P>|z| [95% Conf. Interval]
a
_cons 4.797299 .58835 8.15 0.000 3.644154 5.950444
xb
displacement 10.50887 4.491572 2.34 0.019 1.705552 19.31219
_cons 6672.748 2252.625 2.96 0.003 2257.685 11087.81
lnsigma
_cons 7.82294 .0821995 95.17 0.000 7.661832 7.984048
sigma 2497.237 205.2717 2125.649 2933.784

Now we can use the command nlcom to display the estimate for the coefficient:

. nlcom -exp([a]_cons)

       _nl_1:  -exp([a]_cons)

price Coef. Std. Err. z P>|z| [95% Conf. Interval]
_nl_1 -121.1827 71.29782 -1.70 0.089 -260.9238 18.5585

We may want to help ml a little by providing initial values. For example, if we believe that the variance reported by regress is reasonable, we can include it as an initial value with ml init. Here we used the residual MS reported by regress to set the logarithm of the variance, that is

. display log(sqrt(6499694.33))    
7.8436329

The model above could also be fit by using the following mlexp syntax:

mlexp (ln(normalden(price, -exp({a})*mpg + {b1}*displacement +{b0}, exp({lnsigma}))))

4. Application: Setting a model where parameters are proportions

Now we will apply the previous ideas to a particular example. Let’s assume that we want to fit the model

       y = a1*x1 + a2*x2 + a3*x3 + a4

where a1, a2, and a3 are proportions, i.e., values between zero and one that add up to one. Given this restriction, we need to estimate only two parameters (because a1 = 1 − (a2 + a3)). We need to find a proper parameterization for our restriction. A nice choice is the parameterization used by the command mlogit that fits multinomial logit models.

In this approach, instead of estimating a2 and a3, we estimate t2 and t3, where:

       a2 = exp(t2)/(1+exp(t2)+exp(t3))
       a3 = exp(t3)/(1+exp(t2)+exp(t3))
       a1 = 1/(1+exp(t2)+exp(t3))

Verifying that this parameterization guarantees our restriction is straightforward.

In the following code, the function mynormal_c can be used with the command ml to perform this estimation; the wrapper regprop runs this program and displays the values of interest.

/* begin do file */
cscript
set seed 12345
set obs 1000
generate x1 = invnormal(uniform())
generate x2 = invnormal(uniform())
generate x3 = invnormal(uniform())
generate ep = invnormal(uniform())

generate y = .2*x1 + .5*x2 + .3*x3 + 1 + ep
regress y x1 x2 x3

program mynormal_c
	version 13.0
	args todo b lnf
	tempname t2 t3 a4 lnsigma sigma den a1 a2 a3 
	tempvar mu
	mleval `t2' = `b', eq(1) scalar
	mleval `t3' = `b', eq(2) scalar
	mleval `a4' = `b', eq(3) scalar
	mleval `lnsigma' = `b', eq(4)
	scalar `den' = 1+ exp(`t2') + exp(`t3')
	scalar `a1' =  1/`den'
	scalar `a2' = exp(`t2')/`den'
	scalar `a3'  = exp(`t3')/`den'
	generate double `mu' = `a1'*$x_001 + `a2'*$x_002 + `a3'*$x_003 + `a4'
	scalar `sigma' = exp(`lnsigma')
	mlsum `lnf' = ln( normalden($ML_y1,`mu',`sigma') )
	
end


program regprop
	version 13.0
	syntax varlist
	gettoken lhs rhs: varlist
	global x_001 : word 1 of `rhs'
	global x_002 : word 2 of `rhs'
	global x_003 : word 3 of `rhs'
  

        local f1 func(1/(1+exp(@1)+exp(@2)))
        local der1 der( -exp(@1)/(1+exp(@1)+exp(@2))^2 ///
              -exp(@2)/(1+exp(@1)+exp(@2))^2 ) 
        local f2  func(exp(@1)/(1+exp(@1)+exp(@2))) 
        local der2 der( exp(@1)*(1+exp(@2))/(1+exp(@1)+exp(@2))^2 ///
              -exp(@1)*exp(@2)/(1+exp(@1)+exp(@2))^2 ) 
        local f3 func(exp(@2)/(1+exp(@1)+exp(@2))) 
        local der3 der( -exp(@1)*exp(@2)/(1+exp(@1)+exp(@2))^2 /// 
              exp(@2)*(1+exp(@1))/(1+exp(@1)+exp(@2))^2 ) 
  
        ml model d0 mynormal_c (t2:) (t3:) (c: `lhs' = ) (lnsigma:), ///
        diparm( t2 t3, `f1' `der1' label ($x_001) ci(logit) ) ///
        diparm( t2 t3, `f2' `der2' label ($x_002) ci(logit) ) ///
        diparm( t2 t3, `f3' `der3' label ($x_003) ci(logit) )
        ml maximize 
end
  
regprop y x1 x2 x3
/* end do file */

Here are the results:


                                                  Number of obs   =       1000
                                                  Wald chi2(0)    =          .
Log likelihood = -1439.6632                       Prob > chi2     =          .

y Coef. Std. Err. z P>|z| [95% Conf. Interval]
t2
_cons .9853398 .1686424 5.84 0.000 .6548068 1.315873
t3
_cons .4539756 .1953043 2.32 0.020 .0711861 .8367651
c
_cons 1.03804 .0323145 32.12 0.000 .9747049 1.101376
lnsigma
_cons .0207247 .0223607 0.93 0.354 -.0231014 .0645508
x1 .1903572 .0260333 .1444568 .2466378
x2 .509914 .0264092 .4582315 .5613855
x3 .2997288 .0262758 .2508747 .3536058

Here the option diparm was used to display the actual estimates of the coefficients; besides aesthetic issues, in this case, this option is preferred to the command nlcom because of its flexibility. For example, the option ci(logit) allows us to make sure that the limits of the confidence intervals will also be between zero and one. The syntax for the option diparm is similar to the syntax for the command _diparm.

However, the confidence intervals reported by the option diparm are asymptotically equivalent to confidence intervals reported by nlcom; therefore, you can also use nlcom as in the example in section 1.3.

The model above could also be fit by using the following syntax:

local den  (1+ exp({t2}) + exp({t3}))
local a1 1/`den'
local a2  exp({t2})/`den'
local a3  exp({t3})/`den'
mlexp (ln(normalden(y, `a1'*x1 + `a2'*x2 + `a3'*x3 + {c}, exp({lnsigma}) )))

5. Final remarks

As you probably noticed, the tricky part of this method is finding an appropriate parameterization for the coefficients in your restriction. You may always want to take into account the exponential function to set inequality restrictions, the inverse logit transformation and the hyperbolic tangent transformation to set interval restrictions, and the multinomial logit transformation to set simultaneous interval restrictions, as in the previous section. For example, if you need to set the restrictions:

       a>3
       3<b<10

you can use the transformations:

       a = exp(a1) + 3
       b = 3 + 7*invlogit(b1)

Also, for simplicity, I wrote all the examples by using method d0; however, to achieve convergence, you may need to use methods d1 or d2, in addition to setting good initial values.

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