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
[email protected]
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/
