Entering edit mode

We have (very late) fixed this. The output of cqn now includes a new
component glm.offset, which fixes the issues Michael points out (and
one
other). From the NEWS file
\item THe output object from cqn() now has an additional
component:
glm.offset which is an offset matrix which can directly be used in
a
GLM type model (specifically edgeR). The usage is explained in
the
vignette secion on Import into edgeR. Previously the vignette
recommended using the offset component of the cqn output, which is
wrong, due to a scaling issue. The offset component of cqn is
unchanged. This bug was found by Mike Love
\email{love@molgen.mpg.de} and fixed in CQN 1.5.1.
Kasper
On Wed, Feb 20, 2013 at 4:11 AM, Michael Love <love@molgen.mpg.de>
wrote:
> hi,
>
> re: EDASeq, this was my mistake. The vignette does in fact have the
> multiplication by -1, but my eyes scanned over this:
>
> > library(edgeR)
> > design <- model.matrix(~conditions, data=pData(dataOffset))
> > disp <- estimateGLMCommonDisp(exprs(dataOffset),
> + design, offset=-offst(dataOffset))
> > fit <- glmFit(exprs(dataOffset), design, disp,
> offset=-offst(dataOffset))
>
> Then, I still believe that the cqn offsets should be multiplied by
> -log(2) to be handed to edgeR, which is not currently in the cqn
vignette:
>
> > design <- model.matrix(~ d.mont$sample$group)
> > d.mont.cqn <- estimateGLMCommonDisp(d.mont, design = design,
> + offset = cqn.subset$offset)
>
> But I think it would be desirable for packages on both sides of the
> hand-off to provide an explicit formula for what is meant by
"offset",
> which could be one of:
>
> log(mu) = X * beta + offset
> log2(mu) = X * beta + offset
> log(mu) + offset = X * beta
> log2(mu) + offset = X * beta
>
> best,
>
> Mike
>
>
> On 02/14/2013 06:53 PM, Mike Love wrote:
> > One more thing, another wrinkle is log or log2. For this reason, I
> > shouldn't have used correlation in the example. I believe cqn is
> > giving log2 offsets, EDASeq natural log, and edgeR is accepting
> > natural log.
> >
> > thanks,
> >
> > Mike
> >
> > On 2/14/13 5:21 PM, Michael Love wrote:
> >> hi,
> >>
> >> I could have things totally mixed up here. I was thinking that
GLM
> >> offsets would be positively correlated with log counts.
> >>
> >> for instance in edgeR, mglmLS.R:
> >>
> >> mu <- exp(beta %*% t(X) + offset)
> >>
> >> but it seems offsets from packages EDASeq and cqn are negative
> >> correlated with log counts.
> >>
> >> for example some poisson data:
> >>
> >> > n <- 1000
> >> > covariate <- rnorm(n,6,.5)
> >> > counts <- replicate(2,rpois(n,exp(covariate)))
> >>
> >> > library(EDASeq)
> >> > edaseq.offset <-
> >>
withinLaneNormalization(x=counts,y=covariate,which="full",offset=TRUE)
> >> > cor(edaseq.offset[,1],log(counts[,1]+1))
> >> [1] -0.9765923
> >>
> >> > library(cqn)
> >> > cqn.offset <-
> >>
> cqn(counts,x=covariate,lengths=rep(1,n),lengthMethod="fixed",sizeFac
tors=c(1,1))$offset
> >> > cor(cqn.offset[,1],log(counts[,1]+1))
> >> [1] -0.9950717
> >>
> >>
> >> So EDASeq and cqn are giving similar results here:
> >>
> >> > cor(edaseq.offset[,1],cqn.offset[,1])
> >> [1] 0.9847672
> >>
> >> If I give these to the standard glm function as offsets I would
> >> expect the coefficient for 'covariate' to go from 1 to 0. But it
> >> gives me something like a coefficient of 2, and to get zero I
have to
> >> give it -1 * offset:
> >>
> >> without offset:
> >>
> >> > coef(glm(counts[,1] ~ covariate, family="poisson"))
> >> (Intercept) covariate
> >> -0.02735053 1.00467285
> >>
> >> with offset:
> >>
> >> > coef(glm(counts[,1] ~ covariate, offset=edaseq.offset[,1],
> >> family="poisson"))
> >> (Intercept) covariate
> >> -5.860094 1.954997
> >> > coef(glm(counts[,1] ~ covariate, offset=cqn.offset[,1],
> >> family="poisson"))
> >> (Intercept) covariate
> >> -8.752786 2.457637
> >>
> >> with -1 * offset:
> >>
> >> > coef(glm(counts[,1] ~ covariate, offset=-edaseq.offset[,1],
> >> family="poisson"))
> >> (Intercept) covariate
> >> 5.8253511 0.0498681
> >> > coef(glm(counts[,1] ~ covariate, offset=-cqn.offset[,1],
> >> family="poisson"))
> >> (Intercept) covariate
> >> 8.6975896 -0.4482184
> >>
> >>
> >> thanks for any guidance,
> >>
> >> Mike
> >>
> >
>
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> Bioconductor mailing list
> Bioconductor@r-project.org
> https://stat.ethz.ch/mailman/listinfo/bioconductor
> Search the archives:
> http://news.gmane.org/gmane.science.biology.informatics.conductor
>
[[alternative HTML version deleted]]