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

Question

Я потратил дни на поиск оптимальных моделей, которые соответствовали бы всем стандартным допущениям OLS (нормальное распределение, гомоскедастичность, отсутствие мультиколлинеарности) в R, но с 12 переменными невозможно найти оптимальную комбинацию переменных. Поэтому я пытался создать скрипт, который бы автоматизировал этот процесс.

Вот пример кода для расчетов:

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 <- as.data.frame(cbind(x1,x2,x3,x4,x5))

library(lmtest)
library(car)

model <- lm(x1~x2+x3+x4+x5, data = df)

# check for normal distribution (Shapiro-Wilk-Test)
rs_sd <- rstandard(model)
shapiro.test(rs_sd)

# check for heteroskedasticity (Breusch-Pagan-Test)
bptest(model)

# check for multicollinearity
vif(model)

#-------------------------------------------------------------------------------
# models without outliers
# identify outliers (calculating the Cooks distance, if x > 4/(n-k-1) --> outlier
cooks <- round(cooks.distance(model), digits = 4)
df_no_out <- cbind(df, cooks)
df_no_out <- subset(df_no_out, cooks < 4/(100-4-1))

model_no_out <- lm(x1~x2+x3+x4+x5, data = df_no_out)

# check for normal distribution
rs_sd_no_out<- rstandard(model_no_out)
shapiro.test(rs_sd_no_out)

# check for heteroskedasticity
bptest(model_no_out)

# check for multicollinearity
vif(model_no_out)

Я имею в виду, что нужно пройтись по всем комбинациям var и получить P-VALUES для shapiro.test() и bptest() или VIF-значения для всех созданных моделей, чтобы я мог сравнить значения значимости или мультиколлинеарность соотв. (в моем наборе данных мультиколлинеарность не должна быть проблемой, и поскольку для проверки мультиколлинеарности тест VIF выдает больше значений (для каждого фактора var 1xVIF), которые, вероятно, будет более сложным для реализации в коде), p-значения для shapiro.test + bptest() будет достаточно...).

Я пытался написать несколько сценариев, которые бы автоматизировали процесс, но безуспешно (к сожалению, я не программист). Я знаю, что уже есть некоторые темы, связанные с этой проблемой

Как запустить модели lm, используя все возможные комбинации нескольких переменных и фактора

Нахождение наилучшей комбинации переменных для высоких значений R-квадрата

но я не нашел сценарий, который бы также вычислял просто P-VALUES.

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

Буду очень признателен за любые предложения или помощь в этом.

Answer 1

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

Answer 2

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

Существует одна функция, которая подходит для всех возможных моделей. Затем серия звонков на *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, ]