Tahap 1 : Pendahuluan

1.Business Problem

Tujuan dari project ini adalah untuk memprediksi harga rumah di Boston menggunakan metode regresi linear. Ada beberapa model yang digunakan yaitu pemodelan satu prediktor, pemodelan seluruh prediktor dan pemodelan stepwise. kemudian dilakukan komparasi untuk digunakan satu pemodelan untuk melakukan cross validation dll.

2.Dataset

  • Dataset boston adalah informasi yang dikumpulkan oleh US Census Service terkait dengan harga perumahan dan informasi penunjang di area kawasan perumahan. Ukuran dataset relatif kecil hanya terdiri dari 14 kolom dan 506 baris.

  • Dataset ini dapat diunduh melalui https://www.kaggle.com/datasets/vikrishnan/boston-house-prices

Tahap 2 : Pembacaan dan Pemahaman Data

Tahapan ini adalah melakukan pembacaan data dan memahami data. Pada tahapan ini ditunjukan cara membaca data kemudian di simpan di dalam objek dan ditampilkan secara tabular. Kemudian diketahui informasi tiap kolom data, lalu tiap baris data di dataset.

Pembacaan Data

Tahapan ini adalah melakukan pembacaan data menggunakan fungsi base pada R yaitu read.csv() karena file dataset dalam project ini menggunakan file berformat csv. Kemudian menggunakan fungsi head() untuk melihat cuplikan data. Selanjutnya data dimasukan ke dalam objek yang diberi nama boston.

boston <- read.csv('HousingData.csv')
head(boston, n=3)  # menggunakan 3 data teratas 
##      CRIM ZN INDUS CHAS   NOX    RM  AGE    DIS RAD TAX PTRATIO      B LSTAT
## 1 0.00632 18  2.31    0 0.538 6.575 65.2 4.0900   1 296    15.3 396.90  4.98
## 2 0.02731  0  7.07    0 0.469 6.421 78.9 4.9671   2 242    17.8 396.90  9.14
## 3 0.02729  0  7.07    0 0.469 7.185 61.1 4.9671   2 242    17.8 392.83  4.03
##   MEDV
## 1 24.0
## 2 21.6
## 3 34.7

Pemahaman Data

1. Informasi mengenai setiap kolom

Dalam setiap kolomya diketahui informasi,

  • CRIM = per capita crime rate by town
  • ZN = proportion of residential land zoned for lots over 25,000 sq.ft.
  • INDUS = proportion of non-retail business acres per town.
  • CHAS = Charles River dummy variable (1 if tract bounds river; 0 otherwise)
  • NOX = nitric oxides concentration (parts per 10 million)
  • RM = average number of rooms per dwelling
  • AGE = proportion of owner-occupied units built prior to 1940
  • DIS = weighted distances to five Boston employment centres
  • RAD = index of accessibility to radial highways
  • TAX = full-value property-tax rate per $10,000
  • PTRATIO = pupil-teacher ratio by town
  • B = - 1000(Bk - 0.63)^2 where Bk is the proportion of blacks by town
  • LSTAT = % lower status of the population
  • MEDV = Median value of owner-occupied homes in $1000’s

2. Pengamatan yang diperoleh dari dataset

Dalam setiap barisnya, dataset ini terdiri dari informasi harga rumah di wilayah Boston MEDV, dilengkapi dengan data pelengkap yaitu data sekitar kawasan perumahan, data tarif pajak, data pendidikan dan ras yang tinggal di kawasan perumahan tersebut.

Tahap 3 : Pembersihan Data

1. Cek Kesesuaian Tipe data

Tahap ini melakukan pemeriksaan tipe data. Apakah tipe dataset yang telah disimpan sesuai ? untuk itu lihat struktur data dengan menggunakan fungsi glimpse yang merupakan fungsi dari library dplyr, dimana fungsi ini memiliki kemampuan yang cepat dan tampilan yang lebih rapi dibandingkan fungsi base R.

library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
glimpse(boston)
## Rows: 506
## Columns: 14
## $ CRIM    <dbl> 0.00632, 0.02731, 0.02729, 0.03237, 0.06905, 0.02985, 0.08829,…
## $ ZN      <dbl> 18.0, 0.0, 0.0, 0.0, 0.0, 0.0, 12.5, 12.5, 12.5, 12.5, 12.5, 1…
## $ INDUS   <dbl> 2.31, 7.07, 7.07, 2.18, 2.18, 2.18, 7.87, 7.87, 7.87, 7.87, 7.…
## $ CHAS    <int> 0, 0, 0, 0, 0, 0, NA, 0, 0, NA, 0, 0, 0, 0, NA, 0, 0, 0, 0, 0,…
## $ NOX     <dbl> 0.538, 0.469, 0.469, 0.458, 0.458, 0.458, 0.524, 0.524, 0.524,…
## $ RM      <dbl> 6.575, 6.421, 7.185, 6.998, 7.147, 6.430, 6.012, 6.172, 5.631,…
## $ AGE     <dbl> 65.2, 78.9, 61.1, 45.8, 54.2, 58.7, 66.6, 96.1, 100.0, 85.9, 9…
## $ DIS     <dbl> 4.0900, 4.9671, 4.9671, 6.0622, 6.0622, 6.0622, 5.5605, 5.9505…
## $ RAD     <int> 1, 2, 2, 3, 3, 3, 5, 5, 5, 5, 5, 5, 5, 4, 4, 4, 4, 4, 4, 4, 4,…
## $ TAX     <int> 296, 242, 242, 222, 222, 222, 311, 311, 311, 311, 311, 311, 31…
## $ PTRATIO <dbl> 15.3, 17.8, 17.8, 18.7, 18.7, 18.7, 15.2, 15.2, 15.2, 15.2, 15…
## $ B       <dbl> 396.90, 396.90, 392.83, 394.63, 396.90, 394.12, 395.60, 396.90…
## $ LSTAT   <dbl> 4.98, 9.14, 4.03, 2.94, NA, 5.21, 12.43, 19.15, 29.93, 17.10, …
## $ MEDV    <dbl> 24.0, 21.6, 34.7, 33.4, 36.2, 28.7, 22.9, 27.1, 16.5, 18.9, 15…

