*! sampsi2 version 1.0.1  6/8/1998  by PTS (p.seed@umds.ac.uk)  STB-45 sbe18.1
*! (version 1.0.0 STB-40 sbe18)
*! Adapted version of sampsi to handle repeated measured
*! following Frison & Pocock (1992) Stats in Medicine
*!
*! usage: sampsi2 mean1 mean2, pre(#) post(#) r0(#) r1(#) r01(#) 
*! 		method (post|change|ancova|all) 
*!  		alpha(#) power(#) n1(#) n2(#) ratio(#) sd1(#) sd2(#)
*!
*! Note: continuous outcomes only with repeated measures.

prog define sampsi2

	version 5.0
   	parse "`*'", parse(" ,")
   	local m1 = `1'
   	local m2 = `2'
	confirm number `m1'
	confirm number `m2'
	mac shift 2
	#delimit ;
	local options "pre(real 0) post(real 1) r0(real -1) r1(real -1) 
		r01(real -1) Method(str) SD1(real -1) SD2(real -1) *";
	#delimit cr
	parse "`*'"

	if "`method'" ~= "" & `sd1' == -1 {
 		di in red "method " in ye "`method' " in red "cannot be used without " in ye "sd1"
		exit 198 }
	if (`pre' ~= 0 |`post' ~= 1 ) & `sd1' == -1 {
		di in red "options " in ye "pre " in red "and " in ye "post " in red "cannot be used without " in ye "sd1"
		exit 198 }
	if "`method'" ~= "" & `sd1' == -1 {
		di in red "method " in ye "`method' " in red "cannot be used without " in ye "sd1"
		exit 198 }


	if `r1' ~= -1 & `r0' == -1 {local r0 = `r1' }
	if `r1' ~= -1 & `r01' == -1 {local r01 = `r1' }

	confirm integer number `pre' 
	cap confirm integer number `pre' 
	if _rc { di in red "pre must be a count" 
		  exit 198 }
	cap assert `pre' >= 0
	if _rc { di in red "pre must be a count" 
		  exit 198 }
	cap confirm integer number `post' 
	if _rc { di in red "post must be a count" 
		  exit 198 }
	cap assert `post' >= 0
	if _rc { di in red "post must be a count" 
		  exit 198 }

	if `post' > 1 & `r1' == -1 { 
		di in red "Correlation r1 is needed with more than one post-trial (follow-up) measurement"
		exit 198 }	
	if `pre' ~= 0 & `r01' == -1 { 
		di in red "Correlation r01 is needed with pre-trial (baseline) measurements"
		exit 198 }	
	if `pre' > 1 & `r0' == -1 { 
		di in red "Correlation r0 is needed with more than one pre-trial (baseline) measurement"
		exit 198 }	

	if `r1' >= 1 | `r1' <-1 {
		di in red "Correlation r1 must be between -1 and 1" 
		exit 198 }	
	if `r0' >= 1 | `r0' <-1 {
		di in red "Correlation r0 must be between -1 and 1" 
		exit 198 }	
	if `r01' >= 1 | `r01' <-1 {
		di in red "Correlation r01 must be between -1 and 1" 
		exit 198 }	

	cap assert "`method'" == "" | "`method'" == "post" | "`method'" == "change" | "`method'" == "ancova" | "`method'" == "all"
	if _rc { di in red "method invalid"
		exit 198 }

	if "`method'" == "change" & `pre' == 0 { di in red "Method CHANGE can only be used with baseline measurements"
		exit 198 }
	if "`method'" == "ancova" & `pre' == 0 { di in red "Method ANCOVA can only be used with baseline measurements"
		exit 198 }

	if (`pre' ~= 0 | `post' ~= 1)  & "`method'" == "" { local method "all" }


	if `sd1' == -1 { local sd1 }
	if `sd2' == -1 { local sd2 }

	
/* HEADING */

	if "`method'" == "" {
		if "`sd1'" ~= "" {local sd1 "sd1(`sd1')" }
		if "`sd2'" ~= "" {local sd2 "sd2(`sd2')" }
		_sampsi `m1' `m2', `sd1'  `sd2' `options'
		global S_4 = 1
		exit
		}
	di in gr _n "Power calculations for repeated measures"
	di in gr "(Frison & Pocock 1992) " _n
	if `post' == 1 { di in ye `post' in gr " follow-up measurement" }
	if `post' > 1 { di in ye `post' in gr " follow-up measurements, correlation " in ye %6.3f `r1'}

	if `pre' == 1 {
		di in ye `pre' in gr " baseline measurement
		di in gr "Correlation between baseline & follow-up " in ye %6.3f `r01'
		}
	if `pre' > 1 {
		di in ye `pre' in gr " baseline measurements, correlation " in ye %6.3f `r0'
		di in gr "Correlation between baseline & follow-up " in ye %6.3f `r01'
		}
	if "`sd2'" ~= "" { 
		di in gr "Raw Standard deviations: " in ye %6.3f `sd1' in gr " and " in ye %6.3f `sd2j' in gr " adjustment " in ye %6.3f `sd_r'
		}
	if "`sd2'" == "" { 
		di in gr "Raw Standard deviation: " in ye %6.3f `sd1'
		}


/* POST */
	if "`method'" == "post" | "`method'" == "all" {	

		di _n in gr "Method: POST" 

		local sd_r = ((1 + (`post' - 1) * `r1')/`post')^.5
		local sd1_adj = `sd1'*`sd_r'
		if "`sd2'" ~= "" { 
			local sd2_adj = `sd2'*`sd_r' 
			di in gr "Relative efficiency: " in ye %6.3f `sd_r'^-2 in gr " Adjustment: " in ye %6.3f `sd_r' in gr " Adjusted SDs: " in ye %6.3f `sd1_adj' in gr " and " in ye %6.3f `sd2_adj'
			qui _sampsi `m1' `m2',  sd1(`sd1_adj') sd2(`sd2_adj') `options'
			}
		if "`sd2'" == "" { 
			di in gr "Relative efficiency: " in ye %6.3f `sd_r'^-2 in gr " Adjustment: " in ye %6.3f `sd_r' in gr " Adjusted SD: " in ye %6.3f `sd1_adj' 
			qui _sampsi `m1' `m2',  sd1(`sd1_adj') `options'
			}
		if index("`options'", "onesam") == 0 {di in gr "Sample sizes: n1 = " in ye $S_1 in gr ", n2 = " in ye $S_2}
		if index("`options'", "onesam") ~= 0 {di in gr "Sample size: n = " in ye $S_1}
		di in gr "significance = " in ye %6.3f 1 - $S_level/100 in gr ", power = " in ye %6.3f $S_3
		if "`method'" == "all" & `pre' == 0 { di in gr _n "Methods CHANGE and ANCOVA can only be used with baseline measurements"}
		global S_4 = `sd_r'
		}

/* CHANGE */

	if "`method'" == "change" | ("`method'" == "all" & `pre' >=1 ) {
		/* calculate sd ratio */ 
		di _n in gr "Method: CHANGE" 
		local sd_r = (((1 + (`post' - 1) * `r1')/`post') + ((1 + (`pre'-1) * `r0')/`pre') - 2*`r01')^.5

		local sd1_adj = `sd1'*`sd_r'
		if "`sd2'" ~= "" { 
			local sd2_adj = `sd2'*`sd_r' 
			di in gr "Relative efficiency: " in ye %6.3f `sd_r'^-2 in gr " Adjustment: " in ye %6.3f `sd_r' in gr " Adjusted SDs: " in ye %6.3f `sd1_adj' in gr " and " in ye %6.3f `sd2_adj'
			qui _sampsi `m1' `m2',  sd1(`sd1_adj') sd2(`sd2_adj') `options'
			}
		if "`sd2'" == "" { 
			di in gr "Relative efficiency: " in ye %6.3f `sd_r'^-2 in gr " Adjustment: " in ye %6.3f `sd_r' in gr " Adjusted SD: " in ye %6.3f `sd1_adj' 
			qui _sampsi `m1' `m2',  sd1(`sd1_adj') `options'
			}
		if index("`options'", "onesam") == 0 {di in gr "Sample sizes: n1 = " in ye $S_1 in gr ", n2 = " in ye $S_2}
		if index("`options'", "onesam") ~= 0 {di in gr "Sample size: n = " in ye $S_1}
		di in gr "significance = " in ye %6.3f 1 - $S_level/100 in gr ", power = " in ye %6.3f $S_3
		global S_4 = `sd_r'
		}

/* ANCOVA */
	if "`method'" == "ancova" | ("`method'" == "all" & `pre' >=1 ) {
		/* calculate sd1_adj */ 
		di _n in gr "Method: ANCOVA" 

		local sd_r = ((((1 + (`post' - 1) * `r1')/`post') -(`r01'^2)*`pre'/(1 +(`pre'-1)*`r0')))^.5
		local sd1_adj = `sd1'*`sd_r'

		if "`sd2'" ~= "" { 
			local sd2_adj = `sd2'*`sd_r' 
			di in gr "Relative efficiency: " in ye %6.3f `sd_r'^-2 in gr " Adjustment: " in ye %6.3f `sd_r' in gr " Adjusted SDs: " in ye %6.3f `sd1_adj' in gr " and " in ye %6.3f `sd2_adj'
			qui _sampsi `m1' `m2',  sd1(`sd1_adj') sd2(`sd2_adj') `options'
			}
		if "`sd2'" == "" { 
			di in gr "Relative efficiency: " in ye %6.3f `sd_r'^-2 in gr " Adjustment: " in ye %6.3f `sd_r' in gr " Adjusted SD: " in ye %6.3f `sd1_adj' 
			qui _sampsi `m1' `m2',  sd1(`sd1_adj') `options'
			}
		if index("`options'", "onesam") == 0 {di in gr "Sample sizes: n1 = " in ye $S_1 in gr ", n2 = " in ye $S_2}
		if index("`options'", "onesam") ~= 0 {di in gr "Sample size: n = " in ye $S_1}
		di in gr "significance = " in ye %6.3f 1 - $S_level/100 in gr ", power = " in ye %6.3f $S_3
		global S_4 = `sd_r'
		}

end sampsi2

program define _sampsi
	version 4.0
   	parse "`*'", parse(" ,")
   	local m1 = `1'
   	local m2 = `2'
	confirm number `m1'
	confirm number `m2'
	mac shift 2
	local dalpha = 1 - $S_level/100
	#delimit ;
	local options "Alpha(real `dalpha') Power(real 0.90) N1(int 0) 
		N2(int 0) SD1(real 0) SD2(real 0) Ratio(real 1) 
		ONESAMple ONESIDed" ;
	#delimit cr
	parse "`*'"


	if `alpha'<=0 | `alpha'>=1 { 
		di in red "alpha() out of range"
		exit 198
	}

	if `power'<=0 | `power'>=1 {
		di in red "power() out of range"
		exit 198
	}
	if `sd1'<0 | `sd2'<0 { 
		di in red "sd() out of range"
		exit 198
	}
	if `ratio'<=0 { 
		di in red "ratio() out of range"
		exit 198
	}
	if `n1'<0 | `n2'<0 { 
		di in red "n() out of range"
		exit 198
	}

	tempname diff n0 pbar r1 w w1 w2 za zb
	scalar `diff' = abs(`m1' - `m2')
	scalar `w1' = `m1'*(1 - `m1')
	scalar `w2' = `m2'*(1 - `m2')
	if `sd1' == 0 & `sd2' == 0 & (`w1' < 0 | `w2' < 0) {
		#delimit ;
		di in red "if testing means, then sd(#) must be specified" _n
		"if testing proportions, the proportions must be between 0 and 1" ;
		#delimit cr
		exit 499
	}	
	if `sd1' ~= 0 & `sd2' ~= 0 & "`onesamp'"~="" {
		di in red "for one-sample comparison of mean, only one sd(#) can be specified"
		exit 499
	}	
   	if "`oneside'"~="" { scalar `za' = invnorm(1 - `alpha') }
	else                { scalar `za' = invnorm(1 - `alpha'/2)   }
	if `n1' == 0 & `n2' == 0 {  /* compute sample size */
		scalar `zb' = invnorm(`power')
		if `sd1' ~= 0 | `sd2' ~= 0 { /* means */
			scalar `w' = ((`za'+`zb')/`diff')^2
			if `sd1' == 0 { local sd1 = `sd2' }
			if "`onesamp'"~="" {  /* one sample */
				local n1 = `w'*`sd1'^2
				if `n1' ~= int(`n1') { local n1 = int(`n1' + 1) }
	#delimit ;
	di _n in gr "Estimated sample size for one-sample comparison of mean" _n
	"  to hypothesized value" _n(2)
	"Test Ho: m = " in ye %6.0g `m1' in gr ", where m is the mean in the
population"
	_n(2) "Assumptions:" _n(2) _col(10) "alpha = " in ye %8.4f `alpha' _c ;
	if "`oneside'"~="" { di in gr "  (one-sided)" } ;
	else                { di in gr "  (two-sided)" } ;
	di _col(10) in gr "power = " in ye %8.4f `power' _n in gr " alternative m = "
	in ye %8.0g `m2' _n in gr _col(13) "sd = " in ye %8.0g `sd1' _n(2) in gr
	"Estimated required sample size:" _n(2) _col(14) "n = " in ye %8.0f `n1' ;
	#delimit cr
			}
			else {  /* two sample */
				if `sd2' == 0 { local sd2 = `sd1' }
				local n1 = `w'*(`sd1'^2+`sd2'^2/`ratio')
				if `n1' ~= int(`n1') { local n1 = int(`n1' + 1) }
				local n2 = `ratio'*`n1'
				if `n2' ~= int(`n2') { local n2 = int(`n2' + 1) }
	#delimit ;
	di _n in gr "Estimated sample size for two-sample comparison of means" _n(2)
	"Test Ho: m1 = m2, where m1 is the mean in population 1" _n _col(21)
	"and m2 is the mean in population 2" _n "Assumptions:" _n(2) _col(10)
