Penelitian ITB ( Institut Teknologi Bandung ) Metode Regresi Data Panel dengan R

By . Posted on

Pada Project kali ini, Whitecyber Team mendapatkan Project untuk membantu temen-temen kita di ITB – Institut Teknologi Bandung. Nah pada kali ini saya akan membahas mengenai penelitian tersebut dengan waktu kerja prakteknya. Metode penelitian yang kami gunakan adalah regresi data panel. Yuk, kita mengenal sedikit mengenai regresi data panel.

Regresi Data Panel

Regresi Data Panel adalah gabungan antara data cross section dan data time series, dimana unit cross section yang sama diukur pada waktu yang berbeda. Maka dengan kata lain, data panel merupakan data dari beberapa individu sama yang diamati dalam kurun waktu tertentu. Jika kita memiliki T periode waktu (t =1, 2, …, T) dan N jumlah individu (i = 1, 2, …, N), maka dengan data panel akan memiliki total unit observasi sebanyak N, T. Jika jumlah unit waktu sama untuk setiap individu, maka data disebut balanced panel. Jika sebaliknya, yakni jumlah unit waktu berbeda untuk setiap individu, maka disebut unbalanced panel.

Model Regresi Data Panel

Terdapat tiga pendekatan dalam perhitungan model regresi data panel, yaitu Common Effect Models, Fixed Effect Models, dan Random Effect Models.

  1. Common Effect Models, merupakan pendekatan yang paling sederhana untuk mengestimasi model regresi data panel. Estimasi untuk model ini dapat dilakukan dengan menggunakan metode OLS (Ordinary Least Square).
  2. Fixed Effect Models (FEM), FEM diasumsikan bahwa koefisien slope bernilai konstan tapi intercept bersifat tidak konstan. Metode yang dapat dilakukan untuk estimasi model dalam FEM, yaitu metode Least Square Dummy Variable atau yang sering disebut LSDV. Dalam metode LSDV, estimasi dilakukan dengan memasukkan variabel dummy yang digunakan untuk menjelaskan nilai intersep yang berbeda-beda akibat perbedaan nilai unit.
  3. Random Effect Models (REM), pada FEM atau model efek tetap, perbedaan karakteristik unit dan periode waktu diakomodasikan pada intercept, sehingga intercept dapat berubah antar waktu. Sementara untuk REM atau model efek random, perbedaan karakteristik unit dan periode waktu diakomodasikan pada error atau residual dari model. Dikarenakan ada dua komponen yang berkontribusi pada pembentukan error, yakni unit dan periode waktu, maka random error dalam REM perlu diurai menjadi error gabungan dan error untuk periode waktu

Pemilihan Model Regresi Data Panel

  1. Uji Chow, pengujian yang dilakukan untuk mengetahui apakah metode fixed effect lebih baik digunakan daripada menggunakan metode common effect.
  2. Uji Hausman, pengujian yang dilakukan untuk mengetahui apakah metode random effect merupakan metode yang lebih baik untuk digunakan bila dibandingkan dengan metode fixed effect.
  3. Uji Breusch Pagan, pengujian dilakukan untuk mngetahui apakah terdapat efek individu/waktu (atau keduanya) di dalam panel data.

Nah, sekarang kita langsung menerapkannya menggunakan software R yuk.

Analisis Deskriptif

Mari kita mulai dengan membuat analisa Deskriptif dari tiap-tiap kategori.

  1. GR Gini Rasio Kabupaten / Kota
  2. IPM Indeks Pembangunan Manusia
  3. IP Indeks Pendidikan
  4. IK Indeks Kesehatan
  5. AHH Angka Harapan Hidup
  6. RLS Rata-Rata Lama Sekolah (Tahun)
  7. TPAK Tingkat Partisipasi Angkatan Kerja (Persen)
  8. TPT Tingkat Pengangguran Terbuka Kabupaten / Kota (Persen)
  9. PM Persentase Penduduk Miskin (Persen)
  10. KP Kepadatan Penduduk (Jiwa/km2)
  11. PDRB PDRB per Kapita Atas Dasar Harga Konstan Menurut Kabupaten / Kota (Ribu Rupiah)
  12. PP Pengeluran Pemerintah (Rupiah)
  13. SD Angka Partisipasi Murni SD
  14. SMP Angka Partisipasi Murni SMP
  15. SMA Angka Partisipasi Murni SMA
  16. PK Pengeluaran Per Kapita (Ribu Rupiah / Orang / Tahun)
  17. UMK upah minimum

.

Kita mulai dari membangun Script dan Analisa Data, tampilannya adalah sebagai berikut :

.

SOURCE CODE 1

library(tidyverse)
library(plm)
library(car)
library(gplots)
library(tseries)
library(lmtest)
library(splm)
library(DescTools)
library(rgdal)
library(spdep)
library(spData)
library(raster)
library(sp)
library(nortest)
library(MASS)
library(readxl)
library(ggplot2)
library(grid)
library(gridExtra)
library(RColorBrewer)
library(randtests)
library(PerformanceAnalytics)

#==== Import Data ====
setwd("/Users/ainnamuslidasafira/Documents/data baru 2022")
dt.eksplorasi <- read_xlsx(path = "data baru.xlsx", sheet = "Sheet2 (2)")
colnames(dt.eksplorasi)

KOTKAB <- data.frame(KOTKAB=unlist(dt.eksplorasi[,1]),ntimes=c(rep(5,27)))
KOTKAB <- as.data.frame(lapply(KOTKAB, rep, KOTKAB$ntimes))
KOTKAB$ntimes <- NULL
Tahun <- data.frame(Tahun = rep(c(2018:2022),27))
Long <- data.frame(Long=unlist(dt.eksplorasi[,172]),ntimes=c(rep(5,27)))
Long <- as.data.frame(lapply(Long, rep, Long$ntimes))
Long$ntimes <- NULL
Lat <- data.frame(Lat=unlist(dt.eksplorasi[,173]),ntimes=c(rep(5,27)))
Lat <- as.data.frame(lapply(Lat, rep, Lat$ntimes))
Lat$ntimes <- NULL

