[BioC] Design and not estimable coefficients
Gordon K Smyth
smyth at wehi.EDU.AU
Tue May 1 02:27:04 CEST 2012
Dear Dave,
You are trying to fit an interaction model between exp and mutation, which
is unnecessary and inappropriate because it makes the whole concept of a
mutation effect debatable.
The design I already suggested for you is correct. If you want to adjust
for a batch effect, you need an additive model, not a crossed interaction
model.
You can easily see for yourself, without having to ask me, that your
design matrix is not equivalent, because it has four columns, whereas the
one I suggested has three.
Best wishes
Gordon
---------------------------------------------
Professor Gordon K Smyth,
Bioinformatics Division,
Walter and Eliza Hall Institute of Medical Research,
1G Royal Parade, Parkville, Vic 3052, Australia.
smyth at wehi.edu.au
http://www.wehi.edu.au
http://www.statsci.org/smyth
On Mon, 30 Apr 2012, Dave Canvhet wrote:
> Dear Gordon,
>
> many thanks for your answer !!
>
> Clearly, I need to get some more information about the Intercept terms,
> which is probably a classical parameters in statistics and linear modeling
> (whom I'm not familiar with ...).
>
>
> Also, Could you please give me your opinion using the secodn way of
> designing the design matrix summarized below ?
>
> f <- paste(target$mutation,target$**exp,sep="")
> design <- model.matrix(~0+f)
> colnames(design) <- levels(f)
> using such a design matrix, can I use correctly the following contrast
> matrix to get the gene differentially expressed between mutant and wild
> type ?
> cont.matrix <- makeContrasts(mu_vs_wt = (mu2+mu1) - (wt2+wt1),levels=design)
>
>
> again many thanks for your time
> --
> Dave
>
> 2012/4/26 Gordon K Smyth <smyth at wehi.edu.au>
>
>> Dear Dave,
>>
>> This is a common type of design, with two groups and a batch effect. To
>> find genes DE in mu vs wt adjusting for exp differences:
>>
>> mutation <- factor(target$mutation)
>> mutation <- relevel(mutation, ref="wt")
>> exp <- factor(target$exp)
>> design <- model.matrix(~exp+mutation)
>> fit <- lmFit(x,design)
>> fit <- eBayes(fit)
>> topTable(fit, coef=3)
>>
>> Note that the design matrix has only 3 columns, not four as yours do. The
>> batch effect requires only one extra column.
>>
>> Your design matrix has a redundant column because you've defined the
>> intercept term twice. The first and second columns (wt and mu) add up to
>> the intercept column, and so do the 3rd and 4th columns (exp1 and exp2).
>>
>> You say that limma failed on your design matrix, but really it didn't. It
>> correctly removed the redundant column from your design matrix, and gave
>> you a warning to let you know.
>>
>> Best wishes
>> Gordon
>>
>> ---------------- original message ----------------
>> [BioC] Design and not estimable coefficients
>> Dave Canvhet dcanvhet at gmail.com
>> Wed Apr 25 18:08:59 CEST 2012
>>
>> Dear all,
>>
>> I have the following target
>>
>> target
>>
>> filename mutation exp
>> sample1.cel wt 1
>> sample2.cel wt 1
>> sample3.cel mu 1
>> sample4.cel mu 1
>> sample5.cel wt 2
>> sample6.cel mu 2
>> sample7.cel mu 2
>> sample8.cel mu 2
>>
>>
>> I'm intersting in mu vs wt, but taking into account the exp factor (which
>> seems having most impact on signal as reveal by a PCA : PC1 seems to be
>> assoiated to exp)
>>
>> So I've set the following design so as to model the resulting signal as a
>> combination of these two factor:
>>
>> design = matrix(0, ncol=4, nrow = 8)
>> design[which(target$mutation == "wt"),1] = 1
>> design[which(target$mutation == "mu"),2] = 1
>> design[which(target$exp == 1),3] = 1
>> design[which(target$exp == 2),4] = 1
>> colnames(design) = c("wt","mu","exp1","exp2")
>>
>>
>> design
>> wt mu exp1 exp2
>> [1,] 1 0 1 0
>> [2,] 1 0 1 0
>> [3,] 0 1 1 0
>> [4,] 0 1 1 0
>> [5,] 1 0 0 1
>> [6,] 0 1 0 1
>> [7,] 0 1 0 1
>> [8,] 0 1 0 1
>>
>> fit = lmFit(x, design)
>> but it failed (even before doinf any contrast matrix), raising the
>> following error :
>> Coefficients not estimable: exp2
>>
>> Can you please tell me what is wrong is such design ?
>>
>>
>> SO I've used the approach described in Limma User guide :
>> f <- paste(target$mutation,target$**exp,sep="")
>> f <- factor(f)
>> [1] wt1 wt1 mu1 mu1 wt2 mu2 mu2 mu2
>> Levels: mu1 mu2 wt1 wt2
>>
>>
>> design <- model.matrix(~0+f)
>> colnames(design) <- levels(f)
>> design
>>
>> mu1 mu2 wt1 wt2
>> 1 0 0 1 0
>> 2 0 0 1 0
>> 3 1 0 0 0
>> 4 1 0 0 0
>> 5 0 0 0 1
>> 6 0 1 0 0
>> 7 0 1 0 0
>> 8 0 1 0 0
>> attr(,"assign")
>> [1] 1 1 1 1
>> attr(,"contrasts")
>> attr(,"contrasts")$f
>> [1] "contr.treatment"
>>
>> using such a design matrix, can I use correctly the following contrast
>> matrix to get the gene differentially expressed between mutant and wild
>> type ?
>> cont.matrix <- makeContrasts(mu_vs_wt = (mu2+mu1) -
>> (wt2+wt1),levels=design)
>>
>>
>> many thanks for you answer.
>>
>> ==
>> Dave Canvhet
______________________________________________________________________
The information in this email is confidential and intend...{{dropped:4}}
More information about the Bioconductor
mailing list