TUGAS RESPONSI STA1353 - ANALISIS DATA PANEL (KELOMPOK 4)

2023-04-07

NAMA	NIM
Syifa Maulany	G1401201003
Ervina Dwi Anggrahini	G1401201017
Tantri Gustina Dewi	G1401201031
Fidira Dwi Andari	G1401201059
Irma Wakhidatul Setyowati	G1401201089

Penggunaan Analisis Regresi Data Panel dalam Pemodelan Kemiskinan di Pulau Sumatra Tahun 2008-2020

Pendahuluan

Pulau Sumatera adalah salah satu pulau terbesar di Indonesia, dengan luas wilayah sekitar 473.481 km persegi. Meskipun memiliki sumber daya alam yang melimpah, pulau Sumatera masih menghadapi masalah kemiskinan yang signifikan. Menurut data Badan Pusat Statistik (BPS 2020), tingkat kemiskinan di Sumatera mencapai sekitar 9,41 persen. Angka ini menunjukkan bahwa sekitar 9,41 persen dari total populasi di Sumatera hidup di bawah garis kemiskinan, yang ditetapkan pada Rp 384.258 per kapita per bulan. Masalah kemiskinan di Sumatera dipengaruhi oleh berbagai faktor, seperti akses terbatas terhadap pendidikan, kesehatan, dan lapangan kerja yang layak. Selain itu, infrastruktur yang kurang memadai dan kebijakan pemerintah yang belum efektif juga berdampak pada meningkatnya tingkat kemiskinan di Sumatera. Tingginya tingkat kemiskinan di Sumatera juga berdampak pada masalah sosial, seperti tingginya angka putus sekolah, kekurangan gizi, dan kekurangan fasilitas kesehatan yang memadai. Oleh karena itu, pemerintah dan masyarakat perlu bekerja sama untuk mengatasi masalah kemiskinan di Sumatera dan meningkatkan kesejahteraan masyarakat secara keseluruhan (Cahyani 2021).

Library yang Digunakan

library(readxl)
library(plm)
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:plm':
## 
##     between, lag, lead

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(tidyverse)

## ── Attaching packages
## ───────────────────────────────────────
## tidyverse 1.3.2 ──

## ✔ ggplot2 3.4.0     ✔ purrr   0.3.5
## ✔ tibble  3.1.8     ✔ stringr 1.4.1
## ✔ tidyr   1.2.1     ✔ forcats 0.5.2
## ✔ readr   2.1.3     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::between() masks plm::between()
## ✖ dplyr::filter()  masks stats::filter()
## ✖ dplyr::lag()     masks plm::lag(), stats::lag()
## ✖ dplyr::lead()    masks plm::lead()

library(ggplot2)
library(performance)
library(corrplot)

## corrplot 0.92 loaded

library(gplots)

## 
## Attaching package: 'gplots'
## 
## The following object is masked from 'package:stats':
## 
##     lowess

Data

Data yang digunakan dalam penelitian ini merupakan data sekunder bertipe data panel mengenai kemiskinan di Indonesia tahun 2008-2020. Data ini terdiri dari satu peubah respon berupa KMIS serta empat peubah prediktor, meliputi INFRA, PEKO, INFL, dan PGGR. Unit pengamatan dalam data ini berupa 32 provinsi di Indonesia (DKI Jakarta, Kalimantan Utara, dan tiga provinsi baru Papua tidak dimasukkan dalam data ini).

Input Data

data<-read_excel("C:/SEMESTER 6/ADP/dataresponsiadp.xlsx")
str(data)

## tibble [416 × 7] (S3: tbl_df/tbl/data.frame)
##  $ PROVINSI: chr [1:416] "Bali" "Bali" "Bali" "Bali" ...
##  $ TAHUN   : num [1:416] 2008 2009 2010 2011 2012 ...
##  $ 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 ...
##  $ KMIS    : num [1:416] 6.17 5.13 4.88 4.2 3.95 4.49 4.76 5.25 4.15 4.14 ...

head(data)

## # A tibble: 6 × 7
##   PROVINSI TAHUN  INFRA  PEKO  INFL  PGGR  KMIS
##   <chr>    <dbl>  <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Bali      2008 430814  5.97  9.62  3.31  6.17
## 2 Bali      2009 484340  5.33  4.37  3.13  5.13
## 3 Bali      2010 336352  5.83  8.1   3.06  4.88
## 4 Bali      2011 298867  6.49  3.75  2.95  4.2 
## 5 Bali      2012 348529  6.65  4.71  2.1   3.95
## 6 Bali      2013 398608  6.69  8.16  1.83  4.49

unique(data$PROVINSI)

