Announcement

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

  • Error computing Bootstrap Standard Errors (repost)

    I didn't get any responses to my previous question so I will repost it, hopefully someone is kind enough to respond.

    I am encountering an error when trying to compute robust standard errors for the arhomme command. Since arhomme does not support the svyset extension, I am attempting to manually replicate the survey design to obtain the correct standard errors and coefficients. However, when I bootstrap the coefficients, the standard errors and confidence intervals are empty.

    I have attached the code and output below. Any assistance would be appreciated.


    Code:
    
    * Step 1: Save observed coefficients 
    preserve
    quietly xi: arhomme log_avrg_cost i.inc_d i.endentulism i.race i.age_cat i.male i.education i.veteran i.mothered i.wealth i.smoke_now chronicdisease ///
        [pw=new_weight], select(r11dentst = dentalinsurance_w1 endentulism inc_d race age_cat male education veteran mothered wealth smoke_now chronicdisease) ///
        quantiles(0.10, 0.25, 0.50, 0.75, 0.90) taupoints(29) rhopoints(35) meshsize(0.5) frank nostderrors centergrid(-0.20)
    
    matrix beta = e(b)  
    
    * Store sample size
    global N "`e(N)'"
    global Ns "`e(sN)'"
    
    restore
    
    * Step 2: Initialize bootstrap coefficient matrix
    set seed 12345
    global boot_reps = 50
    mat BOOTmatrix = J(1, 119, .)
    
    * Step 3: Run bootstrap
    preserve
    forvalues i = 1/$boot_reps {
        preserve  // Save original dataset before resampling
        
        gsample [w=new_weight], strata(raestrat) cluster(raehsamp)
        
        quietly xi: arhomme log_avrg_cost i.inc_d i.endentulism i.race i.age_cat i.male i.education i.veteran i.mothered i.wealth i.smoke_now chronicdisease ///
            [pw=new_weight], select(r11dentst = dentalinsurance_w1 endentulism inc_d race age_cat male education veteran mothered wealth smoke_now chronicdisease) ///
            quantiles(0.10, 0.25, 0.50, 0.75, 0.90) taupoints(29) rhopoints(35) meshsize(0.5) frank nostderrors centergrid(-0.20)
    
        if `i' == 1 {
            matrix BOOTmatrix = e(b)  // Initialize BOOTmatrix on first iteration
        }
        else {
            matrix BOOTmatrix = (BOOTmatrix \ e(b))  // Append results
        }
    
        restore  // Restore original dataset
    }
    
    restore
    
    * Convert matrix to variables
    svmat BOOTmatrix
    
    * Step 4: Compute bootstrap statistics
    bootstrap r(mean), reps(100) seed(1234) nodrop: summarize BOOTmatrix*, detail
    
    (running summarize on estimation sample)
    
    Bootstrap replications (100): .........10.........20.........30.........40.........50.........60.........70.........80.........90.........100 done
    
    Bootstrap results                                       Number of obs = 12,711
                                                            Replications  =    100
    
          Command: summarize BOOTmatrix*, detail
            _bs_1: r(mean)
    
    ------------------------------------------------------------------------------
                 |   Observed   Bootstrap                         Normal-based
                 | coefficient  std. err.      z    P>|z|     [95% conf. interval]
    -------------+----------------------------------------------------------------
           _bs_1 |  -5.915026          .        .       .            .           .
    ------------------------------------------------------------------------------

  • #2
    I don't think it's about kindness or lack thereof--it's perhaps about a lack of data to work with. And, someone will have to figure out whether what you are doing is legit. Where did you get this idea for the bootstrap to imitate svy?

    I don't think you want to bootstrap the coefficients but the standard errors. svy and reg,pw give the same coefficients.

    I have no idea if this is of any help or even correct, but the results seem plausible.

    Code:
    clear all
    webuse nhanes2f , clear
    svyset psuid [pweight=finalwgt], strata(stratid)
    
    svy: regress zinc age age2 weight female black orace rural 
    estimates store e1
    
    regress zinc age age2 weight female black orace rural [pw=finalwgt] 
    estimates store e2
    
    matrix S = J(500,7,.)
    matrix P = J(500,7,.)
    
    forvalues i = 1/500 {
        
        quietly {
        
        preserve
        
        gsample [w=finalwgt] , strata(stratid)
    
        regress zinc age age2 weight female black orace rural [pw=finalwgt] 
        
        forv k = 1/7 {
            matrix S[`i',`k'] = r(table)[2,`k']
            matrix P[`i',`k'] = r(table)[4,`k']
        }
        
        nois _dots `i' 0
        
        restore    
        } //quietly
    }
    
    * Convert matrix to variables
    svmat S
    svmat P
    
    tabstat S* , save
    tabstat P* , save
    
    esttab e1 e2 , se

    Comment


    • #3
      might replace gsample with bsample , weight(finalwgt) strata(stratid)

      Comment


      • #4
        my apologies, here are he variables I used

        Code:
        * Example generated by -dataex-. For more info, type help dataex
        clear
        input float log_avrg_cost byte r11dentst float(inc_d dentalinsurance_w1 endentulism race age_cat) byte(male education veteran) float mothered byte wealth float(smoke_now chronicdisease new_weight)
         5.860786 1 0 0 0 1 3 0 5 0 1 4 0 3  4918.557
                . 1 0 1 0 1 1 0 5 0 1 4 0 3  6623.596
         6.478509 1 0 1 0 1 2 0 3 0 0 4 0 0  5815.177
         7.824446 1 0 0 0 1 2 0 4 0 1 4 1 0  5893.801
         6.398595 1 0 1 0 1 3 1 5 1 1 4 0 2  6132.698
         8.699764 1 0 1 0 1 2 0 5 0 0 4 0 0   6568.66
         5.525453 1 0 1 0 4 3 0 5 0 1 4 0 2  6938.598
                . 0 0 1 0 3 3 0 1 0 . 2 0 2 2601.5085
         6.685861 1 0 1 0 3 3 0 1 0 0 1 0 2  3354.929
         6.398595 1 0 1 0 1 2 0 1 0 1 4 0 1   7329.36
         7.378384 1 0 0 0 1 4 1 5 1 0 4 0 2  6378.651
                . 1 0 1 0 3 3 0 1 0 0 1 0 2  8065.042
                . 0 0 1 1 4 3 1 1 0 0 1 0 3  5591.879
                . 1 0 1 1 2 3 0 3 0 0 1 0 1  7468.239
                . 0 0 0 0 2 3 1 1 0 . 1 0 0  6987.121
         7.151485 1 0 1 0 1 3 0 4 0 1 1 0 1  4896.986
                . 0 1 1 1 2 3 0 1 0 . 1 0 1  3817.199
         6.216606 1 0 1 0 1 3 0 4 0 1 4 0 0  5999.094
         5.993961 1 0 1 0 1 3 1 5 1 1 4 0 0  6925.122
          8.03948 1 0 0 0 1 3 0 1 0 1 4 0 1  5896.825
                . 0 0 0 0 1 3 0 3 0 1 1 0 2  8667.065
         5.303305 1 0 1 1 4 3 0 5 0 0 4 0 2  6781.349
         5.303305 0 0 1 0 4 3 1 5 1 0 4 0 3  6224.931
         4.620059 1 0 1 0 1 3 0 3 0 0 1 0 0  6747.077
                . 0 0 1 1 1 3 1 4 1 . 2 0 1  6543.964
         7.313887 1 0 1 0 1 2 0 5 0 1 2 0 0  6165.962
         8.881975 1 0 0 0 1 4 0 5 0 1 4 0 2   3695.68
         5.993961 1 0 0 0 1 3 1 5 1 1 4 0 2  6698.451
         7.824446 1 0 0 0 1 3 0 3 0 1 4 0 0  6189.557
         7.523481 1 0 1 0 1 3 1 5 0 0 4 0 3  6313.635
         6.803505 1 0 1 0 1 3 0 5 0 1 4 0 1  7654.785
                . 1 0 1 0 2 3 0 1 0 0 3 0 2 2294.2178
                . 0 0 0 0 2 3 0 5 0 1 1 0 2 2306.8857
                . 0 0 0 1 2 3 0 1 0 0 4 0 1  2917.955
                . 0 0 0 1 2 2 0 3 0 0 1 0 0 2076.2468
         6.398595 1 0 1 0 2 3 1 3 1 1 3 0 3 3628.3115
                . 1 0 1 0 2 3 0 3 0 1 3 0 0 1960.0643
         6.908755 0 0 0 0 2 3 1 3 0 0 2 . 2  6003.062
         5.786897 1 0 1 0 3 1 0 3 0 1 4 0 0  6675.004
         4.836282 1 0 1 0 4 3 0 1 0 0 4 0 0  5975.286
         5.303305 1 0 1 0 4 4 1 5 0 0 4 0 0  6322.484
         6.908755 1 0 1 0 3 3 1 4 1 1 4 0 1 3176.7256
         3.314186 1 0 1 0 1 2 0 4 0 1 4 0 0  3072.397
          7.31422 1 0 0 0 1 2 0 4 0 0 4 0 1  6713.827
         8.987322 1 0 1 0 2 3 0 5 0 0 1 0 4 3192.3496
                . 0 0 1 1 4 3 0 4 0 0 2 0 0  3568.704
         4.912655 1 1 1 0 1 3 1 1 0 1 2 0 1  6398.812
          5.86221 1 0 1 0 1 3 0 5 0 1 2 0 1  6481.468
         5.239098 1 0 1 0 1 3 1 3 1 0 4 0 2  6808.061
                . 0 0 0 0 2 3 0 4 0 0 3 0 3 11825.906
         6.803505 1 0 1 0 1 4 1 5 1 1 4 0 2  7366.496
         7.438972 1 0 0 0 3 3 0 3 0 0 3 0 1   3320.87
                . 0 0 0 1 3 3 1 2 1 0 3 0 1  3323.894
        4.6151204 0 0 0 0 4 4 1 5 0 0 2 0 1 1206.5812
        3.9318256 1 0 0 0 4 3 0 5 0 0 2 0 2 1571.9827
                . 0 0 0 0 2 3 0 1 0 0 1 1 2  5680.608
         5.351858 1 0 0 0 2 3 0 3 0 0 1 0 1  4061.249
                . 0 0 0 0 2 3 0 1 0 . 1 0 2   3636.16
          6.53814 1 0 0 0 1 3 1 4 1 0 1 0 3 3392.4385
         7.313887 1 1 1 0 4 2 0 1 0 1 2 0 6   3709.96
         8.101981 1 0 0 0 1 4 1 3 1 1 4 0 4  5359.752
          6.29803 1 0 0 0 1 3 0 1 0 1 4 0 1  3885.647
         6.055613 1 0 1 0 1 3 0 5 0 1 3 0 1  4188.258
                . 0 0 0 0 1 3 0 3 0 1 2 0 2  3567.327
          7.91972 1 0 0 0 1 4 0 4 0 0 4 0 3  3992.201
         6.908755 1 0 1 0 1 4 1 5 1 . 4 0 1  3890.897
         7.003974 1 0 0 0 1 4 1 4 1 1 4 0 1 4450.5884
         6.763885 1 0 0 0 1 3 0 3 0 0 4 0 4 4880.6187
                . 1 0 1 0 1 2 0 3 0 1 2 0 1  4865.739
                . 1 0 1 0 1 3 0 2 0 0 3 0 0   4232.61
        1.7917595 1 0 0 0 3 3 0 3 0 . 1 0 3 1945.3035
                . 0 0 0 1 2 3 0 4 0 1 2 0 3  2502.309
                . 0 0 0 0 1 3 0 1 0 0 1 0 2  4858.299
         6.685861 1 0 1 0 1 3 1 4 1 1 2 0 1 4108.6255
         5.463832 1 0 0 0 1 3 1 3 0 . 3 0 2  4545.595
                . 0 0 0 0 1 3 0 3 0 1 4 0 1 4109.1294
         4.620059 1 0 1 0 1 4 1 1 1 0 3 0 0 4699.8203
         6.634634 1 0 0 0 1 3 1 5 1 1 4 0 2  2464.117
         5.525453 0 0 1 0 1 3 1 5 0 0 4 0 2  3485.679
         5.860786 1 0 1 0 1 1 0 5 0 1 4 0 0  4411.027
                . 0 0 0 1 1 3 1 2 0 0 3 0 3  4098.926
         6.246107 1 0 0 0 3 2 0 1 0 0 3 . 2  4772.491
                . 0 0 1 0 1 4 1 1 0 0 2 0 3 4599.3877
         5.673323 1 0 0 0 1 4 0 4 0 0 3 0 0  3694.336
                . 0 0 0 1 2 3 0 3 0 0 1 0 3  2498.341
                . 0 0 0 0 2 3 1 2 0 1 1 1 0  5680.104
                . 0 1 0 1 2 3 1 3 0 1 2 1 0  3492.735
                . 0 0 0 0 2 3 0 4 0 0 2 0 2  2664.004
         3.433987 1 0 1 0 2 3 1 4 1 0 3 0 2  1660.687
        3.9318256 1 1 1 0 2 1 0 4 0 0 3 1 2 1695.4633
         5.170484 1 0 0 0 2 3 1 1 0 1 2 0 3  4255.165
                . 0 0 0 1 1 3 1 1 1 0 3 0 4  5537.976
         7.601402 1 0 0 0 1 3 0 1 0 0 4 0 0  5666.496
                . 0 0 0 1 2 3 1 1 1 . 2 0 3   4951.05
                . 0 0 0 0 2 2 0 3 0 1 3 0 1 1885.2877
                . 0 0 0 1 1 2 0 3 0 . 2 1 2 16301.464
                . 0 1 0 1 1 3 1 4 0 1 4 0 4  5449.032
                . 0 0 0 0 1 2 0 3 0 1 4 0 0  5753.079
          5.01728 1 0 1 0 1 3 0 4 0 1 1 0 3 4752.1553
                . 0 0 0 0 1 3 0 3 0 0 3 0 4 4227.8853
        end
        label values r11dentst YESNO
        label def YESNO 0 "0.no", modify
        label def YESNO 1 "1.yes", modify
        label values inc_d inc_d
        label def inc_d 0 "No", modify
        label def inc_d 1 "Yes", modify
        label values dentalinsurance_w1 dentalinsurance_w1
        label def dentalinsurance_w1 0 "No", modify
        label def dentalinsurance_w1 1 "Yes", modify
        label values endentulism endentulism
        label def endentulism 0 "No", modify
        label def endentulism 1 "Yes", modify
        label values race race
        label def race 1 "White", modify
        label def race 2 "Black", modify
        label def race 3 "Hispanic", modify
        label def race 4 "Other", modify
        label values age_cat age_cat
        label def age_cat 1 "50-59", modify
        label def age_cat 2 "60-69", modify
        label def age_cat 3 "70-79", modify
        label def age_cat 4 "80+", modify
        label values male male
        label def male 0 "Female", modify
        label def male 1 "Male", modify
        label values education EDUC
        label def EDUC 1 "1.lt high-school", modify
        label def EDUC 2 "2.ged", modify
        label def EDUC 3 "3.high-school graduate", modify
        label def EDUC 4 "4.some college", modify
        label def EDUC 5 "5.college and above", modify
        label values veteran veteran
        label def veteran 0 "No", modify
        label def veteran 1 "Yes", modify
        label values mothered mothered
        label def mothered 0 "Less than High School", modify
        label def mothered 1 "High School or Higher", modify
        label values smoke_now smoke_now
        label def smoke_now 0 "Non-Smoker", modify
        label def smoke_now 1 "Currently Smokes", modify

        Comment


        • #5
          I have responded here in this thread (https://www.statalist.org/forums/for...e-with-weights) for a version on how to collect the results easily. I attempted to run the actual bootstrap but had various errors as the example dataset is too small to actually bootstrap (singleton clusters). I still think that writing a program, collecting the results in e(b) and using "simulate" to do the main task should work in theory. I have currently no access to Stata but the last time I tried this I ran into another strange problem. The program would run fine as a standalone but would fail in simulate due to a cryptic error. I am not sure if this might be a bug in the arhomme ado. Maybe I can get back to this at the end of the week.
          Best wishes

          (Stata 16.1 MP)

          Comment


          • #6
            I wrote another program but I don't know if I wrote a proper program

            Comment


            • #7
              Originally posted by Felix Bittmann View Post
              I have responded here in this thread (https://www.statalist.org/forums/for...e-with-weights) for a version on how to collect the results easily. I attempted to run the actual bootstrap but had various errors as the example dataset is too small to actually bootstrap (singleton clusters). I still think that writing a program, collecting the results in e(b) and using "simulate" to do the main task should work in theory. I have currently no access to Stata but the last time I tried this I ran into another strange problem. The program would run fine as a standalone but would fail in simulate due to a cryptic error. I am not sure if this might be a bug in the arhomme ado. Maybe I can get back to this at the end of the week.
              Did you get a chance to run the program on your stata? I can add more observations to the dataex, my dataset has around 12,000 observations

              Comment


              • #8
                What I got to work with a simplified command is this. Maybe this can be expanded to the actual question and dataset:

                Code:
                capture program drop arhomme_bootstrap
                program arhomme_bootstrap, eclass
                preserve
                bsample, strata(raestrat)        //Simplified for testing with small sample, maybe use gsample here
                
                *Percentiles simplified to just the median for testing
                quietly xi:  arhomme log_avrg_cost i.inc_d i.endentulism i.race i.age_cat i.male i.education i.veteran i.mothered i.wealth i.smoke_now chronicdisease[pw=new_weight] ///
                    , select(r11dentst = dentalinsurance_w1 endentulism inc_d race age_cat male education veteran mothered wealth smoke_now chronicdisease) ///
                        quantiles(0.50)  taupoints(29) rhopoints(35)  meshsize(0.5) frank nostderrors centergrid(-0.20)
                
                matrix beta_boot = e(b)
                ereturn post beta_boot
                
                restore
                end
                
                
                *** Test program ***
                arhomme_bootstrap        //Works fine
                ereturn list
                
                *** Simulate ***
                simulate _b, seed(123) reps(5): arhomme_bootstrap
                summarize *
                Best wishes

                (Stata 16.1 MP)

                Comment


                • #9
                  You might get a svy dataset from the Stata catalog and set vce(bootstrap) for a simple regression and compare the results from you program to make sure you're getting something comparable.

                  Comment

                  Working...
                  X