R生成所有可能的交互变量
R generate all possible interaction variables
我有一个带有变量的数据框,比如 a,b,c,d
dat <- data.frame(a=runif(1e5), b=runif(1e5), c=runif(1e5), d=runif(1e5))
并希望在每列之间生成所有可能的双向交互项,即:ab、ac、ad、bc、bd、cd。实际上我的数据框有超过 100 列,所以我不能手动编码。最有效的方法是什么(注意我不想要 ab 和 ba)?
您打算如何处理所有这些交互项?有多种选择,最好取决于您要执行的操作。
如果您想将交互传递给 lm
或 aov
之类的建模函数,则非常简单,只需使用 .^2
语法:
fit <- lm( y ~ .^2, data=mydf )
以上将调用 lm
并告诉它拟合 mydf
中变量的所有主效应和所有 2 种交互作用,不包括 y
.
如果出于某种原因你真的想计算所有的相互作用,那么你可以使用 model.matrix
:
tmp <- model.matrix( ~.^2, data=iris)
这将包括一个用于截距的列和一个用于主效应的列,但如果您不需要它们,您可以删除它们。
如果您需要与建模不同的东西,那么您可以使用 combn
功能,正如@akrun 在评论中提到的那样。
假设预期输出是列名的组合(根据注释应该是a_b
、a_c
等),我们可以在列名上使用combn
数据集并将 m
指定为 2.
combn(colnames(dat), 2, FUN=paste, collapse='_')
#[1] "a_b" "a_c" "a_d" "b_c" "b_d" "c_d"
如果我们需要乘以 'dat' 中的列组合,我们使用 combn
列名输出的每个元素对数据集进行子集化(dat[,x[1]]
,dat[,x[2]]
), 相乘(*
), 转换为'data.frame' (data.frame(
), 通过paste
列的组合设置列名(setNames
)名字。我们在 list
和 cbind
列表元素中创建输出 do.call(cbind
.
do.call(cbind, combn(colnames(dat), 2, FUN= function(x)
list(setNames(data.frame(dat[,x[1]]*dat[,x[2]]),
paste(x, collapse="_")) )))
# a_b a_c a_d b_c b_d c_d
#1 0.26929788 0.17697473 0.26453066 0.55676619 0.83221898 0.54691008
#2 0.06291005 0.08337501 0.04455453 0.10370775 0.05542008 0.07344851
#3 0.53789990 0.47301970 0.03112880 0.51305076 0.03376319 0.02969076
#4 0.41596384 0.34920860 0.25992717 0.53948322 0.40155468 0.33711187
#5 0.16878584 0.21232357 0.09196025 0.08162171 0.03535148 0.04447027
基准
set.seed(494)
dat <- data.frame(a=runif(1e6), b=runif(1e6), c=runif(1e6), d=runif(1e6))
greg <- function()model.matrix( ~.^2, data=dat)
akrun <- function() {do.call(cbind, combn(colnames(dat), 2, FUN= function(x)
list(setNames(data.frame(dat[,x[1]]*dat[,x[2]]),
paste(x, collapse="_")) )))}
system.time(greg())
# user system elapsed
# 1.159 0.024 1.182
system.time(akrun())
# user system elapsed
# 0.013 0.000 0.013
library(microbenchmark)
microbenchmark(greg(), akrun(), times=20L, unit='relative')
# Unit: relative
# expr min lq mean median uq max neval cld
# greg() 39.63122 38.53662 10.23198 18.81274 6.568741 4.642702 20 b
# akrun() 1.00000 1.00000 1.00000 1.00000 1.000000 1.000000 20 a
注意:基准测试因列数和行数而异。在这里,我使用的是 OP post.
中显示的列数
数据
set.seed(24)
dat <- data.frame(a=runif(5), b=runif(5), c=runif(5), d=runif(5))
由于 model.matrix
抱怨只有一个级别的因素,您可能想要使用 stats::terms
labels(terms(~.^2, data = iris[, 1:3]))
# [1] "Sepal.Length" "Sepal.Width" "Petal.Length"
# [4] "Sepal.Length:Sepal.Width" "Sepal.Length:Petal.Length" "Sepal.Width:Petal.Length"
我有一个带有变量的数据框,比如 a,b,c,d
dat <- data.frame(a=runif(1e5), b=runif(1e5), c=runif(1e5), d=runif(1e5))
并希望在每列之间生成所有可能的双向交互项,即:ab、ac、ad、bc、bd、cd。实际上我的数据框有超过 100 列,所以我不能手动编码。最有效的方法是什么(注意我不想要 ab 和 ba)?
您打算如何处理所有这些交互项?有多种选择,最好取决于您要执行的操作。
如果您想将交互传递给 lm
或 aov
之类的建模函数,则非常简单,只需使用 .^2
语法:
fit <- lm( y ~ .^2, data=mydf )
以上将调用 lm
并告诉它拟合 mydf
中变量的所有主效应和所有 2 种交互作用,不包括 y
.
如果出于某种原因你真的想计算所有的相互作用,那么你可以使用 model.matrix
:
tmp <- model.matrix( ~.^2, data=iris)
这将包括一个用于截距的列和一个用于主效应的列,但如果您不需要它们,您可以删除它们。
如果您需要与建模不同的东西,那么您可以使用 combn
功能,正如@akrun 在评论中提到的那样。
假设预期输出是列名的组合(根据注释应该是a_b
、a_c
等),我们可以在列名上使用combn
数据集并将 m
指定为 2.
combn(colnames(dat), 2, FUN=paste, collapse='_')
#[1] "a_b" "a_c" "a_d" "b_c" "b_d" "c_d"
如果我们需要乘以 'dat' 中的列组合,我们使用 combn
列名输出的每个元素对数据集进行子集化(dat[,x[1]]
,dat[,x[2]]
), 相乘(*
), 转换为'data.frame' (data.frame(
), 通过paste
列的组合设置列名(setNames
)名字。我们在 list
和 cbind
列表元素中创建输出 do.call(cbind
.
do.call(cbind, combn(colnames(dat), 2, FUN= function(x)
list(setNames(data.frame(dat[,x[1]]*dat[,x[2]]),
paste(x, collapse="_")) )))
# a_b a_c a_d b_c b_d c_d
#1 0.26929788 0.17697473 0.26453066 0.55676619 0.83221898 0.54691008
#2 0.06291005 0.08337501 0.04455453 0.10370775 0.05542008 0.07344851
#3 0.53789990 0.47301970 0.03112880 0.51305076 0.03376319 0.02969076
#4 0.41596384 0.34920860 0.25992717 0.53948322 0.40155468 0.33711187
#5 0.16878584 0.21232357 0.09196025 0.08162171 0.03535148 0.04447027
基准
set.seed(494)
dat <- data.frame(a=runif(1e6), b=runif(1e6), c=runif(1e6), d=runif(1e6))
greg <- function()model.matrix( ~.^2, data=dat)
akrun <- function() {do.call(cbind, combn(colnames(dat), 2, FUN= function(x)
list(setNames(data.frame(dat[,x[1]]*dat[,x[2]]),
paste(x, collapse="_")) )))}
system.time(greg())
# user system elapsed
# 1.159 0.024 1.182
system.time(akrun())
# user system elapsed
# 0.013 0.000 0.013
library(microbenchmark)
microbenchmark(greg(), akrun(), times=20L, unit='relative')
# Unit: relative
# expr min lq mean median uq max neval cld
# greg() 39.63122 38.53662 10.23198 18.81274 6.568741 4.642702 20 b
# akrun() 1.00000 1.00000 1.00000 1.00000 1.000000 1.000000 20 a
注意:基准测试因列数和行数而异。在这里,我使用的是 OP post.
中显示的列数数据
set.seed(24)
dat <- data.frame(a=runif(5), b=runif(5), c=runif(5), d=runif(5))
由于 model.matrix
抱怨只有一个级别的因素,您可能想要使用 stats::terms
labels(terms(~.^2, data = iris[, 1:3]))
# [1] "Sepal.Length" "Sepal.Width" "Petal.Length"
# [4] "Sepal.Length:Sepal.Width" "Sepal.Length:Petal.Length" "Sepal.Width:Petal.Length"