Comment estimer plusieurs modèles GLM dans Julia?

Question

J'ai un ensemble de données de 5000 variables. Une cible et 4999 covariables. Je veux estimer un glm pour chaque combinaison cible-variable (4999 modèles).

Comment puis-je faire cela sans taper manuellement 4999 formules pour GLM?

Dans RI, il suffit de définir une liste de 4999 chaînes ("cible ~ x1), de convertir chaque chaîne en une formule et d'utiliser map pour estimer plusieurs glm. Y a-t-il quelque chose de similaire qui peut être fait dans Julia? Ou existe-t-il une alternative élégante?

Merci d'avance.

Answer 1

Vous pouvez créer la formule par programmation via des Termobjets. La documentation pour cela peut être trouvée ici , mais considérez l'exemple simple suivant qui devrait répondre à vos besoins:

Commencez avec des données factices

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 │

Maintenant, lorsque vous exécutez un modèle linéaire avec GLM, vous feriez quelque chose comme lm(@formula(y ~ x1), df), qui ne peut en effet pas être facilement utilisé dans une boucle pour construire différentes formules. Nous allons donc suivre la documentation et créer @formuladirectement la sortie de la macro - rappelez-vous que les macros dans Julia transforment simplement la syntaxe en une autre syntaxe, donc elles ne font rien que nous ne pouvons pas écrire nous-mêmes!

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
──────────────────────────────────────────────────────────────────────────

Vous pouvez vérifier par vous-même que ce qui précède est équivalent à lm(@formula(y ~ x1), df).

Nous espérons maintenant que c'est une étape facile pour créer la boucle que vous recherchez (limitée à deux covariables ci-dessous pour limiter la sortie):

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
─────────────────────────────────────────────────────────────────────────

Comme Dave le souligne ci-dessous, il est utile d'utiliser la term()fonction ici pour créer nos termes plutôt que le Term()constructeur directement - c'est parce que names(df)renvoie un vecteur de Strings, alors que le Term()constructeur attend Symbols. term()a une méthode pour Strings qui gère la conversion automatiquement.

Answer 2

Vous pouvez également utiliser l'API de bas niveau et transmettre la variable dépendante en tant que vecteur et la variable indépendante en tant que matrice directement sans même créer de formules. Vous perdrez les noms de coefficient, mais comme vous n'avez qu'une seule variable indépendante dans chaque modèle, c'est probablement OK.

Ceci est documenté dans ?fit. L'appel pour chaque modèle ressemblera à glm([ones(length(x1)) x1], target, dist). La colonne pleine de uns est pour l'interception.