Stata The Stata listserver
[Date Prev][Date Next][Thread Prev][Thread Next][Date index][Thread index]

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


From   "Steichen, Thomas J." <SteichT@rjrt.com>
To   <statalist@hsphsun2.harvard.edu>
Subject   st: RE: RE: RE: Using -byable- and the -by- commands correctly?
Date   Mon, 31 Oct 2005 11:32:25 -0500

> 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/



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