Извлечение опорного уровня от коэффициентов GLM - PullRequest
Новая
       28

Извлечение опорного уровня от коэффициентов GLM

0 голосов
/ 10 июня 2018

Я понимаю, что опорный уровень не включен, но я хотел бы получить возможность взять подобранный объект glm и выяснить, какими были опорные уровни (то есть, не используя знания исходного набора данных).Хранится ли это где-нибудь в glm приспособленном объекте?

Пример данных ниже: 

> btest <- data.frame(var1 = sample(c(1,2,3), 100, replace = T),
+                     var2 = sample(c('a','b','c'), 100, replace = T),
+                     var3 = sample(c('e','f','g'), 100, replace = T),
+                     var4 = rnorm(100, mean = 3, 2),
+                     var5 = sample(c('yes','no'), 100, replace = T))
> summary(glm(var5 ~ var1 + var2 + var3 + var4, data = btest, family = 'binomial'))

Call:
glm(formula = var5 ~ var1 + var2 + var3 + var4, family = "binomial", 
    data = btest)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.6988  -1.0457  -0.6213   1.1224   1.8904  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)  
(Intercept) -0.81827    0.73173  -1.118   0.2635  
var1         0.55923    0.27279   2.050   0.0404 *
var2b       -0.60998    0.53435  -1.142   0.2536  
var2c       -0.60250    0.51706  -1.165   0.2439  
var3f       -0.81899    0.53345  -1.535   0.1247  
var3g        0.21215    0.51907   0.409   0.6828  
var4         0.04429    0.12650   0.350   0.7263  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 137.99  on 99  degrees of freedom
Residual deviance: 128.35  on 93  degrees of freedom
AIC: 142.35

Number of Fisher Scoring iterations: 4

Здесь я хотел бы знать, что var1 и var4 не имеют ссылки, но что контрольный уровень для var2 и var3 являются 'a' и 'e' соответственно.Потому что в конечном итоге я получу таблицу с NA для Estimate для этих переменных на этих эталонных уровнях.

Редактировать: Для людей, которые приходят позже, мне также интересно, в какой степени может помочь подключение к элементу terms приспособленного glm объекта в сочетании с ответом ниже ... 

> btest2 <- glm(var5 ~ var1 + var3 + var2 + var4, data = btest, family = 'binomial')
> btest2$terms
var5 ~ var1 + var3 + var2 + var4
attr(,"variables")
list(var5, var1, var3, var2, var4)
attr(,"factors")
     var1 var3 var2 var4
var5    0    0    0    0
var1    1    0    0    0
var3    0    1    0    0
var2    0    0    1    0
var4    0    0    0    1
attr(,"term.labels")
[1] "var1" "var3" "var2" "var4"
attr(,"order")
[1] 1 1 1 1
attr(,"intercept")
[1] 1
attr(,"response")
[1] 1
attr(,".Environment")
<environment: R_GlobalEnv>
attr(,"predvars")
list(var5, var1, var3, var2, var4)
attr(,"dataClasses")
     var5      var1      var3      var2      var4 
 "factor" "numeric"  "factor"  "factor" "numeric" 
> attr(btest2$terms, 'dataClasses')
     var5      var1      var3      var2      var4 
 "factor" "numeric"  "factor"  "factor" "numeric"

Ответы [ 2 ]

0 голосов
/ 28 февраля 2019

Вот функция, которая извлекает xlevels и использует broom::tidy (с некоторыми дополнительными манипуляциями), так что опорный уровень находится в информационном кадре со всеми другими вашими терминами: 

library(tidyverse)
library(broom)

tidy_coefs_with_ref <- function(mod_obj, sep = "_"){

  tidy_coefs <- tidy(mod_obj) %>% 
    separate(term, c("variable", "level"), sep, remove = FALSE) %>% 
    mutate(variable = paste0(variable, sep))

  xlevels <- mod_obj$xlevels  

  missing_levels <- xlevels %>% 
    enframe() %>% 
    unnest() %>% 
    set_names(c("variable", "level"))

  missing_levels %>% 
    anti_join(tidy_coefs) %>% 
    bind_rows(tidy_coefs) %>% 
    arrange(variable, level)

}

btest <- tibble(var1 = sample(c(1,2,3), 100, replace = T),
                var2 = sample(c('a','b','c'), 100, replace = T),
                var3 = sample(c('e','f','g'), 100, replace = T),
                var4 = rnorm(100, mean = 3, 2),
                var5 = sample(c(TRUE, FALSE), 100, replace = T)) %>% 
  rename_if(is.character, funs(paste0(., "_")))

btest2 <- glm(var5 ~ ., data = btest, family = 'binomial')

tidy_coefs_with_ref(btest2)
#> Warning: Expected 2 pieces. Missing pieces filled with `NA` in 3 rows [1,
#> 2, 7].
#> Joining, by = c("variable", "level")
#> # A tibble: 9 x 7
#>   variable     level term         estimate std.error statistic p.value
#>   <chr>        <chr> <chr>           <dbl>     <dbl>     <dbl>   <dbl>
#> 1 (Intercept)_ <NA>  (Intercept)   0.904       0.835    1.08     0.279
#> 2 var1_        <NA>  var1         -0.126       0.269   -0.468    0.640
#> 3 var2_        a     <NA>         NA          NA       NA       NA    
#> 4 var2_        b     var2_b       -0.719       0.515   -1.40     0.162
#> 5 var2_        c     var2_c       -0.632       0.525   -1.21     0.228
#> 6 var3_        e     <NA>         NA          NA       NA       NA    
#> 7 var3_        f     var3_f       -0.379       0.496   -0.764    0.445
#> 8 var3_        g     var3_g        0.429       0.517    0.829    0.407
#> 9 var4_        <NA>  var4         -0.00833     0.111   -0.0749   0.940

Создано в 2019-02-28 пакетом представлений (v0.2.1)

(Шаг с seperate может быть очищен.)

Также потенциально уместно, вот ссылка на суть, где я использую (расширение) вышеупомянутую функцию с кодированием эффекта , чтобы также извлечь величину эффекта отброшенного уровня: https://gist.github.com/brshallo/f923b9b5c6360ce09beda35c3d1d55e9

0 голосов
/ 10 июня 2018

Если вы сохраняете подгонку, как показано ниже, в переменную my_fit, вы можете сделать my_fit$xlevels.Для всех категориальных переменных вы увидите все их уровни.

Затем вы можете соотнести это с моделью.Например, var1 не в xlevels, поэтому он непрерывный.Var2 имеет 3 уровня (a, b, c,), и у вас есть оценки для b и c.Это означает, что это ссылка.В Var3 есть категории e, f, g, и у вас есть оценки для f и g, поэтому e должно быть ссылкой. 

my_fit <- glm(var5 ~ var1 + var2 + var3 + var4, data = btest, family = 'binomial')


my_fit$xlevels
$var2
[1] "a" "b" "c"

$var3
[1] "e" "f" "g"


> summary(my_fit)

Call:
glm(formula = var5 ~ var1 + var2 + var3 + var4, family = "binomial", 
    data = btest)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.0344  -1.1100   0.5975   0.9605   1.9985  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)   
(Intercept)   0.7321     0.8111   0.903  0.36673   
var1         -0.5061     0.2846  -1.778  0.07533 . 
var2b        -0.1929     0.5904  -0.327  0.74385   
var2c        -0.1968     0.5442  -0.362  0.71764   
var3f        -1.1015     0.5816  -1.894  0.05824 . 
var3g        -0.2004     0.5629  -0.356  0.72187   
var4          0.3945     0.1290   3.058  0.00223 **
---
...