автоматизация lm-тестов со всеми возможными комбинациями переменных и получением значений для: shapiro.test(), bptest(),vif() в R

вы царапаете поверхность того, что сейчас называется статистическим обучением. вступительный текст - «Статистическое обучение с приложениями в R», а текст уровня выпускника - «Элементы статистического обучения». чтобы сделать то, что вам нужно, вы используете функцию regsubsets() из пакета "Leaps". Однако, если вы прочтете хотя бы главу 6 из вводной книги, вы узнаете о перекрестной проверке и начальной загрузке, которые являются современным способом выбора модели.

Следующее автоматизирует подгонку моделей и последующие тесты.

Есть одна функция, которая подходит для всех возможных моделей. Затем серия вызовов функций *apply позволит получить нужные вам значения.

library(lmtest)
library(car)


fitAllModels <- function(data, resp, regr){
  f <- function(M){
    apply(M, 2, function(x){
      fmla <- paste(resp, paste(x, collapse = "+"), sep = "~")
      fmla <- as.formula(fmla)
      lm(fmla, data = data)
    })
  }
  regr <- names(data)[names(data) %in% regr]
  regr_list <- lapply(seq_along(regr), function(n) combn(regr, n))
  models_list <- lapply(regr_list, f)
  unlist(models_list, recursive = FALSE)
}

Теперь данные.

# Make up a data.frame to test the function above.
# Don't forget to set the RNG seed to make the
# results reproducible
set.seed(7646)
x1 <- runif(100, 0, 10)
x2 <- runif(100, 0, 10)
x3 <- runif(100, 0, 10)
x4 <- runif(100, 0, 10)
x5 <- runif(100, 0, 10)

df <- data.frame(x1, x2, x3, x4, x5)

Сначала подберите все модели с "x1" в качестве ответа и другими переменными в качестве возможных регрессоров. Функция может быть вызвана с одним ответом и любым количеством возможных регрессоров, которые вы хотите.

fit_list <- fitAllModels(df, "x1", names(df)[-1])

А теперь последовательность испытаний.

# Normality test, standardized residuals
rs_sd_list <- lapply(fit_list, rstandard)
sw_list <- lapply(rs_sd_list, shapiro.test)
sw_pvalues <- sapply(sw_list, '[[', 'p.value')

# check for heteroskedasticity (Breusch-Pagan-Test)
bp_list <- lapply(fit_list, bptest)
bp_pvalues <- sapply(bp_list, '[[', 'p.value')

# check for multicollinearity, 
# only models with 2 or more regressors
vif_values <- lapply(fit_list, function(fit){
  regr <- attr(terms(fit), "term.labels")
  if(length(regr) < 2) NA else vif(fit)
})

Заметка о расстоянии Кука. В вашем коде вы подмножаете исходный data.frame, создавая новый без выбросов. Это приведет к дублированию данных. Я выбрал только список индексов строк df. Если вы предпочитаете дублированные data.frames, раскомментируйте строку в анонимной функции ниже и закомментируйте последнюю.

# models without outliers
# identify outliers (calculating the 
# Cooks distance, if x > 4/(n - k - 1) --> outlier

df_no_out_list <- lapply(fit_list, function(fit){
  cooks <- cooks.distance(fit)
  regr <- attr(terms(fit), "term.labels")
  k <- length(regr)
  inx <- cooks < 4/(nrow(df) - k - 1)
  #df[inx, ]
  which(inx)
})

# This tells how many rows have the df's without outliers
sapply(df_no_out_list, NROW)

# A data.frame without outliers. This one is the one 
# for model number 8. 
# The two code lines could become a one-liner.
i <- df_no_out_list[[8]]
df[i, ]