在R中将系数名称转换为公式

9

当使用具有因子的公式时,拟合模型会将系数命名为XY,其中X是因子的名称,Y是其特定水平。我希望能够根据这些系数的名称创建一个公式。

原因:如果我将套索拟合到稀疏设计矩阵(如下所示),我希望创建一个仅包含非零系数项的新公式对象。

require("MatrixModels")
require("glmnet")
set.seed(1)
n <- 200
Z <- data.frame(letter=factor(sample(letters,n,replace=T),letters),
                x=sample(1:20,200,replace=T))
f <- ~ letter + x:letter + I(x>5):letter
X <- sparse.model.matrix(f, Z)
beta <- matrix(rnorm(dim(X)[2],0,5),dim(X)[2],1)
y <- X %*% beta + rnorm(n)

myfit <- glmnet(X,as.vector(y),lambda=.05)
fnew <- rownames(myfit$beta)[which(myfit$beta != 0)]
 [1] "letterb"              "letterc"              "lettere"             
 [4] "letterf"              "letterg"              "letterh"             
 [7] "letterj"              "letterm"              "lettern"             
[10] "lettero"              "letterp"              "letterr"             
[13] "letters"              "lettert"              "letteru"             
[16] "letterw"              "lettery"              "letterz"             
[19] "lettera:x"            "letterb:x"            "letterc:x"           
[22] "letterd:x"            "lettere:x"            "letterf:x"           
[25] "letterg:x"            "letterh:x"            "letteri:x"           
[28] "letterj:x"            "letterk:x"            "letterl:x"           
[31] "letterm:x"            "lettern:x"            "lettero:x"           
[34] "letterp:x"            "letterq:x"            "letterr:x"           
[37] "letters:x"            "lettert:x"            "letteru:x"           
[40] "letterv:x"            "letterw:x"            "letterx:x"           
[43] "lettery:x"            "letterz:x"            "letterb:I(x > 5)TRUE"
[46] "letterc:I(x > 5)TRUE" "letterd:I(x > 5)TRUE" "lettere:I(x > 5)TRUE"
[49] "letteri:I(x > 5)TRUE" "letterj:I(x > 5)TRUE" "letterl:I(x > 5)TRUE"
[52] "letterm:I(x > 5)TRUE" "letterp:I(x > 5)TRUE" "letterq:I(x > 5)TRUE"
[55] "letterr:I(x > 5)TRUE" "letteru:I(x > 5)TRUE" "letterv:I(x > 5)TRUE"
[58] "letterx:I(x > 5)TRUE" "lettery:I(x > 5)TRUE" "letterz:I(x > 5)TRUE"

From this I would like to have a formula

~ I(letter=="d") + I(letter=="e") + ...(etc)

我查看了formula()和all.vars(),但都没有找到合适的解决方法。由于出现了不同类型的术语,编写解析函数有点麻烦。例如,当x是数字值而letter是因子时,出现了x:letter这种情况,或者I(x>5):letter这种令人讨厌的情况。

那么,我是否不知道有一些函数可以在公式和其字符表示之间进行转换?


我在 R 中不认识那个公式。 - Gavin Simpson
1
也许我误解了,但您似乎没有完全理解 R 的模型公式。您的公式中没有包括 XY 部分,您只包括了 X,然后 model.matrix()model.frame() 会将 X 的级别扩展到相关的模型矩阵列,即 XY。 - Gavin Simpson
你能解释一下为什么你想要这个公式吗?最终用途是什么? - Gavin Simpson
我有一个由模型f创建的具有200k列的设计矩阵。在拟合套索回归后,我得到了一个带有5k个项的模型g。我想使用g而不是f构建新的设计矩阵。让model.matrix()和model.frame()执行它们的功能后,我想获得与模型矩阵列的某个子集(展开级别后)对应的公式。 - Christopher DuBois
2个回答

4
当我运行代码时,由于没有指定set.seed(),我得到了一些略微不同的结果。我使用"letter_"作为方便的分割标记,而不是使用变量名"letter"。
> fnew <- rownames(myfit$beta)[which(myfit$beta != 0)]

> fnew
 [1] "letter_c" "letter_d" "letter_e" "letter_f" "letter_h" "letter_k" "letter_l"
 [8] "letter_o" "letter_q" "letter_r" "letter_s" "letter_t" "letter_u" "letter_v"
[15] "letter_w"

然后将分割后的内容打包到一个字符矩阵中:
> fnewmtx <- cbind( lapply(sapply(fnew, strsplit, split="_"), "[[", 2),
+ lapply(sapply(fnew, strsplit, split="_"), "[[", 1))

fnewmtx [,1] [,2]
letter_c "c" "letter" letter_d "d" "letter" letter_e "e" "letter" letter_f "f" "letter" 此处省略其余内容

将paste函数的输出结果包装在as.formula()中,这是如何“在公式和其字符表示之间进行转换”的一半答案。另一半是as.character()。

form <- as.formula( paste("~", 
             paste( 
               paste(" I(", fnewmtx[,2], "_ ==", "'",fnewmtx[,1],"') ", sep="") , 
             sep="", collapse="+")
                 ) 
           )  # edit: needed to add back the underscore

