It is not uncommon to believe a variable x predicts a variable y
differently over certain ranges of x. In such instances, you may
wish to fit a piecewise regression model. The simplest scenario
would be fitting two adjoined lines: one line defines the relationship of y
and x for x <= c and the other line defines the
relationship for x > c. For this scenario, we can use the
Stata command
nl to find the value of c that yields the best fitting model.
The nl command in Stata performs nonlinear least-squares estimation and
allows the user to define the function for which it estimates indicated
parameters. It is extremely flexible and very useful, though slightly tricky to
use at first. The page presents a rather simple example. For details
and more examples on nl, see its Stata help page.
For this example, we will use a fictional dataset where the relationship
between x and y is clearly not a single line.
use https://stats.oarc.ucla.edu/wp-content/uploads/2022/07/nl.dta, clear twoway scatter y x
We might look at this plot and
believe that there is a downward trend in y as x increases up to a
certain point in x. After that point, there is an upward trend in
y. Let’s consider the set of parameters we will need to fit.
Our first line will involve a slope and an intercept (a1 and b1);
our second line will also involve a slope (b2) and we can think of the
point at which it meets the first line as its “intercept” defined by the first
intercept, the first slope, and the point at which the lines meet (c). We want to
estimate four total parameters: two slopes, an intercept, and a cut point.
We can indicate these parameters in our nl command and provide starting points for
each parameter based on the plot above.
nl (y = ({a1} + {b1}*x)*(x < {c}) + ///
({a1} + {b1}*{c} + {b2}*(x-{c}))*(x >= {c})), ///
initial(a1 25 b1 -2 c 10 b2 2)
Source | SS df MS
-------------+------------------------------ Number of obs = 200
Model | 8770.59791 3 2923.53264 R-squared = 0.5169
Residual | 8197.31882 196 41.8230552 Adj R-squared = 0.5095
-------------+------------------------------ Root MSE = 6.467075
Total | 16967.9167 199 85.2659132 Res. dev. = 1310.224
------------------------------------------------------------------------------
y | Coef. Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
/a1 | 18.53111 1.382652 13.40 0.000 15.80433 21.2579
/b1 | -1.920463 .2668338 -7.20 0.000 -2.446697 -1.394229
/c | 8.987615 .4400011 20.43 0.000 8.11987 9.855359
/b2 | 2.267615 .1915718 11.84 0.000 1.889808 2.645422
------------------------------------------------------------------------------
Parameter a1 taken as constant term in model & ANOVA table
From the output above, we can see estimates of all four
parameters. We can use the estimate for the cut point c to generate
a new variable, x2, that will allow us to run an ordinary least squares
regression of y on x and x2 that effectively fits a
piecewise function.
gen x2 = x - 8.987615
replace x2 = 0 if x < 8.987615
regress y x x2
Source | SS df MS Number of obs = 200
-------------+------------------------------ F( 2, 197) = 105.39
Model | 8770.59777 2 4385.29889 Prob > F = 0.0000
Residual | 8197.31896 197 41.6107562 R-squared = 0.5169
-------------+------------------------------ Adj R-squared = 0.5120
Total | 16967.9167 199 85.2659132 Root MSE = 6.4506
------------------------------------------------------------------------------
y | Coef. Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
x | -1.920463 .2023035 -9.49 0.000 -2.319422 -1.521505
x2 | 4.188078 .3214448 13.03 0.000 3.554163 4.821993
_cons | 18.53111 1.275748 14.53 0.000 16.01524 21.04699
------------------------------------------------------------------------------
In the regression output, we can see that we have the same sum of
squares we saw in the nl output. We also see that our intercept is
unchanged (a1 in the nl output, _cons in the regress
output), the coefficient for x matches the first slope from nl, and the coefficient for x2 is equal to (b2 – b1).
We can plot the predicted values from the regression above.
predict p graph twoway (scatter y x) (scatter p x)
We have found the optimal point to split our piecewise function in this
scenario. The same process could be used if we wished to fit quadratic or
cubic terms, as long as we carefully described the function and its parameters
in our nl command.