Anggota Kelompok 1
- G1401201004 Alfidhia Rahman Nasa Juhanda
- G1401201050 Hanung Safrizal
- G1401201064 Muhammad Irsyad Robbani
- G1401201081 Zakiyah Azzahra
- G1401201092 Muhammad Zhillan Zakiyyan
lapply(c("kableExtra","tidyverse","plm","tinytex","tseries","lmtest",
"broom", "stargazer", "RColorBrewer"), library, character.only = T)[[1]]## [1] "kableExtra" "stats" "graphics" "grDevices" "utils"
## [6] "datasets" "methods" "base"Data
Data yang digunakan adalah data Tingkat Kemiskinan Dengan 4 Peubah bebas yang diberikan saat responsi. Data ini berjenis data panel dengan individunya adalah 32 Provinsi di Indonesia (34 Provinsi kecuali Provinsi DKI Jakarta dan Kalimantan Utara) selama 13 tahun (Tahun 2008-2020) dengan peubah-peubah sebagai berikut:
KMIS : Kemiskinan (%) (Peubah Respons)
INFR : Infrastrukur (Juta)
PEKO : Pertumbuhan Ekonomi (%)
INFL : Inflasi (%)
PGGR : Pengangguran (%)
Untuk lebih jelasnya, dapat dilihat dari output berikut:
## tibble [416 × 7] (S3: tbl_df/tbl/data.frame)
## $ PROVINSI: Factor w/ 32 levels "Aceh","BABEL",..: 3 3 3 3 3 3 3 3 3 3 ...
## $ TAHUN : Factor w/ 13 levels "2008","2009",..: 1 2 3 4 5 6 7 8 9 10 ...
## $ KMIS : num [1:416] 6.17 5.13 4.88 4.2 3.95 4.49 4.76 5.25 4.15 4.14 ...
## $ INFRA : num [1:416] 430814 484340 336352 298867 348529 ...
## $ PEKO : num [1:416] 5.97 5.33 5.83 6.49 6.65 6.69 6.73 6.03 6.32 5.56 ...
## $ INFL : num [1:416] 9.62 4.37 8.1 3.75 4.71 8.16 8.43 2.75 3.23 3.32 ...
## $ PGGR : num [1:416] 3.31 3.13 3.06 2.95 2.1 1.83 1.9 1.99 1.89 1.48 ...Dapat terlihat bahwa terdapat 416 total observasi, yaitu 32 individu dengan periode 13 tahun. Semua peubah berjenis numerik kontinu. Berikut adalah 13 observasi teratas dari data yang digunakan yang menunjukan semua observasi yang terdapat di Provinsi Bali:
kemis %>% head(13)Metodologi
Metodologi yang digunakan adalah sebagai berikut:
Eksplorasi Data
Menduga parameter model CEM
Menduga parameter model FEM
Melakukan Uji Chow
Menduga Parameter model REM
Melakukan Uji Hausman
Pengujian asumsi terhadap Model terpilih
Jika asumsi tidak terpenuhi, akan dilakukan penanganan asumsi
Interpretasi model terbaik
Berikut adalah dari metode penelitian yang digunakan:
Eksplorasi Data
Plot
ggplot(data=kemis, aes(x=TAHUN, y=KMIS, group = PROVINSI, colour = PROVINSI))+ theme_bw()+
geom_line(size=1.2) +
geom_point( size=3, shape=19, fill="red") +
labs(colour="Provinsi", title = "Tingkat Kemiskinan di Indonesia", subtitle = "Tahun 2008-2020") +
theme(plot.title = element_text(hjust = 0.5,face = "bold"), plot.subtitle = element_text(hjust = 0.5))## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
Interpretasi: Line chart di atas menampilkan sebaran tingkat kemiskinan per provinsi di Indonesia dari tahun 2008 sampai 2020, dapat dilihat bahwa Papua Barat memiliki tingkat kemiskinan paling tinggi, kemudian Provinsi yang memiliki tingkat kemiskinan yang tinggi adalah Maluku dan Papua yang mana ketiga Provinsi tersebut merupakan daerah bagian timur Indonesia. Hal ini mengindikasikan bahwa terdapat perbedaan untuk setiap individu
Map Chart
Map di bawah ini dibuat pada Aplikasi Tableau Public, merupakan sebaran rata-rata kemiskinan di Indonesia pada selama tahun 2008 hingga 2020.
Interpretasi: Sebaran rata-rata kemiskinan di Indonesia yang ditampilkan pada map chart di atas menunjukkan bahwa semakin pekat warna merah daerah tersebut maka semakin tinggi tingkat kemiskinan di provinsi tersebut dan semakin pekat warna hijau daerah tersebut maka semakin rendah tingkat kemiskinan di daerah tersebut. Wilayah Indonesia bagian timur di dominasi oleh warna merah yang artinya masih banyak daerah di bagian timur Indonesia yang memiliki tingkat kemiskinan tinggi, dapat dilihat juga semua provinsi yang terdapat pada Pulau Kalimantan, Pulau Jawa dan Bali memiliki tingkat kemiskinan yang rendah. Hal ini juga mengindikasikan bahwa terdapat perbedaan untuk setiap individu.
Common-Effect Model
Common Effect Model (CEM) atau disebut juga pooled regression model merupakan salah satu model dalam regresi data panel yang mengkombinasikan data cross-section dan time series menjadi satu kesatuan (pooldata) yang kemudian diestimasi modelnya dengan teknik OLS. Akibatnya, perbedaan baik antar individu maupun antar waktu tidak dapat diamati. Dengan kata lain, dalam pendekatan ini dimensi individu maupun waktu tidak diperhatikan. Diasumsikan bahwa perilaku data antar Individu sama dalam berbagai kurun waktu. Persamaan pool model dapat dituliskan sebagai berikut:
\[ \begin{aligned} y_{it} &= \alpha + \beta_{1}x_{1,it}+ \beta_{2}x_{2,it}+...+\beta_{k}x_{k,it}+\varepsilon_{it} \\ \end{aligned} \]
dengan \(y_{ij}\) merupakan variabel dependen individu ke-i dan waktu ke-j. \(\alpha\) adalah intersep, \(x_{p,ij}\) adalah variabel independen ke-p individu ke-i dan waktu ke-j, \(\beta_k\) koefisien variabel bebas, dan \(\varepsilon_{ij}\) merupakan galat individu ke-i dan waktu ke-j.
Selain asumsi kesamaan pengaruh individu dan waktu, asumsi lain yang diperlukan pada model ini adalah asumsi normalitas, homogenitas, dan kebebasan galat.
Model CEM
Program R menyediakan package plm yang dapat memodelkan
data menjadi model CEM.
Dalam pemodelan ini, dapat dibuat hipotesis:
\[ \begin{aligned} H_0 &: \text{CEM adalah model yang tidak layak}\\ H_1 &: \text{CEM adalah model yang layak} \end{aligned} \]
library(plm)
cem <- plm(KMIS ~ INFRA + PEKO + INFL + PGGR, data = kemis, model = "pooling")
summary(cem)## Pooling Model
##
## Call:
## plm(formula = KMIS ~ INFRA + PEKO + INFL + PGGR, data = kemis,
## model = "pooling")
##
## Balanced Panel: n = 32, T = 13, N = 416
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -9.1611 -5.1292 -1.5170 4.5191 19.6026
##
## Coefficients:
## Estimate Std. Error t-value Pr(>|t|)
## (Intercept) 6.6707e+00 1.1961e+00 5.5772 4.438e-08 ***
## INFRA 1.0347e-06 2.9780e-07 3.4746 0.0005661 ***
## PEKO 5.0286e-01 8.2736e-02 6.0779 2.781e-09 ***
## INFL 3.1338e-01 1.0394e-01 3.0151 0.0027284 **
## PGGR 4.8913e-02 1.3531e-01 0.3615 0.7179249
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 17772
## Residual Sum of Squares: 15525
## R-Squared: 0.12645
## Adj. R-Squared: 0.11794
## F-statistic: 14.873 on 4 and 411 DF, p-value: 2.3234e-11Pada dasarnya, pendugaan model yang dilakukan oleh metode CEM menghasilkan model yang sama persis dengan metode Regresi Linear Berganda dengan memandang setiap baris pada data adalah individu yang berbeda-beda. Model regresi yang dihasilkan adalah sebagai berikut:
\[ \begin{aligned} KMIS_{it} &= 6.67 + 1.0347(10^{-6})INFRA_{it}+0.5PEKO_{it}-0.31INFL_{it}+0.049PGGR_{it}+\varepsilon_{it} \\ \end{aligned} \]
Interpretasi:
Setiap kenaikan satu satuan
InfrastrukturmenyebabkanTingkat Kemiskinanmeningkat sebesar1.0347x10^-6persen dengan menganggap peubah lain konstan.Setiap kenaikan satu satuan
Pertumbuhan EkonomimenyebabkanTingkat Kemiskinanmeningkat sebesar0.5persen dengan menganggap peubah lain konstan.Setiap kenaikan satu satuan
InflasimenyebabkanTingkat Kemiskinanmenurun sebesar0.31persen dengan menganggap peubah lain konstan.Setiap kenaikan satu satuan
PengangguranmenyebabkanTingkat Kemiskinanmeningkat sebesar0.049persen dengan menganggap peubah lain konstan.
Untuk semua p-value di bawah 5%, maka terdapat cukup bukti untuk menolak \(H_0\) (kurang dari 5% kemungkinan bahwa \(H_0\) benar).
p-value uji-F < 2.22e-16 << 5%. Artinya, dapat dikatakan bahwa model layak untuk digunakan.
p-value uji-T untuk peubah
INFRA,PEKO, danINFL< 5%. Artinya, dapat dikatakan bahwa peubah tersebut berpengaruh signifikan terhadap model.
\(R^2\) merupakan nilai yang menunjukan kekuatan hubungan antara peubah bebas dan peubah respons. Semakin tinggi nilai \(R^2\), maka semakin baik model tersebut.
- \(R^2\) pada model sebesar 0.12645 atau 12.645% yang mengindikasikan bahwa peubah yang digunakan kurang baik dalam menggambarkan model.
Meskipun sebagian besar parameter model signifikan, tetapi nilai R2 model kurang baik dan fakta bahwa model CEM mengabaikan pengaruh individu dan waktu tidak bisa diabaikan. Kedua faktor tersebut bisa berpengaruh terhadap keheterogenan sisaan pada model. Oleh karena itu, perlu dilakukan perbandingan model CEM dengan FEM dan REM untuk memastikan model terbaik.
Fixed-Effect Model
Fixed-Effect model merupakan salah satu jenis model dalam regresi data panel yang mempertimbangkan adanya pengaruh individu dan waktu. Terdapat asumsi bahwa intersep antar individu dan waktu berbeda tetapi koefisien regresi konstan untuk semua individu dan waktu. Model yang memiliki pengaruh dari salah satu dari individu atau waktu disebut model sisaan satu arah dan model yang memiliki pengaruh dari keduanya disebut model sisaan dua arah. Pendugaan pada model ini biasa dilakukan dengan model Within atau Least Square Dummy Variable (LSDV).
Within vs LSDV (?)
Perlu diperhatikan bahwa model within akan menghasilkan output yang sama seperti model LSDV tanpa referensi atau tanpa intersep. Pembuktiannya dapat dilihat di bawah ini:
fem_within <- plm(KMIS ~ INFRA + PEKO + INFL + PGGR, index = c("PROVINSI", "TAHUN"),
model = "within", effect= "individual", data = kemis)
fem_lsdv_tanpa_referensi <- lm(KMIS ~ INFRA + PEKO + INFL + PGGR + PROVINSI - 1,
data = kemis)
fem_within$coefficients; fem_lsdv_tanpa_referensi$coefficients[1:4]## INFRA PEKO INFL PGGR
## -6.936362e-07 1.829038e-01 1.744263e-01 6.578657e-01## INFRA PEKO INFL PGGR
## -6.936362e-07 1.829038e-01 1.744263e-01 6.578657e-01Dapat terlihat bahwa pengaruh setiap peubah adalah sama. Begitu pun dengan pengaruh setiap individunya:
as.data.frame(cbind(Within=fixef(fem_within, effect = "individual"),
LSDV= fem_lsdv_tanpa_referensi$coefficients[5:36]))Oleh karena itu, kita hanya akan menggunakan satu metode, yaitu model within.
Persamaan Model FEM
Persamaan model Fixed-Effect :
- Pengaruh Individu
\[ \begin{aligned} y_{it} = \alpha + f_i +\beta_1x_{1,it} + \beta_2x_{2,it} + ... + \beta_kx_{k,it} + \varepsilon_{it} \end{aligned} \]
- Pengaruh Waktu
\[ \begin{aligned} y_{it} = \alpha + \lambda_t +\beta_1x_{1,it} + \beta_2x_{2,it} + ... + \beta_kx_{k,it} + \varepsilon_{it} \end{aligned} \]
- Model Dua arah (pengaruh individu dan waktu)
\[ \begin{aligned} y_{it} = \alpha + f_i + \lambda_t +\beta_1x_{1,it} + \beta_2x_{2,it} + ... + \beta_kx_{k,it} + \varepsilon_{it} \end{aligned} \]
FEM INDIVIDUAL
Efek komponen sisaaan satu arah pada individu.
fem <- plm(KMIS ~ INFRA + PEKO + INFL + PGGR, index = c("PROVINSI", "TAHUN"),
model = "within", effect= "individual", data = kemis)
summary(fem)## Oneway (individual) effect Within Model
##
## Call:
## plm(formula = KMIS ~ INFRA + PEKO + INFL + PGGR, data = kemis,
## effect = "individual", model = "within", index = c("PROVINSI",
## "TAHUN"))
##
## Balanced Panel: n = 32, T = 13, N = 416
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -4.76845 -0.95405 -0.11980 0.76873 6.62826
##
## Coefficients:
## Estimate Std. Error t-value Pr(>|t|)
## INFRA -6.9364e-07 1.2802e-07 -5.4184 1.070e-07 ***
## PEKO 1.8290e-01 2.5300e-02 7.2293 2.689e-12 ***
## INFL 1.7443e-01 2.9213e-02 5.9708 5.422e-09 ***
## PGGR 6.5787e-01 7.8248e-02 8.4074 8.533e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 2088
## Residual Sum of Squares: 1025.6
## R-Squared: 0.50881
## Adj. R-Squared: 0.46356
## F-statistic: 98.4064 on 4 and 380 DF, p-value: < 2.22e-16Didapat model sebagai berikut:
\[ \begin{aligned} KMIS_{it} &= f_i-6.9364(10^{-7})INFRA_{it}+0.1829PEKO_{it}-0.17441INFL_{it}+0.6579PGGR_{it}+\varepsilon_{it} \\ \end{aligned} \]
Interpretasi :
Setiap kenaikan satu satuan
InfrastrukturmenyebabkanTingkat Kemiskinanmenurun sebesar6.9364x10^-7persen dengan menganggap peubah lain konstan.Setiap kenaikan satu satuan
Pertumbuhan EkonomimenyebabkanTingkat Kemiskinanmeningkat sebesar0.1829persen dengan menganggap peubah lain konstan.Setiap kenaikan satu satuan
InflasimenyebabkanTingkat Kemiskinanmenurun sebesar0.1744persen dengan menganggap peubah lain konstan.Setiap kenaikan satu satuan
PengangguranmenyebabkanTingkat Kemiskinanmeningkat sebesar0.6579persen dengan menganggap peubah lain konstan.
Untuk semua p-value di bawah 5%, maka terdapat cukup bukti untuk menolak \(H_0\) (kurang dari 5% kemungkinan bahwa \(H_0\) benar).
- Model memiliki
p-valuepada uji F yang lebih kecil dari 0.05 (< 2.22e-16) yang berarti model ini layak untuk digunakan - p-value uji-T untuk peubah
INFRA,PEKO, danINFL< 5%. Artinya, dapat dikatakan bahwa peubah tersebut berpengaruh signifikan terhadap model.
\(R^2\) merupakan nilai yang menunjukan kekuatan hubungan antara peubah bebas dan peubah respons. Semakin tinggi nilai \(R^2\), maka semakin baik model tersebut.
\(R^2\) pada model sebesar 0.50881 atau 50.881% yang mengindikasikan bahwa peubah yang digunakan cukup baik dalam menggambarkan model.
Model dapat menjelaskan sebesar 46.356% dari keragaman data yang dilihat dari \(R^2 Adj\).
Nilai Intersep dari setiap individu
summary(fixef(fem, effect="individual"))## Estimate Std. Error t-value Pr(>|t|)
## Aceh 12.76929 0.90100 14.1724 < 2.2e-16 ***
## BABEL 1.15546 0.63089 1.8315 0.0678101 .
## Bali 1.53673 0.54757 2.8064 0.0052672 **
## Banten -2.43011 1.01291 -2.3991 0.0169149 *
## Bengkulu 13.04266 0.62535 20.8567 < 2.2e-16 ***
## DIY 10.43756 0.59606 17.5110 < 2.2e-16 ***
## Gorontalo 14.13774 0.64218 22.0152 < 2.2e-16 ***
## JABAR 3.27902 1.04650 3.1333 0.0018626 **
## Jambi 3.70584 0.65751 5.6362 3.394e-08 ***
## JATENG 10.57064 0.84614 12.4927 < 2.2e-16 ***
## JATIM 10.36797 0.81561 12.7119 < 2.2e-16 ***
## KALBAR 4.56660 0.69636 6.5578 1.786e-10 ***
## KALSEL 0.50206 0.67719 0.7414 0.4589218
## KALTENG 2.18406 0.64430 3.3898 0.0007727 ***
## KALTIM 0.24362 0.84257 0.2891 0.7726350
## Kepri 0.73340 0.76332 0.9608 0.3372577
## Lampung 11.01549 0.72114 15.2751 < 2.2e-16 ***
## Maluku 14.23684 0.92636 15.3685 < 2.2e-16 ***
## MALUT 3.05684 0.70061 4.3631 1.654e-05 ***
## NTB 13.62690 0.67977 20.0462 < 2.2e-16 ***
## NTT 19.07966 0.67755 28.1597 < 2.2e-16 ***
## PAPUA 24.63917 0.83838 29.3889 < 2.2e-16 ***
## PapuaBarat 21.12064 0.81518 25.9091 < 2.2e-16 ***
## RIAU 3.00880 0.74952 4.0143 7.185e-05 ***
## SULBAR 8.91121 0.59002 15.1033 < 2.2e-16 ***
## SULSEL 5.49197 0.86742 6.3314 6.860e-10 ***
## SULTENG 10.85094 0.69374 15.6411 < 2.2e-16 ***
## SulTenggara 9.98109 0.67618 14.7610 < 2.2e-16 ***
## SULUT 2.02268 0.88163 2.2943 0.0223203 *
## SUMBAR 2.45959 0.79659 3.0876 0.0021656 **
## SUMSEL 9.57356 0.73051 13.1053 < 2.2e-16 ***
## SUMUT 5.63776 0.87574 6.4377 3.662e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Nilai diatas adalah nilai pengaruh konstan dari masing-masing individu yang pada model dapat dituliskan sebagai \(f_i\).
FEM TIME
Efek komponen sisaaan satu arah pada waktu.
fem_time <- plm(KMIS ~ INFRA + PEKO + INFL + PGGR, index = c("PROVINSI", "TAHUN"),
model = "within", effect= "time", data = kemis)
summary(fem_time)## Oneway (time) effect Within Model
##
## Call:
## plm(formula = KMIS ~ INFRA + PEKO + INFL + PGGR, data = kemis,
## effect = "time", model = "within", index = c("PROVINSI",
## "TAHUN"))
##
## Balanced Panel: n = 32, T = 13, N = 416
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -11.1310 -4.3205 -1.7469 4.1329 17.5313
##
## Coefficients:
## Estimate Std. Error t-value Pr(>|t|)
## INFRA 1.9886e-06 3.3696e-07 5.9017 7.693e-09 ***
## PEKO 5.2529e-01 9.4188e-02 5.5771 4.516e-08 ***
## INFL 6.0293e-02 2.0750e-01 0.2906 0.77153
## PGGR -2.5724e-01 1.3818e-01 -1.8616 0.06339 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 16380
## Residual Sum of Squares: 13680
## R-Squared: 0.16481
## Adj. R-Squared: 0.13131
## F-statistic: 19.6832 on 4 and 399 DF, p-value: 8.4451e-15Didapat model sebagai berikut:
\[ \begin{aligned} KMIS_{it} &= \lambda_t+1.989(10^{-6})INFRA_{it}+0.525PEKO_{it}+0.063INFL_{it}-0.257PGGR_{it}+\varepsilon_{it} \\ \end{aligned} \]
Interpretasi :
Setiap kenaikan satu satuan
InfrastrukturmenyebabkanTingkat Kemiskinanmeningkat sebesar1.989x10^-6persen dengan menganggap peubah lain konstan.Setiap kenaikan satu satuan
Pertumbuhan EkonomimenyebabkanTingkat Kemiskinanmeningkat sebesar0.525persen dengan menganggap peubah lain konstan.Setiap kenaikan satu satuan
InflasimenyebabkanTingkat Kemiskinanmeningkat sebesar0.063persen dengan menganggap peubah lain konstan.Setiap kenaikan satu satuan
PengangguranmenyebabkanTingkat Kemiskinanmenurun sebesar0.257persen dengan menganggap peubah lain konstan.
Untuk semua p-value di bawah 5%, maka terdapat cukup bukti untuk menolak \(H_0\) (kurang dari 5% kemungkinan bahwa \(H_0\) benar).
Model memiliki
p-valuepada uji F yang lebih kecil dari 0.05 (< 8.4451e-15) yang berarti model ini layak untuk digunakanp-value uji-T untuk peubah
INFRAdanPEKO< 5% danINFL< 10%. Artinya, dapat dikatakan bahwa peubah tersebut berpengaruh signifikan terhadap model.
\(R^2\) merupakan nilai yang menunjukan kekuatan hubungan antara peubah bebas dan peubah respons. Semakin tinggi nilai \(R^2\), maka semakin baik model tersebut.
\(R^2\) pada model sebesar 0.16481 atau 16.481% yang mengindikasikan bahwa peubah yang digunakan tidak baik dalam menggambarkan model.
Model dapat menjelaskan sebesar 13.131% dari keragaman data yang dilihat dari \(R^2 Adj\) yang mengindikasikan bahwa peubah yang digunakan tidak baik dalam menggambarkan model.
Nilai Intersep dari setiap waktu
summary(fixef(fem_time, effect="time"))## Estimate Std. Error t-value Pr(>|t|)
## 2008 13.2114 2.9719 4.4454 1.138e-05 ***
## 2009 12.4802 1.7951 6.9524 1.474e-11 ***
## 2010 10.5120 2.2926 4.5852 6.079e-06 ***
## 2011 9.3642 1.8790 4.9836 9.318e-07 ***
## 2012 8.2289 1.8189 4.5241 8.011e-06 ***
## 2013 8.1076 2.3144 3.5030 0.0005120 ***
## 2014 8.0184 2.3076 3.4747 0.0005674 ***
## 2015 8.2809 1.7461 4.7426 2.942e-06 ***
## 2016 4.3752 1.8103 2.4169 0.0161026 *
## 2017 5.9215 1.6797 3.5253 0.0004720 ***
## 2018 5.2166 1.6680 3.1274 0.0018928 **
## 2019 4.4697 1.6320 2.7389 0.0064413 **
## 2020 10.1492 1.4282 7.1064 5.523e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Nilai diatas adalah nilai pengaruh konstan dari masing-masing waktu yang pada model dapat dituliskan sebagai \(\lambda_t\).
FEM TWO WAYS
Efek komponen sisaaan satu arah pada individu.
fem_twoways <- plm(KMIS ~ INFRA + PEKO + INFL + PGGR, index = c("PROVINSI", "TAHUN"),
model = "within", effect= "twoways", data = kemis)
summary(fem_twoways)## Twoways effects Within Model
##
## Call:
## plm(formula = KMIS ~ INFRA + PEKO + INFL + PGGR, data = kemis,
## effect = "twoways", model = "within", index = c("PROVINSI",
## "TAHUN"))
##
## Balanced Panel: n = 32, T = 13, N = 416
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -3.718143 -0.713539 0.050615 0.682921 4.957867
##
## Coefficients:
## Estimate Std. Error t-value Pr(>|t|)
## INFRA -6.2411e-07 1.5858e-07 -3.9356 9.922e-05 ***
## PEKO 1.5398e-01 2.4514e-02 6.2813 9.478e-10 ***
## INFL 1.1098e-01 4.7496e-02 2.3367 0.01999 *
## PGGR -6.3412e-02 7.8861e-02 -0.8041 0.42186
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 695.57
## Residual Sum of Squares: 592.21
## R-Squared: 0.14859
## Adj. R-Squared: 0.039851
## F-statistic: 16.0562 on 4 and 368 DF, p-value: 3.9613e-12Didapat model sebagai berikut:
\[ \begin{aligned} KMIS_{it} &= f_i+\lambda_t-6.241(10^{-7})INFRA_{it}+0.154PEKO_{it}+0.111INFL_{it}-0.063PGGR_{it}+\varepsilon_{it} \\ \end{aligned} \]
Interpretasi :
Setiap kenaikan satu satuan
InfrastrukturmenyebabkanTingkat Kemiskinanmenurun sebesar6.241x10^-7persen dengan menganggap peubah lain konstan.Setiap kenaikan satu satuan
Pertumbuhan EkonomimenyebabkanTingkat Kemiskinanmeningkat sebesar0.154persen dengan menganggap peubah lain konstan.Setiap kenaikan satu satuan
InflasimenyebabkanTingkat Kemiskinanmeningkat sebesar0.111persen dengan menganggap peubah lain konstan.Setiap kenaikan satu satuan
PengangguranmenyebabkanTingkat Kemiskinanmenurun sebesar0.063persen dengan menganggap peubah lain konstan.
Untuk semua p-value di bawah 5%, maka terdapat cukup bukti untuk menolak \(H_0\) (kurang dari 5% kemungkinan bahwa \(H_0\) benar).
- Model memiliki
p-valuepada uji F yang lebih kecil dari 0.05 (< 2.22e-16) yang berarti model ini layak untuk digunakan - p-value uji-T untuk peubah
INFRA,PEKO, danINFL< 5%. Artinya, dapat dikatakan bahwa peubah tersebut berpengaruh signifikan terhadap model.
\(R^2\) merupakan nilai yang menunjukan kekuatan hubungan antara peubah bebas dan peubah respons. Semakin tinggi nilai \(R^2\), maka semakin baik model tersebut.
\(R^2\) pada model sebesar 0.14859 atau 14.859% yang mengindikasikan bahwa peubah yang digunakan tidak baik dalam menggambarkan model.
Model dapat menjelaskan sebesar 3.985% dari keragaman data yang dilihat dari \(R^2 Adj\) yang mengindikasikan bahwa peubah yang digunakan tidak baik dalam menggambarkan model.
Nilai Intersep dari setiap individu dan waktu
data.frame(summary(fixef(fem_twoways, effect="twoways")))Nilai diatas adalah nilai pengaruh konstan dari gabungan masing-masing individu dan waktu yang pada model dapat dituliskan sebagai \(f_i+\lambda_t\).
Uji-Chow
Uji Chow adalah uji formal untuk menguji ada atau tidaknya pengaruh spesifik individu dan pengaruh spesifik waktu pada data panel.
Hipotesis pada uji ini yaitu :
- Untuk menguji pengaruh individu
\[ \begin{aligned} H_{01} &: \text{Tidak ada pengaruh spesifik individu} \\ H_{11} &: \text{Ada pengaruh spesifik individu} \end{aligned} \]
- Untuk menguji pengaruh waktu
\[ \begin{aligned} H_{02} &: \text{Tidak ada pengaruh spesifik waktu} \\ H_{12} &: \text{Ada pengaruh spesifik waktu} \end{aligned} \]
Statsitik Uji yang digunakan :
- Pengaruh individu
\[ \begin{aligned} F_{Hitung} = \frac{JKG_{MG} - JKG_{MPTI}}{JKG_{MPTI}} . \frac{NT-N-K}{N-1} \end{aligned} \]
- Pengaruh Waktu
\[ \begin{aligned} F_{Hitung} = \frac{JKG_{MG} - JKG_{MPTW}}{JKG_{MPTW}} . \frac{N-T-K}{T-1} \end{aligned} \]
Uji Chow
\[ \begin{aligned} H_0 &: \text{Model Common Effect} \\ H_1 &: \text{Model Fixed Effect} \end{aligned} \]
\(H_0\) akan ditolak apabila
p-value lebih kecil dari alpha yaitu 0.05
pooltest(cem, fem)##
## F statistic
##
## data: KMIS ~ INFRA + PEKO + INFL + PGGR
## F = 173.29, df1 = 31, df2 = 380, p-value < 2.2e-16
## alternative hypothesis: unstabilityNilai p-value sebesar \(2.2e-16\) berarti keputusan Tolak
H0 atau model fixed effect lebih layak digunakan pada data ini.
Uji pengaruh individu / waktu
Selanjutnya karena yang terpilih adalah model Fixed-Effect akan dilihat komponen yang memiliki pengaruh tetap apakah individu, waktu, atau keduanya.
Efek Individu dan Waktu
plmtest(fem, type = "bp", effect = "twoways")##
## Lagrange Multiplier Test - two-ways effects (Breusch-Pagan)
##
## data: KMIS ~ INFRA + PEKO + INFL + PGGR
## chisq = 1746.8, df = 2, p-value < 2.2e-16
## alternative hypothesis: significant effectsEfek Individu
plmtest(fem, type = "bp", effect = "individual")##
## Lagrange Multiplier Test - (Breusch-Pagan)
##
## data: KMIS ~ INFRA + PEKO + INFL + PGGR
## chisq = 1726.2, df = 1, p-value < 2.2e-16
## alternative hypothesis: significant effectsEfek Waktu
plmtest(fem, type = "bp", effect = "time")##
## Lagrange Multiplier Test - time effects (Breusch-Pagan)
##
## data: KMIS ~ INFRA + PEKO + INFL + PGGR
## chisq = 20.619, df = 1, p-value = 5.605e-06
## alternative hypothesis: significant effectsDari ketiga uji diatas didapat nilai p-value yang lebih kecil dari 0.05 yang berarti keputusannya adalah Tak Tolak \(H_0\) atau terdapat pengaruh yang signifikan dari individu maupun waktu. Sehingga akan diduga parameter pengaruh twoways pada Random Effect Model.
Random-Effect Model
Model pengaruh acak mengasumsikan tidak ada korelasi antara pengaruh spesifik individu dan pengaruh spesifik waktu dengan peubah bebas. Asumsi ini membuat komponen sisaan dari pengaruh spesifik individu dan pengaruh spesifik waktu dimasukan kedalam sisaan. Model sisaan satu arah yaitu:
\[ \begin{aligned} y_{it} &= \alpha + \beta_{1}x_{1,it}+ \beta_{2}x_{2,it}+...+\beta_{k}x_{k,it}+ f_{i} +\varepsilon_{it}\ \text{(pengaruh Individu)} \\ \end{aligned} \]
\[ \begin{aligned} y_{it} &= \alpha + \beta_{1}x_{1,it}+ \beta_{2}x_{2,it}+...+\beta_{k}x_{k,it}+ \lambda_{t} +\varepsilon_{it} \\ \end{aligned} \]
Model sisaan dua arah yaitu:
\[ \begin{aligned} y_{it} &= \alpha + \beta_{1}x_{1,it}+ \beta_{2}x_{2,it}+...+\beta_{k}x_{k,it}+ f_{i}+ \lambda_{t} +\varepsilon_{it} \\ \end{aligned} \] Dengan \(\alpha\) merupakan intersep model dan \(f_i\) dan \(\lambda_t\) merupakan pengaruh acak masing-masing individu dan waktu. Pendugaan parameter model pengaruh acak menggunakan Generalized Least Square (Baltagi 2008).
Generalized Least Square
Analisis data panel tidak sepenuhnya mengikuti prosedur regresi biasa yang tidak dapat menjelaskan keragaman data. Perlunya mempertimbangkan keragaman data tersebut untuk dapat mengetahui dalam penggalian lebih lanjut terhadap informasi yang dikandung data. Salah satu cara untuk melakukan hal tersebut yaitu dengan metode Generalized Least Square (GLS). GLS dipakai sebagai estimator analisis data panel dengan persamaan sistem karena informasi heterogenitas antara individu dan waktu digunakan sebagai informasi yang penting untuk menghasilkan parameter \(\hat{\beta}\) Persamaan dalam pendugaan parameter GLS adalah sebagai berikut.
\[ \begin{aligned} \hat {\beta} &= (X'(\sum)^{-1} X )^{-1} X'(\sum) ^{-1} y \\ \end{aligned} \]
Metode GLS memasukkan struktur matriks ragam peragam residual pada parameter \(\hat{\beta}\) (Ekananda 2016).
Pendugaan Parameter
Jika pada model FEM mengasumsikan bahwa nilai dari pengaruh spesifik individu bersifat fixed pada setiap unit, maka pada model random effect diasumsikan bahwa pengaruh tersebut bersifat acak. Dalam hal ini misalnya model one way individu pada data Grundfeld, maka ke-10 firm memiliki nilai rataan umum intersep yang sama, sedangkan perbedaan intersep antar masing-masing unit firm direfleksikan pada error termnya.
rem_gls <- plm(KMIS ~ INFRA + PEKO + INFL + PGGR, data = kemis,
index = c("PROVINSI", "TAHUN"),
effect = "twoways", model = "random")
summary(rem_gls)## Twoways effects Random Effect Model
## (Swamy-Arora's transformation)
##
## Call:
## plm(formula = KMIS ~ INFRA + PEKO + INFL + PGGR, data = kemis,
## effect = "twoways", model = "random", index = c("PROVINSI",
## "TAHUN"))
##
## Balanced Panel: n = 32, T = 13, N = 416
##
## Effects:
## var std.dev share
## idiosyncratic 1.6093 1.2686 0.047
## individual 32.6514 5.7141 0.947
## time 0.2158 0.4646 0.006
## theta: 0.9385 (id) 0.5653 (time) 0.5648 (total)
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -3.48112 -0.80340 -0.21541 0.46391 6.46679
##
## Coefficients:
## Estimate Std. Error z-value Pr(>|z|)
## (Intercept) 1.0883e+01 1.2462e+00 8.7327 < 2.2e-16 ***
## INFRA -7.8316e-07 1.4391e-07 -5.4421 5.266e-08 ***
## PEKO 1.6972e-01 2.4964e-02 6.7988 1.055e-11 ***
## INFL 1.7252e-01 3.8357e-02 4.4977 6.871e-06 ***
## PGGR 1.6255e-01 7.7831e-02 2.0886 0.03675 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 1018
## Residual Sum of Squares: 773.87
## R-Squared: 0.23978
## Adj. R-Squared: 0.23238
## Chisq: 129.63 on 4 DF, p-value: < 2.22e-16Didapat model sebagai berikut:
\[ \begin{aligned} KMIS_{it} &= 10.883+\lambda_t-7.832(10^{-7})INFRA_{it}+0.1697PEKO_{it}+0.1725INFL_{it}+0.1625PGGR_{it}+f_{i}+ \lambda_{t}+\varepsilon_{it} \\ \end{aligned} \]
Interpretasi :
Setiap kenaikan satu satuan
InfrastrukturmenyebabkanTingkat Kemiskinanmenurun sebesar7.832x10^-7persen dengan menganggap peubah lain konstan.Setiap kenaikan satu satuan
Pertumbuhan EkonomimenyebabkanTingkat Kemiskinanmeningkat sebesar0.1697persen dengan menganggap peubah lain konstan.Setiap kenaikan satu satuan
InflasimenyebabkanTingkat Kemiskinanmeningkat sebesar0.1725persen dengan menganggap peubah lain konstan.Setiap kenaikan satu satuan
PengangguranmenyebabkanTingkat Kemiskinanmeningkat sebesar0.1625persen dengan menganggap peubah lain konstan.
Untuk semua p-value di bawah 5%, maka terdapat cukup bukti untuk menolak \(H_0\) (kurang dari 5% kemungkinan bahwa \(H_0\) benar).
- Model memiliki
p-valuepada uji F yang lebih kecil dari 0.05 (< 2.22e-16) yang berarti model ini layak untuk digunakan - p-value uji-T untuk semua peubah < 5%. Artinya, dapat dikatakan bahwa semua peubah berpengaruh signifikan terhadap model.
\(R^2\) merupakan nilai yang menunjukan kekuatan hubungan antara peubah bebas dan peubah respons. Semakin tinggi nilai \(R^2\), maka semakin baik model tersebut.
\(R^2\) pada model sebesar 0.23978 atau 23.978% yang mengindikasikan bahwa peubah yang digunakan tidak baik dalam menggambarkan model.
Model dapat menjelaskan sebesar 23.238% dari keragaman data yang dilihat dari \(R^2 Adj\) yang mengindikasikan bahwa peubah yang digunakan tidak baik dalam menggambarkan model.
Cek Signifikansi Pengaruh Individu
#efek individu
plmtest(rem_gls,type = "bp", effect="individu")##
## Lagrange Multiplier Test - (Breusch-Pagan)
##
## data: KMIS ~ INFRA + PEKO + INFL + PGGR
## chisq = 1726.2, df = 1, p-value < 2.2e-16
## alternative hypothesis: significant effectsKarena nilai p-value < 0.05, maka tolak \(H_0\) atau pengaruh individu signifikan terhadap model. ## Cek Signifikansi Pengaruh waktu
#efek waktu
plmtest(rem_gls,type = "bp", effect="time")##
## Lagrange Multiplier Test - time effects (Breusch-Pagan)
##
## data: KMIS ~ INFRA + PEKO + INFL + PGGR
## chisq = 20.619, df = 1, p-value = 5.605e-06
## alternative hypothesis: significant effectsKarena nilai p-value < 0.05, maka tolak \(H_0\) atau pengaruh waktu signifikan terhadap model. Artinya, model twoways tepat digunakan untuk metode REM.
tidy_ranef <- tidy(ranef(rem_gls, effect="individual"))
kable(tidy_ranef, digits=3, caption = "Pengaruh Acak Individu",
col.names = c("Firm", "Pengaruh Acak Individu"))| Firm | Pengaruh Acak Individu |
|---|---|
| Aceh | 6.126 |
| BABEL | -7.343 |
| Bali | -7.966 |
| Banten | -7.781 |
| Bengkulu | 4.284 |
| DIY | 1.665 |
| Gorontalo | 5.649 |
| JABAR | -2.633 |
| Jambi | -4.748 |
| JATENG | 2.839 |
| JATIM | 2.102 |
| KALBAR | -3.821 |
| KALSEL | -7.716 |
| KALTENG | -6.500 |
| KALTIM | -6.354 |
| Kepri | -6.480 |
| Lampung | 2.918 |
| Maluku | 7.982 |
| MALUT | -5.020 |
| NTB | 5.297 |
| NTT | 10.094 |
| PAPUA | 15.991 |
| PapuaBarat | 13.705 |
| RIAU | -4.408 |
| SULBAR | -0.234 |
| SULSEL | -1.920 |
| SULTENG | 2.308 |
| SulTenggara | 1.422 |
| SULUT | -4.619 |
| SUMBAR | -4.909 |
| SUMSEL | 1.670 |
| SUMUT | -1.601 |
Output diatas merupakan pengaruh acak dari setiap unit individu, nilai tersebut tersebut menunjukkan seberapa besar perbedaan nilai komponen error acak masing-masing unit indvidu terhadap nilai intersep umum.
tidy_ranef <- tidy(ranef(rem_gls, effect="time"))
kable(tidy_ranef, digits=3, caption = "Pengaruh Acak Waktu",
col.names = c("Firm", "Pengaruh Acak Waktu"))| Firm | Pengaruh Acak Waktu |
|---|---|
| 2008 | 1.605 |
| 2009 | 1.982 |
| 2010 | 0.593 |
| 2011 | 0.169 |
| 2012 | -0.496 |
| 2013 | -0.924 |
| 2014 | -1.205 |
| 2015 | -0.515 |
| 2016 | 0.530 |
| 2017 | -0.476 |
| 2018 | -0.639 |
| 2019 | -0.704 |
| 2020 | 0.080 |
Output diatas merupakan pengaruh acak dari setiap unit waktu, nilai tersebut tersebut menunjukkan seberapa besar perbedaan nilai komponen error acak masing-masing unit waktu terhadap nilai intersep umum.
FEM vs REM
Uji Haussman
Pada bagian sebelumnya, cukup bukti yang menyatakan bahwa model FEM
lebih baik daripada model CEM. Oleh karena itu, tahap selanjutnya akan
dibandingkan kebaikan model FEM jika dibanding model REM. Uji yang
digunakan adalah uji Hausman dengan menggunakan fungsi
phtest(FEM,REM).
Pengujian ini dilakukan untuk memilih model antara model pengaruh tetap atau model pengaruh acak yang sesuai untuk menggambarkan suatu data panel. Uji Hausman didasarkan pada perbedaan penduga model pengaruh tetap \(\hat{\beta}_{MPT}\) dengan penduga model pengaruh acak \(\hat{\beta}_{MPA}\). Kedua penduga konsisten dalam kondisi \(H_0\), tetapi \(\hat{\beta}_{MPA}\) akan bias dan tidak konsisten pada \(H_1\). Pengujian hipotesisnya adalah:
\(H_0\): Model pengaruh acak (REM) adalah model yang tepat
\(H_1\): Model pengaruh tetap (FEM) adalah model yang tepat
#fem dengan rem gls
phtest(fem_twoways, rem_gls)##
## Hausman Test
##
## data: KMIS ~ INFRA + PEKO + INFL + PGGR
## chisq = 40.108, df = 4, p-value = 4.111e-08
## alternative hypothesis: one model is inconsistentKarena nilai p-value yang diperoleh kurang dari dari alpha 5%, maka Tolak H0 yang mengindikasikan bahwa model yang sesuai adalah model FEM.
UJi Diagnostik Sisaan
res <- residuals(fem_twoways)Uji Normalitas
jarque.bera.test(res)##
## Jarque Bera Test
##
## data: res
## X-squared = 65.966, df = 2, p-value = 4.774e-15Histogram
# Histogram
ggplot(as.data.frame(res), aes(x = res)) +
geom_histogram(aes(y = after_stat(density)), color = "white", fill = "steelblue") +
geom_density(color = "red", linewidth = 1) +
theme_minimal()## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Uji Autokorelasi
pbgtest(fem_twoways)##
## Breusch-Godfrey/Wooldridge test for serial correlation in panel models
##
## data: KMIS ~ INFRA + PEKO + INFL + PGGR
## chisq = 194.35, df = 13, p-value < 2.2e-16
## alternative hypothesis: serial correlation in idiosyncratic errorsUji Kehomogenan
bptest(fem_twoways)##
## studentized Breusch-Pagan test
##
## data: fem_twoways
## BP = 16.998, df = 4, p-value = 0.001935Kesimpulan uji diagnostik :
- Sisaan tidak menyebar normal (Asumsi kenormalan tidak terpenuhi)
- Sisaan memiliki autokorelasi (Asumsi kebebasan belum terpenuhi)
- Sisaan tidak homogen (Asumsi kehomogenan belum terpenuhi)
Terlihat bahwa semua asumsi tidak terpenuhi, oleh karena itu, diperlukan penanganan yaitu dengan transformasi menggunakan First Difference. Yaitu sebagai berikut:
Penanganan Asumsi - FEM First Difference
First Difference Estimator (FD-Estimator) merupakan estimator untuk data panel yang digunakan pada model efek tetap satu arah (Oneway-FEM) atau dua arah (twoway-FEM). FD berguna dalam validasi model dimana pada data panel terdapat beberapa nilai amatan yang hilang. Berikut adalah model umum yang digunakan dalah FEM First Difference:
\[ \begin{aligned} \Delta y_{it} &= \beta_{1}\Delta x_{it1} + \beta_{2}\Delta x_{it2} + ... + \beta_{k}\Delta x_{itk} + \varepsilon_{it} \end{aligned} \]
Dalam hal ini, kita akan menggunakan metode FD pada model twoways. Untuk \(\tilde{y}_{it}= \Delta y_{it}\) dan \(\tilde{X}_{it}= \Delta X_{it}\), didefinisikan:
\(\tilde{{y} }_{i,t}={y} _{i,t}-{y} _{i-1,t}-{y} _{i,t-1}+{y} _{i-1,t-1}\)
dan
\(\tilde{{X} }_{s,t}={X} _{s,t}-{X} _{i-1,t}-{X} _{i,t-1}+{X} _{i-1,t-1}\)
Sehingga, kita perlu mentransformasi setiap nilai di peubah respons (Y = KMIS) dan peubah bebas (X) pada data dengan fungsi berikut:
transform_data <- function(df) {
df$dY <- NA # create a new column for the transformed Y values
df$dX1 <- NA # create a new column for the transformed X4 values
df$dX2 <- NA # create a new column for the transformed X2 values
df$dX3 <- NA # create a new column for the transformed X3 values
df$dX4 <- NA # create a new column for the transformed X4 values
df$PROVINSI <- as.numeric(df$PROVINSI)
df$TAHUN <- as.numeric(as.character(df$TAHUN))
for (i in 2:32) {
for (t in 2009:2020) {
y_it <- df$Y[df$PROVINSI == i & df$TAHUN == t]
y_i1t <- df$Y[df$PROVINSI == i-1 & df$TAHUN == t]
y_it1 <- df$Y[df$PROVINSI == i & df$TAHUN == t-1]
y_i1t1 <- df$Y[df$PROVINSI == i-1 & df$TAHUN == t-1]
df$dY[df$PROVINSI == i & df$TAHUN == t] <- y_it - y_i1t - y_it1 + y_i1t1
x1_it <- df$X1[df$PROVINSI == i & df$TAHUN == t]
x1_i1t <- df$X1[df$PROVINSI == i-1 & df$TAHUN == t]
x1_it1 <- df$X1[df$PROVINSI == i & df$TAHUN == t-1]
x1_i1t1 <- df$X1[df$PROVINSI == i-1 & df$TAHUN == t-1]
df$dX1[df$PROVINSI == i & df$TAHUN == t] <- x1_it - x1_i1t - x1_it1 + x1_i1t
x2_it <- df$X2[df$PROVINSI == i & df$TAHUN == t]
x2_i1t <- df$X2[df$PROVINSI == i-1 & df$TAHUN == t]
x2_it1 <- df$X2[df$PROVINSI == i & df$TAHUN == t-1]
x2_i1t1 <- df$X2[df$PROVINSI == i-1 & df$TAHUN == t-1]
df$dX2[df$PROVINSI == i & df$TAHUN == t] <- x2_it - x2_i1t - x2_it1 + x2_i1t
x3_it <- df$X3[df$PROVINSI == i & df$TAHUN == t]
x3_i1t <- df$X3[df$PROVINSI == i-1 & df$TAHUN == t]
x3_it1 <- df$X3[df$PROVINSI == i & df$TAHUN == t-1]
x3_i1t1 <- df$X3[df$PROVINSI == i-1 & df$TAHUN == t-1]
df$dX3[df$PROVINSI == i & df$TAHUN == t] <- x3_it - x3_i1t - x3_it1 + x3_i1t
x4_it <- df$X4[df$PROVINSI == i & df$TAHUN == t]
x4_i1t <- df$X4[df$PROVINSI == i-1 & df$TAHUN == t]
x4_it1 <- df$X4[df$PROVINSI == i & df$TAHUN == t-1]
x4_i1t1 <- df$X4[df$PROVINSI == i-1 & df$TAHUN == t-1]
df$dX4[df$PROVINSI == i & df$TAHUN == t] <- x4_it - x4_i1t - x4_it1 + x4_i1t
}
}
return(df)
}
datapanel <- kemis
colnames(datapanel)[3:7] <- c("Y","X1","X2","X3","X4")
datapanel_transform <- transform_data(datapanel)
datapanel_transform$PROVINSI <- as.character(datapanel$PROVINSI)
datapanel_transform$PROVINSI <- as.factor(datapanel_transform$PROVINSI)
datapanel_transform$TAHUN <- as.factor(datapanel_transform$TAHUN)
datapanel_transform <- datapanel_transform[-c(3:7)]Sehingga dihasilkan data baru:
datapanel_transformMerupakan hasil transformasi dengan metode first difference. Dapat dilihat juga:
datapanel_transform[is.na(datapanel_transform$dY),]Terlihat bahwa untuk setiap provinsi pada tahun 2008, dihasilkan missing value yang disebabkan oleh proses differencing. Selain itu, Individu (Provinsi) Aceh juga menghasilkan missing value untuk semua tahunnya, hal ini karena provinsi Aceh dijadikan sebagai referensi untuk proses differencing.
Selanjutnya, kita lakukan analisis FEM Two ways pada data yang telah ditransformasi
fem_fd <- plm(dY ~ dX1 + dX2 + dX3 + dX4, datapanel_transform, effect = "twoways",
model = "within", index = c("PROVINSI","TAHUN"))
summary(fem_fd)## Twoways effects Within Model
##
## Call:
## plm(formula = dY ~ dX1 + dX2 + dX3 + dX4, data = datapanel_transform,
## effect = "twoways", model = "within", index = c("PROVINSI",
## "TAHUN"))
##
## Balanced Panel: n = 31, T = 12, N = 372
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -3.3121975 -0.3403635 0.0093208 0.3378036 3.8711121
##
## Coefficients:
## Estimate Std. Error t-value Pr(>|t|)
## dX1 -2.4022e-08 1.0283e-07 -0.2336 0.8154
## dX2 2.0223e-02 1.7955e-02 1.1263 0.2609
## dX3 -1.4416e-02 2.4369e-02 -0.5916 0.5546
## dX4 -7.1072e-02 6.0126e-02 -1.1820 0.2380
##
## Total Sum of Squares: 244.66
## Residual Sum of Squares: 242.27
## R-Squared: 0.0097933
## Adj. R-Squared: -0.12689
## F-statistic: 0.806049 on 4 and 326 DF, p-value: 0.522Uji Diagnostik Sisaan ulang
Uji Normalitas
res2 <- residuals(fem_fd)
jarque.bera.test(res2)##
## Jarque Bera Test
##
## data: res2
## X-squared = 498.3, df = 2, p-value < 2.2e-16# Histogram
ggplot(as.data.frame(res2), aes(x = res2)) +
geom_histogram(aes(y = after_stat(density)), color = "white", fill = "steelblue") +
geom_density(color = "red", linewidth = 1) +
theme_minimal()## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Uji Autokorelasi
pbgtest(fem_fd)##
## Breusch-Godfrey/Wooldridge test for serial correlation in panel models
##
## data: dY ~ dX1 + dX2 + dX3 + dX4
## chisq = 83.889, df = 12, p-value = 7.439e-13
## alternative hypothesis: serial correlation in idiosyncratic errorsUji Kehomogenan
bptest(fem_fd)##
## studentized Breusch-Pagan test
##
## data: fem_fd
## BP = 7.0383, df = 4, p-value = 0.1339Dari hasil di atas, dapat terlihat bahwa p-value dari uji
bptest > 0.05 yang mengindikasikan bahwa asumsi
kehomogenan telah terpenuhi dari yang sebelumnya tidak. Artinya,
transformasi ini berhasil mengatasi heterogenitas pada sisaan.
Namun, perlu diperhatikan juga bahwa pada model setelah transformasi, tidak ada parameter yang signifikan dengan nilai p-value uji F sebesar 0.522 > 0.05 yang berarti model tidak layak. Selain itu, nilai R^2-nya sangat kecil yaitu sebesar 0.0097933. Artinya, dapat dikatakan bahwa model sebelum transformasi lebih baik dengan asumsi homogen yang tidak terpenuhi.
Salah satu solusi untuk masalah ini adalah untuk menggunakan transformasi lain. Jika model yang didapat oleh transformasi lain tidak lebih baik, maka solusi alternatifnya adalah dengan menggunakan analisis non-parametrik yang tidak memerlukan asumsi.