Analisis Regresi

Annebel Diestya Clarissa

3/9/2021

Package yang digunakan

library(readxl)     #Membaca file data excel
library(kableExtra) #Tampilan Tabel
library(agricolae)  #Pemeriksaan Asumsi
library(lmtest)     #Untuk pengecekan asumsi
library(car)        #Untuk pengecekan asumsi
library(tseries)    #Untuk pengecekan asumsi
library (ggplot2)   #Visualisasi data
library(corrplot)   #Visualisasi MatrikS Korelasi

Data

Data yang akan digunakan adalah data Body Fat dengan 3 variabel independent dimana terdapat 2 variabel berkorelasi kuat.

3 Variabel independent tersebut adalah :

Triceps Skinfold Thickness
Thigh Circumference
Midarm Circumference

# Import Data
dataBodyFat <- read_xlsx("D:/Data.xlsx")

Berikut adalah data yang akan digunakan (Data terdiri dari 20 observasi dengan 3 variabel)

dataBodyFat %>% kbl(format = "html", caption= "Data Body Fat",align = 'c',longtable = 'T',) %>% kable_material_dark(full_width=F,c("striped", "hover", "condensed", "responsive")) %>% scroll_box(width="100%", height="400px")

Data Body Fat
Subject	Triceps	Thigh	Midarm	BodyFat
1	19.5	43.1	29.1	11.9
2	24.7	49.8	28.2	22.8
3	30.7	51.9	37.0	18.7
4	29.8	54.3	31.1	20.1
5	19.1	42.2	30.9	12.9
6	25.6	53.9	23.7	21.7
7	31.4	58.5	27.6	27.1
8	27.9	52.1	30.6	25.4
9	22.1	49.9	23.2	21.3
10	25.5	53.5	24.8	19.3
11	31.1	56.6	30.0	25.4
12	30.4	56.7	28.3	27.2
13	18.7	46.5	23.0	11.7
14	19.7	44.2	28.6	17.8
15	14.6	42.7	21.3	12.8
16	29.5	54.4	30.1	23.9
17	27.7	55.3	25.7	22.6
18	30.2	58.6	24.6	25.4
19	22.7	48.2	27.1	14.8
20	25.2	51.0	27.5	21.1

Visualisasi Data

Hubungan BodyFat dengan Triceps

scat1=ggplot(dataBodyFat,aes(Triceps,BodyFat))+geom_point()+geom_smooth(method="lm",se=T)

scat1

Berdasarkan pola diatas dimana titik-titik membentuk suatu garis lurus, diduga variabel Triceps memiliki hubungan dengan Variabel BodyFat

Hubungan BodyFat dengan Thigh

scat2=ggplot(dataBodyFat,aes(Thigh,BodyFat))+geom_point()+geom_smooth(method="lm",se=T)

scat2

Berdasarkan pola diatas dimana titik-titik membentuk suatu garis lurus, diduga variabel Thigh memiliki hubungan dengan Variabel BodyFat

Hubungan Bodyfat dengan Midarm

scat3=ggplot(dataBodyFat,aes(Midarm,BodyFat))+geom_point()+geom_smooth(method="lm",se=T)

scat3

Berdasarkan pola diatas, titik-titik menyebar dan tidak membentuk suatu garis lurus, diduga variabel Midarm tidak memiliki hubungan dengan Variabel BodyFat

Model Linier

Pemodelan regresi linier berganda akan dilakukan dengan membuat model yang mencakup keseluruhan variabel independent (Triceps, Thigh, dan Midarm).

# 3 Variabel : Triceps, Thigh, Midarm

lm1 <- lm(formula = BodyFat ~ Triceps + Thigh + Midarm, data=dataBodyFat)

summary(lm1)

## 
## Call:
## lm(formula = BodyFat ~ Triceps + Thigh + Midarm, data = dataBodyFat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.7263 -1.6111  0.3923  1.4656  4.1277 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  117.085     99.782   1.173    0.258
## Triceps        4.334      3.016   1.437    0.170
## Thigh         -2.857      2.582  -1.106    0.285
## Midarm        -2.186      1.595  -1.370    0.190
## 
## Residual standard error: 2.48 on 16 degrees of freedom
## Multiple R-squared:  0.8014, Adjusted R-squared:  0.7641 
## F-statistic: 21.52 on 3 and 16 DF,  p-value: 7.343e-06

