R Automati c Регрессии - PullRequest
Новая
       26

R Automati c Регрессии

0 голосов
/ 27 апреля 2020 
data1=mtcars
CONTROL = c("mpg", "cyl")
X = c("hp","drat","wt","am")
Y = c("vs")

model1= glm(vs ~ hp + mpg + cyl, family = binomial, data = data1)
model2= glm(vs ~ drat + mpg + cyl, family = binomial, data = data1)
model3= glm(vs ~ wt + mpg + cyl, family = binomial, data = data1)
model4= glm(vs ~ am + mpg + cyl, family = binomial, data = data1)

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

i=1:4
model <- vector("list", length(i))
for(1:i){
  model[i]=glm(Y~X[i]+CONTROL,family = binomial, data = data1)
}

Мой желаемый вывод - это кадр данных, такой что: 

VARNAME COEF LL UL
hp
drat
wt
am

Ответы [ 2 ]

1 голос
/ 27 апреля 2020

Мы предполагаем, что вы хотите регрессировать против всех наборов из 3 переменных из 4 в X плюс управляющие переменные.

Кроме того, схема желаемого результата неясна, поэтому мы предположим, что вы хотите, чтобы коэффициенты каждой модели наряду с некоторыми другими неопределенными статистическими данными были использованы, и мы используем логарифмическую вероятность в качестве примера. Каждая выходная строка представляет результаты одного прогона модели с NA в каждой строке, показывающие, какой коэффициент не использовался.

Пакеты не используются. 

# test data is the builtin mtcars as well as CONTROL, X and Y
CONTROL <- c("mpg", "cyl")
X <- c("hp","drat","wt","am")
Y <- "vs"

stats <- function(nm) {
  fo <- reformulate(c(setdiff(X, nm), CONTROL), Y)
  fm <- glm(fo, mtcars, family = binomial)
  coefs <- c(coef(fm), setNames(NA, nm))[c("(Intercept)", X)]
  c(coefs, logLik = logLik(fm))  # add other statistics to this line
}

do.call("rbind", lapply(X, stats))

, давая эту матрицу: 

     (Intercept)          hp        drat        wt         am        logLik
[1,]   150.84593          NA   -4.381111  2.242363  -62.74444 -2.338630e+00
[2,]   137.24734 -0.01079679          NA  1.209804  -60.16064 -2.826391e+00
[3,]  3285.91281 -6.54621462 -344.866970        NA -311.77291 -1.315947e-08
[4,]    84.76241 -0.49550208  -25.801081 24.231253         NA -3.863991e+00
0 голосов
/ 27 апреля 2020

Это даст вам большую часть пути: 

library(broom)
library(dplyr)

char_formula = sprintf("%s ~ %s + %s", Y, paste(CONTROL, collapse = "+"), X)

mods = list()
for(i in seq_along(char_formula)) {
  mods[[i]] = glm(as.formula(char_formula[i]), family = binomial, data = data1)
}

names(mods) = X
lapply(mods, tidy) %>%
  bind_rows(.id = "VARNAME") %>%
  filter(term %in% X)
# # A tibble: 4 x 6
#   VARNAME term  estimate  std.error statistic p.value
#   <chr>   <chr>    <dbl>      <dbl>     <dbl>   <dbl>
# 1 hp      hp     -0.0461     0.0362  -1.27      0.203
# 2 drat    drat   -2.53       1.83    -1.38      0.168
# 3 wt      wt      2.26       1.73     1.30      0.192
# 4 am      am    -61.0    18184.      -0.00335   0.997
...