Introduction

Pada projek kali ini saya akan mengimplementasikan model regresi linear untuk memprediksi berat ikan serta menganalisa hubungan antar variabel dalam kaitannya dengan prediksi berat ikan.

Dataset yang saya gunakan kali ini adalah data Fish Market yang diperoleh di kaggle

Data Preparation

About the Data

Berikut penjelasan singkat dari masing-masing variabel di dataset ini:

Variable	Description
Species	nama spesies ikan
Weight	berat ikan (dalam satuan Gram)
Length1	panjang vertikal (dalam satuan cm)
Length2	panjang diagonal (dalam satuan cm)
Length3	lebar diagonal (dalam cm)
Height	tinggi (dalam cm)
Weight	lebar diagonal (dalam cm)

Karena kita akan memprediksi berat ikan, maka variabel weight sebagai Target Variabel, sedangkan variabel lainnya sebagai variabel prediktor.

Import Library and Data

library(dplyr)
library(tidyr)
library(readr)
library(GGally)
library(car)
library(scales)
library(lmtest)
library(ggplot2)
library(plotly)
library(ggthemes)
library(MLmetrics)
library(performance)

Data Preparation

fish_df <- read.csv("data_input/Fish.csv")

str(fish_df)

## 'data.frame':    159 obs. of  7 variables:
##  $ ï..Species: chr  "Bream" "Bream" "Bream" "Bream" ...
##  $ Weight    : num  242 290 340 363 430 450 500 390 450 500 ...
##  $ Length1   : num  23.2 24 23.9 26.3 26.5 26.8 26.8 27.6 27.6 28.5 ...
##  $ Length2   : num  25.4 26.3 26.5 29 29 29.7 29.7 30 30 30.7 ...
##  $ Length3   : num  30 31.2 31.1 33.5 34 34.7 34.5 35 35.1 36.2 ...
##  $ Height    : num  11.5 12.5 12.4 12.7 12.4 ...
##  $ Width     : num  4.02 4.31 4.7 4.46 5.13 ...

Saya rename kolom pertama dari dataframe menjadi Species

fish_df <- fish_df %>% 
  rename(Species = ï..Species) %>% 
  mutate_if(is.character, as.factor)

str(fish_df)

## 'data.frame':    159 obs. of  7 variables:
##  $ Species: Factor w/ 7 levels "Bream","Parkki",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ Weight : num  242 290 340 363 430 450 500 390 450 500 ...
##  $ Length1: num  23.2 24 23.9 26.3 26.5 26.8 26.8 27.6 27.6 28.5 ...
##  $ Length2: num  25.4 26.3 26.5 29 29 29.7 29.7 30 30 30.7 ...
##  $ Length3: num  30 31.2 31.1 33.5 34 34.7 34.5 35 35.1 36.2 ...
##  $ Height : num  11.5 12.5 12.4 12.7 12.4 ...
##  $ Width  : num  4.02 4.31 4.7 4.46 5.13 ...

Kolom / Informasi pada dataset

names(fish_df)

## [1] "Species" "Weight"  "Length1" "Length2" "Length3" "Height"  "Width"

Cek Missing Value

colSums(is.na(fish_df))

## Species  Weight Length1 Length2 Length3  Height   Width 
##       0       0       0       0       0       0       0

Berikut 10 data awal dari dataframe

head(fish_df,10)

summary(fish_df)

##       Species       Weight          Length1         Length2     
##  Bream    :35   Min.   :   0.0   Min.   : 7.50   Min.   : 8.40  
##  Parkki   :11   1st Qu.: 120.0   1st Qu.:19.05   1st Qu.:21.00  
##  Perch    :56   Median : 273.0   Median :25.20   Median :27.30  
##  Pike     :17   Mean   : 398.3   Mean   :26.25   Mean   :28.42  
##  Roach    :20   3rd Qu.: 650.0   3rd Qu.:32.70   3rd Qu.:35.50  
##  Smelt    :14   Max.   :1650.0   Max.   :59.00   Max.   :63.40  
##  Whitefish: 6                                                   
##     Length3          Height           Width      
##  Min.   : 8.80   Min.   : 1.728   Min.   :1.048  
##  1st Qu.:23.15   1st Qu.: 5.945   1st Qu.:3.386  
##  Median :29.40   Median : 7.786   Median :4.248  
##  Mean   :31.23   Mean   : 8.971   Mean   :4.417  
##  3rd Qu.:39.65   3rd Qu.:12.366   3rd Qu.:5.585  
##  Max.   :68.00   Max.   :18.957   Max.   :8.142  
##

