library(plm)
library(spdep)
library(readxl)
library(tidyverse)
library(dplyr)
library(tidyr)
library(DT)

Dataset

data <- read_excel("C:/Users/Dicky Girsang/Desktop/2023 Kolokium/data_penelitian.xlsx", sheet = "Sheet4")
data <- data %>%
  pivot_longer(cols = -Provinsi, 
               names_to = c(".value", "Tahun"), 
               names_pattern = "(\\w+)_(\\d{4})")

Inputasi nilai Y yang hilang di 2020

dengan rata-rata Y tahun sebelumnya di setiap provinsi

data_filled <- data %>%
  group_by(Provinsi) %>%
  mutate(Y = ifelse(is.na(Y), mean(Y, na.rm = TRUE), Y))
print(data_filled)

## # A tibble: 238 × 7
## # Groups:   Provinsi [34]
##    Provinsi       Tahun     Y    X1    X2    X3    X4
##    <chr>          <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
##  1 ACEH           2016   26.4  38.0  54.4  74.5  848.
##  2 ACEH           2017   23.5  49.8  62.1  70.6  873.
##  3 ACEH           2018   37.1  33.3  70.9  70.7  839.
##  4 ACEH           2019   22    49.8  64    70.0  819.
##  5 ACEH           2020   13.4  65.4  68.0  70.1  815.
##  6 ACEH           2021   12.1  66.7  69.3  74.4  834.
##  7 ACEH           2022    8    65.9  68.0  70.7  807.
##  8 SUMATERA UTARA 2016   24.4  33.5  64.0  69.7 1456.
##  9 SUMATERA UTARA 2017   16    50.9  63.6  67.2 1454.
## 10 SUMATERA UTARA 2018   32.4  25.7  71.9  68.3 1325.
## # ℹ 228 more rows

data <- data_filled
datatable(data, options = list(scrollX = TRUE, pageLength = 10))

summary(data)

##    Provinsi            Tahun                 Y                X1       
##  Length:238         Length:238         Min.   : 0.100   Min.   :21.27  
##  Class :character   Class :character   1st Qu.: 9.325   1st Qu.:50.70  
##  Mode  :character   Mode  :character   Median :18.650   Median :59.27  
##                                        Mean   :18.135   Mean   :58.61  
##                                        3rd Qu.:25.530   3rd Qu.:68.15  
##                                        Max.   :42.700   Max.   :81.46  
##        X2              X3              X4         
##  Min.   :19.63   Min.   :27.44   Min.   :  41.12  
##  1st Qu.:70.38   1st Qu.:55.90   1st Qu.: 195.70  
##  Median :77.52   Median :62.65   Median : 366.16  
##  Mean   :74.49   Mean   :62.25   Mean   : 785.72  
##  3rd Qu.:82.20   3rd Qu.:67.80   3rd Qu.: 838.18  
##  Max.   :92.78   Max.   :90.12   Max.   :4703.30

Deteksi Multikolinearitas

vif_model <- lm(Y ~ X1 + X2 + X3 + X4, data = data)
vif_values <- car::vif(vif_model)
print(vif_values)

##       X1       X2       X3       X4 
## 1.027088 1.287308 1.314296 1.044224

VIF<10, Tidak ada indikasi multikolinearitas yang signifikan dalam model regresi Anda, karena semua nilai VIF berada dalam rentang yang dapat diterima

Estimasi Parameter Regresi Data Panel

cem_model <- plm(Y ~ X1 + X2 + X3 + X4, data = data, model = "within", effect = "individual")
fem_model <- plm(Y ~ X1 + X2 + X3 + X4, data = data, model = "within", effect = "time")
rem_model <- plm(Y ~ X1 + X2 + X3 + X4, data = data, model = "random")
print(summary(cem_model))

