R/fdr.spatial2.R

fdr.spatial2 <- function (object, delta = 2, N = 100, av = "median", edgeNA = FALSE) 
{
   
 if (!(class(object)=="marrayRaw") & !(class(object)=="marrayNorm")){
  stop("Object should be of class marrayRaw or marrayNorm")
}
   
     LL <- list(NULL) 
     FDRpL <- list(NULL)
     FDRnL <- list(NULL)
     ML <- maM(object)   
     index <- c(1:dim(object)[[2]])
  
    

    for (ii in index){
   
      X <- v2m(ML[,ii],Ngr=maNgr(object),Ngc=maNgc(object),Nsr=maNsr(object),Nsc=maNsc(object))
   
    


    XavP <- double(length(X) * N) + NA

    ### GENERATING EMPIRICAL BACKGROUND DISTRIBUTION BASED ON RANDOM PERMUTATION 
    for (i in 1:N) {
        tmp <- ma.matrix(matrix(sample(as.vector(X)), ncol = dim(X)[[2]]), 
            delta = delta, av = av, edgeNA = edgeNA)
        XavP[((i - 1) * length(X) + 1):(i * length(X))] <- as.vector(tmp)
    }
    XavP <- XavP[!is.na(XavP)]
    XavP.l <- length(XavP)

    #### AV. M STATISTICS FOR ORIGINAL DATA 
    Xav <- as.vector(ma.matrix(X, delta = delta, av = av, edgeNA = edgeNA))
    o <- 1:length(Xav)
    ro <- o[rank(Xav)]
    XavS <- sort(Xav)
    XavS.l <- length(XavS)
    XN <- double(length = length(XavS)) + NA

    #### COMPARINING STATISTIC OF ORIGINAL DATA WITH EMPIRICAL BACKGROUND DISTRIBUTION 
    for (i in 1:XavS.l) {
        XN[i] <- sum(XavP >= XavS[i], na.rm = TRUE)
    }
    XN <- XN/(XavP.l/XavS.l)


    #### DETERMINING FALSE DISCOVERY RATES
    pFDR <- double(length = length(Xav)) + NA
    for (i in 1:XavS.l) {
        pFDR[XavS.l - i + 1] <- XN[XavS.l - i + 1]/(XN[XavS.l - 
            i + 1] + i)
    }
    pFDR[pFDR == 0] <- 1/(XavS.l * N)
    nFDR <- double(length = length(Xav)) + NA
    for (i in 1:XavS.l) {
        nFDR[i] <- (XavS.l - XN[i])/((XavS.l - XN[i]) + i)
    }
    nFDR[nFDR == 0] <- 1/(XavS.l * N)
    pFDR <- matrix(pFDR[ro], ncol = dim(X)[[2]])
    nFDR <- matrix(nFDR[ro], ncol = dim(X)[[2]])
      
        FDRpL[[ii]] <-  m2v(pFDR,Ngr=maNgr(object),Ngc=maNgc(object),Nsr=maNsr(object),Nsc=maNsc(object))
        FDRnL[[ii]]  <-  m2v(nFDR,Ngr=maNgr(object),Ngc=maNgc(object),Nsr=maNsr(object),Nsc=maNsc(object))
      
      # FDRpL[[ii]] <-   pFDR; 
      # FDRnL[[ii]] <-   nFDR;   
   }
list(FDRp=FDRpL,FDRn=FDRnL)
}
##########################################################################

Try the OLIN package in your browser

Any scripts or data that you put into this service are public.

OLIN documentation built on Nov. 8, 2020, 7:44 p.m.