##  [1] "Bali"        "NTB"         "NTT"         "JABAR"       "Banten"     
##  [6] "JATENG"      "DIY"         "JATIM"       "KALBAR"      "KALTENG"    
## [11] "KALSEL"      "KALTIM"      "Maluku"      "MALUT"       "PAPUA"      
## [16] "PapuaBarat"  "SULUT"       "Gorontalo"   "SULTENG"     "SULSEL"     
## [21] "SULBAR"      "SulTenggara" "Aceh"        "SUMUT"       "SUMBAR"     
## [26] "RIAU"        "Kepri"       "Jambi"       "SUMSEL"      "BABEL"      
## [31] "Bengkulu"    "Lampung"

unique(data$TAHUN)

##  [1] 2008 2009 2010 2011 2012 2013 2014 2015 2016 2017 2018 2019 2020

data$TAHUN<-as.factor(data$TAHUN)

Split Data By Region

nusa<-data[c(1:39),]
jawa<-data[c(40:104),]
kalimantan<-data[c(105:156),]
maluku<-data[c(157:182),]
papua<-data[c(183:208),]
sulawesi<-data[c(209:286),]
sumatera<-data[c(287:416),]

Berdasarkan data yang tersedia, hanya dipilih Wilayah Pulau Sumatera dengan total 10 provinsi yang akan dilakukan analisis lebih lanjut karena Pulau Sumatera merupakan salah satu pulau terbsesar di Indonesia. Selain itu, kami tidak menggunakan data secara keseluruhan (Indonesia) karena terdapat beberapa provinsi yang datanya tidak teramati dalam data ini.

Membuat Data Panel

pdata <- pdata.frame(sumatera, index=c("PROVINSI","TAHUN"))
head(pdata)

##           PROVINSI TAHUN   INFRA  PEKO  INFL  PGGR  KMIS
## Aceh-2008     Aceh  2008 1005049 -5.27 11.92  9.56 23.53
## Aceh-2009     Aceh  2009 1040297 -5.58  3.72  8.71 21.80
## Aceh-2010     Aceh  2010  812758  2.74  5.86  8.37 20.98
## Aceh-2011     Aceh  2011  912601  5.02  3.43  9.00 19.57
## Aceh-2012     Aceh  2012  932017  5.20  0.22  9.06 18.58
## Aceh-2013     Aceh  2013 1176245  2.83  7.31 10.12 17.72

summary(pdata)

##      PROVINSI      TAHUN        INFRA              PEKO             INFL       
##  Aceh    :13   2008   :10   Min.   : 111046   Min.   :-5.580   Min.   : 0.220  
##  BABEL   :13   2009   :10   1st Qu.: 485142   1st Qu.: 4.173   1st Qu.: 2.518  
##  Bengkulu:13   2010   :10   Median : 909946   Median : 5.190   Median : 3.820  
##  Jambi   :13   2011   :10   Mean   :1107489   Mean   : 4.610   Mean   : 5.032  
##  Kepri   :13   2012   :10   3rd Qu.:1465496   3rd Qu.: 6.125   3rd Qu.: 7.513  
##  Lampung :13   2013   :10   Max.   :4729762   Max.   : 7.190   Max.   :18.400  
##  (Other) :52   (Other):70                                                      
##       PGGR             KMIS       
##  Min.   : 2.600   Min.   : 4.500  
##  1st Qu.: 4.640   1st Qu.: 7.473  
##  Median : 5.675   Median : 9.490  
##  Mean   : 5.953   Mean   :11.179  
##  3rd Qu.: 6.940   3rd Qu.:14.985  
##  Max.   :10.340   Max.   :23.530  
##

glimpse(pdata)

## Rows: 130
## Columns: 7
## $ PROVINSI <fct> Aceh, Aceh, Aceh, Aceh, Aceh, Aceh, Aceh, Aceh, Aceh, Aceh, A…
## $ TAHUN    <fct> 2008, 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2…
## $ INFRA    <pseries> 1005049, 1040297, 812758, 912601, 932017, 1176245, 118570…
## $ PEKO     <pseries> -5.27, -5.58, 2.74, 5.02, 5.20, 2.83, 1.65, -0.72, 3.31, …
## $ INFL     <pseries> 11.92, 3.72, 5.86, 3.43, 0.22, 7.31, 8.09, 1.53, 3.95, 4.…
## $ PGGR     <pseries> 9.56, 8.71, 8.37, 9.00, 9.06, 10.12, 9.02, 9.93, 7.57, 6.…
## $ KMIS     <pseries> 23.53, 21.80, 20.98, 19.57, 18.58, 17.72, 16.98, 17.11, 1…

Eksplorasi Data

Mengecek “Balance”

pdata %>%
  is.pbalanced()

## [1] TRUE

Mengecek Dimensi Waktu

pdata %>%
  is.pconsecutive()

##     Aceh    BABEL Bengkulu    Jambi    Kepri  Lampung     RIAU   SUMBAR 
##     TRUE     TRUE     TRUE     TRUE     TRUE     TRUE     TRUE     TRUE 
##   SUMSEL    SUMUT 
##     TRUE     TRUE

