* Biostat 200 Lab 7

cd "~/work/mb/atcr/200/labs"

capture log close
log using "lab 7", replace

/* In this lab, we will show how to use various power commands in Stata to calculate sample size, power, and minimum detectable effects for use in study planning.  For a comprehensive introduction to this subject, read the intro section of the pdf help for PSS, the suite of Stata power and sample size commands.  In general these functions calculate sample size, power, or minimum detectable effects (MDEs).  By specifying two of these parameters, we get an estimate of the third.

********************************************************************************
One-sample testing against a fixed null hypothesis.  In lab 5, we tested the null hypothesis that mean BMI among men with diabetes does not differ from the population mean among men without diabetes, assumed to be 24 kg/m^2. Now suppose the population mean (SD) among men with diabetes is 25 (2.7).  How large a sample would we need to provide 90% power to reject the null hypothesis that the mean is 24? To do this, we have to specify the MDE (via the two group means) and power. */

power onemean 24 25, sd(2.7) power(.9)

/* Now suppose you have a limited budget and can only afford 50 particpants.  What would your power be? Here we specify the sample size rather than power. */

power onemean 24 25, sd(2.7) n(50)

/* So a little lower than the conventional level of 80%.  It might be reasonable to use a one-sided test in this context, although it can be controversial. */

power onemean 24 25, sd(2.7) n(50) onesided

/* That would fix the problem.  

Stata can also estimate the MDE, fixing the sample size at 50 and power at 80%. */

power onemean 24, sd(2.7) n(50) power(.8)

/* So this indicates the population mean BMI among men with diabetes would have to at least 25.1 for a sample of size 50 to provide 80% power to detect the difference. 

Stata can do any of this for a range of inputs, and output the results as both a table and a graph.  First do this to estimate sample size. Note the syntax to specify numeric lists: Stata automatically expands 25/27 as 25 26 27 28, and likewise interprets .8(.05).95 as .8 .85 .9 .95. */

power onemean 24 (25/28), sd(2.7) power(.8(.05).95) ///
	table graph(name(plot1, replace))

/* This plot shows that as the difference in mean increases, the sample size decreases and increments in power get cheaper.  The parameter delta is the between-group difference in means, standardized by the common SD. 

Now calculate MDEs, varying the sample size and SD. */

power onemean 24, sd(2.5(0.25)3.25) n(50(10)100) power(.8) ///
	table graph(name(plot2, replace))

/* This plot shows that we can detect smaller between-group differences as the sample size increases and the SD decreases.  

Stata can do similar calculations for binary outcomes. Suppose the population prevalence of hypertension (HTN) is 40% and 25% among men with and without diabetes, respectively.  How large a sample of men with diabetes would we need to provide 80% power to reject a null prevalence of 25%? */

power oneproportion .25 .4, power(.8)

* Now calculate for a sample of 50 men with diabetes

power oneproportion .25 .4, n(50)

* Last, get MDEs for a range of sample sizes and power

power oneproportion .25, ///
	n(50(10)100) power(.8(.05).95) table graph(name(plot3, replace))

/* So MDEs decrease with sample size but increase with power.

********************************************************************************
Two-sample testing. 

Now suppose we need to estimate mean BMI in men without diabetes,rather than assuming we know the population mean for this group. You'll see that the pricetag for estimating the control mean is considerable.

By default, for the two-sample calculations the SD is assumed to be the same in both groups, and specified using a single sd() option.  But the command does allow us to specify different SDs. To see how this works, we'll assume that the SDs are 2.4 and 2.7 in men without and with diabetes respectively. */

power twomeans 24 25, sd1(2.4) sd2(2.7) power(.9)

/* Now suppose it is twice as hard to recruit men with diabetes, so you can recruit at most 75 men with diabetes and 150 controls. What would your power be? */

power twomeans 24 25, sd1(2.4) sd2(2.7) n1(150) n2(75)

/* So a little lower than the conventional level of 80%.  It might be reasonable to use a one-sided test in this context. */

power twomeans 24 25, sd1(2.4) sd2(2.7) n1(150) n2(75) onesided

/* We can estimate the MDE with a sample of 225, fixing power at 80%. */

power twomeans 24, sd1(2.4) sd2(2.7) n1(150) n2(75) power(.8)

/* With a continuous outcome, the mean in the reference group doesn't affect the MDE. */

power twomeans 0, sd1(2.4) sd2(2.7) n1(150) n2(75) power(.8)

/* So both calculations indicate that population mean BMI in the two groups would need to differ by about 1.04 kg/m^2 for a sample of size 225 to provide 80% power. 

Now calculate sample size for a range of means among men with diabetes and power from 80 to 95%.  */

power twomeans 24 (25/28), sd1(2.4) sd2(2.7) power(.8(.05).95) ///
	table graph(name(plot4, replace))

/* Now calculate MDEs, varying the sample size and SDs, fixing the ratio of the sample sizes at 2:1. */

power twomeans 24, sd1(2.2 2.6) sd2(2.5 2.9) n(250(25)400) nratio(0.5) power(.8) ///
	table graph(name(plot5, replace))

/* Pretty much the same story as in the one-sample case, except that now the overall sample size is almost 4 times larger!

Now we will do this for binary outcomes, again assuming hypertension prevalence of 25% and 40% among men without and with diabetes. Again, suppose men with diabetes are twice as hard to recruit, so the ratio of the sample sizes is 0.5. */

power twoproportions .25 .4, power(.8) nratio(0.5)

/* Note that in contrast to continuous outcomes, for binary outcomes the reference proportion does matter, because it drives the variance. */

power twoproportions .2 .35, power(.8) nratio(0.5)
power twoproportions .3 .45, power(.8) nratio(0.5)

/* Now calculate power for sample of 75/150 men with/without diabetes. */

power twoproportions .25 .4, n1(150) n2(75)

/* Last, get MDEs for a range of overall sample sizes, fixing the reference mean at 25%, the ratio of men with and without diabetes at 50%, and varying power from 80 to 95%. */

power twoproportions .25, n(200(25)300) nratio(0.5) power(.8(.05).95) ///
	table graph(name(plot6, replace))
	
/* What does this plot indicate about the effects of sample size and power on the MDE? */

log close