如何在Julia中复现R中(y ~ (x1+x2+x3+x4+x5)^2)的模型矩阵？

R中的目标结果 #

以下是R中生成目标模型矩阵的代码及输出：

data = read.csv("https://gist.githubusercontent.com/TJhon/158daa0c2dd06010d01a72dae2af8314/raw/61df065c98ec90b9ea3b8598d1996fb5371a64aa/rnd.csv")

head(model.matrix(y ~ (x1 + x2 + x3 + x4 + x5)^2, data), 3)
#>   (Intercept)         x1         x2         x3        x4 x5       x1:x2
#> 1           1 -0.3007225 -1.3710894  0.3423409  1.322547  2  0.41231744
#> 2           1  0.4674170  0.8728939  0.9534157 -1.007083  1  0.40800548
#> 3           1  0.2085316 -0.3657995 -0.3043694 -1.036938  4 -0.07628076
#>         x1:x3      x1:x4      x1:x5      x2:x3      x2:x4      x2:x5      x3:x4
#> 1 -0.10294961 -0.3977198 -0.6014450 -0.4693799 -1.8133307 -2.7421787  0.4527620
#> 2  0.44564276 -0.4707279  0.4674170  0.8322308 -0.8790769  0.8728939 -0.9601690
#> 3 -0.06347064 -0.2162343  0.8341265  0.1113382  0.3793113 -1.4631979  0.3156121
#>        x3:x5     x4:x5
#> 1  0.6846817  2.645095
#> 2  0.9534157 -1.007083
#> 3 -1.2174775 -4.147751

尝试的Julia代码及问题 #

尝试了三种公式写法，但结果均不符合预期：

• 使用*运算符：生成包含平方项的额外列，矩阵维度为9×31，与R结果不符

  ModelMatrix(ModelFrame(@formula(y ~ (x1 + x2 + x3 + x4 + x5) * (x1 + x2 + x3 + x4 + x5)), data)).m
  

• 使用&运算符：仅生成交互项，缺少主效应列，矩阵维度为9×26，与R结果不符

  ModelMatrix(ModelFrame(@formula(y ~ (x1 + x2 + x3 + x4 + x5) & (x1 + x2 + x3 + x4 + x5)), data)).m
  

• 使用^2运算符：仅生成截距和线性预测值的平方，矩阵维度为9×2，完全不符合需求

  ModelMatrix(ModelFrame(@formula(y ~ (x1 + x2 + x3 + x4 + x5)^2), data)).m
  

需求是获取与R完全一致的数值矩阵和变量名，以便转换为DataFrame。

解决方案 #

核心差异是公式语法：R中(x1+x2+...)^2代表主效应 + 所有两两交互项（不含平方项），而Julia默认^2是对整个表达式平方，*会包含平方项。以下两种方法可实现等价效果：

方法1：显式定义主效应+交互项 #

using CSV, DataFrames, StatsModels

# 读取数据
data = CSV.read(download("https://gist.githubusercontent.com/TJhon/158daa0c2dd06010d01a72dae2af8314/raw/61df065c98ec90b9ea3b8598d1996fb5371a64aa/rnd.csv"), DataFrame)

# 构建等价公式：主效应 + 所有两两交互
f = @formula(y ~ x1 + x2 + x3 + x4 + x5 + (x1 + x2 + x3 + x4 + x5) & (x1 + x2 + x3 + x4 + x5))

# 生成模型框架和矩阵
mf = ModelFrame(f, data)
mm = ModelMatrix(mf).m

# 转换为带变量名的DataFrame
model_df = DataFrame(mm, coefnames(mf))

# 查看前3行
first(model_df, 3)

方法2：用*运算符后剔除平方项 #

# 另一种等价写法：先通过*生成所有项，再减去平方项
f = @formula(y ~ (x1 + x2 + x3 + x4 + x5) * (x1 + x2 + x3 + x4 + x5) - (x1^2 + x2^2 + x3^2 + x4^2 + x5^2))

结果验证 #

生成的DataFrame前3行数值和变量名与R输出完全一致，包含(Intercept)、所有主效应列，以及x1:x2、x1:x3等两两交互项，无多余平方项。

内容的提问来源于stack exchange，提问作者Jhon Kevin Flores Rojas