#code to replicate simulations based on hypertension data set in Sugiyama et al. (2001) 
#run with R version 2.4.1 and R/qtl version 1.07 
library(pLODsims)
data(map10)

#This is the first of eight simulations run in parallel
#For the full set of simulations, replace 01 with the values 02, ..., 08
set.seed(52967665+01)

n.sim <- 500

simres.hyper <-  list()

#specify model parameters
RMSE<-sqrt(11431.61/242)
effects <- c(-2.7384, -2.9568, -6.5385, -2.1764, 3.3048, -2.9130, -7.6507)/RMSE # scaled to have sd=1
qtl.pos <- cbind(c(1,1,4,5,6,15), c(55.1, 100, 41.6, 80.7, 66.9, 37.3))
main.effects <- cbind(1:6, effects[1:6])

intxns <- rbind(c(5,6, effects[7]))

for(i in 1:n.sim) {  
  hypersim <- pLODsims.sim.cross(map10[1:19], qtl.pos=qtl.pos, qtl.main=main.effects, qtl.epis=intxns, n.ind=250, type="bc")
  hypersim <- calc.genoprob(hypersim, step=2.5)

  #pLOD_a with 5% thresholds
  simres.hyper[[i]] <- search.hk(hypersim, T1=2.8040, markers.only=1, refinement=1, 
	max.steps=100, max.df=50, fs2=1, verbose=0, pheno.col=1, restrict.diag=2, incl.intxns=0)
                   
  if(i==floor(i/10)*10){
    cat(i,"\n")
    save(simres.hyper, file="hyper01.RData", compress=T)
  }
}