Model Regresi

Soal Anstat

Data

Data yang digunakan adalah data Body Fat dengan 3 variabel independent

3 Variabel independent tersebut adalah :
- Triceps Skinfold Thickness
- Thigh Circumference
- Midarm Circumference

DataBodyFat <- read.xlsx("BodyFat.xlsx") # reading file from class.
kbl(DataBodyFat)%>%  kable_styling() %>% scroll_box(height = "300px")

Triceps	Thigh	Midarm	BodyFat
19.5	43.1	29.1	11.9
24.7	49.8	28.2	22.8
30.7	51.9	37.0	18.7
29.8	54.3	31.1	20.1
19.1	42.2	30.9	12.9
25.6	53.9	23.7	21.7
31.4	58.5	27.6	27.1
27.9	52.1	30.6	25.4
22.1	49.9	23.2	21.3
25.5	53.5	24.8	19.3
31.1	56.6	30.0	25.4
30.4	56.7	28.3	27.2
18.7	46.5	23.0	11.7
19.7	44.2	28.6	17.8
14.6	42.7	21.3	12.8
29.5	54.4	30.1	23.9
27.7	55.3	25.7	22.6
30.2	58.6	24.6	25.4
22.7	48.2	27.1	14.8
25.2	51.0	27.5	21.1

str(DataBodyFat)

## 'data.frame':    20 obs. of  4 variables:
##  $ Triceps: num  19.5 24.7 30.7 29.8 19.1 25.6 31.4 27.9 22.1 25.5 ...
##  $ Thigh  : num  43.1 49.8 51.9 54.3 42.2 53.9 58.5 52.1 49.9 53.5 ...
##  $ Midarm : num  29.1 28.2 37 31.1 30.9 23.7 27.6 30.6 23.2 24.8 ...
##  $ BodyFat: num  11.9 22.8 18.7 20.1 12.9 21.7 27.1 25.4 21.3 19.3 ...

DataBodyFat<-DataBodyFat

Plot Data

ggpairs(DataBodyFat)

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

Model Linier (3 Variabel)

Model1 <- lm(formula = BodyFat ~ Triceps + Thigh + Midarm, data=DataBodyFat)
summary(Model1)

## 
## 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

Model12 <- lm(Triceps ~  Thigh + Midarm, data=DataBodyFat)
car::vif(Model1)

##  Triceps    Thigh   Midarm 
## 708.8429 564.3434 104.6060

Dari hasil pemodelan tersebut didapatkan nilai Adjusted R-squared sebesar 0.7641. Nilai tersebut menunjukkan adanya hubungan linear dalam model tersebut.

Pengujian parameter koefisien β

H0 : β1 = β2 = β3 = 0
H1 : setidaknya ada satu βi yang tidak sama dengan 0

Dari hasil pengujian tersebut, terlihat bahwa nilai p sangat kecil yaitu 7.343e-06. Sehingga cukup bukti untuk menolah Ho dan dapat disimpulkan bahwa secara simultan model signifikan mempengaruhi

Uji Parsial

Ho : βi = 0 H1 : βi ≠ 0

Hasil uji pada masing masing koefisien menunjukkan nilai t yang kecil sehingga tidak signifikan untuk masing masing koefisien.

Ketika R-Squared tinggi, tetapi koefisien tidak nyata, maka patut dicuragai adanya kolinearitas antara variabel penjelas.

Pemeriksaan Model

VIF

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

vif(Model1)

##  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

BodyFat <- DataBodyFat
mk <- cor(BodyFat)
round(mk,2)

##         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

Jika dilihat dari matriksnya, ketiga variabel ini secara berpasangan saling berkorelasi.

ggcorr(DataBodyFat,label = T, size=3, label_size = 3, hjust=0.95)+
  labs(
    title="Matriks Korelasi Data set"
  )+
  theme_minimal()+
  theme(
    plot.title = element_text(hjust = 0.5),
    axis.title=element_text(size=8,face="bold"), 
    axis.text.y=element_blank()
  )

Linieritas

residuals=resid(Model1)
linearity <- data.frame(residual = residuals, fitted = DataBodyFat$BodyFat)
linearity %>% 
  ggplot(aes(fitted, residual)) + 
  geom_point() + 
  geom_smooth(method = lm) + 
  geom_hline(aes(yintercept = 0)) + 
     theme_minimal()+
  labs(
    title="Residual Vs Fitted", 
    caption = "")

## `geom_smooth()` using formula 'y ~ x'

Normalitas Residual

Kenormalan sisaan akan diuji dengan qqplot dan Shapiro Wilk Test

qq plot

qqPlot(residuals,distribution="norm",main="Normal QQ Plot")

## [1] 14 19