"alpha = " in ye
	%8.4f `alpha' _c ;
	if "`oneside'"~="" { di in gr "  (one-sided)" } ;
	else                { di in gr "  (two-sided)" } ;
	di _col(10) in gr "power = " in ye %8.4f `power' _n in gr _col(13) "m1 = "
	in ye %8.0g `m1' _n in gr _col(13) "m2 = " in ye %8.0g `m2' _n in gr
_col(12) "sd1 = "
	in ye %8.0g `sd1' _n in gr _col(12) "sd2 = " in ye %8.0g `sd2' _n
	in gr _col(10) "n2/n1 = " in ye %8.2f `ratio' _n(2)
	in gr "Estimated required sample sizes:" _n(2) _col(13) "n1 = " in ye %8.0f
`n1' _n
	in gr _col(13) "n2 = " in ye %8.0f `n2' ; 
	#delimit cr
			}
		}
		else { /* proportions */
			if "`onesamp'"~="" {  /* one sample */
				local n1 = ((`za'*sqrt(`w1')+`zb'*sqrt(`w2'))/`diff')^2
				if `n1' ~= int(`n1') { local n1 = int(`n1' + 1) }
	#delimit ;
	di _n in gr "Estimated sample size for one-sample comparison of proportion" _n
	"  to hypothesized value" _n(2)
	"Test Ho: p = " in ye %6.4f `m1' in gr ", where p is the proportion in the
population" _n(2)
	"Assumptions:" _n(2) _col(10) "alpha = " in ye %8.4f `alpha' _c ;
	if "`oneside'"~="" { di in gr "  (one-sided)" } ;
	else                { di in gr "  (two-sided)" } ;
	di _col(10) in gr "power = " in ye %8.4f `power' _n in gr " alternative p =
" in ye %8.4f `m2' 
	_n(2) in gr "Estimated required sample size:" _n(2) _col(14) "n = " in ye
%8.0f `n1' ;
	#delimit cr
			}
			else {  /* two sample */
				scalar `r1' = `ratio' + 1
				scalar `pbar' = (`m1'+`ratio'*`m2')/`r1'
				#delimit ;
				scalar `n0' = ( `za'*sqrt(`r1'*`pbar'*(1-`pbar'))
					       +`zb'*sqrt(`ratio'*`w1'+`w2'))^2
					      /(`ratio'*`diff'^2) ;
				#delimit cr
				local n1 = (`n0'/4)*(1+sqrt(1+2*`r1'/(`n0'*`ratio'*`diff')))^2 
				if `n1' ~= int(`n1') { local n1 = int(`n1' + 1) }
				local n2 = `ratio'*`n1'
				if `n2' ~= int(`n2') { local n2 = int(`n2' + 1) }
	#delimit ;
	di _n in gr "Estimated sample size for two-sample comparison of
