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 model_mat
:
#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 treatment
and 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. treatmentdrug:stageB
).
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
Hello, sorry for replying on an old thread. Did you figure out a solution to your problem? I'm encountering the same with my analysis.