## Oneway (individual) effect Within Model
## 
## Call:
## plm(formula = Y ~ X1 + X2 + X3 + X4, data = data, effect = "individual", 
##     model = "within")
## 
## Balanced Panel: n = 34, T = 7, N = 238
## 
## Residuals:
##      Min.   1st Qu.    Median   3rd Qu.      Max. 
## -14.91971  -3.04976  -0.13102   2.85868  17.69158 
## 
## Coefficients:
##      Estimate Std. Error  t-value Pr(>|t|)    
## X1 -0.5590876  0.0443771 -12.5986  < 2e-16 ***
## X2 -0.1245355  0.0589183  -2.1137  0.03578 *  
## X3  0.0093254  0.1067226   0.0874  0.93046    
## X4  0.0074912  0.0048807   1.5349  0.12640    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    16022
## Residual Sum of Squares: 7319.1
## R-Squared:      0.54319
## Adj. R-Squared: 0.45869
## F-statistic: 59.4557 on 4 and 200 DF, p-value: < 2.22e-16

print(summary(fem_model))

## Oneway (time) effect Within Model
## 
## Call:
## plm(formula = Y ~ X1 + X2 + X3 + X4, data = data, effect = "time", 
##     model = "within")
## 
## Balanced Panel: n = 34, T = 7, N = 238
## 
## Residuals:
##      Min.   1st Qu.    Median   3rd Qu.      Max. 
## -13.46855  -3.37131  -0.23679   3.33873  17.55714 
## 
## Coefficients:
##       Estimate  Std. Error t-value  Pr(>|t|)    
## X1  0.13257638  0.03934249  3.3698 0.0008838 ***
## X2 -0.06986366  0.03669190 -1.9041 0.0581667 .  
## X3 -0.27774167  0.03729299 -7.4476 1.967e-12 ***
## X4 -0.00095365  0.00032503 -2.9340 0.0036893 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    9789
## Residual Sum of Squares: 6598.5
## R-Squared:      0.32593
## Adj. R-Squared: 0.29623
## F-statistic: 27.4398 on 4 and 227 DF, p-value: < 2.22e-16

Terdapat perbedaan dalam hasil antara model dengan efek individu (CEM) dan efek waktu (FEM). - Model CEM memiliki R-squared yang lebih tinggi, menunjukkan bahwa model ini lebih baik dalam menjelaskan variabilitas dalam Y. Namun, dalam model FEM, efek waktu dan hubungan antara peubah X1 dan X3 terbukti signifikan. - Model CEM ini menghasilkan R-squared sebesar 0.54319, yang mengindikasikan bahwa sekitar 54.32% variabilitas dalam Y dapat dijelaskan oleh X1, X2, X3, dan X4 dalam model ini.

Pemilihan Model Regresi Data Panel

# Uji Chow untuk memilih antara CEM dan FEM
chow_test <- pFtest(cem_model, fem_model)
print(chow_test)

## 
##  F test for individual effects
## 
## data:  Y ~ X1 + X2 + X3 + X4
## F = -0.72924, df1 = 27, df2 = 200, p-value = 1
## alternative hypothesis: significant effects

# Uji Chow untuk memilih antara CEM dan FEM
chow_test <- pFtest(cem_model, fem_model)
print(chow_test)

## 
##  F test for individual effects
## 
## data:  Y ~ X1 + X2 + X3 + X4
## F = -0.72924, df1 = 27, df2 = 200, p-value = 1
## alternative hypothesis: significant effects

# Inisialisasi objek selected_model
selected_model <- NULL
selected_model_type <- ""

