Сохранение имен в виде списка из функции, выполняющей многочисленные одномерные регрессионные модели - PullRequest
Новая
       69

Сохранение имен в виде списка из функции, выполняющей многочисленные одномерные регрессионные модели

0 голосов
/ 24 мая 2018

Я строю модель логистической регрессии с набором данных, содержащим около 40 переменных.Первый шаг, который я использую при построении этих типов моделей, - это то, что я запускаю каждую переменную с помощью DV (Hosmer, Lemeshow, & Sturdivant, 2013).Я построил функцию, которая делает это для меня и возвращает значение p для каждого. 

Fit Univariate logistic regression model for each covariate
uni.log2 <- function(x) {
  log.mod2 <- glm(Renewf ~ x, data = dt.train2, family = binomial())
  return(coef(summary(log.mod2))[,4]) #get p-values only
}

Затем я применяю эту функцию к каждому из выбранных столбцов в моем dt 

#apply function to selected IV's 
apply(X = dt.train2[c(3:16)], MARGIN = 2, FUN = uni.log2)
* 1006.* Следующим шагом, который я хотел бы сделать, является проверка этих переменных для значений p с порогом p <0,25 и возвращение списка имен переменных, которые были однозначно значимы при p <0,25. </p>

Кто-нибудь знает, как это можно сделать?

Я могу установить порог и скопировать список имен из многомерной модели, используя этот код: 

threshold <- 0.001
signif_form <- as.formula(paste("Renewf ~ 
",paste(names(which((summary(log.mod2)$coefficients[2: 
(nrow(summary(log.mod2)$coefficients)), 4] < threshold) == TRUE)), collapse 
= "+")))

Но, опять жеЯ не знаю, как вставить имена из серии одномерных регрессионных моделей.Если кто-то знает, как это сделать, я был бы очень признателен за помощь.

Заранее спасибо!

1 Ответ

0 голосов
/ 24 мая 2018

Если вы все еще хотите использовать этот подход после просмотра ссылки, предоставленной @BenBolker (и, возможно, других ресурсов, представляющих опасность ступенчатая регрессия и статистическая значимость ) ...

Следующий код вернет вектор значений p для независимой переменной в каждой регрессии.Я использовал для иллюстрации встроенный фрейм данных mtcars. 

library(tidyverse)
library(broom)

pvals = sapply(names(mtcars)[names(mtcars) != "vs"], function(x) {
  glm(paste("vs ~ ", x), data=mtcars, family=binomial) %>% 
    tidy %>% 
    filter(term==x) %>% pull(p.value)
})

pvals

        mpg         cyl        disp          hp        drat          wt        qsec          am 
0.006590045 0.001917098 0.002453817 0.012340143 0.021777872 0.008672977 0.008813419 0.343628917 
       gear        carb 
0.250981095 0.004157666

Приведенный выше код использует оператор канала (%>%) для объединения функций,После создания модели с glm, tidy возвращает коэффициенты и значения p в виде фрейма данных: 

glm(vs ~ mpg, data=mtcars, family=binomial) %>% 
    tidy

         term   estimate std.error statistic     p.value
1 (Intercept) -8.8330726  3.162274 -2.793266 0.005217877
2         mpg  0.4304135  0.158422  2.716880 0.006590045

Затем filter и pull функции выбирают p-значение для конкретной рассматриваемой переменной: 

glm(vs ~ mpg, data=mtcars, family=binomial) %>% 
  tidy %>% filter(term=="mpg") %>% pull(p.value)

[1] 0.006590045

Оборачивание всего этого в sapply возвращает именованный вектор p-значений, гдеимена являются независимыми переменными в каждой одномерной регрессии.

Чтобы вернуть только элементы ниже порогового значения p: 

pvals[pvals < 0.25]

        mpg         cyl        disp          hp        drat          wt        qsec        carb 
0.006590045 0.001917098 0.002453817 0.012340143 0.021777872 0.008672977 0.008813419 0.004157666

Если вы просто хотитеимена переменных, которые соответствуют критерию порога: 

names(pvals[pvals < 0.25])

Чтобы непосредственно вернуть только элементы ниже порога p-значения: 

pvals = sapply(names(mtcars)[names(mtcars) != "vs"], function(x) {
  glm(paste("vs ~ ", x), data=mtcars, family=binomial) %>% 
    tidy %>% 
    filter(term==x) %>% pull(p.value)
}) %>% .[. < 0.25]

Наконец, упаковав это как функцию длявернуть имена нужных переменных: 

select_vars = function(DV, data, threshold) {
  sapply(names(data)[names(data) != DV], function(x) {
    glm(paste(DV, " ~ ", x), data=data, family=binomial) %>% 
      tidy %>% 
      filter(term==x) %>% pull(p.value)
  }) %>% .[. < threshold] %>% names
}

select_vars("vs", mtcars, 0.25)

[1] "mpg"  "cyl"  "disp" "hp"   "drat" "wt"   "qsec" "carb"

select_vars("Species", iris %>% filter(Species %in% c("versicolor","virginica")), 0.001)

[1] "Sepal.Length" "Petal.Length" "Petal.Width"
...