add functions

[don_microarray.git] / R / power_analysis.R

diff --git a/R/power_analysis.R b/R/power_analysis.R
new file mode 100644 (file)
index 0000000..3598910
--- /dev/null
+++ b/R/power_analysis.R
@@ -0,0 +1,76 @@
+# Power analysis for microarray data
+
+
+# Copyright 2003 By Don Armstrong
+
+# This program is free software; you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation; either version 2 of the License, or
+# (at your option) any later version.
+
+# This program is distributed in the hope that it will be useful, but
+# WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+# General Public License for more details.
+
+# You should have received a copy of the GNU General Public License
+# along with this program; if not, write to the Free Software
+# Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
+# 02111-1307, USA.
+
+
+
+# power.ma.partition
+# calcualate the power for data split into paritions of 200 genes.
+power.ma.partition <- function(dataset,partition=200,n=7,delta.mul=1,sig.level=0.05){
+  power.ma.temp <- dataset[names(sort(apply(dataset,1,mean))),]
+  if (partition>NROW(power.ma.temp)){
+    partition <- NROW(power.ma.temp)
+  }
+  power.ma.p.high <- NROW(power.ma.temp)
+  power.ma.return <- c()
+  for (power.ma.p.low in ceiling(seq(NROW(power.ma.temp),0,length=partition))){
+    if(power.ma.p.low == power.ma.p.high){
+      next;
+    }
+    power.ma.p.sd <- mean(apply(power.ma.temp[power.ma.p.low:power.ma.p.high,],1,sd,na.rm=T))
+    power.ma.p.mean <- mean(apply(power.ma.temp[power.ma.p.low:power.ma.p.high,],1,mean,na.rm=T))
+    power.ma.p.power <- power.t.test(n=n,sd=power.ma.p.sd,delta=power.ma.p.mean*delta.mul,sig.level=0.05)
+    power.ma.p.power <- power.ma.p.power$power
+    power.ma.return <- rbind(cbind(power.ma.p.mean,power.ma.p.power,power.ma.p.sd),power.ma.return)
+    power.ma.p.high <- power.ma.p.low -1;
+  }
+
+  return(power.ma.return)
+}
+
+# power.ma
+# calculate the power of individual genes
+power.ma <- function(dataset,n=7,delta.mul=1,sig.level=0.05){
+  power.ma.p.sd <- apply(dataset,1,sd,na.rm=T)
+  power.ma.p.mean <- apply(dataset,1,mean,na.rm=T)
+  power.ma.p.return <- cbind(power.ma.p.mean,power.ma.p.sd)
+  power.ma.p.power.calc <- function(sdmeanarray,n=7,delta.mul=1,sig.level=0.05){
+    power.ma.p.power.calc.power <- power.t.test(n=n,sd=sdmeanarray[2],delta=sdmeanarray[1]*delta.mul,sig.level=sig.level)
+    return(power.ma.p.power.calc.power$power)
+  }
+  power.ma.p.power <- apply(power.ma.p.return,1,power.ma.p.power.calc,delta.mul=delta.mul,sig.level=sig.level,n=n)
+  power.ma.p.return <- cbind(power.ma.p.return,power.ma.p.power)
+  
+  return(power.ma.p.return)
+}
+
+power.ma.sdmean <- function(dataset,n=7,delta.mul=1,sig.level=0.05){
+  power.ma.p.mean <- dataset[,1]
+  power.ma.p.sd <- dataset[,2]
+  power.ma.p.return <- cbind(power.ma.p.mean,power.ma.p.sd)
+  power.ma.p.power.calc <- function(sdmeanarray,n=7,delta.mul=1,sig.level=0.05){
+    power.ma.p.power.calc.power <- power.t.test(n=n,sd=sdmeanarray[2],delta=sdmeanarray[1]*delta.mul,sig.level=sig.level)
+    return(power.ma.p.power.calc.power$power)
+  }
+  power.ma.p.power <- apply(power.ma.p.return,1,power.ma.p.power.calc,delta.mul=delta.mul,sig.level=sig.level,n=n)
+  power.ma.p.return <- cbind(power.ma.p.return,power.ma.p.power)
+  
+  return(power.ma.p.return)
+}
+