*! version 1.1  20/4/97
pro def ssizebi
version 5.0
parse "`*'", parse(",")
  local pr "`1'"
  local npr: word count `pr'
  if `npr'<2 {
  	di in red "At least two event probabilities required"
  	exit 198
  }
  macro shift
  #delimit ;
  local options "NGroup(integer 2) ARatio(string) COndit noLOcal
                 ALpha(real 0.05) POwer(real 0.8) n(integer 0) TRend
                 DOses(string)";
  #delimit cr
  parse "`*'"
  if max(`alpha',1-`alpha')>=1 { 
  	di in red "alpha() out of range"
  	exit 198
  }
  if max(`power',1-`power')>=1 {
  	di in red "power() out of range"
  	exit 198
  }
  if `n'<0 { 
  	di in red "n() out of range"
  	exit 198
  }
  
  preserve
  if `n'==0 {
	local tit1 " Sample size:"
	local ssize "1"
  }
  else {
  	local tit1 " Power:      "
	local ssize "0"
  }
  drop _all
  local ngroup=max(`ngroup',`npr')
  qui set obs `ngroup'
  tempname PI pibar
  qui gen double `PI'=.
  parse "`pr'",parse(" ")
  local altp "`1'"
  local i "1"
  while `i'<=`npr' {
  	confirm number `1'
  	if max(`1',1-`1')>=1 { 
  		di in red "Event probabilities out of range"
  		exit 198
  	}
  	qui replace `PI'=`1' in `i'
  	local p=round(`1',.01)
	if `i'>1 {
	  	local altp "`altp',`p'"
	}
	macro shift
  	local i=`i'+1
  }
  summ `PI',meanonly
  scalar `pibar'=_result(3)
  if _result(6)<=_result(5) {
  	di in r /*
  	*/ "At least two distinct alternative event probabilities required"
  	exit 198
  }
  /*
    If the number of given alternative probabilities is less than the number of
    treatment groups, assume each of the remaining equals the mean of the given
    event probabilities. This minimises the non-centrality parameter and thus
    gives  conservative estimates of power and sample size.
  */
  if `npr'<`ngroup' {
  	local i=`npr'+1
  	qui replace `PI'=scalar(`pibar') in `i'/l
  	local p=round(scalar(`pibar'),.01)
  	while `i'<=`ngroup' {
  		local altp "`altp',`p'"
  		local i=`i'+1
  	}
  }
  tempname AR sar
  qui gen double `AR'=.
  if "`aratio'"=="" {
  	qui replace `AR'=1/`ngroup'
  	local allocr "Equal group sizes"
  }
  else {
  	scalar `sar'=0
  	parse "`aratio'",parse(" ")
	local allocr "`1'"
  	local i "1"
  	while `i'<=_N {
  		if "`1'"~="" {
  			confirm number `1'
  			if `1'<=0 {
			  	di in r  "Allocation ratio <=0 not alllowed"
			  	exit 198
			}
  			qui replace `AR'=`1' in `i'
  		}
  		else {
  			qui replace `AR'[`i']=`AR'[`i'-1]
  		}
  		scalar `sar'=scalar(`sar')+`AR'[`i']
  		local p=round(`AR'[`i'],.01)
		if `i'>1 {
	  		local allocr "`allocr':`p'"
		}
  		local i=`i'+1
  		macro shift
  	}
  	qui replace `AR'=`AR'/scalar(`sar')
  }
  summ `PI' [w=`AR'],meanonly
  scalar `pibar'=_result(3)

  tempname DOSE
  if "`trend'"~=""|"`doses'"~="" {
	local trtest "Linear trend test: doses = "
  	qui gen double `DOSE'=.
	if "`doses'"=="" {
		qui replace `DOSE'=_n-1
		local doses "1,...,`ngroup'"
	}
	else {
	  	parse "`doses'",parse(" ")
		local score "`1'"
	  	local i "1"
  		while `i'<=_N {
  			if "`1'"~="" {
  				confirm number `1'
	  			if `1'<0 {
				  	di in r  "Dose < 0 not alllowed"
				  	exit 198
				}
	  			qui replace `DOSE'=`1' in `i'
  			}
  			else {
  				qui replace `DOSE'=`DOSE'[`i'-1] in `i'
  			}
  			local p=round(`DOSE'[`i'],.01)
			if `i'>1 {
	  			local score "`score',`p'"
			}
	  		local i=`i'+1
  			macro shift
  		}
		local doses "`score'"
  	}
	summ `DOSE' [w=`AR'],meanonly
	qui replace `DOSE'=`DOSE'-_result(3)
  }

  tempname b MU q0 a D
  local K=`ngroup'-1
  scalar `b'=1-`power'
  if "`condit'"=="" {
	local test0 "Unconditional"
  	local tit2 "Unconditional comparison of `ngroup' binomial proportions"
  	qui gen double `MU'=`PI'-scalar(`pibar')
  	tempname s
  	scalar `s'=scalar(`pibar')*(1-scalar(`pibar'))
	if "`trtest'"=="" {
		local test1 "Chisquare test"
		_sp `MU' `MU' `AR', out(`q0')
		scalar `q0'=scalar(`q0')/scalar(`s')
		scalar `a'=invchi(`K',`alpha')
		if "`local'"=="" {
			if `n'==0 {
			  local n=npnchi(`K',scalar(`a'),`b')/scalar(`q0')
			  local D=`n'*`pibar'
			}
			else {
			  scalar `b'=nchi(`K',`n'*scalar(`q0'),scalar(`a'))
			}
		}
		else {
		  	tempname S sbar W  a0 a1 q1 eta g psi l
		  	qui gen double `S'=`PI'*(1-`PI')
  			summ `S' [w=`AR'],meanonly
		  	scalar `sbar'=_result(3)
			summ `S',meanonly
			scalar `a0'=(_result(18)-scalar(`sbar'))/scalar(`s')
			_sp `MU' `MU' `S' `AR', out(`q1')
			scalar `q1'=scalar(`q1')/scalar(`s')^2
			qui gen double `W'=1-2*`AR'
			_sp `S' `S' `W', out(`a1')
			scalar `a1'=(scalar(`a1')+scalar(`sbar')^2)/scalar(`s')^2
			if `n'==0 {
			  * *******************************
			  * Solve for n iteratvely
			  tempname n0 nl nu b0 sm
			  scalar `sm'=0.001
			  local i "1"
			  scalar `n0'=npnchi(`K',scalar(`a'),`b')/scalar(`q0')
			  _pe2 `a0' `q0' `a1' `q1' `K' `n0' `a' `b0'
			  if abs(scalar(`b0')-scalar(`b'))<=scalar(`sm') {
			    local i "0"
			  }
			  else {
			    if scalar(`b0')<scalar(`b') {
			      scalar `nu'=scalar(`n0')
			      scalar `nl'=scalar(`n0')/2
			    }
			    else {
			      scalar `nl'=scalar(`n0')
			      scalar `nu'=2*scalar(`n0')
			    }
			  }
			  while `i' {
			    scalar `n0'=(scalar(`nl')+scalar(`nu'))/2
			    _pe2 `a0' `q0' `a1' `q1' `K' `n0' `a' `b0'
			    if abs(scalar(`b0')-scalar(`b'))<=scalar(`sm') {
			       local i "0"
			    }
			    else {
			      if scalar(`b0')<scalar(`b') {
			        scalar `nu'=scalar(`n0')
			      }
			      else {
			        scalar `nl'=scalar(`n0')
			      }
			      local i=`i'*((scalar(`nu')-scalar(`nl'))>1)
			    }
			  }
			  local n=scalar(`n0')
			  * *******************************
			}
			else {
			  _pe2 `a0' `q0' `a1' `q1' `K' `n' `a' `b'
			}
		}
	}
	else {
		local test1 "`trtest'`doses'"
		tempname tr q1
		_sp `MU' `DOSE' `AR', out(`tr')
		_sp `DOSE' `DOSE' `AR', out(`q0')
		scalar `q0'=scalar(`q0')*scalar(`s')
		if "`local'"=="" {
			scalar `q1'=scalar(`q0')
		}
		else {
		  	tempname S W 
		  	qui gen double `S'=`PI'*(1-`PI')
			_sp `DOSE' `DOSE' `S' `AR', out(`q1')
		}
		scalar `a'=sqrt(scalar(`q0'))*invnorm(1-`alpha'/2)
		if `n'==0 {
		  scalar `a'=scalar(`a')+sqrt(scalar(`q1'))*invnorm(`power')
		  local n=(scalar(`a')/scalar(`tr'))^2
		}
		else {
		  scalar `a'=abs(scalar(`tr'))*sqrt(`n')-scalar(`a')
		  scalar `b'=1-normprob(scalar(`a')/sqrt(scalar(`q1')))
		}
	}
	local D=`n'*scalar(`pibar')
  }
  else {
	local test0 "Conditional"
  	local tit2 "Conditional test using Peto's approx to the odds ratio"
  	tempname LOR l d v
	scalar `v'=scalar(`pibar')*(1-scalar(`pibar'))
  	qui gen double `LOR'=log(`PI')-log(1-`PI')-log(`PI'[1])+log(1-`PI'[1])
  	qui replace `LOR'=0 in 1
	summ `LOR' [w=`AR'],meanonly
	qui replace `LOR'=`LOR'-_result(3)
	if "`trtest'"=="" {
		local test1 "Chisquare test"
		_sp `LOR' `LOR' `AR', out(`q0')
		scalar `a'=invchi(`K',`alpha')
		if `n'==0 {
		  scalar `l'=npnchi(`K',scalar(`a'),`b')
		  scalar `d'=scalar(`l')
		  scalar `l'=sqrt(scalar(`l')*(scalar(`l')-4*scalar(`q0')*scalar(`v')))
		  scalar `d'=(scalar(`d')+scalar(`l'))/(2*scalar(`q0')*(1-scalar(`pibar')))
		  local n=scalar(`d')/scalar(`pibar')
		}
		else {
		  scalar `d'=`n'*scalar(`pibar')
		  scalar `l'=scalar(`d')*(`n'-scalar(`d'))*scalar(`q0')/(`n'-1)
		  scalar `b'=nchi(`K',scalar(`l'),scalar(`a'))
		}
	}
	else {
  		local test1 "`trtest'`doses'"
		tempname tr
		_sp `DOSE' `LOR' `AR', out(`tr')
		_sp `DOSE' `DOSE' `AR', out(`q0')
		scalar `a'=invnorm(1-`alpha'/2)
		if `n'==0 {
		  scalar `a'=sqrt(scalar(`q0'))*(scalar(`a')+invnorm(`power'))
		  scalar `l'=(scalar(`a')/scalar(`tr'))^2
		  scalar `d'=scalar(`l')
		  scalar `l'=sqrt(scalar(`l')*(scalar(`l')-4*scalar(`v')))
		  scalar `d'=(scalar(`d')+scalar(`l'))/(2*(1-scalar(`pibar')))
		  local n=scalar(`d')/scalar(`pibar')
		}
		else {
		  scalar `d'=`n'*scalar(`pibar')
		  scalar `l'=scalar(`d')*(`n'-scalar(`d'))/(`n'-1)
		  scalar `a'=abs(scalar(`tr'))*sqrt(scalar(`l')/scalar(`q0'))-scalar(`a')
		  scalar `b'=1-normprob(scalar(`a'))
		}
	}
	local D=scalar(`d')
  }

  local n=round(`n',1)+(round(`n',1)<`n')
  local D=int(`D'+1)
  local power=1-`b'
    
  noi dis in gr _col(7) "`tit1' `ngroup'-groups comparison"
  noi dis in gr _col(3) "`tit2'" _n _dup(58) "-"
  if "`trtest'"~="" {
	noi dis in gr "`trtest'" _col(44) in ye "`doses'"
  }
  noi dis in gr "Anticipated event probabilities: " _col(44) in ye "`altp'"
  noi dis in gr "Allocation ratios:               " _col(44) in ye "`allocr'"

  noi dis in gr "Two-sided alpha = "              _col(44) %5.3f in ye `alpha'
  if `ssize'==1 {
  	noi dis in gr "Power = "                _col(44) %5.3f in ye `power' _n
  }
  noi dis in gr "Total sample size = " _col(44) %5.0f in ye `n' 
  noi dis in gr "Expected total number of events = "       _col(44) %5.0f in ye `D'
  if `ssize'==0 {
  	noi dis in gr _n "Power = "             _col(44) %5.2f in ye `power'
  }
  
  global S_1 "`test0' `test1'"
  global S_2 "`altp'"
  global S_3 "`allocr'"
  global S_4=`alpha'
  global S_5=`power'
  global S_6=`n'
  global S_7=`D'
  end
