###################################################### # Chap 5: Whole genome scans: the significance level # ###################################################### OU.sim <- function(beta,markers,n.iter=10000) { V <- outer(markers,markers,function(x,y) exp(-beta*abs(x-y))) mu <- rep(0,length(markers)) Z <- mvrnorm(n.iter,mu,V) return(Z) } markers <- seq(0,80,by=20) beta <- 0.02 outer(markers,markers,function(x,y) exp(-beta*abs(x-y))) OU.sim(beta,markers,10) library(MASS) OU.sim(beta,markers,10) markers <- seq(0,80,by=10) Z <- OU.sim(beta,markers,30) plot(range(markers),c(-3,3),type="n",xlab="cM",ylab="Z") for (i in 1:30) lines(markers,Z[i,],col=gray(0.7*(i-1)/30)) abline(h=c(-1.96,1.96),lty=2) Z <- NULL for (i in 1:20) Z <- cbind(Z,OU.sim(beta,markers)) dim(Z) mean(abs(Z[,1]) >= 1.96) mean(apply(abs(Z),1,max) >= 1.96) plot(c(0,5),c(0,1),type="n",xlab="|Z|",ylab="Density") lines(density(Z[,1],from=0),lwd=1.5) lines(density(apply(abs(Z),1,max)),col=gray(0.4)) legend(0,1,legend=c("|Z_1|","max_t |Z_t|"),lty=c(1,1), col=gray(c(0,0.4)),lwd=c(1.5,1)) mean(apply(abs(Z),1,max) >= 3.50) Z.max <- apply(abs(Z),1,max) z.val <- seq(3.52,3.58,by=0.01) p <- NULL for(z in z.val) p <- c(p, mean(Z.max >= z)) names(p) <- z.val p n.tests <- c(1,2,3,5,9,17,81)*20 bon <- qnorm(1-0.025/n.tests) names(bon) <- paste("cM", c(80,40,25,20,10,5,1),sep="") bon markers <- 0:80 Z <- NULL for (i in 1:20) Z <- cbind(Z,OU.sim(beta,markers)) marker.set <- list(cM80=c(41),cM40=c(21,61),cM25=c(16,41,66), cM20=c(1,21,41,61,81),cM10=seq(1,81,by=10),cM5=seq(1,81,by=5)) marker.set <- lapply(marker.set, function(m) as.vector(outer(m,(0:19)*81,'+'))) Z.max <- apply(abs(Z),1,max) for(i in seq(length(marker.set),1)) Z.max <- cbind(apply(abs(Z[,marker.set[[i]]]),1,max),Z.max) colnames(Z.max) <- names(bon) p.bon <- NULL for(i in 1:length(bon)) p.bon <- c(p.bon, mean(Z.max[,i] >= bon[i])) names(p.bon) <-colnames(Z.max) p.bon Nu <- function(x) { y <- x/2 (1/y)*(pnorm(y)-0.5)/(y*pnorm(y) + dnorm(y)) } OU.approx <- function(z,beta,Delta, length=1600,chr=20,center=0,test="two-sided") { d <- switch(test,"two-sided"=2,"one-sided"=1) p <- 1-exp(-d*chr*(1-pnorm(z)) -d*beta*length*z*dnorm(z)*Nu(z*sqrt(2*beta*Delta))) return(p-center) } app <- bon app[1] <- uniroot(OU.approx,c(2,4), beta=0.02,Delta=80,length=0,center=0.05)$root app[2] <- uniroot(OU.approx,c(2,4), beta=0.02,Delta=40,length=20*(60-20),center=0.05)$root app[3] <- uniroot(OU.approx,c(2,4), beta=0.02,Delta=25,length=20*(65-15),center=0.05)$root app[4] <- uniroot(OU.approx,c(2,4),beta=0.02,Delta=20,center=0.05)$root app[5] <- uniroot(OU.approx,c(2,4),beta=0.02,Delta=10,center=0.05)$root app[6] <- uniroot(OU.approx,c(2,4),beta=0.02,Delta=5,center=0.05)$root app[7] <- uniroot(OU.approx,c(2,4),beta=0.02,Delta=1,center=0.05)$root app p.app <- NULL for(i in 1:length(app)) p.app <- c(p.app,mean(Z.max[,i] >= app[i])) rbind(p.bon,p.app)