Question: Design matrix not of full rank.
0
gravatar for rmf
2.7 years ago by
rmf560
rmf560 wrote:

I have a dataset like so

df <- data.frame(condition=rep(c("c","l"),each=12),family=rep(c(20,21,27,30),each=3))
df
   condition family
1          c     20
2          c     20
3          c     20
4          c     21
5          c     21
6          c     21
7          c     27
8          c     27
9          c     27
10         c     30
11         c     30
12         c     30
13         l     20
14         l     20
15         l     20
16         l     21
17         l     21
18         l     21
19         l     27
20         l     27
21         l     27
22         l     30
23         l     30
24         l     30

I would like to compare condition c and l within families. I created a model like so:

design.mat <- model.matrix(~condition+family,df)

creating this model matrix

   (Intercept) conditionl family20 family21 family27 family30
1            1          0        1        0        0        0
2            1          0        1        0        0        0
3            1          0        1        0        0        0
4            1          0        0        1        0        0
5            1          0        0        1        0        0
6            1          0        0        1        0        0
7            1          0        0        0        1        0
8            1          0        0        0        1        0
9            1          0        0        0        1        0
10           1          0        0        0        0        1
11           1          0        0        0        0        1
12           1          0        0        0        0        1
13           1          1        1        0        0        0
14           1          1        1        0        0        0
15           1          1        1        0        0        0
16           1          1        0        1        0        0
17           1          1        0        1        0        0
18           1          1        0        1        0        0
19           1          1        0        0        1        0
20           1          1        0        0        1        0
21           1          1        0        0        1        0
22           1          1        0        0        0        1
23           1          1        0        0        0        1
24           1          1        0        0        0        1
attr(,"assign")
[1] 0 1 2 2 2 2
attr(,"contrasts")
attr(,"contrasts")$condition
[1] "contr.treatment"

attr(,"contrasts")$family
[1] "contr.treatment"

I then ran

eb <- DGEList(counts)
ebg <- estimateGLMCommonDisp(eb,design.mat)

but I get this error:

Error in glmFit.default(y, design = design, dispersion = dispersion, offset = offset,  : 
  Design matrix not of full rank.  The following coefficients not estimable:
 family30

I also tried

design.mat <- model.matrix(~0+condition+family,df)

as suggested in another post, but it still does not work. Any ideas on how this can be tackled is highly appreciated. Thanks.

edger rna-seq R • 2.6k views
ADD COMMENTlink written 2.7 years ago by rmf560

Just drop the column 'family 20' from the design matrix. In this setting condition:c / family:20 is your baseline class

ADD REPLYlink modified 2.7 years ago • written 2.7 years ago by russhh4.2k

Yes. That did fix the issue. Thanks a lot. Would you like to add that as an answer?

ADD REPLYlink written 2.7 years ago by rmf560

Not really, it's pretty elementary. Can you see why the matrix you posted has rank 5?

ADD REPLYlink written 2.7 years ago by russhh4.2k

When calling design.mat <- model.matrix(~condition+family,df), the resulting object has no conditionc column. It was automatically removed as the first level of the condition factor. But then why was family20 not removed in a similar way?

ADD REPLYlink written 2.7 years ago by rmf560
1

It's a trap ... to keep us all in business. Plus I don't think model.matrix could really be automated for all the possible designs in the world

ADD REPLYlink written 2.7 years ago by russhh4.2k

Can I just check that your real code uses

family=factor(rep(c(20,21,27,30),each=3))

rather than family=rep(c(20,21,27,30),each=3)

ADD REPLYlink written 2.7 years ago by russhh4.2k

hmmm... that might be it. I just tried it again and family20 does seem to be removed atutomatically.

ADD REPLYlink written 2.7 years ago by rmf560
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1377 users visited in the last hour