如何强制 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 对字符向量进行排序的方式。例如，family 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