################################ # Chap 12: Association studies # ################################ p <- seq(0.1,0.9,by=0.2) R <- seq(1,5,by=0.1) plot(range(R),c(0,1),type="n",xlab="GRR", ylab="noncentrality") for (i in 1:length(p)) { ncp <- (R-1)*sqrt(p[i]*(1-p[i]))/(p[i]*(R-1)+1) lines(R,ncp,lty=i,lwd=1.5) } legend(1,1,legend=paste("p = ",p, sep=""),lty=1:length(p),lwd=rep(1.5,length(p))) CR <- read.table("CaseRandom.csv",header=TRUE,sep=",") summary(CR) table(CR$group,CR$snp1) chisq.test(table(CR$group,CR$snp1)) allele <- matrix(c(2,1,0,0,1,2),3,2) table(CR$group,CR$snp1)%*%allele chisq.test(table(CR$group,CR$snp1)%*%allele) p <- matrix(1,2,3) colnames(p) <- colnames(CR)[3:5] rownames(p) <- c("genotype","allele") for(m in 1:3) { tab <- table(CR$group,CR[,m+2]) p[1,m] <- chisq.test(tab)$p.value p[2,m] <- chisq.test(tab%*%allele)$p.value } round(p,6) n <- 500 pheno <- c(rep(1,n),rep(0,n)) pop <- c(rbinom(n,1,0.6),rbinom(n,1,0.4))+1 gam <- matrix(nrow=length(pheno),ncol=2) gam[pop==1,] <- rbinom(2*sum(pop==1),1,0.3) gam[pop==2,] <- rbinom(2*sum(pop==2),1,0.7) y <- rep(pheno,2) x <- as.vector(gam) z <- rep(pop,2) chisq.test(y,x) mantelhaen.test(y,x,z) glm(y~x,family=binomial()) summary(glm(y~x,family=binomial()))$coef summary(glm(y~x+z,family=binomial()))$coef na <- nr <- 500 n.mark <- 10^5 x <- qchisq(1-0.05/n.mark,df=1) F <- seq(0,4*10^(-4),length=100) v0 <- 1/(2*na) + 1/(2*nr) v1 <- F*(2*na-1)/(2*na) + v0 sig <- 1-exp(-(1-pchisq(x*v0/v1,df=1))*n.mark) plot(range(F),range(c(0,sig)),type="n",xlab="F", ylab="Sig. level") lines(F,sig,lty=2) abline(h=0.05) legend(0,max(sig),legend=c("Assumed significance", "Actual significance"),lty=1:2) (1/2^4)*100/choose(500,2) (1/2^6)*100/choose(500,2)