pbsize {gap} | R Documentation |
Power for population-based association design
pbsize(kp, gamma = 4.5, p = 0.15, alpha = 5e-08, beta = 0.2)
kp |
population disease prevalence. |
gamma |
genotype relative risk assuming multiplicative model. |
p |
frequency of disease allele. |
alpha |
type I error rate. |
beta |
type II error rate. |
This function implements Scott et al. (1997) statistics for population-based association design. This is based on a contingency table test and accurate level of significance can be obtained by Fisher's exact test.
The returned value is scaler containing the required sample size.
Jing Hua Zhao extracted from rm.c.
Scott WK, Pericak-Vance MA, Haines JL (1997). “Genetic analysis of complex diseases.” Science, 275(5304), 1327; author reply 1329-30.
Long AD, Grote MN, Langley CH (1997). Genetic analysis of complex traits. Science 275: 1328.
Rosner B (2000). Fundamentals of biostatistics, 5 edition. Duxbury, Pacific Grove, CA. ISBN 9780534370688.
Armitage P, Colton T (eds.) (2005). Encyclopedia of biostatistics, 2 edition. John Wiley, Chichester, West Sussex, England ; Hoboken, NJ. ISBN 9780470849071.
kp <- c(0.01,0.05,0.10,0.2)
models <- matrix(c(
4.0, 0.01,
4.0, 0.10,
4.0, 0.50,
4.0, 0.80,
2.0, 0.01,
2.0, 0.10,
2.0, 0.50,
2.0, 0.80,
1.5, 0.01,
1.5, 0.10,
1.5, 0.50,
1.5, 0.80), ncol=2, byrow=TRUE)
outfile <- "pbsize.txt"
cat("gamma","p","p1","p5","p10","p20\n",sep="\t",file=outfile)
for(i in 1:dim(models)[1])
{
g <- models[i,1]
p <- models[i,2]
n <- vector()
for(k in kp) n <- c(n,ceiling(pbsize(k,g,p)))
cat(models[i,1:2],n,sep="\t",file=outfile,append=TRUE)
cat("\n",file=outfile,append=TRUE)
}
table5 <- read.table(outfile,header=TRUE,sep="\t")
unlink(outfile)
# Alzheimer's disease
g <- 4.5
p <- 0.15
alpha <- 5e-8
beta <- 0.2
z1alpha <- qnorm(1-alpha/2) # 5.45
z1beta <- qnorm(1-beta)
q <- 1-p
pi <- 0.065 # 0.07 and zbeta generate 163
k <- pi*(g*p+q)^2
s <- (1-pi*g^2)*p^2+(1-pi*g)*2*p*q+(1-pi)*q^2
# LGL formula
lambda <- pi*(g^2*p+q-(g*p+q)^2)/(1-pi*(g*p+q)^2)
# mine
lambda <- pi*p*q*(g-1)^2/(1-k)
n <- (z1alpha+z1beta)^2/lambda
cat("\nPopulation-based result: Kp =",k, "Kq =",s, "n =",ceiling(n),"\n")