# Jika hasil uji Chow signifikan, maka pilih model dengan efek individu (CEM)
if (chow_test$p.value < 0.05) {
  selected_model <- cem_model
  selected_model_type <- "CEM"
} else {
  # Jika uji Chow tidak signifikan, lakukan uji Breusch Pagan
  bp_test <- plmtest(fem_model, type = c("bp"))
  print(bp_test)
  
  # Jika hasil uji Breusch Pagan signifikan, pilih model dengan efek waktu (FEM)
  if (bp_test$p.value < 0.05) {
    selected_model <- fem_model
    selected_model_type <- "FEM"
  } else {
    # Jika uji Breusch Pagan tidak signifikan, lakukan uji Hausman
    hausman_test <- phtest(fem_model, cem_model)
    print(hausman_test)
    
    # Jika hasil uji Hausman menunjukkan perbedaan signifikan antara model CEM dan FEM, pilih model CEM
    if (hausman_test$p.value < 0.05) {
      selected_model <- cem_model
      selected_model_type <- "CEM"
    } else {
      # Jika uji Hausman tidak signifikan, pilih model dengan efek waktu (FEM)
      selected_model <- fem_model
      selected_model_type <- "FEM"
    }
  }
}

## 
##  Lagrange Multiplier Test - (Breusch-Pagan)
## 
## data:  Y ~ X1 + X2 + X3 + X4
## chisq = 61.471, df = 1, p-value = 4.494e-15
## alternative hypothesis: significant effects

print(summary(selected_model))

## Oneway (time) effect Within Model
## 
## Call:
## plm(formula = Y ~ X1 + X2 + X3 + X4, data = data, effect = "time", 
##     model = "within")
## 
## Balanced Panel: n = 34, T = 7, N = 238
## 
## Residuals:
##      Min.   1st Qu.    Median   3rd Qu.      Max. 
## -13.46855  -3.37131  -0.23679   3.33873  17.55714 
## 
## Coefficients:
##       Estimate  Std. Error t-value  Pr(>|t|)    
## X1  0.13257638  0.03934249  3.3698 0.0008838 ***
## X2 -0.06986366  0.03669190 -1.9041 0.0581667 .  
## X3 -0.27774167  0.03729299 -7.4476 1.967e-12 ***
## X4 -0.00095365  0.00032503 -2.9340 0.0036893 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    9789
## Residual Sum of Squares: 6598.5
## R-Squared:      0.32593
## Adj. R-Squared: 0.29623
## F-statistic: 27.4398 on 4 and 227 DF, p-value: < 2.22e-16

print(paste("Model terpilih:", selected_model_type))

## [1] "Model terpilih: FEM"

Jadi, hasil ini mengindikasikan bahwa model FEM adalah yang terbaik untuk data panel ini. Peubah X1, X2, X3, dan X4 masing-masing memiliki pengaruh pada Y, dengan X3 memiliki pengaruh negatif yang sangat signifikan. R-squared yang lebih rendah menunjukkan bahwa sebagian besar variasi dalam Y tetap tidak dijelaskan oleh model ini, dan faktor lain mungkin perlu dipertimbangkan

Uji asumsi klasik pada data panel terpilih

Uji normalitas dengan Jarque-Bera

library(tseries)

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

jarque_bera_test <- jarque.bera.test(fem_model$residuals)
print(jarque_bera_test)

## 
##  Jarque Bera Test
## 
## data:  fem_model$residuals
## X-squared = 3.6229, df = 2, p-value = 0.1634

Uji heteroskedastisitas dengan Bptest

library(lmtest)

## Loading required package: zoo

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

bptest(fem_model)

## 
##  studentized Breusch-Pagan test
## 
## data:  fem_model
## BP = 21.357, df = 4, p-value = 0.000269

Uji autokorelasi dengan Durbin Watson

dw_test <- plm::pwartest(fem_model, alternative = "two.sided")
print(dw_test)

## 
##  Wooldridge's test for serial correlation in FE panels
## 
## data:  fem_model
## F = 138.13, df1 = 1, df2 = 202, p-value < 2.2e-16
## alternative hypothesis: serial correlation

Secara keseluruhan, model mungkin memiliki residu yang berdistribusi normal, tetapi menunjukkan heteroskedastisitas dan korelasi serial. Mungkin perlu menangani masalah-masalah ini saat menginterpretasikan hasil atau mempertimbangkan model alternatif untuk memperbaiki masalah tersebut

TRI-U Spatial Panel Analysis

Dicky Girsang

2023-10-30