Praktikum 8 Analisis Regresi

Package

library(readxl)

## Warning: package 'readxl' was built under R version 4.5.2

library(dplyr) # untuk manipulasi data

## Warning: package 'dplyr' was built under R version 4.5.3

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

library(car) # untuk multikolinearitas

## Warning: package 'car' was built under R version 4.5.3

## Loading required package: carData

## Warning: package 'carData' was built under R version 4.5.2

## 
## Attaching package: 'car'

## The following object is masked from 'package:dplyr':
## 
##     recode

library(olsrr) # untuk best subset

## Warning: package 'olsrr' was built under R version 4.5.2

## 
## Attaching package: 'olsrr'

## The following object is masked from 'package:datasets':
## 
##     rivers

library(MASS) # untuk forward, backward, dan stepwise

## Warning: package 'MASS' was built under R version 4.5.3

## 
## Attaching package: 'MASS'

## The following object is masked from 'package:olsrr':
## 
##     cement

## The following object is masked from 'package:dplyr':
## 
##     select

library(glmnet) # untuk ridge dan lasso

## Warning: package 'glmnet' was built under R version 4.5.3

## Loading required package: Matrix

## Loaded glmnet 5.0

select <- dplyr::select

Data

Data yang digunakan memiliki 38 amatan berupa Provinsi di Indonesia dengan 1 peubah respon yaitu Angka Harapan Hidup (AHH) dan 5 peubah penjelas, yaitu Rata-rata Lama Sekolah (RLS), PDRB per kapita ADHB (PDRB), Jumlah Penduduk (JP), Tingkat Pengangguran Terbuka (TPT), dan Persentase Penduduk Miskin (PPM).

data <- read_xlsx("C:/Users/hp/Downloads/dataset regresi.xlsx")
data

## # A tibble: 38 × 9
##    Provinsi Angka Harapan Hidup …¹ Persentase Angkatan …² Rata-rata Lama Sekol…³
##    <chr>                     <dbl>                  <dbl>                  <dbl>
##  1 ACEH                       141.                   65.1                   9.64
##  2 SUMATER…                   141.                   71.4                   9.93
##  3 SUMATER…                   141.                   70.3                   9.44
##  4 RIAU                       145.                   66.3                   9.43
##  5 JAMBI                      144.                   68.9                   8.9 
##  6 SUMATER…                   142.                   70.8                   8.57
##  7 BENGKULU                   140.                   71.7                   9.04
##  8 LAMPUNG                    143.                   70.4                   8.36
##  9 KEP. BA…                   143.                   68.9                   8.33
## 10 KEP. RI…                   143.                   69.2                  10.5 
## # ℹ 28 more rows
## # ℹ abbreviated names: ¹`Angka Harapan Hidup (tahun)`,
## #   ²`Persentase Angkatan Kerja Terhadap Penduduk Usia Kerja (TPAK) (%)`,
## #   ³`Rata-rata Lama Sekolah (tahun)`
## # ℹ 5 more variables: `PDRB per kapita ADHB (Ribu Rp)` <dbl>,
## #   `Jumlah Penduduk (Ribu Jiwa)` <dbl>,
## #   `Tingkat Pengangguran Terbuka (%)` <dbl>, …

data <- select(data, -Provinsi, -`Persentase Angkatan Kerja Terhadap Penduduk Usia Kerja (TPAK) (%)`, -`Rumah Tangga Yang Memiliki Akses Air Minum Layak (%)`) 
# menghapus peubah yang tidak akan digunakan

data <- rename(data,
               Y  = `Angka Harapan Hidup (tahun)`,
               X1 = `Rata-rata Lama Sekolah (tahun)`,
               X2 = `PDRB per kapita ADHB (Ribu Rp)`,
               X3 = `Jumlah Penduduk (Ribu Jiwa)`,
               X4 = `Tingkat Pengangguran Terbuka (%)`,
               X5 = `Persentase Penduduk Miskin (%)`)
# mengganti nama peubah agar lebih mudah digunakan

Uji F Parsial dan Sekuensial

