Lawrence Joseph
Bayesian statistics
|
|
Please Note: Prof. Joseph has retired. These pages are left up in case they prove useful, but the pages and software will no longer be updated. All material and software is "as is" with no guarantees of functionality or correctness. |
|
bhpd1
Sample size calculations for binomial proportions via highest posterior density intervals
Summary
bhpd1 is a FORTRAN 77 program to calculate sample size requirements
in the context of a binomial experiment. The program was written by
Lawrence Joseph, David B. Wolfson and Roxane du Berger and
developed in a SunOS 4.1.2 operating environment on a SPARC
station ELC. bhpd1 uses three different Bayesian approaches to
sample size calculations based on highest posterior density
intervals (HPD). These criteria are given by:
The successful programming of the algorithms is largely due to the following sources for which we are thankful: Abramowitz, M., and Stegun, I.A. (ed.) (1972) Handbook of Mathematical Functions New York: Dover Publications Inc. Documentation
Detailed information about the implementation of the Bayesian
approaches to sample size requirements are found in our manuscript:
Sample size calculations for binomial proportions via highest posterior density intervals by Lawrence Joseph, David B. Wolfson and Roxane du Berger. The Statistician, 1995;44(2):143-154 Questions or comments regarding the program bhpd1 or requests for a copy of the manuscript may be e-mailed to Lawrence Joseph at lawrence.joseph@mcgill.ca. Files
Once you open the package bhpd1.doc (with sh), you will note the presence
of the following files:
Ex.input README.BHPD1 bhpd.f frqsmp.f readpar.f util.f Makefile baysmp.f cover.f headers.f search.f Makefile
For convenience, we have provided a Makefile to compile and link the FORTRAN code into one executable file, called bhpd1. The command make works on a Sparc station. If you encounter unnecessary difficulties, compile and link with the following command: f77 -o bhpd1 *.f Example of execution
bhpd1 is primarily designed for interactive use. bhpd1 usually executes
in a reasonable amount of time (a few minutes). When the coverage
requested is large (close to one), or the interval length
small (close to zero), or both, bhpd1 will finish execution in a much
longer time. In the latter case, it is wise to run bhpd1 in background.
We have not yet experienced strange program behaviours in the extreme
coverage and length limits, but this occurence is not impossible.
To illustrate interactive use of bhpd1, type % bhpd1 *************************************************************************** This program calculates the sample size requirements for a Binomial Parameter using four different methods: 1. Standard Normal Theory Based Frequentist Method 2. Bayesian Average Coverage Criterion (ACC) 3. Bayesian Average Length Criterion (ALC) 4. Bayesian Worst Outcome Criterion (WOC) This program also calculates minimum possible HPD coverage probability and maximum HPD length for each given sample size, as well as the pre-posterior probability that these events will occur, according to the prior input. These criteria are described in Joseph,L., Wolfson, D. and du Berger, R. Sample size calculations for binomial proportions via highest posterior density intervals, The Statistician, to appear in 1995, 44(2). The following parameters must be input: - Beta prior parameters (a,b) for the binomial parameter. The mean of this distribution is a/(a+b), and the variance is a*b/((a+b)^2*(a+b+1)); - The desired HPD coverage probability; - The desired HPD interval length; - The point estimate for the binomial parameter (used for frequentist method only). *************************************************************************** The default values for the prior parameters are a = 1.00000 and b = 1.00000, which give a uniform distribution Type <y><enter> to accept these values, other <key> to MODIFY: y The default value for HPD coverage is = 0.95000 Type <y><enter> to accept this value, other <key> to MODIFY: n Enter HPD coverage .96 The default value for HPD length is = 0.10000 Type <y><enter> to accept this value, other <key> to MODIFY: n Enter HPD length .09 The default value for the point estimate of the binomial parameter is a/(a+b) = 0.50000, the mean of the prior distribution Type <y><enter> to accept this value, other <key> to MODIFY: y *************************************************************************** You have chosen the following parameters: Beta parameters: a = 0.10000000000D+01 b = 0.10000000000D+01 HPD coverage probability c = 0.96000000000D+00 HPD interval length w = 0.90000000000D-01 Point estimate p = 0.50000000000D+00 Type <y><enter> to CONFIRM, other <key> to further MODIFY any parameter: n Modify which parameter? Type [<a><b><c><w><p>or<all>]<enter> a Enter Beta a and Beta b 2 1.5 *************************************************************************** You have chosen the following parameters: Beta parameters: a = 0.20000000000D+01 b = 0.15000000000D+01 HPD coverage probability c = 0.96000000000D+00 HPD interval length w = 0.90000000000D-01 Point estimate p = 0.57142857143D+00 Type <y><enter> to CONFIRM, other <key> to further MODIFY any parameter: y *************************************************************************** Do you wish results output to Screen or File? Type <y><enter> for SCREEN, other <key> for FILE: y *************************************************************************** The program may take a few minutes to run. Running time varies with the input parameters. Please be patient... *************************************************************************** Given point estimate = 0.57142857143 Frequentist sample size = 511 *************************************************************************** Average Coverage Criterion Step 0 511 Step 1 255 Step 2 383 Step 3 447 Step 4 415 Step 5 399 Step 6 407 Step 7 411 Step 8 413 Step 9 412 ACC sample size = 412 For a sample size of 412, the minimum HPD coverage of fixed length 0.90000D-01 over all possible data values is 0.93381538084D+00. This interval is (0.45560132480D+00, 0.54560132480D+00). For a sample size of 412, the maximum length required for 0.96000000000D+00 equal-tailed (non HPD) coverage is 0.10055879113D+00. This interval is (0.45031918322D+00, 0.55087797435D+00). For a sample size of 412, the maximum HPD length required for 0.96000000000D+00 HPD coverage is 0.10055879082D+00. This interval is (0.45032111542D+00, 0.55087990624D+00). The event of 205 successes in 412 trials occurs with prior probability of 0.31976102838D-02. *************************************************************************** Average Length Criterion Step 0 412 Step 1 206 Step 2 309 Step 3 360 Step 4 386 Step 5 373 Step 6 379 Step 7 376 Step 8 377 Step 9 378 ALC sample size = 378 For a sample size of 378, the maximum HPD length of fixed coverage 0.96000D+00 over all possible data values is 0.10492601361D+00. This interval is (0.44819090894D+00, 0.55311692255D+00). For a sample size of 378, the minimum HPD coverage required for 0.90000000000D-01 HPD length is 0.92164407237D+00. This interval is (0.45565519843D+00, 0.54565519843D+00). The event of 188 successes in 378 trials occurs with prior probability of 0.34832345729D-02. *************************************************************************** Worst Outcome Criterion Step 0 1022 Step 1 511 Step 2 766 Step 3 638 Step 4 574 Step 5 542 Step 6 526 Step 7 518 Step 8 514 Step 9 516 Step 10 515 WOC sample size = 516 For a sample size of 516, the minimum HPD coverage of fixed length 0.90000D-01 over all possible data values is 0.96007401092D+00. This interval is (0.45548047885D+00, 0.54548047885D+00). For a sample size of 516, the maximum length required for 0.96000000000D+00 equal-tailed (non HPD) coverage is 0.89966629412D-01. This interval is (0.45549593033D+00, 0.54546255974D+00). For a sample size of 516, the maximum HPD length required for 0.96000000000D+00 HPD coverage is 0.89966629269D-01. This interval is (0.45549716615D+00, 0.54546379542D+00). The event of 257 successes in 516 trials occurs with prior probability of 0.25564038885D-02. *************************************************************************** We have included in this package a parameter file called Ex.input to illustrate background execution of bhpd1. This is especially useful when a job is suspected to require much time. % cat Ex.input n #Default Beta parameters a,b where mean=a/(a+b) 2,1.5 #Beta parameters n #Default Coverage 0.99 #Coverage n #Default interval length 0.04 #Length n #Default binomial parameter estimate 0.5 #Point estimate of binomial parameter y #Choice confirmation n #Screen output b-2-1.5-.99-.04-.5 #File name For nice background execution of bhpd1, type % nice bhpd1 < Ex.input > out & bhpd1 reads the contents of Ex.input and dumps to out all the successfully read parameter information, including the file name b-2-1.5-.99-.04-.5 . This file will contain sample size calculations for each of the Bayesian criteria. % cat -n b-2-1.5-.99-.04-.5 1 Given the following user input: 2 3 Beta parameters: a = 0.20000000000D+01 4 5 b = 0.15000000000D+01 6 7 HPD coverage probability c = 0.99000000000D+00 8 9 HPD interval length w = 0.40000000000D-01 10 11 12 *************************************************************************** 13 14 Given point estimate = 0.50000000000 15 Frequentist sample size = 4147 16 17 *************************************************************************** 18 19 20 Average Coverage Criterion 21 22 Step 0 4147 23 Step 1 2073 24 Step 2 3110 25 Step 3 3628 26 Step 4 3369 27 Step 5 3498 28 Step 6 3433 29 Step 7 3465 30 Step 8 3449 31 Step 9 3441 32 Step 10 3437 33 Step 11 3435 34 Step 12 3434 35 36 37 ACC sample size = 3435 38 39 For a sample size of 3435, 40 the minimum HPD coverage of fixed length 0.40000D-01 41 over all possible data values is 0.98104004914D+00. 42 This interval is (0.47992732919D+00, 0.51992732919D+00). 43 44 For a sample size of 3435, 45 the maximum length 46 required for 0.99000000000D+00 equal-tailed (non HPD) 47 coverage is 0.43909109763D-01. 48 This interval is (0.47797281841D+00, 0.52188192817D+00). 49 50 For a sample size of 3435, 51 the maximum HPD length 52 required for 0.99000000000D+00 HPD 53 coverage is 0.43909109763D-01. 54 This interval is (0.47797279023D+00, 0.52188189999D+00). 55 56 The event of 1717 successes 57 in 3435 trials occurs with 58 prior probability of 0.38573669212D-03. 59 60 *************************************************************************** 61 62 63 Average Length Criterion 64 65 Step 0 3435 66 Step 1 1717 67 Step 2 2576 68 Step 3 3005 69 Step 4 3220 70 Step 5 3112 71 Step 6 3058 72 Step 7 3031 73 Step 8 3044 74 Step 9 3037 75 Step 10 3040 76 Step 11 3042 77 Step 12 3041 78 79 80 ALC sample size = 3041 81 82 For a sample size of 3041, 83 the maximum HPD length of fixed coverage 0.99000D+00 84 over all possible data values is 0.46661439137D-01. 85 This interval is (0.47658723049D+00, 0.52324866962D+00). 86 87 For a sample size of 3041, 88 the minimum HPD coverage 89 required for 0.40000000000D-01 HPD 90 length is 0.97274051750D+00. 91 This interval is (0.47991791841D+00, 0.51991791841D+00). 92 93 The event of 1520 successes 94 in 3041 trials occurs with 95 prior probability of 0.43567886022D-03. 96 97 *************************************************************************** 98 99 100 Worst Outcome Criterion 101 102 Step 0 4147 103 Step 1 2073 104 Step 2 3110 105 Step 3 3628 106 Step 4 3887 107 Step 5 4017 108 Step 6 4082 109 Step 7 4114 110 Step 8 4130 111 Step 9 4138 112 Step 10 4142 113 Step 11 4140 114 Step 12 4141 115 116 117 WOC sample size = 4141 118 119 For a sample size of 4141, 120 the minimum HPD coverage of fixed length 0.40000D-01 121 over all possible data values is 0.99000456637D+00. 122 This interval is (0.47993971436D+00, 0.51993971436D+00). 123 124 For a sample size of 4141, 125 the maximum length 126 required for 0.99000000000D+00 equal-tailed (non HPD) 127 coverage is 0.39997549597D-01. 128 This interval is (0.47994095896D+00, 0.51993850855D+00). 129 130 For a sample size of 4141, 131 the maximum HPD length 132 required for 0.99000000000D+00 HPD 133 coverage is 0.39997549597D-01. 134 This interval is (0.47994093956D+00, 0.51993848915D+00). 135 136 The event of 2070 successes 137 in 4141 trials occurs with 138 prior probability of 0.32000608811D-03. 139 140 *************************************************************************** 141 142 143 Elapsed time 1301.940 seconds 144 Output description
In this section, we describe in detail the output
file b-2-1.5-.99-.04-.5 produced above.
Lines 102 to 139: (Worst Outcome Criterion)
Copyright
© Copyright Lawrence Joseph, David Wolfson and Roxane du Berger 1994
bhpd1 is a program written by Lawrence Joseph, David Wolfson and Roxane du Berger at the Division of Clinical Epidemiology, Department of Medicine, Montreal General Hospital. This program is an implementation of the manuscript Sample size calculations for binomial proportions via highest posterior density intervals by Lawrence Joseph, David B. Wolfson and Roxane du Berger. The Statistician, 1995;44(2):143-154 You are free to use this program, for non-commercial purposes only, under two conditions:
Please e-mail questions, comments, responses or requests to lawrence.joseph@mcgill.ca or write to Lawrence Joseph Division of Clinical Epidemiology McGill University Health Centre Royal Victoria Hospital 687 Pine Avenue West V Building, Room V2.10 Montreal, Quebec Canada, H3A 1A1 |