现在的输出是一个合适的类对象:

> class(form)
[1] "formula"
> form
~I(letter_ == "c") + I(letter_ == "d") + I(letter_ == "e") + 
    I(letter_ == "f") + I(letter_ == "h") + I(letter_ == "k") + 
    I(letter_ == "l") + I(letter_ == "o") + I(letter_ == "q") + 
    I(letter_ == "r") + I(letter_ == "s") + I(letter_ == "t") + 
    I(letter_ == "u") + I(letter_ == "v") + I(letter_ == "w")

我觉得有趣的是,as.formula 转换将字母周围的单引号变成了双引号。

编辑:现在问题增加了一两个维度,我的建议是跳过重新创建公式的步骤。请注意,myfit$beta 的行名与 X 的列名完全相同,因此可以使用非零行名作为索引,在 X 矩阵中选择列:

> str(X[ , which( colnames(X) %in% rownames(myfit$beta)[which(myfit$beta != 0)] )] )
Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
  ..@ i       : int [1:429] 9 54 91 157 166 37 55 68 117 131 ...
  ..@ p       : int [1:61] 0 5 13 20 28 36 42 50 60 68 ...
  ..@ Dim     : int [1:2] 200 60
  ..@ Dimnames:List of 2
  .. ..$ : chr [1:200] "1" "2" "3" "4" ...
  .. ..$ : chr [1:60] "letter_b" "letter_c" "letter_e" "letter_f" ...
  ..@ x       : num [1:429] 1 1 1 1 1 1 1 1 1 1 ...
  ..@ factors : list()

> myfit2 <- glmnet(X[ , which( colnames(X) %in% rownames(myfit$beta)[which(myfit$beta != 0)] )] ,as.vector(y),lambda=.05)
> myfit2

Call:  glmnet(x = X[, which(colnames(X) %in% rownames(myfit$beta)[
                                           which(myfit$beta != 0)])], 
              y = as.vector(y), lambda = 0.05) 

     Df   %Dev Lambda
[1,] 60 0.9996   0.05

谢谢DWin!这很接近,但并不是一个通用的解决方案。我已经更新了我的示例。是否有一个快速修复可以扩展您的解决方案? - Christopher DuBois
另外,在创建fnewmtx时,strsplit不应该使用split =“”而应该使用split =“_”吧? - Christopher DuBois
非常正确。我会进行编辑。我正在尝试处理交互公式的复杂性,但不要期望我能很快回复......更新:尝试编辑,但下划线在原始代码段中,未在SO页面上显示。重新编码,现在可能更好了。 - IRTFM
1
@Christopher 不确定你上面提到的X是哪个。DWin的解决方案是对你的稀疏模型矩阵X进行索引。从你的问题中,我得到了一个明显的印象,你可以形成该矩阵以便使用glmnet()拟合模型。DWin的子集生成一个新的稀疏矩阵,其中仅包含系数非零的列。你想要从这个减少的稀疏设计矩阵创建完整的设计矩阵(而不是稀疏版本)吗?如果是这样,请尝试:newX <- as.matrix(X[ , which( colnames(X) %in% rownames(myfit$beta)[which(myfit$beta != 0)] )])newX有60列。 - Gavin Simpson
1
抱歉没有表达清楚。假设我创建了我的第一个稀疏模型矩阵X_train并使用glmnet()进行拟合。(对于我的特定问题,这需要5-15分钟。)现在我想使用新数据X_test创建一个新的稀疏矩阵,但只包含系数非零的列。出于计算原因,我想避免为具有所有列的X_test构建稀疏矩阵,因此上面的子集解决方案对我不起作用(除非我漏掉了什么)。 - Christopher DuBois
显示剩余3条评论

2
克里斯托弗,经过对sparse.model.matrix等的考虑和检查,你所要求的似乎有些复杂。你没有解释为什么不想为X_test形成完整的稀疏模型矩阵,因此很难提供除下面两个选项之外的建议。
如果X_test中有大量观测值,因此不希望为计算原因生成完整的稀疏矩阵以在predict()中使用,则将X_test拆分为两个或多个样本块,并逐个形成每个块的稀疏模型矩阵,在使用后丢弃即可更加迅速。
否则,您需要详细研究Matrix软件包中的代码。从sparse.model.matrix开始,并注意它然后调用Matrix:::model.spmatrix并在该函数中定位对Matrix:::fac2Sparse的调用。您可能需要利用这些函数的代码,但使用修改后的fac2Sparse来实现您想要实现的内容。
很抱歉我无法提供现成的脚本来完成此任务,但这是一个重要的编码任务。如果您选择这条路,请查看Matrix软件包中的“稀疏模型矩阵”vignette,并获取软件包源代码(来自CRAN),以查看我提到的函数是否在源代码中有更好的文档记录(例如,fac2Sparse没有Rd文件)。您还可以向Matrix的作者(Martin Maechler和Doug Bates)寻求建议,尽管请注意,这两位都有特别繁重的教学任务。
祝你好运!

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接