F Sekuensial

# Model 1: tanpa prediktor
m1 <- lm( Y ~ 1, data = data)

# Model 2: tambah X1
m2 <- lm(Y ~ X1, data = data)

# Model 3: tambah X2
m3 <- lm(Y ~ X1 + X2, data = data)

# Model 4 : tambah X3
m4 <- lm(Y ~ X1 + X2 + X3, data = data)

# Model 5 : tambah X4
m5 <- lm(Y ~ X1 + X2 + X3 + X4, data = data)

# Model 6 : tambah X5
m6 <- lm(Y ~ X1 + X2 + X3 + X4 + X5, data = data)

anova(m1, m2, m3, m4, m5, m6)

## Analysis of Variance Table
## 
## Model 1: Y ~ 1
## Model 2: Y ~ X1
## Model 3: Y ~ X1 + X2
## Model 4: Y ~ X1 + X2 + X3
## Model 5: Y ~ X1 + X2 + X3 + X4
## Model 6: Y ~ X1 + X2 + X3 + X4 + X5
##   Res.Df     RSS Df Sum of Sq       F    Pr(>F)    
## 1     37 1050.18                                   
## 2     36  742.02  1   308.156 28.6166 7.187e-06 ***
## 3     35  718.07  1    23.952  2.2243  0.145651    
## 4     34  482.73  1   235.341 21.8546 5.100e-05 ***
## 5     33  447.43  1    35.299  3.2780  0.079612 .  
## 6     32  344.59  1   102.838  9.5499  0.004119 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Variabel X1 berpengaruh signifikan terhadap Y (p < 0.001). Setelah variabel X1 dimasukkan, variabel X2 tidak memberikan tambahan pengaruh yang signifikan (p = 0.146). Variabel X3 kemudian memberikan pengaruh tambahan yang signifikan setelah variabel X1 dan X2 dimasukkan (p < 0.001). Akan tetapi, variabel X4 tidak memberikan pengaruh tambahan yang signifikan (p = 0.080). Terakhir, variabel X5 memberikan pengaruh tambahan yang signifikan setelah variabel X1, X2, X3, dan X4 dimasukkan ke dalam model (p < 0.01).

Berdasarkan hasil tersebut, model yang cukup baik adalah Model 2 (Y ~ X1) dan Model 4 (Y ~ X1 + X2 + X3), karena keduanya menunjukkan peningkatan kemampuan model yang signifikan. Pemilihan model akhir dapat disesuaikan dengan tujuan analisis. Model 2 lebih sederhana dan mudah diinterpretasikan, sedangkan Model 4 memberikan informasi yang lebih lengkap dengan melibatkan lebih banyak variabel penjelas.

F Parsial

Akan diuji apakah variabel X3 signifikan setelah semua variabel lain dikontrol

# full model
full <- lm(Y ~ X1 + X2 + X3 + X4 + X5, data = data)

# reduced (tanpa X3)
red <- lm(Y ~ X1 + X2 + X4 + X5, data = data)

anova(red, full)

## Analysis of Variance Table
## 
## Model 1: Y ~ X1 + X2 + X4 + X5
## Model 2: Y ~ X1 + X2 + X3 + X4 + X5
##   Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
## 1     33 496.79                                  
## 2     32 344.59  1     152.2 14.133 0.0006849 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Variabel X3 berpengaruh signifikan terhadap Y setelah mengontrol X1, X2, X4, dan X5 (p < 0.05).

Pemodelan Awal

model1 <- lm(Y ~ ., data=data)
summary(model1)

## 
## Call:
## lm(formula = Y ~ ., data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.9260 -2.0078  0.1767  1.7943  9.4791 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.361e+02  6.879e+00  19.780  < 2e-16 ***
## X1           9.977e-01  7.748e-01   1.288 0.207062    
## X2           1.790e-05  9.770e-06   1.832 0.076320 .  
## X3           2.010e-04  5.346e-05   3.759 0.000685 ***
## X4          -6.226e-01  4.994e-01  -1.247 0.221545    
## X5          -3.599e-01  1.165e-01  -3.090 0.004119 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.282 on 32 degrees of freedom
## Multiple R-squared:  0.6719, Adjusted R-squared:  0.6206 
## F-statistic:  13.1 on 5 and 32 DF,  p-value: 5.592e-07

