Announcement

Collapse
No announcement yet.
X
  • Filter
  • Time
  • Show
Clear All
new posts

  • Reporting comparable estimates and boottest result for OLS and Poisson regression

    Dear community,


    I would like to populate the following table
    naloxone drug rehab
    t-stat P value t-stat P value
    Linear regression
    Beta
    CV1
    WBS
    CV3
    Poisson regression
    Beta
    CV1
    WBS

    The table reports the points estimates (Beta) computed with OLS or Poisson regression. Both point estimates should be reported as average marginal effects. The table also reports cluster robust errors (CV1), wild bootstrap clustered errors (WBS) and clustered errors jackknife errors (CV3).


    Consider this hypothetical data.

    Code:
    * Example generated by -dataex-. For more info, type help dataex
    clear
    input byte(Drugrehab Naloxone) float HN long st byte(_est_est1 _est_est2 _est_est3 _est_est4)
     68 96 0 12 1 1 1 1
     59 47 0  9 1 1 1 1
     59 69 1 43 1 1 1 1
     43 21 1 27 1 1 1 1
     83 32 0 43 1 1 1 1
     41 25 0 42 1 1 1 1
     51  0 0 17 1 1 1 1
     43 25 0 20 1 1 1 1
     57 15 0  2 1 1 1 1
     45 40 0 29 1 1 1 1
     65 74 1 10 1 1 1 1
     58  0 0 34 1 1 1 1
     31  0 0 47 1 1 1 1
     47  0 0 10 1 1 1 1
    100  0 0 16 1 1 1 1
     39 80 0 21 1 1 1 1
     48 60 0 35 1 1 1 1
     16 57 1 35 1 1 1 1
     78  0 0 46 1 1 1 1
     40  0 0 28 1 1 1 1
    end
    label values st st
    label def st 2 "AL", modify
    label def st 9 "FL", modify
    label def st 10 "GA", modify
    label def st 12 "IA", modify
    label def st 16 "KS", modify
    label def st 17 "KY", modify
    label def st 20 "MD", modify
    label def st 21 "ME", modify
    label def st 27 "NC", modify
    label def st 28 "ND", modify
    label def st 29 "NE", modify
    label def st 34 "NY", modify
    label def st 35 "OH", modify
    label def st 42 "TX", modify
    label def st 43 "VA", modify
    label def st 46 "WI", modify
    label def st 47 "WV", modify

    This is my code

    **************************
    eststo clear
    **************************
    cls
    foreach g in ///
    Naloxone ///
    Drugrehab ///
    {
    foreach h in ///
    gaussian ///
    poisson ///
    {
    eststo: quietly glm ///
    `g' ///
    HN ///
    , ///
    family(`h') ///
    vce(cl st ) //
    estimates store results
    quietly boottest HN, nograph
    estadd scalar WBS_p = r(p)
    estadd scalar WBS_z = r(z)
    estimates restore results
    eststo: margins, dydx(HN) post
    }
    **************************
    }
    **************************
    esttab, ///
    cells(b p t) ///
    star(+ 0.10 * 0.05 ** 0.01 *** 0.001) ///
    stats(WBS_p WBS_z) ///
    keep(HN) //
    **************************

    This is the output. The problem is that scalars WBS_p and WBS_z are not reported.

    Naloxone Naloxone Drugrehab Drugrehab
    b/p/t b/p/t b/p/t b/p/t b/p/t b/p/t b/p/t b/p/t
    main
    HN 16.94118 16.94118 0.4037365 14.13078 5.631579 5.631579 0.0985611 5.386363
    0.0033685 0.0033685 0.0173224 0.0001545 0.2831716 0.2831716 0.3071589 0.2601239
    2.93194 2.93194 2.379793 3.783774 1.073222 1.073222 1.021202 1.126098
    WBS_p
    WBS_z

    If I do not store results or use margins, as here

    eststo clear
    **************************
    cls
    foreach g in ///
    Naloxone ///
    Drugrehab ///
    {
    foreach h in ///
    gaussian ///
    poisson ///
    {
    eststo: quietly glm ///
    `g' ///
    HN ///
    , ///
    family(`h') ///
    vce(cl st ) //
    // estimates store results
    quietly boottest HN, nograph
    estadd scalar WBS_p = r(p)
    estadd scalar WBS_z = r(z)
    // estimates restore results
    // eststo: margins, dydx(HN) post
    }
    **************************
    }
    **************************
    esttab, ///
    cells( b p t) ///
    star(+ 0.10 * 0.05 ** 0.01 *** 0.001) ///
    stats(WBS_p WBS_z) ///
    keep(HN) //

    I have bootstrapped p and z (t for OLS) values, but these estimates are not directly comparable since Poisson regression follows a nonlinear routine.
    Naloxone Naloxone Drugrehab Drugrehab
    b/p/t b/p/t b/p/t b/p/t
    main
    HN 16.94118 0.4037365 5.631579 0.0985611
    0.0033685 0.0173224 0.2831716 0.3071589
    2.93194 2.379793 1.073222 1.021202
    WBS_p 0.4574575 0.4354354 0.4984985 0.5005005
    WBS_z 1.027171 1.027171 0.7620507 0.7620507

    I suspect that `estimates store' rewrite something. Please help.

    PS. I would like to compute CV3 for Poisson too, but doing it with vce(jackknife, mse cluster(st)) would take too long (actual data is much larger), and boottest does not support jk and score at the same time. Is this correct? If you could advise on CV3 for Poisson, too, I'd be grateful.

    C1 and CV3 follow the language of
    https://www.sciencedirect.com/scienc...04407622000781


    Kind regards,
    Sergey Alexeev | ​Senior Research Fellow - Health Economist​ ​
    NHMRC Clinical Trials Centre
    The University of Sydney, Faculty of Medicine and Health
    Kind regards,
    Sergey Alexeev | ​The University of Sydney
    https://alexeev.pw/

  • #2
    Do you have the same problem if you use test instead of boottest? If so then this, is a question about proper use of eststo and esttab.

    Comment


    • #3
      Dear David,

      Thank you for getting back to me. My apologies for an absolutely terrible explanation of my issue in the original post.

      My issue persists. Here is a clear illustration of the problem. I am running this code.


      Code:
      **************************
      sysuse auto.dta, clear
      **************************
      eststo clear
      glm ///
      mpg ///
      weight ///
      foreign ///
      , ///
      family(poisson) ///
      irls //
      **************************
      estimates store results
      boottest weight, nograph
      estadd scalar WBS_p = r(p)
      estadd scalar WBS_z = r(z)
      estimates restore results
      eststo: margins, dydx(weight) post
      **************************
      esttab , ///
      cells( b(star fmt(3)) se( fmt(3) par) p( fmt(3) par) t( fmt(3) par) ) ///
      star(+ 0.10 * 0.05 ** 0.01 *** 0.001) ///
      stats(N WBS_p WBS_z) ///
      keep(weight) //
      **************************

      I get this table.

      Code:
      ----------------------------
                            (1)   
                                  
                       b/se/p/t   
      ----------------------------
      weight             -0.007***
                        (0.001)   
                        (0.000)   
                       (-7.626)   
      ----------------------------
      N                      74   
      WBS_p                       
      WBS_z                       
      ----------------------------
      With WBS_p and WBS_z emply, but would like to have results from boottest there. I spend so much time trying to figure this out and nothing helps.
      Kind regards,
      Sergey Alexeev | ​The University of Sydney
      https://alexeev.pw/

      Comment


      • #4
        After the boottest command, does "display r(p), r(z)" show the results? If so, then I don't think this is a question about boottest, but maybe about estadd.

        Comment


        • #5
          boottest and estout are from SSC, as you are asked to explain in FAQ Advice #12. This is a bit of a mess. With the -post- option of margins, you overwrite the previously stored results. So you need to hold the boottest statistics in locals or scalars if you want to include them with the margins output.

          Code:
          sysuse auto.dta, clear
          **************************
          eststo clear
          glm ///
          mpg ///
          weight ///
          foreign ///
          , ///
          family(poisson) ///
          irls //
          **************************
          estimates store results
          boottest weight, nograph
          local WBS_p = r(p)
          local WBS_z = r(z)
          estimates restore results
          eststo m1: margins, dydx(weight) post
          estadd local WBS_p= cond(`WBS_p'<0.001, "<0.001", "`:di %4.3f `WBS_p''"): m1
          estadd scalar WBS_z = `WBS_z': m1
          **************************
          esttab m1, ///
          cells( b(star fmt(3)) se( fmt(3) par) p( fmt(3) par) t( fmt(3) par) ) ///
          star(+ 0.10 * 0.05 ** 0.01 *** 0.001) ///
          stats(N WBS_p WBS_z) ///
          keep(weight) //
          ******
          Res.:

          Code:
          . esttab m1, ///
          > cells( b(star fmt(3)) se( fmt(3) par) p( fmt(3) par) t( fmt(3) par) ) ///
          > star(+ 0.10 * 0.05 ** 0.01 *** 0.001) ///
          > stats(N WBS_p WBS_z) ///
          > keep(weight) //
          
          ----------------------------
                                (1)  
                                      
                           b/se/p/t  
          ----------------------------
          weight             -0.007***
                            (0.001)  
                            (0.000)  
                           (-7.626)  
          ----------------------------
          N                  74.000  
          WBS_p              <0.001  
          WBS_z             -10.556  
          ----------------------------
          Last edited by Andrew Musau; 18 Dec 2023, 01:43.

          Comment


          • #6
            Dear David Roodman,

            Thank you very much again for your time. This has always been an estout issue; I named this topic incorrectly. I'm sorry. Your work on boottest is outstanding. Thank you for your contributions and care.

            Dear Andrew Musau,

            Your responses are so epically good, as usual. Everything works as expected. I should somehow buy you a cup of coffee (a bottle of whisky?) for all you work here. You helped me so much throughout the years.

            Warm regards, and have a great holiday season!!!
            Sergey Alexeev
            Kind regards,
            Sergey Alexeev | ​The University of Sydney
            https://alexeev.pw/

            Comment


            • #7
              Originally posted by Sergey Alexeev View Post
              I should somehow buy you a cup of coffee (a bottle of whisky?)
              Thanks for the offer. If I am rational at all, I should take the bottle of whisky. For one, I could sell it and buy several cups of coffee. But, I am not known to be a coffee drinker, I prefer tea like the Brits!

              Comment

              Working...
              X