retrieve the lower and upper levels of the 95%CI after simulations

[Federica Fregonese]

Dear Statalist community,
I am working with Stata 13 on Mac. I am using a program to run simulations of Poisson regression (1,000 different samples drawn from a larger population):

program define program1, rclass
 drop _all
 use "/International Health/DBs/bdd/An2013sub2.dta"
  sample 1,c  by(m)
         poisson SMc Mc, offset(n_s)
  return scalar b1=exp(_b[Mc])
  return scalar b2=_se[Mc]
 end
  simulate b_1=r(b1) b_2=r(b2) , reps (1000): program1
  sum b_1 b_2, d
 program drop _all

After the simulation I have the median of incidence-rate ratio [sum b_1], but not its 95%CI. How can I retrieve the lower and upper levels of the 95%CI for the 1,000 simulations (so that I can have their summary statistics and present a 95%CI of my median IRR)?
Thanks a lot
federica

[Mike Lacy]

Most of the information in the results table of -poisson- and other regression commands is stored in a matrix called r(table).
(This is not documented in the help for -poisson- or even -regress-. Indeed, -search r(table)- returns nothing.) Individual elements of r(table) can be retrieved with el().

Note, though, as you want IRRs, you might use the irr option of y x, IRR-. So, your -program1- might contain something like this:

poisson SMc Mc, offset(n_s) irr
return scalar ll_b1 = el(r(table), 2,5)  // ll of Mc b1 should be row 5, column 2; check this
... etc

.
And, as a side note, do you want exp(_b[Mc]) but not exp for _se[Mc]?
Regards, Mike

[Mike Lacy]

Whoops, please ignore. I see that I misread your question, and it appears you want something other than the CIs of each simulation rep. - Mike

[Maarten Buis]

centile b_1, c(2.5 97.5) will give the 2.5th and 97.5th percentiles, which could be interpreted as your lower and upper bounds for the 95% CI (assuming the simulation is correct, which without further information I cannot judge).

[Federica Fregonese]

Thank you Mike and Maarten for very useful suggestions.

Just as follow-up note, the median of lower and upper levels of 95%CI (found using

return scalar llb1 = el(r(table), 5,1)
 return scalar ulb1 = el(r(table), 6,1)

(5 and 6 being the columns and 1 the row for the 95%CI),

are more ‘conservative’ (e.g. wider) than the distribution of irr
(with:

centile b_1, c(2.5 97.5)

,

Thanks again
federica

[Maarten Buis]

Whcih one you should choose depends on what you exactly want.

[Federica Fregonese]

I'd like to present median irr, with its 95%CI, in order to show if it is significant (or not). Therefore I think I should use the median of lower and upper limits...
Thanks

[Maarten Buis]

That is not enough to track down exactly what you want. Why did you use simulation? What is the sampling uncertainty you are trying to capture with your CI (the one from the simulation or from sampling the data from the population)? I sugest you find someone local with whom you can discuss this issue face to face.