For a particular tool I need a model matrix which allows me to build desired contrasts and which is a full rank, i.e. no columns are linear combinations of other columns.
The experimental design involves four factors:
- type (T) with two levels, A and B
- group (G) with two levels, C and T
- time point (TP) with three levels
- subject id (SID)
Both A and B samples were taken from each subject. Subjects either belong to group C or T (control vs treatment). Multiple samples were taken from each subject at different time points.
The comparisons I desire to make are always within a type (no comparisons between types). I want to test interaction between time points and groups, so for example (T.TP1-C.TP1)-(T.TP0-C.TP0)
.
The only problem is that the model matrix must be full rank, and I don't know how to achieve it.
Here is a mock example:
mock <- data.frame(
ID=paste0("ID", 1:16),
type=rep(c("A", "B"), each=8),
treatment=rep(c("C", "T"), each=4),
tp=c("T1", "T2"),
PID=rep(paste0("P.", 1:8), each=2))
which gives
ID type treatment tp PID
1 ID1 A C T1 P.1
2 ID2 A C T2 P.1
3 ID3 A C T1 P.2
4 ID4 A C T2 P.2
5 ID5 A T T1 P.3
6 ID6 A T T2 P.3
7 ID7 A T T1 P.4
8 ID8 A T T2 P.4
9 ID9 B C T1 P.5
10 ID10 B C T2 P.5
11 ID11 B C T1 P.6
12 ID12 B C T2 P.6
13 ID13 B T T1 P.7
14 ID14 B T T2 P.7
15 ID15 B T T1 P.8
16 ID16 B T T2 P.8
Normally, without the repeated measures, I would do something like
mock$ttt <- with(mock, paste(type, treatment, tp, sep="_"))
mm <- model.matrix(~ 0 + ttt, mock)
...and then define contrasts (B_T_T2-B_C_T2)-(B_T_T1-B_C_T1)
to test an interaction between time points and treatment within type B.
However, I am at loss how to do it with the repeated measures. I tried the following:
mock$type_pid <- paste0(mock$type, "_", mock$PID)
mm <- model.matrix(~ 0 + type_pid + type:tp:treatment, mock)
I get a matrix which is not fully ranked, however I have the coefficients I need for my contrasts. How can I get a fully ranked matrix with the necessary coefficients?
Please note that I am not trying to fit a mixed random model (despite the repeated measures), because my particular setup does not allow for this.