INTRODUCTION

    Data yang digunakan saat ini berkaitan dengan rumah-rumah yang ditemukan di distrik California tertentu dan beberapa ringkasan statistik tentang rumah-rumah tersebut berdasarkan data sensus tahun 1990. Dari data ini juga kita ingin mengetahui hubungan antar variable populasi (Jumlah total orang yang tinggal dalam satu blok) terhadap variable yag lainnya.

link : https://www.kaggle.com/datasets/camnugent/california-housing-prices

IMPORT LIBRARY

library(GGally)
library(dplyr)
library(lubridate)
library(MLmetrics)
library(data.table)
library(performance)
library(scales)
library(ggplot2)
library(lmtest)
library(plotly)
library(car)
library(caret)
library(tidyr)
library(tidyverse)
library(nortest)

READ DATA

Hal pertama yang harus kita lakukan yaitu import data kita terlebih dahulu read.csv.

housing <- read.csv("housing.csv")

DATA INSPECTIONS

Kita gunakan head() untuk melihat urutan data teratas dan gunakan tail() untuk melihat urutan data paling bawah.

head(housing)
tail(housing)

Cek struktur tipe data yang kita miliki menggunakan glimpse()

glimpse(housing)
#> Rows: 20,640
#> Columns: 10
#> $ longitude          <dbl> -122.23, -122.22, -122.24, -122.25, -122.25, -122.2…
#> $ latitude           <dbl> 37.88, 37.86, 37.85, 37.85, 37.85, 37.85, 37.84, 37…
#> $ housing_median_age <int> 41, 21, 52, 52, 52, 52, 52, 52, 42, 52, 52, 52, 52,…
#> $ total_rooms        <int> 880, 7099, 1467, 1274, 1627, 919, 2535, 3104, 2555,…
#> $ total_bedrooms     <int> 129, 1106, 190, 235, 280, 213, 489, 687, 665, 707, …
#> $ population         <int> 322, 2401, 496, 558, 565, 413, 1094, 1157, 1206, 15…
#> $ households         <int> 126, 1138, 177, 219, 259, 193, 514, 647, 595, 714, …
#> $ median_income      <dbl> 8.3252, 8.3014, 7.2574, 5.6431, 3.8462, 4.0368, 3.6…
#> $ ocean_proximity    <chr> "NEAR BAY", "NEAR BAY", "NEAR BAY", "NEAR BAY", "NE…
#> $ median_house_value <int> 452600, 358500, 352100, 341300, 342200, 269700, 299…

Dari fungsi glimps() yang telah kita jalankan diatas, kita mendapatkan informasi bahwa terdapat 20,640 rows dan 10 columns.

`longitude`             : Ukuran seberapa jauh letak sebuah rumah di arah barat; nilai yang lebih tinggi terletak lebih jauh ke barat
`latitude`              : Ukuran seberapa jauh letak rumah di utara; nilai yang lebih tinggi terletak lebih jauh ke utara
`housing_median_age`    : Usia rata-rata sebuah rumah dalam satu blok; angka yang lebih rendah adalah bangunan yang lebih baru
`total_rooms`           : Jumlah total kamar dalam satu blok
`total_bedrooms`        : Jumlah total kamar tidur dalam satu blok
`population`            : Jumlah total orang yang tinggal dalam satu blok
`households`            : Jumlah total rumah tangga, sekelompok orang yang tinggal dalam satu unit rumah, untuk satu blok
`median_income`         : Pendapatan median rumah tangga dalam satu blok rumah (diukur dalam puluhan ribu Dolar AS)
`ocean_proximity`       : Nilai rumah median untuk rumah tangga dalam satu blok (diukur dalam Dolar AS)
`median_house_value`    : Lokasi rumah w.r.t samudera/laut

CLEANSING DATA & PRE-PROCESSING

Pada proses ini kita akan mengecheck apakah data kita terdapat missing value atau tidak, kita check dengan menggunakan fungsi anyNA() atau colSums(is.na()) atau bisa menggunakan keduanya.