EDA

I. Persebaran data `Species`

ggplot(fish_df, aes(x = Species, y=Weight)) +
  geom_boxplot(aes(fill = Species)) +
  coord_flip() +
  geom_hline(yintercept = mean(fish_df$Weight), linetype = 3) +
  theme(legend.position = 'none') + 
  labs(
    title = "Weight Distribution",
    subtitle = "group by Species",
    x = "Species", y = "Weights (in Gram)"
  )

fish_df %>% group_by(Species) %>% 
  summarise(freq=n()) %>% 
  arrange(-freq) %>% 
  ggplot(aes(x=reorder(Species, -freq), 
             y=freq, fill = Species)) +
  geom_col() +
  theme(legend.position = 'none') + 
  labs(
    title = "Species Frequencies",
    x = "Species", y = "Frequencies"
  )

Species paling banyak di dataset ini secara berurutan, yakni Perch, Bream, dan Roach. Meskipun saya tidak familiar dengan nama species ikan di dataset ini, namun dengan menganalisa data Weight per species saya dapat gambaran bahwa kelompok ikan species Smelt, Roach, dan Perkki cenderung memiliki berat yang lebih ringan dibandingkan kelompok ikan dari species Bream, Whitefish,Perch, maupun Pike

II. Distribusi data numerik

boxplot(fish_df)

Rentang data Weight sangat jauh dibandingkan dengan data lainnya. Sebagai informasi data berat ikan di dataset ini dalam satuan gram, sedangkan dimensi ikan (Length, Height, dan Width) diukur dalam satuan centimeter

Weight <- fish_df$Weight
h <- hist(Weight)
xfit<-seq(min(Weight),max(Weight))
yfit<-dnorm(xfit,mean=mean(Weight),sd=sd(Weight))
yfit <- yfit*diff(h$mids[1:2])*length(Weight)
lines(xfit, yfit, col="blue", lwd=2)

plot(density(Weight))

fish_df %>% select(-Weight) %>% 
  pivot_longer(-Species, names_to = "variable", values_to = "values") %>% 
  group_by(variable) %>% 
  #filter(variable == "Height") %>% 
  ggplot(aes(x = values, fill = variable)) +
  geom_density(alpha = 0.4)

Jika dilihat sekilas dari plot, data variabel target (Weight) maupun data prediktor tidak berdistribusi normal

III. Species manakah yang paling berat? (secara rataan)

fish_df %>% group_by(Species) %>% 
  summarise(avg_Weight = mean(Weight),
            avg_Height = mean(Height),
            avg_Width = mean(Width)) %>% 
  arrange(-avg_Weight) %>% 
  ggplot(aes(x=avg_Weight,
             y=reorder(Species,avg_Weight),
             fill = Species)) +
  geom_col() +
  geom_label(aes(label = round(avg_Weight,2))) +
  geom_vline(xintercept = mean(fish_df$Weight), linetype = 3) +
  theme(legend.position = 'none') + 
  labs(
    title = "Average Weight Rank",
    subtitle = "per species",
    x = "Weights (in Gram)", y = "Species"
  )

## `summarise()` ungrouping output (override with `.groups` argument)

IV. Korelasi antar data

pairs(fish_df)

ggcorr(fish_df, label=T )

cor(fish_df %>% select_if(is.numeric))

##            Weight   Length1   Length2   Length3    Height     Width
## Weight  1.0000000 0.9157117 0.9186177 0.9230436 0.7243453 0.8865066
## Length1 0.9157117 1.0000000 0.9995173 0.9920310 0.6253779 0.8670497
## Length2 0.9186177 0.9995173 1.0000000 0.9941026 0.6404408 0.8735467
## Length3 0.9230436 0.9920310 0.9941026 1.0000000 0.7034089 0.8785202
## Height  0.7243453 0.6253779 0.6404408 0.7034089 1.0000000 0.7928810
## Width   0.8865066 0.8670497 0.8735467 0.8785202 0.7928810 1.0000000

Jika dilihat dari scatter plot per masing-masing variabel numerik mayoritas memiliki korelasi kuat dan positif.

V. Weight vs Height per `Species`

Berdasarkan nilai korelasi saya coba mencari tahu hubungan variabel target terhadap berdasarkan variabel prediktor kategoriknya (Species) untuk mengetahui apakah ada hubungan atau pola tertentu

plot_weightvsheight <- fish_df %>% select(c(Height,Weight,Species)) %>% 
ggplot(aes(x=Height, y=Weight)) +
  geom_point(aes(col = Species)) +
  geom_smooth(method = "lm")

