Announcement

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

  • Bootstrap loop.

    Dear All,

    I am estimating a simple predictive model for the number of opioid related deaths. I want to check how predictions improve as the sample size increases. For that, I am trying to bootstrap the prediction, within a loop, that estimates the predictions as sample size increases from n=500 to n=2500.

    Here is a dataex of count 100 of the sample:

    Code:
    * Example generated by -dataex-. To install: ssc install dataex
    clear
    input double(ICD_Opioids Sex T509) byte(Black age030 Home)
    0 1 1 0 0 1
    0 1 1 0 0 1
    1 1 0 0 0 1
    0 1 1 0 0 1
    0 1 1 0 0 0
    0 1 1 0 0 0
    0 1 1 0 0 1
    0 1 1 0 0 0
    1 1 0 0 0 1
    1 1 0 0 1 1
    1 1 0 0 0 1
    0 1 1 1 0 1
    0 1 1 0 0 1
    0 0 1 0 1 0
    0 1 1 0 0 1
    1 1 0 0 1 0
    0 0 1 0 0 1
    0 0 1 0 1 0
    0 1 1 0 0 0
    0 1 1 0 0 0
    0 0 1 0 0 1
    0 0 1 1 0 1
    0 1 1 0 0 0
    1 1 0 0 0 1
    1 0 0 0 0 0
    0 0 1 0 0 0
    0 1 1 0 0 1
    0 0 1 0 0 1
    1 1 0 0 1 0
    1 0 0 0 0 1
    1 0 0 0 0 0
    1 1 0 0 1 0
    0 0 0 1 0 0
    0 1 1 0 0 1
    1 1 0 0 0 1
    0 1 1 0 0 0
    1 0 0 0 0 1
    1 1 0 0 0 0
    0 0 1 0 0 1
    1 1 0 0 0 1
    0 0 1 0 1 0
    0 0 1 0 0 1
    1 1 0 0 1 0
    0 1 1 0 0 1
    0 0 0 0 0 1
    1 0 0 0 1 1
    0 1 1 1 0 1
    0 1 1 0 0 0
    1 0 0 0 1 0
    0 1 1 0 1 0
    1 1 0 0 1 1
    0 1 1 0 0 1
    1 1 0 0 0 0
    1 0 0 1 0 0
    1 1 0 0 1 1
    0 1 1 0 1 1
    0 1 0 1 0 0
    0 1 1 0 0 0
    1 0 0 0 0 1
    0 1 1 0 0 1
    0 1 1 0 0 1
    1 1 0 0 0 0
    0 1 1 0 1 0
    1 1 0 0 0 1
    0 0 1 0 0 1
    0 1 1 0 0 1
    0 1 1 0 0 1
    1 1 0 0 1 0
    0 0 1 0 1 1
    0 1 1 0 1 0
    0 1 1 0 0 1
    0 1 1 0 0 1
    1 1 0 0 1 1
    0 0 1 0 0 1
    0 1 1 0 0 1
    0 1 1 0 1 0
    0 1 0 0 0 1
    0 1 1 0 1 1
    0 0 1 0 0 1
    1 0 0 0 0 1
    0 0 1 0 0 1
    1 0 0 0 1 0
    0 1 1 0 0 0
    0 1 1 0 0 1
    1 1 0 0 0 0
    0 1 1 0 1 1
    1 1 0 0 0 1
    0 1 1 0 1 1
    1 1 0 0 0 1
    1 0 0 0 0 1
    1 0 0 0 1 0
    1 0 0 0 0 0
    1 1 0 1 0 1
    0 1 1 0 0 0
    0 0 1 0 0 1
    1 1 0 1 0 1
    0 1 1 0 0 0
    1 1 0 0 1 1
    0 0 1 0 0 1
    0 1 1 0 1 0
    end
    label values ICD_Opioids ICD_Opioids
    label def ICD_Opioids 0 "Absent", modify
    label def ICD_Opioids 1 "Present", modify
    label values Sex Sex
    label def Sex 0 "Female", modify
    label def Sex 1 "Male", modify
    label values T509 T509
    label def T509 0 "Absent", modify
    label def T509 1 "Present", modify
    And the loop I have written is as follows:

    g Male = 1 if Sex == 1
    replace Male = 0 if Sex == 0
    g MalexBlack = Male*Black


    forvalues i=500/2500 {
    capture program drop myboot
    program define myboot, rclass
    preserve
    bsample

    probit ICD_Opioids Sex Black Home age030 if (T509~=1)

    predict death_hat if T509==1, pr
    gen deathpredicted=1 if death_hat>=.5
    recode deathpredicted .=0
    replace deathpredicted = . if T509 !=1

    egen numdeaths_ruhm_tmp = sum(deathpredicted) // Calculate estimated number of opioid deaths with Ruhm approach
    sum numdeaths_ruhm_tmp // Summarize this--will extract mean in next line
    return scalar numdeaths_ruhm = r(mean)

    egen numdeaths_tox_tmp = sum(T509_TOXOPIOID) // Calculate estimated number of opioid deaths from toxicology data
    sum numdeaths_tox_tmp // Summarize this--will extract mean in next line
    return scalar numdeaths_tox = r(mean)

    return scalar numdeaths_diff = numdeaths_ruhm - numdeaths_tox // Main statistic of interest. Is it different from 0?

    restore
    end

    bootstrap numdeaths_diff = r(numdeaths_diff) numdeaths_ruhm = r(numdeaths_ruhm) numdeaths_tox = r(numdeaths_tox), size(`i') /*saving(bootstrap10, replace)*/ reps(500) seed(1234): myboot
    save "bootstrap.dta", append
    }

    It doesn't run at all. I am a somewhat experienced user of Stata but with no experience with loops. So I suspect this has multiple problems but I don't even know where to begin since it does not run at all. I will be grateful for all/ any help please.
    Sincerely,
    Sumedha.

  • #2
    Well, there are several clear problems with this code. Given the limited data example, I can't run my corrected code to test it, but I think it will work with a real data set.

    First, the block of code from -capture program drop myboot- through -end-, which defines the program myboot, should not be inside the loop. It should appear before the loop.

    Second, the code in that program refers to a non-existent variable T509_TOXOPIOID. If your real data set has such a variable, then no change is needed. If it does not, you need to replace this with a reference to an actual variable in your data. Perhaps you meant just T509?

    The -save- command does not have an -append- option. There is a separate -append- command for that. See -help append-. What you need to do, to accumulate the results of each bootstrap, is to start with an empty results file saved, and then inside your loop append the latest results to it, and re-save the results file. Then when you are all done, the results file will have everything in one place, and you can -use- it to do further work with those results.

    Since your -probit- inside -myboot- is resricted to observations where T509 = 1, you should, before entering your loop, drop all the observations where T509 != 1. The reason is that if you fail to do that, your bootstrap samples will include a random mixture of observations with T509 == 1 and T509 != 1. Since the latter will be excluded from the -probit-, your actual sample size will not be the value of `i' from your loop, but will be some random number less than or equal to that and will vary in each of your bootstrap samples. -bootstrap- is likely to get confused about that and send you disturbing messages. And some of your samples may even contain so few T509 == 1 observations that you can't even do the -probit- estimation.

    Even if you get this running correctly, looping -forvalues i = 500/2500- will do this once for every integer i between 500 and 2500 inclusive. That's 2000 bootstraps, each of which entails 500 replications, for a total of 1,000,000 probit regressions. Do you really care about the difference between sample size 500 and sample size 501? If you do, then proceed, but prepare to have your children, or perhaps your grandchildren, analyze the results. I would suggest something more like -forvalues i = 500(100)2500-, which will look at 500, 600, 700,..., 2400, 2500, or perhaps even fewer values of i.

    Since you are looping, do not use the -seed()- option of the bootstrap command: by resetting the same seed with each bootstrap you are not getting independent results. So use a -set seed- command near the top of the code, and omit -seed()- from your -bootstrap- command.

    Your real data set will need to contain at least 2,500 T509 == 1 observations in order to allow bootstrap with a size as large as 2500. The bootstrap resample size must be less than or equal to the number of observations in the data.

    Although it is perfectly acceptable to -preserve- and -restore- your data set each time you execute -myboot-, I don't see any need for that here. Your code in -myboot- does not destroy any information; all it does is generate some new variables. If you just -drop- those newly created variables just before the -end- statement, then you can eliminate the -preserve-/-restore- dance, and by doing this you will greatly speed up execution. (Preserve writes the file to disk, and restore reads it. So, with them, you will be spending a lot of time doing file input-output operations, which are painfully slow.)

    In developing and testing this code, I recommend that you reduce -reps(500)- with something fairly small like -reps(10)-. The results of a 10-replication bootstrap aren't really useful, but you will be able to run the program much faster to find out if it is working or to find where it breaks. Then once it works, going back to -reps(500)- it almost certainly will still work. Similarly, I would reduce the -forvalues- loop to something like -forvalues i = 500(1000)2500- during development and testing. If it runs properly for those values of i, it will run with more values in between filled-in. So change the loop back to the actual values of i you want only after you have code that runs without difficulty.

    There may well be other problems with the code. But at a minimum you will need to fix these things. And if it does not perform as expected after you do, it may be easier for you to identify and fix whatever problems remain. If you need more help to do that, when you post back, please do not say it "doesn't work." Be sure to show the new code you are using (N.B. exactly!) and be very explicit about what precisely you got and why it isn't what you were hoping for. Ideally, show the actual Stata output as copied from the Results window, and an example of the data from your results data set.

    Comment


    • #3
      Thank you so much for your guidance, Prof. Schechter. I have the code running now but I need some help getting the right information into my new data. This is the code that I have running now:


      capture program drop myboot
      program define myboot, rclass
      preserve
      bsample

      probit ICD_Opioids Sex Black Home age030 if (T509~=1)
      //Do not bootstrap SEs here, since we don't need them and will slow down code.

      predict death_hat if T509==1, pr
      gen deathpredicted=1 if death_hat>=.5
      recode deathpredicted .=0
      replace deathpredicted = . if T509 !=1

      egen numdeaths_ruhm_tmp = sum(deathpredicted) // Calculate estimated number of opioid deaths with Ruhm approach
      sum numdeaths_ruhm_tmp // Summarize this--will extract mean in next line
      return scalar numdeaths_ruhm = r(mean)

      egen numdeaths_tox_tmp = sum(T509_TOXOPIOID) // Calculate estimated number of opioid deaths from toxicology data
      sum numdeaths_tox_tmp // Summarize this--will extract mean in next line
      return scalar numdeaths_tox = r(mean)

      return scalar numdeaths_diff = numdeaths_ruhm - numdeaths_tox // Main statistic of interest. Is it different from 0?
      sum numdeaths_diff

      restore
      end


      forvalues i=500(10)1169 {
      local rand_seed = 1234 + `i'
      set seed `rand_seed'
      bootstrap numdeaths_diff = r(numdeaths_diff) `i', size(`i') /*saving(bootstrap10, replace)*/ reps(500): myboot
      append using "C:\Users\Sumedha\Documents\OPIOIDS DiD\RUHMvsTOX\RuhmVsTox_MicroMacromerged.dta"
      save "C:\Users\Sumedha\Documents\OPIOIDS DiD\RUHMvsTOX\RuhmVsTox_MicroMacromerged.dta", replace
      }



      Right now, what I get from running the code is the three variables -
      numdeaths_diff,
      numdeaths_ruhm , and numdeaths_tox - one for each of the 500 reps
      in the new data. Addtionally, I get the following output table analyzing these estimates :

      Bootstrap results Number of obs = 1,169
      Replications = 499

      command: myboot
      numdeaths_d~f: r(numdeaths_diff)
      numdeaths_r~m: r(numdeaths_ruhm)
      numdeaths_tox: r(numdeaths_tox)


      Observed Bootstrap Normal-based
      Coef. Std. Err. z P>z [95% Conf. Interval]

      numdeaths_diff 72 22.35286 3.22 0.001 28.1892 115.8108
      numdeaths_ruhm 624 29.72488 20.99 0.000 565.7403 682.2597
      numdeaths_tox 552 22.97391 24.03 0.000 506.972 597.028

      Note: One or more parameters could not be estimated in 1 bootstrap replicate;
      standard-error estimates include only complete replications.

      Instead of the new data with
      numdeaths_diff,
      numdeaths_ruhm , and numdeaths_tox for each rep
      , I only want the new data to record for each `i' the following 5 variables from the output table above -
      numdeaths_diff, numdeaths_ruhm , numdeaths_tox, `i' and z score of numdeaths_diff. Is there a way to return those values to the new dataset instead, please?

      I am very grateful for your input.
      Sincerely,
      Sumedha.

      Comment


      • #4
        It is easy enough to include `i' in the output: just add -gen i = `i'- right before the -append- command.

        But I don't quite get what you are asking about the z score of numdeaths_diff. First of all, there are no z-scores lying around here. Do you mean the z-statistic that Stata gives you in the summary output? Or do you actually want to standardize the numdeaths_diff variable and calculate z-scores for them. The z-statistic would be a single number that applies to the entire 500 (or 499 in this case) observations, whereas z-scores would take on (usually) distinct values for each observation.

        Comment


        • #5
          Thank you for your response, Prof. Schechter. In the new dataset I create, I ONLY want to include the information from the summary output. I do not want the entire 499 observations. I want the
          numdeaths_diff, numdeaths_ruhm , numdeaths_tox, `i' and z statistic (not score) of numdeaths_diff from the summary output. So, if my loop is for values of `i' between 500(10)2000 then I should have 150 observations in all. Sorry for not writing it correctly.

          Comment


          • #6
            I have also made numerous modifications to your program -myboot-. I believe it will produce the results you are looking for. Much of the calculation in the original version is just wasted effort, and you can get the same results much quicker with less code, or faster versions of it. Also, while I have removed all the -egen- commands anyway as all they do is find the average value of a constant, for future reference the -egen, sum()- function has been renamed -egen, total()-. This was done to eliminate confusion with the -gen, sum()- function which does something rather different. While -egen, sum()- is still recognized, it is best to avoid using it. I have removed the -preserve- and -restore- and just dropped the new variables created.

            I have moved the -set seed- command outside the loop: setting the random number seed should not be re-done each time through the loop; that weaken the randomness of the sequence you generate. Random number seeds: set 'em (once) and forget 'em.

            The following code shows the approach to getting what you want, but is not tested. It does not run with your -dataex- example for the same reasons stated earlier in the thread: there are variable name incompatibilities and the sample is too small for the analysis.

            Code:
            clear*
            tempfile results
            save `results', emptyok
            
            //   HERE YOU MUST READ IN THE DATA: use whatever
            
            capture program drop myboot
            program define myboot, rclass
            
                probit ICD_Opioids Sex Black Home age030 if (T509~=1)
                //Do not bootstrap SEs here, since we don't need them and will slow down code.
            
                predict death_hat if T509==1, pr
                gen deathpredicted=1 if death_hat>=.5
                recode deathpredicted .=0
                replace deathpredicted = . if T509 !=1
            
                sum deathpredicted, meanonly // Summarize this--will extract total in next line
                scalar numdeaths_ruhm = r(sum)
                return scalar numdeaths_ruhm = numdeaths_ruhm
            
                sum T509, meanonly // Summarize this--will extract total in next line
                scalar numdeaths_tox = r(sum)
                return scalar numdeaths_tox = numdeaths_tox
            
                return scalar numdeaths_diff = numdeaths_ruhm - numdeaths_tox // Main statistic of interest. Is it different from 0?
            
                drop deathpredicted death_hat
            
            end
            
            set seed 1234 // OUTSIDE LOOP; SET IT ONCE AND FORGET IT!!!
            
            forvalues i=500(10)1169 {
                bootstrap numdeaths_diff = r(numdeaths_diff), size(`i')  reps(500) noisily: myboot
                matrix M = r(table)
                gen z_stat = M[3, 1]
                gen i = `i'
                append using `results'
                save `"`results'"', replace
            } 
            
            use `results', clear
            The temporary file `results' will now contain, one observation for each value of i, that observation containing numdeaths_diff, `i', and the associated z-statistic. You can save it as "C:\Users\Sumedha\Documents\OPIOIDS DiD\RUHMvsTOX\RuhmVsTox_MicroMacromerged.dta" at the end.

            If you need to post back and show your updated code, please remember to use code delimiters, as I have done here, around it so that the indentation (and readability) is preserved. If you are not familiar with code delimiters, read Forum FAQ #12. (If you are not indenting your code, that is a habit you should take up immediately. Indentation makes code much more readable and understandable, which means fewer mistakes, and easier debugging, testing, and maintenance.)




            Comment


            • #7
              Dear Prof. Schechter,

              You are so kind! Thank you for the code and detailed explanation. I will keep in mind the correct way to include code.

              Sorry for continuing to request help but I have one last additional question. How can I also include some sensitivity and specificity measures in this bootstrap program? I know a simple estat classification will do the job after the probit but I am not sure how to include that in the bootstrap program. With my older version of the boostrap, one prior to your much cleaner code I got the following error:

              g Male = 1 if Sex == 1
              (416 missing values generated)

              . replace Male = 0 if Sex == 0
              (416 real changes made)

              . g MalexBlack = Male*Black

              . capture program drop myboot

              . program define myboot, rclass
              1. preserve
              2. bsample
              3.
              . probit ICD_Opioids Sex Black Home age030 if (T509~=1)
              4. //Do not bootstrap SEs here, since we don't need them and wi
              > ll slow down code.
              .
              .
              . predict death_hat if T509==1, pr
              5. gen deathpredicted=1 if death_hat>=.5
              6. recode deathpredicted .=0
              7. replace deathpredicted = . if T509 !=1
              8.
              . egen numdeaths_ruhm_tmp = sum(deathpredicted) // Calculate esti
              > mated number of opioid deaths with Ruhm approach
              9. sum numdeaths_ruhm_tmp // Summarize this--will extract mean
              > in next line
              10. return scalar numdeaths_ruhm = r(mean)
              11.
              . egen numdeaths_tox_tmp = sum(T509_TOXOPIOID) // Calculate estim
              > ated number of opioid deaths from toxicology data
              12. sum numdeaths_tox_tmp // Summarize this--will extract mean i
              > n next line
              13. return scalar numdeaths_tox = r(mean)
              14.
              . return scalar numdeaths_diff = numdeaths_ruhm - numdeaths_tox /
              > / Main statistic of interest. Is it different from 0?
              15.
              .
              . estat classification
              16. return scalar corr=r(P_corr)
              . return scalar sensitivity=r(P_p1)
              17. return scalar specificity=r(P_n0)
              18. restore
              19. end

              .
              . bootstrap numdeaths_diff = r(numdeaths_diff) numdeaths_ruhm = r(numdeaths_ruhm)
              > numdeaths_tox = r(numdeaths_tox) corr=r(P_corr) sensitivity=r(P_p1) specificity=r(P_n0
              > ), saving(bootstrap10, replace) reps(500) seed(1234): myboot
              (running myboot on estimation sample)
              'r(P_corr)' evaluated to missing in full sample
              r(322);
              Not sure why it will not work since it does work if I run it after the probit without the bootstrapping.

              I am sorry to keep bothering you ... but I am very grateful for your help and advise.
              Sincerely,
              Sumedha.

              Comment


              • #8
                I'm really not sure. It may be that you need to run -estat classification- and the -return scalar- commands associated with it closer to the -probit- command itself. Perhaps the intervening commands are wiping out something that -estat classification- needs to run properly. But I can't actually point to a command that looks like it will do that.

                I see that you have put back most of the things into your program myboot that I took out. You are, of course, free to do that, but I don't see why you are doing that. Everything I removed from the code is completely unnecessary and just gums up the works. It makes the code longer, harder to understand, and slows things down (inside a loop, where speed actually matters). Also, including -bsample- within program myboot is actually an error: -bootstrap- itself calls -bsample-. You not only do not need to do it explicitly within your program, you will get incorrect results from doing that because you are causing Stata to now do another re-sampling from the original re-sample, so the distributions of results will not be right.

                I suggest that you run the code down through the definition of program myboot, and then run

                Code:
                myboot
                Stata will then attempt to run your program on the full data set. The difference from running it under -bootstrap-, however, is that you will get to see the output from all of the commands along the way, and you may be able to see why -estat ic- failed to produce results for you.

                Comment


                • #9
                  Dear Prof. Schechter,

                  I had not had a chance to run your code yet and hence I asked about the estat classification using the older code. I just ran your code, with estat classification right after running the probit, but am still out of luck
                  clear*

                  . tempfile results

                  . save `results', emptyok
                  (note: dataset contains 0 observations)
                  file C:\Users\Sumedha\AppData\Local\Temp\ST_01000001.tm p saved

                  .
                  . use "C:\Users\Sumedha\Documents\OPIOIDS DiD\RUHMvsTOX\RuhmVsTox_MicroMacromerged.dta",
                  > clear

                  .
                  . capture program drop myboot

                  . program define myboot, rclass
                  1.
                  . probit ICD_Opioids Sex Black Home age030 if (T509~=1)
                  2. //Do not bootstrap SEs here, since we don't need them and will slow down code.
                  . estat classification
                  3. return scalar corr=r(P_corr)
                  4. return scalar sensitivity=r(P_p1)
                  5. return scalar specificity=r(P_n0)
                  6.
                  . predict death_hat if T509==1, pr
                  7. gen deathpredicted=1 if death_hat>=.5
                  8. recode deathpredicted .=0
                  9. replace deathpredicted = . if T509 !=1
                  10.
                  . sum deathpredicted, meanonly // Summarize this--will extract total in next line
                  11. scalar numdeaths_ruhm = r(sum)
                  12. return scalar numdeaths_ruhm = numdeaths_ruhm
                  13.
                  . sum T509, meanonly // Summarize this--will extract total in next line
                  14. scalar numdeaths_tox = r(sum)
                  15. return scalar numdeaths_tox = numdeaths_tox
                  16.
                  . return scalar numdeaths_diff = numdeaths_ruhm - numdeaths_tox // Main statistic of
                  > interest. Is it different from 0?
                  17.
                  . drop deathpredicted death_hat
                  18.
                  . end

                  .
                  . set seed 1234 // OUTSIDE LOOP; SET IT ONCE AND FORGET IT!!!

                  .
                  . forvalues i=500(10)1169 {
                  2. bootstrap numdeaths_diff = r(numdeaths_diff) corr=r(P_corr) sensitivity=r(P_p1)
                  > specificity=r(P_n0), size(`i') reps(500) noisily: myboot
                  3. matrix M = r(table)
                  4. gen z_stat = M[3, 1]
                  5. gen i = `i'
                  6. append using `results'
                  7. save `"`results'"', replace
                  8. }
                  bootstrap: First call to myboot with data as is:

                  . myboot

                  Iteration 0: log likelihood = -239.56588
                  Iteration 1: log likelihood = -218.40066
                  Iteration 2: log likelihood = -218.07161
                  Iteration 3: log likelihood = -218.0713
                  Iteration 4: log likelihood = -218.0713

                  Probit regression Number of obs = 497
                  LR chi2(4) = 42.99
                  Prob > chi2 = 0.0000
                  Log likelihood = -218.0713 Pseudo R2 = 0.0897


                  ICD_Opioids Coef. Std. Err. z P>z [95% Conf. Interval]

                  Sex .4084474 .1424403 2.87 0.004 .1292695 .6876253
                  Black -.4821259 .1597527 -3.02 0.003 -.7952355 -.1690164
                  Home .3079318 .1387389 2.22 0.026 .0360086 .579855
                  age030 .7664434 .1762885 4.35 0.000 .4209242 1.111963
                  _cons .3812713 .1486202 2.57 0.010 .0899811 .6725615


                  Probit model for ICD_Opioids

                  True --------
                  Classified D ~D Total

                  401 86 487
                  3 7 10

                  Total 404 93 497

                  Classified + if predicted Pr(D) >= .5
                  True D defined as ICD_Opioids != 0

                  Sensitivity Pr( + D) 99.26%
                  Specificity Pr( -~D) 7.53%
                  Positive predictive value Pr( D +) 82.34%
                  Negative predictive value Pr(~D -) 70.00%

                  False + rate for true ~D Pr( +~D) 92.47%
                  False - rate for true D Pr( - D) 0.74%
                  False + rate for classified + Pr(~D +) 17.66%
                  False - rate for classified - Pr( D -) 30.00%

                  Correctly classified 82.09%

                  (497 missing values generated)
                  (6 missing values generated)
                  (deathpredicted: 6 changes made)
                  (497 real changes made, 497 to missing)
                  'r(P_corr)' evaluated to missing in full sample
                  r(322);

                  end of do-file
                  Any further insights? It is giving the sensitivity, specificity etc. but having trouble returning those estimates to scalar?

                  Thank you, again.
                  Sincerely,
                  Sumedha.

                  Comment


                  • #10
                    Try adding the -noisily- option to the -bootstrap- command so that we can see more of what's going on there.

                    Comment


                    • #11
                      Its already there....

                      Comment


                      • #12
                        Sorry, you're right. Oh, it's clear. -myboot- does not return a parameter r(P_corr). Rather, it returns a parameter called r(corr) which is calculated as equal to the r(P_corr) that -estat classification- returns. There is a similar problem with your sensitivity and specificity. -bootstrap- does not reach in to the scalars returned by -estat classification-, it only reaches for scalars returned by -myboot-, and you gave them different names. So your -bootstrap- command needs to read:

                        Code:
                        bootstrap numdeaths_diff = r(numdeaths_diff) corr=r(corr) sensitivity=r(sensitivity) ///
                            specificity=r(specificity), size(`i') reps(500) noisily: myboot

                        Comment


                        • #13
                          You are the best! Thank you, Prof. Schechter.

                          Comment

                          Working...
                          X