What happens to ProbCE object that you create here? That might be what
you tried to rbind, rather than INB in your last line of the cycle. I
would assume this is what you want to do.
As always, there are multiple ways of getting what you need there. You
were given a Mata solution, which looks an overkill to me, and hardly
is directly usable in further analysis. If you want a Stata variable
as a result, there are two ways of going there.
1. A somewhat clumsy way of trying to have everyting in the same data
set: assume that you have -eff- and -cost- variables for each
observation. Then you could do this:
capture noisily set obs 1000
* might give you an error message if you already have more than 1000
observations;
* we might want to capture that error message
gen INBsign = .
gen ceas = .
* this is what we will fill in
gen n = _n in 1/1000
* _n is the notation for the current observation number
* the range restriction limits the values to be from 1 to 1000
forvalues n=1/1000 {
* n is what is known as a local macro in Stata jargon;
* has a special reference format with opening and closing quotes
replace INBsign = (eff*`n' - cost>0)
sum INBsign, meanonly
replace ceas = r(mean) in `n'
* uses the return value from -summarize- command
* replaces a single value in `n'-th line of data
}
2. What I think is a more elegant solution is to create a separate
Stata data set. In your R code, you create a separate ceas object;
unfortunately Stata can only think of one data set in memory:
postfile topost n probce using ceas, replace
* declares a new file ceas.dta to be created that is to contain two
variables, n and probce
count if ~missing(eff)
local denom = r(N)
* again, it is using the return value of the previous command
forvalues n=1/1000 {
count if eff*`n'-cost > 0
post topost (`n') ( r(N)/`denom' )
* writes a new line to the ceas data set
}
postclose topost
use ceas, clear
HTH, Stas
On Thu, Aug 27, 2009 at 12:44 AM, Frank Peter<[email protected]> wrote:
> Dear All,
>
> This is to thank everybody for all their answers, I have really learn a lot from the answers.
>
> I was trying to generate cost-effectiveness acceptability curve. The R codes I was referring to were:
> ========================================
> ceac = c()
> denom = length(effectiveness)
> for (n in 1000) {
> INB = effectiveness*n - cost
> ProbCE = length(subset(INB, INB>=0))/denom
> ceac = rbind(ceac, c(n, INB))
> }
> ========================================
>
> Thanks
>
> Frank
> reference:
> Fenwick E, Marshall DA, Levy AR, Nichol G. Using and interpreting cost-effectiveness acceptability curves: an example using data from a trial of management strategies for atrial fibrillation. BMC Health Serv Res. 2006 Apr 19;6:52. (http://www.biomedcentral.com/1472-6963/6/52)
>
--
Stas Kolenikov, also found at http://stas.kolenikov.name
Small print: I use this email account for mailing lists only.
*
* For searches and help try:
* http://www.stata.com/help.cgi?search
* http://www.stata.com/support/statalist/faq
* http://www.ats.ucla.edu/stat/stata/