data <- data.frame(KOTKAB,Tahun,Long,Lat,
                   GR = unlist(data.frame(t(dt.eksplorasi[,2:6]))),
                   IP= unlist(data.frame(t(dt.eksplorasi[,7:11]))),
                   IK= unlist(data.frame(t(dt.eksplorasi[,12:16]))),
                   TPAK= unlist(data.frame(t(dt.eksplorasi[,17:21]))),
                   TPT= unlist(data.frame(t(dt.eksplorasi[,22:26]))),
                   PM = unlist(data.frame(t(dt.eksplorasi[,27:31]))),
                   KP = unlist(data.frame(t(dt.eksplorasi[,32:36]))),
                   PDRB = unlist(data.frame(t(dt.eksplorasi[,37:41]))),
                   PP = unlist(data.frame(t(dt.eksplorasi[,42:46]))),
                   SD = unlist(data.frame(t(dt.eksplorasi[,47:51]))),
                   SMP = unlist(data.frame(t(dt.eksplorasi[,52:56]))),
                   SMA = unlist(data.frame(t(dt.eksplorasi[,57:61]))),
                   gr = unlist(data.frame(t(dt.eksplorasi[,62:66]))),
                   ip= unlist(data.frame(t(dt.eksplorasi[,67:71]))),
                   ik= unlist(data.frame(t(dt.eksplorasi[,72:76]))),
                   tpak= unlist(data.frame(t(dt.eksplorasi[,77:81]))),
                   tpt= unlist(data.frame(t(dt.eksplorasi[,82:86]))),
                   pm = unlist(data.frame(t(dt.eksplorasi[,87:91]))),
                   kp = unlist(data.frame(t(dt.eksplorasi[,92:96]))),
                   pdrb = unlist(data.frame(t(dt.eksplorasi[,97:101]))),
                   pp = unlist(data.frame(t(dt.eksplorasi[,102:106]))),
                   sd = unlist(data.frame(t(dt.eksplorasi[,107:111]))),
                   smp = unlist(data.frame(t(dt.eksplorasi[,112:116]))),
                   sma = unlist(data.frame(t(dt.eksplorasi[,117:121]))),
                   AHH = unlist(data.frame(t(dt.eksplorasi[,122:126]))),
                   RLS = unlist(data.frame(t(dt.eksplorasi[,127:131]))),
                   PK = unlist(data.frame(t(dt.eksplorasi[,132:136]))),
                   UMK = unlist(data.frame(t(dt.eksplorasi[,137:141]))),
                   ahh = unlist(data.frame(t(dt.eksplorasi[,142:146]))),
                   rls = unlist(data.frame(t(dt.eksplorasi[,147:151]))),
                   pk = unlist(data.frame(t(dt.eksplorasi[,152:156]))),
                   umk = unlist(data.frame(t(dt.eksplorasi[,157:161]))),
                   IPM = unlist(data.frame(t(dt.eksplorasi[,162:166]))),
                   ipm = unlist(data.frame(t(dt.eksplorasi[,167:171]))),
                   row.names = NULL)
View(data)


#=====Grafik Pergerakan Rasio Gini=====
rasiogini<-read_excel(path = "data baru.xlsx", sheet="Sheet6")
Tahun<-rasiogini$Tahun
GR<-rasiogini$Y
KabKota<-rasiogini$KOTKAB
ggplot(rasiogini,aes(Tahun, GR, colour=KabKota)) + geom_line()

#=====Plot Korelasi Sebelum Transformasi=====
datakorelasi<-data.frame(data$GR,data$IPM,data$TPT,data$PM,data$PP,data$KP,data$PDRB)
datakorelasi
chart.Correlation(datakorelasi, histogram=TRUE, pch=1)

#=======cek vif sebelum transformasi===========
#Tahun 2018
as.data.frame(vif(lm(GR18~IPM18+PP18+PM18+PDRB18+KP18+TPT18, data=dt.eksplorasi))) 
#Tahun 2019
as.data.frame(vif(lm(GR19~IPM19+PP19+PM19+PDRB19+KP19+TPT19, data=dt.eksplorasi)))
#Tahun 2020
as.data.frame(vif(lm(GR20~IPM20+PP20+PM20+PDRB20+KP20+TPT20, data=dt.eksplorasi)))
#Tahun 2021
as.data.frame(vif(lm(GR21~IPM21+PP21+PM21+PDRB21+KP21+TPT21, data=dt.eksplorasi)))
#Tahun 2022
as.data.frame(vif(lm(GR22~IPM22+PP22+PM22+PDRB22+KP22+TPT22, data=dt.eksplorasi)))
#Tahun 2018-2022
as.data.frame(vif(lm(GR~IPM+PP+PM+PDRB+KP+TPT, data=data)))


#====chart korelasi transformasi logaritma natural====
datakorelasi1<-data.frame(data$lnGR,data$lnIPM,data$TPT,data$PM,data$lnPP,data$lnKP,data$lnPDRB)
datakorelasi1
chart.Correlation(datakorelasi1, histogram=TRUE, pch=1)