Mengecek Sebaran Data

Boxplot

boxplot(sumatera[,-c(1,2)])

Pada boxplot di atas, diketahui bahwa sebaran data sumatera tidak sama, karena peubah INFRA menggunakan satuan juta sedangkan peubah lainnya menggunakan satuan persen.

Histogram

hist(sumatera$INFRA)

hist(sumatera$PEKO)

hist(sumatera$INFL)

hist(sumatera$PGGR)

hist(sumatera$KMIS)

Dari histogram-histogram di atas, diketahui bahwa sebaran data peubah INFRA, INFL, dan KMIS cenderung menjulur ke kanan. Dapat diketahui pula bahwa sebaran data peubah PEKO menjulur ke kiri, sedangkan sebaran data peubah PGGR cenderung hampir simetris sehingga nilai modus, median, dan rataannya cenderung sama.

Plot Korelasi

corrplot(cor(sumatera[,-c(1,2)]), type = "upper")

corrplot(cor(sumatera[c(1:13),-c(1,2)]),type = "upper")

corrplot(cor(sumatera[c(117:130),-c(1,2)]),type = "upper")

Dari corrplot di atas, dapat diketahui bahwa korelasi dari data keseluruhan sumatera terlihat kecil. Namun, jika difokuskan pada provinsi tertentu, terdapat pasangan peubah dengan korelasi tinggi. Contohnya yaitu ketika dilihat korelasi antarpeubah pada provinsi Aceh, diketahui bahwa peubah-peubah yang memiliki korelasi yang cukup tinggi yaitu PGGR dengan INFRA, KMIS dengan INFRA, dan PEKO dengan KMIS. Contoh selanjutnya yaitu pada provinsi Lampung, diketahui peubah-peubah yang memiliki korelasi yang cukup tinggi yaitu INFRA dan PGGR, INFRA dan KMIS, serta PGGR dan KMIS.

Multiple Linechart

##peubah respon dibedakan berdasarkan individu
ggplot(sumatera, aes(x=TAHUN, y=KMIS, group= PROVINSI)) +
  geom_line(aes(color = PROVINSI), size = 1.15)+
  geom_point(size = 0.9)+
  theme_classic()

par(mfrow = c(3, 2))

##peubah respon dibedakan berdasarkan tahun
par(mfrow = c(3, 2))

ggplot(sumatera, aes(x=PROVINSI, y=KMIS, group= TAHUN)) +
  geom_line(aes(color = TAHUN), size = 1.15)+
  geom_point(size = 0.9)+
  theme_classic()

ggplot(sumatera, aes(x=TAHUN, y=INFRA, group= PROVINSI)) +
  geom_line(aes(color = PROVINSI), size = 1.15)+
  geom_point(size = 0.9)+
  theme_classic()

ggplot(sumatera, aes(x=TAHUN, y=PEKO, group= PROVINSI)) +
  geom_line(aes(color = PROVINSI), size = 1.15)+
  geom_point(size = 0.9)+
  theme_classic()

ggplot(sumatera, aes(x=TAHUN, y=INFL, group= PROVINSI)) +
  geom_line(aes(color = PROVINSI), size = 1.15)+
  geom_point(size = 0.9)+
  theme_classic()

ggplot(sumatera, aes(x=TAHUN, y=PGGR, group= PROVINSI)) +
  geom_line(aes(color = PROVINSI), size = 1.15)+
  geom_point(size = 0.9)+
  theme_classic()

##Y~X
par(mfrow = c(2, 2))

ggplot(sumatera, aes(x=INFRA, y = KMIS))+
  geom_point(aes(color= TAHUN))

ggplot(sumatera, aes(x=INFRA, y = KMIS))+
  geom_point(aes(color= PROVINSI))

ggplot(sumatera, aes(x=PEKO, y = KMIS))+
  geom_point(aes(color= TAHUN))

ggplot(sumatera, aes(x=PEKO, y = KMIS))+
  geom_point(aes(color= PROVINSI))

ggplot(sumatera, aes(x=INFL, y = KMIS))+
  geom_point(aes(color= TAHUN))

ggplot(sumatera, aes(x=INFL, y = KMIS))+
  geom_point(aes(color= PROVINSI))

ggplot(sumatera, aes(x=PGGR, y = KMIS))+
  geom_point(aes(color= TAHUN))

ggplot(sumatera, aes(x=PGGR, y = KMIS))+
  geom_point(aes(color= PROVINSI))

plotmeans(KMIS ~ PROVINSI, main="Keragaman Antar Individu",data)
plotmeans(KMIS ~ TAHUN, main="Keragaman Antar Waktu",data)

#terdapat keragaman baik pada individu, dan sepertinya tidak terdapat keragaman antar waktu
coplot( KMIS ~ TAHUN|PROVINSI, type="b", data=pdata)

