I have a proteomics experiment I want to analyse with limma. It is a crossover design. Each subject received a treatment over the two periods, separated by a washout interval. Some subjects received the drug in the first period and placebo in the second, some the other way round. Samples were collected at the start and end of each period, so there are 4 samples from each subject. For example, subject 1 received drug during the first period and placebo during the second period and the four samples from this subject represent the effect of the baseline, drug, baseline, and placebo.

The purpose of the experiment is to see the effect of the drug, as compared to the placebo.

This is the treatment allocation for the first 4 subjects (there are more in the actual dataset):

> meta
   subject period course order    treatment
   <fct>   <fct>  <fct>  <chr>    <chr>    
 1 1       P1     start  DrugPlac baseline 
 2 1       P1     end    DrugPlac drug     
 3 1       P2     start  DrugPlac baseline 
 4 1       P2     end    DrugPlac placebo  
 5 2       P1     start  PlacDrug baseline 
 6 2       P1     end    PlacDrug placebo  
 7 2       P2     start  PlacDrug baseline 
 8 2       P2     end    PlacDrug drug     
 9 3       P1     start  DrugPlac baseline 
10 3       P1     end    DrugPlac drug     
11 3       P2     start  DrugPlac baseline 
12 3       P2     end    DrugPlac placebo  
13 4       P1     start  PlacDrug baseline 
14 4       P1     end    PlacDrug placebo  
15 4       P2     start  PlacDrug baseline 
16 4       P2     end    PlacDrug drug

Following this answer from Aaron Lun I would do the following limma approach:

design_mat <- model.matrix(~ 0 + subject + treatment:order, data = meta)
design_mat <- design_mat[,!grepl("baseline", colnames(design_mat))] # get to full rank.
colnames(design_mat) <- colnames(design_mat) |>
    str_replace(":", "_")

colnames(design_mat)
[1] "subject1"                       "subject2"                       "subject3"                      
[4] "subject4"                       "subject5"                       "treatmentdrug_orderDrugPlac"   
[7] "treatmentplacebo_orderDrugPlac" "treatmentdrug_orderPlacDrug"    "treatmentplacebo_orderPlacDrug"

con <- makeContrasts(
  treatmentdrug_orderDrugPlac - treatmentdrug_orderPlacDrug,        # effect of order
  treatmentplacebo_orderDrugPlac - treatmentplacebo_orderPlacDrug,  # effect of order
  treatmentdrug_orderDrugPlac - treatmentplacebo_orderDrugPlac,     # ?
  treatmentdrug_orderPlacDrug - treatmentplacebo_orderPlacDrug,     # ?
  levels = design_mat
)

fit <- tab |>
    lmFit(design_mat) |> 
    contrasts.fit(contrasts = con) |> 
    eBayes()

If I understand this correctly, the first two contrasts would tell me if there is any effect of order (there should not be any). However, after that, I'm simply interested in the effect of the drug, while there are two contrast, each under different order.

For comparison, this answer from Aaron Lun would suggest a different approach:

design_mat <- model.matrix(~ 0 + subject:treatment + treatment:course, data = meta)
colnames(design_mat)
 [1] "subject1:treatmentbaseline"  "subject2:treatmentbaseline"  "subject3:treatmentbaseline" 
 [4] "subject4:treatmentbaseline"  "subject5:treatmentbaseline"  "subject1:treatmentdrug"     
 [7] "subject2:treatmentdrug"      "subject3:treatmentdrug"      "subject4:treatmentdrug"     
[10] "subject5:treatmentdrug"      "subject1:treatmentplacebo"   "subject2:treatmentplacebo"  
[13] "subject3:treatmentplacebo"   "subject4:treatmentplacebo"   "subject5:treatmentplacebo"  
[16] "treatmentbaseline:courseend" "treatmentdrug:courseend"     "treatmentplacebo:courseend"

This is even more difficult to understand.

I'm confused here. In "normal" experiments I only use interaction terms to see if there is an emergent effect of two different variables, but these variables are always included in the model. The approaches suggested above use interaction terms only, which is new to me. Yes, I know I should check for the effect of the placebo-drug order. But how do I do the drug-placebo comparison? The simplest approach would be:

design_mat <- model.matrix(~ subject + treatment, data = meta)
colnames(design_mat)
[1] "(Intercept)"      "subject2"         "subject3"         "subject4"         "subject5"         "treatmentdrug"   
[7] "treatmentplacebo"

which takes subjects into account and gives me directly the effect of the drug and placebo over the baseline. But I feel like I'm missing some of the experiment's design complexity here.

Any other ideas how to approach this design?