Power computations for differential expression
Description
powfindgenes
evaluates the posterior expected number of true positives
(e.g. true gene discoveries) if one were to obtain an additional batch
of data. It uses either a GaGa or a normalnormal model fit on a pilot
data set.
Usage
1 2  powfindgenes(fit, x, groups, batchSize = 1, fdrmax = 0.05, genelimit,
v0thre = 1, B = 1000, mc.cores=1)

Arguments
fit 
GaGa/MiGaGa or normalnormal model fit using pilot data

x 

groups 
If 
batchSize 
Number of additional samples to obtain per group. 
fdrmax 
Upper bound on FDR. 
.
genelimit 
Only the 
v0thre 
Only genes with posterior probability of being equally
expressed < 
B 
Number of simulations from the GaGa predictive distribution to be used to estimate the posterior expected number of true positives. 
mc.cores 
If 
Details
The routine simulates data from the posterior predictive distribution
of a GaGa or normalnormal model. That is, first it simulates parameter values (differential
expression status, mean expression levels etc.) from the posterior
distribution. Then it simulates data using the
parameter values drawn from the posterior.
Finally the simulated data is used to determine the differential status
of each gene, controlling the Bayesian FDR at the fdrmax
level,
as implemented in findgenes
.
As the differential expression status is known for each gene, one can
evaluate the number of true discoveries in the reported gene list.
In order to improve speed, hyperparameters are not reestimated when computing posterior probabilities for the posterior predictive simulated data.
Value
m 
Posterior expected number of true positives (as estimated by
the sample mean of 
s 
Standard error of the estimate i.e. SD of the simulations/sqrt(B) 
Author(s)
David Rossell
References
Rossell D. GaGa: a simple and flexible hierarchical model for microarray data analysis. http://rosselldavid.googlepages.com.
See Also
findgenes
, fitGG
, fitNN
,
parest
. See powclasspred
for
power calculations for sample classification.
Examples
1 2 3 4 5 6 7 8 9 10 11  #Simulate data and fit GaGa model
set.seed(1)
x < simGG(n=20,m=2,p.de=.5,a0=3,nu=.5,balpha=.5,nualpha=25)
gg1 < fitGG(x,groups=1:2,method='EM')
gg1 < parest(gg1,x=x,groups=1:2)
#Expected nb of TP for 1 more sample per group
powfindgenes(gg1,x=x,groups=1:2,batchSize=1,fdrmax=.05)$m
#Expected nb of TP for 10 more samples per group
powfindgenes(gg1,x=x,groups=1:2,batchSize=10,fdrmax=.05)$m