ggplotly(plot_weightvsheight)

plot_weightvsheight2 <- plot_weightvsheight +
facet_wrap(~Species, ncol = 3)
plot_weightvsheight2

Variabel Weight berkorelasi positif dengan Height untuk setiap Species. Hal yang menarik dari pola pengamatan tersebut adalah kenaikan nilai Weight pada Species Perch dan Pike lebih besar dibandingkan Species lainnya untuk setiap penambahan satu cm Height

VI. Weight vs Width per `Species`

plot_weightvswidth <- fish_df %>% select(c(Width,Weight,Species)) %>% 
ggplot(aes(x=Width, y=Weight)) +
  geom_point(aes(col = Species)) +
  geom_smooth(method = "lm")
ggplotly(plot_weightvswidth)

plot_weightvwidth2 <- plot_weightvswidth +
facet_wrap(~Species, ncol = 3)
plot_weightvwidth2

Variabel Weight berkorelasi positif dengan Width untuk setiap Species. Hal yang menarik dari pola pengamatan tersebut, penambahan satu cm Width pada Species Smelt tidak berdampak signifikan terhadap kenaikan nilai Weight

Modelling

Train-Test Split (Holdout Validation)

Dari 159 data observasi saya bagi 80% sebagai data-train dan 20% sebagai data-test. Data-train digunakan untuk melakukan pemodelan dan data-test dianggap sebagai unseen data yang digunakan untuk menguji seberapa baik model yang dibuat.

Untuk pemodelan saya akan membuang kolom Species dari dataset karena berdasarkan analisis di bagian EDA tidak memiliki pengaruh yang signifikan terhadap variabel target.

fish_df_new <- fish_df %>% select(-Species)

set.seed(123)
row_data <- nrow(fish_df_new)
index <- sample(row_data, row_data*0.8)

data_train <- fish_df_new[ index, ]
data_test <- fish_df_new[ -index, ]

Model1 : Full predictor

model1 <- lm(Weight ~ .,data_train)
summary(model1)

## 
## Call:
## lm(formula = Weight ~ ., data = data_train)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -262.28  -73.51  -23.75   68.23  413.26 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -503.86      32.84 -15.343   <2e-16 ***
## Length1        64.49      47.97   1.344   0.1813    
## Length2       -23.55      49.75  -0.473   0.6368    
## Length3       -12.86      20.40  -0.630   0.5296    
## Height         24.06      10.17   2.365   0.0196 *  
## Width          13.76      24.72   0.557   0.5788    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 126.3 on 121 degrees of freedom
## Multiple R-squared:  0.8832, Adjusted R-squared:  0.8783 
## F-statistic: 182.9 on 5 and 121 DF,  p-value: < 2.2e-16

\[ y = -503.86 + 64.49 Length_1 - 23.55 Length_2 - 12.86 Length_3 + 24.06 Height + 13.76 Width \]

AIC(model1)

## [1] 1597.294

Model2: Stepwise - Backward

Selanjutnya akan dicoba membuat model dengan penentuan variabel prediktor secara otomatis menggunakan metode step-wise backward elimination.

