Tugas Mandiri PSD
Aida Darajati G1401211016
2023-09-28
Packages
lapply(c("glmnet","lmridge"),library,character.only=T)[[1]]## Warning: package 'glmnet' was built under R version 4.2.3## Loading required package: Matrix## Loaded glmnet 4.1-8## [1] "glmnet" "Matrix" "stats" "graphics" "grDevices" "utils"
## [7] "datasets" "methods" "base"Input data Chiken Weight dari R
ChickWeight Dataset:
Dataset ChickWeight berisi data berat badan ayam selama beberapa periode waktu. Dataset ini sering digunakan untuk mengilustrasikan konsep analisis data longitudinal.
Variabel:
weight: Berat badan ayam (dalam gram).
Time: Waktu (hari).
Chick: Nomor identifikasi ayam.
Diet: Diet yang diberikan kepada ayam (1-4).
dataayam <- ChickWeight
head(dataayam)## weight Time Chick Diet
## 1 42 0 1 1
## 2 51 2 1 1
## 3 59 4 1 1
## 4 64 6 1 1
## 5 76 8 1 1
## 6 93 10 1 1Seleksi Peubah
library(leaps)forward <- regsubsets(weight~Time+Chick+Diet,data=dataayam, method="forward")## Warning in leaps.setup(x, y, wt = wt, nbest = nbest, nvmax = nvmax, force.in =
## force.in, : 3 linear dependencies foundsf<-summary(forward)
sf## Subset selection object
## Call: regsubsets.formula(weight ~ Time + Chick + Diet, data = dataayam,
## method = "forward")
## 53 Variables (and intercept)
## Forced in Forced out
## Time FALSE FALSE
## Chick.L FALSE FALSE
## Chick.Q FALSE FALSE
## Chick.C FALSE FALSE
## Chick^4 FALSE FALSE
## Chick^5 FALSE FALSE
## Chick^6 FALSE FALSE
## Chick^7 FALSE FALSE
## Chick^8 FALSE FALSE
## Chick^9 FALSE FALSE
## Chick^10 FALSE FALSE
## Chick^11 FALSE FALSE
## Chick^12 FALSE FALSE
## Chick^13 FALSE FALSE
## Chick^14 FALSE FALSE
## Chick^15 FALSE FALSE
## Chick^16 FALSE FALSE
## Chick^17 FALSE FALSE
## Chick^18 FALSE FALSE
## Chick^19 FALSE FALSE
## Chick^20 FALSE FALSE
## Chick^21 FALSE FALSE
## Chick^22 FALSE FALSE
## Chick^23 FALSE FALSE
## Chick^24 FALSE FALSE
## Chick^25 FALSE FALSE
## Chick^26 FALSE FALSE
## Chick^27 FALSE FALSE
## Chick^28 FALSE FALSE
## Chick^29 FALSE FALSE
## Chick^30 FALSE FALSE
## Chick^31 FALSE FALSE
## Chick^32 FALSE FALSE
## Chick^33 FALSE FALSE
## Chick^34 FALSE FALSE
## Chick^35 FALSE FALSE
## Chick^36 FALSE FALSE
## Chick^37 FALSE FALSE
## Chick^38 FALSE FALSE
## Chick^39 FALSE FALSE
## Chick^40 FALSE FALSE
## Chick^41 FALSE FALSE
## Chick^42 FALSE FALSE
## Chick^43 FALSE FALSE
## Chick^44 FALSE FALSE
## Chick^45 FALSE FALSE
## Chick^46 FALSE FALSE
## Chick^47 FALSE FALSE
## Chick^48 FALSE FALSE
## Chick^49 FALSE FALSE
## Diet2 FALSE FALSE
## Diet3 FALSE FALSE
## Diet4 FALSE FALSE
## 1 subsets of each size up to 8
## Selection Algorithm: forward
## Time Chick.L Chick.Q Chick.C Chick^4 Chick^5 Chick^6 Chick^7 Chick^8
## 1 ( 1 ) "*" " " " " " " " " " " " " " " " "
## 2 ( 1 ) "*" "*" " " " " " " " " " " " " " "
## 3 ( 1 ) "*" "*" " " " " " " " " " " " " " "
## 4 ( 1 ) "*" "*" " " " " " " " " " " " " " "
## 5 ( 1 ) "*" "*" " " " " " " " " "*" " " " "
## 6 ( 1 ) "*" "*" " " " " " " " " "*" " " " "
## 7 ( 1 ) "*" "*" " " " " " " " " "*" " " " "
## 8 ( 1 ) "*" "*" " " "*" " " " " "*" " " " "
## Chick^9 Chick^10 Chick^11 Chick^12 Chick^13 Chick^14 Chick^15 Chick^16
## 1 ( 1 ) " " " " " " " " " " " " " " " "
## 2 ( 1 ) " " " " " " " " " " " " " " " "
## 3 ( 1 ) " " " " " " " " " " " " " " " "
## 4 ( 1 ) " " " " " " " " "*" " " " " " "
## 5 ( 1 ) " " " " " " " " "*" " " " " " "
## 6 ( 1 ) " " " " " " " " "*" " " " " " "
## 7 ( 1 ) " " " " " " " " "*" " " " " " "
## 8 ( 1 ) " " " " " " " " "*" " " " " " "
## Chick^17 Chick^18 Chick^19 Chick^20 Chick^21 Chick^22 Chick^23
## 1 ( 1 ) " " " " " " " " " " " " " "
## 2 ( 1 ) " " " " " " " " " " " " " "
## 3 ( 1 ) " " " " " " " " " " " " " "
## 4 ( 1 ) " " " " " " " " " " " " " "
## 5 ( 1 ) " " " " " " " " " " " " " "
## 6 ( 1 ) " " " " " " " " " " " " " "
## 7 ( 1 ) " " " " " " " " " " " " " "
## 8 ( 1 ) " " " " " " " " " " " " " "
## Chick^24 Chick^25 Chick^26 Chick^27 Chick^28 Chick^29 Chick^30
## 1 ( 1 ) " " " " " " " " " " " " " "
## 2 ( 1 ) " " " " " " " " " " " " " "
## 3 ( 1 ) " " " " " " " " " " " " " "
## 4 ( 1 ) " " " " " " " " " " " " " "
## 5 ( 1 ) " " " " " " " " " " " " " "
## 6 ( 1 ) " " " " " " " " " " " " " "
## 7 ( 1 ) " " " " " " " " " " " " " "
## 8 ( 1 ) " " " " " " " " " " " " " "
## Chick^31 Chick^32 Chick^33 Chick^34 Chick^35 Chick^36 Chick^37
## 1 ( 1 ) " " " " " " " " " " " " " "
## 2 ( 1 ) " " " " " " " " " " " " " "
## 3 ( 1 ) " " " " " " " " " " " " " "
## 4 ( 1 ) " " " " " " " " " " " " " "
## 5 ( 1 ) " " " " " " " " " " " " " "
## 6 ( 1 ) " " " " " " " " " " " " " "
## 7 ( 1 ) " " " " " " " " " " " " " "
## 8 ( 1 ) " " " " " " " " " " " " " "
## Chick^38 Chick^39 Chick^40 Chick^41 Chick^42 Chick^43 Chick^44
## 1 ( 1 ) " " " " " " " " " " " " " "
## 2 ( 1 ) " " " " " " " " " " " " " "
## 3 ( 1 ) " " " " " " " " " " " " " "
## 4 ( 1 ) " " " " " " " " " " " " " "
## 5 ( 1 ) " " " " " " " " " " " " " "
## 6 ( 1 ) " " " " " " " " " " " " " "
## 7 ( 1 ) " " " " " " " " " " " " " "
## 8 ( 1 ) " " " " " " " " " " " " " "
## Chick^45 Chick^46 Chick^47 Chick^48 Chick^49 Diet2 Diet3 Diet4
## 1 ( 1 ) " " " " " " " " " " " " " " " "
## 2 ( 1 ) " " " " " " " " " " " " " " " "
## 3 ( 1 ) " " " " " " " " " " " " " " "*"
## 4 ( 1 ) " " " " " " " " " " " " " " "*"
## 5 ( 1 ) " " " " " " " " " " " " " " "*"
## 6 ( 1 ) " " " " " " " " " " " " "*" "*"
## 7 ( 1 ) " " " " " " " " " " "*" "*" "*"
## 8 ( 1 ) " " " " " " " " " " "*" "*" "*"par(mfrow=c(2,2))
plot(sf$rsq, xlab="Jumlah Variabel", ylab="R2",type="l");points(1:8,sf$rsq,col=ifelse(sf$rsq==sf$rsq[which.max(sf$rsq)],"coral","steelblue"),cex=1.5,pch=16)
plot(sf$adjr2, xlab="Jumlah Variabel", ylab="Adj R2",type="l");points(1:8,sf$adjr2,col=ifelse(sf$adjr2==sf$adjr2[which.max(sf$adjr2)],"coral","steelblue"),cex=1.5,pch=16)
plot(sf$bic, xlab="Jumlah Variabel", ylab="BIC",type="l");points(1:8,sf$bic,col=ifelse(sf$bic==sf$bic[which.min(sf$bic)],"coral","steelblue"),cex=1.5,pch=16)
plot(sf$cp, xlab="Jumlah Variabel", ylab="CP",type="l");points(1:8,sf$cp,col=ifelse(sf$cp==sf$cp[which.min(sf$cp)],"coral","steelblue"),cex=1.5,pch=16)
Hasil dari seleksi peubah menunjukkan bahwa peubah Time dan semua peubah dummy Diet signifikan terhadap model sehingga peubah tersebut akan dimasukkan ke dalam model. Sementara itu, peubah dummy dari Chick banyak yang tidak signifikan terhadap model dan peubah tersebut berisi nomor identifikasi ayam (dengan kata lain nama2 ayam yang diubah menjadi numerik) sehingga peubah Chick tidak akan dimasukkan ke dalam model.
Regresi Klasik
modelklasik <- lm(weight~Time+Diet, data=dataayam)
summary(modelklasik)##
## Call:
## lm(formula = weight ~ Time + Diet, data = dataayam)
##
## Residuals:
## Min 1Q Median 3Q Max
## -136.851 -17.151 -2.595 15.033 141.816
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 10.9244 3.3607 3.251 0.00122 **
## Time 8.7505 0.2218 39.451 < 2e-16 ***
## Diet2 16.1661 4.0858 3.957 8.56e-05 ***
## Diet3 36.4994 4.0858 8.933 < 2e-16 ***
## Diet4 30.2335 4.1075 7.361 6.39e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 35.99 on 573 degrees of freedom
## Multiple R-squared: 0.7453, Adjusted R-squared: 0.7435
## F-statistic: 419.2 on 4 and 573 DF, p-value: < 2.2e-16Didapatkan model regresi klasik dengan R-squared sebesar 74.53%. Model ini memiliki p-value < 0.05 artinya cukup bukti untuk menyatakan bahwa minimal ada satu peubah penjelas yang berpengaruh signifikan terhadap peubah respon.
Berdasarkan hasil uji t, peubah penjelas yang berpengaruh yaitu peubah time, diet2, diet3, dan diet4. Dugaan model klasik yang diperoleh sebagai berikut.
\[Weight=10.9244+8.7505Time+16.1661Diet2+36.4994Diet3+30.2335Diet4\]
Interpretasi
Intersep: nilai rataan peubah respon (weight/berat ayam) ketika seluruh peubah penjelas bernilai 0 sebesar 10.9244 gram (signifikan)
Time: 8.7505, artinya jika waktu percobaan diet ayam bertambah satu hari, maka beratnya naik dengan rata-rata sebesar 8.7505 gram (signifikan)
Diet2: 16.1661, artinya ayam dengan diet 2 memiliki berat rata-rata 16.1661 gram lebih tinggi dari ayam dengan diet 1 (signifikan)
Diet3: 36.4994, artinya ayam dengan diet 3 memiliki berat rata-rata 36.4994 gram lebih tinggi dari ayam dengan diet 1 (signifikan)
Diet4: 30.2335, artinya ayam dengan diet 4 memiliki berat rata-rata 30.2335 gram lebih tinggi dari ayam dengan diet 1 (signifikan)
##cek multikol
library(car)## Loading required package: carData##
## Attaching package: 'car'## The following object is masked from 'package:lmridge':
##
## vifvif_values <- car::vif(modelklasik)
vif_values## GVIF Df GVIF^(1/(2*Df))
## Time 1.000832 1 1.000416
## Diet 1.000832 3 1.000139Dapat diketahui bahwa nilai VIF pada kedua peubah penjelas relatif rendah dan mendekati 1. Hal ini menunjukkan bahwa tidak ada multikolinieritas yang signifikan antara kedua peubah tersebut sehingga tidak perlu dilakukan penanganan. Akan tetapi, sebagai latihan akan dicoba metode penyusutan dengan regresi ridge dan lasso.
Regresi Ridge (glmnet)
Peubah
x <- data.matrix(dataayam[, c('Time', 'Diet')])
y <- dataayam$weightCV
cv.r<-cv.glmnet(x,y,alpha=0);plot(cv.r)
Best Model
best.lr<-cv.r$lambda.min
bestridge<-glmnet(x,y,alpha=0,lambda=best.lr);coef(bestridge)## 3 x 1 sparse Matrix of class "dgCMatrix"
## s0
## (Intercept) 10.656231
## Time 8.090758
## Diet 10.936109Dari hasil penyusutan menggunakan regresi Ridge dengan fungsi glmnet didapatkan dugaan best model sebagai berikut. \[Weight=10.656231+8.090758Time+10.936109Diet\] ## Interpretasi
Intersep: nilai rataan peubah respon (weight/berat ayam) ketika seluruh peubah penjelas bernilai 0 sebesar 10.656231 gram
Time: 8.090758, artinya jika waktu percobaan diet ayam bertambah satu hari, maka beratnya naik dengan rata-rata sebesar 8.090758 gram
Diet: 10.936109, artinya jika diet yang diberikan kepada ayam bertambah satu satuan, maka beratnya naik dengan rata-rata sebesar 10.936109 gram
##Fungsi R-Square
rsq<-function(bestmodel,bestlambda,x,y){
#y duga
y.duga <- predict(bestmodel, s = bestlambda, newx = x)
#JKG dan JKT
jkt <- sum((y - mean(y))^2)
jkg <- sum((y.duga- y)^2)
#find R-Squared
rsq <- 1 - jkg/jkt
return(rsq)
}R-Square Ridge
rsq(bestridge,best.lr,x,y)## [1] 0.7335174Model regresi Ridge tersebut memiliki nilai R-squared sebesar 73.35174%.
Regresi Ridge (lmridge)
library(lmridge)
lmr <- lmridge(dataayam$weight~dataayam$Time+dataayam$Diet, scaling="centered", data=dataayam)
plot(lmr)
Summary
summary(lmr)##
## Call:
## lmridge.default(formula = dataayam$weight ~ dataayam$Time + dataayam$Diet,
## data = dataayam, scaling = "centered")
##
##
## Coefficients: for Ridge parameter K= 0
## Estimate Estimate (Sc) StdErr (Sc) t-value (Sc) Pr(>|t|)
## Intercept 10.9244 10.9244 4.0644 2.6878 0.0074 **
## dataayam$Time 8.7505 8.7505 0.2216 39.4857 <2e-16 ***
## dataayam$Diet2 16.1661 16.1661 4.0823 3.9601 0.0001 ***
## dataayam$Diet3 36.4994 36.4994 4.0823 8.9409 <2e-16 ***
## dataayam$Diet4 30.2335 30.2335 4.1039 7.3670 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Ridge Summary
## R2 adj-R2 DF ridge F AIC BIC
## 0.74530 0.74400 4.00009 419.90889 4145.31467 7838.58705
## Ridge minimum MSE= 50.22119 at K= 0
## P-value for F-test ( 4.00009 , 574 ) = 7.250106e-169
## -------------------------------------------------------------------Didapatkan model regresi Ridge dengan R-squared sebesar 74.53%. Dari hasil penyusutan menggunakan regresi Ridge dengan fungsi lmridge didapatkan dugaan model sebagai berikut.
\[Weight=10.9244+8.7505Time+16.1661Diet2+36.4994Diet3+30.2335Diet4\]
Interpretasi
Intersep: nilai rataan peubah respon (weight/berat ayam) ketika seluruh peubah penjelas bernilai 0 sebesar 10.9244 gram (signifikan)
Time: 8.7505, artinya jika waktu percobaan diet ayam bertambah satu hari, maka beratnya naik dengan rata-rata sebesar 8.7505 gram (signifikan)
Diet2: 16.1661, artinya ayam dengan diet 2 memiliki berat rata-rata 16.1661 gram lebih tinggi dari ayam dengan diet 1 (signifikan)
Diet3: 36.4994, artinya ayam dengan diet 3 memiliki berat rata-rata 36.4994 gram lebih tinggi dari ayam dengan diet 1 (signifikan)
Diet4: 30.2335, artinya ayam dengan diet 4 memiliki berat rata-rata 30.2335 gram lebih tinggi dari ayam dengan diet 1 (signifikan)
Regresi Lasso (glmnet)
CV
cv.l<-cv.glmnet(x,y,alpha=1);plot(cv.l)
Best Model
best.ll<-cv.l$lambda.min
bestlasso<-glmnet(x,y,alpha=1,lambda=best.ll);coef(bestlasso)## 3 x 1 sparse Matrix of class "dgCMatrix"
## s0
## (Intercept) 2.562690
## Time 8.722327
## Diet 11.528596Dari hasil penyusutan menggunakan regresi Lasso dengan fungsi lmridge didapatkan dugaan model sebagai berikut.
\[Weight=2.562690+8.722327Time+11.528596Diet\] ## Interpretasi
Intersep: nilai rataan peubah respon (weight/berat ayam) ketika seluruh peubah penjelas bernilai 0 sebesar 2.562690 gram
Time: 8.722327, artinya jika waktu percobaan diet ayam bertambah satu hari, maka beratnya naik dengan rata-rata sebesar 8.722327 gram
Diet: 11.528596, artinya jika diet yang diberikan kepada ayam bertambah satu satuan, maka beratnya naik dengan rata-rata sebesar 11.528596 gram
R-Square Lasso
rsq(bestlasso,best.ll,x,y)## [1] 0.7378211Model regresi Lasso tersebut memiliki nilai R-squared sebesar 73.78211%.
Perbandingan Model
R-squared regresi klasik: 74.53%
R-squared regresi Ridge (glmnet): 73.35174%
R-squared regresi Ridge (lmridge): 74.53%
R-squared regresi Lasso (glmnet): 73.78211%
Model dari dataayam yang dilakukan penanganan dengan
metode penyusutan Ridge dan Lasso sejak awal tidak terindikasi adanya
multikolinearitas, sehingga nilai koefisien dan R-squared pada model
setelah dilakukan penanganan hanya mengalami sedikit perubahan apabila
dibandingkan dengan model sebelum penanganan (model klasik). Berdasarkan
nilai R-squared yang didapatkan, model dari metode penyusutan regresi
Ridge(lmridge) memiliki nilai yang paling tinggi yaitu sebesar 74.53%.
Nilai R-squared yang lebih tinggi menunjukkan model yang lebih baik
dalam menjelaskan keragaman dalam data.
Interpretasi Model Terbaik
Dugaan model terbaik dari regresi(lmridge) diperoleh sebagai berikut.
\[Weight=10.9244+8.7505Time+16.1661Diet2+36.4994Diet3+30.2335Diet4\]
Intersep: nilai rataan peubah respon (weight/berat ayam) ketika seluruh peubah penjelas bernilai 0 sebesar 10.9244 gram (signifikan)
Time: 8.7505, artinya jika waktu percobaan diet ayam bertambah satu hari, maka beratnya naik dengan rata-rata sebesar 8.7505 gram (signifikan)
Diet2: 16.1661, artinya ayam dengan diet 2 memiliki berat rata-rata 16.1661 gram lebih tinggi dari ayam dengan diet 1 (signifikan)
Diet3: 36.4994, artinya ayam dengan diet 3 memiliki berat rata-rata 36.4994 gram lebih tinggi dari ayam dengan diet 1 (signifikan)
Diet4: 30.2335, artinya ayam dengan diet 4 memiliki berat rata-rata 30.2335 gram lebih tinggi dari ayam dengan diet 1 (signifikan)