Pengecekan Model

Dalam membangun model regresi linier, metode estimasi yang digunakan adalah dengan metode estimasi kuadrat terkecil (ordinary least square), yaitu dengan meminimumkan jumlah kuadrat residual.

Terdapat asumsi yang harus diuji dalam membangun model regresi linier tersebut. Asumsi model tersebut sering juga disebut sebagai asumsi klasik yang terdiri atas :

Uji normalitas residual
Uji asumsi variasi galat yang bersifat konstan (homoskedastisitas)
Uji asumsi tidak adanya serial korelasi dari galat (autokorelasi)
Uji multikolinieritas antarvariabel independen.

Uji kenormalan

Kenormalan sisaan akan diuji dengan Shapiro Wilk Test , Uji Jarque Bera, dan qqplot.

Uji Shapiro Wilk

shapiro.test(residuals(lm1))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(lm1)
## W = 0.96603, p-value = 0.6698

Hipotesis

H0 : Residual berdistribusi normal
H1 : Residual tidak berdistribusi normal

Tingkat signifikansi = 0.05

Statistik Uji : p-value = 0.6698

Keputusan : Karena p-value (0.6698) > alpha (0.05), maka gagal tolak H0

Kesimpulan : Dengan tingkat keyakinan 95%, kita yakin bahwa residual berdistribusi normal.

Uji Jarque Bera

# analisis residual
residual=resid(lm1) #menampilkan residual dari model 2

jarque.bera.test(residual) #normalitas

## 
##  Jarque Bera Test
## 
## data:  residual
## X-squared = 0.83048, df = 2, p-value = 0.6602

Hipotesis

H0 : Residual berdistribusi normal
H1 : Residual tidak berdistribusi normal

Tingkat signifikansi = 0.05

Statistik Uji : p-value = 0.6602

Keputusan : Karena p-value (0.6602) > alpha (0.05), maka gagal tolak H0

Kesimpulan : Dengan tingkat keyakinan 95%, kita yakin bahwa residual berdistribusi normal.

qq plot

qqPlot(residual,distribution="norm",main="Normal QQ Plot") #normalitas dengan plot

## [1] 14 19

qqplot menunjukkan residual ada di sepanjang garis 45 derajat, sehingga mengindikasikan bahwa residual menyebar normal.

Dilihat dari uji Shapiro Wilk, Jarque Bera, dan qqplot, maka dapat disimpulkan bahwa model yang dibuat memenuhi asumsi kenormalan sisaan.

Kehomogenan Ragam (Homoskedastisitas)

Untuk menguji homoskedastisitas akan digunakan uji Breusch Pagan

bptest(lm1, data=dataBodyFat)

## 
##  studentized Breusch-Pagan test
## 
## data:  lm1
## BP = 5.1452, df = 3, p-value = 0.1615

Hipotesis

H0 : Ragam Homogen
H1 : Ragam tidak homogen

Tingkat Signifikansi : 0.05

Statistik Uji : pvalue=0.1615

Keputusan : Karena p-value (0.1615) > alpha (0.05) maka gagal tolak H0.

Kesimpulan : Dengan tingkat keyakinan 95%, kita yakin bahwa ragam homogen (asumsi homoskedastisitas terpenuhi)

Autokorelasi

Uji autokorelasi akan diuji dengan uji Durbin Watson

dwtest(lm1)

## 
##  Durbin-Watson test
## 
## data:  lm1
## DW = 2.2429, p-value = 0.6698
## alternative hypothesis: true autocorrelation is greater than 0

Hipotesis

H0 : Tidak ada auto korelasi dari galat
H1 : Ada auto korelasi dari galat

Tingkat signifikansi : 0.05

Statistik uji : pvalue=0.6698

Keputusan : Karena nilai pvalue (0.6698) > alpha (0.05) maka gagal tolak H0

Kesimpulan : Dengan tingkat kepercayaan 95%, kita yakin bahwa tidak ada auto korelasi pada galat.

Pencilan

Pencilan dapat dideteksi dengan menggunakan boxplot dan scatter plot

# boxplot 
boxplot((dataBodyFat[,5]), main="Body Fat")

# scatterplot 
plot(dataBodyFat$BodyFat)

# Mengeluarkan nilai yang outlier
nilai_outlier = boxplot.stats(dataBodyFat$BodyFat)$out
nilai_outlier

