Re: st: capturing scalars in simulation

[Richard Hiscock <richardjhiscock@gmail.com>]

Nick,
once again problem solved & understanding increased
thanks Richard



On 23 Feb 2013, at 20:51, Nick Cox <njcoxstata@gmail.com> wrote:

> In your program the first scalar you pick up from -concord- is
> r(LOA_ll) and you chose to call it r(ll). Similarly in other (although
> not all) cases you gave scalars different names from those used by
> -concord-.
> 
> You can do that, but then you are committed to using the names you
> gave elsewhere, not the names used by -concord-. This is because your
> program as an r-class program overwrites previous r-class results.
> 
> it would be simpler just to run -concord- within your program and let
> its r-class results be those you pick up. Otherwise what you are doing
> is taking something out of one bag and putting it in another bag,
> which is not needed. That way round, -mysim1- is not rclass.
> -simulate- would just be picking up the results emitted by -concord-
> without the need for re-bagging.
> 
> A quite different comment is that within -mysim1- it seems that a lot
> of data management is done again and again, each time round the
> simulation. I realise that you are taking a different sample each
> time, but the rest of the management looks as if it could and should
> be taken outside the program, not least because doing this for real
> would, as you realise, require many more than 10 repetitions.
> 
> capture program drop mysim1
> *! mysim1 Bland-Altman LOA from random pairs RJH 23021013
> program define mysim1
> version 12.1
> drop _all
> 
> use  "file location",
> ** data in long format with id meth replicate
> 
> *obtain random replicate for each of methods 2 & 3 (meth 1 single value)
> sort id meth
> by id meth : sample 1,count
> 
> ** obtain requred estimates using concord program
> quietly reshape wide hb repl, i(id) j(meth)
> rename hb1 meth1
> rename hb2 meth2
> rename hb3 meth3
> 
> keep id meth*
> concord  meth2 meth1,
> end
> 
> simulate ll =r(LOA_ll) ul =r(LOA_ul) diff =r(diff) sd =r(sd_diff) n
> =r(N) , reps(10) noisily : mysim1
> 
> 
> On Sat, Feb 23, 2013 at 8:48 AM, Richard Hiscock
> <richardjhiscock@gmail.com> wrote:
>> Hi all,
>> I have method comparison data with haemoglobin level measured by three methods (pathology (1 reading), meth1 (3 readings) & meth2 (3 readings)) at the same time on the same subjects. Replicates are considered exchangeable within subject/method.
>> 
>> As part of the analysis I am trying to create a dataset with a single reading for meth1 & meth2 which have been randomly sampled from within replicates for the subject & method.
>> Then using concord[SJ 10-4 by T Steichen & N Cox] to obtain scalars for the mean difference & its SD + upper & lower limits of agreement
>> 
>> I then wish to simulate this so as to obtain percentile based CI for these scalars (say 1000 reps)
>> I have written the following program which successfully creates the datasets and when using noisily the analysis of each dataset occurs but  I am unable to capture the scalars produced by concord except for one (the mean difference [listed as r(diff)] in concord program help file).
>> 
>> Any guidance would be appreciated.
>> 
>> thanks Richard Hiscock
>> 
>> **TRY RANDOM selection from replicates
>> 
>> capture program drop mysim1
>> *! mysim1 Bland-Altman LOA from random pairs RJH 23021013
>> program define mysim1, rclass
>> version 12.1
>> drop _all
>> 
>> use  "file location",
>> ** data in long format with id meth replicate
>> 
>> *obtain random replicate for each of methods 2 & 3 (meth 1 single value)
>> sort id meth
>> by id meth : sample 1,count
>> 
>> ** obtain requred estimates using concord program
>> quietly reshape wide hb repl, i(id) j(meth)
>> rename hb1 meth1
>> rename hb2 meth2
>> rename hb3 meth3
>> 
>> keep id meth*
>> concord  meth2 meth1,
>> return scalar ll =r(LOA_ll)
>> return scalar ul =r(LOA_ul)
>> return scalar diff =r(diff)
>> return scalar sd =r(sd_diff)
>> return scalar n =r(N)
>> end
>> 
>> simulate ll =r(LOA_ll) ul =r(LOA_ul) diff =r(diff) sd =r(sd_diff) n =r(N) , reps(10) noisily : mysim1
> 
> *
> *   For searches and help try:
> *   http://www.stata.com/help.cgi?search
> *   http://www.stata.com/support/faqs/resources/statalist-faq/
> *   http://www.ats.ucla.edu/stat/stata/


*
*   For searches and help try:
*   http://www.stata.com/help.cgi?search
*   http://www.stata.com/support/faqs/resources/statalist-faq/
*   http://www.ats.ucla.edu/stat/stata/