julia - 如何在 Julia 中估计多个 GLM 模型？

Question

我有一个包含 5000 个变量的数据集。一个目标和 4999 个协变量。我想为每个目标变量组合(4999 个模型)估计一个 glm。

如果不为 GLM 手动输入 4999 公式，我该如何做到这一点？

在 R 中，我会简单地定义一个包含 4999 个字符串的列表 ("target ~ x1)，将每个字符串转换为一个公式并使用 map 来估计多个 glm。在 Julia 中是否可以做类似的事情？或者是否有一个优雅的替代方案？

提前致谢。

Answer 1

您可以通过 Term 对象以编程方式创建公式。可以找到相关文档 here ，但请考虑以下应该满足您需求的简单示例:

从虚拟数据开始

julia> using DataFrames, GLM

julia> df = hcat(DataFrame(y = rand(10)), DataFrame(rand(10, 5)))
10×6 DataFrame
│ Row │ y         │ x1        │ x2       │ x3        │ x4         │ x5       │
│     │ Float64   │ Float64   │ Float64  │ Float64   │ Float64    │ Float64  │
├─────┼───────────┼───────────┼──────────┼───────────┼────────────┼──────────┤
│ 1   │ 0.0200963 │ 0.924856  │ 0.947904 │ 0.429068  │ 0.00833488 │ 0.547378 │
│ 2   │ 0.169498  │ 0.0915296 │ 0.375369 │ 0.0341015 │ 0.390461   │ 0.835634 │
│ 3   │ 0.900145  │ 0.502495  │ 0.38106  │ 0.47253   │ 0.637731   │ 0.814095 │
│ 4   │ 0.255163  │ 0.865253  │ 0.791909 │ 0.0833828 │ 0.741899   │ 0.961041 │
│ 5   │ 0.651996  │ 0.29538   │ 0.161443 │ 0.23427   │ 0.23132    │ 0.947486 │
│ 6   │ 0.305908  │ 0.170662  │ 0.569827 │ 0.178898  │ 0.314841   │ 0.237354 │
│ 7   │ 0.308431  │ 0.835606  │ 0.114943 │ 0.19743   │ 0.344216   │ 0.97108  │
│ 8   │ 0.344968  │ 0.452961  │ 0.595219 │ 0.313425  │ 0.102282   │ 0.456764 │
│ 9   │ 0.126244  │ 0.593456  │ 0.818383 │ 0.485622  │ 0.151394   │ 0.043125 │
│ 10  │ 0.60174   │ 0.8977    │ 0.643095 │ 0.0865611 │ 0.482014   │ 0.858999 │

现在，当您使用 GLM 运行线性模型时，您会执行类似 lm(@formula(y ~ x1), df) 的操作，这确实不容易在循环中使用构造不同的公式。因此，我们将按照文档直接创建 @formula 宏的输出——记住 Julia 中的宏只是将语法转换为其他语法，所以它们不会做我们自己无法编写的任何事情!

julia> lm(Term(:y) ~ Term(:x1), df)
StatsModels.TableRegressionModel{LinearModel{GLM.LmResp{Array{Float64,1}},GLM.DensePredChol{Float64,LinearAlgebra.Cholesky{Float64,Array{Float64,2}}}},Array{Float64,2}}

y ~ 1 + x1

Coefficients:
──────────────────────────────────────────────────────────────────────────
                 Coef.  Std. Error      t  Pr(>|t|)   Lower 95%  Upper 95%
──────────────────────────────────────────────────────────────────────────
(Intercept)   0.428436    0.193671   2.21    0.0579  -0.0181696   0.875041
x1           -0.106603    0.304597  -0.35    0.7354  -0.809005    0.595799
──────────────────────────────────────────────────────────────────────────

你可以自己验证上面的等价于lm(@formula(y ~ x1), df)。

现在，希望它是构建您正在寻找的循环的简单步骤(限于以下两个协变量以限制输出):

julia> for x ∈ names(df[:, Not(:y)])[1:2]
           @show lm(term(:y) ~ term(x), df)
       end
lm(term(:y) ~ term(x), df) = StatsModels.TableRegressionModel{LinearModel{GLM.LmResp{Array{Float64,1}},GLM.DensePredChol{Float64,LinearAlgebra.Cholesky{Float64,Array{Float64,2}}}},Array{Float64,2}}

y ~ 1 + x1

Coefficients:
──────────────────────────────────────────────────────────────────────────
                 Coef.  Std. Error      t  Pr(>|t|)   Lower 95%  Upper 95%
──────────────────────────────────────────────────────────────────────────
(Intercept)   0.428436    0.193671   2.21    0.0579  -0.0181696   0.875041
x1           -0.106603    0.304597  -0.35    0.7354  -0.809005    0.595799
──────────────────────────────────────────────────────────────────────────
lm(Term(:y) ~ Term(x), df) = StatsModels.TableRegressionModel{LinearModel{GLM.LmResp{Array{Float64,1}},GLM.DensePredChol{Float64,LinearAlgebra.Cholesky{Float64,Array{Float64,2}}}},Array{Float64,2}}

y ~ 1 + x2

Coefficients:
─────────────────────────────────────────────────────────────────────────
                 Coef.  Std. Error      t  Pr(>|t|)  Lower 95%  Upper 95%
─────────────────────────────────────────────────────────────────────────
(Intercept)   0.639633    0.176542   3.62    0.0068   0.232527    1.04674
x2           -0.502327    0.293693  -1.71    0.1256  -1.17958     0.17493
─────────────────────────────────────────────────────────────────────────

正如 Dave 在下面指出的，在这里使用 term() 函数来创建我们的术语比直接使用 Term() 构造函数更有帮助 - 这是因为 names(df) 返回 String 的向量，而 Term() 构造函数需要 Symbol。 term() 有一个用于 String 的方法，可以自动处理转换。