P02_STA1353_Common Effect Model

Data

library(readxl)
dataa <- read_excel("C:/Users/LENOVO/Documents/DIAH/datapanel.xlsx") 
knitr::kable(head(dataa))

dplyr::glimpse(dataa)

## Rows: 1,480
## Columns: 8
## $ KabKota  <chr> "Jembrana", "Jembrana", "Jembrana", "Jembrana", "Jembrana", "…
## $ Tahun    <dbl> 2014, 2015, 2016, 2017, 2018, 2019, 2020, 2021, 2014, 2015, 2…
## $ JPM      <dbl> 15800, 15830, 14530, 14780, 14350, 13550, 12600, 14240, 24400…
## $ PDRB     <dbl> 9020, 10198, 11168, 12169, 13206, 14141, 13465, 13510, 15066,…
## $ RLS      <dbl> 7.30, 7.54, 7.59, 7.62, 7.95, 8.22, 8.23, 8.35, 7.91, 8.07, 8…
## $ PPD      <dbl> 10944, 11168, 11343, 11468, 11666, 11902, 11790, 11675, 13492…
## $ HLS      <dbl> 11.48, 11.88, 12.27, 12.40, 12.61, 12.63, 12.65, 12.92, 12.04…
## $ Provinsi <chr> "Bali", "Bali", "Bali", "Bali", "Bali", "Bali", "Bali", "Bali…

summary(dataa)

##    KabKota              Tahun           JPM              PDRB       
##  Length:1480        Min.   :2014   Min.   :  3550   Min.   :   134  
##  Class :character   1st Qu.:2016   1st Qu.: 14138   1st Qu.:  2417  
##  Mode  :character   Median :2018   Median : 24295   Median :  5087  
##                     Mean   :2018   Mean   : 31015   Mean   :  9116  
##                     3rd Qu.:2019   3rd Qu.: 37365   3rd Qu.: 11448  
##                     Max.   :2021   Max.   :222190   Max.   :190318  
##       RLS              PPD             HLS          Provinsi        
##  Min.   : 0.630   Min.   : 3607   Min.   : 2.16   Length:1480       
##  1st Qu.: 6.860   1st Qu.: 7053   1st Qu.:11.89   Class :character  
##  Median : 7.740   Median : 8358   Median :12.57   Mode  :character  
##  Mean   : 7.668   Mean   : 8783   Mean   :12.37                     
##  3rd Qu.: 8.652   3rd Qu.:10279   3rd Qu.:13.22                     
##  Max.   :12.510   Max.   :19992   Max.   :16.89

Data terdiri dari 1480 unit pengamatan pada level kabupaten/kota dengan series amatan dari 2014 hingga 2021.

Eksplorasi Data

Jumlah Penduduk Miskin di wilayah timur sepanjang tahun 2014 - 2018

library(dplyr)

## 
## Attaching package: 'dplyr'

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

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

dataprovtimur <- dataa %>% 
  dplyr::select(c(Provinsi, Tahun, JPM))%>%
  group_by(Tahun) %>%
 summarise(across(where(is.numeric),sum))

dataprovtimur

## # A tibble: 8 × 2
##   Tahun     JPM
##   <dbl>   <dbl>
## 1  2014 5909092
## 2  2015 5919737
## 3  2016 5785240
## 4  2017 5764370
## 5  2018 5648650
## 6  2019 5596240
## 7  2020 5564170
## 8  2021 5714440

library(ggplot2)
library(hrbrthemes)

## NOTE: Either Arial Narrow or Roboto Condensed fonts are required to use these themes.

##       Please use hrbrthemes::import_roboto_condensed() to install Roboto Condensed and

##       if Arial Narrow is not on your system, please see https://bit.ly/arialnarrow

