Логистическая регрессия на группу

Question

У меня есть такой же фрейм данных, что и у набора данных iris:

library(datasets)
df <- iris
head(df)

Я пытаюсь увидеть, какой столбец / переменная лучше всего предсказывает каждый вид.

Я пробовал логистическую регрессию:

model <- glm(Species ~
                  Petal.Width+
                  Petal.Length+
                  Sepal.Width+
                  Sepal.Length,
                data=df,
                family = binomial(logit))

model 
summary(model)
anova(model)

Это показывает мне, какая переменная является наиболее значимой.Тем не менее, я хотел бы наблюдать это для всего набора данных и также для каждого вида (то есть, какая переменная лучше всего предсказывает каждый вид).

Любая помощь будет очень признательна.

Answer 1

Вы можете group_by(Species) перед запуском регрессий.Рассмотрим следующий пример:

library(datasets)
df <- iris
### models for each species
library(dplyr)
models <-  df %>% 
        group_by(Species) %>% 
        do(model = glm(Species ~
                               Petal.Width+
                               Petal.Length+
                               Sepal.Width+
                               Sepal.Length,
                       data = .,
                       family = binomial(logit)))
models$model # list of the three models with their estimates
# output
[[1]]

Call:  glm(formula = Species ~ Petal.Width + Petal.Length + Sepal.Width + 
    Sepal.Length, family = binomial(logit), data = .)

Coefficients:
 (Intercept)   Petal.Width  Petal.Length   Sepal.Width  Sepal.Length  
  -2.657e+01    -1.391e-14    -2.777e-15     3.349e-15     1.244e-15  

Degrees of Freedom: 49 Total (i.e. Null);  45 Residual
Null Deviance:      0 
Residual Deviance: 2.901e-10    AIC: 10

[[2]]

Call:  glm(formula = Species ~ Petal.Width + Petal.Length + Sepal.Width + 
    Sepal.Length, family = binomial(logit), data = .)

Coefficients:
 (Intercept)   Petal.Width  Petal.Length   Sepal.Width  Sepal.Length  
  -2.657e+01    -6.410e-14    -1.630e-14     3.634e-14     4.943e-14  

Degrees of Freedom: 49 Total (i.e. Null);  45 Residual
Null Deviance:      0 
Residual Deviance: 2.901e-10    AIC: 10

[[3]]

Call:  glm(formula = Species ~ Petal.Width + Petal.Length + Sepal.Width + 
    Sepal.Length, family = binomial(logit), data = .)

Coefficients:
 (Intercept)   Petal.Width  Petal.Length   Sepal.Width  Sepal.Length  
  -2.657e+01     1.135e-15     1.749e-14    -6.585e-15    -9.822e-15  

Degrees of Freedom: 49 Total (i.e. Null);  45 Residual
Null Deviance:      0 
Residual Deviance: 2.901e-10    AIC: 10

### to get std.error, p.value, etc.
library(broom)
models %>% tidy(model)
# output
# A tibble: 15 x 6
# Groups:   Species [3]
      Species         term      estimate std.error     statistic   p.value
       <fctr>        <chr>         <dbl>     <dbl>         <dbl>     <dbl>
 1     setosa  (Intercept) -2.656607e+01  786965.7 -3.375760e-05 0.9999731
 2     setosa  Petal.Width -1.390512e-14  523902.1 -2.654145e-20 1.0000000
 3     setosa Petal.Length -2.776770e-15  316800.5 -8.765044e-21 1.0000000
 4     setosa  Sepal.Width  3.349353e-15  200749.8  1.668422e-20 1.0000000
 5     setosa Sepal.Length  1.243936e-15  221428.0  5.617790e-21 1.0000000
 6 versicolor  (Intercept) -2.656607e+01  615841.5 -4.313783e-05 0.9999656
 7 versicolor  Petal.Width -6.410448e-14  475196.2 -1.349011e-19 1.0000000
 8 versicolor Petal.Length -1.629537e-14  226400.3 -7.197591e-20 1.0000000
 9 versicolor  Sepal.Width  3.634246e-14  225869.3  1.609004e-19 1.0000000
10 versicolor Sepal.Length  4.943463e-14  156829.6  3.152123e-19 1.0000000
11  virginica  (Intercept) -2.656607e+01  608647.0 -4.364774e-05 0.9999652
12  virginica  Petal.Width  1.135375e-15  223583.6  5.078078e-21 1.0000000
13  virginica Petal.Length  1.748522e-14  186273.0  9.386880e-20 1.0000000
14  virginica  Sepal.Width -6.584951e-15  202705.8 -3.248527e-20 1.0000000
15  virginica Sepal.Length -9.821934e-15  165119.0 -5.948396e-20 1.0000000