r - 如何强制 model.matrix/model.Matrix() 保持因子水平的原始顺序......？

Question

我有一个简单的假设血统

> 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))

所以，就data.frame而言，我的谱系如下:

> 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

基于 dam 和 sire 变量，我创建了一个 family 字段(family=damxsire，仅当 dam 和 sire 不为零时，NA 在另一种情况下)，即

> 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

从我的family变量中，我想得到一个设计矩阵(Zfam)，即

> 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

结果矩阵是:

> 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

出于某种未知的原因，model.Matrix 对系列级别进行了重新排序，以便系列 1000x7 首先出现。问题是，对于后期分析，我需要根据数据的原始顺序(家庭级别的原始顺序)构建 Zfam 矩阵。

预期输出

> 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 

另一方面，我认为还有一个问题与 R 对字符向量进行排序的方式有关。例如，家庭1000x7从最后一个位置移动到第一个位置(这里开始我的头痛)

> 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>

我还尝试了另一种不切实际的方法。例如，使用 model.matrix 函数(来自 stats 包)，请参阅以下代码

> 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

但是，我得到了与之前相同的结果......

目前，我的解决方法是以 model.Matrix(来自 MatrixModel 或 stats 包)排序的方式对数据进行排序。然而，这种人为的数据排序会给我的其余分析带来问题(这个阶段只是广泛分析的开始)。我确信有更好/更有效的方法来完成这项任务......

任何帮助将不胜感激。

Answer 1

矩阵中列的顺序取决于family列中的因子水平。因此，您可以按照您想要的顺序分配因子水平。

在这种情况下，您需要按照它们出现的顺序，以便您可以使用unique。

library(dplyr)

datafam <- pedigree %>% 
  mutate(family=ifelse((sire==0 | dam==0), NA, paste(dam, sire, sep="x")),
         family = factor(family, levels = unique(family)))

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
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