Plotting confidence intervals around quadratic fitted regression lines

[Dany Shakeel]

Hi,

I am running the code below and it generates and plots four lines. Instead, I want to plot bands of the confidence interval for each line (l4, l8, l12 and lboth). How to do that?

reg smean byr byrsq test1 test2 test3 test4 if ES & math
gen l4 =_b[_cons]+_b[byr]*byr + _b[byrsq]*byrsq if ES & math

reg smean byr byrsq test1 test2 test3 test4 if MS & math
gen l8 =_b[_cons]+_b[byr]*byr + _b[byrsq]*byrsq if MS & math

reg smean byr byrsq test1 test2 test3 test4 if HS & math
gen l12 =_b[_cons]+_b[byr]*byr + _b[byrsq]*byrsq if HS & math

reg smean byr byrsq test1 test2 test3 test4 ES MS HS if math
gen lboth =_b[_cons]+_b[byr]*byr + _b[byrsq]*byrsq

sort birthyr
twoway (scatter l4 birthyr if ES, msize(vsmall) connect(l) mcolor(red) lcolor(red)) ///
(scatter l8 birthyr if MS, msize(vsmall) connect(l) mcolor(navy) lcolor(navy)) ///
(scatter l12 birthyr if HS, msize(vsmall)connect(l) mcolor(dkgreen) lcolor(dkgreen)) ///
(scatter lboth birthyr, msize(vsmall) connect(l) mcolor(maroon) lcolor(maroon) ///
title("Levels for four categories", size(small)) xtitle("Year", size(small)) legend(off) xlabel(1900(10)2000, labsize(small)) ///
ytitle("Mean value", size(small)) ylab(-.8(.1).5, labsize(small) angle(horizontal)) ///
name(fig5, replace) yscale(titlegap(*-20)) ///
text(.1 1980 "first case", size(2) color(red)) ///
text(.13 1986 "second case", size(2) color(navy)) ///
text(-.02 1990 "third case", size(2) color(dkgreen)) ///
text(.13 1997 "Combined", size(2) color(maroon)))

[Bruce Weaver]

Dany, I suspect you'll have an easier time getting the plots you want via -margins- and -marginsplot-. See the examples on this UCLA webpage, for example. In order for -margins- to work properly, you'll need to use factor variable notation when specifying your regression models. In other words, replace byr byrsq with c.byr##c.byr. For more info about factor variables, type -help fvvarlist-. HTH.

[Dany Shakeel]

Bruce, I tried your suggestion for the first line l4 (code below). However, it gives just one data point. I know there is a qfit command for plotting quadratic fitted regression lines. The problem is that qfit does not allow controls test1-test4. Any suggestions?

reg smean byr c.byr##c.byr test1-test4 if ES & math
margins, dydx(byr)
marginsplot

[Bruce Weaver]

See if this example helps you get started.

sysuse auto, clear
summarize weight
local wtmin = r(min)
local wtmax = r(max)

regress mpg c.weight##c.weight foreign length turn
quietly margins, at(weight=(`wtmin'(10)`wtmax')) vsquish
marginsplot, recast(line) recastci(rarea) ciopt(color(*.3))

[Dany Shakeel]

Thanks Bruce. When I use the margins command and compare the start and end points of the graph with the ones I had generated for quadratic fit

(gen l4 =_b[_cons]+_b[byr]*byr + _b[byrsq]*byrsq if ES & math)

there seems to be discrepancies between them. It seems that margins command does not generate exactly same values as the quadractic fit.

[Dany Shakeel]

I successfully ran the following command to fit variable smean on a quadratic fit of variable byr with controls ES-HS

sum byr if !missing(smean), meanonly
local byrmin=r(min)
local byrmax=r(max)
reg smean c.byr##c.byr ES-HS [w=ORIGWT]
quietly margins if !missing(smean), at(byr=(`byrmin'(1)`byrmax'))
marginsplot, ylab(-1(.5)1) recast(line) recastci(rarea) ciopt(color(*.3))

Now I want to calculate if the difference between the margins of byrmin and byrmax is statistically significant. This would allow me to test if the change in start and end points of the plotted line is statically significant.

Can anyone help with this calculation?

[Bruce Weaver]

Hi Dany. If you generate an example of what you are trying to do with one of the datasets that comes with Stata (e.g., the auto dataset I used in #4), people will be much more inclined to try to help. Please also enclose and code with code delimiters. To do so, highlight the code (or output) and then click on the # symbol in the toolbar.

Meanwhile, I wonder if -lincom- will give you what you are looking for in #6.

[Dany Shakeel]

Hi Bruce. Thanks for your response. and I will use your suggestions for posting in future. In the example in #4, I want to calculate if the margin at wtmin is statistically different from the margin at wtmax. If you look at the graph generated by marginsplot in #4, the begin and end points of the line have confidence intervals. I want to test if the change from beginning of the line to its end is statically significant. Thus a difference in magnitudes of margins with a p-value is all what I need.

I appreciate your help.

[Bruce Weaver]

Here's an example showing what I had in mind when I suggested using -lincom-. Notice that I had to use -nlcom- rather than -lincom- to get it to work.

sysuse auto, clear

regress mpg c.weight##c.weight foreign length turn
*regress, coefleg noheader

quietly summarize weight if e(sample)
local wtmin = r(min)
local wtmax = r(max)

* Let yhat1 = yhat when weight = min(weight)
* Let yhat2 = yhat when weight = max(weight)

local yhat1  (_b[weight]*`wtmin' + _b[c.weight#c.weight]*`wtmin'^2)
local yhat2  (_b[weight]*`wtmax' + _b[c.weight#c.weight]*`wtmax'^2)

* Plot yhat as a function of weight
quietly margins, at(weight=(`wtmin'(10)`wtmax')) vsquish
marginsplot, recast(line) recastci(rarea) ciopt(color(*.3)) ///
xlab(2000(500)5000, grid) ylab(5(5)35, grid)
* Display fitted values for weight = min and weight = max
margins, at(weight=(`wtmin' `wtmax')) vsquish
display 16.266-29.413

* I had in mind using -lincom- as follows:
* lincom `yhat2' - `yhat1'
* But this generates an error:  "not possible with test"
* Therefore, I then tried -nlcom- in place of -lincom-.
nlcom `yhat2' - `yhat1'
* As you can see, -nlcom- works.  But it yields a z-test rather than a t-test.

[Dany Shakeel]

Thank you very much Bruce. Much appreciated!

[Dany Shakeel]

Hi Bruce,

Could you help with two more issues related to above regression:

1) I need to run multiple regressions and producing an output in CSV or Excel would immensely help.

How to output the the following in Excel or CSV:

a) # of observations from the regression,
b) the value at the start of margin i.e. 'yhat1' and its standard error,
c) the value at the end of margin i.e. 'yhat2' and its standard error and
d) the value of the difference from nlcom `yhat2' - `yhat1' and its standard error

Having asterisks for p-values for b), c and d) would be helpful

2) I am using combomarginsplot command to plot 4 or 5 margins from such regressions. I'd like the resulting STATA graph to display the values of 'yhat2' at the end and 'yhat1' at the start for each margin line. The displayed values in the graph can be rounded to two decimal places for ease. Is there a way to do that?

Thanks.