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

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


From   "David Winter" <d.l.winter@bham.ac.uk>
To   <statalist@hsphsun2.harvard.edu>
Subject   st: RE: RE: Using -byable- and the -by- commands correctly?
Date   Mon, 31 Oct 2005 12:17:36 -0000

Certainly Nick, 

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	

-----Original Message-----
From: owner-statalist@hsphsun2.harvard.edu
[mailto:owner-statalist@hsphsun2.harvard.edu] On Behalf Of Nick Cox
Sent: 31 October 2005 10:39
To: statalist@hsphsun2.harvard.edu
Subject: st: RE: Using -byable- and the -by- commands correctly?

Please show us the whole program. 

Nick 
n.j.cox@durham.ac.uk 

David Winter
 
> Firstly, thanks to Svend Juul for the helpful comments on survival
> analysis.
> 
> 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/

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