#==VIF data setelah transformasi logaritma, multikol di IPM dan KP==
#Tahun 2018
vif(lm(lnGR18~lnIPM18+lnPP18+PM18+lnPDRB18+TPT18+lnKP18, data = dt.eksplorasi))
#Tahun 2019
vif(lm(lnGR19~lnIPM19+lnPP19+PM19+lnPDRB19+TPT19+lnKP19, data = dt.eksplorasi))
#Tahun 2020
vif(lm(lnGR20~lnIPM20+lnPP20+PM20+lnPDRB20+TPT20+lnKP20, data = dt.eksplorasi))
#Tahun 2021
vif(lm(lnGR21~lnIPM21+lnPP21+PM21+lnPDRB21+TPT21+lnKP21, data = dt.eksplorasi)) 
#Tahun 2022
vif(lm(lnGR22~lnIPM22+lnPP22+PM22+lnPDRB22+TPT22+lnKP22, data = dt.eksplorasi)) 
#Tahun 2018-2022
vif(lm(lnGR~lnIPM+lnPP+PM+lnPDRB+TPT+lnKP, data = data))

#==VIF data setelah transformasi logaritma natural TANPA IPM==
#Tahun 2018
vif(lm(lnGR18~lnPP18+PM18+lnPDRB18+TPT18+lnKP18, data = dt.eksplorasi))
#Tahun 2019
vif(lm(lnGR19~lnPP19+PM19+lnPDRB19+TPT19+lnKP19, data = dt.eksplorasi))
#Tahun 2020
vif(lm(lnGR20~lnPP20+PM20+lnPDRB20+TPT20+lnKP20, data = dt.eksplorasi))
#Tahun 2021
vif(lm(lnGR21~lnPP21+PM21+lnPDRB21+TPT21+lnKP21, data = dt.eksplorasi)) 
#Tahun 2022
vif(lm(lnGR22~lnPP22+PM22+lnPDRB22+TPT22+lnKP22, data = dt.eksplorasi)) 
#Tahun 2018-2022
vif(lm(lnGR~lnPP+PM+lnPDRB+TPT+lnKP, data = data))


#==== Regresi data panel =====
datapanel <- pdata.frame(data, index = c("KOTKAB","Tahun"))
head(datapanel)
model1 <- lnGR~lnIPM+TPT+PM+lnPP+lnPDRB
model2 <- lnGR~TPT+PM+lnPP+lnKP+lnPDRB

#Model CEM
mod.cem <- plm(model2, data=datapanel, model = "pooling")
summary(mod.cem)
model_cem<-summary(mod.cem)
mean(model_cem$residuals^2)

#Model FEM
mod.fem <- plm(model2, data=datapanel, model = "within")
summary(mod.fem)
fixef(mod.fem) #Melihat efek individu setiap intercept

#Model REM
mod.rem <- plm(model2, data=datapanel, model = "random")
summary(mod.rem)
ranef(mod.rem) #Melihat efek individu

#Model CEM VS MODEL FEM
#Uji Chow
pFtest(mod.fem, mod.cem) #p-value < 0.05 = dipilih model FEM

#Model FEM VS Model REM
#Uji Haussman
phtest(mod.rem, mod.fem) #p-value < 0.05 = pilih model FEM

#Efek dua arah
plmtest(mod.fem, effect = "twoways", type = "bp") #p-value < 0.05; Minimal terdapat salah satu antara efek individu dan waktu yang signifikan
plmtest(mod.fem, effect = "individual", type = "bp") #p-value < 0.05; efek individu signifikan
plmtest(mod.fem, effect = "time", type = "bp") #p-value > 0.05;  efek waktu tidak signifikan

#Uji Pesaran
pcdtest(mod.fem, test = c("cd")) #p-value > 0.05; tidak ada korelasi spasial antar individu

#pemeriksaan asumsi
#kenormalan sisaan, kolmogorov-smirnov, p-value>0.05 = tak tolak H0 = sisaan menyebar normal
galat<- mod.fem$residuals
ks.test(galat, "pnorm", mean=mean(galat), sd=sd(galat))
galat
hist(galat)
qqnorm(galat)
qqline(galat)
#heteroskedastisitas/kehomogenan ragam sisaan, breusch-pagan, p-value>0.05 = tak tolak H0 = ragam sisaan homogen
galaat <- abs(galat)
regresi <- lm(galaat~TPT+PM+lnPP+lnKP+lnPDRB, data=datapanel)
summary(regresi)
bptest(mod.fem)
#uji autokorelasi/kebebasan sisaan, runs test, p-value>0.05, tak tolak H0 = sisaan saling bebas
RunsTest(galat)
err.sar<-residuals(mod.fem)
RunsTest(err.sar)
runs.test(err.sar)

#selain pakai R-squared, kalau pake MAPE, model dikatakan baik bila nilai MAPE <10%

.

SOURCE CODE 2

library(readxl)
library(plm)
library(performance)
library(normtest)
library(nortest)
library(pcse)

#memasukkan DATABASE
analisadata <- read_excel("Tugas/R/IPBUniversity/csv/analisadata.xlsx")
head(analisadata)
View(analisadata)

#melihat struktur datanya
str(analisadata)

#membuat analisa deskriptif terlebih dahulu agar ada gambaran umumnya
summary(analisadata)

# menggunakan index mines karena kolom 1 (kotkab) sampai kolom 2 (tahun) tidak kita gunakan
cor(analisadata[,-c(1:2)])

# membuat data panel
paneldata <- pdata.frame(analisadata, index=c("TAHUN"))

# membuat model
# karena data angkanya ada yang sampai jutaan (uang)
# maka MODEL kita bagi menjadi 2 model (model angka biasa dan model jutaan)

### --------------- POOLING MODEL -----------------------------
#MODEL 1 
model1 <- Y~IP+IK+TPAK+TPT+PM+KP+PDRB+SD+SMP+SMA+AHH+RLS+PK+IPM
## Pooled
pooled1 <- plm(model1, paneldata, model = "pooling")
summary(pooled1)
## setiap model kita cek asumsinya
check_collinearity(pooled1)
## menyimpan residnya
residpooled1 <- pooled1$residuals
## kita uji normalitasnya
jb.norm.test(residpooled1)
## kita uji autocorelasinya
check_autocorrelation(pooled1)
## kita uji heteroskedsitas
check_heteroscedasticity(pooled1)