Dari plot di atas, dapat diketahui bahwa terdapat keragaman yang cukup baik pada individu, tetapi tidak terlihat bahwa terdapat keragaman antarwaktu.

Common Effect Model (CEM)

Common Effect Model (CEM) merupakan salah satu model data panel yang menggabungkan data time series dan cross-section. Data gabungan ini dianggap sebagai suatu kesatuan pengamatan sehingga untuk mengestimasi parameter model ini dapat digunakan ordinary least square (OLS). CEM berfokus pada konsep “efek umum” atau “efek total”, yaitu bagaimana suatu faktor mempengaruhi hasil secara keseluruhan, tidak memperhitungkan interaksi antara faktor-faktor lain (Setiawan dan Kusrini 2010).

cem<-plm(scale(sumatera$KMIS)~INFRA+PEKO+INFL+PGGR, data = pdata, model = "pooling")
summary(cem)

## Pooling Model
## 
## Call:
## plm(formula = scale(sumatera$KMIS) ~ INFRA + PEKO + INFL + PGGR, 
##     data = pdata, model = "pooling")
## 
## Balanced Panel: n = 10, T = 13, N = 130
## 
## Residuals:
##      Min.   1st Qu.    Median   3rd Qu.      Max. 
## -2.276831 -0.461508  0.022039  0.482178  1.699348 
## 
## Coefficients:
##                Estimate  Std. Error t-value  Pr(>|t|)    
## (Intercept) -1.5386e+00  4.0645e-01 -3.7854 0.0002371 ***
## INFRA        2.4350e-07  9.3871e-08  2.5940 0.0106207 *  
## PEKO        -5.9160e-02  3.3516e-02 -1.7651 0.0799852 .  
## INFL         4.2735e-02  2.4328e-02  1.7567 0.0814246 .  
## PGGR         2.2284e-01  4.9075e-02  4.5407 1.301e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    129
## Residual Sum of Squares: 95.541
## R-Squared:      0.25937
## Adj. R-Squared: 0.23567
## F-statistic: 10.944 on 4 and 125 DF, p-value: 1.2183e-07

multicollinearity(cem)

## # Check for Multicollinearity
## 
## Low Correlation
## 
##   Term  VIF   VIF 95% CI Increased SE Tolerance Tolerance 95% CI
##  INFRA 1.10 [1.01, 1.72]         1.05      0.91     [0.58, 0.99]
##   PEKO 1.24 [1.09, 1.65]         1.11      0.81     [0.61, 0.92]
##   INFL 1.16 [1.04, 1.61]         1.08      0.86     [0.62, 0.96]
##   PGGR 1.19 [1.06, 1.62]         1.09      0.84     [0.62, 0.95]

model_performance(cem)

## # Indices of model performance
## 
## AIC     |    AICc |     BIC |    R2 | R2 (adj.) |  RMSE | Sigma
## ---------------------------------------------------------------
## 340.886 | 341.569 | 358.092 | 0.259 |     0.236 | 0.857 | 0.874

res.cem <- residuals(cem)

Berdasarkan output di atas, dapat dilihat bahwa pada model CEM seluruh peubah memiliki nilai VIF < 10, sehingga tidak terdeteksi adanya multikolinearitas antar peubah-peubahnya. Sementara itu, nilai \(R^2\) _adjective_ yang diperoleh pada model ini terbilang cukup kecil, yaitu hanya sebesar 23.6% dan RMSE sebesar 0.857.

Diagnostik Sisaan

##normality (H0: normal)
library(tseries)

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

jarque.bera.test(res.cem)

## 
##  Jarque Bera Test
## 
## data:  res.cem
## X-squared = 1.0689, df = 2, p-value = 0.586

Berdasarkan Jarque Bera test, didapatkan p-value > α = 5% maka terima H0. Artinya sisaan model menyebar normal.

##histogram
hist(res.cem, 
     xlab = "sisaan",
     col = "#27D3D3", 
     breaks=30,  
     prob = TRUE) 
lines(density(res.cem), # density plot
      lwd = 2, # thickness of line
      col = "chocolate3")

#plotqqnorm
set.seed(369)
res.cem1 <- as.numeric(res.cem)
qqnorm(res.cem1,datax=T, col="blue")
qqline(rnorm(length(res.cem1),mean(res.cem1),sd(res.cem1)),datax=T, col="red")

Secara eksploratif, dapat dilihat bahwa sisaan menyebar normal. Hal ini sesuai dengan hasil uji Jarque Bera.

##autocorelastion (H0: sisaan saling bebas)
pbgtest(cem)

## 
##  Breusch-Godfrey/Wooldridge test for serial correlation in panel models
## 
## data:  scale(sumatera$KMIS) ~ INFRA + PEKO + INFL + PGGR
## chisq = 95.502, df = 13, p-value = 1.228e-14
## alternative hypothesis: serial correlation in idiosyncratic errors

