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



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