Сравнение двух разных методов регрессии
1. Введение
COVID-19 — это заболевание, вызываемое коронавирусом тяжелого острого респираторного синдрома 2 (SARS-CoV-2). COVID-19 может вызывать нарушения дыхательной системы, начиная от легких симптомов, таких как грипп, и заканчивая легочными инфекциями, такими как пневмония (ВОЗ, 2020).
У большинства инфицированных вирусом разовьется респираторное заболевание легкой или средней степени тяжести, и они выздоровеют без необходимости специального лечения. Некоторые, однако, станут очень больными и потребуют медицинской помощи. Когда зараженный человек кашляет, чихает, говорит, поет или дышит, вирус может распространяться через микроскопические частицы жидкости изо рта или носа. Размер этих частиц варьируется от крупных респираторных капель до крошечных аэрозолей. Если вы чувствуете себя плохо, очень важно соблюдать респираторный этикет, например, кашлять в согнутый локоть, оставаться дома и самоизолироваться до выздоровления.
В этой статье мы обсудим развитие COVID-19 в Индонезии, начиная с января 2020 года, когда произошел первый случай COVID-19, до июля 2021 года, когда PPKM закончился и перешел на Emergency PPKM. В этой статье данные будут обрабатываться с использованием языка программирования R в качестве среды обработки данных.
В этой статье мы рассмотрим, какие факторы способствовали росту числа новых смертей от COVID-19 в Индонезии в период с января 2020 года по июль 2021 года. Этот набор данных будет использоваться в этой статье. Набор данных, использованный в этой статье, собран открыто с сайтов covid19.go.id, Kemendagri.go.id, bps.go.id и bnpb-inacovid19.hub.arcgis.com. А в качестве среды обработки данных для обсуждения мы будем использовать язык программирования R.
Используя штрафную регрессию, в частности, регрессию Риджа и Лассо, мы проанализируем, какие факторы влияют на увеличение числа новых случаев смерти в Индонезии. где есть этапы выполнения моделирования, а именно:
- В начале мы подготовим данные;
- Затем мы разделяем данные;
- Кроме того, мы создали модель для обучения данных;
- Кроме того, выбор лучшей модели с проверкой данных;
- После этого мы проверили данные;
- И последний шаг — сравнение созданной нами модели с другими моделями и ее интерпретация.
Прежде чем мы начнем раздел кодирования, мы должны импортировать набор данных, который мы будем использовать позже.
#Import The Dataframe
#I already stored the dataset on my local
df <- read_csv("C:/Users/User/Documents/R/covid_19_indonesia_time_series_all.csv")
df
Набор данных, который мы хотим использовать, называется «covid_19_indonesia_time_series_all». Этот набор данных был собран с 8 января 2020 г. по 9 июля 2021 г. Он содержит 16 283 строки и 41 столбец.
2. Подготовка данных
Следующим этапом будет манипулирование и оценка набора данных, который у нас есть. Чтобы избежать посторонних данных и неправильного анализа наших данных, мы должны сначала очистить данные.
#Data Cleaning
drop_cols <- c('Name',
'Item',
'Kind',
'Hidden',
'City or Regency',
'Special Status',
'Location ISO Code',
'Province',
'Country',
'Continent',
'Total Regencies',
'Total Cities',
'Total Districts',
'Total Cases per Million',
'New Deaths per Million',
'Total Deaths per Million',
'Case Fatality Rate',
'Case Recovered Rate',
'Growth Factor of New Cases',
'Growth Factor of New Deaths',
'Time Zone',
'Total Urban Villages',
'Total Rural Villages',
'New Cases per Million'
)
df <- df %>% select(-drop_cols)
df %>% str()
Первым шагом в очистке данных в этом наборе данных является удаление столбцов, которые создают неопределенность и больше не нужны в нашем подходе.
Следующим шагом является замена любых неоднозначных значений, например, в столбцах «Местоположение» и «Остров». В столбце «Местоположение» есть значение «Индонезия», которое несколько сбивает с толку с точки зрения имени, но значение, которое идет с ним, все же может быть принято во внимание, поэтому мы пытаемся преобразовать значение «Индонезия» в «Неизвестно». в таком случае. Затем в столбце «Остров» есть значение NA. Поскольку значения NA в пределах «Остров» потенциально представляют неизвестный остров, мы решили заполнить значения NA значением «Неизвестно». Поскольку эти два столбца и их значения все еще можно исследовать. , мы решили не удалять строки «Индонезия» или «NA», а вместо этого изменить их на «Неизвестно», чтобы мы могли продолжить их оценку.
df$Location<-replace(df$Location,df$Location=='Indonesia','Unknown') df$Island <- ifelse(is.na(df$Island),'Unknown', df$Island) df %>% view()
Затем мы проверяем, осталось ли какое-либо значение NA в нашем наборе данных.
indx <- apply(df, 2, function(x) any(is.na(x) | is.infinite(x))) indx
Оказывается, в нашем наборе данных нет пропущенного значения. Затем мы проверяем, есть ли в нашем наборе данных дубликаты данных.
duplicated(df) sum(duplicated(df)) #There is no duplicates data in our dataframe
В результате в нашем наборе данных нет дубликатов данных. С этого момента наши данные прямо сейчас готовы к использованию для манипуляций, потому что теперь они чистые.
3. Исследовательский анализ данных
Какова цель ЭДА? EDA используется для изучения идей, которые существуют в наборах данных, которыми мы владеем. С сильным EDA вы можете лучше понять проблему, выполнить соответствующую подготовку данных и моделирование, а также узнать, как оценить ваше решение для данных, и все это может помочь вам решить вашу проблему. Вот почему нам нужно сначала выполнить EDA, прежде чем мы создадим модель. Для получения дополнительных разъяснений относительно EDA об этом наборе данных вы можете проверить эту статью здесь.
Во-первых, мы хотим знать, сколько всего случаев Covid-19 было зарегистрировано в Индонезии до настоящего времени (9 июля 2021 г.). Таким образом, мы могли бы получить первоначальный обзор набора данных, который у нас есть.
Оказывается, в Джакарте больше всего случаев по сравнению с любым другим местом. Возможно, это связано с тем, что Джакарта находится недалеко от Депока, места первого случая заболевания COVID-19 в Индонезии (Горбиано, 2020). И Джакарта имеет наибольшую мобилизацию в регионе ДЖАБОДЕТАБЕК (Джакарта, Богор, Депок, Тангеранг, Бекаси), потому что Джакарта является центром бизнеса, а также правительства в Индонезии.
BREAKING: Джокови объявляет о первых двух подтвержденных случаях заболевания COVID-19 в Индонезии
Президент Джоко «Джокови Видодо объявил, что два индонезийца дали положительный результат на новый коронавирус… www.thejakartapost.com»
Далее, мы хотели бы знать, как продвигается дело в каждой области, в частности, в четырех местах с самым высоким уровнем заболеваемости в Индонезии: «DKI Джакарта», «Джава Барат (Западная Ява)», «Джава Тенгах (Центральная Ява)» и «Джава Тимур (Восточная Ява)».
Когда мы смотрим на график четырех основных мест с точки зрения общего числа новых случаев Covid-19 в районе Явы, мы видим, что все они имели примерно одинаковое пиковое время, а именно в начале года. в декабре и январе, а затем наибольшие в июне-июле. Однако, если мы присмотримся, то увидим, что общее число новых случаев на Западной Яве выше, чем в Джакарте в начале 2021 года, и по мере приближения праздника Курбан-байрам (Ид — большой праздник для мусульман), в конце мая — начале июня общее число случаев заболевания на Западной Яве является самым высоким в Индонезии, хотя к концу Джакарта превзошла общее число новых случаев (Nurita, 2021).
Мы замечаем закономерность после просмотра данных по новым случаям, но не можем понять, связано ли это с политикой правительства, например, с политикой «ППКМ», которая привела к уменьшению числа случаев в конце 2020 года, или с чем-то еще. Нам нужно исследовать глубже, и в этом случае мы должны различать даты, когда применялся ПКМ.
Мы пытаемся идентифицировать между PPKM, Emergency PPKM, PSBB и No Policy (отсутствие политики происходит, когда правительство вообще не проводило никакой политики). Мы можем видеть из кода ниже:
df$PPKM = ifelse(df$Date >= "2021-01-11" & df$Date <= "2021-07-02", "PPKM", ifelse(df$Date >= "2021-07-03","PPKM Darurat", ifelse(df$Date >= "2020-04-17" & df$Date <= "2020-12-07","PSBB","No Policy"))) df
- Политика не применяется с января 2020 г. по 16 апреля 2020 г., а затем повторяется с 8 декабря 2020 г. по 10 января 2021 г.;
- Политика PSBB действует с 17 апреля 2020 г. по 7 декабря 2020 г.;
- Политика PPKM действует с 11 января 2021 г. по 2 июля 2021 г.;
- Аварийный PPKM происходит с 3 июля 2021 г. по 20 июля 2021 г. (но он не будет полностью отображаться в этом наборе данных, потому что этот набор данных заканчивается 9 июля 2021 г., проверьте здесь).
Затем мы создаем новый столбец для нашего набора данных и показываем график:
Согласно линейному графику выше, когда правительство Индонезии приняло политику по внедрению PPKM, которая была реализована 3 января 2021 г., PPKM удалось сократить количество случаев из-за отсутствия политики PSBB, которая закончилась 7 декабря 2020 г., но Успех PPKM не удался из-за того, что многие решили взять праздник Ид. Администрация внедрила экстренный PPKM в результате увеличения числа новых случаев во время каникул Ид (но в этом наборе данных еще не видно, удалось ли экстренному PPKM сократить количество случаев Covid-19).
Из этого EDA мы можем сделать следующие выводы.
- Общее количество новых случаев на Яве, особенно в четырех самых высоких случаях (Джакарта, Западная Ява, Центральная Ява и Восточная Ява), является самым большим из всех мест;
- В Джакарте зарегистрировано наибольшее количество новых случаев среди четырех самых высоких случаев, возможно, потому, что она находится недалеко от того места, где произошел первый случай заболевания COVID-19, а также имеет самый высокий уровень мобильности в Индонезии, потому что это центр бизнеса и правительства в Индонезия;
- PPKM — наиболее эффективная политика, которую правительство внедрило в связи с пандемией COVID-19, даже несмотря на то, что когда она приближается к празднику Курбан-Байрам, количество новых случаев резко увеличилось (что также требует дальнейшего изучения);
- Мы можем предположить, что внедрение PSBB и PPKM в Индонезии привело к сокращению числа случаев, особенно PPKM, и кажется, что правительству Индонезии следует проводить ограничительную политику.
4. Новое прогнозирование смертей с использованием штрафной регрессии
Теперь мы хотели бы предсказать новые случаи смерти от COVID-19 в Индонезии, используя регрессию со штрафом, в частности регрессию хребта и регрессию Лассо. Как упоминалось выше, для этого нужно выполнить несколько шагов, начиная с подготовки данных и заканчивая сравнением и интерпретацией созданной нами модели.
4.1 Подготовка данных
Во-первых, нам нужно подготовить данные, которые мы хотим использовать, выбрав необходимые функции, которые мы хотим использовать.
drop_cols_Prediction <- c('Date',
'Location',
'Location Level',
'Island',
'Area (km2)',
'Population Density',
'Longitude',
'Latitude',
'year',
'month',
'day'
)
new_df <- df %>% select(-drop_cols_Prediction)
new_df %>% view() %>% str()
Почему мы должны отказаться от всех этих столбцов? Это потому, что мы хотим избежать ненужных данных для манипуляций, а также R не может запускать множество функций, что приведет к медленной работе программы.
Следующим шагом является однократное кодирование. Какой смысл использовать однократное кодирование? Это довольно просто, поэтому мы можем наблюдать, какая политика правительства Индонезии вызывает увеличение числа новых случаев смерти.
library(caret)
dummy <- dummyVars(" ~ .", data=new_df)
new_df <- data.frame(predict(dummy, newdata = new_df))
new_df
colnames(new_df) <- c('New_Cases','New_Deaths','New_Recovered','New_Active_Cases','Total_Cases','Total_Deaths','Total_Recovered','Total_Active_Cases','Population','PPKM_No_Policy', 'PPKM_PPKM', 'PPKM_PPKM_Darurat', 'PPKM_PSBB')
new_df %>% str()
И результат:
'data.frame': 16283 obs. of 13 variables: $ New_Cases : num 9 0 0 0 0 0 0 0 0 0 ... $ New_Deaths : num 3 1 1 3 3 0 0 0 0 0 ... $ New_Recovered : num 1 23 14 8 1 0 0 0 0 4 ... $ New_Active_Cases : num 5 -24 -15 -11 -4 0 0 0 0 -4 ... $ Total_Cases : num 9 9 9 9 9 9 9 9 9 9 ... $ Total_Deaths : num 5 6 7 10 13 13 13 13 13 13 ... $ Total_Recovered : num 10 33 47 55 56 56 56 56 56 60 ... $ Total_Active_Cases: num -6 -30 -45 -56 -60 -60 -60 -60 -60 -64 ... $ Population : num 40479023 40479023 40479023 40479023 40479023 ... $ PPKM_No_Policy : num 1 1 1 1 1 1 1 1 1 1 ... $ PPKM_PPKM : num 0 0 0 0 0 0 0 0 0 0 ... $ PPKM_PPKM_Darurat : num 0 0 0 0 0 0 0 0 0 0 ... $ PPKM_PSBB : num 0 0 0 0 0 0 0 0 0 0 ...
4.2 Разделить данные
Нам нужно разделить данные на три части: обучающие, проверочные и тестовые данные. Используя эти коды:
#Prediction #Split the data to Train, Validation, and Test data set.seed(123) sample <- sample.split(new_df$New_Deaths, SplitRatio = .80) pre_train_df <- subset(new_df, sample == TRUE) sample_train <- sample.split(pre_train_df$New_Deaths, SplitRatio = .80) train_df <- subset(pre_train_df, sample_train == TRUE) validation_df <- subset(pre_train_df, sample_train == FALSE) test_df <- subset(new_df, sample == FALSE)
4.3 Выбор функций
Чтобы предотвратить мультиколлинеарность, мы можем нарисовать расширенный парный график поверх обучающих данных. С помощью функции pairs.panels() из библиотеки psych.
#Feature Selection
#Correlation study
library(psych)
pairs.panels(pre_train_df,
method = "pearson", # correlation method
hist.col = "#00AFBB",
density = TRUE, # show density plots
ellipses = TRUE # show correlation ellipses
)
Мы обнаружили, что корреляция между «новые случаи»-«новые выздоровевшие», «новые случаи»-«всего случаев», «новые случаи»-«всего смертей», «новые случаи»-«всего выздоровевших» и «новые случаи» -Total Active Cases» очень высоки (чтобы определить, является ли функция высокой корреляцией или нет, увидев значение абсолютной корреляции. Если есть значение более 80%, мы можем просто сказать, что корреляция высока). Мы решили включить только «Новые случаи» для представления функций, а также потому, что «Новые случаи» наиболее коррелируют с целевой переменной, которой является «Новые смерти».
На следующем шаге нам нужно удалить столбцы «Новые выздоровевшие», «Всего случаев», «Всего смертей», «Всего выздоровевших» и «Всего активных случаев».
# drop correlated columns
library(dplyr)
drop_cols_sample <- c('New_Recovered',
'Total_Cases',
'Total_Deaths',
'Total_Recovered',
'Total_Active_Cases')
train_df <- train_df %>% select(-drop_cols_sample)
validation_df <- validation_df %>% select(-drop_cols_sample)
test_df <- test_df %>% select(-drop_cols_sample)
4.4 Ридж-регрессия и регрессия Лассо
Что такое регрессия Риджа и Лассо (оператор наименьшего абсолютного сокращения и выбора)? Ридж-регрессия и регрессия Лассо — это линейные регрессии с измененной функцией потерь.
Если лямбда является неотрицательным параметром регуляризации, чем больше лямбда, тем сильнее регуляризуется регрессия. Если лямбда равна 0, регрессия является обычной линейной регрессией. Но, разница между Ридж-регрессиями и ЛАССО-регрессиями, коэффициенты в Ридж-регрессии никогда не бывают равны нулю, однако в Лассо некоторые из них будут точно равны нулю (другими словами, исключены из модели).
4.4.1 Регрессия хребта
# feature preprocessing # to ensure we handle categorical features x <- model.matrix(New_Deaths ~ ., train_df)[,-1] y <- train_df$New_Deaths # ridge regression # fit multiple ridge regression with different lambda # lambda = [0.01, 0.1, 1, 10] ridge_reg_pointzeroone <- glmnet(x, y, alpha = 0, lambda = 0.01) coef(ridge_reg_pointzeroone) ridge_reg_pointone <- glmnet(x, y, alpha = 0, lambda = 0.1) coef(ridge_reg_pointone) ridge_reg_one <- glmnet(x, y, alpha = 0, lambda = 1) coef(ridge_reg_pointone) ridge_reg_ten <- glmnet(x, y, alpha = 0, lambda = 10) coef(ridge_reg_ten)
# comparison on validation data # to choose the best lambda # Make predictions on the validation data x_validation <- model.matrix(New_Deaths ~ ., validation_df)[,-1] y_validation <- validation_df$New_Deaths RMSE_ridge_pointzeroone <- sqrt(mean((y_validation - predict(ridge_reg_pointzeroone, x_validation))^2)) RMSE_ridge_pointzeroone #8.252002 RMSE_ridge_pointone <- sqrt(mean((y_validation - predict(ridge_reg_pointone, x_validation))^2)) RMSE_ridge_pointone #8.242678 RMSE_ridge_one <- sqrt(mean((y_validation - predict(ridge_reg_one, x_validation))^2)) RMSE_ridge_one #8.223498 <- best RMSE_ridge_ten <- sqrt(mean((y_validation - predict(ridge_reg_ten, x_validation))^2)) RMSE_ridge_ten #8.682836
# true evaluation on test data # using the best model --> RMSE_ridge_one x_test <- model.matrix(New_Deaths ~ ., test_df)[,-1] y_test <- test_df$New_Deaths # RMSE RMSE_ridge_best <- sqrt(mean((y_test - predict(ridge_reg_one, x_test))^2)) RMSE_ridge_best # 0.06423617 # MAE MAE_ridge_best <- mean(abs(y_test-predict(ridge_reg_one, x_test))) MAE_ridge_best # MAPE MAPE_ridge_best <- mean(abs((predict(ridge_reg_one, x_test) - y_test))/y_test) MAPE_ridge_best
Для MAE это означает, что в среднем наш прогноз отклоняется от истинного accept_prob на 3,3, а также для RMSE. Стандартное отклонение ошибок прогноза составляет 8,7, т.е. от линии регрессии, остатки в основном отклоняются в пределах +- 8,7. К сожалению, мы не можем определить MAPE, потому что в этих данных есть абсолютное нулевое значение.
4.4.2 Лассо-регрессия
# feature preprocessing # to ensure we handle categorical features x <- model.matrix(New_Deaths ~ ., train_df)[,-1] y <- train_df$New_Deaths # lasso regression # fit multiple lasso regression with different lambda # lambda = [0.01, 0.1, 1, 10] lasso_reg_pointzeroone <- glmnet(x, y, alpha = 1, lambda = 0.01) coef(lasso_reg_pointzeroone) lasso_reg_pointone <- glmnet(x, y, alpha = 1, lambda = 0.1) coef(lasso_reg_pointone) lasso_reg_one <- glmnet(x, y, alpha = 1, lambda = 1) coef(lasso_reg_pointone) lasso_reg_ten <- glmnet(x, y, alpha = 1, lambda = 10) coef(lasso_reg_ten)
# comparison on validation data # to choose the best lambda # Make predictions on the validation data RMSE_lasso_pointzeroone <- sqrt(mean((y_validation - predict(lasso_reg_pointzeroone, x_validation))^2)) RMSE_lasso_pointzeroone # 8.249907 RMSE_lasso_pointone <- sqrt(mean((y_validation - predict(lasso_reg_pointone, x_validation))^2)) RMSE_lasso_pointone # 8.223708 RMSE_lasso_one <- sqrt(mean((y_validation - predict(lasso_reg_one, x_validation))^2)) RMSE_lasso_one # 8.153449 <- Best RMSE_lasso_ten <- sqrt(mean((y_validation - predict(lasso_reg_ten, x_validation))^2)) RMSE_lasso_ten # 9.12958
# true evaluation on test data # using the best model --> RMSE_lasso_pointone RMSE_lasso_best <- sqrt(mean((y_test - predict(lasso_reg_one, x_test))^2)) RMSE_lasso_best # 8.62048 # MAE MAE_lasso_best <- mean(abs(y_test-predict(lasso_reg_one, x_test))) MAE_lasso_best # MAPE MAPE_lasso_best <- mean(abs((predict(lasso_reg_one, x_test) - y_test))/y_test) MAPE_laso_best
Для MAE это означает, что в среднем наш прогноз отклоняется от истинного accept_prob на 3,3, а также для RMSE. Стандартное отклонение ошибок прогноза составляет 8,6, т.е. от линии регрессии, остатки в основном отклоняются в пределах +- 8,6. К сожалению, мы можем определить MAPE, потому что эти данные имеют нулевое значение.
5. В заключение
Поскольку в данных есть нулевое значение, значение MAPE для Ridge и Lasso соответственно не может быть определено.
Основываясь на этом прогнозе модели с регрессией регуляризации, модель, которую мы используем для прогнозирования новых смертей, выглядит следующим образом:
Итак, модель, которую мы выбираем для предсказания новых смертей, можно интерпретировать следующим образом:
Новые случаи смерти = -0,21 + (0,023 новых случая) + (-0,0074 новых активных случая) + (0,000000101 населения) + (-0,325 без политики) + (0,242 млн. млн) + (7,2 млн. млн Дарурат)
Если мы хотим определить прогноз, какой фактор оказывает наибольшее влияние, вы можете выбрать между Ridge и Lasso в качестве наиболее подходящего метода для прогнозирования того, какой фактор имеет наибольшее влияние. Результат сравнения между Ridge и Lasso выше заключается в том, что мы выбираем регрессию гребня, потому что ее значение RMSE больше, чем значение RMSE регрессии Lasso. Это означает, что регрессия хребта является наиболее подходящим методом для прогнозирования в этом наборе данных.
И из выбранной модели можно интерпретировать, что увеличение новых случаев смерти - для случая в Индонезии - PPKM увеличит вероятность увеличения числа смертей на 0,242, если другие характеристики будут зафиксированы. Почему бы не взять Аварийный ППКМ? Здесь мы не делаем такого вывода, поскольку данные по аварийному ПКМ слишком ранние и неполные.
Об авторе
Галих Триас Фароби верит в тяжелую и умную работу. Это человек, который постоянно стремится узнавать что-то новое. В настоящее время он расширяет свои знания в области науки о данных и рад поделиться со всеми как техническими, так и нетехническими аспектами науки о данных.
Он твердо верит в развитие мышления роста и продолжение учебы на протяжении всей жизни. Пожалуйста, не стесняйтесь общаться; давайте учиться вместе!
Ссылки
[1] Farobi, G. T. EDA и визуализация случаев COVID-19 в Индонезии в период с января 2020 г. по июль 2021 г. с использованием R (2021). https://medium.com/@farobiobi/eda-and-visualization-of-indonesias-covid-19-cases-in-period-january-2020-to-july-2021-using-r-30c9d1b5ea62
[2] Горбиано. BREAKING: Jokowi объявляет о первых двух подтвержденных случаях COVID-19 в Индонезии (2020 г.). https://www.thejakartapost.com/news/2020/03/02/breaking-jokowi-announces-indonesias-first-two-confirmed-covid-19-cases.html
[3] Хендратно. Набор данных по COVID-19 в Индонезии (2021 г.). https://www.kaggle.com/hendratno/covid19-indonesia
[4] Нурита, Деви. Случаи Covid-19 в Кудусе резко выросли на 7,594% после праздника Ид, 53,4% по стране (2021 г.). https://en.tempo.co/read/1470785/covid-19-cases-in-kudus-skyrocket-7-594-after-eid-holiday-53-4-national
[5] Парамадита, Анджарсари. PPKM Darurat необходим для устойчивого восстановления в Индонезии (2021 г.). https://www.thejakartapost.com/news/2021/07/06/ppkm-darurat-necessary-for-sustainable-recovery-in-indonesia.html
[6] Сухенда, Дио. Индонезия сталкивается с всплеском COVID-19 после Idul Fitri из Центральной Явы, Восточная Ява (2021 г.). https://www.thejakartapost.com/news/2021/06/08/indonesia-faces-post-idul-fitri-covid-19-surge-from-central-java-east-java.html
[7] ВОЗ. Коронавирусная болезнь (COVID-19) (2020 г.). https://www.who.int/health-topics/coronavirus#tab=tab_1
