How edgeR estimate the common dispersion?
1
1
Entering edit mode
Dan • 0
@63648ae2
Last seen 11 weeks ago
United States

In the paper Small-sample estimation of negative binomial dispersion, with applications to SAGE data, section 4.1 Conditional maximum likelihood:

For RNA sequencing data, assume the counts for a single tag across n libraries is a negative binomial random variable. Consider Y1,,Yn as independent and NB(μi=miλ,ϕ), where mi is the library size (i.e. total number of tags sequenced for library i) and λ represents the proportion of the library that is a particular tag. The probability mass function is:

(1)f(yi;μ,ϕ)=P(Y=yi)=Γ(ϕ1+yi)Γ(ϕ1)Γ(yi+1)(1ϕ1μ1+1)yi(1ϕμ+1)ϕ1

E(Y)=μ and Var(Y)=μ+ϕμ2

If all libraries are the same size (i.e. mi≡m), the sum Z=Y1++YnNB(nmλ,ϕn1)

If we write out the log-likelihood, dropping terms that don't involve ϕ, we get:

logΓ(yi+ϕ1)nlogΓ(ϕ1)+yilog(ϕμ1+ϕμ)+nϕ1log(11+ϕμ)

Rewriting the log terms and rearranging gives us:

logΓ(yi+ϕ1)nlogΓ(ϕ1)+yilogϕμ(yi+nϕ1)log(1+ϕμ)

Some more rearranging:

logΓ(yi+ϕ1)nlogΓ(ϕ1)+zlogϕμ(z+nϕ1)log(1+ϕμ)

Substituting y¯=(1/n)yi for μ, as it is the MLE for μ and is what enables us to get rid of λ (which is hidden inside μ) gives us the one-dimensional optimization problem:

maxϕ[logΓ(yi+ϕ1)nlogΓ(ϕ1)+ny¯logϕy¯n(y¯+ϕ1)log(1+ϕy¯)]

How to derive conditional maximum likelihood for ϕ that is in the paper?

lY|Z=z(ϕ)=[i=1nlogΓ(yi+ϕ1)]+logΓ(nϕ1)logΓ(z+nϕ1)nlogΓ(ϕ1)

edgeR • 712 views
ADD COMMENT
2
Entering edit mode
Dan • 0
@63648ae2
Last seen 11 weeks ago
United States

I try to answer this question. Based on the Conditional Likelihood defined in https://rss.onlinelibrary.wiley.com/doi/epdf/10.1111/j.2517-6161.1996.tb02101.x, the conditional log-likelihood for ϕ of Y conditioned on Z, dropping terms that don't involve ϕ, is: lY|Z=z(y;ϕ)=logf(y;μ,ϕ)logfz(z;nμ,n1ϕ)

=logΓ(yi+ϕ1)nlogΓ(ϕ1)+yilog(ϕμ1+ϕμ)+nϕ1log(11+ϕμ)

logΓ(z+nϕ1)+logΓ(nϕ1)zlog(ϕμ1+ϕμ)nϕ1log(11+ϕμ)

=[i=1nlogΓ(yi+ϕ1)]+logΓ(nϕ1)logΓ(z+nϕ1)nlogΓ(ϕ1)

ADD COMMENT

Login before adding your answer.

Traffic: 855 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6