# source("Rc2.pl")
# R code for generating power curves for Association (Case-Control) studies using SNPs
# inspired from "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.
# and Freidlin, Zheng, Li, Gastwirth (2002) Hum Hered 53:146-152
alpha=0.05; beta=0.2 # beta is prob Type II error
no.tests=1 # Can adjust no. tests to number of SNPs to be tested if
alpha=alpha/no.tests
# Bonferronni correction is desired
xs=c(0,1,1) # (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.8; # This is population frequency of dis.allele
K=0.01 # K is prevalence
N=3925 # N is sample size, including cases and controls
theta=0.2 # theta is desired ratio of number of cases to sample size.
step.size=0.1 # If you can't look at the code and see the effect of
# changing step.size, don't play w/ step.size. 
# FOR A PROTECTIVE "disease allele": Can use negative step size.
#################################################################
# Build the power curve
gamma2=1-step.size/2
power=0.05
g2.vec=NULL; pow.vec=NULL
while(power<0.98){ 
# gamma's are rel risks
gamma2=gamma2+step.size
if(xs[2]==0){ gamma1=1; model="recessive" }
if(xs[2]==0.5){ gamma1=(1+gamma2)/2; model="additive (on rel risk scale)" }
if(xs[2]==1){ gamma1=gamma2; model="dominant" }
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, respectively

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))
stopifnot(mu1!=0)
if(mu1>0){
	power=1-pnorm((qnorm(1-alpha/2)*sqrt(sigma1.2(theta, xs, p, q)+mu1^2)-N^(1/2)*mu1)/
		sqrt(sigma1.2(theta, xs, p, q)))
}else{
	power=pnorm((-qnorm(1-alpha/2)*sqrt(sigma1.2(theta, xs, p, q)+mu1^2)-N^(1/2)*mu1)/
		sqrt(sigma1.2(theta, xs, p, q)))
}
pow.vec=c(pow.vec, power)
g2.vec=c(g2.vec, gamma2)
} # end: while(power<0.90){ 
par(oma=c(1,1,3,0))
plot(g2.vec, pow.vec, xlab="Relative risk given 2 copies of disease allele",
ylab="Power", cex.lab=1.3, main="")
mtext(paste("Power curve for", model, "model"), line=1.5, outer=T, cex=1.4)
mtext(paste("Sample size =", N, " Prevalence=", K), line=-.5, outer=T)
no.cases=round(N*theta); no.controls=N-no.cases;
mtext(paste("No. cases=", no.cases, " No. controls=", no.controls),line=-2, outer=T)