[BioC] Desing matrix in limma

[Dario Beraldi]

Hi Jim,

Many thanks for the timely reply!
I see your point and I didn't realize that duplicateCorrelation() can  
be used to model random effects.
If I understand it correctly, anim_id as random effect would make more  
sense since we are interested in the difference between these specific  
time points, but not in the difference between these specific animals.  
Also, I guess anim_id as random would consume 1 df instead of 2 as  
fixed.

Anyway, glad to hear my reasoning was fine!

All the best
Dario


Quoting "James W. MacDonald" <jmacdon at med.umich.edu> on Wed, 13 Apr  
2011 09:59:54 -0400:

> Hi Dario,
>
> On 4/13/2011 6:46 AM, Dario Beraldi wrote:
>> Hello,
>>
>> I'm using limma to analyze affymetrix arrays from a time course
>> experiment (time points 0, 2, 7, 24). We are interested in comparing 0
>> vs 2, 0 vs 7, 0 vs 24.
>>
>> The experimental design looks like this:
>>
>> anim_id<- as.factor(rep(c('L1', 'L2', 'L3'), each= 4))
>> time_point<- rep(c(0, 2, 7, 24), 3)
>> slide_id<- paste(anim_id, time_point, sep= '.')
>> (targets<- data.frame(slide_id, anim_id, time_point))
>>
>> # slide_id anim_id time_point
>> #1 L1.0 L1 0
>> #2 L1.2 L1 2
>> #3 L1.7 L1 7
>> #4 L1.24 L1 24
>> #5 L2.0 L2 0
>> #6 L2.2 L2 2
>> #7 L2.7 L2 7
>> #8 L2.24 L2 24
>> #9 L3.0 L3 0
>> #10 L3.2 L3 2
>> #11 L3.7 L3 7
>> #12 L3.24 L3 24
>>
>> I wanted to analyze all the arrays at once to fully take advantage of
>> the eBayes function. Also I wanted to capitalize on the fact that the
>> same animal is sampled at each time point (0, 2, 7, 24 hours).
>>
>> So I set up the design matrix like this:
>>
>> design<- model.matrix(~0+anim_id)
>> rownames(design)<- paste(targets$anim_id, targets$time_point, sep= '.')
>> design<- cbind(design, tp.2vs0= c(0,1,0,0, 0,1,0,0, 0,1,0,0)) ## 0 vs 2
>> design<- cbind(design, tp.7vs0= c(0,0,1,0, 0,0,1,0, 0,0,1,0)) ## 0 vs 7
>> design<- cbind(design, tp.24vs0= c(0,0,0,1, 0,0,0,1, 0,0,0,1)) ## 0 vs 24
>> design
>> # anim_idL1 anim_idL2 anim_idL3 tp.2vs0 tp.7vs0 tp.24vs0
>> # L1.0 1 0 0 0 0 0
>> # L1.2 1 0 0 1 0 0
>> # L1.7 1 0 0 0 1 0
>> # L1.24 1 0 0 0 0 1
>> # L2.0 0 1 0 0 0 0
>> # L2.2 0 1 0 1 0 0
>> # L2.7 0 1 0 0 1 0
>> # L2.24 0 1 0 0 0 1
>> # L3.0 0 0 1 0 0 0
>> # L3.2 0 0 1 1 0 0
>> # L3.7 0 0 1 0 1 0
>> # L3.24 0 0 1 0 0 1
>>
>> Now I can procede 'as usual':
>> fit <- lmFit(limma_affy, design)
>> fit <- eBayes(fit)
>>
>> topTable(fit, coef= 4) ## 0 vs 2
>> topTable(fit, coef= 5) ## 0 vs 7
>> topTable(fit, coef= 6) ## 0 vs 24
>>
>> I verified that the logFC I get from topTable is the same as the logFC
>> calculated manually, which tells me that the design is correct.
>>
>> However, I would like to make sure that this procedure is correct, Would
>> anyone be so kind to confirm it?
>
> It is correct, in the strict sense that you are doing what you think  
> you are doing.
>
> However, one could argue that it would be a better idea to fit the  
> animals as a random effect rather than a fixed effect, using  
> duplicateCorrelation() to account for the intra-animal covariance  
> structure.
>
> There are multiple examples of using duplicateCorrelation() in the  
> Limma User's Guide.
>
> Best,
>
> Jim
>>
>> All the best
>>
>> Dario
>>
>>
>
> -- 
> James W. MacDonald, M.S.
> Biostatistician
> Douglas Lab
> University of Michigan
> Department of Human Genetics
> 5912 Buhl
> 1241 E. Catherine St.
> Ann Arbor MI 48109-5618
> 734-615-7826
> **********************************************************
> Electronic Mail is not secure, may not be read every day, and should  
> not be used for urgent or sensitive issues



-- 

Dr. Dario Beraldi
Institute of Evolutionary Biology
University of Edinburgh
West Mains Road
Edinburgh EH9 3JT
Scotland, UK


-- 
The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.