## numeric(0)

Dilihat dari hasil boxplot, scatter, dan hasil perhitungan, tidak ditemukan adanya pencilan

Multikolinieritas

Kolinieritas akan diuji dengan VIF (Variance Inflation Factor). Jika, nilai VIF > 5, maka terjadi multikolinieritas.

vif(lm1)

##  Triceps    Thigh   Midarm 
## 708.8429 564.3434 104.6060

Ketiga variabel (Triceps, Thigh, dan Midarm) memiliki nilai VIF > 5 dan nilainya sangat besar sehingga mengindikasikan bahwa ada multikolinieritas

Matriks Korelasi

Uji multikolinieritas dengan VIF > 5 menunjukkan adanya multikolinieritas. Pada data juga sudah diketahui bahwa ada 2 variabel yang berkorelasi kuat. Akan diuji variabel mana yang berkorelasi dengan menggunakan matrik korelasi

Perhitungan Matriks Korelasi

bodyfat <- dataBodyFat[c(2,3,4,5)]

mk <- cor(bodyfat)

round(mk,2) #Membulatkan 2 angka di belakang koma

##         Triceps Thigh Midarm BodyFat
## Triceps    1.00  0.92   0.46    0.84
## Thigh      0.92  1.00   0.08    0.88
## Midarm     0.46  0.08   1.00    0.14
## BodyFat    0.84  0.88   0.14    1.00

Visualisasi Matriks Korelasi

corrplot(mk, type="lower",
         order = "hclust", # mengurutkan berdasarkan hierarchical clustering
         tl.col= "black", # warna tulisan
         addCoef.col = "black", # tambahkan koefisien korelasi
         diag=FALSE, #menyembunyikan koefisien pada diagonal
         tl.srt= 45, # kemiringan tulisan 45 derajat
         method = "circle") # Bentuk Visualisasi

Terlihat dalam matriks korelasi bahwa terdapat korelasi yang tinggi antara Triceps dan Thigh sebesar 0.92 yang menyebabkan adanya kolinieritas.

Solusi

Untuk mengatasi masalah multikolinieritas dalam model maka perlu dibuat model yang baru yang tidak mengandung variabel yang berkorelasi tinggi.

Cara yang dilakukan yaitu dengan mengeluarkan variabel yang memiliki p-value terbesar serta tidak signifikan. Begitu seterusnya hingga semua variabel telah signifikan.

Model yang kita punya sebelumnya adalah

Model 1 - Model Awal

# 3 Variabel : Triceps, Thigh, Midarm

lm1 <- lm(formula = BodyFat ~ Triceps + Thigh + Midarm, data=dataBodyFat)

summary(lm1)

## 
## Call:
## lm(formula = BodyFat ~ Triceps + Thigh + Midarm, data = dataBodyFat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.7263 -1.6111  0.3923  1.4656  4.1277 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  117.085     99.782   1.173    0.258
## Triceps        4.334      3.016   1.437    0.170
## Thigh         -2.857      2.582  -1.106    0.285
## Midarm        -2.186      1.595  -1.370    0.190
## 
## Residual standard error: 2.48 on 16 degrees of freedom
## Multiple R-squared:  0.8014, Adjusted R-squared:  0.7641 
## F-statistic: 21.52 on 3 and 16 DF,  p-value: 7.343e-06

Jika dilihat dari summary model 1, p-value terbesar adalah variabel Thigh yaitu 0.285 sehingga akan dikeluarkan dari model dan dibuat Model ke-2 dengan 2 variabel yaitu Triceps dan Midarm.

Model 2 - Triceps & Midarm

# 2 Variabel : Triceps, Midarm

lm2 <- lm(formula = BodyFat ~ Triceps + Midarm, data=dataBodyFat)

summary(lm2)

## 
## Call:
## lm(formula = BodyFat ~ Triceps + Midarm, data = dataBodyFat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.8794 -1.9627  0.3811  1.2688  3.8942 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   6.7916     4.4883   1.513   0.1486    
## Triceps       1.0006     0.1282   7.803 5.12e-07 ***
## Midarm       -0.4314     0.1766  -2.443   0.0258 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.496 on 17 degrees of freedom
## Multiple R-squared:  0.7862, Adjusted R-squared:  0.761 
## F-statistic: 31.25 on 2 and 17 DF,  p-value: 2.022e-06