proportions" _n(2)
	"Test Ho: p1 = p2, where p1 is the proportion in population 1" _n
	_col(21) "and p2 is the proportion in population 2" _n
	"Assumptions:" _n(2) _col(10) "alpha = " in ye %8.4f `alpha' _c ;
	if "`oneside'"~="" { di in gr "  (one-sided)" } ;
	else                { di in gr "  (two-sided)" } ;
	di _col(10) in gr "power = " in ye %8.4f `power' _n in gr _col(13)
	"p1 = " in ye %8.4f `m1' _n in gr _col(13)
	"p2 = " in ye %8.4f `m2' _n in gr _col(10) "n2/n1 = " in ye %6.2f `ratio' _n(2)
	in gr "Estimated required sample sizes:" _n(2) _col(13) "n1 = " in ye %8.0f
`n1' _n
	in gr _col(13) "n2 = " in ye %8.0f `n2' ;
	#delimit cr
			}
		}
	}
/* Compute power. */
	else {
		if `n2' == 0 {  /* determine n2 from n1 and ratio */
			local n2 = `ratio'*`n1'
			if `n2' ~= int(`n2') { local n2 = int(`n2' + 1) }
		}
		else if `n1' == 0 {  /* determine n1 from n2 and ratio */
			local n1 = `n2'/`ratio'
			if `n1' ~= int(`n1') { local n1 = int(`n1' + 1) }
		}
		if `sd1' ~= 0 | `sd2' ~= 0 { /* means */
			if `sd1' == 0 { local sd1 = `sd2' }
			if "`onesamp'"~="" {  /* one sample */
				local power = normprob(`diff'*sqrt(`n1')/`sd1' - `za')
	#delimit ;
	di _n in gr "Estimated power for one-sample comparison of mean" _n
	"  to hypothesized value" _n(2)
	"Test Ho: m = " in ye %6.0g `m1' in gr ", where m is the mean in the
population"
	_n(2) "Assumptions:" _n(2) _col(10) "alpha = " in ye %8.4f `alpha' _c ;
	if "`oneside'"~="" { di in gr "  (one-sided)" } ;
	else                { di in gr "  (two-sided)" } ;
	di in gr " alternative m = " in ye %8.0g `m2' _n in gr _col(13) "sd = "
	in ye %8.0g `sd1' _n in gr " sample size n = " in ye %8.0f `n1' _n(2) 
	in gr "Estimated power:" _n(2) _col(10) "power = " in ye %8.4f `power' ;
	#delimit cr
			}
			else {  /* two sample */
				if `sd2' == 0 { local sd2 = `sd1' }
				scalar `w' = sqrt(`sd1'^2/`n1' + `sd2'^2/`n2')
				local power = normprob(`diff'/`w' - `za')
	#delimit ;
	di _n in gr "Estimated power for two-sample comparison of means" _n(2)
	"Test Ho: m1 = m2, where m1 is the mean in population 1" _n _col(21)
	"and m2 is the mean in population 2" _n "Assumptions:" _n(2) _col(10)
	"alpha = " in ye %8.4f `alpha' _c ;
	if "`oneside'"~="" { di in gr "  (one-sided)" } ;
	else                { di in gr "  (two-sided)" } ;
	di in gr _col(13) "m1 = " in ye %8.0g `m1' _n in gr _col(13) "m2 = "
	in ye %8.0g `m2' _n in gr _col(12) "sd1 = " in ye %8.0g `sd1' _n
	in gr _col(12) "sd2 = " in ye %8.0g `sd2' _n in gr "sample size n1 = " 
	in ye %8.0f `n1' _n in gr _col(13) "n2 = " in ye %8.0f `n2' _n
	in gr _col(10) "n2/n1 = " in ye %8.2f `n2'/`n1' _n(2)
	in gr "Estimated power:" _n(2) _col(10) "power = " in ye %8.4f `power' ;
	#delimit cr
			}
		}
		else { /* proportions */
			if "`onesamp'"~="" {  /* one sample */
				local power = normprob((`diff'*sqrt(`n1') - `za'*sqrt(`w1'))/sqrt(`w2'))
	#delimit ;
	di _n in gr "Estimated power for one-sample comparison of proportion" _n
	"  to hypothesized value" _n(2)
	"Test Ho: p = " in ye %6.4f `m1' in gr ", where p is the proportion in the
population"
	_n(2) "Assumptions:" _n(2) _col(10) "alpha = " in ye %8.4f `alpha' _c ;
	if "`oneside'"~="" { di in gr "  (one-sided)" } ;
	else                { di in gr "  (two-sided)" } ;
	di in gr " alternative p = " in ye %8.4f `m2' _n in gr " sample size n = "
in ye %8.0f
	`n1' _n(2) in gr "Estimated power:" _n(2) _col(10) "power = " in ye %8.4f
`power' ;
	#delimit cr
			}
			else {  /* two sample */
				local ratio = `n2'/`n1'
				scalar `r1' = `ratio' + 1
				scalar `pbar' = (`m1'+`ratio'*`m2')/`r1'
				scalar `n0' = (`n1' - `r1'/(2*`ratio'*`diff'))^2/`n1'
				#delimit ;
				scalar `zb' = (`diff'*sqrt(`ratio'*`n0') -
`za'*sqrt(`r1'*`pbar'*(1-`pbar')))
					      /sqrt(`ratio'*`w1'+`w2') ;
				local power = normprob(`zb') ;
	di _n in gr "Estimated power for two-sample comparison of proportions" _n(2)
	"Test Ho: p1 = p2, where p1 is the proportion in population 1" _n
	_col(21) "and p2 is the proportion in population 2" _n
	"Assumptions:" _n(2) _col(10) "alpha = " in ye %8.4f `alpha' _c ;
	if "`oneside'"~="" { di in gr "  (one-sided)" } ;
	else                { di in gr "  (two-sided)" } ;
	di in gr _col(13) "p1 = " in ye %8.4f `m1' _n in gr _col(13)
	"p2 = " in ye %8.4f `m2' _n in gr "sample size n1 = " 
	in ye %8.0f `n1' _n in gr _col(13) "n2 = " in ye %8.0f `n2' _n
	in gr _col(10) "n2/n1 = " in ye %8.2f `n2'/`n1' _n(2)
	in gr "Estimated power:" _n(2) _col(10) "power = " in ye %8.4f `power' ;
	#delimit cr
			}
		}
	}
	global S_1 = `n1'
	global S_2 = `n2'
	global S_3 = `power'
end