I have a simple hypothetical pedigree
> dam <- c(0, 0, 0, 0, 2, 4, 5, 6, 9, 1000)
> sire <- c(0, 0, 0, 0, 1, 3, 1, 3, 8, 7)
> ID <- c(1:length(dam))
So, in terms of data.frame, my pedigree is as follow:
> pedigree <- data.frame(ID, dam, sire)
> pedigree
ID dam sire
1 1 0 0
2 2 0 0
3 3 0 0
4 4 0 0
5 5 2 1
6 6 4 3
7 7 5 1
8 8 6 3
9 9 9 8
10 10 1000 7
Based on dam and sire variables, I have created a family field (family=damxsire, only when both dam and sire are different to zero, NA in another case), i.e.
> datafam <- pedigree %>%
mutate(family=ifelse((sire==0 | dam==0), NA, as.vector(paste(dam, sire, sep="x")))) %>%
mutate_at(vars(family), as.factor)
> datafam
ID dam sire family
1 1 0 0 <NA>
2 2 0 0 <NA>
3 3 0 0 <NA>
4 4 0 0 <NA>
5 5 2 1 2x1
6 6 4 3 4x3
7 7 5 1 5x1
8 8 6 3 6x3
9 9 9 8 9x8
10 10 1000 7 1000x7
and from my family variable, I want to get a design matrix (Zfam), that is
> form1 <- formula(~ family -1)
> termsf1 <- terms(form1, keep.order = TRUE)
> mf1 <- model.frame(termsf1, data=datafam, na.action= na.pass)
> Zfam <- as.matrix(MatrixModels::model.Matrix(form1, mf1, sparse=FALSE))
> Zfam[is.na(Zfam)] <- 0 # replaces any missing values in Z by zeros
The resulting matrix is:
> Zfam
family1000x7 family2x1 family4x3 family5x1 family6x3 family9x8
1 0 0 0 0 0 0
2 0 0 0 0 0 0
3 0 0 0 0 0 0
4 0 0 0 0 0 0
5 0 1 0 0 0 0
6 0 0 1 0 0 0
7 0 0 0 1 0 0
8 0 0 0 0 1 0
9 0 0 0 0 0 1
10 1 0 0 0 0 0
For some unknown reason, model.Matrix reorders the family levels, so that the family 1000x7 appears first. The problem is that for a post-analysis, I need the Zfam matrix to be built based on the original order of the data (original order of family levels).
Expected output
> Zfam
family2x1 family4x3 family5x1 family6x3 family9x8 family1000x7
1 0 0 0 0 0 0
2 0 0 0 0 0 0
3 0 0 0 0 0 0
4 0 0 0 0 0 0
5 1 0 0 0 0 0
6 0 1 0 0 0 0
7 0 0 1 0 0 0
8 0 0 0 1 0 0
9 0 0 0 0 1 0
10 0 0 0 0 0 1
On the other hand, I think that an additional issue is concerned to the way to R sorts character vectors. For example, family 1000x7 was moved from the last to the first position (and here begins my headache)
> datafam[with(datafam, order(family)), ]
ID dam sire family
10 10 1000 7 1000x7
5 5 2 1 2x1
6 6 4 3 4x3
7 7 5 1 5x1
8 8 6 3 6x3
9 9 9 8 9x8
1 1 0 0 <NA>
2 2 0 0 <NA>
3 3 0 0 <NA>
4 4 0 0 <NA>
I also tried with another impractical approach. For example, using model.matrix function (from stats package), see the following code
> form1 <- formula(~ family -1)
> termsf1 <- terms(form1, keep.order = TRUE)
> mf1 <- model.frame(termsf1, data=datafam, na.action= na.pass)
> Zfam <- as.matrix(stats::model.matrix(form1, mf1, sparse=FALSE))
> Zfam[is.na(Zfam)] <- 0 # replaces any missing values in Z by zeros
However, I got the same previous result…
For now, my workaround has been to sort the data in the same way it is sorted by the model.Matrix (from MatrixModel or stats packages). However, this artificial ordering of the data will give me problems for the rest of the analysis (this stage is just the beginning of an extensive analysis). I'm sure there is a better/efficient way to do this task ...
Any help would be greatly appreciated.