#MODEL 2 
model2 <- Y~PP+UMK
## Pooled
pooled2 <- plm(model2, paneldata, model = "pooling")
summary(pooled2)
## setiap model kita cek asumsinya
check_collinearity(pooled2)
## menyimpan residnya
residpooled2 <- pooled1$residuals
## kita uji normalitasnya
jb.norm.test(residpooled2)
## kita uji autocorelasinya
check_autocorrelation(pooled2)
## kita uji heteroskedsitas
check_heteroscedasticity(pooled2)

### --------------- FIXED EFFECT MODEL -----------------------------
#MODEL 1
fixed1 <-plm(model1, paneldata,model="within", effect = "individual")
summary(fixed1)
## setiap model kita cek asumsinya
## menyimpan residnya
residfixed1 <- fixed1$residuals
## kita uji normalitasnya
jb.norm.test(residfixed1)
## kita uji autocorelasinya
check_autocorrelation(fixed1)
## kita uji heteroskedsitas
check_heteroscedasticity(fixed1)
pwartest(model1, data=paneldata)

#MODEL 2
fixed2 <-plm(model2, paneldata,model="within", effect = "individual")
summary(fixed2)
## setiap model kita cek asumsinya
## menyimpan residnya
residfixed2 <- fixed2$residuals
## kita uji normalitasnya
jb.norm.test(residfixed2)
## kita uji autocorelasinya
check_autocorrelation(fixed2)
## kita uji heteroskedsitas
check_heteroscedasticity(fixed2)
pwartest(model2, data=paneldata)

### --------------- RANDOM EFFECTS -----------------------------
#MODEL 1 tidak bisa dites karena RANDOM hanya bisa max 5 variable
#MODEL 2
random2 <-plm(model2, paneldata,model="random", effect = "individual")
summary(random2)
## setiap model kita cek asumsinya
## menyimpan residnya
residrandom2 <- random2$residuals
## kita uji normalitasnya
jb.norm.test(residrandom2)
## kita uji autocorelasinya
check_autocorrelation(random2)
## kita uji heteroskedsitas
check_heteroscedasticity(random2)

##chow test : pilih antara pooled vs fixed
## kalau hasilnya signifikan berarti yang menang fixed
chow_panel1 <- pFtest(fixed1,pooled1)
chow_panel1
chow_panel2 <- pFtest(fixed2,pooled2)
chow_panel2
## hasilnya signifikan berarti yang manang fixed

##hausman test: pilih antara random vs fixed
hausman_panel <- phtest(fixed2, random2)
hausman_panel

##breuschpagan test : pilih antara pooled vs random
bp_panel<-plmtest(pooled1, type=c("bp"))
bp_panel

model2 <- Y~PP+UMK
gls <- pggls(Y~PP+UMK+lag(Y), data=paneldata, model = "pooling")
summary(gls)
residgls<-gls$residual
jb.norm.test(residgls)

#PCSE
ols <- lm(model1, analisadata)
pcse <- pcse(ols, groupN = paneldata$Y, groupT = paneldata$TAHUN)
summary(pcse)

.

RESULT / HASILNYA


# A tibble: 6 × 19
  KOTKAB   TAHUN     Y    IP    IK  TPAK   TPT    PM    KP  PDRB            PP    SD   SMP   SMA   AHH   RLS    PK      UMK   IPM
  <chr>    <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>         <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>    <dbl> <dbl>
1 Bogor     2018 0.416  60.8  78.2  62.8  9.83  7.14  2155 25373 7377983448292  99.2  84.0  52.8  70.9  7.88 10323 3483667.  69.7
2 NA        2019 0.404  62.3  78.5  65.4  9.11  6.66  2201 26709 7879058217777  99.3  83.4  53.3  71.0  8.29 10683 3763406.  70.6
3 NA        2020 0.381  62.3  78.7  62.6 14.3   7.69  2246 28476 8405493380579  99.1  84.3  53.2  71.2  8.3  10317 4083670   70.4
4 NA        2021 0.396  62.4  79.0  62.6 12.2   8.13  2025 29072 8708210000000  99.2  85.3  54.7  71.4  8.31 10410 4217206   70.6
5 NA        2022 0.4    62.5  79.5  63.8 10.6   7.73  1987 30173 7738408726077  95.9  84.4  48.4  71.6  8.34 10860 4217206   71.2
6 Sukabumi  2018 0.329  56.6  77.7  62.7  7.84  6.76   594 17938 3842321046647  99.4  81.9  49.7  70.5  6.8   8618 2583557.  66.0
> 
> str(analisadata)
tibble [135 × 19] (S3: tbl_df/tbl/data.frame)
 $ KOTKAB: chr [1:135] "Bogor" NA NA NA ...
 $ TAHUN : num [1:135] 2018 2019 2020 2021 2022 ...
 $ Y     : num [1:135] 0.416 0.404 0.381 0.396 0.4 0.329 0.347 0.334 0.343 0.309 ...
 $ IP    : num [1:135] 60.8 62.3 62.3 62.4 62.5 ...
 $ IK    : num [1:135] 78.2 78.5 78.7 79 79.5 ...
 $ TPAK  : num [1:135] 62.8 65.4 62.6 62.5 63.8 ...
 $ TPT   : num [1:135] 9.83 9.11 14.29 12.22 10.64 ...
 $ PM    : num [1:135] 7.14 6.66 7.69 8.13 7.73 6.76 6.22 7.09 7.7 7.34 ...
 $ KP    : num [1:135] 2155 2201 2246 2025 1987 ...
 $ PDRB  : num [1:135] 25373 26709 28476 29072 30173 ...
 $ PP    : num [1:135] 7.38e+12 7.88e+12 8.41e+12 8.71e+12 7.74e+12 ...
 $ SD    : num [1:135] 99.2 99.3 99.1 99.2 95.9 ...
 $ SMP   : num [1:135] 84 83.5 84.3 85.3 84.4 ...
 $ SMA   : num [1:135] 52.8 53.3 53.2 54.7 48.4 ...
 $ AHH   : num [1:135] 70.9 71 71.2 71.4 71.7 ...
 $ RLS   : num [1:135] 7.88 8.29 8.3 8.31 8.34 6.8 7.02 7.07 7.1 7.11 ...
 $ PK    : num [1:135] 10323 10683 10317 10410 10860 ...
 $ UMK   : num [1:135] 3483667 3763406 4083670 4217206 4217206 ...
 $ IPM   : num [1:135] 69.7 70.7 70.4 70.6 71.2 ...
