# 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)
}