Эта статья является (надеюсь) последней в серии статей, которые я опубликовал здесь, на Medium, в которых подробно объясняется, как воспроизвести мои результаты для тестового раздела предстоящей статьи, посвященной оценочной исчерпывающей регрессии, алгоритму выбора оптимальной переменной, который является новым. вычислительно осуществимый побочный продукт исчерпывающей регрессии, также известной как выбор лучшего подмножества или регрессия всех подмножеств. Предыдущие статьи можно найти здесь и здесь соответственно, а здесь — ссылка на репозиторий GitHub.
Контрольный показатель № 2: обратная пошаговая регрессия исключения
Я пробежался назад по 260 000 синтетических выборочных наборов данных, используя функцию step() из пакета stats в R. Прежде всего, просто как быстрый и простой способ проверить, правильно ли я пометил свои рабочие пространства. по мере их загрузки я включаю следующие три базовые команды R ненадолго после загрузки WorkSpace в верхнюю часть скрипта, который запускает наш алгоритм выбора оптимальной переменной 2nd Benchmark, и снова в самом конце:
length(datasets) head(DS_names_list) tail(DS_names_list)
Таким образом, если бы я назвал WorkSpace от «0,75–4–1–1 до 0,75–4–10–500», они должны вернуть: 5000, «0,75–4–1–1», «0,75–4–1–2». «0,75–4–1–3» «0,75–4–1–4» «0,75–4–1–5» «0,75–4–1–6» и «0,75–4–1–495» «0,75 –4–1–496» «0,75–4–1–497» «0,75–4–1–498» «0,75–4–1–499» «0,75–4–1–500» соответственно. Итак, если вместо этого он возвращает:
> length(datasets) [1] 10000 > head(DS_names_list) [1] "0.75-3-1-1" "0.75-3-1-2" "0.75-3-1-3" "0.75-3-1-4" "0.75-3-1-5" "0.75-3-1-6" > tail(DS_names_list) [1] "0.75-4-10-495" "0.75-4-10-496" "0.75-4-10-497" "0.75-4-10-498" "0.75-4-10-499" "0.75-4-10-500"
Затем я знаю, что допустил простую ошибку при присвоении имени этому WorkSpace ранее, и я могу исправить ее, когда я сохраняю показатели производительности для ступенчатой регрессии обратного исключения, выполняемой для этих наборов данных, в файл csv.
Код R для запуска второго теста
Вот как получился мой окончательный код после улучшения моего кода после первоначального использования неэффективного цикла for (большое нет-нет в R, если вы заботитесь о времени выполнения):
# load all necessary packages
library(plyr)
library(dplyr)
library(stringi)
library(purrr)
library(stats)
library(leaps)
library(lars)
library(data.table)
library(parallel)
### Step 2: Run a Backward Elimination Stepwise Regression
### function on each of the 260,000 datasets.
### Assign the full models to their corresponding datasets and
### store these in the object "full_models"
CL <- makeCluster(detectCores() - 3L)
system.time( clusterExport(CL, c('datasets')) )
set.seed(11) # for reproducibility
time_taken_to_fit_BE <- system.time( BE.fits <- parLapply(cl = CL, X = datasets, \(ds_i) {
full_models <- lm(ds_i$Y ~ ., data = ds_i)
back <- stats::step(object = full_models, scope = formula(full_models),
direction = 'back', trace = FALSE) }) )
time_taken_to_fit_BE
# Extract the elapsed time
BE_elapsed_time <- time_taken_to_fit_BE["elapsed"]
stopCluster(CL)
rm(CL)
Здесь много всего происходит, поэтому я постараюсь разбить это по пунктам шаг за шагом. Набор предварительно загруженных библиотек очень похож на случаи, когда я использовал функции enet() и lars() для запуска, а затем реплицировал 1-й тест, только с Библиотека purr добавлена к библиотеке эластичных сетей и без нее по очевидным причинам.
Функция makeCluster() и функция detectCores() взяты из параллельной библиотеки. Функция detectCores() просто расшифровывает, сколько ядер вам нужно выбрать, у моего ноутбука 8, поэтому 8 - 3 означает, что будет выделено 5 для любых задач, которые запускаются/выполняются с использованием CL после запуска этого строка кода.
Что касается лассо ранее, я установил случайное начальное значение 11 для всех регрессий BE, которые я запускаю, но вместо того, чтобы снова использовать lapply(), как я делал на каждом другом этапе пути до сих пор , я использую его версию из параллельной библиотеки, функцию parLapply(). Единственное отличие синтаксиса функции parLapply() от стандартной lapply() заключается в том, что вы должны явно указать имя списка установленных вами ядер. что для меня просто CL, сокращение от основного списка.
Во всех важных разделах этого конкретного скрипта, где оценивается/соответствует алгоритм эталонного оптимального выбора переменных для каждого из 260 тыс. наборов данных, я использовал parLapply(datasets), \(i)\n (где \n здесь просто псевдокод означает «остальная часть продолжается на следующей строке») синтаксис внутри итерационной части, как и в предыдущих сценариях, но на этот раз я использовал ds_i, сокращение от dataset_i, и, что более важно, Я также явно включил волнистые скобки {} функции R, чтобы более четко разграничить, где эта часть начинается и заканчивается, потому что эта намного сложнее, чем любая из предыдущих функций, запускаемых на каждом наборе данных, определенном в их функциях lapply.
Первый шаг, если вы извините за каламбур, заключается в оценке классического решения OLS для модели со всеми 30 регрессорами-кандидатами, включенными для каждого набора данных, с использованием
full_models <- lm(ds_i$Y ~ ., data = ds_i)
Поскольку функция lm() в R выполняет стандартные регрессии OLS, а синтаксис Y ~ . указывает на включение всех возможных независимых переменных. Если бы это не было установлено заранее, в аргумент object в функции step() прямо под ним не было бы ничего помещено, и, что важно, не было бы условия остановки/ верхний предел аргумента scope:
back <- stats::step(object = full_models, scope = formula(full_models),
direction = 'back', trace = FALSE) }) )
Одна вещь, которую следует упомянуть здесь с точки зрения гибкости разрешенного синтаксиса в step(), заключается в том, что вы можете установить direction равным «назад» или «назад». , это не имеет значения.
Еще одна вещь, время выполнения этого значительно выше, чем у любого из методов подгонки лассо, даже у cv.glmnet, вот сколько времени требуется только для 50 наборов данных:
> time_taken_to_fit_BE user system elapsed 0.03 0.05 12.14
Когда вы запускаете 260 000 и вам также нужно запустить Forward Selection, который занимает не так много времени, но занимает, может быть, 60–80% времени, это займет весь день и, возможно, всю ночь.
Оттуда, чтобы извлечь коэффициенты, а затем и имена выбранных по очереди переменных, используется следующий код:
# Extract the names of all IVs selected (besides their intercepts) by each
# of the 260k Backward Elimination Stepwise Regressions estimated above
# and store those estimates in the list object 'IVs_Selected_by_BE'
set.seed(11) # for reproducibility
IVs_Selected_by_BE <- lapply(seq_along(BE.fits),
\(i) names(coef(BE.fits[[i]])[-1]))
BE_IV_Candidates <- function(var_names, IVs_Selected_by_BE) {
Candidate_Vars <- setdiff(var_names, IVs_Selected_by_BE)
return(Candidate_Vars)}
# now do the same but for the candidate factors NOT selected
IVs_Not_Selected_by_BE <- lapply(seq_along(BE.fits), \(j)
j <- (BE_IV_Candidates(var_names, IVs_Selected_by_BE[[j]])))
После этого рассчитываются показатели производительности, и хотя логика остается той же, что и при запуске алгоритма 1st Benchmark, я включу сюда код, чтобы эта статья могла быть автономной еще немного. Первые две части, часть, которая вычисляет TP, TN, FP, FN, а также общее количество P и общее количество N для каждого набора данных, включены в один длинный фрагмент ниже:
### Count up how many Variables Selected match the true
### structural equation variables for that dataset in order
### to measure Backward Stepwise's performance.
# the True Positive Rate
num_of_Positives <- lapply(Structural_Variables, function(i) { length(i) })
BE_TPs_list <- lapply(seq_along(datasets), \(i)
sum(IVs_Selected_by_BE[[i]] %in%
Structural_Variables[[i]]))
BE_TPRs_list = lapply(seq_along(datasets), \(j)
j <- (BE_TPs_list[[j]]/num_of_Positives[[j]]) )
# the number of False Positives & True Negatives for each Regression
num_of_Negatives <- lapply(Structural_Variables, function(i) {30 - length(i)})
BE_FPs_list <- lapply(seq_along(datasets), \(i)
sum(!(IVs_Selected_by_BE[[i]] %in%
Structural_Variables[[i]])))
BE_TNs_list <- lapply(seq_along(datasets), \(k)
sum(IVs_Not_Selected_by_BE[[k]] %in%
Nonstructural_Variables[[k]]))
# the number of False Negatives Selected by each BE
BE_FNs_list <- lapply(seq_along(datasets), \(i)
sum(IVs_Not_Selected_by_BE[[i]] %in%
Structural_Variables[[i]]))
# the False Positive Rate = FP/(FP + TN)
BE_FPRs_list = lapply(seq_along(datasets), \(j)
j <- (BE_FPs_list[[j]])/(BE_FPs_list[[j]] + BE_TNs_list[[j]]))
BE_FPRs_list2 = lapply(seq_along(datasets), \(i)
i <- (BE_FPs_list[[i]])/(num_of_Negatives[[i]]))
# the True Negative Rate
BE_TNRs_list <- lapply(seq_along(datasets), \(j)
j <- (BE_TNs_list[[j]])/((num_of_Negatives[[j]])))
BE_TNRs_list2 <- lapply(BE_FPRs_list, \(i)
i <- (1 - i))
# the False Negative Rate = FN/(FN + TP)
BE_FNRs_list = lapply(seq_along(datasets), \(j)
j <- (BE_FNs_list[[j]])/(BE_FNs_list[[j]] + BE_TPs_list[[j]]))
## calculate the accuracy and F1 score with help from GPT 4
BE_Accuracy_list <- lapply(seq_along(datasets), function(i)
(BE_TPs_list[[i]] + BE_TNs_list[[i]])/(BE_TPs_list[[i]] + BE_TNs_list[[i]] + BE_FPs_list[[i]] + BE_FNs_list[[i]]))
# First calculate precision and TPR for each dataset
BE_Precision_list <- lapply(seq_along(datasets), function(i)
BE_TPs_list[[i]]/(BE_TPs_list[[i]] + BE_FPs_list[[i]]))
# Then calculate F1 score for each dataset
BE_F1_Score_list <- lapply(seq_along(datasets), function(i)
2 * (BE_Precision_list[[i]] * BE_TPRs_list[[i]])/(BE_Precision_list[[i]] + BE_TPRs_list[[i]]))
Я включил несколько способов расчета нескольких показателей производительности, чтобы перепроверить их достоверность с помощью альтернативных формул.
Одна вещь, которую важно понимать в отношении каждого из приведенных выше вычислений, заключается в том, что с точки зрения их классов объектов в R, в данный момент все они являются списками, поэтому мне действительно нужно привести их в более полезную форму, чтобы я мог агрегировать их, чем это. Как и в предыдущем примере в предыдущей статье, я использую функцию unlist() для достижения этой цели, я привожу только первый такой пример явно, потому что остальные идентичны:
# the True Positive Rates as a vector rather than a list BE_TPRs <- unlist(BE_TPRs_list) BE_mean_TPR <- round(mean(BE_TPRs), 3)
Насколько хорошо работало обратное исключение?
Это результаты при использовании функции step() из пакета статистики R, выбранные переменные и, следовательно, производительность с использованием функции stepAIC() из пакета MASS идентична, а время выполнения нет. Что касается различий во времени выполнения в среднем, версия пакета статистики занимает больше времени на основе сравнения времени выполнения для 6 подмножеств наборов данных по 5 или 10 000. Я также запускал их с помощью stepAIC(), чтобы убедиться, что результаты были одинаковыми, а также сравните время выполнения.
Время выполнения
Общее время работы составило 1533 минуты или 25,5 часов, что просто чудовищно!
Показатели производительности для задач классификации округлены до 3 знаков после запятой
Средняя точность = 0,859
Среднее положительное прогностическое значение, также известное как точность = 0,670.
Средний балл F1 = 0,789.
Средняя истинная положительная частота, также известная как чувствительность = 1,00.
Средний истинный отрицательный показатель, также известный как специфичность = 0,799.
Средняя частота ложноположительных результатов = 0,201.
Средняя частота ложноотрицательных результатов = 0
Оценка эффективности с точки зрения эконометрики моделирования структурными уравнениями
Количество выбранных недоопределенных спецификаций регрессии: 0
Количество выбранных спецификаций регрессии с завышенной спецификацией: 256 090.
Количество выбранных правильно указанных спецификаций: 3 910
Контрольный показатель № 3: пошаговая регрессия прямого выбора
Как вы можете себе представить, большая часть кода, используемого для запуска BE, может быть повторно использована синтаксически без изменений (просто с другими префиксами в именах), но некоторые из них требуют тонкой и точной настройки.
Самая важная часть и часть с наиболее легко запутанными различиями (как я узнал из первых рук) — это часть кода, которая подходит для всех N прямых пошаговых регрессий, она включена ниже:
### Step 3: Run a Forward Selection Stepwise Regression
### function on each of the 260,000 datasets.
### Assign the null models to their corresponding datasets
### and store these in the object "nulls".
CL <- makeCluster(detectCores() - 3L)
system.time( clusterExport(CL, c('datasets')) )
set.seed(11) # for reproducibility
time_taken_to_fit_FS <- system.time( FS.fits <- parLapply(cl = CL, X = datasets, \(ds_i) {
null_models <- lm(ds_i$Y ~ 1, data = ds_i)
full_models <- lm(ds_i$Y ~ ., data = ds_i)
forward <- stats::step(object = null_models, scope = formula(full_models),
direction = 'forward', trace = FALSE) }) )
time_taken_to_fit_FS
# Extract the elapsed time
FS_elapsed_time <- time_taken_to_fit_FS["elapsed"]
stopCluster(CL)
rm(CL)
Основными отличиями здесь в функции подгонки являются добавление еще одной модели для предварительной подгонки N раз, а именно, нулевой модели, в которую включен только перехват, но не один из 30 регрессоров-кандидатов, а затем тот факт, что эти null_models занимают место в аргументе object в функции step(), которую full_models взяли во 2-м тесте, это абсолютно важно. Я считаю полезным думать об аргументах object и scope здесь как о нижней и верхней границах пространства поиска для оптимальных спецификаций регрессии. Еще одно отличие состоит в том, что аргумент direction теперь установлен равным «вперед», но это отличие тривиально и очевидно.
После этой части все действительно то же самое, что и для пошаговой регрессии с обратным исключением, только со всеми префиксами BE, замененными префиксами FS.
Репликация результатов выбора переменных с помощью функции stepAIC() из пакета R MASS
Код, используемый для повторной подгонки обратных ступенчатых регрессий, чтобы убедиться, что результаты, с точки зрения которых предикторы-кандидаты выбраны 2-м эталоном для каждого из 260 000 синтетических выборочных наборов данных, воспроизводимы, будет показан, а некоторые различия объяснены в эта секция. На этот раз я явно включил обе версии, которые распараллеливают и не распараллеливают выполнение этих регрессий, чего я не делал раньше, несмотря на то, что для каждого теста было две версии. Я делаю это здесь, потому что я использовал оба в разное время, обычно в зависимости от того, делал ли я что-то еще на своем ноутбуке во время их запуска, тем более, когда я запускал 2-й и 3-й тесты из-за гораздо более высокого времени выполнения, чем когда я запускал 1-й тест.
Ниже приведен код, используемый для воспроизведения исходных значений, сделанных с помощью функции step(), с функцией stepAIC():
# load all necessary packages
library(plyr)
library(dplyr)
library(readr)
library(stringi)
library(purrr)
library(stats)
library(MASS)
library(leaps)
library(lars)
library(parallel)
library(data.table)
### Step 2: Run a Backward Elimination Stepwise Regression
### function on each of the 260,000 datasets again,
### but this time via the stepAIC function from the MASS package.
CL <- makeCluster(detectCores() - 3L)
system.time( clusterExport(CL, c('datasets')) )
set.seed(11) # for reproducibility
time_to_fit_BE <- system.time( BE.fits <- parLapply(cl = CL, X = datasets, \(ds_i) {
null_models <- lm(ds_i$Y ~ 1, ds_i)
full_models <- lm(ds_i$Y ~ ., data = ds_i)
back <- MASS::stepAIC(full_models, direction = 'backward',
scope = list(upper = full_models, lower = null_models),
trace = FALSE) }) )
time_to_fit_BE
stopCluster(CLs)
rm(CL)
#set.seed(11)
#time_to_fit_BE <- system.time( BE.fits <- lapply(X = datasets, \(ds_i) {
# null_models <- lm(ds_i$Y ~ 1, ds_i)
# full_models <- lm(ds_i$Y ~ ., data = ds_i)
# back <- MASS::stepAIC(full_models, direction = 'backward',
# scope = list(upper = full_models, lower = null_models),
# trace = FALSE) }) )
#time_to_fit_BE
# Extract the elapsed time
BE_elapsed_time <- time_to_fit_BE["elapsed"]
Первое, на что следует обратить внимание, это то, что в то время как с функцией step() только Forward Stepwise требуется дополнительный аргумент null_models для функции подбора Backward Stepwise, с stepAIC( ), оба обязательны. Единственное другое существенное различие в синтаксисе на этот раз — это разница между
scope = Formula(full_models) и
scope = list(upper = full_models, Lower = null_models)
Следующие строки кода, используемые для извлечения и присвоения оценочных коэффициентов списку, за которыми следует список независимых переменных, выбранных при запуске BE с использованием этого альтернативного пути в R, одинаковы, поэтому их не нужно здесь использовать. То же самое касается того, как рассчитывались показатели производительности классификации и как они использовались для подсчета количества выбранных регрессий с недостаточным, избыточным и правильным указанием.
На этом мой подробный обзор и объяснение моего кода R, используемого для запуска алгоритмов Benchmark в этом исследовательском проекте, завершаются.
Насколько хорошо работает предварительный отбор?
Время выполнения
Общее время выполнения составило 1178 минут или 19,6 часа, что опять же, столько времени, что оно граничит с моральным отвращением.
Показатели производительности для задач классификации округлены до 3 знаков после запятой
Средняя точность = 0,874
Средняя положительная прогностическая ценность, также известная как точность = 0,680.
Средний балл F1 = 0,798.
Средняя истинная положительная частота, также известная как чувствительность = 1,00.
Средний истинный отрицательный показатель, также известный как специфичность = 0,823.
Средняя частота ложноположительных результатов = 0,177.
Средняя частота ложноотрицательных результатов = 0
Оценка эффективности с точки зрения эконометрики моделирования структурными уравнениями
Количество выбранных недоопределенных спецификаций регрессии: 0
Количество выбранных спецификаций регрессии с завышенной спецификацией: 226 027.
Количество выбранных правильно указанных спецификаций: 3 973
Количество регрессий, в которых пропущен хотя бы один структурный фактор: 0
Количество включенных регрессий хотя бы с одним посторонним фактором: 226 027.
p.s.
Если вы на самом деле пытаетесь полностью воспроизвести какой-либо из тестовых методов на всех 260 000 используемых исходных наборах данных, отправьте мне электронное письмо по адресу [email protected], и я отвечу со ссылкой на папку в моем OneDrive.