Uji Shapiro Wilk

shapiro.test(residuals)

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

Kehomogenan Ragam

Untuk menguji Homoskedastisitas akan digunakan uji Breusch Pagan

bptest(Model1, data=DataBodyFat)

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

homogen <- data.frame(residual = residuals, fitted = DataBodyFat$BodyFat)
homogen %>% 
  ggplot(aes(fitted, residual)) + 
  geom_point() + 
  geom_smooth(method = "h") + 
  geom_hline(aes(yintercept = 0)) + 
     theme_minimal()+
    labs(
    title="Residual Vs Fitted", 
    caption = "")

## `geom_smooth()` using formula 'y ~ x'

## Warning: Computation failed in `stat_smooth()`:
## object 'h' of mode 'function' was not found

Outlier

Hat Values

hatvalues(Model1)

##          1          2          3          4          5          6          7 
## 0.34120920 0.15653638 0.44042770 0.11242972 0.36109984 0.13151364 0.19433721 
##          8          9         10         11         12         13         14 
## 0.16418081 0.19278940 0.24051819 0.13935816 0.10929380 0.21357666 0.18808377 
##         15         16         17         18         19         20 
## 0.34830629 0.11439069 0.12532943 0.22828343 0.13235798 0.06597771

hatval = data.frame(hatv=hatvalues(Model1)) %>% mutate(idx=1:20)
hatval

##          hatv idx
## 1  0.34120920   1
## 2  0.15653638   2
## 3  0.44042770   3
## 4  0.11242972   4
## 5  0.36109984   5
## 6  0.13151364   6
## 7  0.19433721   7
## 8  0.16418081   8
## 9  0.19278940   9
## 10 0.24051819  10
## 11 0.13935816  11
## 12 0.10929380  12
## 13 0.21357666  13
## 14 0.18808377  14
## 15 0.34830629  15
## 16 0.11439069  16
## 17 0.12532943  17
## 18 0.22828343  18
## 19 0.13235798  19
## 20 0.06597771  20

plot = ggplot(data=hatval,
  aes(x=idx,y=hatv))+ylim(0,0.8)+
  geom_point()+ 
  labs(x ="Observation",y="Hat Value ke-i", title = "Hat Values") + 
  geom_abline(aes(intercept =12/20,slope=0))
plot

Berdasarkan nilai HatValues, tidak ditemukan adanya observation yang melebihi dari 12/20 sehingga tidak mengandung pencilan

Residual Plot

residuals=residuals(Model1)
par(mfrow=c(2,2))
plot(Model1)

Perbaikan Model

Dengan menggunakan regsubsets maka model yang diajukan ada dua model, yaitu model dengan 1 variabel thigh dan model dengan 2 variabel prediktor yaitu triceps dan midarm

models <- regsubsets(BodyFat~.,data=DataBodyFat, nvmax = 2)
summary(models)

## Subset selection object
## Call: regsubsets.formula(BodyFat ~ ., data = DataBodyFat, nvmax = 2)
## 3 Variables  (and intercept)
##         Forced in Forced out
## Triceps     FALSE      FALSE
## Thigh       FALSE      FALSE
## Midarm      FALSE      FALSE
## 1 subsets of each size up to 2
## Selection Algorithm: exhaustive
##          Triceps Thigh Midarm
## 1  ( 1 ) " "     "*"   " "   
## 2  ( 1 ) "*"     " "   "*"

res.sum <- summary(models)
data.frame(
  Adj.R2 = which.max(res.sum$adjr2),
  CP = which.min(res.sum$cp),
  BIC = which.min(res.sum$bic)
)

##   Adj.R2 CP BIC
## 1      2  1   1

# Pemodelan regresi
Model2 <- lm(formula = BodyFat ~ Triceps + Midarm, data=DataBodyFat)
Model3<- lm(formula = BodyFat ~ Thigh , data=DataBodyFat)

# kriteria pemilihan model
AIC=c(AIC(Model1),AIC(Model2),AIC(Model3))
BIC=c(BIC(Model1),BIC(Model2),BIC(Model3))
AdjRsquared=c(summary(Model1)[[9]],summary(Model2)[[9]],summary(Model3)[[9]])
RSE=c(summary(Model1)[[6]],summary(Model2)[[6]],summary(Model3)[[6]])
Model=c("Model1","Model2","Model3")

Kriteria=data.frame(Model,AdjRsquared,RSE,AIC,BIC)
Kriteria

##    Model AdjRsquared      RSE      AIC      BIC
## 1 Model1   0.7641133 2.479981 98.62471 103.6034
## 2 Model2   0.7610022 2.496282 98.09925 102.0822
## 3 Model3   0.7583215 2.510242 97.46550 100.4527

