library(readxl)   
## Warning: package 'readxl' was built under R version 4.4.2
library(ggplot2)   
## Warning: package 'ggplot2' was built under R version 4.4.3
library(car)       
## Warning: package 'car' was built under R version 4.4.3
## Loading required package: carData
## Warning: package 'carData' was built under R version 4.4.3
library(lmtest)
## Warning: package 'lmtest' was built under R version 4.4.3
## Loading required package: zoo
## Warning: package 'zoo' was built under R version 4.4.3
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(olsrr)
## Warning: package 'olsrr' was built under R version 4.4.3
## 
## Attaching package: 'olsrr'
## The following object is masked from 'package:datasets':
## 
##     rivers
library(MASS)
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:olsrr':
## 
##     cement
library(nlme)
## Warning: package 'nlme' was built under R version 4.4.3
data <- read_excel("C:/Users/Resea/Documents/Data Anreg Ulya Fatimah.xlsx")
data
## # A tibble: 34 × 7
##    Provinsi                PDRB    IDN    ILN   TPT  PLKI   PPK
##    <chr>                  <dbl>  <dbl>  <dbl> <dbl> <dbl> <dbl>
##  1 ACEH                  41408.  8883.  249.   6.03  59.9 10334
##  2 SUMATERA UTARA        68306. 21574  1181.   5.89  58.5 11049
##  3 SUMATERA BARAT        54327.  4488.  121.   5.94  64.2 11380
##  4 RIAU                 154522. 48243. 2042.   4.23  52.1 11448
##  5 JAMBI                 79850.  8939    45.1  4.53  59.7 11160
##  6 SUMATERA SELATAN      71958. 25602. 1479.   4.11  63.0 11472
##  7 BENGKULU              46300.  7219.   76.1  3.42  67.8 11172
##  8 LAMPUNG               48191.  7626.  221.   4.23  70.7 10769
##  9 KEP. BANGKA BELITUNG  67813.  7961.   72.5  4.56  51.1 13589
## 10 KEP. RIAU            154065.  8857.  764.   6.8   33.7 14998
## # ℹ 24 more rows
str(data)
## tibble [34 × 7] (S3: tbl_df/tbl/data.frame)
##  $ Provinsi: chr [1:34] "ACEH" "SUMATERA UTARA" "SUMATERA BARAT" "RIAU" ...
##  $ PDRB    : num [1:34] 41408 68306 54327 154522 79850 ...
##  $ IDN     : num [1:34] 8883 21574 4488 48243 8939 ...
##  $ ILN     : num [1:34] 248.6 1181.3 120.7 2042.3 45.1 ...
##  $ TPT     : num [1:34] 6.03 5.89 5.94 4.23 4.53 4.11 3.42 4.23 4.56 6.8 ...
##  $ PLKI    : num [1:34] 59.9 58.5 64.2 52.1 59.7 ...
##  $ PPK     : num [1:34] 10334 11049 11380 11448 11160 ...
summary(data)
##    Provinsi              PDRB             IDN             ILN        
##  Length:34          Min.   : 23078   Min.   : 1174   Min.   :   8.3  
##  Class :character   1st Qu.: 48233   1st Qu.: 5487   1st Qu.: 109.9  
##  Mode  :character   Median : 64110   Median : 8489   Median : 458.4  
##                     Mean   : 81942   Mean   :19779   Mean   :1444.7  
##                     3rd Qu.: 77359   3rd Qu.:24595   3rd Qu.:1442.1  
##                     Max.   :322619   Max.   :95202   Max.   :8283.7  
##       TPT             PLKI            PPK       
##  Min.   :2.270   Min.   :33.67   Min.   : 7562  
##  1st Qu.:3.487   1st Qu.:52.71   1st Qu.:10125  
##  Median :4.320   Median :59.80   Median :11276  
##  Mean   :4.614   Mean   :59.35   Mean   :11470  
##  3rd Qu.:5.763   3rd Qu.:65.11   3rd Qu.:12285  
##  Max.   :7.520   Max.   :84.43   Max.   :19373
# Eksplorasi Data
# Melihat korelasi antar variabel
cor(data[, sapply(data, is.numeric)], use = "complete.obs")
##            PDRB        IDN        ILN        TPT       PLKI        PPK
## PDRB  1.0000000  0.4844513  0.2804425  0.2465797 -0.6425002  0.5220616
## IDN   0.4844513  1.0000000  0.6536220  0.4551561 -0.4124790  0.5184135
## ILN   0.2804425  0.6536220  1.0000000  0.3228454 -0.2135567  0.1858715
## TPT   0.2465797  0.4551561  0.3228454  1.0000000 -0.5772479  0.3221756
## PLKI -0.6425002 -0.4124790 -0.2135567 -0.5772479  1.0000000 -0.7610313
## PPK   0.5220616  0.5184135  0.1858715  0.3221756 -0.7610313  1.0000000
# Visualisasi pasangan variabel
data_num <- data[sapply(data, is.numeric)]
pairs(data_num)