Selanjutnya akan dibuat 2 model lagi dengan 1 variabel yaitu Triceps dan Thigh. Pertimbangannya adalah dengan melihat matriks korelasi, kedua variabel ini memiliki korelasi yang tinggi dengan variabel respon.

Model 3 - Triceps

Akan dibentuk model dengan 1 variabel Triceps karena memiliki nilai korelasi 0.84 dengan variabel respon Bodyfat.

# 1 Variabel : Triceps

lm3 <- lm(formula = BodyFat ~ Triceps , data=dataBodyFat)

summary(lm3)

## 
## Call:
## lm(formula = BodyFat ~ Triceps, data = dataBodyFat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.1195 -2.1904  0.6735  1.9383  3.8523 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -1.4961     3.3192  -0.451    0.658    
## Triceps       0.8572     0.1288   6.656 3.02e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.82 on 18 degrees of freedom
## Multiple R-squared:  0.7111, Adjusted R-squared:  0.695 
## F-statistic:  44.3 on 1 and 18 DF,  p-value: 3.024e-06

Model 4 - Thigh

Akan dibentuk model dengan 1 variabel Thigh karena memiliki nilai korelasi 0.88 dengan variabel respon Bodyfat.

# 1 Variabel : Thigh

lm4 <- lm(formula = BodyFat ~ Thigh , data=dataBodyFat)

summary(lm4)

## 
## Call:
## lm(formula = BodyFat ~ Thigh, data = dataBodyFat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.4949 -1.5671  0.1241  1.3362  4.4084 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -23.6345     5.6574  -4.178 0.000566 ***
## Thigh         0.8565     0.1100   7.786  3.6e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.51 on 18 degrees of freedom
## Multiple R-squared:  0.771,  Adjusted R-squared:  0.7583 
## F-statistic: 60.62 on 1 and 18 DF,  p-value: 3.6e-07

Setelah terbentuk 4 model regresi, dilanjutkan dengan memilih model yang terbaik dengan melihat beberapa kriterianya.

Kriteria yang digunakan adalah AIC, BIC, Adjusted Rsquared, dan Residual Standard Error.

Metode AIC dan BIC