> #membuat analisa deskriptif terlebih dahulu agar ada gambaran umumnya
> summary(analisadata)
    KOTKAB              TAHUN            Y                IP              IK             TPAK            TPT        
 Length:135         Min.   :2018   Min.   :0.2840   Min.   :53.88   Min.   :75.32   Min.   :55.74   Min.   : 1.560  
 Class :character   1st Qu.:2019   1st Qu.:0.3415   1st Qu.:58.38   1st Qu.:78.79   1st Qu.:62.65   1st Qu.: 7.120  
 Mode  :character   Median :2020   Median :0.3630   Median :62.33   Median :80.23   Median :64.71   Median : 8.750  
                    Mean   :2020   Mean   :0.3691   Mean   :64.00   Mean   :80.40   Mean   :65.05   Mean   : 8.592  
                    3rd Qu.:2021   3rd Qu.:0.3970   3rd Qu.:69.25   3rd Qu.:82.19   3rd Qu.:66.96   3rd Qu.:10.015  
                    Max.   :2022   Max.   :0.4890   Max.   :77.33   Max.   :85.35   Max.   :79.92   Max.   :14.290  
       PM               KP               PDRB             PP                  SD              SMP             SMA       
 Min.   : 2.070   Min.   :  393.0   Min.   :12869   Min.   :7.603e+11   Min.   : 94.59   Min.   :69.93   Min.   :43.45  
 1st Qu.: 6.655   1st Qu.:  848.5   1st Qu.:17285   1st Qu.:2.283e+12   1st Qu.: 96.22   1st Qu.:79.55   1st Qu.:52.60  
 Median : 8.260   Median : 1435.0   Median :21601   Median :2.951e+12   Median : 98.25   Median :81.28   Median :60.31  
 Mean   : 8.279   Mean   : 3984.6   Mean   :29202   Mean   :3.340e+12   Mean   : 97.85   Mean   :81.31   Mean   :60.41  
 3rd Qu.:10.350   3rd Qu.: 6787.0   3rd Qu.:32532   3rd Qu.:4.153e+12   3rd Qu.: 99.11   3rd Qu.:84.22   3rd Qu.:67.96  
 Max.   :13.130   Max.   :15798.0   Max.   :85820   Max.   :8.708e+12   Max.   :100.00   Max.   :90.07   Max.   :77.18  
      AHH             RLS               PK             UMK               IPM       
 Min.   :68.96   Min.   : 5.980   Min.   : 7597   Min.   :1558794   Min.   :64.62  
 1st Qu.:71.22   1st Qu.: 7.390   1st Qu.: 9412   1st Qu.:2016580   1st Qu.:68.24  
 Median :72.15   Median : 8.170   Median :10414   Median :2678029   Median :70.86  
 Mean   :72.26   Mean   : 8.545   Mean   :10836   Mean   :2866903   Mean   :71.76  
 3rd Qu.:73.42   3rd Qu.: 9.555   3rd Qu.:11612   3rd Qu.:3604240   3rd Qu.:74.33  
 Max.   :75.48   Max.   :11.470   Max.   :17639   Max.   :4816921   Max.   :82.50  
> # menggunakan index mines karena kolom 1 (kotkab) sampai kolom 2 (tahun) tidak kita gunakan
> cor(analisadata[,-c(1:2)])
               Y          IP          IK        TPAK         TPT          PM          KP         PDRB          PP            SD
Y     1.00000000  0.54411737  0.37292731 -0.47472596  0.37609569 -0.19335493  0.51513932  0.344179113  0.01756584 -0.2632787576
IP    0.54411737  1.00000000  0.73557768 -0.19901775  0.17163682 -0.68487845  0.89069495  0.368166066 -0.04678506 -0.1027328519
IK    0.37292731  0.73557768  1.00000000 -0.22304017  0.34324589 -0.52053673  0.70862606  0.382594557  0.13128853 -0.0597562267
TPAK -0.47472596 -0.19901775 -0.22304017  1.00000000 -0.56747008  0.13486944 -0.21277889 -0.093037651 -0.12069408  0.0452407575
TPT   0.37609569  0.17163682  0.34324589 -0.56747008  1.00000000 -0.09747395  0.30205178  0.361568675  0.28088197 -0.1075479306
PM   -0.19335493 -0.68487845 -0.52053673  0.13486944 -0.09747395  1.00000000 -0.63507070 -0.405814125 -0.24318689 -0.2596669474
KP    0.51513932  0.89069495  0.70862606 -0.21277889  0.30205178 -0.63507070  1.00000000  0.403105994  0.01095925 -0.0976565077
PDRB  0.34417911  0.36816607  0.38259456 -0.09303765  0.36156868 -0.40581412  0.40310599  1.000000000  0.27074866  0.0082376178
PP    0.01756584 -0.04678506  0.13128853 -0.12069408  0.28088197 -0.24318689  0.01095925  0.270748660  1.00000000  0.3748156530
SD   -0.26327876 -0.10273285 -0.05975623  0.04524076 -0.10754793 -0.25966695 -0.09765651  0.008237618  0.37481565  1.0000000000
SMP  -0.28681015 -0.01959959 -0.09367943  0.29130027 -0.36772600 -0.05904069 -0.20317286 -0.306159640  0.20420023  0.3403163670
SMA   0.25895612  0.64946320  0.43137138 -0.02293289 -0.06090017 -0.30834836  0.63606426  0.159482195 -0.37027663 -0.2098399434
AHH   0.37277826  0.73558381  0.99999914 -0.22294782  0.34327974 -0.52052833  0.70859683  0.382745854  0.13133939 -0.0597486001
RLS   0.57890523  0.98380481  0.75730926 -0.23348160  0.23431249 -0.67640794  0.89573393  0.380187845 -0.03836342 -0.1191404301
PK    0.40037341  0.77444408  0.74072400 -0.12692133  0.25792519 -0.68852960  0.82098746  0.560016112  0.18927131 -0.0007443967
UMK   0.27552184  0.42412008  0.53614761 -0.23939005  0.54083097 -0.56399614  0.39566197  0.563300310  0.48507163  0.1840636440
IPM   0.51066444  0.95046685  0.83623579 -0.19425062  0.25060540 -0.71405331  0.89925810  0.485648045  0.05616496 -0.0789936791
             SMP           SMA         AHH         RLS            PK           UMK         IPM
