# R code for calculating sample sizes for Association (Case-Control) studies using SNPs
# Could be used to reproduce columns "n" of Table II 
# of "On estimation of the variance in Cochran-Armitage trend tests for genetic 
# association using case-control studies" by Zheng and Gastwirth; 
# Statistics in Medicine 2006, 25 :3150-3157.

# An innovation with this code is that the Bonferronni correction can be applied.

# source("sample-size.R")
alpha=0.05; beta=0.2 # beta is prob Type II error
number.of.tests=1 # This can be adjusted for the Bonferronni adjustment
alpha=alpha/number.of.tests
xs=c(0,0,1) # c(0,0,1) is recessive model.  
# Put c(0,0.5,1) for additive on rel risk scale or
# c(0, 1, 1) for dominant model
dis.allele=0.1; 
theta=0.5 # theta is desired ratio of number of cases to sample size.
# In this example, case-control ratio is 1
K=0.1 # K is prevalence

gamma1=1; gamma2=3; # gamma's are rel risks
g=c((1-dis.allele)^2, 2*dis.allele*(1-dis.allele), dis.allele^2)

f=rep(NA, 3) # f's are penetrances
f[1]=K/(g[1]+gamma1*g[2]+gamma2*g[3])
f[2]=f[1]*gamma1; f[3]=f[1]*gamma2


p=f*g/K; q=(1-f)*g/(1-K) # p's and q's are cond probs given that one is a case or cont

sigma1.2=function(theta, xs, p, q){
sumx2p=sum(xs^2*p)-(sum(xs*p))^2
sumx2q=sum(xs^2*q)-(sum(xs*q))^2
return(theta*(1-theta)*((1-theta)*sumx2p+theta*sumx2q))
}
mu1=theta*(1-theta)*sum(xs*(p-q))
n=(qnorm(1-alpha/2)*(sigma1.2(1-theta, xs, p, q)+mu1^2)^0.5+qnorm(1-beta)*
	sqrt(sigma1.2(theta, xs, p, q)))^2/mu1^2