Metode AIC dan BIC adalah metode yang dapat digunakan untuk memilih model regresi terbaik yang ditemukan oleh Akaike dan Schwarz (Grasa,1989). Kedua metode ini didasarkan pada metode maximum likelihood estimation (MLE). Menurut metode AIC dan BIC, model regresi terbaik adalahmodel regresi yang mempunyai nilai AIC dan BIC terkecil`.

# kriteria pemilihan model
AIC=c(AIC(lm1),AIC(lm2),AIC(lm3),AIC(lm4))
BIC=c(BIC(lm1),BIC(lm2),BIC(lm3),BIC(lm4))
AdjustedRsquared=c(summary(lm1)[[9]],summary(lm2)[[9]],
                   summary(lm3)[[9]],summary(lm4)[[9]])
ResidualStandardError=c(summary(lm1)[[6]],summary(lm2)[[6]],
                        summary(lm3)[[6]],summary(lm4)[[6]])
Model=c("LM1","LM2","LM3","LM4")

Kriteria=data.frame(Model,AdjustedRsquared,
                    ResidualStandardError,AIC,BIC)
Kriteria

##   Model AdjustedRsquared ResidualStandardError       AIC      BIC
## 1   LM1        0.7641133              2.479981  98.62471 103.6034
## 2   LM2        0.7610022              2.496282  98.09925 102.0822
## 3   LM3        0.6950464              2.819769 102.11652 105.1037
## 4   LM4        0.7583215              2.510242  97.46550 100.4527

Tabel di atas menunjukkan nilai masing-masing kriteria untuk setiap model. Dapat disimpulkan bahwa model 4 adalah model yang terbaik karena memiliki kriteria yang lebih baik dari model lainnya.

Sehingga, dipilih model 4 dengan 1 variabel independent yaitu Thigh yang akan dilakukan pengecekan asumsi klasik.

Pengecekan Asumsi Model Baru

par(mfrow=c(2,2))
plot(lm4)

par(mfrow=c(1,1))

Uji kenormalan

Kenormalan sisaan akan diuji dengan Shapiro Wilk Test , Uji Jarque Bera, dan qqplot.

Uji Shapiro Wilk

shapiro.test(residuals(lm4))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(lm4)
## W = 0.97401, p-value = 0.8363

Hipotesis

H0 : Residual berdistribusi normal
H1 : Residual tidak berdistribusi normal

Tingkat signifikansi = 0.05

Statistik Uji : p-value = 0.8363

Keputusan : Karena p-value (08363) > alpha (0.05), maka gagal tolak H0

Kesimpulan : Dengan tingkat keyakinan 95%, kita yakin bahwa residual berdistribusi normal.

Uji Jarque Bera

# analisis residual
residual=resid(lm4) #menampilkan residual dari model 2

jarque.bera.test(residual) #normalitas

## 
##  Jarque Bera Test
## 
## data:  residual
## X-squared = 0.56188, df = 2, p-value = 0.7551

Hipotesis

H0 : Residual berdistribusi normal
H1 : Residual tidak berdistribusi normal

Tingkat signifikansi = 0.05

Statistik Uji : p-value = 0.7551

Keputusan : Karena p-value (0.7551) > alpha (0.05), maka gagal tolak H0

Kesimpulan : Dengan tingkat keyakinan 95%, kita yakin bahwa residual berdistribusi normal.

qq plot

qqPlot(residual,distribution="norm",main="Normal QQ Plot") #normalitas dengan plot

## [1] 13  8

qqplot menunjukkan error ada di sepanjang garis 45 derajat, sehingga mengindikasikan bahwa residual menyebar normal.

Dilihat dari uji Shapiro Wilk, Jarque Bera, dan qqplot, maka dapat disimpulkan bahwa model yang dibuat memenuhi asumsi kenormalan sisaan.

Kehomogenan Ragam (Homoskedastisitas)

Untuk menguji homoskedastisitas akan digunakan uji Breusch Pagan

bptest(lm4, data=dataBodyFat)

## 
##  studentized Breusch-Pagan test
## 
## data:  lm4
## BP = 0.85834, df = 1, p-value = 0.3542

Hipotesis

H0 : Ragam Homogen
H1 : Ragam tidak homogen

Tingkat Signifikansi : 0.05

Statistik Uji : pvalue=0.3542

Keputusan : Karena p-value (0.3542) > alpha (0.05) maka gagal tolak H0.

Kesimpulan : Dengan tingkat keyakinan 95%, kita yakin bahwa ragam homogen (asumsi homoskedastisitas terpenuhi)

Autokorelasi

Uji autokorelasi akan diuji dengan uji Durbin Watson

dwtest(lm4)

## 
##  Durbin-Watson test
## 
## data:  lm4
## DW = 2.5001, p-value = 0.8604
## alternative hypothesis: true autocorrelation is greater than 0

Hipotesis

H0 : Tidak ada auto korelasi dari galat
H1 : Ada auto korelasi dari galat

Tingkat signifikansi : 0.05

Statistik uji : pvalue=0.8604

Keputusan : Karena nilai pvalue (0.7851) > alpha (0.05) maka gagal tolak H0

Kesimpulan : Dengan tingkat kepercayaan 95%, kita yakin bahwa tidak ada auto korelasi pada galat.

Uji-uji asumsi menunjukkan bahwa solusi menghilangkan variabel Triceps dan Midarm pada model cukup berhasil untuk menghasilkan model yang baik.

Persamaan Regresi

lm4 <- lm(formula = BodyFat ~ Thigh, data=dataBodyFat)

summary(lm4)

## 
## Call:
## lm(formula = BodyFat ~ Thigh, data = dataBodyFat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.4949 -1.5671  0.1241  1.3362  4.4084 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -23.6345     5.6574  -4.178 0.000566 ***
## Thigh         0.8565     0.1100   7.786  3.6e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.51 on 18 degrees of freedom
## Multiple R-squared:  0.771,  Adjusted R-squared:  0.7583 
## F-statistic: 60.62 on 1 and 18 DF,  p-value: 3.6e-07

Persamaan Regresi

BodyFat = -23.6345 + (0.8565) Thigh

Interpretasi :

Setiap kenaikan satu satuan Thigh menyebabkan Body Fat naik sebesar 0.8565 dengan mengganggap variabel lain konstan.

Multiple R-squared = 0.771 yang artinya BodyFat dapat dijelaskan oleh Thigh sebesar 77.1% sedangkan 22.9% dijelaskan oleh faktor lain diluar model.