# Plot
dataprovtimur %>%
  tail(10) %>%
  ggplot( aes(x=Tahun, y=JPM)) +
    geom_line( color="Blue") +
    geom_point(shape=21, color="black", fill="#69b3a2", size=3) +
  scale_x_continuous(limits = c(2014, 2021), breaks = seq(2014, 2021, 1))+
    theme_ipsum() +
  theme(plot.title=element_text(size=15))+
    labs(title="Jumlah Penduduk Miskin di Wilayah Timur", subtitle =" Tahun 2014 - 2021")

Berdasarkan grafik di atas dapat dapat dilihat bahwa jumlah penduduk miskin pada wilayah timur cenderung mengalami penurunan pada rentang waktu 2014 hingga 2020. Namun jumlah penduduk miskin terlihat meningkat di tahun 2021 yang bersamaan dengan kejadian pandemi Covid 19.

Selanjutnya akan dilihat pola dari jumlah penduduk miskin setiap tahun pada setiap unit individu pada level provinsi.

dataProv <- dataa %>% 
  dplyr::select(c(Provinsi, Tahun, JPM))%>%
  group_by(Provinsi,Tahun) %>%
 summarise(across(where(is.numeric),sum))

## `summarise()` has grouped output by 'Provinsi'. You can override using the
## `.groups` argument.

knitr::kable(dataProv)

dataProv %>% 
  mutate(provinsi = forcats::fct_reorder2(Provinsi, Tahun, JPM)) %>%
ggplot(mapping = aes(y = JPM, x = Tahun ,group = provinsi, colour=provinsi)) +
   geom_line(size=1.2) +
  geom_point(shape=21, size=2) +
  scale_x_continuous(limits = c(2014, 2021), breaks = seq(2014, 2021, 1))+
  theme(plot.title=element_text(size=12))+
    labs(title="Jumlah Penduduk Miskin Per Provinsi", subtitle =" Tahun 2014-2021")+
  theme_minimal()

Berdasarkan grafik di atas dapat dilihat jumlah penduduk miskin di setiap provinsi pada wilayah timur cenderung mengalami penurunan pada rentang waktu 2014 hingga 2020. Namun jumlah penduduk miskin terlihat meningkat di tahun 2021 yang bersamaan dengan kejadian pandemi Covid 19. Gambaran ini sejalan dengan pola yang terjadi pada plot pola pergerakan jumlah penduduk miskin sebelumnya. Provinsi NTT, Papua, Sulawesi Selatan, dan NTB merupakan empat provinsi yang memiliki jumlah penduduk miskin tertinggi pada kurun waktu tersebut. Sedangkan provinsi Maluku Utara, Sulawesi Barat dan Bali merupakan provinsi dengan jumlah penduduk miskin yang lebih rendah dibandingkan sembilan provinsi lainnya di wilayah timur.

Untuk melihat rataan jumlah penduduk antar provinsi pada kurun waktu tersebut dapat diamati melalui grafik dibawah ini.

dataProv1 <- dataProv %>% 
  dplyr::select(-Tahun) %>%
  group_by(Provinsi) %>%
 summarise(across(where(is.numeric),mean)) %>%
  arrange(JPM)
dataProv1

## # A tibble: 13 × 2
##    Provinsi               JPM
##    <chr>                <dbl>
##  1 Maluku Utara        85908.
##  2 Sulawesi Barat     157736.
##  3 Bali               181726.
##  4 Gorontalo          195782.
##  5 Sulawesi Utara     197711.
##  6 Papua Barat        227839 
##  7 Maluku             320172.
##  8 Sulawesi Tenggara  347554 
##  9 Sulawesi Tengah    418094 
## 10 NTB                771582.
## 11 Sulawesi Selatan   793312.
## 12 Papua              905200 
## 13 NTT               1135124.

dataProv1 %>%
  mutate(Prov = forcats::fct_reorder(Provinsi, JPM)) %>%