insight :
1. seluruh data bertipe numerik
2. dataset terdiri dari 506 baris dan 14 kolom

Agar semakin representatif, maka kita perlu melihat masing-masing nilai unique(berbeda) dalam tiap kolomnya. Gunakan fungsi n_distinct untuk melihat keseluruhan data.

sapply(boston, n_distinct)
##    CRIM      ZN   INDUS    CHAS     NOX      RM     AGE     DIS     RAD     TAX 
##     485      27      77       3      81     446     349     412       9      66 
## PTRATIO       B   LSTAT    MEDV 
##      46     357     439     229

insight :
1. terdapat dua kolom yang memiliki unique sedikit yaitu kolom CHAS dan kolom RAD
2. mayoritas kolom umumnnya memiliki nilai unique puluhan hingga ratusan

Dari insight yang diperoleh diatas, tipe data kolom CHAS akan di ubah menjadi tipe data factor karena kolom tersebut jika dilihat dari informasi kolom yang ada merupakan data biner. Sementara untuk data kolom RAD tidak diubah karena untuk kepentingan pengoalahan lanjutan (ingin melihat korelasi). Ubah tipe data dapat dilakukan dengan menggunakan fungsi mutate() yang merupakan fungsi untuk memanipulasi data salah satunya ubah tipe data yang berasal dari library dplyr. Jangan lupa masukan ke dalam variable setelah ubah tipe data.

boston <- boston |> 
  mutate(
    
    CHAS = as.factor(CHAS)
    
  )

2. Mengatasi Missing Value

Tahapan ini adalah melakukan cek missing value atau data NA. Dimana data tersebut akan dihapus jika terdapat pada dataset dikarenakan data missing value mampu menjadikan data tidak sesuai/ tidak representatif. Untuk mendeteksi apakah suatu dataset terdapat missing value gunakan fungsi anyNA(), dimana apabila hasilnya TRUE menunjukan bahwa dalam suatu dataset terdapat missing value.

anyNA(boston)
## [1] TRUE

Hasil diatas menunjukan adanya missing value pada dataset boston.Untuk itu kita lakukan inspeksi lebih lanjut untuk mengetahui seberapa banyak data misssing value ini.

colSums(is.na(boston)) # digunakan untuk melihat sebaran missing value pada masing-masing kolom
##    CRIM      ZN   INDUS    CHAS     NOX      RM     AGE     DIS     RAD     TAX 
##      20      20      20      20       0       0      20       0       0       0 
## PTRATIO       B   LSTAT    MEDV 
##       0       0      20       0
dim(na.omit(boston)) #melihat dataset apabila missing value dihilangkan
## [1] 394  14

dari nilai diatas diketahui, apabila missing value dihapus maka akan menghapus lebih dari 20% data. (506-394)/506. Karena terlalu banyak data yang hilang jika dihapus maka, akan dilakukan imputasi atau memasukan nilai missing value menggunakan nilai mean. Nilai mean dipilih karena dapat mendeskripsikan keseluruhan nilai dalam dataset.

#Proses Imputasi mean
boston$CRIM <- ifelse(is.na(boston$CRIM),
                      ave(boston$CRIM, FUN = function(x) mean(x, na.rm = T)), boston$CRIM)

boston$ZN <- ifelse(is.na(boston$ZN),
                      ave(boston$ZN, FUN = function(x) mean(x, na.rm = T)), boston$ZN)

boston$INDUS <- ifelse(is.na(boston$INDUS),
                      ave(boston$INDUS, FUN = function(x) mean(x, na.rm = T)), boston$INDUS)

boston$AGE <- ifelse(is.na(boston$AGE),
                      ave(boston$AGE, FUN = function(x) mean(x, na.rm = T)), boston$AGE)

boston$LSTAT <- ifelse(is.na(boston$LSTAT),
                      ave(boston$LSTAT, FUN = function(x) mean(x, na.rm = T)), boston$LSTAT)

Lakukan kembali periksa data missing value

colSums(is.na(boston))
##    CRIM      ZN   INDUS    CHAS     NOX      RM     AGE     DIS     RAD     TAX 
##       0       0       0      20       0       0       0       0       0       0 
## PTRATIO       B   LSTAT    MEDV 
##       0       0       0       0

Dari proses imputasi yang sudah dilakukan, diketahui bahwa data factor tidak dapat dilakukan imputasi. Untuk itu missing value yang ada di kolom CHAS akan dihapuskan karena memiliki persentase yang kecil dibandingkan total data

boston <- na.omit(boston)

3. Menghapus Data duplikat

melakukan penghapusan data yang sama jika ada. data dihapus karena akan memperkecil variasi data. Menghapus data duplikat dapat menggunakan fungsi distinct() yang terdapat pada library dplyr

boston <-  boston |>
  distinct()

4. Cek Kembali Structure Data

memastikan kembali kesesuaian tipe data. Hal ini penting untuk menghindari code yang belum dimasukan ke dalam objek sehingga tidak terjadi pernyimpanan hasil pengoalahan.

