## S3 method for class 'DGEList':
glmFit(y, design=NULL, dispersion=NULL, prior.count=0.125, start=NULL, \dots)
## S3 method for class 'default':
glmFit(y, design=NULL, dispersion=NULL, offset=NULL, lib.size=NULL, weights=NULL,
prior.count=0.125, start=NULL, \dots)
glmLRT(glmfit, coef=ncol(glmfit$design), contrast=NULL)
DGEList
object with (at least) elements counts
(table of unadjusted counts) and samples
(data frame containing information about experimental group, library size and normalization factor for the library size)NULL
will be extracted from y
, with order of precedence: genewise dispersion, trended dispersions, common dispersion.y
giving offsets for the log-linear models. Can be a scalor or a vector of length ncol(y)
, in which case it is expanded out to a matrix.ncol(y)
giving library sizes. Only used if offset=NULL
, in which case offset
is set to log(lib.size)
. Defaults to colSums(y)
.DGEGLM
object, usually output from glmFit
.design
. Defaults to the last coefficient. Ignored if contrast
is specified.design
. If specified, then takes precedence over coef
.glmFit
produces an object of class DGEGLM
containing components counts
, samples
, genes
and abundance
from y
plus the following new components:nrow(y)
by ncol(design)
.nrow(y)
by ncol(design)
. It exists only when prior.count
is not 0.y
.glmLRT
produces objects of class DGELRT
with the same components as for glmfit
plus the following:y
containing the log2-fold-changes, likelhood ratio statistics and p-values, ready to be displayed by topTags
.table
contains the following columns:y
.glmFit
and glmLRT
implement generalized linear model (glm) methods developed by McCarthy et al (2012).glmFit
fits genewise negative binomial glms, all with the same design matrix but possibly different dispersions, offsets and weights.
When the design matrix defines a one-way layout, or can be re-parametrized to a one-way layout, the glms are fitting very quickly using mglmOneGroup
.
Otherwise the default fitting method, implemented in mglmLevenberg
, uses a Fisher scoring algorithm with Levenberg-style damping.
Positive prior.count
cause the returned coefficients to be shrunk in such a way that fold-changes between the treatment conditions are decreased.
In particular, infinite fold-changes are avoided.
Larger values cause more shrinkage.
The returned coefficients are affected but not the likelihood ratio tests or p-values.
glmLRT
conducts likelihood ratio tests for one or more coefficients in the linear model.
If coef
is used, the null hypothesis is that all the coefficients indicated by coef
are equal to zero.
If contrast
is non-null, then the null hypothesis is that the specified contrasts of the coefficients are equal to zero.
For example, a contrast of c(0,1,-1)
, assuming there are three coefficients, would test the hypothesis that the second and third coefficients are equal.
mglmOneGroup
or mglmLevenberg
.topTags
displays results from glmLRT
.
nlibs <- 3
ngenes <- 100
dispersion.true <- 0.1
# Make first gene respond to covariate x
x <- 0:2
design <- model.matrix(~x)
beta.true <- cbind(Beta1=2,Beta2=c(2,rep(0,ngenes-1)))
mu.true <- 2^(beta.true %*% t(design))
# Generate count data
y <- rnbinom(ngenes*nlibs,mu=mu.true,size=1/dispersion.true)
y <- matrix(y,ngenes,nlibs)
colnames(y) <- c("x0","x1","x2")
rownames(y) <- paste("gene",1:ngenes,sep=".")
d <- DGEList(y)
# Normalize
d <- calcNormFactors(d)
# Fit the NB GLMs
fit <- glmFit(d, design, dispersion=dispersion.true)
# Likelihood ratio tests for trend
results <- glmLRT(fit, coef=2)
topTags(results)
# Estimate the dispersion (may be unreliable with so few genes)
d <- estimateGLMCommonDisp(d, design, verbose=TRUE)
Run the code above in your browser using DataLab