ggplot(aes(x = Prov, y = JPM)) +
 geom_bar(stat="identity", fill="#6D67E4", alpha=.6, width=.6) +
  geom_text(aes(label =round(JPM)), hjust= -0.1, size=3)+
  ylim(c(0, 1200000))+
  coord_flip()+
  theme_minimal()+ 
   labs(title = 'Rata-Rata Jumlah Penduduk Miskin Per Provinsi')

Diagnostik Sisaan

res.ols <- residuals(ols)
yhat.ols <- fitted(ols)


#normality
library(tseries)

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

jarque.bera.test(res.ols)

## 
##  Jarque Bera Test
## 
## data:  res.ols
## X-squared = 6600.7, df = 2, p-value < 2.2e-16

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

#plotqqnorm
qqnorm(res.ols,datax=T, col="blue")
qqline(rnorm(length(res.ols),mean(res.ols),sd(res.ols)),datax=T, col="red")

library(lmtest)

## Loading required package: zoo

## 
## Attaching package: 'zoo'

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

dwtest(ols)

## 
##  Durbin-Watson test
## 
## data:  ols
## DW = 0.23289, p-value < 2.2e-16
## alternative hypothesis: true autocorrelation is greater than 0

bptest(ols)

## 
##  studentized Breusch-Pagan test
## 
## data:  ols
## BP = 134.57, df = 4, p-value < 2.2e-16

#Hasil plot menunjukkan ragam tidak homogen
plot(yhat.ols,res.ols, 
     xlab = "fitted value", 
     ylab = "residuals")

Diagnostik Sisaan

res.cem <- residuals(common)
yhat.cem <- fitted(common)

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

## 
##  Jarque Bera Test
## 
## data:  res.cem
## X-squared = 6600.7, df = 2, p-value < 2.2e-16

shapiro.test(res.cem)

## 
##  Shapiro-Wilk normality test
## 
## data:  res.cem
## W = 0.85831, p-value < 2.2e-16

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

#plotqqnorm
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")

pbgtest(common)

## 
##  Breusch-Godfrey/Wooldridge test for serial correlation in panel models
## 
## data:  JPM ~ PDRB + RLS + PPD + HLS
## chisq = 1153.6, df = 8, p-value < 2.2e-16
## alternative hypothesis: serial correlation in idiosyncratic errors

bptest(common)

## 
##  studentized Breusch-Pagan test
## 
## data:  common
## BP = 134.57, df = 4, p-value < 2.2e-16

Transformasi Logaritma Natural pada Seluruh Peubah

common1 <- plm(log(JPM)~log(PDRB)+log(RLS)+log(PPD)+log(HLS), data=dataa, model="pooling")

summary(common1)

## Pooling Model
## 
## Call:
## plm(formula = log(JPM) ~ log(PDRB) + log(RLS) + log(PPD) + log(HLS), 
##     data = dataa, model = "pooling")
## 
## Balanced Panel: n = 185, T = 8, N = 1480
## 
## Residuals:
##      Min.   1st Qu.    Median   3rd Qu.      Max. 
## -2.386016 -0.418585 -0.018865  0.453462  1.710436 
## 
## Coefficients:
##              Estimate Std. Error  t-value  Pr(>|t|)    
## (Intercept) 11.610960   0.716536  16.2043 < 2.2e-16 ***
## log(PDRB)    0.476109   0.022684  20.9888 < 2.2e-16 ***
## log(RLS)    -1.567425   0.104153 -15.0493 < 2.2e-16 ***
## log(PPD)    -0.755473   0.095486  -7.9119 4.938e-15 ***
## log(HLS)     1.731967   0.166756  10.3862 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    848.55
## Residual Sum of Squares: 561.59
## R-Squared:      0.33818
## Adj. R-Squared: 0.33639
## F-statistic: 188.428 on 4 and 1475 DF, p-value: < 2.22e-16

multicollinearity(common1)

