*! 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')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