R/ma_loess_norm.R

   1 library(stats)
   2
   3
   4 # the appropriate method is to average over the \hat{M}'s to find the
   5 # average fit, and then use that to correct the data.
   6
   7 # ma.loessnorm <- function(array,iterations=2){
   8 #   n <- ncol(array)
   9 #   for (q in 1:iterations){
  10 #     bar.m.hat.array <- array(0,dim=dim(array))
  11 #     for (i in 1:(n-1)){
  12 #       for (j in (i+1):n){
  13 #         print(paste(i,j))
  14 #         # fit lowess
  15 #         # [we have to use loess so we can predict using predict.loess
  16 #         mva <- cbind(log2(array[,i]/array[,j]),0.5*log2(array[,i]*array[,j]))
  17 #         colnames(mva) <- c("m","a")
  18 #         mva <- data.frame(na.omit(mva))
  19 #         temp.loess <- loess(mva[,"m"]~mva[,"a"],control = loess.control(surface = "direct"))
  20 #         # predict for i, predict for j (\hat{M})
  21 #         # subtract out \hat{M} from M
  22 #         m.hat <- predict(temp.loess,mva[,"a"])
  23 #         # weight correction by the inverse of 2(n-1) . [so that there will be n adjustments]
  24 #         bar.m.hat.array[,i] <- bar.m.hat.array[,i]+m.hat/(n-1)
  25 #         bar.m.hat.array[,j] <- bar.m.hat.array[,j]+m.hat/(n-1)
  26 #         print(paste(i,j))
  27 #       }
  28 #     }
  29 #     # correct arrays
  30 #     result.array <- array(0,dim=dim(array))
  31 #     i <- 1
  32 #     for (j in (i+1):n){
  33 #       mva <- cbind(log2(array[,i]/array[,j]),0.5*log2(array[,i]*array[,j]))
  34 #       colnames(mva) <- c("m","a")
  35 #       if (j == 2){
  36 #         # correct Mi the first time through
  37 #         result.array[,i] <- 2^(mva[,"a"]+0.5*(mva[,"m"]-bar.m.hat.array[,i]))
  38 #       }
  39 #       result.array[,j] <- 2^(mva[,"a"]-0.5*(mva[,"m"]-bar.m.hat.array[,i]))
  40 #     }
  41 #     rownames(result.array) <- rownames(array)
  42 #     colnames(result.array) <- colnames(array)
  43 #     array <- result.array
  44 #   }
  45 #   return(array)
  46 # }
  47
  48
  49 ma.loessnorm2 <- function(array,iterations=2,small.positive=0.1,print.status=F){
  50   n <- ncol(array)
  51   # strip na values.
  52   array <- na.omit(array)
  53   for (q in 1:iterations){
  54     if (print.status){
  55       print(q)
  56     }
  57     result.array <- array(0,dim=dim(array))
  58     for (i in 1:(n-1)){
  59       for (j in (i+1):n){
  60         #print(paste(i,j))
  61         # replace 0's or negative values with a small positive value
  62         array[array[,i]<=0,i] <- small.positive
  63         array[array[,j]<=0,j] <- small.positive
  64         # fit lowess
  65         # [we have to use loess so we can predict using predict.loess
  66         mva <- cbind(log2(array[,i]/array[,j]),0.5*log2(array[,i]*array[,j]))
  67         colnames(mva) <- c("m","a")
  68         mva <- data.frame(mva)
  69         #temp.loess <- loess(m~a,mva,control = loess.control(surface = "direct"))
  70         temp.loess <- loess(m~a,mva,control=loess.control(surface="interpolate",statistics="approximate",
  71                                                           trace.hat="approximate",iterations=2))
  72         # predict for i, predict for j (\hat{M})
  73         # subtract out \hat{M} from M
  74         m.hat <- predict(temp.loess,mva[,"a"])
  75         # weight correction by the inverse of 2(n-1) . [so that there will be n adjustments]
  76         result.array[,i] <- result.array[,i]+2^(mva[,"a"]+0.5*(mva[,"m"]-m.hat))/(n-1)
  77         result.array[,j] <- result.array[,j]+2^(mva[,"a"]-0.5*(mva[,"m"]-m.hat))/(n-1)
  78       }
  79     }
  80     # correct arrays
  81     rownames(result.array) <- rownames(array)
  82     colnames(result.array) <- colnames(array)
  83     array <- result.array
  84   }
  85   return(array)
  86 }