glimpse(boston)
## Rows: 486
## Columns: 14
## $ CRIM    <dbl> 0.00632, 0.02731, 0.02729, 0.03237, 0.06905, 0.02985, 0.14455,…
## $ ZN      <dbl> 18.0, 0.0, 0.0, 0.0, 0.0, 0.0, 12.5, 12.5, 12.5, 12.5, 12.5, 0…
## $ INDUS   <dbl> 2.31, 7.07, 7.07, 2.18, 2.18, 2.18, 7.87, 7.87, 7.87, 7.87, 7.…
## $ CHAS    <fct> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ NOX     <dbl> 0.538, 0.469, 0.469, 0.458, 0.458, 0.458, 0.524, 0.524, 0.524,…
## $ RM      <dbl> 6.575, 6.421, 7.185, 6.998, 7.147, 6.430, 6.172, 5.631, 6.377,…
## $ AGE     <dbl> 65.2, 78.9, 61.1, 45.8, 54.2, 58.7, 96.1, 100.0, 94.3, 82.9, 3…
## $ DIS     <dbl> 4.0900, 4.9671, 4.9671, 6.0622, 6.0622, 6.0622, 5.9505, 6.0821…
## $ RAD     <int> 1, 2, 2, 3, 3, 3, 5, 5, 5, 5, 5, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,…
## $ TAX     <int> 296, 242, 242, 222, 222, 222, 311, 311, 311, 311, 311, 307, 30…
## $ PTRATIO <dbl> 15.3, 17.8, 17.8, 18.7, 18.7, 18.7, 15.2, 15.2, 15.2, 15.2, 15…
## $ B       <dbl> 396.90, 396.90, 392.83, 394.63, 396.90, 394.12, 396.90, 386.63…
## $ LSTAT   <dbl> 4.98000, 9.14000, 4.03000, 2.94000, 12.71543, 5.21000, 19.1500…
## $ MEDV    <dbl> 24.0, 21.6, 34.7, 33.4, 36.2, 28.7, 27.1, 16.5, 15.0, 18.9, 21…

Tahap 4 : Analisis Data Eksploratif

Dikarenakan metode yang digunakan adalah regresi. maka tahapan ini akan melakukan eksplorasi data target MEDV dengan seluruh prediktornya sehingga diharapkan dapat diketahui hubungan antara variable. Untuk menemukan hal tersebut dalam tahapan ini dilakukan cek data sebaran dan cek korelasi antar data.

1. Cek Sebaran Data antar Variable

Mengetahui sebaran data penting dilakukan dalam regresi linear, hal tersebut dikarenakan untuk melakukan regresi linear kita harus menggunakan data yang memiliki pola linear. sehingga mengatahui pola/sebaran penting untuk memastikan bahwa data yang digunakan dengan model yang diimplementasikan sesuai. Untuk cek sebaran data kita dapat lakukan secara visual karena mampu memberikan pengambilan insight secara cepat. Gunakan fungsi pairs() untuk membuat visualisasi hubungan antar variable. fungsi tersebut sangat cocok digunakan untuk dataset yang memiliki prediktor yang banyak karena mampu menggabungkan keseluruhan plot dalam satu plot.

pairs(~MEDV + PTRATIO + B + LSTAT + DIS + RM + CRIM, data = boston, main = 'Pola Antar-Variabel pada Dataset Boston')

insight =
1. dari hasil plot diatas diketahui bahwa, variable MEDV cenderung linear dengan variable RM
2. variable MEDV cenderung linear dengan variable LSTAT

2. Cek Korelasi antar Variable

Melihat korelasi antar variable dalam regresi linear adalah cara yang tepat untuk mendapatkan variable prediktor potensial. Hal ini dikarenakan korelasi menunjukan seberapa kuat hubungan linear antar dua variable. Korelasi memiliki hubungan dengan regresi yang menunjukan apakah masing-masing variable saling mempengaruhi. Untuk melakukan korelasi dapat menggunakan fungsi cor(), Namun karena kita mengharapkan interpretasi yang mudah memberikan insight secara cepat dan melakukan pengoalahan terhadap banyak variable sekaligus maka gunakan fungsi ggcorr() di library GGally untuk memberikan visualisasi korelasi antar banyak variable.

library(GGally)
## Loading required package: ggplot2
## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2
ggcorr(boston, label = T, hjust = 0.9)
## Warning in ggcorr(boston, label = T, hjust = 0.9): data in column(s) 'CHAS' are
## not numeric and were ignored

insight :
1. korelasi hubungan linear paling kuat diperoleh oleh Varibel RM sebesar 0.7 terhadap Variable target MEDV
2. korelasi yang paling tidak berhubungan diperoleh oleh Variable DIS sebesar 0.2 terhadap Variable target MEDV
3. Variable RM akan digunakan untuk pemodelan satu prediktor

Tahap 5 : Pemodelan Regresi dan Interpretasi

Dalam pengerjaan ini dilakukan penggunaan 3 buah model regresi liniear yaitu model_ols (model 1 prediktor), model_all (model seluruh prediktor) dan model_step (model menggunakan metode stepwise) menggunakan data boston. Kemudian ketiga model tersebut dibandingkan menurut R-Squared dan dipilih nilai tertinggi untuk melakukan tahapan selanjutnya

1. Pemodelan Regresi Linear 1 Prediktor

Pemodelan ini menggunakan prediktor RM yang merupakan prediktor dengan korelasi paling tinggi dengan target.

model_ols <- lm(MEDV ~ RM, data = boston)
summary(model_ols)
## 
## Call:
## lm(formula = MEDV ~ RM, data = boston)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -24.255  -2.487   0.108   3.117  39.867 
## 
## Coefficients:
##             Estimate Std. Error t value            Pr(>|t|)    
## (Intercept) -36.8574     2.6983  -13.66 <0.0000000000000002 ***
## RM            9.4547     0.4259   22.20 <0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.558 on 484 degrees of freedom
## Multiple R-squared:  0.5045, Adjusted R-squared:  0.5035 
## F-statistic: 492.8 on 1 and 484 DF,  p-value: < 0.00000000000000022

insight :
- Prediktor variable RM merupakan prediktor yang signifikan
- Dari model_ols data diperoleh nilai Intercept = -36.8574 dan koefisien RM = 9.4547 Sehingga diperoleh formula model regresinya adalah: MEDV = -36.8574 + 9.4547*RM
- variansi variabel target MEDV dijelaskan sebesar 50.35% oleh variable RM. Sisanya dijelaskan oleh variabel lain yang belum diketahui.

2. Pemodelan Regresi Linear Seluruh Prediktor

Pemodelan ini menggunakan seluruh prediktor yang tersedia. Tujuan dibuat model ini untuk mengetahui apakah model_all mampu menghasilkan pemodelan yang lebih efisien dibandingkan model_ols atau sebaliknya.