Persamaan yang terbentuk dari seluruh peubah tersebut menghasilkan nilai koefisien determinasi sebesar 62.06% dengan 2 parameter yang tidak signifikan yaitu peubah X1 dan X4. Selanjutnya akan diperiksa ada tidaknya multikolinearitas pada model

\(\hat{AHH} = 136.1+0.9977RLS+0.00001790PDRB+0.000201JP-0.6226TPT-0.3599PPM\)

Pemeriksaan Multikolinieritas

vif(model1)

##       X1       X2       X3       X4       X5 
## 3.109331 1.268760 1.260777 1.713793 2.121326

Hasil pemeriksaan nilai VIF menunjukkan semua peubah penjelas memiliki VIF<10 sehingga tidak ada multikolinearitas

Penyeleksian Peubah

Best Subset

bs <- ols_step_best_subset(model1)
bs

##    Best Subsets Regression   
## -----------------------------
## Model Index    Predictors
## -----------------------------
##      1         X5             
##      2         X3 X5          
##      3         X2 X3 X5       
##      4         X1 X2 X3 X5    
##      5         X1 X2 X3 X4 X5 
## -----------------------------
## 
##                                                     Subsets Regression Summary                                                     
## -----------------------------------------------------------------------------------------------------------------------------------
##                        Adj.        Pred                                                                                             
## Model    R-Square    R-Square    R-Square     C(p)        AIC         SBIC        SBC         MSEP        FPE       HSP       APC  
## -----------------------------------------------------------------------------------------------------------------------------------
##   1        0.4930      0.4790      0.4494    15.4412    214.1524    105.2803    219.0652    562.0274    15.5674    0.4225    0.5633 
##   2        0.5982      0.5752      0.5426     7.1852    207.3180     99.2954    213.8684    458.5432    13.0079    0.3546    0.4707 
##   3        0.6503      0.6194      0.5804     4.1063    204.0430     97.0909    212.2309    411.2041    11.9392    0.3273    0.4320 
##   4        0.6559      0.6142      0.5457     5.5543    205.4230     98.8734    215.2485    417.1915    12.3900    0.3422    0.4483 
##   5        0.6719      0.6206      0.5367     6.0000    205.6206     99.9610    217.0837    410.7003    12.4687    0.3474    0.4512 
## -----------------------------------------------------------------------------------------------------------------------------------
## AIC: Akaike Information Criteria 
##  SBIC: Sawa's Bayesian Information Criteria 
##  SBC: Schwarz Bayesian Criteria 
##  MSEP: Estimated error of prediction, assuming multivariate normality 
##  FPE: Final Prediction Error 
##  HSP: Hocking's Sp 
##  APC: Amemiya Prediction Criteria

Backward, Forward, dan Stepwise

null<-lm(Y ~ 1, data=data) 
full<-lm(Y ~ ., data=data)
 
step(full, scope=list(lower=null, upper=full),data=data, direction='backward', trace=0)

## 
## Call:
## lm(formula = Y ~ X2 + X3 + X5, data = data)
## 
## Coefficients:
## (Intercept)           X2           X3           X5  
##   1.431e+02    2.064e-05    1.704e-04   -4.415e-01

step(null, scope=list(lower=null, upper=full),data=data, direction='forward', trace=0)

## 
## Call:
## lm(formula = Y ~ X5 + X3 + X2, data = data)
## 
## Coefficients:
## (Intercept)           X5           X3           X2  
##   1.431e+02   -4.415e-01    1.704e-04    2.064e-05

step(null, scope=list(lower=null, upper=full),data=data, direction='both', trace=0)

## 
## Call:
## lm(formula = Y ~ X5 + X3 + X2, data = data)
## 
## Coefficients:
## (Intercept)           X5           X3           X2  
##   1.431e+02   -4.415e-01    1.704e-04    2.064e-05

