如何在 Julia 中估计许多 GLM 模型？

Question

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

如何在不为 GLM 手动输入 4999 公式的情况下做到这一点？

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

提前致谢。

Answer 1

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

从虚拟数据开始

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)返回一个Strings的向量，而Term()构造函数需要Symbols。sterm()有一个String自动处理转换的方法。