R/power_analysis.R

   1 # Power analysis for microarray data
   2
   3
   4 # Copyright 2003 By Don Armstrong
   5
   6 # This program is free software; you can redistribute it and/or modify
   7 # it under the terms of the GNU General Public License as published by
   8 # the Free Software Foundation; either version 2 of the License, or
   9 # (at your option) any later version.
  10
  11 # This program is distributed in the hope that it will be useful, but
  12 # WITHOUT ANY WARRANTY; without even the implied warranty of
  13 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14 # General Public License for more details.
  15
  16 # You should have received a copy of the GNU General Public License
  17 # along with this program; if not, write to the Free Software
  18 # Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
  19 # 02111-1307, USA.
  20
  21
  22
  23 # power.ma.partition
  24 # calcualate the power for data split into paritions of 200 genes.
  25 power.ma.partition <- function(dataset,partition=200,n=7,delta.mul=1,sig.level=0.05){
  26   power.ma.temp <- dataset[names(sort(apply(dataset,1,mean))),]
  27   if (partition>NROW(power.ma.temp)){
  28     partition <- NROW(power.ma.temp)
  29   }
  30   power.ma.p.high <- NROW(power.ma.temp)
  31   power.ma.return <- c()
  32   for (power.ma.p.low in ceiling(seq(NROW(power.ma.temp),0,length=partition))){
  33     if(power.ma.p.low == power.ma.p.high){
  34       next;
  35     }
  36     power.ma.p.sd <- mean(apply(power.ma.temp[power.ma.p.low:power.ma.p.high,],1,sd,na.rm=T))
  37     power.ma.p.mean <- mean(apply(power.ma.temp[power.ma.p.low:power.ma.p.high,],1,mean,na.rm=T))
  38     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)
  39     power.ma.p.power <- power.ma.p.power$power
  40     power.ma.return <- rbind(cbind(power.ma.p.mean,power.ma.p.power,power.ma.p.sd),power.ma.return)
  41     power.ma.p.high <- power.ma.p.low -1;
  42   }
  43
  44   return(power.ma.return)
  45 }
  46
  47 # power.ma
  48 # calculate the power of individual genes
  49 power.ma <- function(dataset,n=7,delta.mul=1,sig.level=0.05){
  50   power.ma.p.sd <- apply(dataset,1,sd,na.rm=T)
  51   power.ma.p.mean <- apply(dataset,1,mean,na.rm=T)
  52   power.ma.p.return <- cbind(power.ma.p.mean,power.ma.p.sd)
  53   power.ma.p.power.calc <- function(sdmeanarray,n=7,delta.mul=1,sig.level=0.05){
  54     power.ma.p.power.calc.power <- power.t.test(n=n,sd=sdmeanarray[2],delta=sdmeanarray[1]*delta.mul,sig.level=sig.level)
  55     return(power.ma.p.power.calc.power$power)
  56   }
  57   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)
  58   power.ma.p.return <- cbind(power.ma.p.return,power.ma.p.power)
  59
  60   return(power.ma.p.return)
  61 }
  62
  63 power.ma.sdmean <- function(dataset,n=7,delta.mul=1,sig.level=0.05){
  64   power.ma.p.mean <- dataset[,1]
  65   power.ma.p.sd <- dataset[,2]
  66   power.ma.p.return <- cbind(power.ma.p.mean,power.ma.p.sd)
  67   power.ma.p.power.calc <- function(sdmeanarray,n=7,delta.mul=1,sig.level=0.05){
  68     power.ma.p.power.calc.power <- power.t.test(n=n,sd=sdmeanarray[2],delta=sdmeanarray[1]*delta.mul,sig.level=sig.level)
  69     return(power.ma.p.power.calc.power$power)
  70   }
  71   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)
  72   power.ma.p.return <- cbind(power.ma.p.return,power.ma.p.power)
  73
  74   return(power.ma.p.return)
  75 }
  76