Hasil penyeleksian melalui best subset, backward, forward, dan stepwise menghasilkan model yang sama, yaitu dengan 3 peubah penjelas (PDRB, JP, dan PPM)

model2 <- lm(Y ~ X2 + X3 + X5, data=data)
summary(model2)

## 
## Call:
## lm(formula = Y ~ X2 + X3 + X5, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.6192 -1.9644  0.2728  1.7115 11.0680 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.431e+02  1.624e+00  88.099  < 2e-16 ***
## X2           2.064e-05  9.171e-06   2.250  0.03102 *  
## X3           1.704e-04  4.910e-05   3.469  0.00144 ** 
## X5          -4.415e-01  8.609e-02  -5.129 1.17e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.287 on 34 degrees of freedom
## Multiple R-squared:  0.6503, Adjusted R-squared:  0.6194 
## F-statistic: 21.07 on 3 and 34 DF,  p-value: 6.813e-08

AIC(model2); BIC(model2)

## [1] 204.043

## [1] 212.2309

Model yang dihasilkan memiliki 3 peubah penjelas yang signifikan pada taraf nyata 5% dengan nilai adjusted r square sebesar 61.94% yang lebih kecil daripada model dengan semua peubah. Nilai AIC pada model ini adalah 204.043 dan BIC sebesar 212.2309

\[\hat{AHH} = 143.1+0.00002064PDRB+0.0001704JP-0.4415PPM\]

Jika pada metode seleksi peubah yang digunakan menghasilkan model yang berbeda, maka bisa dibandingkan nilai MSE, adjusted R-square, AIC, BIC, dan kriteria lainnya.

Shrinkage Methods

Ridge

x <- as.matrix(data[, -1])
y <- data$Y

rid <- cv.glmnet(x, y, alpha = 0, nfolds = 5)  # Ridge regression
# secara default, nfolds = 10 sehingga akan buruk ketika n<30
rid

## 
## Call:  cv.glmnet(x = x, y = y, nfolds = 5, alpha = 0) 
## 
## Measure: Mean-Squared Error 
## 
##     Lambda Index Measure    SE Nonzero
## min  1.237    87   11.47 3.312       5
## 1se  7.953    67   14.43 3.967       5

coef(rid, s = "lambda.min")

## 6 x 1 sparse Matrix of class "dgCMatrix"
##                lambda.min
## (Intercept)  1.358409e+02
## X1           8.386232e-01
## X2           1.514394e-05
## X3           1.572976e-04
## X4          -2.644206e-01
## X5          -3.039155e-01

Model yang terbentuk menggunakan semua peubah penjelas, pemodelan ridge akan lebih terlihat berguna jika digunakan peubah yang terdapat multikolinearitas sehingga tidak perlu disisihkan ketika ada yang memiliki VIF > 10

\(\hat{AHH} = 135.8123+0.8299393RLS+0.00001492887PDRB+0.0001542184JP-0.2428861TPT-0.2992546PPM\)

Lasso

lass <- cv.glmnet(x, y, alpha=1, nfolds = 5)
lass

## 
## Call:  cv.glmnet(x = x, y = y, nfolds = 5, alpha = 1) 
## 
## Measure: Mean-Squared Error 
## 
##     Lambda Index Measure    SE Nonzero
## min 0.3606    26   12.29 2.457       4
## 1se 1.0035    15   14.39 2.535       4

coef(lass,s="lambda.min")

## 6 x 1 sparse Matrix of class "dgCMatrix"
##                lambda.min
## (Intercept)  1.400811e+02
## X1           3.410312e-01
## X2           1.356544e-05
## X3           1.439335e-04
## X4           .           
## X5          -3.727143e-01

Lasso menghasilkan pemodelan regresi dengan penyusutan/penyeleksian peubah yang memiliki koefisien yang sangat mendekati atau bahkan tepat nol

\[\hat{AHH}=139.0028+0.4912041RLS+0.00001521624PDRB+0.0001606378JP−1.171122TPT−0.3726701PPM\]