rm(list=ls()) setwd(getwd()) source("subfunction.R") pc0=read.table(file="pc2a.txt",header=FALSE, skip=1) pc0=as.matrix(pc0) sample=800 K=20 M=0.01 L=100 rep=100 rep1=1000 model=2 RR=3 mm=4 N=L*rep pc=0 xx=matrix(0,sample,rep*L) pvper=matrix(0,rep1,mm) power=matrix(0,8,mm+3) nn0=0 alpha=0.05 for(model in c(0)) { if(model==0) npc=5 else npc=10 for(herit in c(0.01,0.015,0.02,0.025)) { for(i in 1:rep) { print(paste("model=",model,"herit=",herit,"rep=",i)) n1=(i-1)*L+1 n2=i*L xx[,n1:n2]=generate_genotype(sample, M, K, L) } yy=trait_nulla(sample,K,model,RR,2,4) lambda=uncorrect_lambda(yy,xx) h1=new_choose_h(yy, xx, pc0, npc) for(kk in c(10)) { nn0=nn0+1 ncau=kk/2 nprot=0 for(i in 1:rep1) { print(paste("model=",model,"herit=",herit,"kk=",kk,"rep=",i)) n1=(i-1)*kk+1 n2=i*kk x=xx[,n1:n2] y=gen_power(sample, x, herit, nprot, 0.01, ncau)+yy pvper[i,1]=uncorrect_region(y,x) pvper[i,2]=GC_lambda(pvper[i,1],lambda) pvper[i,3]=pc_linear_region(y, x, pc0[,1:10]) pvper[i,4]=pc_nonp_region(y, x, pc0, npc, h1) } power[nn0,1]=model power[nn0,2]=herit power[nn0,3]=kk for(k in 1:mm) { power[nn0,(k+3)]=sum(pvper[,k]