saving scalars in appended file after running simulations

[Mark Oyama]

Hi all, newbie to this listserve--hoping to figure out how to append scalar values to a data file after running multiple time to event simulations.

Currently, following stcox, I generate a new variable from the scalar matrix that I call coeff. Would like to save this value to a separate data file and then when I run the simulation again, append the new scalar value to the separate data file. Objective is to run the simulation ~100 times and end up with a data set of 100 scalar values that I can then average.

set obs 100
gen t25 = -1/.99 * log(1-uniform())
gen t99 = -1/.99 * log(1-uniform())
gen time = min(t99,t25)
gen byte fail = (t25<t99) + 2*(t25>=t99)
replace fail = 0 if time>=2
replace time = 2 if time>2
set obs 200
gen treat = _n>100
*change hazard for treated cohort
*change hazard for first event from 0.25 to 0.1875 (HR=0.75)
replace t25 = -1/.5 * log(1-uniform()) in 100/l
*change HR for second (competing event) from 0.99 to 0.495 (HR=0.50)
replace t99 = -1/.99 * log(1-uniform()) in 100/l
replace time = min(t99,t25)
replace fail = (t25<t99) + 2*(t25>=t99)
replace fail = 0 if time>=2
stset time, f(fail==1) noshow
stcox treat
matrix b=e(b)
svmat b, name(coeff)

**Would like to save the scalar matrix value e(b) to dataset at variable coeff, then send this value to a different file (postfile??)
**Then run simulation again and save next coeff values to the new file etc...

replace t25 = -1/.99 * log(1-uniform()) in 1/100
replace t99 = -1/.99 * log(1-uniform()) in 1/100
replace time = min(t99,t25) in 1/100
replace fail = (t25<t99) + 2*(t25>=t99) in 1/100
replace fail = 0 if time>=2 in 1/100
replace time = 2 if time>2 in 1/100

*change hazard for treated cohort
*change hazard for first event from 0.25 to 0.1875 (HR=0.75)
replace t25 = -1/.5 * log(1-uniform()) in 100/l
*change HR for second (competing event) from 0.99 to 0.495 (HR=0.50)
replace t99 = -1/.99 * log(1-uniform()) in 100/l
replace time = min(t99,t25)
replace fail = (t25<t99) + 2*(t25>=t99)
replace fail = 0 if time>=2
stset time, f(fail==1) noshow
stcox treat
matrix b=e(b)

[Charlie Joyez]

Mark,
Try the following (it uses some Mata, perhaps it would have been possible without, but learning to use some Mata is also great.)
There is certainly proper way to do this, but I didn't wanted to modify much the core of your code, so you can understand what happens (what I added is in bold):

local nbloop=100 /*Choose your nb of iterations*/
mata matcoeff=J(`nbloop',1,.) /*defines an empty Mata matrix (vector) of the same dimension*/


forvalues n=1/`nbloop' {
clear
set obs 100
gen t25 = -1/.99 * log(1-uniform())
gen t99 = -1/.99 * log(1-uniform())
gen time = min(t99,t25)
gen byte fail = (t25<t99) + 2*(t25>=t99)
replace fail = 0 if time>=2
replace time = 2 if time>2
set obs 200
gen treat = _n>100
*change hazard for treated cohort
*change hazard for first event from 0.25 to 0.1875 (HR=0.75)
replace t25 = -1/.5 * log(1-uniform()) in 100/l
*change HR for second (competing event) from 0.99 to 0.495 (HR=0.50)
replace t99 = -1/.99 * log(1-uniform()) in 100/l
replace time = min(t99,t25)
replace fail = (t25<t99) + 2*(t25>=t99)
replace fail = 0 if time>=2
stset time, f(fail==1) noshow
stcox treat
return list
matrix b=e(b)
svmat b, name(coeff)

local coeff=coeff[1]
mata matcoeff[`n',]=`coeff'

}
clear
set obs `nbloop'
mata st_addvar("double","est_coeff")
mata st_store(.,"est_coeff",matcoeff)

I simply store the coefficient into a local, then this local into a mata matrix, and loop this for 100 (but you can easily change it) iterations.
At the end I convert the mata matrix into a Stata variable, with the 100 estimated coefficient.

Hope this helps,
Charlie

[Mark Oyama]

Hi Charlie,
This is awesome. Thanks for the help and intro to Mata!!
Mark

[Charlie Joyez]

You're welcome,
I'm glad this correspond to what you wanted, I wasn't sure at first.

I don't know what you're looking at, but if you want to see how this coefficient should be when randomly drawn, I first advice you to increase the number of iteration to 1000 at least (it doesn't take much time to run, so it should be ok), it should lower the variance, and give you a mean closer to the "real" mean of a null model.

Don't hesitate to post again,
Best,
Charlie