1 # Power analysis for microarray data
4 # Copyright 2003 By Don Armstrong
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.
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.
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
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)
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){
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;
44 return(power.ma.return)
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)
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)
60 return(power.ma.p.return)
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)
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)
74 return(power.ma.p.return)