Hasil Breusch-Godfrey/Wooldridge test menunjukkan bahwa p-value < α = 5% maka tolak H0. Artinya, terdapat autokorelasi atau sisaan tidak saling bebas.

##homogen (H0: ragam homogen)
library(lmtest)

## Loading required package: zoo

## 
## Attaching package: 'zoo'

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

bptest(cem)

## 
##  studentized Breusch-Pagan test
## 
## data:  cem
## BP = 21.619, df = 4, p-value = 0.0002386

Hasil Breusch-Pagan test menunjukkan bahwa p-value < α = 5% maka tolak H0. Artinya, ragam sisaan tidak homogen.

Fixed Effect Model (FEM)

Fixed Effect Model mempertimbangkan adanya pengaruh individu dan waktu. Model ini memiliki asumsi bahwa intersep antar individu dan waktu berbeda, namun koefisien regresi konstan untuk semua individu dan waktu. Jika ada pengaruh waktu atau individu dalam fixed effect model, maka model disebut model sisaan satu arah. Sementara itu, jika ada pengaruh waktu dan individu, maka model disebut model sisaan dua arah.

#individual effect
fem.ind<-plm(scale(sumatera$KMIS)~INFRA+PEKO+INFL+PGGR, data = pdata, model = "within",
             effect = "individual", index = c("PROVINSI, TAHUN"))
summary(fem.ind)

## Oneway (individual) effect Within Model
## 
## Call:
## plm(formula = scale(sumatera$KMIS) ~ INFRA + PEKO + INFL + PGGR, 
##     data = pdata, effect = "individual", model = "within", index = c("PROVINSI, TAHUN"))
## 
## Balanced Panel: n = 10, T = 13, N = 130
## 
## Residuals:
##     Min.  1st Qu.   Median  3rd Qu.     Max. 
## -0.50039 -0.13958 -0.04659  0.15065  0.82360 
## 
## Coefficients:
##          Estimate  Std. Error t-value  Pr(>|t|)    
## INFRA -1.7377e-07  4.2393e-08 -4.0990 7.724e-05 ***
## PEKO   1.3401e-02  1.1278e-02  1.1883  0.237134    
## INFL   2.1814e-02  7.4449e-03  2.9301  0.004081 ** 
## PGGR   1.1817e-01  2.4681e-02  4.7881 5.014e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    15.68
## Residual Sum of Squares: 7.5958
## R-Squared:      0.51556
## Adj. R-Squared: 0.46127
## F-statistic: 30.8627 on 4 and 116 DF, p-value: < 2.22e-16

Pendugaan dengan model FEM Individu memberikan nilai \(R^2\) sebesar 0.51556 atau 51.56% dari keragaman data yang dapat dijelaskan oleh model. Hal ini menunjukkan model FEM Individu tidak terlalu baik digunakan untuk data persentase kemiskinan di Pulau Sumatra. Didapatkan bahwa nilai p value < α = 0.05. Artinya, peubah penjelas yang digunakan memiliki pengaruh yang signifikan terhadap persentase kemiskinan pada taraf nyata 5%. Persamaan model dapat dinyatakan sebagai berikut:

\(KMIS = -0.00000017 INFRA + 0.0134 PEKO + 0.0218 INFL + 0.1182 PGGR\)

Setiap pertambahan infrastruktur sebesar 1 juta unit, persentase kemiskinan akan menurun sebesar 0.0000017 persen dengan asumsi peubah lain tetap. Peubah lainnya berhubungan positif sehingga peningkatan pada masing-masing peubah pertumbuhan ekonomi, inflasi, dan pengangguran akan meningkatkan kemiskinan sebesar 0.0134 persen, 0.0218 persen, dan 0.1182 persen dengan asumsi peubah lainnya tetap.

summary(fixef(fem.ind, effect = "individual"))

##          Estimate Std. Error t-value  Pr(>|t|)    
## Aceh      0.69066    0.26189  2.6372  0.009505 ** 
## BABEL    -0.82136    0.15908 -5.1633 1.017e-06 ***
## Bengkulu -1.27415    0.16005 -7.9608 1.284e-12 ***
## Jambi    -1.25542    0.18040 -6.9592 2.174e-10 ***
## Kepri    -1.87565    0.21573 -8.6943 2.667e-14 ***
## Lampung  -1.17164    0.19937 -5.8766 4.098e-08 ***
## RIAU     -0.16670    0.20557 -0.8109  0.419082    
## SUMBAR   -1.93660    0.22866 -8.4693 8.826e-14 ***
## SUMSEL    0.66754    0.20440  3.2659  0.001435 ** 
## SUMUT     0.31749    0.25948  1.2236  0.223589    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

res.fem.ind <- residuals(fem.ind)

Berdasarkan hasil di atas, terlihat bahwa masing-masing provinsi di Pulau Sumatra memiliki pengaruh tetap yang berbeda-beda terhadap model.