colSums(is.na(housing))
#>          longitude           latitude housing_median_age        total_rooms 
#>                  0                  0                  0                  0 
#>     total_bedrooms         population         households      median_income 
#>                207                  0                  0                  0 
#>    ocean_proximity median_house_value 
#>                  0                  0
housingNA <- housing %>% 
  select(-ocean_proximity)

Kita menemukan missing value pada columns Price dan Rating. Pada columns Price terdapat 20 missing value, sedangkan columns Rating terdapat 133 missing value. Maka dari itu kita remove untuk columns yang terdapat missing value dengan menggunakan fungsi na.omit().

housingNA <- na.omit(housingNA)

anyNA(housingNA)
#> [1] FALSE

EXPLORATORY DATA ANALYSIS (EDA)

Analisis persebaran data meggunakan str() dan summary()

str(housingNA)
#> 'data.frame':    20433 obs. of  9 variables:
#>  $ longitude         : num  -122 -122 -122 -122 -122 ...
#>  $ latitude          : num  37.9 37.9 37.9 37.9 37.9 ...
#>  $ housing_median_age: int  41 21 52 52 52 52 52 52 42 52 ...
#>  $ total_rooms       : int  880 7099 1467 1274 1627 919 2535 3104 2555 3549 ...
#>  $ total_bedrooms    : int  129 1106 190 235 280 213 489 687 665 707 ...
#>  $ population        : int  322 2401 496 558 565 413 1094 1157 1206 1551 ...
#>  $ households        : int  126 1138 177 219 259 193 514 647 595 714 ...
#>  $ median_income     : num  8.33 8.3 7.26 5.64 3.85 ...
#>  $ median_house_value: int  452600 358500 352100 341300 342200 269700 299200 241400 226700 261100 ...
#>  - attr(*, "na.action")= 'omit' Named int [1:207] 291 342 539 564 697 739 1098 1351 1457 1494 ...
#>   ..- attr(*, "names")= chr [1:207] "291" "342" "539" "564" ...
summary(housingNA)
#>    longitude         latitude     housing_median_age  total_rooms   
#>  Min.   :-124.3   Min.   :32.54   Min.   : 1.00      Min.   :    2  
#>  1st Qu.:-121.8   1st Qu.:33.93   1st Qu.:18.00      1st Qu.: 1450  
#>  Median :-118.5   Median :34.26   Median :29.00      Median : 2127  
#>  Mean   :-119.6   Mean   :35.63   Mean   :28.63      Mean   : 2636  
#>  3rd Qu.:-118.0   3rd Qu.:37.72   3rd Qu.:37.00      3rd Qu.: 3143  
#>  Max.   :-114.3   Max.   :41.95   Max.   :52.00      Max.   :39320  
#>  total_bedrooms     population      households     median_income    
#>  Min.   :   1.0   Min.   :    3   Min.   :   1.0   Min.   : 0.4999  
#>  1st Qu.: 296.0   1st Qu.:  787   1st Qu.: 280.0   1st Qu.: 2.5637  
#>  Median : 435.0   Median : 1166   Median : 409.0   Median : 3.5365  
#>  Mean   : 537.9   Mean   : 1425   Mean   : 499.4   Mean   : 3.8712  
#>  3rd Qu.: 647.0   3rd Qu.: 1722   3rd Qu.: 604.0   3rd Qu.: 4.7440  
#>  Max.   :6445.0   Max.   :35682   Max.   :6082.0   Max.   :15.0001  
#>  median_house_value
#>  Min.   : 14999    
#>  1st Qu.:119500    
#>  Median :179700    
#>  Mean   :206864    
#>  3rd Qu.:264700    
#>  Max.   :500001

Selanjutnya kita akan menemukan hubungan antar variabel

ggcorr(housingNA, label = TRUE, label_size = 2.9, hjust = 1, layout.exp = 2)

MODELING & INTERPRETATION

Sebelum melakukan proses modeling lebih jauh, kita akan akan membagi data menjadi 2 yaitu data train dan data test.