Apabila menginginkan model dengan R-squared yang lebih besar maka dapat menggunakan 2 variabel prediktor, namun jika mencari nilai CP dan BIC yang lebih kecil dapat menggunakan satu variabel prediktor.

Pemilihan Model Baru

Model 2 Variabel

Dengan menggunakan regsubsets salah satu alternatif yang diberikan adalah model dengan 2 variabel prediktor yaitu triceps dan midarms

# Pemodelan regresi
Model2 <- lm(formula = BodyFat ~ Triceps + Midarm, data=DataBodyFat)
summary(Model2)

## 
## 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

Model tersebut dapat dituliskan dengan persamaan:

Y = β0 + β1X1 + β2X2

Body fat = 6.7916 + 1.006 * Triceps Skinfold Thickness - 0.4314 * Midarm Circumference

Model memiliki nilai Adjusted R-Squared sebesar 0.761 dan nilai BIC sebesar 102.0822

Residual Plot

residuals=residuals(Model2)
par(mfrow=c(2,2))
plot(Model2)

Model 1 Variabel

Rekomendasi lain yang adalah memasukkan satu variabel prediktor yaitu Thigh. Sehingga modelnya menjadi:

# Pemodelan regresi
Model3<- lm(formula = BodyFat ~ Thigh , data=DataBodyFat)
summary(Model3)

## 
## 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

Model tersebut dapat dituliskan dengan persamaan:

Y = β0 + β1X1

Body fat = -23.6345 + 0.8565 * Thigh Circumference

Model memiliki nilai Adjusted R-Squared sebesar 0.7582 dan nilai BIC sebesar 100.4527

Residual Plot

residuals=residuals(Model3)
par(mfrow=c(2,2))
plot(Model3)

Linieritas

residuals=resid(Model3)
linearity <- data.frame(residual = residuals, fitted = DataBodyFat$BodyFat)
linearity %>% 
  ggplot(aes(fitted, residual)) + 
  geom_point() + 
  geom_smooth(method = lm) + 
  geom_hline(aes(yintercept = 0)) + 
     theme_minimal()+
  labs(
    title="Residual Vs Fitted", 
    caption = "")

## `geom_smooth()` using formula 'y ~ x'

Normalitas

Kenormalan sisaan akan diuji dengan Shapiro Wilk Test dan qqplot.

qq plot

qqPlot(residuals,distribution="norm",main="Normal QQ Plot")

## [1] 13  8

Uji Shapiro Wilk

shapiro.test(residuals)

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

Kehomogenan Ragam

Untuk menguji Homoskedastisitas akan digunakan uji Breusch Pagan

bptest(Model3, data=DataBodyFat)

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

homogen <- data.frame(residual = residuals, fitted = DataBodyFat$BodyFat)
homogen %>% 
  ggplot(aes(fitted, residual)) + 
  geom_point() + 
  geom_smooth(method = "h") + 
  geom_hline(aes(yintercept = 0)) + 
     theme_minimal()+
    labs(
    title="Residual Vs Fitted", 
    caption = "")

## `geom_smooth()` using formula 'y ~ x'

## Warning: Computation failed in `stat_smooth()`:
## object 'h' of mode 'function' was not found

Outlier

Hat Values

hatvalues(Model3)

##          1          2          3          4          5          6          7 
## 0.17509056 0.05360511 0.05102358 0.06881768 0.20454764 0.06431538 0.15320136 
##          8          9         10         11         12         13         14 
## 0.05166128 0.05309803 0.06042772 0.10663399 0.10873916 0.09189009 0.14331319 
##         15         16         17         18         19         20 
## 0.18779844 0.07003930 0.08276254 0.15603643 0.06694300 0.05005551

hatva = data.frame(hatv=hatvalues(Model3)) %>% mutate(idx=1:20)
hatva

##          hatv idx
## 1  0.17509056   1
## 2  0.05360511   2
## 3  0.05102358   3
## 4  0.06881768   4
## 5  0.20454764   5
## 6  0.06431538   6
## 7  0.15320136   7
## 8  0.05166128   8
## 9  0.05309803   9
## 10 0.06042772  10
## 11 0.10663399  11
## 12 0.10873916  12
## 13 0.09189009  13
## 14 0.14331319  14
## 15 0.18779844  15
## 16 0.07003930  16
## 17 0.08276254  17
## 18 0.15603643  18
## 19 0.06694300  19
## 20 0.05005551  20

plot = ggplot(data=hatva,
  aes(x=idx,y=hatv))+ylim(0,0.4)+
  geom_point()+ 
  labs(x ="Observation",y="Hat Value ke-i", title = "Hat Values") + 
  geom_abline(aes(intercept =6/20,slope=0))
plot