Diagnostik Sisaan

##normality (H0: normal)
library(tseries)
jarque.bera.test(res.fem.ind)

## 
##  Jarque Bera Test
## 
## data:  res.fem.ind
## X-squared = 18.907, df = 2, p-value = 7.842e-05

Berdasarkan Jarque Bera test, didapatkan p-value < α = 5% maka tolak H0. Artinya sisaan model tidak menyebar normal.

##histogram
hist(res.fem.ind, 
     xlab = "sisaan",
     col = "#27D3D3", 
     breaks=30,  
     prob = TRUE) 
lines(density(res.fem.ind), # density plot
      lwd = 2, # thickness of line
      col = "chocolate3")

##plotqqnorm
set.seed(369)
res.fem.ind1 <- as.numeric(res.fem.ind)
qqnorm(res.fem.ind1,datax=T, col="blue")
qqline(rnorm(length(res.fem.ind1),mean(res.fem.ind1),sd(res.fem.ind1)),datax=T, col="red")

Secara eksploratif, dapat dilihat bahwa sisaan tidak menyebar normal. Hal ini sesuai dengan hasil uji Jarque Bera.

##autocorelastion (H0: sisaan saling bebas)
pbgtest(fem.ind)

## 
##  Breusch-Godfrey/Wooldridge test for serial correlation in panel models
## 
## data:  scale(sumatera$KMIS) ~ INFRA + PEKO + INFL + PGGR
## chisq = 58.621, df = 13, p-value = 9.258e-08
## alternative hypothesis: serial correlation in idiosyncratic errors

Hasil Breusch-Godfrey/Wooldridge test menunjukkan bahwa p-value < α = 5% maka tolak H0. Artinya, terdapat autokorelasi atau sisaan tidak saling bebas.

##homogen (H0: ragam homogen)
library(lmtest)
bptest(fem.ind)

## 
##  studentized Breusch-Pagan test
## 
## data:  fem.ind
## BP = 21.619, df = 4, p-value = 0.0002386

Hasil Breusch-Pagan test menunjukkan bahwa p-value < α = 5% maka tolak H0. Artinya, ragam sisaan tidak homogen.

FEM dengan Pengaruh Spesifik Individu dan atau Waktu

#Pengaruh spesifik Individu dan atau waktu

#efek individu
plmtest(fem.ind,type = "bp", effect="individual" )

## 
##  Lagrange Multiplier Test - (Breusch-Pagan)
## 
## data:  scale(sumatera$KMIS) ~ INFRA + PEKO + INFL + PGGR
## chisq = 497.99, df = 1, p-value < 2.2e-16
## alternative hypothesis: significant effects

#efek waktu
plmtest(fem.ind,type = "bp", effect="time" )

## 
##  Lagrange Multiplier Test - time effects (Breusch-Pagan)
## 
## data:  scale(sumatera$KMIS) ~ INFRA + PEKO + INFL + PGGR
## chisq = 0.00032526, df = 1, p-value = 0.9856
## alternative hypothesis: significant effects

#efek individu dan waktu
plmtest(fem.ind,type = "bp", effect="twoways" )

## 
##  Lagrange Multiplier Test - two-ways effects (Breusch-Pagan)
## 
## data:  scale(sumatera$KMIS) ~ INFRA + PEKO + INFL + PGGR
## chisq = 497.99, df = 2, p-value < 2.2e-16
## alternative hypothesis: significant effects

Uji Lagrange Multiplier (LM) Breusch-Pagan menguji signifikansi pengaruh individu dan/atau waktu pada model. Hasil uji menunjukkan bahwa pengaruh individu dan pengaruh individu dan waktu masing-masing signifikan pada taraf nyata 5%. Sedangkan, pengaruh waktu tidak signifikan pada taraf nyata 5%.

Uji Chouw

#Uji Chouw
#memilih antara CEM dan FEM
pooltest(cem,fem.ind)

## 
##  F statistic
## 
## data:  scale(sumatera$KMIS) ~ INFRA + PEKO + INFL + PGGR
## F = 149.23, df1 = 9, df2 = 116, p-value < 2.2e-16
## alternative hypothesis: unstability

#better goes with FEM

Diperoleh p-value < α = 5%, maka tolak H0. Disimpulkan bahwa Fixed Effect Model Pengaruh Individu lebih baik digunakan daripada Common Effect Model pada taraf nyata 5%.

Uji Kebebasan Pada Unit Cross Section

#Uji kebebasan pada unit cross section
pcdtest(fem.ind, test = c("lm"), index=NULL,w =NULL )

## 
##  Breusch-Pagan LM test for cross-sectional dependence in panels
## 
## data:  scale(sumatera$KMIS) ~ INFRA + PEKO + INFL + PGGR
## chisq = 164.55, df = 45, p-value = 1.582e-15
## alternative hypothesis: cross-sectional dependence