model_all <- lm(MEDV ~., data = boston)
summary(model_all)
## 
## Call:
## lm(formula = MEDV ~ ., data = boston)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -16.7626  -2.8534  -0.5878   1.7728  28.0120 
## 
## Coefficients:
##               Estimate Std. Error t value             Pr(>|t|)    
## (Intercept)  30.310808   5.314841   5.703      0.0000000207828 ***
## CRIM         -0.107352   0.033081  -3.245             0.001257 ** 
## ZN            0.038322   0.013991   2.739             0.006396 ** 
## INDUS        -0.033622   0.061754  -0.544             0.586394    
## CHAS1         2.969330   0.885392   3.354             0.000862 ***
## NOX         -16.004708   3.889637  -4.115      0.0000457476721 ***
## RM            4.487547   0.433634  10.349 < 0.0000000000000002 ***
## AGE          -0.013400   0.013246  -1.012             0.312212    
## DIS          -1.458987   0.202662  -7.199      0.0000000000024 ***
## RAD           0.261734   0.067651   3.869             0.000125 ***
## TAX          -0.010153   0.003813  -2.663             0.008012 ** 
## PTRATIO      -0.921395   0.136109  -6.770      0.0000000000385 ***
## B             0.010677   0.002818   3.789             0.000171 ***
## LSTAT        -0.443971   0.051303  -8.654 < 0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.818 on 472 degrees of freedom
## Multiple R-squared:  0.7392, Adjusted R-squared:  0.732 
## F-statistic: 102.9 on 13 and 472 DF,  p-value: < 0.00000000000000022

insight =
- Berdasarkan Significant Predictornya, mayoritas prediktor adalah prediktor yang signifikan, terlihat dari hasil nilai Pr(>|t|) banyak yang dibawah 0.05

  • Berdasarkan nilai slope :
    • Setiap kenaikan 1 satuan pada prediktor CRIM, maka nilai target MEDV variabel akan menurun sebanyak 0.107352 dengan catatan variabel lain tetap
    • Setiap kenaikan 1 satuan pada prediktor ZN, maka nilai target MEDV variabel akan naik sebanyak 0.038322 dengan catatan variabel lain tetap
    • Variable CHAS1 meningkatkan nilai MEDV sebesar 2.969330 poin dibandingkan jika bukan CHAS1
    • Setiap kenaikan 1 satuan pada prediktor NOX, maka nilai target MEDV variabel akan menurun sebanyak 16.004708 dengan catatan variabel lain tetap
    • Setiap kenaikan 1 satuan pada prediktor RM, maka nilai target MEDV variabel akan naik sebanyak 4.487547 dengan catatan variabel lain tetap
    • Setiap kenaikan 1 satuan pada prediktor DIS, maka nilai target MEDV variabel akan menurun sebanyak 1.458987 dengan catatan variabel lain tetap
    • Setiap kenaikan 1 satuan pada prediktor RAD, maka nilai target MEDV variabel akan naik sebanyak 0.261734 dengan catatan variabel lain tetap
    • Setiap kenaikan 1 satuan pada prediktor TAX, maka nilai target MEDV variabel akan turun sebanyak 0.010153 dengan catatan variabel lain tetap
    • Setiap kenaikan 1 satuan pada prediktor PTRATIO, maka nilai target MEDV variabel akan turun sebanyak 0.921395 dengan catatan variabel lain tetap
    • Setiap kenaikan 1 satuan pada prediktor B, maka nilai target MEDV variabel akan naik sebanyak 0.010677 dengan catatan variabel lain tetap
    • Setiap kenaikan 1 satuan pada prediktor LSTAT, maka nilai target MEDV variabel akan turun sebanyak 0.443971 dengan catatan variabel lain tetap
  • variansi variable target MEDV dapat dijelaskan oleh variansi CRIM sampai LSTAT sebesar 73.20% . sisanya dijelaskan oleh variable lain yang belum diketahui. model ini lebih efisien karena hasil yang diperoleh lebih baik dari pemodelan satu prediktor model_ols.

3. Pemodelan Regresi Linear Menggunakan Stepwise

Pemodelan ini menggunakan metode stepwise backward. Tujuan dibuat model ini untuk mengetahui apakah model_step mampu menghasilkan pemodelan yang lebih efisien dibandingkan model_all atau sebaliknya.

#stepwise regression: backward elimination
model_step <- step(object = model_all, 
                       direction = 'backward'
                       ) #tidak mengeluarkan output trace = FALSE
