phenoSimu=function(X,h2,alpha,NQTN,distibution){ n=nrow(X) m=ncol(X) #Sampling QTN QTN.position=sample(m,NQTN,replace=F) SNPQ=as.matrix(X[,QTN.position]) QTN.position SNPQ[1:5,1:3] #Simulate phenotype if(distribution=="norm") {addeffect=rnorm(NQTN,0,1) }else {addeffect=alpha^(1:NQTN)} effect=SNPQ%*%addeffect effectvar=var(effect) residualvar=(effectvar-h2*effectvar)/h2 residual=rnorm(n,0,residualvar) y=as.data.frame(effect+residual) #View addeffect effect[1:3] effectvar residualvar residual[1:3] y[1:3,] return(list(addeffect = addeffect, y=y, add = effect, residual = residual,QTN.position=QTN.position,SNPQ=SNPQ)) }