*! Date        : 25 Oct 2005
*! Version     : 1.00
*! Authors     : Adrian Mander
*! Email       : adrian.mander@mrc-hnr.cam.ac.uk
*!
*! Sample size calculations for matched case/control study

prog def sampsi_mcc, rclass
version 9.0
syntax [varlist] [, P0(real 0.5) ALT(real 1.5) PHI(real 0.2) M(integer 1) N1(integer 100) Alpha(real 0.05) POWER(real 0.9) ////
Solve(string) ]

/*
 Null is the probability of Control being exposed
 Alt is the odds ratio alternative
 phi is the correlation between case and control exposure status

 n is the number of cases
 m is the number of matched controls

 p0 is probability that control is exposed
 p1 is the probability that a case is exposed
*/

if "`solve'"=="" local solve "n"

local or = `alt'

local q0 = 1-`p0'
local o1 = `or'*(`p0'/`q0')
local p1 = `o1'/(1+`o1')
local q1 = 1-`p1'

local p11 = `p1'*`p0'+`phi'* sqrt(`p0'*`p1'*`q0'*`q1')
local p10 = `p1'*`q0'-`phi'* sqrt(`p0'*`p1'*`q0'*`q1')
local p01 = `q1'*`p0'-`phi'* sqrt(`p0'*`p1'*`q0'*`q1')
local p00 = `q1'*`q0'+`phi'* sqrt(`p0'*`p1'*`q0'*`q1')

local p0p = `p11'/`p1'
local p0m = `p01'/`q1'

/* Not sure what to do with negatives?? it is possible with this phi parameterisation*/
/*
if `p0m'<0 local p0m 0
else if `p0m'>1 local p0m 1
*/

local q0p = 1-`p0p'
local q0m = 1-`p0m'

/*
local testphi = (`p11'*`p00'-`p10'*`p01')/sqrt(`p1'*`p0'*`q1'*`q0')
local testor = ( `p1'/`q1' )/( `p0'/`q0' )
local testp0m = `p0' - `phi'* sqrt(`p0'*`p1'*`q0'/`q1')
*/

local eor 0
local vor 0
local e1 0
local v1 0
forv i=1/`m' {
  qui _tk `i' `m' `p1' `p0p' `q0p' `q1' `p0m' `q0m'
  local eor = `eor'+ ( `i'*`r(tk)'*`or' ) / (`i'*`or'+`m'-`i'+1 )
  local vor = `vor'+ ( `i'*`r(tk)'*`or'*(`m'-`i'+1) ) / ((`i'*`or'+`m'-`i'+1 )^2 )
  local e1 = `e1'+ ( `i'*`r(tk)' ) / (`i'+`m'-`i'+1 )
  local v1 = `v1'+ ( `i'*`r(tk)'*(`m'-`i'+1) ) / ((`i'+`m'-`i'+1 )^2 )
}

if "`solve'"=="power" {
  local zalpha = invnorm(1-`alpha'/2)
  local prob = 1 + normprob(  ( sqrt(`n1')*(`e1'-`eor')-`zalpha'*sqrt(`v1') )/( sqrt(`vor') )  )  - ////
  normprob( ( sqrt(`n1')*(`e1'-`eor')+`zalpha'*sqrt(`v1') )/( sqrt(`vor') ) )

  di
  di as text "Estimate power for Matched Case Control Study"
  di "Test Ho: Odds ratio=1 Ha: Odds ratio = alt. OR
  di
  di "Assumptions:"
  di
  di as text "                Alpha = " as res %9.4f `alpha'
  di as text "(number of controls)M = " as res %9.4f `m'
  di as text "  Prob. Exp. Controls = " as res %9.4f `p0'
  di as text "               Alt OR = " as res %9.4f `alt'
  di as text "                    N = " as res %9.4f `n1'
  di
  di
  di as text "Estimated Power:"
  di
  di "       Power = " as res `prob'
  return local N_1 = `n1'
  return local N_2 = `n1'
  return local power= `prob'
}
if "`solve'"=="n" {
  local zalpha = invnorm(1-`alpha'/2)
  local zbeta = invnorm(`power')

  local n1 = ( ( `zbeta'*sqrt(`vor')+`zalpha'*sqrt(`v1') )^2 )/( (`e1'-`eor')^2 )
  local n = int(`n1')+1

  di
  di as text "Estimate Sample Size for Matched Case Control Study"
  di "Test Ho: Odds ratio=1 Ha: Odds ratio = alt. OR
  di
  di "Assumptions:"
  di
  di as text "                Alpha = " as res %9.4f `alpha'
  di as text "(number of controls)M = " as res %9.4f `m'
  di as text "  Prob. Exp. Controls = " as res %9.4f `p0'
  di as text "            Alt OR    = " as res %9.4f `alt'
  di as text "                Power = " as res %9.4f `power'
  di
  di
  di as text "Estimated Number of Cases:"
  di
  di "       N = " as res `n'

  return local N_1 = `n'
  return local N_2 = `n'
  return local power= `power'
}

end

/* Calculate the probability of observing k exposed subjects among a case and its m controls */

prog def _tk ,rclass
args k m p1 p0p q0p q1 p0m q0m

local temp = (`p1')* comb(`m',`k'-1)*(`p0p')^(`k'-1)*(`q0p')^(`m'-`k'+1)+(`q1')*comb(`m',`k')*(`p0m')^(`k')*(`q0m')^(`m'-`k')
return local tk = (`temp')

end