## Start:  AIC=1542.16
## MEDV ~ CRIM + ZN + INDUS + CHAS + NOX + RM + AGE + DIS + RAD + 
##     TAX + PTRATIO + B + LSTAT
## 
##           Df Sum of Sq   RSS    AIC
## - INDUS    1      6.88 10964 1540.5
## - AGE      1     23.76 10981 1541.2
## <none>                 10957 1542.2
## - TAX      1    164.61 11122 1547.4
## - ZN       1    174.15 11131 1547.8
## - CRIM     1    244.47 11202 1550.9
## - CHAS     1    261.10 11218 1551.6
## - B        1    333.26 11290 1554.7
## - RAD      1    347.48 11305 1555.3
## - NOX      1    393.04 11350 1557.3
## - PTRATIO  1   1063.84 12021 1585.2
## - DIS      1   1203.14 12160 1590.8
## - LSTAT    1   1738.51 12696 1611.7
## - RM       1   2486.17 13443 1639.5
## 
## Step:  AIC=1540.46
## MEDV ~ CRIM + ZN + CHAS + NOX + RM + AGE + DIS + RAD + TAX + 
##     PTRATIO + B + LSTAT
## 
##           Df Sum of Sq   RSS    AIC
## - AGE      1     23.81 10988 1539.5
## <none>                 10964 1540.5
## - ZN       1    180.78 11145 1546.4
## - TAX      1    238.05 11202 1548.9
## - CRIM     1    242.29 11206 1549.1
## - CHAS     1    255.08 11219 1549.6
## - B        1    338.53 11303 1553.2
## - RAD      1    402.62 11367 1556.0
## - NOX      1    437.71 11402 1557.5
## - PTRATIO  1   1107.58 12072 1585.2
## - DIS      1   1234.65 12199 1590.3
## - LSTAT    1   1747.80 12712 1610.3
## - RM       1   2542.30 13506 1639.8
## 
## Step:  AIC=1539.52
## MEDV ~ CRIM + ZN + CHAS + NOX + RM + DIS + RAD + TAX + PTRATIO + 
##     B + LSTAT
## 
##           Df Sum of Sq   RSS    AIC
## <none>                 10988 1539.5
## - ZN       1    199.33 11187 1546.2
## - TAX      1    244.64 11233 1548.2
## - CRIM     1    246.48 11234 1548.3
## - CHAS     1    249.32 11237 1548.4
## - B        1    327.12 11315 1551.8
## - RAD      1    420.78 11409 1555.8
## - NOX      1    529.82 11518 1560.4
## - PTRATIO  1   1139.76 12128 1585.5
## - DIS      1   1240.16 12228 1589.5
## - LSTAT    1   2044.71 13033 1620.5
## - RM       1   2534.96 13523 1638.4
summary(model_step)
## 
## Call:
## lm(formula = MEDV ~ CRIM + ZN + CHAS + NOX + RM + DIS + RAD + 
##     TAX + PTRATIO + B + LSTAT, data = boston)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -16.5944  -2.8705  -0.5188   1.8591  27.9794 
## 
## Coefficients:
##               Estimate Std. Error t value             Pr(>|t|)    
## (Intercept)  30.744417   5.297174   5.804     0.00000001186813 ***
## CRIM         -0.107709   0.033031  -3.261             0.001191 ** 
## ZN            0.040587   0.013841   2.932             0.003527 ** 
## CHAS1         2.883220   0.879164   3.279             0.001116 ** 
## NOX         -17.483763   3.657140  -4.781     0.00000233459563 ***
## RM            4.435654   0.424171  10.457 < 0.0000000000000002 ***
## DIS          -1.374601   0.187935  -7.314     0.00000000000111 ***
## RAD           0.276743   0.064956   4.260     0.00002461776653 ***
## TAX          -0.011182   0.003442  -3.249             0.001242 ** 
## PTRATIO      -0.942039   0.134348  -7.012     0.00000000000814 ***
## B             0.010537   0.002805   3.757             0.000194 ***
## LSTAT        -0.460114   0.048991  -9.392 < 0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.815 on 474 degrees of freedom
## Multiple R-squared:  0.7384, Adjusted R-squared:  0.7324 
## F-statistic: 121.7 on 11 and 474 DF,  p-value: < 0.00000000000000022

insight:
- nilai AIC terendah adalah 1539.52
- variansi variable target MEDV dapat dijelaskan menggunakan model stepwise backward sebesar 73.24% . sisanya dijelaskan oleh variable lain yang belum diketahui. model ini sedikit lebih efisien karena hasil yang diperoleh sedikit lebih baik dari pemodelan keseluruhan prediktor model_all

4. Kesimpulan

Untuk tahapan selajutnya akan menggunakan model_step karena memberikan nilai R-Squared yang paling tinggi. sebagaimana rangkuman hasil dibawah ini.

summary(model_ols)$r.squared       #menggunakan R-squared karena 1 variable prediktor
## [1] 0.5045159
summary(model_all)$adj.r.squared   #menggunakan Adjusted R-squared karena lebih dari 1 variable prediktor
## [1] 0.7319824
summary(model_step)$adj.r.squared  #menggunakan Adjusted R-squared karena lebih dari 1 variable prediktor
## [1] 0.7323658

Tahap 6 : Uji Silang

Uji silang atau cross validation adalah membagi dataset boston menjadi data test boston_test dan data train boston_train. Bertujuan untuk mengetahui seberapa baik model memprediksi data yang belum dilihat. Untuk membagi data menjadi dua gunakan fungsi sample() untuk melakukan sampling, gunakan set.seed untuk mengunci data. Proporsi perbandingan yang digunakan adalah 80:20 mengacu hal yang umum digunakan dalam uji silang.

1. Subsetting Index

# digunakan untuk mengunci nilai random
set.seed(100)

# index sampling
index <- sample(nrow(boston), nrow(boston)*0.8) #nilai 0.8 atau 80% tidak mutlak tergantung kebutuhan

# splitting
boston_train <- boston[index,] #digunakan untuk membuat model
boston_test <- boston[-index,] #digunakan untuk pengujian

2. Melakukan Uji Silang pada Model terpilih (model_step)

Model Train

model_all_train <- lm(MEDV ~., data = boston_train)
model_all_test <- lm(MEDV ~., data = boston_test)

#stepwise regression: backward elimination
model_step_train <- step(object = model_all_train, 
                       direction = 'backward'
                       ) #tidak mengeluarkan output trace = FALSE
