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.
7 # ma.loessnorm <- function(array,iterations=2){
9 # for (q in 1:iterations){
10 # bar.m.hat.array <- array(0,dim=dim(array))
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)
30 # result.array <- array(0,dim=dim(array))
33 # mva <- cbind(log2(array[,i]/array[,j]),0.5*log2(array[,i]*array[,j]))
34 # colnames(mva) <- c("m","a")
36 # # correct Mi the first time through
37 # result.array[,i] <- 2^(mva[,"a"]+0.5*(mva[,"m"]-bar.m.hat.array[,i]))
39 # result.array[,j] <- 2^(mva[,"a"]-0.5*(mva[,"m"]-bar.m.hat.array[,i]))
41 # rownames(result.array) <- rownames(array)
42 # colnames(result.array) <- colnames(array)
43 # array <- result.array
49 ma.loessnorm2 <- function(array,iterations=2,small.positive=0.1,print.status=F){
52 array <- na.omit(array)
53 for (q in 1:iterations){
57 result.array <- array(0,dim=dim(array))
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
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)
81 rownames(result.array) <- rownames(array)
82 colnames(result.array) <- colnames(array)