e <- 0.05 # the prevalence of exposure
g <- 0.1 # the prevalence of genotype
Re <- 1 # the relative risk for exposure alone
Rg <- 1 # the relative risk for genotype alone
Ri <- 2 # the effect of the gene-environment interaction
alpha <- 0.05
power <- 0.8
k <- 2 # how many time more control than case?
### Make a 2*2*2 table
table <- rep(0,8)
dim(table) <- c(2,2,2)
a <- rep(0,2)
exp <- table
table[1,1,1] <- g*e*Rg*Re*Ri
table[1,1,2] <- g*e
table[1,2,1] <- (1-e)*g*Rg
table[1,2,2] <- (1-e)*g
table[2,1,1] <- (1-g)*e*Re
table[2,1,2] <- (1-g)*e
table[2,2,1] <- (1-g)*(1-e)
table[2,2,2] <- (1-g)*(1-e)
table[,,1] <- table[,,1]/sum(table[,,1])
table[,,2] <- k*table[,,2]
### Mantel-Haenszel common odds ratio
Rmh <- (table[1,1,1]*table[1,2,2]/sum(table[1,,])
+ table[2,1,1]*table[2,2,2]/sum(table[2,,]))/
(table[1,1,2]*table[1,2,1]/sum(table[1,,])
+ table[2,1,2]*table[2,2,1]/sum(table[2,,]));
MH <- function(x) Rmh*(sum(table[i,1,])-x)*(sum(table[i,,1])-x)
- x*(sum(table[i,,])-sum(table[i,1,])-sum(table[i,,1])+x)
OR <- function(x) x[1,1]*x[2,2]/x[1,2]/x[2,1]
for(i in 1:2){
a[i]<-uniroot(MH, c(0,min(sum(table[i,1,]),sum(table[i,,1]))),tol = 1e-10)$root
exp[i,,] <- c(a[i], sum(table[i,,1])-a[i], sum(table[i,1,])-a[i],
sum(table[i,,])-sum(table[i,1,])-sum(table[i,,1])+a[i])
}
OR(exp[1,,])
OR(exp[2,,])
vn <- sum(1/exp)
va <- sum(1/table)
n <- (qnorm(1-alpha/2)*sqrt(vn) + qnorm(power)*sqrt(va))^2/log(Ri)^2
