1 item has been added to your cart.
| Title | Fitting a regression with interval constraints | |
| Author | Isabel Canette, StataCorp |
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:
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 4.
If you also need to set linear constraints, you can use the option constraint of the command ml.
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.
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.
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 = 74F(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 | Coefficient 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 17
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 | Coefficient 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:
mleval `mu' = `b', eq(1)
args todo b lnf
mleval `mu' = `b', eq(1)
replace `mu' = `b'[1,1]*mpg + `b'[1,2]*displacement + `b'[1,3]
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}))))
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 17
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 | Coefficient 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}))))
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 17
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 | Coefficient 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 | Coefficient Std. err. z P>|z| [95% conf. interval] | |
| _nl_1 | -121.1827 71.29782 -1.70 0.089 -260.9238 18.55855 | |
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}))))
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 17
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 17
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 = 1,000
Wald chi2(0) = .
Log likelihood = -1414.503 Prob > chi2 = .
| y | Coefficient Std. err. z P>|z| [95% conf. interval] | |
| t2 | ||
| _cons | .7524279 .1504271 5.00 0.000 .4575961 1.04726 | |
| t3 | ||
| _cons | .2617819 .1745923 1.50 0.134 -.0804128 .6039765 | |
| c | ||
| _cons | .9913625 .0314883 31.48 0.000 .9296466 1.053078 | |
| lnsigma | ||
| _cons | -.0044355 .0223607 -0.20 0.843 -.0482616 .0393906 | |
| x1 | .2261733 .0259554 .1793569 .2810251 | |
| x2 | .4799727 .0262131 .4289861 .5313799 | |
| x3 | .2938541 .0256475 .2461986 .3464938 | |
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}) )))
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.
