DeBUG with GLM BLUE

jiabowang · Apr 1, 2024 · 7fe3120 · 7fe3120
1 parent d3df1a5
commit 7fe3120
Show file tree

Hide file tree

Showing 3 changed files with 44 additions and 14 deletions.
diff --git a/R/GAPIT.EMMAxP3D.R b/R/GAPIT.EMMAxP3D.R
@@ -212,8 +212,25 @@
     X.beta <- X%*%beta
 
     beta.cv=beta
-    BLUE=X.beta
-
+    # BLUE=X.beta
+    # if(var(CVI[,2])!=0)
+    # {
+    #   XCV=cbind(1,as.matrix(CVI[,-1]))
+    # }else{
+    #   XCV=as.matrix(CVI[,-1])                               
+    # }
+    #CV.Extragenetic specified
+    # beta.Extragenetic=beta
+    XCVI=X[,-c(1:(1+CV.Extragenetic)),drop=FALSE]
+    XCVN=X[,c(1:(1+CV.Extragenetic)),drop=FALSE]
+    beta.I=as.numeric(beta)[-c(1:(1+CV.Extragenetic))]
+    beta.N=as.numeric(beta)[c(1:(1+CV.Extragenetic))]
+    BLUE.N=XCVN%*%beta.N
+    BLUE.I=rep(0,length(BLUE.N))
+    if(CV.Extragenetic!=0)BLUE.I=XCVI%*%beta.I
+    #Interception only
+    # if(length(beta)==1)XCV=X
+    BLUE=cbind(BLUE.N,BLUE.I)
     if(!is.null(K))
     {
       U.times.yv.minus.X.beta <- crossprod(U,(yv-X.beta))
@@ -858,17 +875,24 @@
                               XCVN=XCV[,c(1:(1+CV.Extragenetic)),drop=FALSE]
                               beta.I=as.numeric(beta)[-c(1:(1+CV.Extragenetic))]
                               beta.N=as.numeric(beta)[c(1:(1+CV.Extragenetic))]
-                              BLUE.I=try(XCVI%*%beta.I,silent=TRUE)
-                              BLUE.N=try(XCVN%*%beta.N,silent=TRUE)
+                              # print(is.null(beta.I))
+                              BLUE.I=rep(0,nrow(XCVI))
+                              BLUE.N=rep(0,nrow(XCVN))
+                              if(length(beta.I)>0)BLUE.I=try(XCVI%*%beta.I,silent=TRUE)
+                              if(length(beta.N)>0)BLUE.N=try(XCVN%*%beta.N,silent=TRUE)
                             }else{
                               XCVI=as.matrix(cbind(1,data.frame(CVI[,-1])))
                               beta.I=beta
-                              BLUE.I=try(XCVI%*%beta.I,silent=TRUE)
+                              BLUE.I=rep(0,length(BLUE.I))
+                              if(length(beta.I)>0)BLUE.I=try(XCVI%*%beta.I,silent=TRUE)
                               BLUE.N=rep(0,length(BLUE.I))
                             }
     #Interception only
                             # if(length(beta)==1)XCV=X
+                            # print(head(BLUE.N))
+                            # print(head(BLUE.I))
                             BLUE=cbind(BLUE.N,BLUE.I)
+                            # print(head(BLUE))
 #                           if(!is.null(CV.Extragenetic))
 #                           {
 #                             XCV=XCV[,-c(1:(1+CV.Extragenetic))]
@@ -987,7 +1011,7 @@
   Timmer=GAPIT.Timmer(Timmer=Timmer,Infor="GWAS done for this Trait")
   Memory=GAPIT.Memory(Memory=Memory,Infor="GWAS done for this Trait")
   #print("GAPIT.EMMAxP3D accomplished successfully!")
-
+  # print(head(BLUE))
     return(list(ps = ps, REMLs = -2*REMLs, stats = stats, effect.est = effect.est, rsquare_base = rsquare_base, rsquare = rsquare, dfs = dfs, df = df, tvalue = tvalue, stderr = stderr,maf=maf,nobs = nobs,Timmer=Timmer,Memory=Memory,
         vgs = vgs, ves = ves, BLUP = BLUP, BLUP_Plus_Mean = BLUP_Plus_Mean,
         PEV = PEV, BLUE=BLUE, logLM = logLM,effect.snp=effect.est,effect.cv=beta.cv))