* *****************************************************************************

prog def _sp
* Calculate sum of products.
	local varlist "req ex min(1)"
	local options "OUt(string)"
	parse "`*'"
	parse "`varlist'",parse(" ")
	tempname SP
	qui gen double `SP'=`1'
	macro shift
	while "`1'"~="" {
		qui replace `SP'=`SP'*`1'
		macro shift
	}
	summ `SP',meanonly
	scalar `out'=_result(18)
end
* *****************************************************************************

prog def _pe2
* Calculate beta=P(type II error)
	local a0 "`1'"
	local q0 "`2'"
	local a1 "`3'"
	local q1 "`4'"
	local k  "`5'"
	local n  "`6'"
	local a  "`7'"
	local b  "`8'"

	tempname b0 b1 f l
	scalar `b0'=scalar(`a0')+scalar(`n')*scalar(`q0')
	scalar `b1'=scalar(`a1')+2*scalar(`n')*scalar(`q1')
	scalar `l'=scalar(`b0')^2-`k'*scalar(`b1')
	scalar `f'=sqrt(scalar(`l')*(scalar(`l')+`k'*scalar(`b1')))
	scalar `l'=(scalar(`l')+scalar(`f'))/scalar(`b1')
	scalar `f'=scalar(`a')*(`k'+scalar(`l'))/scalar(`b0')
	scalar `b'=nchi(`k',`l',`f')
end