Hello,

This is yet again a question about limma design matrices. I would like to
perform analysis on single channel arrays (Affy). The experimental design
is: paired samples before and after treatment, and I would also like to take
into account batch effect.

I am having problems with the batch effect, by reading through mails in this
mailing list, I've gotten the impression that you can include batch to your
design matrix, but when I try to do that, I get an error message. Here is a
simulated example of what I am trying to do:

Create example date, patient (5 patients, 2 samples), treatment (A or B),
and batch (three batches, 1-3)

> patient <- factor(rep(1:5,2))
> patient
 [1] 1 2 3 4 5 1 2 3 4 5
Levels: 1 2 3 4 5
> treatment <- factor(c(rep("A",5), rep("B",5)))
> treatment
 [1] A A A A A B B B B B
Levels: A B
> batch <- factor(c(rep(1,3), rep(2,3), rep(3,4)))
> batch
 [1] 1 1 1 2 2 2 3 3 3 3
Levels: 1 2 3

If I would like to perform a basic analysis between the paired samples and
not taking into account batch effect, I could just do:

> design <- model.matrix(~patient+treatment)
> design
   (Intercept) patient2 patient3 patient4 patient5 treatmentB
1            1        0        0        0        0          0
2            1        1        0        0        0          0
3            1        0        1        0        0          0
4            1        0        0        1        0          0
5            1        0        0        0        1          0
6            1        0        0        0        0          1
7            1        1        0        0        0          1
8            1        0        1        0        0          1
9            1        0        0        1        0          1
10           1        0        0        0        1          1

And then use treatmentB to obtain the differences between treatments A and
B.

But, how could I include batch effect to the model? My inital thought was:

> design <- model.matrix(~patient+treatment+batch)
> design
   (Intercept) patient2 patient3 patient4 patient5 treatmentB batch2 batch3
1            1        0        0        0        0          0      0      0
2            1        1        0        0        0          0      0      0
3            1        0        1        0        0          0      0      0
4            1        0        0        1        0          0      1      0
5            1        0        0        0        1          0      1      0
6            1        0        0        0        0          1      1      0
7            1        1        0        0        0          1      0      1
8            1        0        1        0        0          1      0      1
9            1        0        0        1        0          1      0      1
10           1        0        0        0        1          1      0      1

But when I do that with actual data, I get the following error message:

Coefficients not estimable: batch2 batch3
Warning message:
Partial NA coefficients for 10163 probe(s)

Any ideas what I am doing wrong, and more importantly, what I should do to
take that batch effect into account?

Thanks,
John

	[[alternative HTML version deleted]]