## Start:  AIC=1212.43
## MEDV ~ CRIM + ZN + INDUS + CHAS + NOX + RM + AGE + DIS + RAD + 
##     TAX + PTRATIO + B + LSTAT
## 
##           Df Sum of Sq     RSS    AIC
## - INDUS    1      3.85  8218.3 1210.6
## - AGE      1      8.56  8223.0 1210.8
## <none>                  8214.5 1212.4
## - TAX      1    140.89  8355.4 1217.0
## - ZN       1    194.86  8409.3 1219.5
## - CRIM     1    212.72  8427.2 1220.3
## - RAD      1    256.41  8470.9 1222.3
## - B        1    319.61  8534.1 1225.2
## - NOX      1    322.20  8536.7 1225.3
## - CHAS     1    423.99  8638.5 1230.0
## - PTRATIO  1    666.38  8880.8 1240.7
## - DIS      1    803.12  9017.6 1246.6
## - LSTAT    1   1185.11  9399.6 1262.7
## - RM       1   2017.49 10232.0 1295.6
## 
## Step:  AIC=1210.61
## MEDV ~ CRIM + ZN + CHAS + NOX + RM + AGE + DIS + RAD + TAX + 
##     PTRATIO + B + LSTAT
## 
##           Df Sum of Sq     RSS    AIC
## - AGE      1      8.53  8226.8 1209.0
## <none>                  8218.3 1210.6
## - TAX      1    197.04  8415.4 1217.8
## - ZN       1    198.98  8417.3 1217.9
## - CRIM     1    211.20  8429.5 1218.5
## - RAD      1    296.07  8514.4 1222.3
## - B        1    322.26  8540.6 1223.5
## - NOX      1    356.36  8574.7 1225.1
## - CHAS     1    420.43  8638.7 1228.0
## - PTRATIO  1    705.43  8923.7 1240.6
## - DIS      1    827.96  9046.3 1245.8
## - LSTAT    1   1200.93  9419.3 1261.5
## - RM       1   2042.18 10260.5 1294.7
## 
## Step:  AIC=1209.01
## MEDV ~ CRIM + ZN + CHAS + NOX + RM + DIS + RAD + TAX + PTRATIO + 
##     B + LSTAT
## 
##           Df Sum of Sq     RSS    AIC
## <none>                  8226.8 1209.0
## - TAX      1    200.31  8427.2 1216.3
## - ZN       1    211.44  8438.3 1216.9
## - CRIM     1    214.27  8441.1 1217.0
## - RAD      1    306.51  8533.4 1221.2
## - B        1    316.57  8543.4 1221.7
## - CHAS     1    418.51  8645.4 1226.3
## - NOX      1    424.69  8651.5 1226.5
## - PTRATIO  1    735.75  8962.6 1240.2
## - DIS      1    854.67  9081.5 1245.4
## - LSTAT    1   1322.02  9548.9 1264.8
## - RM       1   2045.10 10271.9 1293.2
summary(model_step_train)
## 
## Call:
## lm(formula = MEDV ~ CRIM + ZN + CHAS + NOX + RM + DIS + RAD + 
##     TAX + PTRATIO + B + LSTAT, data = boston_train)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -11.0614  -2.7806  -0.5588   1.7029  26.2418 
## 
## Coefficients:
##               Estimate Std. Error t value             Pr(>|t|)    
## (Intercept)  26.688349   5.933924   4.498    0.000009170482794 ***
## CRIM         -0.112440   0.035930  -3.129             0.001888 ** 
## ZN            0.047524   0.015288   3.109             0.002022 ** 
## CHAS1         4.519986   1.033496   4.373    0.000015852753374 ***
## NOX         -17.074205   3.875487  -4.406    0.000013770317823 ***
## RM            4.675047   0.483562   9.668 < 0.0000000000000002 ***
## DIS          -1.293183   0.206910  -6.250    0.000000001113947 ***
## RAD           0.258100   0.068958   3.743             0.000210 ***
## TAX          -0.010856   0.003588  -3.026             0.002651 ** 
## PTRATIO      -0.873920   0.150706  -5.799    0.000000014164362 ***
## B             0.011109   0.002921   3.804             0.000166 ***
## LSTAT        -0.424466   0.054607  -7.773    0.000000000000074 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.678 on 376 degrees of freedom
## Multiple R-squared:  0.747,  Adjusted R-squared:  0.7396 
## F-statistic: 100.9 on 11 and 376 DF,  p-value: < 0.00000000000000022

Model Test

#stepwise regression: backward elimination
model_step_test <- step(object = model_all_test, 
                       direction = 'backward'
                       ) #tidak mengeluarkan output trace = FALSE