set.seed(100)

# index sampling
index <- sample(nrow(housingNA), nrow(housingNA)*0.75)

# splitting
housing_train <- housingNA[index,]
housing_test <- housingNA[-index,]

Model Fitting

Membuat model untuk keseluruhan prediktor dengan menggunakan fungsi lm().

housing_lm1 <- lm(formula = population ~.,
               data = housing_train)

summary(housing_lm1)
#> 
#> Call:
#> lm(formula = population ~ ., data = housing_train)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -4144.7  -183.1   -38.2   132.9 21944.6 
#> 
#> Coefficients:
#>                          Estimate     Std. Error t value             Pr(>|t|)
#> (Intercept)        -3762.51335423   499.85585903  -7.527   0.0000000000000547
#> longitude            -62.40858854     5.72052485 -10.910 < 0.0000000000000002
#> latitude             -94.06912686     5.40516285 -17.404 < 0.0000000000000002
#> housing_median_age     0.50452136     0.32572202   1.549                0.121
#> total_rooms            0.18316836     0.00564698  32.436 < 0.0000000000000002
#> total_bedrooms        -1.46708913     0.05038131 -29.120 < 0.0000000000000002
#> households             3.34710807     0.04915800  68.089 < 0.0000000000000002
#> median_income          3.94039060     3.22578783   1.222                0.222
#> median_house_value    -0.00152173     0.00005021 -30.305 < 0.0000000000000002
#>                       
#> (Intercept)        ***
#> longitude          ***
#> latitude           ***
#> housing_median_age    
#> total_rooms        ***
#> total_bedrooms     ***
#> households         ***
#> median_income         
#> median_house_value ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 444.8 on 15315 degrees of freedom
#> Multiple R-squared:  0.8506, Adjusted R-squared:  0.8505 
#> F-statistic: 1.09e+04 on 8 and 15315 DF,  p-value: < 0.00000000000000022

Hasil interpretasi dari summary(model_lm):

Significant Predictor : - Predictor yang berperan signifikan yaitu : longitude, latitude, housing_median_age, total_rooms, total_bedrooms, households, ocean_proximityINLAND, ocean_proximityNEAR BAY, ocean_proximityNEAR OCEAN, median_house_value. - Untuk variabel yang tidak signifikan dapat kita buang ketika membuat model. Interpretasi Coefficient : - Ketika longitude naik 1 nilai maka population turun sebesar 70.17730852 dengan catatan variabel lain bernilai tetap. - Ketika latitude naik 1 nilai maka population turun sebesar 91.13704249 dengan catatan variabel lain bernilai tetap. - Ketika total_rooms naik 1 nilai maka population naik sebesar 0.18509604 dengan catatan variabel lain bernilai tetap. - Ketika total_bedrooms naik 1 nilai maka population turun sebesar 1.39794410 dengan catatan variabel lain bernilai tetap. - Ketika households naik 1 nilai maka population naik sebesar 3.24173544 dengan catatan variabel lain bernilai tetap. - Ketika median_house_value naik 1 nilai maka population turun sebesar 0.00151332 dengan catatan variabel lain bernilai tetap. Goodness of Fit (adj r-squared) :

summary(housing_lm1)$r.squared
#> [1] 0.8505642

Adjusted R-squared: 0.2665256 - Prediktor yang kita gunakan mampu menjelaskan sebesar 26.62% target variabel kita. - Model kita mampu menangkap informasi sebesar 26.62% dari seluruh prediktor yang digunakan untuk memprediksi variabel target.

Kita akan membuat model baru dari hasil model diatas yang mempunya signifikansi.

housing_lm2 <- lm(housing_median_age ~ longitude + latitude + population + total_rooms + total_bedrooms + households + median_house_value,
                data = housingNA)

