[BioC] edgeR: GLM & residuals and model fitting & hypothesis testing

[Gordon K Smyth]

On Thu, 16 Feb 2012, Susanne Franssen wrote:

> Dear Gordon,
>
> thanks, a lot for you answer, however I am a little more uncertain now
> what to do.
>
>> 1) You may know that there are many different definitions of "residuals"
>> from a generalized linear models, and none of them have the properties of
>> residuals from normal linear models.  I have not found residuals very useful
>> myself in a genomic differential expression context.
>>
>> What sort of residuals did you want, and what were you planning to do with
>> them?
>
> Indeed I am not aware of this difference!
> I was thinking of the residual errors in a normal linear model. The
> problem in my analysis is that I only have 3 df (4 libs) to start
> with, so I cannot model the complete model
> individual+treatment+individual*treatment. So an idea was to model
> individual+treatment and afterwards do a model on the interaction only
> on the remaining residual errors as new observed vector.
> Is something like this at all possible or am I completely wrong here?

No, it is not possible.

Please read the section in the edgeR User's Guide "What to do if you have 
no replicates".

>> 2) You must always fit the model with all relevant factors:
>> individual*treatment.  You cannot do meaningful inference from a reduced
>> model until you have determined that the other factors are not important.
>> Hence you have to deal with the interaction first.
>
> It makes sense to have to put in all relevent factors to the model.
> However, concerning this I don't understand your comment - I have to
> deal with the interaction first, wouldn't the full model be:
> individual+treatment+individual*treatment (which I can't do because of
> a lack of df)

In R, the interaction is represented by individual:treatment.  The 
notation individual*treatment is a shorthand for interaction plus main 
effects, i.e., individual*treatment = 
individual+treatment+individual:treatment.

Please read the "Introduction to R" that comes with your R installation, 
in particular Section 11 on statistical models.

Best wishes
Gordon

> Cosidering my samples & factors:
>>> I have a fully crossed design with 2 factors and 2 factor levels each:
>>> individual <- as.factor(c("indA","indA","indB","indB"))
>>> treatment <- as.factor(c("treat1","treat2","treat1","treat2"))
>
> how would the cmd for the design matrix look like:
> design <- model.matrix(~??)
> and with which cmd would I fit the model and test for it:
> fit <- glmFit(d.GLM,design)
> lrt <- glmLRT(d.GLM,fit, coeff=?)
>
> Thanks a lot,
> Susanne
>
>
>
>>
>> Best wishes
>> Godon
>>
>>> Date: Tue, 14 Feb 2012 17:44:31 +0100
>>> From: Susanne Franssen <s.franssen at uni-muenster.de>
>>> To: bioconductor at r-project.org
>>> Subject: [BioC] edgeR: GLM & residuals and model fitting & hypothesis
>>>        testing
>>>
>>> Dear all,
>>>
>>>
>>> 1) GLM & residuals:
>>>
>>> I have a question concerning the use of GLMs in edgeR and the analysis
>>> of the residuals after model fitting.
>>>
>>> I have followed all the steps until model fitting, e.g.:
>>> glmfit.D <- glmFit(D, design, dispersion = D$tagwise.dispersion)
>>>
>>> The results I obtain from the fitting are the following catgories:
>>>>
>>>> names(glmfit.D)
>>>
>>> [1] "coefficients"  "fitted.values" "fail"          "not.converged"
>>> [5] "deviance"      "df.residual"   "abundance"     "design"
>>> [9] "offset"        "dispersion"    "method"        "counts"
>>> [13] "samples"
>>>
>>> What would be the best way to obtain the residuals for the "genewise"
>>> GLMs?
>>>
>>>
>>>
>>> 2) model fitting & hypothesis testing:
>>>
>>> I have a fully crossed design with 2 factors and 2 factor levels each:
>>> individual <- as.factor(c("indA","indA","indB","indB"))
>>> treatment <- as.factor(c("treat1","treat2","treat1","treat2"))
>>>
>>> in general I would be interested in 3 different aspects:
>>> a) effect of individual
>>> b) effect of treatment
>>> c) interaction between individual and treatment
>>>
>>> What would be the best way to test for those effects, would I rather
>>> test for all three aspects individually, i.e.:
>>> a) design <- model.matrix(~individual)
>>> b) design <- model.matrix(~treatment)
>>> c) design <- model.matrix(~individual*treatment)
>>>
>>> or doesn't it also make sense to model
>>> design <- model.matrix(~individual+treatment)
>>> and test for
>>> a) lrt.cd_ind <- glmLRT(D, glmfit.D, coef=2)
>>> b) lrt.cd_treat <- glmLRT(D, glmfit.D, coef=3)
>>> ... this way the effect of both factors could be accounted for in the
>>> model?!
>>>
>>>
>>> Thanks a lot!
>>> Susanne

______________________________________________________________________
The information in this email is confidential and intend...{{dropped:5}}