Notice: On April 23, 2014, Statalist moved from an email list to a forum, based at statalist.org.
From | Nick Cox <njcoxstata@gmail.com> |
To | statalist@hsphsun2.harvard.edu |
Subject | Re: st: capturing scalars in simulation |
Date | Sat, 23 Feb 2013 09:51:54 +0000 |
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/