# Membuat Model Regresi Linear Berganda
model_awal <- lm(PDRB ~ IDN + ILN + TPT + PLKI + PPK, data = data)
summary(model_awal)
## 
## Call:
## lm(formula = PDRB ~ IDN + ILN + TPT + PLKI + PPK, data = data)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -73691 -30306  -1828  11598 124120 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  5.073e+05  1.652e+05   3.071  0.00470 **
## IDN          1.082e+00  5.277e-01   2.051  0.04977 * 
## ILN         -7.200e-01  4.975e+00  -0.145  0.88597   
## TPT         -1.571e+04  7.660e+03  -2.051  0.04970 * 
## PLKI        -5.071e+03  1.432e+03  -3.540  0.00142 **
## PPK         -6.298e+00  6.339e+00  -0.993  0.32899   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 45520 on 28 degrees of freedom
## Multiple R-squared:  0.5413, Adjusted R-squared:  0.4593 
## F-statistic: 6.607 on 5 and 28 DF,  p-value: 0.0003543

Maka, model awalnya adalah PDRB = 5.073+1.082IDN-7.200ILN-1.57TPT-5.071PLKI–6.298PPK.

# Cek Asumsi

# Normalitas Residual
qqnorm(residuals(model_awal))
qqline(residuals(model_awal))

shapiro.test(residuals(model_awal))  # Uji Shapiro-Wilk
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(model_awal)
## W = 0.94024, p-value = 0.06269

Residual berdistribusi normal (p-value > 0.05).

# Multikolinearitas
vif(model_awal)
##      IDN      ILN      TPT     PLKI      PPK 
## 2.717815 1.877402 1.881736 3.540932 3.284688

Tidak ada multikolinearitas.

# Homoskedastisitas
bptest(model_awal) # Uji Breusch-Pagan
## 
##  studentized Breusch-Pagan test
## 
## data:  model_awal
## BP = 12.656, df = 5, p-value = 0.02683

Homoskedastisitas tidak terpenuhi (seharusnya p-value > 0.05).

# Independence of Residuals
dwtest(model_awal)  # Uji Durbin-Watson
## 
##  Durbin-Watson test
## 
## data:  model_awal
## DW = 1.4211, p-value = 0.03338
## alternative hypothesis: true autocorrelation is greater than 0

Adanya autokorelasi. Seharusnya nilai DW mendekati 2 → Residual independen.

library(nlme)

# Model Generalized Least Squares (GLS) untuk koreksi autokorelasi AR(1) dan penanganan heteroskedastisitas
model_gls_hetero <- gls(PDRB ~ IDN + ILN + TPT + PLKI + PPK,
                        data = data,
                        correlation = corAR1(form = ~ 1),
                        weights = varIdent(form = ~ 1))

