Question: Limma, model with several factors

0

Jenny Drnevich •

**1.9k**wrote:Hi Gordon,
I'm confused by your response to Ingrid last May as to how to
construct a design matrix for her 2 factor experiment with repeated
measures on 1 factor. I came across this post in searching the mailing
list because I have a similar experimental design, but I was getting
the dreaded "Coefficients not estimable" warning message with what I
thought was the same design matrix you suggested to Ingrid. Your
example below actually produces a singular design matrix:
#reproducible example:
> Subject <- gl(4,2)
> treatment <- gl(2,2,8)
> timepoint <- gl(2,1,8)
> Treat.Time <- factor(paste(treatment,timepoint,sep="."))
> design <- model.matrix(~0+Treat.Time+Subject)
> colnames(design)[1:4] <- levels(Treat.Time)
> design
1.1 1.2 2.1 2.2 Subject2 Subject3 Subject4
1 1 0 0 0 0 0 0
2 0 1 0 0 0 0 0
3 0 0 1 0 1 0 0
4 0 0 0 1 1 0 0
5 1 0 0 0 0 1 0
6 0 1 0 0 0 1 0
7 0 0 1 0 0 0 1
8 0 0 0 1 0 0 1
attr(,"assign")
[1] 1 1 1 1 2 2 2
attr(,"contrasts")
attr(,"contrasts")$Treat.Time
[1] "contr.treatment"
attr(,"contrasts")$Subject
[1] "contr.treatment"
>
> is.fullrank(design)
[1] FALSE
>
> nonEstimable(design)
[1] "Subject4"
My searches of the mailing list led to a previous question also about
the same experimental design, in which your response seems to indicate
that you can't treat Subject as a fixed effect:
https://stat.ethz.ch/pipermail/bioconductor/2010-March/032523.html
Is this really true? I found this page
(http://ww2.coastal.edu/kingw/statistics/R-tutorials/repeated.html)
about modeling repeated measures factors in R, and it seems you could
model a single gene like this:
#not-reproducible
> aov.out = aov(expression.level ~ treatment * timepoint +
Error(Subject/timepoint), data=one.gene)
But this model formulation doesn't work with model.matrix(). So I'm
confused about the best way to model my experiment: I have 11 subjects
divided into 4 treatments with 3 timepoint measurements per subject
and we're interested in treatment effects, time effects and
interaction effects. I'm worried about simply treating Subject as a
random effect because when I look separately within each treatment,
only 1 of the 4 treatments actually has any substantial Subject
effect! (Subject correlation from duplicateCorrelation() for the
entire experiment is 0.038, but within each treatment the values are
-0.036, -0.042, -0.041 and 0.110). If none of them had an observable
Subject effect, I could just ignore it and do a standard 2-factor
ANOVA design, but because only one treatment has an effect I'm worried
the overall correlation value over-estimates the Subject effect in 3
treatments and under-estimates in the 1 treatment. Can I possibly drop
some of the Subject coefficients from
model.matrix(~0+Treat.Time+Subject) and only keep those pertaining to
the treatment that does show Subject effects?
I would greatly appreciate any advice or suggestions!
Jenny
-----Original Message-----
From: bioconductor-bounces@r-project.org [mailto:bioconductor-
bounces@r-project.org] On Behalf Of Gordon K Smyth
Sent: Monday, May 20, 2013 7:29 PM
To: Ingrid Dahlman [guest]
Cc: Bioconductor mailing list
Subject: [BioC] Limma, model with several factors
Dear Ingrid,
If you block on subject (as in Section 8.4 of the User's Guide), then
you have automatically blocked on center as well, because the subjects
are different in the two centers.
You experiment can be analyzed as suggested in Section 8.4.2 of the
User's Guide.
First combine treatment and timepoint into one factor (following
Section
8.5.2):
Treat.Time <- factor(paste(treatment,timepoint,sep="."))
Then make a design matrix including the blocking factor:
design <- model.matrix(~0+Treat.Time+Subject)
colnames(design)[1:4] <- levels(Treat.Time)
cont.matrix <- makeContrasts((F.A-F.B)-(P.A-P.B),levels=design)
etc.
The rule is that the experimental factors of interest are combined
into the treatment factors, whereas the nuisance factors are blocked
on.
Best wishes
Gordon
> Date: Mon, 20 May 2013 01:01:08 -0700 (PDT)
> From: "Ingrid Dahlman [guest]" <guest at="" bioconductor.org="">
> To: bioconductor at r-project.org, ingrid.dahlman at ki.se
> Subject: [BioC] Limma, model with several factors
>
>
> I have carefully read the Limma users guide but still have not
sorted
> out how I design the contrasts for the following Project.
>
> In want to compare the effect (After vs before) of two different
> treatments, F and P. The study was carried our in two different
centers.
> This can be illustrated as follows:
> Subject center treatment timepoint
> 1 1 F B
> 1 1 F A
> 2 1 P B
> 2 1 P A
> 3 2 F B
> 3 2 F A
> 4 2 P B
> 4 2 P A
> I want to compare (F.A-F.B)-(P.A-P.B), blocking for subject.
However,
> in addition, I would like to block for center. I.e. the center is
like
> a batch effect.
> Is it possible to block for two factors, suject and center, in the
> same test in Limma?
> /Ingrid
>
> -- output of sessionInfo():
>
> See section 8.7 in the user guide.
>
> --
> Sent via the guest posting facility at bioconductor.org.
>
______________________________________________________________________
The information in this email is confidential and
intend...{{dropped:9}}

ADD COMMENT
• link
•
modified 6.1 years ago
by
Gordon Smyth ♦

**38k**• written 6.1 years ago by Jenny Drnevich •**1.9k**