Plot qreg2 + question about comparing confidence intervals in quantile regression with panel data

[Verena Balke]

For my thesis, I am conducting quantile regression for panel data using the qreg2 command in Stata 13.0. The Parente-Santos Silva test shows intra-cluster correlation.
To visualize my results I would like to plot the confidence interval obtained from the quantile regression, but I am also not sure what to plot it against.
I have two questions:

1. To plot the results from qreg2 I have tried to use the grqreg command, but that does not work. Is there a similar command for plotting my results from qreg2?

2. Also, usually the OLS confidence interval is included in the plot to visualize differences between the confidence intervals. Would it not be better to include the confidence interval from a fixed effects regression (I have conducted both LSDV and within estimator). However, in the papers about panel data quantile regression I have read so far OLS confidence intervals are used to make a comparison with the confidence intervals obtained from quantile regression. Could someone explain to me why that is?

[Scott Merryman]

Regarding question one, here is one way:

sysuse auto,clear
tempname memhold
tempfile results
postfile `memhold' q b cilow cihigh using `results'
forv i = .1(.1).9 {
    qui qreg2 price mpg, q(`i')
    local b = _b[mpg]
    local cilow  = _b[mpg]-_se[mpg]*invttail(e(df_r),.025)    
    local cihigh  = _b[mpg]+_se[mpg]*invttail(e(df_r),.025)    
    post `memhold' (`i') (`b') (`cilow') (`cihigh')
}
postclose `memhold'
use `results'
twoway rarea cih cil q, color(gs14)  || line b q,  lc(blue) || , legend(off) xlabel(.1(.1).9)

[Michael Foreman]

Regarding question one, here is one way:

sysuse auto,clear
tempname memhold
tempfile results
postfile `memhold' q b cilow cihigh using `results'
forv i = .1(.1).9 {
qui qreg2 price mpg, q(`i')
local b = _b[mpg]
local cilow = _b[mpg]-_se[mpg]*invttail(e(df_r),.025)
local cihigh = _b[mpg]+_se[mpg]*invttail(e(df_r),.025)
post `memhold' (`i') (`b') (`cilow') (`cihigh')
}
postclose `memhold'
use `results'
twoway rarea cih cil q, color(gs14) || line b q, lc(blue) || , legend(off) xlabel(.1(.1).9)

Dear Scott,

I've been having the same question, so thank you for the solution! When I tested the code, however, the graph has no confidence intervals (please see below). Do you know what the problem might be?
It looks like the confidence intervals are not generated properly

. su

    Variable |        Obs        Mean    Std. Dev.       Min        Max
-------------+---------------------------------------------------------
           q |          9          .5    .2738613         .1         .9
           b |          9   -168.4324    118.5083  -391.5217        -64
       cilow |          0
      cihigh |          0

Thank you once again for your help!

[Joao Santos Silva]

Dear All,

Thanks Scott Merryman for providing the code; I am often asked this question and now I can just point people to your code.
Michael Foreman, I have tested the code in Stata 16 and in Stata 11 and it worked like a charm; maybe you need to update Stata?
Finally, Verena Balke, if you are using panel data, you may want to consider using xtqreg, which allows for fixed effects.

Best wishes,

Joao

[Michael Foreman]

Dear Joao,

Thank you for your reply. My Stata 13.1 is up to date, relatively speaking (please see the code below, it produces the graphs without CIs as I posted earlier). Is there something else I might be doing wrong?

. update