diff --git a/R/GAPIT.Main.R b/R/GAPIT.Main.R
@@ -733,7 +733,7 @@ function(Y,
 
                 Timmer=p3d$Timmer
                 Memory=p3d$Memory
-
+                # print(head(p3d$BLUE))
                 Timmer=GAPIT.Timmer(Timmer=Timmer,Infor="Post PreP3D")
                 Memory=GAPIT.Memory(Memory=Memory,Infor="Post PreP3D")
 
@@ -1311,6 +1311,7 @@ BLUP=rep(0,nrow(CVI))
 PEV=rep(0,nrow(CVI))
 gs=NULL
 gs$BLUP=cbind(as.data.frame(CVI[,1]),Group,RefInf,ID,BLUP,PEV)
+# print(head(gs$BLUP))
 colnames(gs$BLUP)=c("Taxa","Group","RefInf","ID","BLUP","PEV")
 if(!byPass) 
 {
@@ -1406,6 +1407,7 @@ if((!byPass)&(!Model.selection)){
  print("GAPIT before BLUP and BLUE")
  #print(dim(p3d$BLUE))
  BLUE=data.frame(cbind(data.frame(CV.taxa),data.frame(p3d$BLUE)))
+ print(head(BLUE))
  colnames(BLUE)=c("Taxa","BLUE.N","BLUE.I")
  QTNs=rep(0,nrow(BLUE))
  #Initial BLUP as BLUe and add additional columns

diff --git a/R/GAPIT.Power.compare.R b/R/GAPIT.Power.compare.R
@@ -1,4 +1,4 @@
-`GAPIT.Power.compare` <-function(G=NULL,GD=NULL,GM=NULL,PCA.total=3,KI=NULL,CV=NULL,nrep=10,h2=0.85,NQTN=5,maxOut=100,model=c("GLM","FarmCPU"),seed=99163,WS=c(1e0,1e3,1e4,1e5,1e6,1e7),
+`GAPIT.Power.compare` <-function(G=NULL,GD=NULL,GM=NULL,PCA.total=3,KI=NULL,CV=NULL,nrep=10,h2=0.85,NQTN=5,maxOut=100,model=c("GLM","FarmCPU"),seed=99163,WS=c(1e0),
 myalpha=c(.01,.05,.1,.2,.3,.4,.5,.6,.7,.8,.9,1)){
 # Object: compare to Power against FDR for GLM,MLM,CMLM,ECMLM,SUPER
 # nrep:repetition times
@@ -70,40 +70,44 @@ utils::write.csv(cbind(myalpha,rep.Power.Alpha.store[[j]]),paste(h2,"_",NQTN,".P
 
 }
 
-grDevices::pdf(paste("GAPIT.Power.compare to multiple models", ".pdf", sep = ""), width = 4.5, height = 4.5,pointsize=9)
+
+for(k in 1:length(WS))
+{
+grDevices::pdf(paste("GAPIT.Power.compare to multiple models ",WS[k], ".pdf", sep = ""), width = 4.5, height = 4.5,pointsize=9)
 graphics::par(mar = c(5,6,5,3))
 	#win.graph(width=6, height=4, pointsize=9)
 	#palette(c("blue","red","green4","brown4","orange",rainbow(5)))
 	plot.color=rainbow(length(all.method))
 # print(head(rep.FDR.store[[j]]))
 kkt=NULL
+# print(head(rep.FDR.store))
+# print(head(rep.power.store))
 for(j in 1:length(all.method))
 {	
 	if(j==1)plot(as.numeric(as.matrix(rep.FDR.store[[j]])),as.numeric(as.matrix(rep.power.store[[j]])),bg="lightgray",xlab="FDR",ylab="Power",ylim=c(0,1),xlim=c(0,1),main="Power against FDR",type="o",pch=20,col=plot.color[j],cex=1.0,cex.lab=1.3, cex.axis=1, lwd=2,las=1)
    if(j!=1) graphics::lines(as.numeric(as.matrix(rep.power.store[[j]]))~as.numeric(as.matrix(rep.FDR.store[[j]])), lwd=2,type="o",pch=20,col=plot.color[j])
 	kkt<-cbind(kkt,as.numeric(as.matrix(rep.Power.Alpha.store[[j]])))
-
 }
 	graphics::legend("bottomright",c(all.method), pch = 20, lty =1,col=plot.color,lwd=2,cex=1.0,bty="n")
 grDevices::dev.off()
-
+colnames(kkt)=c(all.method)
 ###add type I error and power###
 utils::write.csv(cbind(myalpha,kkt),paste(h2,"_",NQTN,".Type I error.Power.by.multiple models.",nrep,".csv",sep=""))
 
 myalpha1<-myalpha/10
-grDevices::pdf(paste("GAPIT.Type I error_Power.compare to multiple models", ".pdf", sep = ""), width = 6, height = 4.5,pointsize=9)
+grDevices::pdf(paste("GAPIT.Type I error_Power.compare to multiple models ",WS[k], ".pdf", sep = ""), width = 6, height = 4.5,pointsize=9)
 graphics::par(mar = c(5,6,5,3))
 for(j in 1:length(all.method))
 {	
-
 	 if(j==1)plot(as.numeric(myalpha1),as.numeric(rep.Power.Alpha.store[[j]][,1]),log="x",bg="lightgray",xlab="Type I error",ylab="Power",main="Power against Type I error",type="o",pch=20,col=plot.color[j],cex=1.0,cex.lab=1.3, cex.axis=1, lwd=2,las=1,ylim=c(min(kkt),max(kkt)))
    if(j!=1) graphics::lines(rep.Power.Alpha.store[[j]][,1]~myalpha1, lwd=2,type="o",pch=20,col=plot.color[j])
 	 kkt<-cbind(kkt,rep.Power.Alpha.store[[j]][,1])
 
 }
 graphics::legend("bottomright",c(all.method), pch = 20, lty =1,col=plot.color,lwd=2,cex=1.0,bty="n")
 grDevices::dev.off()
-
+print("Power.Compare with multiple methods have been done!!!")
+}
 }#end compare to GLM,MLM,CMLM,ECMLM,SUPER
 #=============================================================================================