summary(housing_lm2)
#> 
#> Call:
#> lm(formula = housing_median_age ~ longitude + latitude + population + 
#>     total_rooms + total_bedrooms + households + median_house_value, 
#>     data = housingNA)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -34.974  -8.307   0.341   7.782  57.860 
#> 
#> Coefficients:
#>                           Estimate      Std. Error t value            Pr(>|t|)
#> (Intercept)        -276.1813029842   10.8107198691 -25.547 <0.0000000000000002
#> longitude            -3.4429234965    0.1237569964 -27.820 <0.0000000000000002
#> latitude             -2.8982819927    0.1182913639 -24.501 <0.0000000000000002
#> population            0.0001510235    0.0001812806   0.833              0.4048
#> total_rooms          -0.0029046380    0.0001098744 -26.436 <0.0000000000000002
#> total_bedrooms        0.0013922287    0.0010661913   1.306              0.1916
#> households            0.0030849362    0.0012253319   2.518              0.0118
#> median_house_value    0.0000075901    0.0000008697   8.727 <0.0000000000000002
#>                       
#> (Intercept)        ***
#> longitude          ***
#> latitude           ***
#> population            
#> total_rooms        ***
#> total_bedrooms        
#> households         *  
#> median_house_value ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 11.29 on 20425 degrees of freedom
#> Multiple R-squared:  0.1959, Adjusted R-squared:  0.1957 
#> F-statistic: 711.1 on 7 and 20425 DF,  p-value: < 0.00000000000000022

Goodness of Fit (adj r-squared) :

summary(housing_lm2)$adj.r.squared
#> [1] 0.1956694

Adjusted R-squared: 0.1956694 - Prediktor yang kita gunakan mampu menjelaskan sebesar 19.56% target variabel kita. - Model kita mampu menangkap informasi sebesar 19.56% dari seluruh prediktor yang digunakan untuk memprediksi variabel target.

Kita akan membuat model step-weise regression

housing_backward <- step(object = housing_lm1, 
                       direction = 'backward',
                       trace = 0)

summary(housing_backward)
#> 
#> Call:
#> lm(formula = population ~ longitude + latitude + total_rooms + 
#>     total_bedrooms + households + median_house_value, data = housing_train)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -4111.2  -182.2   -39.0   132.7 21952.3 
#> 
#> Coefficients:
#>                          Estimate     Std. Error t value             Pr(>|t|)
#> (Intercept)        -3854.03350959   490.55556499  -7.856  0.00000000000000421
#> longitude            -63.71823283     5.60120564 -11.376 < 0.0000000000000002
#> latitude             -95.23655192     5.31648457 -17.913 < 0.0000000000000002
#> total_rooms            0.18549709     0.00474102  39.126 < 0.0000000000000002
#> total_bedrooms        -1.48744946     0.04729940 -31.448 < 0.0000000000000002
#> households             3.35129386     0.04909586  68.260 < 0.0000000000000002
#> median_house_value    -0.00147840     0.00003776 -39.154 < 0.0000000000000002
#>                       
#> (Intercept)        ***
#> longitude          ***
#> latitude           ***
#> total_rooms        ***
#> total_bedrooms     ***
#> households         ***
#> median_house_value ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 444.8 on 15317 degrees of freedom
#> Multiple R-squared:  0.8505, Adjusted R-squared:  0.8505 
#> F-statistic: 1.453e+04 on 6 and 15317 DF,  p-value: < 0.00000000000000022

Goodness of Fit (adj r-squared) :

summary(housing_backward)$adj.r.squared
#> [1] 0.8504742

Adjusted R-squared: 0.2659507 - Prediktor yang kita gunakan mampu menjelaskan sebesar 26.59% target variabel kita. - Model kita mampu menangkap informasi sebesar 26.59% dari seluruh prediktor yang digunakan untuk memprediksi variabel target.

Model Comparisson

Tujuan dari model comparisson yaitu mendapatkan model terbaik untuk prediksi variabel target.

Coba kita lihat Goodness of Fit (R-Squared) dari ketiga model berdasarkan nilai R-squarednya.

summary(housing_lm1)$r.squared
#> [1] 0.8505642
summary(housing_lm2)$adj.r.squared
#> [1] 0.1956694
summary(housing_backward)$adj.r.squared
#> [1] 0.8504742

