# We use the following functions: OU.sim <- function(beta,markers,n.iter=10^4) { require(MASS) 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) } add.qtl <- function(Z,beta,markers,q,xi) { d <- dim(Z) if (length(markers) != d[2]) stop("Number of columns of simulated matrix does not match the number of markers") mu <- xi*exp(-beta*abs(markers-q)) Z <- sweep(Z,2,mu,"+") return(Z) } # Determine the threshold markers <- seq(0,80,by=10) beta <- 0.02 Z <- NULL for (i in 1:20) Z <- cbind(Z,OU.sim(beta,markers)) Z.max <- apply(abs(Z),1,max) quantile(Z.max,0.95) # Find the power q <- 35 xi <- 4.5 Z <- OU.sim(beta,markers) Z.1 <- add.qtl(Z,beta,markers,q,xi) mean(apply(abs(Z.1),1,max)> quantile(Z.max,0.95))