Bootstrap own program

[Jonas Weinert]

Dear Community,
I need to bootstrap a set of commands. Namely,
1. Run a regression
2. Obtain fitted values
3. Obtain mean from fitted values

I am doing thi in the context of medical cost estimation. Unfortunately, my code does not work:
Running the simulate command returns only empty values. I am guessing that somewhere in Step 3 (my program), the storing of the mean from my fitted values does not work.

Code:

* Step 1

quietly glm cost time if death==1, link(log) family(gamma)
predict GLM2
sum GLM2
return list

matrix GLMmean = (r(mean))
matrix list GLMmean



* Step 2: program
capture program drop GLMboot
program define GLMboot, rclass
preserve
    bsample
        glm cost time if death==1, link(log) family(gamma)
        predict GLMbootfit2
        sum GLMbootfit2
        return scalar GLMmean_1 = r(mean)
restore        
end

* Step 3: Monte Carlo 
preserve
simulate GLMmean_1=r(mean), reps(10) seed(12345): GLMboot, rclass
bstat,  n(2000)
estat bootstrap, all

Thank you!

[Ken Chui]

Welcome to Statalist!

I can find two mistakes if that helps. But also read the whole thing before doing any revision.

1. I have never seen specifying a simulation command with options like ", rclass bstat, n(2000)". I tried that and got an error. Generally, within one level of command, if you're typing two or more commas, it's usually against the grammar. I'd suggest moving that 2000 up to the line 4 of the program instead, as "bsample 2000". And take away all the extra options after the "GLMboot" in the simulate statement.
2. In your program, you named the rclass object "r(GLMmean_1)". But in the simulate statement, you asked for a non-existing "r(mean)". Try covert that back to r(GLMmean_1).

First, a disclaimer is that I am not very well-versed in the bootstrap command, but are you working off someone else's code that has been tested? I'm a bit confused by the approach. You ended with an estat boostrap, but the command before that was not bootstrap. It feels somewhat like this was duct taped together. I'd suggest check out help bootstrap, which allows both setting repetition (rep()) and sample size (size()), and it's very close to the concept of simulate that it takes user-written program as well.

clear
sysuse auto, clear
keep price weight foreign
expand 1000

capture program drop GLMboot
program define GLMboot, rclass
    glm price weight if foreign==1, link(log) family(gamma)
    capture drop GLMbootfit2
    predict GLMbootfit2
    sum GLMbootfit2
    return scalar GLMmean_1 = r(mean)
end

bootstrap testvar = r(GLMmean_1), size(2000) rep(10) saving(output, replace): GLMboot

[Jonas Weinert]

Thank you very much Ken!
Your code works perfectly and your annotations make a lot of sense!