## # Check for Multicollinearity
## 
## Low Correlation
## 
##       Term  VIF   VIF 95% CI Increased SE Tolerance Tolerance 95% CI
##  log(PDRB) 2.28 [2.11, 2.47]         1.51      0.44     [0.40, 0.47]
##   log(RLS) 4.72 [4.32, 5.18]         2.17      0.21     [0.19, 0.23]
##   log(PPD) 3.20 [2.94, 3.49]         1.79      0.31     [0.29, 0.34]
##   log(HLS) 4.13 [3.78, 4.52]         2.03      0.24     [0.22, 0.26]

model_performance(common1,metrics = "all")

## # Indices of model performance
## 
## AIC       |      AICc |       BIC |    R2 | R2 (adj.) |  RMSE | Sigma
## ---------------------------------------------------------------------
## 32545.875 | 32545.932 | 32577.674 | 0.338 |     0.336 | 0.616 | 0.617

res.cem <- residuals(common1)
yhat.cem <- fitted(common1)


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

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

shapiro.test(res.cem)

## 
##  Shapiro-Wilk normality test
## 
## data:  res.cem
## W = 0.99439, p-value = 2.263e-05

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

#plotqqnorm
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")

pbgtest(common1)

## 
##  Breusch-Godfrey/Wooldridge test for serial correlation in panel models
## 
## data:  log(JPM) ~ log(PDRB) + log(RLS) + log(PPD) + log(HLS)
## chisq = 991.42, df = 8, p-value < 2.2e-16
## alternative hypothesis: serial correlation in idiosyncratic errors

bptest(common1)

## 
##  studentized Breusch-Pagan test
## 
## data:  common1
## BP = 121.77, df = 4, p-value < 2.2e-16

Diagnostik Sisaan

res.gls <- residuals(common.gls)

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

## 
##  Jarque Bera Test
## 
## data:  res.gls
## X-squared = 4.5408, df = 2, p-value = 0.1033

shapiro.test(res.gls)

## 
##  Shapiro-Wilk normality test
## 
## data:  res.gls
## W = 0.99564, p-value = 0.0002855

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

#plotqqnorm
res.gls1 <- as.numeric(res.gls)
qqnorm(res.gls1,datax=T, col="blue")
qqline(rnorm(length(res.gls1),mean(res.gls1),sd(res.gls1)),datax=T, col="red")

bptest(common.gls)

## 
##  studentized Breusch-Pagan test
## 
## data:  common.gls
## BP = 121.77, df = 4, p-value < 2.2e-16

Informasi Tambahan (silakan dieksplorasi):

Dalam mengestimasi parameter model CEM, terdapat 4 metode estimasi yang dapat digunakan (Srihardianti et al. 2016), yaitu

• Ordinary Least Square (OLS), jika bersifat homoskedastik dan tidak ada cross-sectional correlation;

• Weighted Least Square (WLS), jika bersifat heteroskedastik dan tidak ada cross-sectional correlation;

• Seemingly Uncorrelated Regression (SUR), jika bersifat heteroskedastik dan ada cross-sectional correlation; dan

• Feasible Generalized Least Square (FGLS) dengan proses autoregressive (AR), jika bersifat heteroskedastik dan ada korelasi antar waktu pada residualnya. Pada regresi data panel, estimator yang digunakan adalah FGLS yang memiliki 2 tahap estimasi yang berbeda dengan OLS. Sehingga pengujian asumsi klasik, seperti homoskedastisitas dan non-autokorelasi, tidak dapat digunakan dalam estimator ini karena tidak lagi relevan dengan konsep pengujian asumsi klasik pada regresi OLS. Namun, pengujian multikolinieritas dan normalitas akan tetap dilakukan.

Referensi

1. Gujarati DN. 2004. Basic Econometrics 4th ed. New York:McGraw-Hill.
2. Srihardianti M, Mustafid, Prahutama A. 2016. Metode regresi data panel untuk peramalan konsumsi energi di Indonesia. Jurnal Gaussian. 5(3):475-485.
3. Pustaka lain yang relevan