2.5 years ago by
Cambridge, United Kingdom
For your specific contrast of interest, I would set up the design matrix to block on the sibling factor:
data$SibShip <- factor(data$SibShip) # Make sure it is a factor
design <- model.matrix(~0 + SibShip + Treatment:Schedule, data)
design <- design[,-grep("TreatmentC", colnames(design))] # Get to full rank
colnames(design) <- sub(":", "_", colnames(design)) # Rename for convenience
Now, if we have a look at the column names of the design matrix, we get:
 "SibShip1" "SibShip2" "SibShip3"
 "SibShip4" "SibShip5" "SibShip6"
 "TreatmentT_ScheduleA" "TreatmentT_ScheduleB"
The first 6 coefficients are the sibling-specific effects, which are not of interest to us. The last two coefficients represent the effect (i.e., log-fold change in expression) of treatment over control in schedules A and B, respectively. If you want to find differences in the treatment effect between schedules, you can set up a contrast like this:
con <- makeContrasts(TreatmentT_ScheduleA - TreatmentT_ScheduleB, levels=design)
Obviously, you can also drop either coefficient separately to determine if there is a non-zero treatment effect for the corresponding schedule. This is often a good idea to help in interpreting the log-fold change from the above contrast, e.g., is a positive log-fold change due to an increase in the treatment effect for schedule A or a decrease in the effect for schedule B?
You could also use
duplicateCorrelation here. However, this is unnecessary as the calculation of the treatment effect is done within the blocking factor, such that the sibling effect cancels out naturally. If you wanted to compare, say, treatment samples directly between schedules, then you would have to block on sibling in
duplicateCorrelation instead of in the design matrix.
modified 2.5 years ago
2.5 years ago by
Aaron Lun • 25k