Metode step-wise regression ini akan menghasilkan formula optimum berdasarkan nilai `AIC yang terendah, dimana semakin rendah nilai Akaike Information Criterion (AIC) tersebut, maka model yang dihasilkan semakin baik.

model2 <- step(model1, direction = "backward")

## Start:  AIC=1234.88
## Weight ~ Length1 + Length2 + Length3 + Height + Width
## 
##           Df Sum of Sq     RSS    AIC
## - Length2  1      3575 1933890 1233.1
## - Width    1      4943 1935258 1233.2
## - Length3  1      6340 1936656 1233.3
## - Length1  1     28836 1959151 1234.8
## <none>                 1930316 1234.9
## - Height   1     89252 2019567 1238.6
## 
## Step:  AIC=1233.12
## Weight ~ Length1 + Length3 + Height + Width
## 
##           Df Sum of Sq     RSS    AIC
## - Width    1      2783 1936674 1231.3
## - Length3  1     12860 1946750 1232.0
## <none>                 1933890 1233.1
## - Length1  1     68894 2002785 1235.6
## - Height   1     95062 2028952 1237.2
## 
## Step:  AIC=1231.3
## Weight ~ Length1 + Length3 + Height
## 
##           Df Sum of Sq     RSS    AIC
## <none>                 1936674 1231.3
## - Length3  1     31263 1967937 1231.3
## - Length1  1    151406 2088080 1238.9
## - Height   1    299135 2235808 1247.5

summary(model2)

## 
## Call:
## lm(formula = Weight ~ Length1 + Length3 + Height, data = data_train)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -267.50  -68.26  -24.57   66.06  402.36 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -500.765     31.375 -15.961  < 2e-16 ***
## Length1       49.868     16.082   3.101  0.00239 ** 
## Length3      -21.282     15.103  -1.409  0.16133    
## Height        27.898      6.401   4.359 2.73e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 125.5 on 123 degrees of freedom
## Multiple R-squared:  0.8828, Adjusted R-squared:  0.8799 
## F-statistic: 308.7 on 3 and 123 DF,  p-value: < 2.2e-16

AIC(model2)

## [1] 1593.712

\[ y = -503.77 + 49.87 Length_1 - 21.28 Length_2 + 27.9 Height \]

Model3: Stepwise - Forward

Jika metode backward mengurangi prediktor satu persatu (dari full predictor) sampai mendapatkan nilai AIC yang rendah, metode Forward akan mencari prediktor dengan mencoba menambahkan prediktor satu persatu hingga mendapatkan AIC yg optimal

model_intercept <- lm(Weight~1, data_train)
model3 <- step(model_intercept,
               scope = list(lower = model_intercept, upper = model1),
               direction = "forward")

## Start:  AIC=1497.54
## Weight ~ 1
## 
##           Df Sum of Sq      RSS    AIC
## + Length3  1  14270135  2249958 1246.3
## + Length2  1  14038129  2481964 1258.8
## + Length1  1  13948048  2572045 1263.3
## + Width    1  12653908  3866185 1315.1
## + Height   1   8328269  8191824 1410.5
## <none>                 16520093 1497.5
## 
## Step:  AIC=1246.34
## Weight ~ Length3
## 
##           Df Sum of Sq     RSS    AIC
## + Width    1    219626 2030332 1235.3
## + Height   1    161878 2088080 1238.9
## <none>                 2249958 1246.3
## + Length1  1     14150 2235808 1247.5
## + Length2  1      9299 2240659 1247.8
## 
## Step:  AIC=1235.3
## Weight ~ Length3 + Width
## 
##           Df Sum of Sq     RSS    AIC
## <none>                 2030332 1235.3
## + Height   1   27546.8 2002785 1235.6
## + Length2  1    3095.6 2027236 1237.1
## + Length1  1    1379.7 2028952 1237.2

summary(model3)

## 
## Call:
## lm(formula = Weight ~ Length3 + Width, data = data_train)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -283.69  -74.76  -28.23   84.91  403.11 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -518.105     32.628 -15.879  < 2e-16 ***
## Length3       21.614      2.041  10.589  < 2e-16 ***
## Width         53.984     14.740   3.662 0.000368 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 128 on 124 degrees of freedom
## Multiple R-squared:  0.8771, Adjusted R-squared:  0.8751 
## F-statistic: 442.5 on 2 and 124 DF,  p-value: < 2.2e-16

AIC(model3)

## [1] 1597.71

\[ y = -518.11 + 21.614 Length_3 + 53.98 Width \]

Model4: Stepwise - Both Direction

Sendangkan metode Stepwise Both Direction adalah kombinasi dari Forward dan Backward hingga diperoleh nilai AIC yang optimum

model4 <- step(model_intercept, 
               scope = list(lower = model_intercept, upper = model1),
               direction = "both")

## Start:  AIC=1497.54
## Weight ~ 1
## 
##           Df Sum of Sq      RSS    AIC
## + Length3  1  14270135  2249958 1246.3
## + Length2  1  14038129  2481964 1258.8
## + Length1  1  13948048  2572045 1263.3
## + Width    1  12653908  3866185 1315.1
## + Height   1   8328269  8191824 1410.5
## <none>                 16520093 1497.5
## 
## Step:  AIC=1246.34
## Weight ~ Length3
## 
##           Df Sum of Sq      RSS    AIC
## + Width    1    219626  2030332 1235.3
## + Height   1    161878  2088080 1238.9
## <none>                  2249958 1246.3
## + Length1  1     14150  2235808 1247.5
## + Length2  1      9299  2240659 1247.8
## - Length3  1  14270135 16520093 1497.5
## 
## Step:  AIC=1235.3
## Weight ~ Length3 + Width
## 
##           Df Sum of Sq     RSS    AIC
## <none>                 2030332 1235.3
## + Height   1     27547 2002785 1235.6
## + Length2  1      3096 2027236 1237.1
## + Length1  1      1380 2028952 1237.2
## - Width    1    219626 2249958 1246.3
## - Length3  1   1835853 3866185 1315.1

summary(model3)

## 
## Call:
## lm(formula = Weight ~ Length3 + Width, data = data_train)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -283.69  -74.76  -28.23   84.91  403.11 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -518.105     32.628 -15.879  < 2e-16 ***
## Length3       21.614      2.041  10.589  < 2e-16 ***
## Width         53.984     14.740   3.662 0.000368 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 128 on 124 degrees of freedom
## Multiple R-squared:  0.8771, Adjusted R-squared:  0.8751 
## F-statistic: 442.5 on 2 and 124 DF,  p-value: < 2.2e-16

AIC(model4)

## [1] 1597.71

\[ y = -518.11 - 21.61 Length_3 + 53.98 Width \]

Evaluation

Performance

Selanjutnya saya coba komparasi performa ankeempat model regresi linear yang telah dibuat. Kali ini saya menggunakan method compare_performance() dari package performance

compare_performance(model1, model2, model3, model4, rank=T, verbose=F, metrics = "common")

Dari komparasi empat model, model2 (stepwise- Backward) memiliki performa paling baik dengan nilai R2-adjusted palign tinggi (0.8799 ) dan skor AIC terendah (1593.71)

Untuk mengukur performance model2 saya menggunakan metric Root Mean Squared Error (RMSE). \[ RMSE = \sqrt{\frac{1}{n} \sum (\hat y - y)^2} \]

model2_predict_dtest <- predict(model2, data_test)

rmse_dtrain_m2 <- RMSE(y_pred = model2$fitted.values, 
                       y_true = data_train$Weight)
rmse_dtest_m2 <- RMSE(y_pred = model2_predict_dtest, 
                      y_true = data_test$Weight)

print(paste("RMSE data-train= ", round(rmse_dtrain_m2,2)))

## [1] "RMSE data-train=  123.49"

print(paste("RMSE data-test= ", round(rmse_dtest_m2,2)))

## [1] "RMSE data-test=  118.16"

Selisih skor RMSE antara data-train dan data test, dimana RMSE data-train sedikit lebih tinggi dibandingkan data-test. Sehingga bisa disimpulkan bahwa model1 tidak overfit maupun underfit (meskipun memiliki kecenderungan underfit)

Model Assumption

Sebagai model parametrik, least square harus memenuhi empat asumsi: linearity, normality, Homocesdasticity, dan Multicolinearity. Aspek asumsi tersebut perlu dipenuhi agar interpretasi model yang diperoleh tidak bersifat bias.

Asumsi ini hanya perlu dipenuhi jika tujuan membuat model regresi linear adalah menginginkan interpretasi atau melihat efek dari setiap prediktor terhadap nilai target variabel. Jika hanya ingin menggunakan regresi linear untuk melakukan prediksi, maka asumsi model tidak wajib dipenuhi.

Kita cek keempat asumsi tersebut terhadap model2

1. Linearity

df_res_vs_act <- data.frame(residual = model2$residuals, 
                            fitted = model2$fitted.values)

df_res_vs_act %>% ggplot(aes(fitted, residual)) + 
  geom_point() + geom_smooth() + 
  geom_hline(aes(yintercept = 0)) + 
  theme(panel.grid = element_blank(), panel.background = element_blank())

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

Berdasarkan residual plot vs prediksi ini , tidak terlihat pola yang signifikan sehingga dapat dikatakan bahwa model ini linear.

2. Normality

Residual atau error dari data perlu berdistribusi normal. Ketika kita membuat sebuah model regresi linear, persamaan yang akan kita dapat adalah sebagai berikut:

\[ y = \alpha + \beta\times x + \epsilon \\ \epsilon \sim Normal(0, \sigma) \]

dimana (\(\epsilon\)) adalah error atau residual yang tidak ditangkap model.

Model berasumsi bahwa error ini memiliki distribusi normal yang berpusat di angka 0. Untuk memeriksa distribusi dari nilai residual apakah berdistribusi normal (membentuk kurva lonceng) atau tidak, kita bisa menggunakan histogram.

hist(model2$residuals)

plot(density(model2$residuals))

Terlihat dari plot histogram, distribusi residual banyak disekitar angka nol.

Alternatif lainnya adalah melakukan uji statistik menggunakan uji Shapiro-Wilk dengan hipotesis:

\(H_0\) = Residu berdistribusi normal
\(H_1\) = Residu tidak berdistribusi normal

shapiro.test(model2$residuals)

## 
##  Shapiro-Wilk normality test
## 
## data:  model2$residuals
## W = 0.96106, p-value = 0.001047

Hasil uji Shapiro-Wilk, p-value < 0.05 sehingga kesimpulan secara statistik residual model tidak berdistribusi normal (tolak \(H_0\))

Menggunakan QQ-plot, jika data berdistribusi normal maka semua titik ada di dalam area biru dan mengikuti pola garis biru

qqPlot(model2$residuals)

## 143 145 
##  35  99

Dari hasil QQ-Plot, terdapat mayoritas residual berada di sekitar garis biru atau bisa dikatakan residual berdistribusi normal

3. Homocesdasticity

Homocesdasticity menunjukkan bahwa residual atau error bersifat konstan atau tidak membentuk pola tertentu. Jika error membentuk pola tertentu seperti garis linear atau mengerucut, maka kita sebut dengan Heterocesdasticity dan akan berpengaruh pada nilai standard error pada estimate/koefisien prediktor yan bias (terlalu sempit atau terlalu lebar).

Untuk pengujian secara statistik kali ini menggunakan uji Bresuch-Pagan dengan hipotesis:

\(H_0\): error bersifat konstan (Homocesdasticity)
\(H_1\): error bersifat tidak konstan (Heteroscesdasticity)

pengujian Bresuch-Pagan menggunakan method bptest dari package lmtestdocumentation

bptest(model2)

## 
##  studentized Breusch-Pagan test
## 
## data:  model2
## BP = 60.461, df = 3, p-value = 4.684e-13

Hasil uji Bresuch-Pagan, p-value < 0.05 sehingga kesimpulan secara statistik tolak \(H_0\) atau error tidak bersifat constant (heterocesdasticity)

4. Multicollinearity

Multikolinearitas terjadi ketika antar variabel prediktor yang digunakan pada model memiliki hubungan yang kuat. Ada atau tidak multikolinearitas dapat dilihat dari nilai Variance Inflation Factor (VIF). VIF merupakan ukuran yang menjelaskan seberapa besar variansi koefisien yang meningkat karena multikolinieritas

Ketika nilai VIF lebih dari 10 artinya terjadi multikolinearitas. Jika hal tersebut terjadi, maka bisa dipilih salah satu variabel yang dibuang dari model yang memiliki VIF > 10.

Asumsi: - Antara target variabel dan prediktor ada korelasi - Antara prediktor dengan prediktor lain tidak ada korelasi pengujian menggunakan fungsi vif() dari package car

vif(model2)

##    Length1    Length3     Height 
## 217.499860 259.609031   6.147768

Terlihat bahwa terdapat dua variabel yang memiliki nilai VIF diatas 10, yakni Length1 dan Length3. Dari pengujian diatas bisa disimpulkan terdapat multikolinaritas di model ini.

Model Improvement

Model Tuning

Berdasarkan hasil pengujian terdapat beberapa bagian yang belum memenuhi aspek yakni residu tidak berdistribusi normal, heterocesdasticity, dan Multicolinearity.

Jika dilihat kembali di bagian EDA, range nilai Weight sangat besar dibandingkan dengan variabel lain. Perbedaan range yang besar antar variabel akan berpengaruh terhadap model dan interpretasi yang dihasilkan¹. Oleh karena saya coba standarisasi satuan pengukuran yang dipakai di dataset (rescaling) agar range nilai antara variabel target dan variabel prediktor.

Berdasarkan hasil penelusuran terkait “Standard weight in fish”² (khususnya ikan air tawar), panjang dan berat ikan masing-masing diukur dalam satuan milimeter dan gram

note:

pada referensi yang sama, ternyata dijelaskan bahwa hubungan antara length (L) dan Weight(W) adalah Non-Linear, \[W = aL^b \].

Oleh karena itu saya akan coba tuning data dengan mentransformasikan data terkait dimensi ikan, yang awalnya diukur dalam centimeter menjadi milimeter serta drop variabel Length3 (salah satu variabel yang memiliki skor VIF > 10)

fish_df_tr <- fish_df %>%
  select(-c(Species, Length3)) %>% 
  mutate(Length1 = Length1*100,
         Length2 = Length2*100,
         Height = Height*100,
         Width = Width*100)

boxplot(fish_df_tr)

summary(fish_df_tr)

##      Weight          Length1        Length2         Height      
##  Min.   :   0.0   Min.   : 750   Min.   : 840   Min.   : 172.8  
##  1st Qu.: 120.0   1st Qu.:1905   1st Qu.:2100   1st Qu.: 594.5  
##  Median : 273.0   Median :2520   Median :2730   Median : 778.6  
##  Mean   : 398.3   Mean   :2625   Mean   :2842   Mean   : 897.1  
##  3rd Qu.: 650.0   3rd Qu.:3270   3rd Qu.:3550   3rd Qu.:1236.6  
##  Max.   :1650.0   Max.   :5900   Max.   :6340   Max.   :1895.7  
##      Width      
##  Min.   :104.8  
##  1st Qu.:338.6  
##  Median :424.9  
##  Mean   :441.7  
##  3rd Qu.:558.5  
##  Max.   :814.2

Dari hasil plot maupun summary, range data masing-masing variabel jadi lebih seimbang dibandingkan data awal sebelum di rescaling

ggcorr(fish_df_tr, label = T)

Buat model baru

# split the data
set.seed(123)
data_train2 <- fish_df_tr[index, ]
data_test2 <- fish_df_tr[-index, ]

model tune 1

# train the model - using Full Predictor
model_tune1 <- lm(Weight ~ ., data = data_train2)
summary(model_tune1)

## 
## Call:
## lm(formula = Weight ~ ., data = data_train2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -264.71  -74.16  -27.71   72.24  414.37 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -507.56081   32.23076 -15.748  < 2e-16 ***
## Length1        0.63255    0.47807   1.323 0.188269    
## Length2       -0.36212    0.45409  -0.797 0.426738    
## Height         0.18539    0.05155   3.597 0.000466 ***
## Width          0.23705    0.18981   1.249 0.214103    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 126 on 122 degrees of freedom
## Multiple R-squared:  0.8828, Adjusted R-squared:  0.8789 
## F-statistic: 229.7 on 4 and 122 DF,  p-value: < 2.2e-16

AIC(model_tune1)

## [1] 1595.711

model tune 2

# Create new model using stepwise backward
model_tune2 <- step(model_tune1, direction = "backward")

## Start:  AIC=1233.3
## Weight ~ Length1 + Length2 + Height + Width
## 
##           Df Sum of Sq     RSS    AIC
## - Length2  1     10095 1946750 1232.0
## - Width    1     24759 1961414 1232.9
## - Length1  1     27790 1964446 1233.1
## <none>                 1936656 1233.3
## - Height   1    205346 2142001 1244.1
## 
## Step:  AIC=1231.96
## Weight ~ Length1 + Height + Width
## 
##           Df Sum of Sq     RSS    AIC
## - Width    1     21186 1967937 1231.3
## <none>                 1946750 1232.0
## - Height   1    214068 2160818 1243.2
## - Length1  1   1914619 3861370 1316.9
## 
## Step:  AIC=1231.34
## Weight ~ Length1 + Height
## 
##           Df Sum of Sq     RSS    AIC
## <none>                 1967937 1231.3
## - Height   1    604109 2572045 1263.3
## - Length1  1   6223888 8191824 1410.5

summary(model_tune2)

## 
## Call:
## lm(formula = Weight ~ Length1 + Height, data = data_train2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -282.40  -71.00  -24.24   79.14  382.39 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -503.90320   31.41984  -16.04  < 2e-16 ***
## Length1        0.27291    0.01378   19.80  < 2e-16 ***
## Height         0.20128    0.03262    6.17 8.91e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 126 on 124 degrees of freedom
## Multiple R-squared:  0.8809, Adjusted R-squared:  0.879 
## F-statistic: 458.5 on 2 and 124 DF,  p-value: < 2.2e-16

AIC(model_tune2)

## [1] 1593.746

compare_performance(model_tune1, model_tune2, rank=T, verbose=F, metrics = "common")

Performance

modelt2_predict_dtest <- predict(model_tune2, data_test2)

rmse_dtraintune_m2 <- RMSE(y_pred = model_tune2$fitted.values, 
                       y_true = data_train2$Weight)

rmse_dtesttune_m2 <- RMSE(y_pred = modelt2_predict_dtest, 
                      y_true = data_test2$Weight)

print(paste("RMSE data-train (after rescaling)= ", round(rmse_dtraintune_m2,2)))

## [1] "RMSE data-train (after rescaling)=  124.48"

print(paste("RMSE data-test (after rescaling)= ", round(rmse_dtesttune_m2,2)))

## [1] "RMSE data-test (after rescaling)=  130.14"

skor RMSE di model tuning lebih besar dibandingkan model sebelum dituning

Evaluation

Kita cek keempat asumsi tersebut terhadap model_tune2

1. Linearity

df_res_vs_act <- data.frame(residual = model_tune2$residuals, 
                            fitted = model_tune2$fitted.values)

df_res_vs_act %>% ggplot(aes(fitted, residual)) + 
  geom_point() + geom_smooth() + 
  geom_hline(aes(yintercept = 0)) + 
  theme(panel.grid = element_blank(), panel.background = element_blank())

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

2. Normality

uji Shapiro-Wilk

\(H_0\) = Residu berdistribusi normal \(H_1\) = Residu tidak berdistribusi normal

shapiro.test(model_tune2$residuals)

## 
##  Shapiro-Wilk normality test
## 
## data:  model_tune2$residuals
## W = 0.96039, p-value = 0.0009212

Hasil uji Shapiro-Wilk, p-value < 0.05 sehingga kesimpulan secara statistik residual model tidak berdistribusi normal (tolak \(H_0\))

qqPlot(model_tune2$residuals)

## 143 144 
##  35  26

3. Homocesdasticity

uji Bresuch-Pagan \[ H_0 = error\ bersifat\ konstan\ (Homocesdasticity)\\ H_1 = error\ bersifat\ tidak\ konstan\ (Heteroscesdasticity) \]

bptest(model_tune2)

## 
##  studentized Breusch-Pagan test
## 
## data:  model_tune2
## BP = 56.115, df = 2, p-value = 6.528e-13

Hasil uji Bresuch-Pagan, p-value < 0.05 sehingga kesimpulan secara statistik tolak \(H_0\) atau error tidak bersifat constant (heterocesdasticity)

4. Multicolinearity

Cek nilai VIF

vif(model_tune2)

##  Length1   Height 
## 1.584595 1.584595

model hasil rescaling sudah tidak bersifat multicolinearity.

Conclusion

Berdasarkan modelling dan hasil analisis regresi linear yang telah dilakukan dapat disimpulkan sebagai berikut:

Dari empat model regresi linear yang telah dibuat, model dari proses stepwise backward (model2) yang memiliki performa yang lebih baik dibandingkan ketiga model lainnya.
Hasil tuning dengan rescale data masih belum memenuhi uji asumsi model linear.
Beberapa point uji asumsi tidak terpenuhi (residu tidak berdistribusi normal, heterosscesdacity, dan multikoliearity) sehingga interpretasi dari model bisa bersifat bias, namun secara umum model ini masih bisa digunakan untuk prediksi.
Permasalahan prediksi ini menarik untuk dicoba menggunakan model prediksi lain, seperti logistic regression, svm, ann

LBB: Regression Model

Husain A. Hidayat

Januay 29, 2021

Introduction

Data Preparation

About the Data

Import Library and Data

Data Preparation

EDA

I. Persebaran data `Species`

II. Distribusi data numerik

III. Species manakah yang paling berat? (secara rataan)

IV. Korelasi antar data

V. Weight vs Height per `Species`

VI. Weight vs Width per `Species`

Modelling

Train-Test Split (Holdout Validation)

Model1 : Full predictor

Model2: Stepwise - Backward

Model3: Stepwise - Forward

Model4: Stepwise - Both Direction

Evaluation

Performance

Model Assumption

1. Linearity

2. Normality

3. Homocesdasticity

4. Multicollinearity

Model Improvement

Model Tuning

model tune 1

model tune 2

Performance

Evaluation

1. Linearity

2. Normality

3. Homocesdasticity

4. Multicolinearity

Conclusion

Annotations

LBB: Regression Model

Husain A. Hidayat

Januay 29, 2021

Introduction

Data Preparation

About the Data

Import Library and Data

Data Preparation

EDA

I. Persebaran data Species

II. Distribusi data numerik

III. Species manakah yang paling berat? (secara rataan)

IV. Korelasi antar data

V. Weight vs Height per Species

VI. Weight vs Width per Species

Modelling

Train-Test Split (Holdout Validation)

Model1 : Full predictor

Model2: Stepwise - Backward

Model3: Stepwise - Forward

Model4: Stepwise - Both Direction

Evaluation

Performance

Model Assumption

1. Linearity

2. Normality

3. Homocesdasticity

4. Multicollinearity

Model Improvement

Model Tuning

model tune 1

model tune 2

Performance

Evaluation

1. Linearity

2. Normality

3. Homocesdasticity

4. Multicolinearity

Conclusion

Annotations

I. Persebaran data `Species`

V. Weight vs Height per `Species`

VI. Weight vs Width per `Species`