## Start:  AIC=343.66
## MEDV ~ CRIM + ZN + INDUS + CHAS + NOX + RM + AGE + DIS + RAD + 
##     TAX + PTRATIO + B + LSTAT
## 
##           Df Sum of Sq    RSS    AIC
## - ZN       1      0.66 2455.8 341.68
## - CHAS     1      2.14 2457.3 341.74
## - AGE      1     14.09 2469.2 342.22
## - INDUS    1     14.75 2469.9 342.24
## - TAX      1     19.45 2474.6 342.43
## - B        1     32.53 2487.7 342.95
## - CRIM     1     39.88 2495.0 343.23
## <none>                 2455.2 343.66
## - NOX      1     64.46 2519.6 344.20
## - RAD      1     82.37 2537.5 344.89
## - LSTAT    1    296.70 2751.9 352.84
## - DIS      1    361.56 2816.7 355.12
## - PTRATIO  1    415.31 2870.5 356.97
## - RM       1    468.37 2923.5 358.77
## 
## Step:  AIC=341.68
## MEDV ~ CRIM + INDUS + CHAS + NOX + RM + AGE + DIS + RAD + TAX + 
##     PTRATIO + B + LSTAT
## 
##           Df Sum of Sq    RSS    AIC
## - CHAS     1      1.96 2457.8 339.76
## - AGE      1     14.83 2470.6 340.27
## - INDUS    1     16.57 2472.4 340.34
## - TAX      1     19.40 2475.2 340.45
## - B        1     32.47 2488.3 340.97
## - CRIM     1     39.41 2495.2 341.24
## <none>                 2455.8 341.68
## - NOX      1     65.98 2521.8 342.28
## - RAD      1     83.07 2538.9 342.94
## - LSTAT    1    300.28 2756.1 350.99
## - DIS      1    388.75 2844.6 354.08
## - RM       1    467.73 2923.5 356.77
## - PTRATIO  1    529.39 2985.2 358.81
## 
## Step:  AIC=339.76
## MEDV ~ CRIM + INDUS + NOX + RM + AGE + DIS + RAD + TAX + PTRATIO + 
##     B + LSTAT
## 
##           Df Sum of Sq    RSS    AIC
## - AGE      1     15.02 2472.8 338.36
## - TAX      1     17.54 2475.3 338.46
## - INDUS    1     18.24 2476.0 338.48
## - B        1     30.93 2488.7 338.99
## - CRIM     1     38.60 2496.4 339.29
## <none>                 2457.8 339.76
## - NOX      1     66.80 2524.6 340.39
## - RAD      1     81.47 2539.2 340.96
## - LSTAT    1    327.98 2785.8 350.04
## - DIS      1    388.35 2846.1 352.14
## - RM       1    481.52 2939.3 355.29
## - PTRATIO  1    527.79 2985.6 356.82
## 
## Step:  AIC=338.36
## MEDV ~ CRIM + INDUS + NOX + RM + DIS + RAD + TAX + PTRATIO + 
##     B + LSTAT
## 
##           Df Sum of Sq    RSS    AIC
## - TAX      1     17.65 2490.4 337.05
## - INDUS    1     19.75 2492.6 337.14
## - B        1     28.82 2501.6 337.49
## - CRIM     1     38.46 2511.2 337.87
## <none>                 2472.8 338.36
## - NOX      1     81.45 2554.2 339.53
## - RAD      1     85.34 2558.1 339.68
## - DIS      1    381.93 2854.7 350.43
## - RM       1    480.35 2953.1 353.75
## - LSTAT    1    516.81 2989.6 354.96
## - PTRATIO  1    527.02 2999.8 355.29
## 
## Step:  AIC=337.05
## MEDV ~ CRIM + INDUS + NOX + RM + DIS + RAD + PTRATIO + B + LSTAT
## 
##           Df Sum of Sq    RSS    AIC
## - B        1     31.87 2522.3 336.30
## - CRIM     1     38.69 2529.1 336.57
## - INDUS    1     46.02 2536.5 336.85
## <none>                 2490.4 337.05
## - NOX      1     88.17 2578.6 338.46
## - RAD      1     98.98 2589.4 338.87
## - DIS      1    418.90 2909.3 350.29
## - RM       1    485.78 2976.2 352.52
## - LSTAT    1    531.63 3022.1 354.02
## - PTRATIO  1    550.37 3040.8 354.62
## 
## Step:  AIC=336.3
## MEDV ~ CRIM + INDUS + NOX + RM + DIS + RAD + PTRATIO + LSTAT
## 
##           Df Sum of Sq    RSS    AIC
## <none>                 2522.3 336.30
## - INDUS    1     55.68 2578.0 336.44
## - CRIM     1     80.72 2603.0 337.39
## - RAD      1     90.25 2612.6 337.75
## - NOX      1     95.98 2618.3 337.96
## - DIS      1    428.69 2951.0 349.68
## - RM       1    473.22 2995.5 351.15
## - LSTAT    1    536.22 3058.5 353.19
## - PTRATIO  1    538.12 3060.4 353.25
summary(model_step_test)
## 
## Call:
## lm(formula = MEDV ~ CRIM + INDUS + NOX + RM + DIS + RAD + PTRATIO + 
##     LSTAT, data = boston_test)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -13.6358  -3.1889  -0.6427   2.1088  25.6825 
## 
## Coefficients:
##              Estimate Std. Error t value  Pr(>|t|)    
## (Intercept)  45.61423   11.51794   3.960  0.000151 ***
## CRIM         -0.13401    0.07941  -1.688  0.094974 .  
## INDUS        -0.21219    0.15139  -1.402  0.164496    
## NOX         -19.61079   10.65615  -1.840  0.069053 .  
## RM            3.86863    0.94674   4.086 0.0000958 ***
## DIS          -1.70026    0.43717  -3.889  0.000194 ***
## RAD           0.19550    0.10955   1.785  0.077744 .  
## PTRATIO      -1.20889    0.27743  -4.357 0.0000352 ***
## LSTAT        -0.50887    0.11699  -4.350 0.0000362 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.324 on 89 degrees of freedom
## Multiple R-squared:  0.7285, Adjusted R-squared:  0.7041 
## F-statistic: 29.85 on 8 and 89 DF,  p-value: < 0.00000000000000022

insight:
- dari perbandingan nilai tersebut diperoleh bahwa nilai R-Squared dari model_step_train sebesar 74.85% dan nilai R-Squared dari model_step_test sebesar 67.19%
- nilai AIC terendah yang diperoleh dari model_step_train sebesar 1217.48 dan nilai AIC terendah yang diperoleh dari model_step_test sebesar 327.53

Tahap 7 : Pemodelan Prediksi

Tahapan ini adalah membandingkan nilai MEDV dari train dataset dengan nilai MEDV dari test dataset (Hasil Prediksi)

library(ggplot2)
head(predict(object = model_step_train,
        newdata = boston_test, interval = 'confidence'), 5)
##         fit      lwr      upr
## 5  25.46751 23.93732 26.99770
## 8  12.99321 10.58237 15.40406
## 9  20.22645 18.37575 22.07716
## 10 21.76962 20.35242 23.18682
## 18 13.09865 11.90591 14.29140
head(predict(object = model_step_train,
        newdata = boston_test, interval = 'prediction'), 5)
##         fit       lwr      upr
## 5  25.46751 16.143558 34.79145
## 8  12.99321  3.484969 22.50146
## 9  20.22645 10.844576 29.60833
## 10 21.76962 12.463545 31.07569
## 18 13.09865  3.824108 22.37320
model_predict <- predict(model_step_train, boston_test)
test.boston <- cbind(boston_test, pred = model_predict)
ggplot() +
  geom_density(data=test.boston,
               mapping = aes(x = MEDV, color ='real'))+
  geom_density(data=test.boston,
               mapping = aes(x = pred, color ='predicted'))+
  scale_color_discrete(name = 'MEDV Distribution')

insight:
- bahwa antara model(real) dengan hasil prediksi tidak 100% cocok
- untuk perbandingan anatara data train dan data test akan bahas bagian berikutnya.

Tahap 8 : Pemodelan Komparasi

Dari tiga buah model yang telah dibuat selanjutnya, dilakukan perbandingan antar ketiga model tersebut. Perbandingan yang dilakukan meliputi Perbandingan informasi model R-Squared dan Metrics Evaluasi.

1. Komparasi Goodness of Fit

Goodness of Fit yang dinyatakan dalam nilai R-Squared merupakan nilai yang menjelaskan seberapa baik prediktor dapat menjelaskan keberagaman kelas target.Semakin mendekati 1, semakin baik prediktor dalam menjelaskan target.

