PCA to Reduce Multicollinearity
1 Setup
library(dplyr) # data manipulation
library(GGally) # cek korelasi
library(performance) # cek asumsi model2 Business Question
- Tujuan: Membuat model regresi untuk prediksi crime rate.
- Fungsi PCA: Pemenuhan asumsi no multicollinearity antar prediktor.
Catatan: pada contoh berikut tidak dilakukan train-test splitting. Pada kenyataannya PCA hanya dilakukan pada data train saja, kemudian data test akan ditransformasi menggunakan eigen vector yang telah dipelajari dari data train.
3 Read Data
Data crime merupakan data sosio-ekonomi dari 47 negara
bagian US di tahun 1960. Nilai pada variable tampaknya telah
diskalakan.
crime <- read.csv("crime.csv")
glimpse(crime)## Rows: 47
## Columns: 16
## $ percent_m <int> 151, 143, 142, 136, 141, 121, 127, 131, 157, 140,~
## $ is_south <int> 1, 0, 1, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 0~
## $ mean_education <int> 91, 113, 89, 121, 121, 110, 111, 109, 90, 118, 10~
## $ police_exp60 <int> 58, 103, 45, 149, 109, 118, 82, 115, 65, 71, 121,~
## $ police_exp59 <int> 56, 95, 44, 141, 101, 115, 79, 109, 62, 68, 116, ~
## $ labour_participation <int> 510, 583, 533, 577, 591, 547, 519, 542, 553, 632,~
## $ m_per1000f <int> 950, 1012, 969, 994, 985, 964, 982, 969, 955, 102~
## $ state_pop <int> 33, 13, 18, 157, 18, 25, 4, 50, 39, 7, 101, 47, 2~
## $ nonwhites_per1000 <int> 301, 102, 219, 80, 30, 44, 139, 179, 286, 15, 106~
## $ unemploy_m24 <int> 108, 96, 94, 102, 91, 84, 97, 79, 81, 100, 77, 83~
## $ unemploy_m39 <int> 41, 36, 33, 39, 20, 29, 38, 35, 28, 24, 35, 31, 2~
## $ gdp <int> 394, 557, 318, 673, 578, 689, 620, 472, 421, 526,~
## $ inequality <int> 261, 194, 250, 167, 174, 126, 168, 206, 239, 174,~
## $ prob_prison <dbl> 0.084602, 0.029599, 0.083401, 0.015801, 0.041399,~
## $ time_prison <dbl> 26.2011, 25.2999, 24.3006, 29.9012, 21.2998, 20.9~
## $ crime_rate <int> 791, 1635, 578, 1969, 1234, 682, 963, 1555, 856, ~
Deskripsi kolom:
percent_m: percentage of males aged 14-24is_south: whether it is in a Southern state. 1 for Yes, 0 for No.mean_education: mean years of schoolingpolice_exp60: police expenditure in 1960police_exp59: police expenditure in 1959labour_participation: labour force participation ratem_per1000f: number of males per 1000 femalesstate_pop: state populationnonwhites_per1000: number of non-whites resident per 1000 peopleunemploy_m24: unemployment rate of urban males aged 14-24unemploy_m39: unemployment rate of urban males aged 35-39gdp: gross domestic product per headinequality: income inequalityprob_prison: probability of imprisonmenttime_prison: average time served in prisonscrime_rate: crime rate in an unspecified category
4 Exploratory Data Analysis
Cek korelasi antar kolom numerik.
ggcorr(crime, label = T, layout.exp = 3, hjust = 1)Terdapat korelasi yang cukup kuat antar prediktor (selain kolom
crime_rateyang akan menjadi kolom target) terutama padapolice_exp59danpolice_exp60.
5 Data Pre-processing: PCA
Melakukan PCA hanya pada kolom prediktor.
crime_pca <- prcomp(select(crime, -crime_rate), # hapus kolom target
scale = T) # melakukan scaling data
summary(crime_pca)$importance %>%
as.data.frame()6 Modelling
Membandingkan kedua model:
model_ori: Model dengan menggunakan semua variabel predictor dari datacrimeawalmodel_pca: Model dengan menggunakan semua PC sebagai predictor
6.1
model_ori
model_ori <- lm(crime_rate ~ ., data = crime)
summary(model_ori)##
## Call:
## lm(formula = crime_rate ~ ., data = crime)
##
## Residuals:
## Min 1Q Median 3Q Max
## -395.74 -98.09 -6.69 112.99 512.67
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -5984.2876 1628.3184 -3.675 0.000893 ***
## percent_m 8.7830 4.1714 2.106 0.043443 *
## is_south -3.8035 148.7551 -0.026 0.979765
## mean_education 18.8324 6.2088 3.033 0.004861 **
## police_exp60 19.2804 10.6110 1.817 0.078892 .
## police_exp59 -10.9422 11.7478 -0.931 0.358830
## labour_participation -0.6638 1.4697 -0.452 0.654654
## m_per1000f 1.7407 2.0354 0.855 0.398995
## state_pop -0.7330 1.2896 -0.568 0.573845
## nonwhites_per1000 0.4204 0.6481 0.649 0.521279
## unemploy_m24 -5.8271 4.2103 -1.384 0.176238
## unemploy_m39 16.7800 8.2336 2.038 0.050161 .
## gdp 0.9617 1.0367 0.928 0.360754
## inequality 7.0672 2.2717 3.111 0.003983 **
## prob_prison -4855.2658 2272.3746 -2.137 0.040627 *
## time_prison -3.4790 7.1653 -0.486 0.630708
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 209.1 on 31 degrees of freedom
## Multiple R-squared: 0.8031, Adjusted R-squared: 0.7078
## F-statistic: 8.429 on 15 and 31 DF, p-value: 0.0000003539
Setiap nilai Estimate yang ditampilkan dapat diinterpretasi karena merepresentasikan nilai aslinya.
6.2
model_pca
# menggunakan nilai PC, lalu gabung dengan kolom target
crime_data_pc <- crime_pca$x %>%
as.data.frame() %>%
mutate(crime_rate = crime$crime_rate)
head(crime_data_pc)# buat model
model_pca <- lm(crime_rate ~ ., data = crime_data_pc)
summary(model_pca)##
## Call:
## lm(formula = crime_rate ~ ., data = crime_data_pc)
##
## Residuals:
## Min 1Q Median 3Q Max
## -395.74 -98.09 -6.69 112.99 512.67
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 905.09 30.50 29.680 < 0.0000000000000002 ***
## PC1 65.22 12.56 5.191 0.0000124368 ***
## PC2 -70.08 18.42 -3.806 0.000625 ***
## PC3 25.19 21.77 1.157 0.255987
## PC4 69.45 28.59 2.429 0.021143 *
## PC5 -229.04 31.49 -7.274 0.0000000349 ***
## PC6 -60.21 41.44 -1.453 0.156305
## PC7 117.26 54.34 2.158 0.038794 *
## PC8 28.72 55.60 0.517 0.609159
## PC9 -37.18 63.57 -0.585 0.562890
## PC10 56.32 68.95 0.817 0.420261
## PC11 30.59 73.54 0.416 0.680272
## PC12 289.61 86.09 3.364 0.002059 **
## PC13 81.79 117.06 0.699 0.489962
## PC14 219.19 127.48 1.719 0.095517 .
## PC15 -622.21 453.79 -1.371 0.180174
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 209.1 on 31 degrees of freedom
## Multiple R-squared: 0.8031, Adjusted R-squared: 0.7078
## F-statistic: 8.429 on 15 and 31 DF, p-value: 0.0000003539
Setiap nilai Estimate yang ditampilkan tidak dapat diinterpretasi karena nilai PC adalah nilai data asli yang sudah tercampur dengan formula tertentu (menggunakan eigen vector).
7 Evaluasi Model
compare_performance(model_ori, model_pca) %>%
as.data.frame()Performa model
model_oridanmodel_pcasama persis. Hal ini karena padamodel_pcakita menggunakan keseluruhan informasi pada data awalcrime, sehingga hasilnya sama saja.
8 Uji Asumsi
Asumsi linear model agar dikatakan Best Linear Unbiased Estimate:
- Linearity of Predictor
- Normality of Residual
- Homoscedasticity of Residual
- No multicollinearity
8.1
model_ori
performance::check_model(model_ori)performance::check_normality(model_ori)## OK: residuals appear as normally distributed (p = 0.779).
performance::check_heteroscedasticity(model_ori)## OK: Error variance appears to be homoscedastic (p = 0.346).
performance::check_collinearity(model_ori) %>%
as.data.frame() %>%
arrange(-VIF)Pada
model_ori, asumsi dari multicolinearity tidak terpenuhi karena ada nilai VIF yang diatas 10. Kolompolice_exp59danpolice_exp60memiliki korelasi yang sangat kuat.
8.2
model_pca
performance::check_model(model_pca)performance::check_normality(model_pca)## OK: residuals appear as normally distributed (p = 0.779).
performance::check_heteroscedasticity(model_pca)## OK: Error variance appears to be homoscedastic (p = 0.346).
performance::check_collinearity(model_pca) %>%
as.data.frame() %>%
arrange(-VIF)Pada
model_pca, asumsi no multicollinearity menjadi terpenuhi. Nilai VIF 1 menandakan semua PC yang digunakan tidak memiliki korelasi sama sekali.
9 Conclusion
Metode PCA dapat mengatasi permasalahan multicollinearity pada model linear. Namun model regresi linear kehilangan kemampuan untuk diinterpretasi, karena nilai PC sudah bercampur.