summary(model_gls_hetero)
## Generalized least squares fit by REML
##   Model: PDRB ~ IDN + ILN + TPT + PLKI + PPK 
##   Data: data 
##        AIC      BIC    logLik
##   767.1937 777.8513 -375.5969
## 
## Correlation Structure: AR(1)
##  Formula: ~1 
##  Parameter estimate(s):
##       Phi 
## 0.4585625 
## 
## Coefficients:
##                Value Std.Error   t-value p-value
## (Intercept) 384459.0 145944.38  2.634284  0.0136
## IDN              1.2      0.44  2.772778  0.0098
## ILN             -1.1      4.29 -0.253088  0.8020
## TPT         -18304.5   6467.68 -2.830146  0.0085
## PLKI         -4128.1   1251.18 -3.299342  0.0026
## PPK              0.5      5.78  0.093855  0.9259
## 
##  Correlation: 
##      (Intr) IDN    ILN    TPT    PLKI  
## IDN   0.379                            
## ILN  -0.318 -0.642                     
## TPT  -0.644 -0.352  0.203              
## PLKI -0.931 -0.261  0.210  0.550       
## PPK  -0.870 -0.457  0.367  0.378  0.692
## 
## Standardized residuals:
##         Min          Q1         Med          Q3         Max 
## -1.53114016 -0.78570886 -0.06965709  0.45106898  1.84263912 
## 
## Residual standard error: 48184.54 
## Degrees of freedom: 34 total; 28 residual
# Ambil residual dari model_gls
resid_gls_hetero <- residuals(model_gls_hetero, type = "normalized")
acf(resid_gls_hetero)

Garis biru putus-putus itu adalah batas signifikansi. Kalau batang ACF jatuh di antara batas biru, artinya tidak signifikan → tidak ada autokorelasi di lag tersebut. Kalau batang keluar dari batas biru, baru dianggap ada autokorelasi.

Dapat dilihat hanya di lag ke-0 (nilai 1, ini normal karena residual dengan dirinya sendiri = 1). Untuk lag 1, 2, 3, dst, semua batangnya kecil dan masih di dalam batas biru. Kesimpulannya tidak ada pola autokorelasi yang kuat di residual.

plot(fitted(model_gls_hetero), residuals(model_gls_hetero, type = "normalized"))
abline(h = 0, col = "red")

Titik-titik residual tersebar cukup acak di sekitar garis merah (y=0). Tidak terlihat pola melebar (seperti kipas), menyempit, atau berbentuk sistematis tertentu. Variasi residual di sepanjang fitted values terlihat cukup konsisten. Kesimpulannya heteroskedastisitas sudah tertangani dengan cukup baik.

# Seleksi Model
# stepwise pada model
model_terbaik <- ols_step_both_p(model_awal)
summary(model_terbaik$model)
## 
## Call:
## lm(formula = paste(response, "~", paste(preds, collapse = " + ")), 
##     data = l)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -67804 -33733  -2616  15388 108435 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.626e+05  7.753e+04   4.677 5.79e-05 ***
## PLKI        -4.004e+03  9.365e+02  -4.275 0.000178 ***
## IDN          8.561e-01  3.612e-01   2.370 0.024397 *  
## TPT         -1.299e+04  7.028e+03  -1.849 0.074380 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 44750 on 30 degrees of freedom
## Multiple R-squared:  0.5249, Adjusted R-squared:  0.4774 
## F-statistic: 11.05 on 3 and 30 DF,  p-value: 4.722e-05
# Pakai stepAIC untuk BIC (menggunakan k = log(n))
n <- nrow(data)
model_bic <- stepAIC(model_awal, direction = "both", k = log(n), trace = FALSE)
summary(model_bic)
## 
## Call:
## lm(formula = PDRB ~ IDN + TPT + PLKI, data = data)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -67804 -33733  -2616  15388 108435 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.626e+05  7.753e+04   4.677 5.79e-05 ***
## IDN          8.561e-01  3.612e-01   2.370 0.024397 *  
## TPT         -1.299e+04  7.028e+03  -1.849 0.074380 .  
## PLKI        -4.004e+03  9.365e+02  -4.275 0.000178 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 44750 on 30 degrees of freedom
## Multiple R-squared:  0.5249, Adjusted R-squared:  0.4774 
## F-statistic: 11.05 on 3 and 30 DF,  p-value: 4.722e-05

