[BioC] mixed-effects analysis within DEseq

[Mike Love]

Hi Akiko,

I believe that DESeq treats all effects as "fixed", i.e. parameters to 
be estimated.

Though with such small number of groups and replicates, it should not 
make a big difference in the estimates or the testing of the treatment 
effect.  For example, below I fit a generalized linear model and a 
generalized linear mixed model to Poisson counts of the design that you 
have (2 x 3 x 2).  The coefficient for the fixed effect, x, is almost 
identical, as is the change in deviance between the full and reduced 
model (this 435.33 which is the basis for the Chi-squared test).  I did 
not explore what might happen with interaction terms.

 > n <- 12
 > x <- factor(rep(0:1,each=n/2))
 > z <- factor(rep(c(0:2),times=n/3))
 > y <- rpois(n, exp(1 + as.integer(x) + as.integer(z)))
 > glm.fit1 <- glm(y ~ x + z, family=poisson)
 > coef(glm.fit1)
(Intercept)          x1          z1          z2
    2.919019    1.105920    0.949297    2.009069

 > library(lme4)
 > glmix.fit1 <- glmer(y ~ x + (1|z), family=poisson)
 > coef(glmix.fit1)
$z
   (Intercept)       x1
0    2.928748 1.105916
1    3.868466 1.105916
2    4.926723 1.105916

 > glm.fit0 <- glm(y ~ z, family=poisson)
 > glmix.fit0 <- glmer(y ~ (1|z), family=poisson)
 > anova(glm.fit0, glm.fit1, test="Chi")
Analysis of Deviance Table

Model 1: y ~ z
Model 2: y ~ x + z
   Resid. Df Resid. Dev Df Deviance  Pr(>Chi)
1         9     448.91
2         8      13.58  1   435.33 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

 > anova(glmix.fit0, glmix.fit1)
Data:
Models:
glmix.fit0: y ~ (1 | z)
glmix.fit1: y ~ x + (1 | z)
            Df    AIC    BIC   logLik  Chisq Chi Df Pr(>Chisq)
glmix.fit0  2 472.70 473.67 -234.351
glmix.fit1  3  39.37  40.82  -16.685 435.33      1  < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


best,

Mike


On 10/6/12 12:00 PM, bioconductor-request at r-project.org wrote:
>
> Message: 4
> Date: Fri, 05 Oct 2012 16:15:18 +0200
> From: Akiko Sugio <akiko.sugio at rennes.inra.fr>
> To: bioconductor at r-project.org
> Subject: [BioC] mixed-effects analysis within DEseq
> Message-ID: <506EEB76.1020102 at rennes.inra.fr>
> Content-Type: text/plain
>
> Dear all,
>
> We are currently analyzing RNA-seq data using DEseq and would like to
> ask for your advice regarding the possibility to treat a variable as a
> "random effect" rather than a "fixed effect".
>
> In our experiment, we are working with clonal individuals (aphids). We
> are working with three different clones. Each of these clones is fed
> with two different diets (treatment).We have two biological replicates
> for each condition. Hence, we have a total of 12 RNAseq samples (3
> different clones x two diet treatment x 2 replicates). We are mainly
> interested in identifying genes that are differentially expressed in
> response to the diet.
>
> Currently, we are running GLM models in DESeq:
>
> *Model0 <- fitNbinomGLMs( cds, count ~diet * cloneID )*
>
> We then test for the significance of each variable (i.e. cloneID, diet,
> and the interaction diet:cloneID) by comparing the model including the
> variable of interest to the model estimated without this variable.
>
> To test for the effect of the diet:
>
> *Model1 <- fitNbinomGLMs( cds, count ~cloneID)*
>
> *Model2<-fitNbinomGLMs( cds, count ~cloneID + diet )*
>
> *Diet_effect <- nbinomGLMTest(Model2, Model1)*
>
> *Diet_effect_adjusted< p.adjust(Diet_effect, method="BH" )*
>
> To test for the effect of the interaction:
>
> *Model0 <- fitNbinomGLMs( cds, count ~diet * cloneID ).*
>
> *Model3 <- fitNbinomGLMs( cds, count ~diet + cloneID ).*
>
> *Interaction_effect<- nbinomGLMTest(Model0, Model3)*
>
> *Interaction_effect _adjusted <- p.adjust(Interaction_effect, method="BH" )*
>
> For genes with a significant interaction between diet and CloneID, we do
> not interpret the effect of each variable (i.e. diet and CloneID).
>
> We are however not entirely satisfied by our model, because we think it
> would be better to treat the variable cloneID as a random effect. Is it
> possible to conduct such mixed-effects analysis within DEseq?
>
> Thank you very much in advance.
> Akiko
>