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

如何解决如何在Julia中估算许多GLM模型？

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

如果不为GLM手动键入4999公式怎么办？

在RI中，只需定义4999个字符串的列表（“ target〜x1”），将每个字符串转换为公式，然后使用map来估算多个glm。在Julia中可以做类似的事情吗？替代方案？

谢谢。

解决方法

您可以通过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}}}},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,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),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)返回了{{ 1}}，而String构造函数期望Term()。 Symbol具有term()的方法，可以自动处理转换。

,

您还可以使用低级API并直接将自变量作为向量传递，将自变量作为矩阵直接传递，甚至无需建立公式。您将丢失系数名称，但是由于每个模型中只有一个自变量，因此可能没问题。

此内容记录在?fit中。每个模型的调用看起来像glm([ones(length(x1)) x1],target,dist)。满一的列用于拦截。