# 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) } power.qtl <- function(xi,p.power,Z,beta,markers,q,z.critical) { p <- mean(apply(abs(add.qtl(Z,beta,markers,q,xi)),1,max)> z.critical) return(p-p.power) } mouse <- 30 primer <- 70 genotype <- 2 sig.level <- 0.95 p.power <- 0.85 beta <- 0.02 minimal.effect <- function(n.mark) { # Determine the threshold markers <- seq(0,80,length=n.mark) Z <- NULL for (i in 1:20) Z <- cbind(Z,OU.sim(beta,markers)) Z.max <- apply(abs(Z),1,max) z.critical <- quantile(Z.max,sig.level) # Find the ncp q <- markers[2]/2 Z <- OU.sim(beta,markers) ncp <- uniroot(power.qtl,c(1,10),p.power=p.power, Z=Z,beta=beta,markers=markers,q=q,z.critical=z.critical)$root # Determine the number of mice and minimal effect n.mice <- (50000 - n.mark*20*primer)/(n.mark*20*genotype + mouse) genetic.effect <- 2*ncp/sqrt(n.mice) results <- c(genetic.effect, n.mice) names(results) <- c("genetic.effect","n.mice") return(results) } # Running minimal.effect(2) minimal.effect(3) minimal.effect(4) minimal.effect(5)