Y    -0.28681015  0.2589561229  0.37277826  0.57890523  0.4003734079  0.2755218380  0.51066444
IP   -0.01959959  0.6494632024  0.73558381  0.98380481  0.7744440845  0.4241200751  0.95046685
IK   -0.09367943  0.4313713833  0.99999914  0.75730926  0.7407240032  0.5361476122  0.83623579
TPAK  0.29130027 -0.0229328894 -0.22294782 -0.23348160 -0.1269213280 -0.2393900489 -0.19425062
TPT  -0.36772600 -0.0609001652  0.34327974  0.23431249  0.2579251880  0.5408309703  0.25060540
PM   -0.05904069 -0.3083483585 -0.52052833 -0.67640794 -0.6885296013 -0.5639961377 -0.71405331
KP   -0.20317286  0.6360642628  0.70859683  0.89573393  0.8209874635  0.3956619729  0.89925810
PDRB -0.30615964  0.1594821954  0.38274585  0.38018785  0.5600161125  0.5633003100  0.48564804
PP    0.20420023 -0.3702766304  0.13133939 -0.03836342  0.1892713074  0.4850716334  0.05616496
SD    0.34031637 -0.2098399434 -0.05974860 -0.11914043 -0.0007443967  0.1840636440 -0.07899368
SMP   1.00000000 -0.0715105512 -0.09377835 -0.09623895 -0.1575585733 -0.1050651356 -0.09647659
SMA  -0.07151055  1.0000000000  0.43133541  0.64397905  0.5500878373  0.0005546047  0.63366069
AHH  -0.09377835  0.4313354071  1.00000000  0.75729382  0.7407709224  0.5362003310  0.83626051
RLS  -0.09623895  0.6439790454  0.75729382  1.00000000  0.7744610381  0.4801706539  0.94532303
PK   -0.15755857  0.5500878373  0.74077092  0.77446104  1.0000000000  0.6034661397  0.92011965
UMK  -0.10506514  0.0005546047  0.53620033  0.48017065  0.6034661397  1.0000000000  0.55280748
IPM  -0.09647659  0.6336606856  0.83626051  0.94532303  0.9201196482  0.5528074778  1.00000000
> # membuat data panel
> paneldata <- pdata.frame(analisadata, index=c("TAHUN"))
> ### --------------- POOLING MODEL -----------------------------
> #MODEL 1 
> model1 <- Y~IP+IK+TPAK+TPT+PM+KP+PDRB+SD+SMP+SMA+AHH+RLS+PK+IPM
> ## Pooled
> pooled1 <- plm(model1, paneldata, model = "pooling")
> summary(pooled1)
Pooling Model

Call:
plm(formula = model1, data = paneldata, model = "pooling")

Balanced Panel: n = 5, T = 27, N = 135

Residuals:
       Min.     1st Qu.      Median     3rd Qu.        Max. 
-0.05825408 -0.01655490  0.00050312  0.01410330  0.08217751 

Coefficients:
               Estimate  Std. Error t-value  Pr(>|t|)    
(Intercept)  4.9017e+01  2.4619e+01  1.9910 0.0487533 *  
IP          -4.6593e-03  4.8155e-03 -0.9676 0.3352129    
IK           1.5547e+00  8.0112e-01  1.9407 0.0546427 .  
TPAK        -4.0314e-03  8.1451e-04 -4.9495 2.456e-06 ***
TPT         -1.6033e-03  1.5227e-03 -1.0529 0.2944989    
PM           7.0927e-03  1.4189e-03  4.9988 1.988e-06 ***
KP           8.9953e-07  1.5027e-06  0.5986 0.5505693    
PDRB         4.6191e-07  1.6586e-07  2.7849 0.0062237 ** 
SD          -1.9517e-03  1.7524e-03 -1.1137 0.2676233    
SMP          1.2427e-04  7.3621e-04  0.1688 0.8662358    
SMA         -1.3785e-03  3.9824e-04 -3.4616 0.0007446 ***
AHH         -2.4042e+00  1.2322e+00 -1.9512 0.0533689 .  
RLS          3.3755e-02  1.0466e-02  3.2252 0.0016221 ** 
PK          -3.8616e-06  8.8256e-06 -0.4375 0.6625028    
IPM          8.2882e-03  1.0207e-02  0.8120 0.4183816    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Total Sum of Squares:    0.22415
Residual Sum of Squares: 0.081053
R-Squared:      0.63839
Adj. R-Squared: 0.59621
F-statistic: 15.1323 on 14 and 120 DF, p-value: < 2.22e-16
> ## setiap model kita cek asumsinya
> check_collinearity(pooled1)
# Check for Multicollinearity