pcdtest(fem.ind, test = c("cd"), index=NULL,w =NULL )

## 
##  Pesaran CD test for cross-sectional dependence in panels
## 
## data:  scale(sumatera$KMIS) ~ INFRA + PEKO + INFL + PGGR
## z = 9.017, p-value < 2.2e-16
## alternative hypothesis: cross-sectional dependence

Berdasarkan hasil uji Breusch-Pagan didapatkan p-value < α = 5%, maka tolak H0. Artinya, terdapat dependensi antar unit individu atau antar unit individu tidak saling bebas.

Random Effect Model (REM)

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 dimasukkan kedalam sisaan.

#rem dengan metode gls, efek individu
rem_gls_ind <- plm(scale(sumatera$KMIS)~INFRA+PEKO+INFL+PGGR, data = pdata, 
                   index = c("PROVINSI, TAHUN"), 
                   effect = "individual", model = "random")
summary(rem_gls_ind)

## Oneway (individual) effect Random Effect Model 
##    (Swamy-Arora's transformation)
## 
## Call:
## plm(formula = scale(sumatera$KMIS) ~ INFRA + PEKO + INFL + PGGR, 
##     data = pdata, effect = "individual", model = "random", index = c("PROVINSI, TAHUN"))
## 
## Balanced Panel: n = 10, T = 13, N = 130
## 
## Effects:
##                   var std.dev share
## idiosyncratic 0.06548 0.25589 0.157
## individual    0.35184 0.59316 0.843
## theta: 0.8812
## 
## Residuals:
##      Min.   1st Qu.    Median   3rd Qu.      Max. 
## -0.485532 -0.178232 -0.010677  0.113255  0.972466 
## 
## Coefficients:
##                Estimate  Std. Error z-value  Pr(>|z|)    
## (Intercept) -7.5724e-01  2.8373e-01 -2.6689 0.0076106 ** 
## INFRA       -1.5487e-07  4.4321e-08 -3.4942 0.0004754 ***
## PEKO         1.2265e-02  1.1936e-02  1.0276 0.3041553    
## INFL         2.2233e-02  7.8915e-03  2.8173 0.0048426 ** 
## PGGR         1.2772e-01  2.5703e-02  4.9692 6.722e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    17.279
## Residual Sum of Squares: 9.2112
## R-Squared:      0.46691
## Adj. R-Squared: 0.44985
## Chisq: 109.482 on 4 DF, p-value: < 2.22e-16

ranef(rem_gls_ind, effect = "individual")

##       Aceh      BABEL   Bengkulu      Jambi      Kepri    Lampung       RIAU 
##  1.3176670 -0.1127668 -0.5576312 -0.5412018 -1.1742615 -0.4798091  0.5055431 
##     SUMBAR     SUMSEL      SUMUT 
## -1.2464764  1.3318781  0.9570588

Nilai intersep pada model pengaruh acak sebesar -7.5724e-01, nilai ini menggambarkan rataan umum intersep model. Koefisien INFRA sebesar -1.5487e-07; PEKO sebesar 1.2265e-02; INFL sebesar 2.2233e-02;PGGR sebesar 1.2772e-01 dengan nilai R-Squared sebesar 46.7%. Hal ini menunjukkan bahwa model REM dengan efek individu tidak terlalu baik digunakan untuk data persentase kemiskinan di Pulau Sumatra. Dapat dilihat pula berdasarkan nilai p-value, dengan taraf signifikansi 5%, hanya terdapat satu peubah penjelas yang signifikan yaitu peubah PEKO.

Pemeriksaan Efek dalam Model

Selanjutnya dapat melakukan pengecekan pada model pengaruh acak yang dibangun, apakah terdapat pengaruh satu arah atau dua arah. Pada prinsipnya pengecekan pengaruh sama halnya seperti yang dilakukan pada model pengaruh tetap dengan fungsi plmtest().

#pemeriksaan efek dalam model
##efek individu
plmtest(rem_gls_ind,type = "bp", effect="individual")

## 
##  Lagrange Multiplier Test - (Breusch-Pagan)
## 
## data:  scale(sumatera$KMIS) ~ INFRA + PEKO + INFL + PGGR
## chisq = 497.99, df = 1, p-value < 2.2e-16
## alternative hypothesis: significant effects

##efek waktu
plmtest(rem_gls_ind,type = "bp", effect="time")

## 
##  Lagrange Multiplier Test - time effects (Breusch-Pagan)
## 
## data:  scale(sumatera$KMIS) ~ INFRA + PEKO + INFL + PGGR
## chisq = 0.00032526, df = 1, p-value = 0.9856
## alternative hypothesis: significant effects

##efek individu dan waktu
plmtest(rem_gls_ind,type = "bp", effect="twoways")

