Very many thanks to Tom (Steichen), Nick (Cox) and Svend (Juul) for their help on resolving this issue.
Dave
-----Original Message-----
From: [email protected]
[mailto:[email protected]]On Behalf Of Steichen,
Thomas J.
Sent: 31 October 2005 16:32
To: [email protected]
Subject: st: RE: RE: RE: Using -byable- and the -by- commands correctly?
> David Winter writes (in part):
> > I have written a program to calculate a two-tailed poisson statistic,
> >
> > 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?
Yes, you failed to use mark or marksample in the program.
As the help for byable indicates:
"the sample is automatically restricted to the appropriate subsample
when myprog uses the mark or marksample commands; see help mark."
With 3 changes, the program becomes:
program poissontt, sortpreserve byable(recall, noheader)
/* make program usable with 'by' command - see manual [P] page 12 */
version 8
marksample touse /* added to allow subsetting data */
quietly
gsort - `touse' /* <- added to bring selected record to top of file so */
local poissval = 0 /* observed = `1', which is really observed = `1'[1], */
local minprob = 0 /* loads the correct O (and expected loads correct E) */
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' if `touse' /* changed to replace only in current record */
noisily
end
And you get:
.. list
+-------------------------------------------------+
| sex ageband observed expected PoissS~t |
|-------------------------------------------------|
1. | MALE 5 5 .2585 .0000155 |
2. | MALE 10 11 .6054 1.15e-10 |
3. | MALE 15 6 .6516 .000122 |
4. | MALE 20 5 .6435 .0010803 |
+-------------------------------------------------+
Tom
-----------------------------------------
CONFIDENTIALITY NOTE: This e-mail message, including any attachment(s),
contains information that may be confidential, protected by the
attorney-client or other legal privileges, and/or proprietary non-
public information. If you are not an intended recipient of this
message or an authorized assistant to an intended recipient, please
notify the sender by replying to this message and then delete it from
your system. Use, dissemination, distribution, or reproduction of this
message and/or any of its attachments (if any) by unintended recipients
is not authorized and may be unlawful.
*
* 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/
*
* 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/
