[Date Prev][Date Next][Thread Prev][Thread Next][Date index][Thread index]

From |
"Nick Cox" <n.j.cox@durham.ac.uk> |

To |
<statalist@hsphsun2.harvard.edu> |

Subject |
st: RE: RE: RE: Using -byable- and the -by- commands correctly? |

Date |
Mon, 31 Oct 2005 15:48:01 -0000 |

With one thing and another I think I can reduce your program to this, although it's not tested: program poissontt, byable(recall, noheader) version 9 quietly { args observed expected toreplace if `observed' == 0 local maxprob = 1 else if `expected' < 0.000001 local maxprob = 0.000001 else { forval I = 0 /`= `observed' - 1' { local tempcalc = /// ((exp(-`expected')) * (`expected'^`I')) / exp(lnfactorial(`I')) local sumcalc = `sumcalc' + `tempcalc' } local maxprob = 1 -`sumcalc' } forval I = 0/`observed' { local tempcalc = /// ((exp(-`expected')) * (`expected'^`I')) / exp(lnfactorial(`I')) local minprob = `minprob' + `tempcalc' } if `observed' == 0 & `expected' == 0 local poissval = 0 else if `minprob' <= `maxprob' local poissval = `minprob' else local poissval = `maxprob' replace `toreplace' = max(2 * `poissval', 1) } end I suspect there's scope for further simplification, but the dodgiest details are that there is nothing in your internal syntax that looks like a standard byable operation and that your last -replace- is for the whole of the variable. I doubt that's what you want. A broader comment is that your approach of making this byable and then looping over the data does not look the best strategy to me, as calculations look as if they can and should be recast as calculations on variables. Finally, this is not my usual territory, but biostat people may recognise the problem as one already solved. Nick n.j.cox@durham.ac.uk David Winter > please see below as requested. This program is based on Fortran > functions incorporated in a Person-Years analysis program. > > capture program drop poissontt > program poissontt, sortpreserve byable(recall, noheader) > /* make program usable with 'by' command - see manual [P] page 12 */ > version 9 > quietly > local poissval = 0 > local minprob = 0 > local maxprob = 0 > local observed = `1' /* O */ > local expected = `2' /* E */ > local replacevar = "replace `3' = " > /* probability Observed no >= O */ > if `observed' == 0 { > local maxprob = 1 > } > else if `expected' < 0.000001 { > local maxprob = 0.000001 > } > else { > local I = 0 > local sumcalc = 0 > local tempcalc = 0 > local endloop = `observed' - 1 > while `I' <= `endloop' { > local tempcalc = ((exp(-1 * `expected')) * > (`expected' ^ `I')) / exp(lnfactorial(`I')) > local sumcalc = `sumcalc' + `tempcalc' > local I = `I' + 1 > } > local maxprob = 1 -`sumcalc' > } > /* probability Observed no <= O */ > local I = 0 > local sumcalc = 0 > local tempcalc = 0 > local endloop = `observed' > while `I' <= `endloop' { > local tempcalc = ((exp(-1 * `expected')) * (`expected' ^ > `I')) / exp(lnfactorial(`I')) > local sumcalc = `sumcalc' + `tempcalc' > local I = `I' + 1 > } > local minprob = `sumcalc' > /* evaluate probsbility to calculate statistic */ > if `observed' == 0 & `expected' == 0 { > local poissval = 0 /* empty */ > } > else { > if `minprob' <= `maxprob' { > local poissval = `minprob' > } > else { > local poissval = `maxprob' > } > local poissval = 2 * `poissval' > if `poissval' > 1 { > local poissval = 1 > } > } > `replacevar' `poissval' > noisily > end > > The following 4 records are part of a 27 record dataset sorted by > > ageband (5,10,15,20..85) within sex (MALE,FEMALE). > > > > Sex Ageband Observed Expected PoissStat > > MALE 5 5 .2585 > > MALE 10 11 .6054 > > MALE 15 6 .6516 > > MALE 20 5 .6435 > > > > I have written a program to calculate a two-tailed poisson > statistic, > > that has three arguements: > > > > `1' = Observed number > > `2' = Expected number > > `3' = variable to hold the poisson statistic > > > > I have made the program -byable- : > > > > Program poissontt, sortpreserve byable(recall, no header) > > .. > > .. > > local poissval = 0 /* macro to store the result */ > > .. > > local observed = `1' > > local expected =`2' > > local replacevar = "replace `3' = " > > .. > > <define poissval> > > .. > > `replacevar' `poissval' > > end > > > > In a -.do- file, I run the commands > > > > gen PoissStat = . > > by sex ageband: poissontt Observed Expected PoissStat > > > > Stata diligently loops through the data and using -set trace- I have > > tracked the iterations. It does not appear to pick up the values of > > Observed and Expected for each record, rather those of the > > first record > > only. The calculation is correct for record 1 (Obs=5, Expect=.258, > > poisson(two-tailed=.0000153) and the PoissStat variable is > updated in > > each record, but with the value for record 1 only. > > > > Is there something I am missing either in the program or in > using the > > -by- command? * * For searches and help try: * http://www.stata.com/support/faqs/res/findit.html * http://www.stata.com/support/statalist/faq * http://www.ats.ucla.edu/stat/stata/

- Prev by Date:
**Re: st: Access to results after using -adjust-** - Next by Date:
**st: Help about Ranksum - porder option** - Previous by thread:
**st: RE: RE: Using -byable- and the -by- commands correctly?** - Next by thread:
**st: RE: RE: RE: Using -byable- and the -by- commands correctly?** - Index(es):

© Copyright 1996–2014 StataCorp LP | Terms of use | Privacy | Contact us | What's new | Site index |