I have a design given by treatment and categorical stages, each with 5 replicates that I want to model in limma, that looks like this:
treatment stage 1 ctrl A 2 ctrl B 3 ctrl C 4 ctrl D 1 drug A 2 drug B 3 drug C 4 drug D
I want to model the individual effect of treatment, stage, as well as their interaction. To do so I'm considering
design = ~0+ treatment*stage. I then specify the contrasts that I'm interested, namely:
drug - ctrl D - B (drug_D - drug_B) - (ctrl_D - ctrl_B)
so I specify
model_mat = model.matrix( eval( parse ( text = "~0+ treatment + stage" ) ), data=some_df) .
This generates the model matrix that then I want to use to build my contrasts, for the 3 different contrasts above. However, some fail because not all columns are present in
#works makeContrasts(contrasts = 'treatmentdrug - treatmentctrl', levels = model_mat) #works makeContrasts(contrasts = 'stageD - stageB', levels = model_mat) # no columns in model_mat contain treatmentdrg:stage # only columns with treatmentctrl:stage are present # (fails because treatmentdrug:stageD and treatmentdrug:stageB are not present) # makeContrasts( contrasts = '(treatmentdrug:stageD - treatmentdrug:stageB) - (treatmentctrl:stageD - treatmentctrl:stageB)', levels=model_mat)
I understand that some of the columns are not there as they can be given by combinations of the other factors, but this becomes trickier in more complex designs such as the last contrast that I want to make. Is there a way to overcome this issue?
I see that the alternative given the tests I want to make could be to specify
design1=~0 + treatment + stage which would let me model the first 2 contrats, and
design1=~0 + treatment:stage which would let me model the last. But in the case of
design2, then the model does not take into account the individual effect of each
stage, but only their combination, which I don't think it is correct mathematically.
In other words, what I am asking is how to express all contrasts that are not present? (e.g.
From what I understand,
makeContrasts gives everything with reference to the first alphabetical group, so no treatmentdrug will show up for the interactions if I take a different design such as those above. In this case how would you compute `treatmentdrug:stageD - treatmentdrug:stageB) - (treatmentctrl:stageD - treatmentctrl:stageB)
Any suggestions on how to solve these issues? Thanks in advance