Update status
    Last check for updates:  23 Mar 2020
    New update available:    none         (as of 23 Mar 2020)
    Current update level:    16 Dec 2016  (what's new)

Possible actions

    Check for available updates           (or type -update query-)

. sysuse auto,clear
(1978 Automobile Data)

. 
. tempname memhold

. 
. tempfile results

. 
. postfile `memhold' q b cilow cihigh using `results'

. 
. forv i = .1(.1).9 {
  2. 
.     qui qreg2 price mpg, q(`i')
  3. 
.     local b = _b[mpg]
  4. 
.     local cilow  = _b[mpg]-_se[mpg]*invttail(e(df_r),.025)    
  5. 
.     local cihigh  = _b[mpg]+_se[mpg]*invttail(e(df_r),.025)    
  6. 
.     post `memhold' (`i') (`b') (`cilow') (`cihigh')
  7. 
. }

. postclose `memhold'
. use `results'
. twoway rarea cih cil q, color(gs14)  || line b q,  lc(blue) || , legend(off) xlabel(.1(.1).9)

[Joao Santos Silva]

Dear Michael Foreman,

Does the code generate the cih, cil, and q variables? If so, maybe you can try different plots or different colours.

Best wishes,

Joao

[Michael Foreman]

Dear Joao,

Thank you for your advice. It looks like the code does not generate confidence intervals (please see below). I am really perplexed because I just copied and pasted exactly the same code. It has something to do with " *invttail(e(df_r),.025) " part. if I replace "e(df_r)" with a value, say "73", the code generates confidence intervals, but if i keep "e(df_r)" as in the original code, i get missing values. Would you have any idea why this might be the case with my Stata? Are there some setting that could have been manually and now cause this issue? Apprently, "e(df_r)" does not work in my case. I am really perplexed by why it works on other machines but not on mine.

. sysuse auto,clear
(1978 Automobile Data)

. tempname memhold

. tempfile results

. postfile `memhold' q b cilow cihigh using `results'

. forv i = .1(.1).9 {
  2.     qui qreg2 price mpg, q(`i')
  3.     local b = _b[mpg]
  4.     local cilow  = _b[mpg]-_se[mpg]*invttail(e(df_r),.025)    
  5.     local cihigh  = _b[mpg]+_se[mpg]*invttail(e(df_r),.025)    
  6.     post `memhold' (`i') (`b') (`cilow') (`cihigh')
  7. }

. postclose `memhold'

. use `results'

. list

     +---------------------------------+
     |  q           b   cilow   cihigh |
     |---------------------------------|
  1. | .1   -75.71429       .        . |
  2. | .2         -64       .        . |
  3. | .3         -84       .        . |
  4. | .4        -107       .        . |
  5. | .5   -135.6667       .        . |
     |---------------------------------|
  6. | .6   -120.6667       .        . |
  7. | .7   -202.5833       .        . |
  8. | .8   -334.7391       .        . |
  9. | .9   -391.5217       .        . |
     +---------------------------------+

[Joao Santos Silva]

Dear Michael Foreman,

I should have asked this before: do you have the latest version of qreg2?

Best wishes,

Joao

[Michael Foreman]

Dear Joao,

Following your advice, I reinstalled qreg2 and now I am getting the confidence intervals! It is a mystery to me what the problem was.. Thank you for your help!

[Joao Santos Silva]

Dear Michael Foreman,

An earlier version of qreg2 did not save the degrees of freedom.

Best wishes,

Joao

[Michael Foreman]

Thank you, Joao, this clarifies it all!

[Suna Marasli]

Regarding question one, here is one way:

sysuse auto,clear
tempname memhold
tempfile results
postfile `memhold' q b cilow cihigh using `results'
forv i = .1(.1).9 {
qui qreg2 price mpg, q(`i')
local b = _b[mpg]
local cilow = _b[mpg]-_se[mpg]*invttail(e(df_r),.025)
local cihigh = _b[mpg]+_se[mpg]*invttail(e(df_r),.025)
post `memhold' (`i') (`b') (`cilow') (`cihigh')
}
postclose `memhold'
use `results'
twoway rarea cih cil q, color(gs14) || line b q, lc(blue) || , legend(off) xlabel(.1(.1).9)

Dear all,
I tried applying this suggestion to visualize my estimation, like below:

tempname memhold
tempfile results
postfile `memhold' q b cilow cihigh using `results'
forv i = .1(.1).9 {
    xi: qreg2 lavprod zip $general i.country i.yy logdeg_t logsz, cluster(new_account_id) quantile(.75) wlsiter(10) q(`i')
    local b = _b[zip]
    local cilow  = _b[zip]-_se[zip]*invttail(e(df_r),.025)    
    local cihigh  = _b[zip]+_se[zip]*invttail(e(df_r),.025)    
    post `memhold' (`i') (`b') (`cilow') (`cihigh')
}
postclose `memhold'
use `results'
twoway rarea cih cil q, color(gs14)  || line b q,  lc(blue) || , legend(off) xlabel(.1(.1).9)

But I get an error about the q option, as "option q() not allowed". I thought it may be related to using dummy variables and "xi", but I couldn't solve the problem.

I am using Stata 16 on a Windows pc, and installed qreg2 just yesterday, I guess it is up to date.

Best regards,
Suna

[Scott Merryman]

You have specified the quantile(#) twice:

...quantile(.75) wlsiter(10) q(`i')

[Suna Marasli]

Thank you Scott Merryman! Such a beginners' mistake.

Yet, I still cannot make it work. My code is as below:

tempname memhold
tempfile results
postfile `memhold' q b cilow cihigh using `results'
forv i = .1(.1).9 {
    qreg2 lavprod zip $general i.country i.yy logdeg_t logsz, cluster(new_account_id) wlsiter(10) q(.50)
    local b = _b[zip]
    local cilow  = _b[zip]-_se[zip]*invttail(e(df_r),.025)  
    local cihigh  = _b[zip]+_se[zip]*invttail(e(df_r),.025)  
    post `memhold' (0.50) (`b') (`cilow') (`cihigh')
}
postclose `memhold'
use `results'
twoway rarea cih cil q, color(gs14)  || line b q,  lc(blue) || , legend(off) xlabel(.1(.1).9)

After the line " use `results' " I get an error like:
"no; dataset in memory has changed since last saved"

I am not experienced with coding graphics in Stata. I cannot make sense of this error. Help greatly appreciated.

Best regards,
Suna

[FernandoRios]

Hi Suna,
perhaps the following code may help with what yo need to do

 clear all
 use http://fmwww.bc.edu/RePEc/bocode/o/oaxaca.dta, clear
 
 forvalues i=.05(.05).95 {
  qui:qreg lnwage educ exper tenure female, q(`i')
  matrix aux=r(table)
  matrix bqreg=nullmat(bqreg)\aux["b",....]
  matrix llqreg=nullmat(llqreg)\aux["ll",....]
  matrix uuqreg=nullmat(uuqreg)\aux["ul",....]
  qui:qreg2 lnwage educ exper tenure female, q(`i')
  matrix aux=r(table)
  matrix bqreg2=nullmat(bqreg2)\aux["b",....]
  matrix llqreg2=nullmat(llqreg2)\aux["ll",....]
  matrix uuqreg2=nullmat(uuqreg2)\aux["ul",....]
 }
 ssc install lbsvmat
 lbsvmat bqreg
 lbsvmat llqreg
 lbsvmat uuqreg

 lbsvmat llqreg2
 lbsvmat uuqreg2
 
 gen q=_n*5 if  _n*5<=95
 
 two rarea llqreg1 uuqreg1 q, color(green%50)  ||  ///
    rarea llqreg21 uuqreg21 q, color(blue%50) || line bqreg1 q, ///
    legend(order( 1 "qreg CI" 2 "qreg2 CI" 3 "Point estimate"))

 two rarea llqreg2 uuqreg2 q, color(green%50)  ||  ///
    rarea llqreg22 uuqreg22 q, color(blue%50) || line bqreg2 q, ///
    legend(order( 1 "qreg CI" 2 "qreg2 CI" 3 "Point estimate"))

HTH