Paired Experimental Design
Entering edit mode
Last seen 4 months ago
United Kingdom


Hypothetical question I was thinking about. If I have an experiment that arrayed 2 paired experiments, is it possible to model them together?


Treatments <- factor(c(rep(c("ControlA", "TreatA"), 3),
                       rep(c("ControlB", "TreatB"), 3)))
Pairings   <- factor(paste0("pair",
                            rep(1:6, each=2)))
design     <- model.matrix(~Pairings + Treatments)
"(Intercept)"        "Pairingspair2"      "Pairingspair3"      "Pairingspair4"      "Pairingspair5"      "Pairingspair6"      "TreatmentsControlB" "TreatmentsTreatA"   "TreatmentsTreatB"

If I wanted the Contrast that denotes the difference between ControlB and TreatB, or the difference between ControlA and TreatA, would I need to make a contrast matrix, or is that accessible through one of the columns in the design matrix? 


limma model.matrix • 1.5k views
Entering edit mode
Last seen 2 days ago
United States

In that design, the last three coefficients are comparisons to the baseline, which is TreatmentsControlA (e.g., TreatmentsControlB is actually TreatmentsControlB - TreatmentsControlA), so you already get the comparison of TreatmentsTreatA and TreatmentsControlA, but any other comparisons (TreatmentsTreatB - TreatmentsControlB) have to be done using a contrasts matrix. Since you need a contrasts matrix anyway, you could eliminate the intercept and just make all the comparisons you need.

So yes, you can model the experiments together, and contingent upon the intra-group variability being relatively similar between the different groups, this will increase power. If one of the experiments is more variable than the other, then you may reduce power in the less variable experiment because you may be introducing extra variability. Or you may be correcting for artificially low variability, depending on your interpretation.

Entering edit mode

Well, if I fit a none intercept model (~0 + Pairings + Treatments), then as expected, TreatmentsControlA gets dropped as the baseline. Would it make more sense to use the model ~0 + Treatments + Pairings so that all terms are present for the contrast matrix? Then it'd be a case of:

design     <- model.matrix(~0 + Treatments + Pairings)
cont_mat   <- makeContrasts(TreatmentsTreatA-TreatmentsControlA,
                            levels = colnames(design))
"TreatmentsControlA" "TreatmentsControlB" "TreatmentsTreatA"   "TreatmentsTreatB"   "Pairingspair2"      "Pairingspair3"      "Pairingspair4"      "Pairingspair5"      "Pairingspair6"

That should still utilise the information between the pairs wouldn't it? 

Entering edit mode

In fact, the design matrix isn't of full rank; you need to drop TreatmentsControlB for the coefficients to be estimable. The TreatmentsTreatA term represents the treatment effect in experiment A, while the TreatmentsTreatB term represents the effect in experiment B. You can drop them individually to test for the treatment effect in each experiment, or you can do fancier stuff, e.g., compare them to each other, drop both of them, test the average of them against zero.

Also check out Gordon's comment (A: Correct assumptions of using limma moderated t-test) regarding unequal variances.

Edit: Note that I'm referring to the design matrix in your original post. The design matrix in your comment would require you to drop, e.g., Pairingspair4 to get to full column rank. It's a lot harder to interpret the coefficients in the second design compared to the first, because you don't have a blocking term for each pairing.

Entering edit mode

Cheers Aaron, really useful. I agree that the first model is far easier to interpret after dropping those terms to get full rank. I've seen the post with Gordon's answer before, but useful to go over again! 


Login before adding your answer.

Traffic: 514 users visited in the last hour
Help About
Access RSS

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6