Model stepwise dan model BIC hasilnya sama.

# Cek Asumsi
# Normalitas Residual
qqnorm(residuals(model_bic))
qqline(residuals(model_bic))

shapiro.test(residuals(model_bic))  # Uji Shapiro-Wilk
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(model_bic)
## W = 0.94039, p-value = 0.06332

Residual berdistribusi normal (p-value > 0.05).

# Multikolinearitas
vif(model_bic)
##      IDN      TPT     PLKI 
## 1.317167 1.639308 1.566162

Tidak ada multikolinearitas.

# Homoskedastisitas
bptest(model_bic) # Uji Breusch-Pagan
## 
##  studentized Breusch-Pagan test
## 
## data:  model_bic
## BP = 8.9416, df = 3, p-value = 0.03008

Homoskedastisitas tidak terpenuhi (p-value < 0.05).

# Independence of Residuals
dwtest(model_bic)  # Uji Durbin-Watson
## 
##  Durbin-Watson test
## 
## data:  model_bic
## DW = 1.2633, p-value = 0.009609
## alternative hypothesis: true autocorrelation is greater than 0

Nilai DW mendekati 2 → Residual independen. Tetapi terjadi autokorelasi.

library(nlme)

# Menggunakan data asli dari model lm
# Asumsi: model_bic adalah objek lm yang sudah dibuat sebelumnya
data_for_gls <- model_bic$model  # Ini mengambil data yang digunakan dalam model lm

# Model GLS dengan koreksi AR(1) dan heteroskedastisitas
model_gls_bic <- gls(PDRB ~ IDN + TPT + PLKI,
                    data = data_for_gls,
                    correlation = corAR1(form = ~ 1),
                    weights = varIdent(form = ~ 1),
                    method = "ML")  # ML untuk komparasi model

summary(model_gls_bic)
## Generalized least squares fit by maximum likelihood
##   Model: PDRB ~ IDN + TPT + PLKI 
##   Data: data_for_gls 
##        AIC     BIC    logLik
##   827.4198 836.578 -407.7099
## 
## Correlation Structure: AR(1)
##  Formula: ~1 
##  Parameter estimate(s):
##       Phi 
## 0.4270663 
## 
## Coefficients:
##                Value Std.Error   t-value p-value
## (Intercept) 393888.9  69887.51  5.636041  0.0000
## IDN              1.2      0.32  3.685315  0.0009
## TPT         -18169.8   5851.54 -3.105126  0.0041
## PLKI         -4200.6    875.75 -4.796524  0.0000
## 
##  Correlation: 
##      (Intr) IDN    TPT   
## IDN  -0.056              
## TPT  -0.699 -0.218       
## PLKI -0.927  0.071  0.442
## 
## Standardized residuals:
##         Min          Q1         Med          Q3         Max 
## -1.78668034 -0.85842603 -0.07622174  0.43561528  2.05910848 
## 
## Residual standard error: 43055.74 
## Degrees of freedom: 34 total; 30 residual
BIC(model_bic)
## [1] 838.0654
BIC(model_gls_bic)
## [1] 836.578

Penurunan BIC sebesar 1.4874 menunjukkan model GLS lebih baik. Aturan praktis: Selisih BIC > 2 dianggap signifikan, jadi perbaikan ini cukup baik meskipun tidak dramatis.

# Interpretasi Model
# Coefficients interpretation
coefficients(model_gls_bic$model)
## corStruct 
##  0.912606

Model GLS BIC atau stepwsise (karena hasilnya sama) dengan koreksi autokorelasi AR(1) dan heteroskedastisitas menunjukkan perbaikan dibanding model OLS, dengan penurunan BIC dari 838.07 menjadi 836.58. Variabel construct memiliki pengaruh positif signifikan terhadap PDRB (β=0.913, p<0.05). Diagnostic plot menunjukkan residual sudah memenuhi asumsi normalitas dan homoskedastisitas setelah koreksi.

Maka model terbaik adalah model GLS BIC