Проверка комбинаций регрессоров - PullRequest
Новая
       8

Проверка комбинаций регрессоров

0 голосов
/ 08 ноября 2018

У меня есть фрейм данных, который содержит 100 значений для 6 регрессоров, x1-x6 и 100 значений для независимой переменной, y. Моя цель - оценить множественную линейную регрессию y на 2 регрессорах и выбрать модель с наивысшим квадратом R. Мне нужно проверить все возможные комбинации х. Например, оцените модели y по x1, x2; у на х1, х3; у на х2, х3 и тд. Как я могу проверить все эти возможные комбинации, а затем запустить все регрессии? Возможно, мне нужно как-то использовать функцию combn (), но я не знаю, как это сделать вместе с оценкой регрессии

1 Ответ

0 голосов
/ 08 ноября 2018

Сначала я создам воспроизводимый набор данных, поскольку вы не опубликовали ни одного. 

set.seed(1)    # make the results reproducible
n <- 100
X <- matrix(rnorm(6*n), ncol = 6)
y <- X %*% sample(6) + rnorm(n)

X <- as.data.frame(X)
names(X) <- paste0("x", 1:6)

data <- cbind(X, y)
head(data)

Следующее решает проблему в последовательных циклах *apply. Это происходит в следующих шагах.

  1. Получите полные имена регрессоров, зная только, что они начинаются с "x".
  2. Создайте матрицу из всех возможных комбинаций два на два.
  3. Запустить все возможные регрессии в первом цикле lapply. Формула соединяется с paste и затем приводится к классу "formula".
  4. Получите summary, где значение R ^ 2 вычисляется с другим lapply.
  5. Вот и все. Теперь в этих двух списках есть вся необходимая информация, и стандартные операторы поднабора, такие как '[[', могут извлекать все, что нам нужно.
    1. R ^ 2 сначала извлекаются с sapply, поскольку они образуют вектор значений.
    2. Извлечение коэффициентов модели с максимальным значением R ^ 2 и резюме модели.

Итак, давайте сделаем это. 

regress <- grep("^x", names(data), value = TRUE)

regress_mat <- combn(regress, 2)

lm_list <- apply(regress_mat, 2, function(reg){
  fmla <- paste("y", paste(reg, collapse = "+"), sep = "~")
  fmla <- as.formula(fmla)
  lm(fmla, data)
})

smry_list <- lapply(lm_list, summary)
rsq <- sapply(smry_list, '[[', "r.squared")

coef(lm_list[[which.max(rsq)]])
#(Intercept)          x3          x4 
# -0.1400852   5.0556334   5.7352798

smry_list[[which.max(rsq)]]

#Call:
#  lm(formula = fmla, data = data)
#
#Residuals:
#     Min       1Q   Median       3Q      Max 
#-11.6529  -4.6230   0.1127   3.7821  11.9342 
#
#Coefficients:
#              Estimate Std. Error t value Pr(>|t|)    
#  (Intercept)  -0.1401     0.5763  -0.243    0.808    
#  x3            5.0556     0.5626   8.987 2.07e-14 ***
#  x4            5.7353     0.5867   9.775 4.11e-16 ***
#  ---
#  Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
#Residual standard error: 5.754 on 97 degrees of freedom
#Multiple R-squared:  0.6714,   Adjusted R-squared:  0.6646 
#F-statistic:  99.1 on 2 and 97 DF,  p-value: < 2.2e-16
...