Low Correlation

 Term  VIF           VIF 95% CI Increased SE Tolerance Tolerance 95% CI
 TPAK 1.77 [    1.47,     2.25]         1.33      0.57     [0.45, 0.68]
  TPT 2.47 [    1.99,     3.18]         1.57      0.41     [0.31, 0.50]
   PM 3.08 [    2.44,     3.99]         1.75      0.33     [0.25, 0.41]
 PDRB 1.92 [    1.58,     2.45]         1.39      0.52     [0.41, 0.63]
   SD 1.57 [    1.32,     1.98]         1.25      0.64     [0.50, 0.75]
  SMP 2.05 [    1.68,     2.62]         1.43      0.49     [0.38, 0.59]
  SMA 2.46 [    1.98,     3.16]         1.57      0.41     [0.32, 0.50]

High Correlation

 Term      VIF           VIF 95% CI Increased SE Tolerance Tolerance 95% CI
   IP   206.36 [  153.92,   276.79]        14.37  4.85e-03     [0.00, 0.01]
   IK 6.19e+05 [4.61e+05, 8.30e+05]       786.56  1.62e-06     [0.00, 0.00]
   KP    10.08 [    7.65,    13.38]         3.17      0.10     [0.07, 0.13]
  AHH 6.19e+05 [4.61e+05, 8.30e+05]       786.47  1.62e-06     [0.00, 0.00]
  RLS    45.75 [   34.24,    61.25]         6.76      0.02     [0.02, 0.03]
   PK    78.87 [   58.91,   105.69]         8.88      0.01     [0.01, 0.02]
  IPM   460.23 [  343.09,   617.48]        21.45  2.17e-03     [0.00, 0.00]
> ## menyimpan residnya
> residpooled1 <- pooled1$residuals
> ## kita uji normalitasnya
> jb.norm.test(residpooled1)

	Jarque-Bera test for normality

data:  residpooled1
JB = 6.9199, p-value = 0.0335

> ## kita uji autocorelasinya
> check_autocorrelation(pooled1)
Warning: Autocorrelated residuals detected (p = 0.012).
> ## kita uji heteroskedsitas
> check_heteroscedasticity(pooled1)
Warning: Heteroscedasticity (non-constant error variance) detected (p = 0.022).
> #MODEL 2 
> model2 <- Y~PP+UMK
> ## Pooled
> pooled2 <- plm(model2, paneldata, model = "pooling")
> summary(pooled2)
Pooling Model

Call:
plm(formula = model2, data = paneldata, model = "pooling")

Balanced Panel: n = 5, T = 27, N = 135

Residuals:
      Min.    1st Qu.     Median    3rd Qu.       Max. 
-0.0720909 -0.0266696 -0.0040942  0.0241757  0.1213524 

Coefficients:
               Estimate  Std. Error t-value  Pr(>|t|)    
(Intercept)  3.3834e-01  1.1042e-02 30.6416 < 2.2e-16 ***
PP          -3.7674e-15  2.3518e-15 -1.6019 0.1115735    
UMK          1.5125e-08  4.1050e-09  3.6845 0.0003333 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Total Sum of Squares:    0.22415
Residual Sum of Squares: 0.20318
R-Squared:      0.093533
Adj. R-Squared: 0.079799
F-statistic: 6.81019 on 2 and 132 DF, p-value: 0.0015319
> ## setiap model kita cek asumsinya
> check_collinearity(pooled2)
# Check for Multicollinearity

Low Correlation

 Term  VIF   VIF 95% CI Increased SE Tolerance Tolerance 95% CI
   PP 1.31 [1.13, 1.72]         1.14      0.76     [0.58, 0.88]
  UMK 1.31 [1.13, 1.72]         1.14      0.76     [0.58, 0.88]
> ## menyimpan residnya
> residpooled2 <- pooled1$residuals
> ## kita uji normalitasnya
> jb.norm.test(residpooled2)

	Jarque-Bera test for normality

data:  residpooled2
JB = 6.9199, p-value = 0.0325

> ## kita uji autocorelasinya
> check_autocorrelation(pooled2)
Warning: Autocorrelated residuals detected (p = 0.008).
> ## kita uji heteroskedsitas
> check_heteroscedasticity(pooled2)
OK: Error variance appears to be homoscedastic (p = 0.763).
> ### --------------- FIXED EFFECT MODEL -----------------------------
> #MODEL 1
> fixed1 <-plm(model1, paneldata,model="within", effect = "individual")
> summary(fixed1)
Oneway (individual) effect Within Model

Call:
plm(formula = model1, data = paneldata, effect = "individual", 
    model = "within")

Balanced Panel: n = 5, T = 27, N = 135

Residuals:
      Min.    1st Qu.     Median    3rd Qu.       Max. 
-0.0575176 -0.0181776 -0.0012779  0.0140140  0.0817065 

Coefficients:
        Estimate  Std. Error t-value  Pr(>|t|)    