Dari di ke tiga model yang telah kita buat, kita telah mengetahui bahwa R-squared paling tinggi yaitu model_lm1 0.8537841.

Model Evaluation

Kita akan melakukan model evaluasi berdasarkan hasil prediksi model menggunakan nilai RMSE. Apakah nilai RMSE cukup kecil dibandingkan dengan range target variabel.

Setelah kita menemukan R-Squared yang paling tinggi, kita buat Model Evaluation nya.

#RMSE
RMSE(housing_lm1$fitted.values, housing_train$population)
#> [1] 444.6968
range(housing_train$population)
#> [1]     3 35682

Nilai RMSE dari model model_lm1 sudah cukup. Tetapi perlu dicek apakah terdapat nilai outlier atau tidak, jika ada nilai outlier maka mungkin menjadi RMSE ini nilainya besar dan sebaiknya modelnya kita perbaiki.

ASUMSI LINEAR REGRESSION

Asumsi model linear regression:

  1. Linearity
  2. Normality of Residuals
  3. Homoscedasticity of Residuals
  4. No Multicolinearity

Pada asumsi ini, kita hanya perlu cek asumsi kepada model yang akan digunakan. TidaK perlu semua model di cek asumsinya. Kita akan cek asumsi untuk model_lm1.

Linearity

Plot ini merupakan scatter plot dengan sumbu x adalah nilai fitted values (hasil prediksi variabel respon) dan sumbu y adalah nilai residual.

plot(housing_lm1, which = 1)

Nilai residual “bounce randomly” di sekitar nilai 0. Model backward dapat dikatakan punya hubungan linear antara prediktor dengan targetnya.

Normality of Residuals

Model linear regression diharapkan menghasilkan error yang berdistribusi normal. Dengan begitu, error lebih banyak berkumpul di sekitar angka nol.

#visualisasi histogram
hist(housing_lm1$residuals)

shapiro.test(housing_lm1$residuals[0:5000])
#> 
#>  Shapiro-Wilk normality test
#> 
#> data:  housing_lm1$residuals[0:5000]
#> W = 0.51096, p-value < 0.00000000000000022

p-value < alpha (0.05) sehingga error tidak berdistrubi normal.

library(olsrr)

ols_plot_resid_hist(housing_lm1)

ols_plot_resid_qq(housing_lm1)

Homoscedasticity of Residuals

Diharapkan error yang dihasilkan oleh model menyebar secara acak atau dengan variasi konstan. Apabila divisualisasikan maka error tidak berpola.

plot(x = housing_lm1$fitted.values, y = housing_lm1$residuals)
abline(h = 0, col = "red")

bptest(housing_lm1)
#> 
#>  studentized Breusch-Pagan test
#> 
#> data:  housing_lm1
#> BP = 433.07, df = 8, p-value < 0.00000000000000022

p-value < alpha (0.05) error menyebar TIDAK konstan atau heteroscedasticity.

No Multicollinearity

Multicollinearity adalah kondisi adanya korelasi antar prediktor yang kuat. Hal ini tidak diinginkan karena menandakan prediktor redundan pada model, yang seharusnya dapat dipilih salah satu saja dari variable yang hubungannya amat kuat tersebut. Harapannya tidak terjadi multicollinearity.

vif(housing_lm1)
#>          longitude           latitude housing_median_age        total_rooms 
#>          10.165889          10.300575           1.306404          11.889046 
#>     total_bedrooms         households      median_income median_house_value 
#>          35.351083          27.753258           2.922782           2.595673

ada beberapa prediktor yang nilai vif > 10, artinya terjadi multicollinearity pada model.

SUMMARY

Variable yang kita gunakan untuk model yaitu housing_lm1 dengan uji coba semua variabel. Dengan R-squared tertinggi adalah 0.8447789 (84.47%). Accuracy model menggunakan RMSE mendapatkan nilai sebesar 442.6317. Dan untuk assumsi model hanya model Linearity yang memenuhi uji asumsi.