G2P=function(X,h2,alpha,NQTN,distribution,a2=0){ n=nrow(X) m=ncol(X) #Sampling QTN QTN.position=sample(m,NQTN,replace=F) SNPQ=as.matrix(X[,QTN.position]) QTN.position #QTN effects if(distribution=="normal") {addeffect=rnorm(NQTN,0,1) }else {addeffect=alpha^(1:NQTN)} #Simulate phenotype effect=SNPQ%*%addeffect effectvar=var(effect) #Interaction cp=0*effect nint=4 if(a2>0&NQTN>=nint){ for(i in nint:nint){ Int.position=sample(NQTN,i,replace=F) cp=apply(SNPQ[,Int.position],1,prod) } print(dim(cp)) cpvar=var(cp) intvar=(effectvar-a2*effectvar)/a2 if(is.na(cp[1]))stop("something wrong in simulating interaction") if(cpvar!=0){ print(c(effectvar,intvar,cpvar,var(cp),a2)) print(dim(cp)) cp=cp/sqrt(cpvar) cp=cp*sqrt(intvar) effectvar=effectvar+intvar }else{cp=0*effect} } residualvar=(effectvar-h2*effectvar)/h2 residual=rnorm(n,0,sqrt(residualvar)) y=effect+residual+cp return(list(addeffect = addeffect, y=y, add = effect, residual = residual, QTN.position=QTN.position, SNPQ=SNPQ,int=cp)) }