Question: Paired design in limma
gravatar for rasmus.agren
7 months ago by
rasmus.agren0 wrote:


I wonder if someone could double check that this design is fine. I have a design with paired samples and I'd like to test for the effect of the treatment and the interaction between treatment and group.

This is the experiment: Overweight and normal patients are investigated for their response to a treatment. Samples are taken from each patient before and after the treatment. I'm interested in identifying genes which respond to the treatment differently in obese vs normal. This is to be tested with a fixed effects model, not using the duplicateCorrelation option.

group <- paste0(obese_normal, ".", before_after)
design <- design.matrix(~ 0 + subject + group)

This design matrix isn't full rank (1 dependent column). I therefore remove the "obese.before" column. "normal.before" was the first level so it's already removed. After that it's full rank. Am I correct in that I can now test for the interaction effect with:

contrast <- makeContrasts(obese.after - normal.after, levels=design)


limma paired • 165 views
ADD COMMENTlink modified 7 months ago by Aaron Lun25k • written 7 months ago by rasmus.agren0
Answer: Paired design in limma
gravatar for Aaron Lun
7 months ago by
Aaron Lun25k
Cambridge, United Kingdom
Aaron Lun25k wrote:

It's a lot easier to explain if you make up a little example for us:

subject <- factor(rep(1:4, each=2))
obese_normal <- rep(c("obese", "normal"), each=4)
before_after <- rep(c("before", "after"), 4)

Your design matrix is pretty much the same as what I would recommend, though this requires a bit of care with the set-up. There's several ways to get the right result, but this is what I would do:

# Renaming to make things simpler.
D <- obese_normal
T <- before_after

design <- model.matrix(~0 + subject + D:T)
design <- design[,!grepl("before", colnames(design))]
colnames(design) <- make.names(colnames(design))

The first four columns are the patient blocking factors, which are not scientifically interesting. The last two columns are the average treatment effects in normal and obese patients respectively. These can be dropped individually to determine the treatment effect in each condition, or compared to each other to determine the differences in treatment effects, i.e., the interaction between obesity and treatment, which is what you wanted.

I've used D:T rather than group because group will automatically drop the first level, which in this case is "normal.after" (it depends on how you named the levels). However, I want to retain the "after" terms, as it makes more sense to think of the treatment effect as after - before. This requires me to get all levels in the output from model.matrix in order to explicitly drop the terms I don't want.

ADD COMMENTlink written 7 months ago by Aaron Lun25k

Great, thanks! I didn't know that D:T would retain all levels. That's useful so you don't have to fiddle with factor level order.

ADD REPLYlink written 7 months ago by rasmus.agren0
Please log in to add an answer.


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