IP   -3.6874e-03  4.9122e-03 -0.7507 0.4543700    
IK    1.5525e+00  8.0410e-01  1.9307 0.0559632 .  
TPAK -4.3871e-03  9.5261e-04 -4.6054 1.062e-05 ***
TPT  -1.2442e-03  1.9230e-03 -0.6470 0.5188964    
PM    5.6756e-03  1.6587e-03  3.4218 0.0008599 ***
KP    1.5001e-06  1.6191e-06  0.9265 0.3561223    
PDRB  4.4535e-07  1.6740e-07  2.6605 0.0089085 ** 
SD   -3.6799e-03  2.0317e-03 -1.8113 0.0726843 .  
SMP   1.8798e-04  7.3961e-04  0.2542 0.7998206    
SMA  -1.3297e-03  4.0437e-04 -3.2884 0.0013341 ** 
AHH  -2.4009e+00  1.2368e+00 -1.9412 0.0546624 .  
RLS   2.9070e-02  1.0698e-02  2.7174 0.0075890 ** 
PK   -2.8743e-06  8.8959e-06 -0.3231 0.7471981    
IPM   6.5157e-03  1.0319e-02  0.6315 0.5289881    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Total Sum of Squares:    0.22069
Residual Sum of Squares: 0.078163
R-Squared:      0.64582
Adj. R-Squared: 0.59086
F-statistic: 15.1083 on 14 and 116 DF, p-value: < 2.22e-16
> ## setiap model kita cek asumsinya
> ## menyimpan residnya
> residfixed1 <- fixed1$residuals
> ## kita uji normalitasnya
> jb.norm.test(residfixed1)

	Jarque-Bera test for normality

data:  residfixed1
JB = 7.6861, p-value = 0.027

> ## kita uji autocorelasinya
> check_autocorrelation(fixed1)
Warning: Autocorrelated residuals detected (p = 0.030).
> ## kita uji heteroskedsitas
> check_heteroscedasticity(fixed1)
OK: Error variance appears to be homoscedastic (p = 0.078).
> pwartest(model1, data=paneldata)

	Wooldridge's test for serial correlation in FE panels

data:  plm.model
F = 25.373, df1 = 1, df2 = 128, p-value = 1.571e-06
alternative hypothesis: serial correlation

> #MODEL 2
> fixed2 <-plm(model2, paneldata,model="within", effect = "individual")
> summary(fixed2)
Oneway (individual) effect Within Model

Call:
plm(formula = model2, data = paneldata, effect = "individual", 
    model = "within")

Balanced Panel: n = 5, T = 27, N = 135

Residuals:
      Min.    1st Qu.     Median    3rd Qu.       Max. 
-0.0698032 -0.0272597 -0.0041819  0.0234140  0.1189913 

Coefficients:
       Estimate  Std. Error t-value  Pr(>|t|)    
PP  -3.7202e-15  2.4058e-15 -1.5463 0.1244901    
UMK  1.5343e-08  4.2904e-09  3.5760 0.0004931 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Total Sum of Squares:    0.22069
Residual Sum of Squares: 0.20053
R-Squared:      0.09135
Adj. R-Squared: 0.048757
F-statistic: 6.43413 on 2 and 128 DF, p-value: 0.0021747
> ## setiap model kita cek asumsinya
> ## menyimpan residnya
> residfixed2 <- fixed2$residuals
> ## kita uji normalitasnya
> jb.norm.test(residfixed2)

	Jarque-Bera test for normality

data:  residfixed2
JB = 7.5362, p-value = 0.0275

> ## kita uji autocorelasinya
> check_autocorrelation(fixed2)
Warning: Autocorrelated residuals detected (p = 0.010).
> ## kita uji heteroskedsitas
> check_heteroscedasticity(fixed2)
OK: Error variance appears to be homoscedastic (p = 0.603).
> pwartest(model2, data=paneldata)

	Wooldridge's test for serial correlation in FE panels

data:  plm.model
F = 10.469, df1 = 1, df2 = 128, p-value = 0.001545
alternative hypothesis: serial correlation

> ### --------------- RANDOM EFFECTS -----------------------------
> #MODEL 1 tidak bisa dites karena RANDOM hanya bisa max 5 variable
> #MODEL 2
> random2 <-plm(model2, paneldata,model="random", effect = "individual")
> summary(random2)
Oneway (individual) effect Random Effect Model 
   (Swamy-Arora's transformation)

Call:
plm(formula = model2, data = paneldata, effect = "individual", 
    model = "random")

Balanced Panel: n = 5, T = 27, N = 135

Effects:
                   var  std.dev share
idiosyncratic 0.001567 0.039580     1
individual    0.000000 0.000000     0
theta: 0

Residuals:
      Min.    1st Qu.     Median    3rd Qu.       Max. 
-0.0720909 -0.0266696 -0.0040942  0.0241757  0.1213524 

Coefficients:
               Estimate  Std. Error z-value  Pr(>|z|)    
(Intercept)  3.3834e-01  1.1042e-02 30.6416 < 2.2e-16 ***
PP          -3.7674e-15  2.3518e-15 -1.6019 0.1091831    
UMK          1.5125e-08  4.1050e-09  3.6845 0.0002292 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Total Sum of Squares:    0.22415
Residual Sum of Squares: 0.20318
R-Squared:      0.093533
Adj. R-Squared: 0.079799
Chisq: 13.6204 on 2 DF, p-value: 0.0011025
> ## setiap model kita cek asumsinya
> ## menyimpan residnya
> residrandom2 <- random2$residuals
> ## kita uji normalitasnya
> jb.norm.test(residrandom2)

	Jarque-Bera test for normality

data:  residrandom2
JB = 8.8147, p-value = 0.018

> ## kita uji autocorelasinya
> check_autocorrelation(random2)
Warning: Autocorrelated residuals detected (p < .001).
> ## kita uji heteroskedsitas
> check_heteroscedasticity(random2)
OK: Error variance appears to be homoscedastic (p = 0.763).
> ##chow test : pilih antara pooled vs fixed
> ## kalau hasilnya signifikan berarti yang menang fixed
> chow_panel1 <- pFtest(fixed1,pooled1)
> chow_panel1

	F test for individual effects

data:  model1
F = 1.0723, df1 = 4, df2 = 116, p-value = 0.3735
alternative hypothesis: significant effects

Terimakasih

COMMENTS (1)

Leave a Reply

Your email address will not be published. Required fields are marked *