## 
##  Lagrange Multiplier Test - two-ways effects (Breusch-Pagan)
## 
## data:  scale(sumatera$KMIS) ~ INFRA + PEKO + INFL + PGGR
## chisq = 497.99, df = 2, p-value < 2.2e-16
## alternative hypothesis: significant effects

Berdasarkan hasil pengujian efek, maka disimpulkan bahwa:

Pengaruh Individu: Tolak H0
Pengaruh Waktu: Terima H0
Pengaruh Two Way: Tolak H0

sehingga model pengaruh acak pengaruh individu yang sesuai.

Diagnostik Sisaan Model REM

##diagnostik sisaan model rem ind
res.rem_ind <- residuals(rem_gls_ind)

#normality
library(tseries)
jarque.bera.test(res.rem_ind)

## 
##  Jarque Bera Test
## 
## data:  res.rem_ind
## X-squared = 45.371, df = 2, p-value = 1.405e-10

Berdasarkan Jarque Bera Test, diperoleh nilai p-value < taraf signifikansi 5%, sehingga dapat diketahui bahwa sisaan model tidak menyebar normal.

#kolmogorov smirnov
ks.test(res.rem_ind, "pnorm")

## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  res.rem_ind
## D = 0.31724, p-value = 8.646e-12
## alternative hypothesis: two-sided

Berdasarkan uji Kolmogorov-Smirnov didapatkan p-value < taraf signifikansi 5%, sehingga Tolak H0 yang artinya sisaan model tidak berdistribusi normal.

#histogram
hist(res.rem_ind, 
     xlab = "sisaan",
     col = "#27D3D3", 
     breaks=30,  
     prob = TRUE) 
lines(density(res.rem_ind), # density plot
      lwd = 2, # thickness of line
      col = "chocolate3")

#plotqqnorm
set.seed(1353)
res.rem_ind1 <- as.numeric(res.rem_ind)
qqnorm(res.rem_ind1,datax=T, col="blue")
qqline(rnorm(length(res.rem_ind1),mean(res.rem_ind1),sd(res.rem_ind1)),datax=T, col="red")

Secara eksploratif, dapat dilihat bahwa sisaan tidak menyebar normal.

#autocorelastion
adf.test(res.rem_ind, k=2) #alternatif : Terdapat autokorelasi

## 
##  Augmented Dickey-Fuller Test
## 
## data:  res.rem_ind
## Dickey-Fuller = -5.5371, Lag order = 2, p-value = 0.01
## alternative hypothesis: stationary

Hasil uji Dicky-Fuller menunjukkan nilai p-value < taraf signifikansi 5%, sehingga dapat disimpulkan bahwa terdapat autokorelasi.

#homogen
bptest(rem_gls_ind)

## 
##  studentized Breusch-Pagan test
## 
## data:  rem_gls_ind
## BP = 21.619, df = 4, p-value = 0.0002386

Hasil uji Breusch Pagan menunjukkan nilai p-value < taraf signifikansi 5%, sehingga dapat disimpulkan bahwa ragam sisaan tidak homogen.

Pemilihan Model Terbaik antara CEM dan FEM

Hipotesis yang digunakan adalah:

H0 : Random Effect Model

H1 : Fixed Effect Model

phtest(fem.ind, rem_gls_ind)#pvalue >5%, artinya model yang sesuai adalah rem_gls_ind

## 
##  Hausman Test
## 
## data:  scale(sumatera$KMIS) ~ INFRA + PEKO + INFL + PGGR
## chisq = 3.5313, df = 4, p-value = 0.4731
## alternative hypothesis: one model is inconsistent

Dari hausman Test di atas, didapatkan nilai p-value sebesar 0.47 yang mana lebih besar dari 5%. Artinya, model yang sesuai adalah model rem_gls_ind.

Akan diperiksa antara model pengaruh acak dengan common effect model dengan uji pengganda laggrang.

#memeriksa antara model pengaruh acak dengan common effect model dengan uji pengganda laggrang.
cem <- plm(scale(sumatera$KMIS)~INFRA+PEKO+INFL+PGGR, data = pdata, model="pooling")
c.ind <- plmtest(cem, effect = "individual", type = c("bp"))
c.ind

## 
##  Lagrange Multiplier Test - (Breusch-Pagan)
## 
## data:  scale(sumatera$KMIS) ~ INFRA + PEKO + INFL + PGGR
## chisq = 497.99, df = 1, p-value < 2.2e-16
## alternative hypothesis: significant effects

Setelah diperiksa, dapat diketahui bahwa pooling model tidak cukup. Hal ini terjadi karena terapat efek panel.

Daftar Pustaka

Cahyani RD. 2021. Analisis Kemiskinan di Pulau Sumatra Tahun 2017-2020[skripsi]. Surakarta: Universitas Muhammadiyah Surakarta.

Setiawan, Kusrini. 2010. Ekonometrika. Yogyakarta: Penerbit ANDI.