Fit GaGa hierarchical model
Description
fitGG fits GaGa/MiGaGa hierarchical models, either via a fully Bayesian approach or via maximum likelihood.
fitNN fits a normalnormal hierarchical model (wrapper to call emfit in package EBarrays with the LNNMV model). fitNNSingleHyp is the same as fitNN but only considers the pattern that all groups are different from each other.
adjustfitNN corrects a small samplesize bias in the fitNN estimation procedure.
Usage
1 2 3 4 5 6 7 8 
Arguments
x 

groups 
If 
patterns 
Matrix indicating which groups are put together under
each pattern, i.e. the hypotheses to consider for each
gene. 
equalcv 

nclust 
Number of clusters in the MiGaGa model. 
method 

B 
Number of iterations to fit the model. For 
priorpar 
List with prior parameter values. It must have
components 
parini 
list with components 
trace 
For 
fit 

pitrue 
Grid of true 
nsim 
Number of datasets to simulate for each 
mc.cores 
If package 
Details
For GaGa/MiGaGa models, an approximation is used to sample faster from the
posterior distribution of the gamma shape parameters and to compute
the normalization constants (needed to evaluate the likelihood). These
approximations are implemented in rcgamma
and mcgamma
.
The cooling scheme in method=='SA'
uses a temperature equal to
1/log(1+i)
, where i
is the iteration number.
The EM implementation in method=='quickEM'
is a quick EM
algorithm that usually delivers hyperparameter estimates very similar
to those obtained via the slower method=='EM'
. Additionally,
the GaGa model inference has been seen to be robust to moderate
changes in the hyperparameter estimates in most datasets.
fitNN
is a wrapper to emfit
in package EBarrays with the LNNMV model.
This procedure estimates hyperparameters using the method of moments
(MoM), which typically results in overestimating the proportion of
differentially expressed genes, which we denote by pi1.
adjustfitNN
corrects this bias by repeatedly simulating from
the prior predictive of a normalnormal model. Simulations are
performed for a grid of pi1 values, so that the expected bias can be
evaluated for each of them. The bias is then modeled as a smooth
function of pi1 using function gam
from package mgcv
.
Finally, the value of pi1 is biasadjusted and the posterior
probabilities are recomputed using the updated pi1 value.
Value
fitGG
returns an object of class gagafit
, with components
parest 
Hyperparameter estimates. Only returned if 
mcmc 
Object of class 
lhood 
For 
nclust 
Same as input argument. 
patterns 
Same as input argument, converted to object of class

fitNN
returns an analogous object of class nnfit
. The
component nn.fit
is the object returned by emfit
.
Author(s)
David Rossell
References
Rossell D. (2009) GaGa: a Parsimonious and Flexible Model for Differential Expression Analysis. Annals of Applied Statistics, 3, 10351051.
Yuan, M. and Kendziorski, C. (2006). A unified approach for simultaneous gene clustering and differential expression identification. Biometrics 62(4): 10891098.
See Also
parest
to estimate hyperparameters and compute
posterior probabilities after a GaGa or MiGaGa
fit. findgenes
to find differentially expressed
genes. classpred
to predict the group that a new sample
belongs to.
Examples
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16  library(gaga)
set.seed(10)
n < 100; m < c(6,6)
a0 < 25.5; nu < 0.109
balpha < 1.183; nualpha < 1683
probpat < c(.95,.05)
xsim < simGG(n,m,p.de=probpat[2],a0,nu,balpha,nualpha,equalcv=TRUE)
x < exprs(xsim)
#Frequentist fit: EM algorithm to obtain MLE
groups < pData(xsim)$group[c(6,12)]
patterns < matrix(c(0,0,0,1),2,2)
colnames(patterns) < c('group 1','group 2')
gg1 < fitGG(x[,c(6,12)],groups,patterns=patterns,method='EM',trace=FALSE)
gg1 < parest(gg1,x=x[,c(6,12)],groups)
gg1