summary(model_step_test)$adj.r.squared   #menggunakan Adjusted R-squared karena lebih dari 1 variable prediktor
## [1] 0.7040534
summary(model_step_train)$adj.r.squared  #menggunakan Adjusted R-squared karena lebih dari 1 variable prediktor
## [1] 0.7395868

insight :
- model_step_test memiliki adj.r.squared 0.704, sedikit lebih rendah dari model_step_train yang memiliki adj.r.squared 0.7395

2. Komparasi Metrics Evaluasi

Dalam hal ini adalah membandingkan nilai error melalui nilai MAE, MAPE, dan RMSE.

library(MLmetrics)
## Warning: package 'MLmetrics' was built under R version 4.3.1
## 
## Attaching package: 'MLmetrics'
## The following object is masked from 'package:base':
## 
##     Recall
MAE_train <- MAE(y_pred = predict(model_step, boston_train), 
    y_true = boston_train$MEDV)

MAE_test <- MAE(y_pred = predict(model_step, boston_test), 
    y_true = boston_test$MEDV)


MAPE_train <- MAPE(y_pred = predict(model_step, boston_train), 
    y_true = boston_train$MEDV)

MAPE_test <- MAPE(y_pred = predict(model_step, boston_test), 
    y_true = boston_test$MEDV)


RMSE_train <- RMSE(y_pred = predict(model_step, boston_train), 
    y_true = boston_train$MEDV)

RMSE_test <- RMSE(y_pred = predict(model_step, boston_test), 
    y_true = boston_test$MEDV)

data.frame('MAE_train' = MAE_train,'MAE_test' = MAE_test,'MAPE_train' = MAPE_train, 'MAPE_test' = MAPE_test, 'RMSE_train' = RMSE_train,'RMSE_test' = RMSE_test)
##   MAE_train MAE_test MAPE_train MAPE_test RMSE_train RMSE_test
## 1  3.238351 3.607121  0.1663126 0.1623774   4.626884  5.231005

Insight:
- nilai MAE_train yang diperoleh sebesar 3.309764 dan nilai MAE_test yang diperoleh sebesar 3.324385
- nilai MAPE_train yang diperoleh sebesar 0.1678104 dan nilai MAPE_test yang diperoleh sebesar 0.1564473
- nilai RMSE_train yang diperoleh sebesar 4.680828 dan nilai MAPE_test yang diperoleh sebesar 5.037413

Tahap 9 : Pemodelan Asumsi Regresi

1. Linearitas

Linearity artinya target variabel dengan prediktornya memiliki hubungan yang linear atau hubungannya bersifat garis lurus.

plot(model_step, which = 1)

insight:
- model_backward memenuhi asumsi linearity sebagaimana kondisi yang diharapkan Nilai residual “bounce randomly” di sekitar nilai 0.

2. Normalitas

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

1. Visualisasi histogram residual menggunakan fungsi hist()

hist(model_step$residuals)

2. Uji statistik dengan shapiro.test()

Shapiro-Wilk hypothesis test:

  • H0: error berdistribusi normal
  • H1: error TIDAK berdistribusi normal
  • Kondisi yang diharapkan: H0
shapiro.test(model_step$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  model_step$residuals
## W = 0.89403, p-value < 0.00000000000000022

insight:
- p-value < 0.05 sehingga tolak H0 –> terima H1 –> model_backward mempunyai error TIDAK berdistribusi normal –> asumsi Normality of Residuals tidak terpenuhi

3. Homoscedasticity

Diharapkan error yang dihasilkan oleh model menyebar secara acak atau dengan variasi konstan. Apabila divisualisasikan maka error tidak berpola. Kondisi ini disebut juga sebagai homoscedasticity

1. Visualisasi scatter plot

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

2. Uji statistik dengan bptest() dari package lmtest

Breusch-Pagan hypothesis test:

  • H0: error menyebar konstan atau homoscedasticity
  • H1: error menyebar TIDAK konstan atau heteroscedasticity
  • Kondisi yang diharapkan: H0
library(lmtest)
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
bptest(model_step)
## 
##  studentized Breusch-Pagan test
## 
## data:  model_step
## BP = 58.104, df = 11, p-value = 0.00000002082

insight:
p-value < 0.05 sehingga tolak H0 (menerima H1) –> model mempunyai error yang menyebar tidak konstan atau heteroscedasticity –> asumsinya tidak terpenuhi

4. Multicolinearity

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

Uji VIF (Variance Inflation Factor) dengan fungsi vif() dari package car:

  • nilai VIF > 10: terjadi multicollinearity pada model
  • nilai VIF < 10: tidak terjadi multicollinearity pada model
  • Kondisi yang diharapkan: VIF < 10
# vif dari model
library(car)
## Loading required package: carData
## 
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
## 
##     recode
vif(model_step)
##     CRIM       ZN     CHAS      NOX       RM      DIS      RAD      TAX 
## 1.701916 2.133868 1.054347 3.725858 1.840166 3.247334 6.569621 6.963005 
##  PTRATIO        B    LSTAT 
## 1.777505 1.361923 2.510595

insight:
- tidak ada multikolinearity pada model_step

Tahap 10: Kesimpulan

Model_step menjadi model yang paling baik karena memiliki nilai R-squared paling tinggi senilai 0.7323658 dibandingkan model lain, secara intrepretasi menunjukan bahwa model model_step dapat memprediksi 71.3%. ketika model step dilakukan cross validation, prediksi dan evaluasi diketahui bahwa model step train senilai 0.7040534 dan model step test senilai 0.7395868. Kemudian Hasil pengujian error menggunakan metode MAE menunjukan error kurang lebih sekitar 3.288 USD. Pada pengujian asumsi, model ini hanya berhasil melewati pengujian linearity dan Multicollinearity sedangkan pada gagal dalam pengujian Normality dan Heteroscedasticity. Kesimpulannya, model ini masih kurang baik jika hendak digunakan untuk memprediksi harga rumah terkait dataset ini.