PENDAHULUAN
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. ### 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} \\ \end{aligned} \]
atau
\[ \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} \]
Pendugaan parameter model pengaruh acak menggunakan Generalized Least Square (Baltagi 2008).
Pendugaan Parameter
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'(\Sigma)-1 \times X)-1\times X'(\Sigma)-1\times y \\ \end{aligned} \]
Metode GLS memasukkan struktur matriks ragam peragam residual pada parameter \(\hat\beta\) (Ekananda 2016).
Spesifikasi Model
Spesifikasi model dilakukan untuk memilih model yang sesuai di antara pendekatan dugaan model yang digunakan antara model pengaruh tetap atau model pengaruh acak. Berikut beberapa uji dalam menentukan spesifikasi model.
Uji Haussman
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 hipotesis yang digunakan yaitu:
\(H_{0}\) = Model pengaruh acak adalah model yang tepat.
\(H_{1}\) = Model pengaruh tetap adalah model yang tepat.
Statistik uji Hausmann yaitu: \[ \begin{aligned} X^2_{hitung} &= (\hat\beta MPA-\hat\beta MPT)'(VMPT-VMPA)-1(\hat\beta MPA-\hat\beta MPT) \\ \end{aligned} \] Keputusan tolak \(H_{0}\) jika \(χ^2_{hitung}\) > \(χ^2_{k,α}\) atau jika nilai p-value < α (Baltagi,2005).
PACKAGES
library(readxl)
library(performance)
library(lmtest)## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericlibrary(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, unionlibrary(kableExtra) #untuk tampilan tabel##
## Attaching package: 'kableExtra'## The following object is masked from 'package:dplyr':
##
## group_rowslibrary(RMySQL)## Loading required package: DBIlibrary(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats 1.0.0 ✔ readr 2.1.4
## ✔ ggplot2 3.4.1 ✔ stringr 1.5.0
## ✔ lubridate 1.9.2 ✔ tibble 3.2.0
## ✔ purrr 1.0.1 ✔ tidyr 1.3.0## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::between() masks plm::between()
## ✖ dplyr::filter() masks stats::filter()
## ✖ kableExtra::group_rows() masks dplyr::group_rows()
## ✖ dplyr::lag() masks plm::lag(), stats::lag()
## ✖ dplyr::lead() masks plm::lead()
## ℹ Use the ]8;;http://conflicted.r-lib.org/conflicted package]8;; to force all conflicts to become errorslibrary(quantmod)## Loading required package: xts
##
## ################################### WARNING ###################################
## # We noticed you have dplyr installed. The dplyr lag() function breaks how #
## # base R's lag() function is supposed to work, which breaks lag(my_xts). #
## # #
## # Calls to lag(my_xts) that you enter or source() into this session won't #
## # work correctly. #
## # #
## # All package code is unaffected because it is protected by the R namespace #
## # mechanism. #
## # #
## # Set `options(xts.warn_dplyr_breaks_lag = FALSE)` to suppress this warning. #
## # #
## # You can use stats::lag() to make sure you're not using dplyr::lag(), or you #
## # can add conflictRules('dplyr', exclude = 'lag') to your .Rprofile to stop #
## # dplyr from breaking base R's lag() function. #
## ################################### WARNING ###################################
##
## Attaching package: 'xts'
##
## The following objects are masked from 'package:dplyr':
##
## first, last
##
## Loading required package: TTR
## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zooDATA
Inisiasi Data
data_adp <- read_excel("C:/Users/HP/Downloads/Kelompok 6 ADP/dataadptahun.xlsx")
head(data_adp)## # A tibble: 6 × 6
## Provinsi Tahun IPM PDRB `Rata-rata_lama_sekolah` Jml_Penduduk_Miskin
## <chr> <dbl> <chr> <chr> <chr> <chr>
## 1 ACEH 2019 71.9 4.14 9.18 814.60
## 2 ACEH 2020 71.99 -0.37 9.33 824.41
## 3 ACEH 2021 72.18 2.79 9.37 842.25
## 4 BALI 2019 75.38 5.6 8.84 160.38
## 5 BALI 2020 75.5 -9.33 8.95 181.06
## 6 BALI 2021 75.69 -2.47 9.06 206.72data_adp$IPM <- as.numeric(data_adp$IPM)
data_adp$PDRB <- as.numeric(data_adp$PDRB)
data_adp$`Rata-rata_lama_sekolah` <- as.numeric(data_adp$`Rata-rata_lama_sekolah`)
data_adp$Jml_Penduduk_Miskin<- as.numeric(data_adp$Jml_Penduduk_Miskin)
str(data_adp)## tibble [102 × 6] (S3: tbl_df/tbl/data.frame)
## $ Provinsi : chr [1:102] "ACEH" "ACEH" "ACEH" "BALI" ...
## $ Tahun : num [1:102] 2019 2020 2021 2019 2020 ...
## $ IPM : num [1:102] 71.9 72 72.2 75.4 75.5 ...
## $ PDRB : num [1:102] 4.14 -0.37 2.79 5.6 -9.33 -2.47 5.26 -3.39 4.44 4.94 ...
## $ Rata-rata_lama_sekolah: num [1:102] 9.18 9.33 9.37 8.84 8.95 9.06 8.74 8.89 8.93 8.73 ...
## $ Jml_Penduduk_Miskin : num [1:102] 815 824 842 160 181 ...Data Panel
datapanel <- pdata.frame(data_adp)
pdim(datapanel)## Balanced Panel: n = 34, T = 3, N = 102Check Balancing Data
datapanel %>%
is.pbalanced()## [1] TRUECek Dimensi Waktu
datapanel %>%
is.pconsecutive()## ACEH BALI BANTEN
## TRUE TRUE TRUE
## BENGKULU DI YOGYAKARTA DKI JAKARTA
## TRUE TRUE TRUE
## GORONTALO JAMBI JAWA BARAT
## TRUE TRUE TRUE
## JAWA TENGAH JAWA TIMUR KALIMANTAN BARAT
## TRUE TRUE TRUE
## KALIMANTAN SELATAN KALIMANTAN TENGAH KALIMANTAN TIMUR
## TRUE TRUE TRUE
## KALIMANTAN UTARA KEP. BANGKA BELITUNG KEP. RIAU
## TRUE TRUE TRUE
## LAMPUNG MALUKU MALUKU UTARA
## TRUE TRUE TRUE
## NUSA TENGGARA BARAT NUSA TENGGARA TIMUR PAPUA
## TRUE TRUE TRUE
## PAPUA BARAT RIAU SULAWESI BARAT
## TRUE TRUE TRUE
## SULAWESI SELATAN SULAWESI TENGAH SULAWESI TENGGARA
## TRUE TRUE TRUE
## SULAWESI UTARA SUMATERA BARAT SUMATERA SELATAN
## TRUE TRUE TRUE
## SUMATERA UTARA
## TRUEEksplorasi Data
Grafik Data
ggplot(data = datapanel, aes(x = Tahun, y = Jml_Penduduk_Miskin)) +
geom_line() +
labs(x = "Tahun", y = "Jumlah Penduduk Miskin") +
theme(legend.position = "none")+ theme_bw()
Scatter Plot IPM Vs Jumlah Penduduk Miskin
ggplot() + geom_point(aes(datapanel$IPM, datapanel$Jml_Penduduk_Miskin),colour = 'red') + ggtitle('IPM Vs Jumlah Penduduk Miskin') +
xlab('IPM') + ylab('Jumlah Penduduk Miskin')
Scatter Plot Rata-Rata Lama Sekolah Vs Jumlah Penduduk Miskin
ggplot() + geom_point(aes(datapanel$Rata.rata_lama_sekolah, datapanel$Jml_Penduduk_Miskin),colour = 'red') + ggtitle('Rata-Rata Lama Sekolah Vs Jumlah Penduduk Miskin') +
xlab('Rata-Rata Lama Sekolah') + ylab('Jumlah Penduduk Miskin')
Scatter Plot PDRB Vs Jumlah Penduduk Miskin
ggplot() + geom_point(aes(datapanel$PDRB, datapanel$Jml_Penduduk_Miskin),colour = 'red') + ggtitle('PDRB Vs Jumlah Penduduk Miskin') +
xlab('PDRB') + ylab('Jumlah Penduduk Miskin')
Boxplot Semua Peubah
data_adp$Provinsi <- as.factor(data_adp$Provinsi)
I <- ggplot(data_adp, aes(x = Provinsi, y = IPM, color=Provinsi)) +
geom_boxplot() + ggtitle("IPM") + theme(legend.position="none")
P <- ggplot(data_adp, aes(x = Provinsi, y = PDRB, color=Provinsi)) +
geom_boxplot() + ggtitle("PDRB") + theme(legend.position="none")
R <- ggplot(data_adp, aes(x = Provinsi, y = `Rata-rata_lama_sekolah`, color=Provinsi)) +
geom_boxplot() + ggtitle("Rata-Rata Lama Sekolah") + theme(legend.position="none")
J <- ggplot(data_adp, aes(x = Provinsi, y = Jml_Penduduk_Miskin, color=Provinsi)) +
geom_boxplot() + ggtitle("Jumlah Penduduk Miskin") + theme(legend.position="none")
library(gridExtra)##
## Attaching package: 'gridExtra'## The following object is masked from 'package:dplyr':
##
## combinegrid.arrange(I,P,R,J, ncol=2)
Plot Keragaman Antar Individu
library(gplots)##
## Attaching package: 'gplots'## The following object is masked from 'package:stats':
##
## lowessplotmeans(Jml_Penduduk_Miskin ~ `Rata-rata_lama_sekolah`, main="Keragaman Antar Individu",data_adp)## Warning in qt((1 + p)/2, ns - 1): NaNs produced## Warning in arrows(x, li, x, pmax(y - gap, li), col = barcol, lwd = lwd, :
## zero-length arrow is of indeterminate angle and so skipped
## Warning in arrows(x, li, x, pmax(y - gap, li), col = barcol, lwd = lwd, :
## zero-length arrow is of indeterminate angle and so skipped
## Warning in arrows(x, li, x, pmax(y - gap, li), col = barcol, lwd = lwd, :
## zero-length arrow is of indeterminate angle and so skipped## Warning in arrows(x, ui, x, pmin(y + gap, ui), col = barcol, lwd = lwd, :
## zero-length arrow is of indeterminate angle and so skipped
## Warning in arrows(x, ui, x, pmin(y + gap, ui), col = barcol, lwd = lwd, :
## zero-length arrow is of indeterminate angle and so skipped
## Warning in arrows(x, ui, x, pmin(y + gap, ui), col = barcol, lwd = lwd, :
## zero-length arrow is of indeterminate angle and so skipped
Plot Keragaman Antar Waktu
plotmeans(Jml_Penduduk_Miskin ~ Tahun, main="Keragaman Antar Waktu",data_adp)
Plot Keragaman Antar Waktu Per-Provinsi
coplot(Jml_Penduduk_Miskin ~ Tahun|Provinsi, type="b", data=datapanel)
FEM
fem <- plm(Jml_Penduduk_Miskin ~ IPM+PDRB+Rata.rata_lama_sekolah, datapanel, model = "within",effect= "individual",index = c("Provinsi","Tahun"))
summary(fem)## Oneway (individual) effect Within Model
##
## Call:
## plm(formula = Jml_Penduduk_Miskin ~ IPM + PDRB + Rata.rata_lama_sekolah,
## data = datapanel, effect = "individual", model = "within",
## index = c("Provinsi", "Tahun"))
##
## Balanced Panel: n = 34, T = 3, N = 102
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -388.2785 -22.6733 -7.5712 34.2204 217.8140
##
## Coefficients:
## Estimate Std. Error t-value Pr(>|t|)
## IPM -61.7904 71.0678 -0.8695 0.387796
## PDRB -1.8985 2.1848 -0.8690 0.388071
## Rata.rata_lama_sekolah 541.6641 169.0449 3.2043 0.002099 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 530900
## Residual Sum of Squares: 401060
## R-Squared: 0.24457
## Adj. R-Squared: -0.17382
## F-statistic: 7.01462 on 3 and 65 DF, p-value: 0.00036953Pendugaan dengan model FEM memberikan nilai R^2-Adj sebesar -0.17382. Nilai dugaan koefisien IPM Sebesar -61.7904, PDRB sebesar -1.8985, dan Rata-rata lama sekolah sebesar 541.6641. Nilai dari p-value < α (5%) menunjukkan bahwa x memiliki pengaruh yang signifikan terhadap Jml_Penduduk_Miskin pada taraf 5%.
Melihat Fixed Effect dari FEM , Nilai Intersep yang Berbeda-beda pada Masing-masing Individu
fixef(fem, type="level")## ACEH BALI BANTEN
## 247.7234 -2.4805 465.3698
## BENGKULU DI YOGYAKARTA DKI JAKARTA
## -54.7248 269.3254 -570.3106
## GORONTALO JAMBI JAWA BARAT
## 210.5492 75.2327 3702.3728
## JAWA TENGAH JAWA TIMUR KALIMANTAN BARAT
## 4228.3132 4576.1882 562.9377
## KALIMANTAN SELATAN KALIMANTAN TENGAH KALIMANTAN TIMUR
## 102.9734 -113.9621 -325.9614
## KALIMANTAN UTARA KEP. BANGKA BELITUNG KEP. RIAU
## -440.6052 135.8215 -658.0862
## LAMPUNG MALUKU MALUKU UTARA
## 1025.8168 -757.0660 -553.1488
## NUSA TENGGARA BARAT NUSA TENGGARA TIMUR PAPUA
## 992.3077 1059.0001 1037.8321
## PAPUA BARAT RIAU SULAWESI BARAT
## 128.1536 57.8264 -15.3178
## SULAWESI SELATAN SULAWESI TENGAH SULAWESI TENGGARA
## 696.1671 -61.3563 -158.4145
## SULAWESI UTARA SUMATERA BARAT SUMATERA SELATAN
## -440.3159 -37.1591 968.7935
## SUMATERA UTARA
## 584.5630Diagnostik Sisaan
res.fem <- residuals(fem)Normality
library(tseries)
jarque.bera.test(res.fem)##
## Jarque Bera Test
##
## data: res.fem
## X-squared = 981.51, df = 2, p-value < 2.2e-16Histogram
hist(res.fem,
xlab = "Sisaan",
col = "#27D3D3",
breaks=30,
prob = TRUE)
lines(density(res.fem), # density plot
lwd = 2, # thickness of line
col = "chocolate3")
Plot QQ-Norm
set.seed(1353)
res.fem1 <- as.numeric(res.fem)
qqnorm(res.fem1,datax=T, col="blue")
qqline(rnorm(length(res.fem1),mean(res.fem1),sd(res.fem1)),datax=T, col="red")
Auto-Corelation
pbgtest(fem) #alternatif : Terdapat autokorelasi##
## Breusch-Godfrey/Wooldridge test for serial correlation in panel models
##
## data: Jml_Penduduk_Miskin ~ IPM + PDRB + Rata.rata_lama_sekolah
## chisq = 30.244, df = 3, p-value = 1.226e-06
## alternative hypothesis: serial correlation in idiosyncratic errorsHomogen
library(lmtest)
bptest(fem)##
## studentized Breusch-Pagan test
##
## data: fem
## BP = 19.783, df = 3, p-value = 0.0001882Kesimpulan :
Kenormalan sisaan Belum Terpenuhi
Kebebasan sisaan Belum Terpenuhi
Kehomogenan ragam Belum Terpenuhi
GLS : One Way Random Effect Model
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 Jumlah Penduduk Miskin, maka ke-34 Provinsi memiliki nilai rataan umum intersep yang sama, sedangkan perbedaan intersep antar masing-masing unit Provinsi direfleksikan pada error term-nya.
One Way Individual
model <- `Jml_Penduduk_Miskin` ~ `IPM` + `PDRB` + `Rata.rata_lama_sekolah`
rem_gls_ind <- plm(model, data = datapanel,
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 = model, data = datapanel, effect = "individual",
## model = "random", index = c("Provinsi", "Tahun"))
##
## Balanced Panel: n = 34, T = 3, N = 102
##
## Effects:
## var std.dev share
## idiosyncratic 6.170e+03 7.855e+01 0.006
## individual 1.039e+06 1.019e+03 0.994
## theta: 0.9556
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -275.291 -36.312 -15.820 11.984 372.314
##
## Coefficients:
## Estimate Std. Error z-value Pr(>|z|)
## (Intercept) -646.3765 2828.0785 -0.2286 0.81921
## IPM -13.6978 51.4274 -0.2664 0.78997
## PDRB -3.0735 2.2875 -1.3436 0.17908
## Rata.rata_lama_sekolah 278.2680 150.2211 1.8524 0.06397 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 762090
## Residual Sum of Squares: 691390
## R-Squared: 0.092764
## Adj. R-Squared: 0.064992
## Chisq: 10.0204 on 3 DF, p-value: 0.018393Selanjutnya output dirapihkan dengan menggunakan library
stargazer,kableExtra, dan
broom.
library(stargazer)##
## Please cite as:## Hlavac, Marek (2022). stargazer: Well-Formatted Regression and Summary Statistics Tables.## R package version 5.2.3. https://CRAN.R-project.org/package=stargazerstargazer(rem_gls_ind, type='text')##
## ==================================================
## Dependent variable:
## ---------------------------
## Jml_Penduduk_Miskin
## --------------------------------------------------
## IPM -13.698
## (51.427)
##
## PDRB -3.073
## (2.287)
##
## Rata.rata_lama_sekolah 278.268*
## (150.221)
##
## Constant -646.376
## (2,828.078)
##
## --------------------------------------------------
## Observations 102
## R2 0.093
## Adjusted R2 0.065
## F Statistic 10.020**
## ==================================================
## Note: *p<0.1; **p<0.05; ***p<0.01library(kableExtra)library(broom)
kable(tidy(rem_gls_ind), digits=3,
caption="GLS_Random Effect Model_Individu")| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | -646.376 | 2828.078 | -0.229 | 0.819 |
| IPM | -13.698 | 51.427 | -0.266 | 0.790 |
| PDRB | -3.073 | 2.287 | -1.344 | 0.179 |
| Rata.rata_lama_sekolah | 278.268 | 150.221 | 1.852 | 0.064 |
Pengecekan Pengaruh Spesifik Model
Selanjutnya, dilakukan pemeriksaan mengenai model random effect yang dibangun, apakah terdapat pengaruh efek individu, waktu, atau keduanya.
Efek Individu dan Waktu
plmtest(rem_gls_ind, type = "bp", effect = "twoways" )##
## Lagrange Multiplier Test - two-ways effects (Breusch-Pagan)
##
## data: model
## chisq = 99.864, df = 2, p-value < 2.2e-16
## alternative hypothesis: significant effectsEfek Individu
plmtest(rem_gls_ind,type = "bp", effect = "individu" )##
## Lagrange Multiplier Test - (Breusch-Pagan)
##
## data: model
## chisq = 98.776, df = 1, p-value < 2.2e-16
## alternative hypothesis: significant effectsEfek Waktu
plmtest(rem_gls_ind,type = "bp", effect = "time" )##
## Lagrange Multiplier Test - time effects (Breusch-Pagan)
##
## data: model
## chisq = 1.0883, df = 1, p-value = 0.2968
## alternative hypothesis: significant effectsBerdasarkan pengujian pengaruh, didapatkan kesimpulan bahwa \(H_{0}\) ditolak pada model dengan pengaruh two-way dan pengaruh individu. Oleh karena itu, model random effect yang dibangun memiliki pengaruh individu.
Maka dilihatlah pengaruh individunya sebagai berikut :
tidy_ranef_ind <- tidy(ranef(rem_gls_ind, effect="individual"))## Warning: 'tidy.numeric' is deprecated.
## See help("Deprecated")kable(tidy_ranef_ind, digits=3,caption = "Pengaruh Acak Individu",
col.names = c("Provinsi", "Pengaruh Acak Individu"))| Provinsi | Pengaruh Acak Individu |
|---|---|
| ACEH | -119.059 |
| BALI | -632.002 |
| BANTEN | -42.244 |
| BENGKULU | -517.346 |
| DI YOGYAKARTA | -421.653 |
| DKI JAKARTA | -880.222 |
| GORONTALO | -388.581 |
| JAMBI | -460.361 |
| JAWA BARAT | 3114.722 |
| JAWA TENGAH | 3428.180 |
| JAWA TIMUR | 3806.751 |
| KALIMANTAN BARAT | -101.980 |
| KALIMANTAN SELATAN | -480.539 |
| KALIMANTAN TENGAH | -621.185 |
| KALIMANTAN TIMUR | -785.735 |
| KALIMANTAN UTARA | -827.901 |
| KEP. BANGKA BELITUNG | -534.644 |
| KEP. RIAU | -985.355 |
| LAMPUNG | 432.377 |
| MALUKU | -837.159 |
| MALUKU UTARA | -813.776 |
| NUSA TENGGARA BARAT | 281.401 |
| NUSA TENGGARA TIMUR | 577.657 |
| PAPUA | 532.581 |
| PAPUA BARAT | -355.380 |
| RIAU | -396.018 |
| SULAWESI BARAT | -472.792 |
| SULAWESI SELATAN | 90.417 |
| SULAWESI TENGAH | -428.045 |
| SULAWESI TENGGARA | -565.011 |
| SULAWESI UTARA | -797.421 |
| SUMATERA BARAT | -503.794 |
| SUMATERA SELATAN | 417.502 |
| SUMATERA UTARA | 286.614 |
Output diatas menggambarkan pengaruh acak dari setiap unit individu dimana masing-masing nilai menunjukkan seberapa besar perbedaan nilai komponen error acak masing-masing unit indvidu terhadap nilai intersep umum.
Diagnostik Model One Way REM
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 = 406.66, df = 2, p-value < 2.2e-16Kolmogorov-Smirnov
ks.test(res.rem_ind, "pnorm")##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: res.rem_ind
## D = 0.59419, p-value < 2.2e-16
## alternative hypothesis: two-sidedHistogram
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")
Plot QQ-Norm
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")
Auto-Corelation
adf.test(res.rem_ind, k=2) #alternatif : Terdapat autokorelasi##
## Augmented Dickey-Fuller Test
##
## data: res.rem_ind
## Dickey-Fuller = -3.0967, Lag order = 2, p-value = 0.122
## alternative hypothesis: stationaryHomogen
bptest(rem_gls_ind)##
## studentized Breusch-Pagan test
##
## data: rem_gls_ind
## BP = 19.783, df = 3, p-value = 0.0001882Output 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.
Misalnya, nilai 320.978552 pada Provinsi 1 menunjukkan besarnya perbedaan komponen error acak Provinsi 1 terhadap intersep umum. Artinya Provinsi 1 memiliki intersep 320.978552 yang lebih tinggi dari intersep umum. Begitu pula interpretasi untuk lainnya.
Pada package R, pendugaan model REM disediakan beberapa pilihan metode random yang dapat digunakan yaitu :
walhus: Wallace and Hussain (1969)swar: Swamy Arora (1972)amemiya1: Amemiya (1971)ht: Hausman and Taylor (1981)nerlove: Nerlove (1971)
rem.walhus <- update(rem_gls_ind, random.method = "walhus")
rem.amemiya <- update(rem_gls_ind, random.method = "amemiya")
rem.nerlove <- update(rem_gls_ind, random.method = "nerlove")
rem.models_ind <- list(swar = rem_gls_ind, walhus = rem.walhus,amemiya = rem.amemiya, nerlove = rem.nerlove)
sapply(rem.models_ind, coef)## swar walhus amemiya nerlove
## (Intercept) -646.37648 -1326.135886 -378.251420 -154.581639
## IPM -13.69776 28.398240 -27.324754 -37.458430
## PDRB -3.07347 -4.254152 -2.713043 -2.456763
## Rata.rata_lama_sekolah 278.26799 10.464332 359.398637 416.923384One Way Time
Pada ilustrasi ini, pemodelan satu arah random effect model tetap dilakukan, meski dalam pengujian spesifikasi efek model, pengaruh waktu tidak signifikan (Terima \(H_{0}\))
library(plm)
model <- `Jml_Penduduk_Miskin` ~ `IPM` + `PDRB` + `Rata.rata_lama_sekolah`
rem_gls_time <- plm(model, data = datapanel,
index = c("Provinsi","Tahun"),
effect = "time",
model = "random", random.method = "walhus")
summary(rem_gls_time)## Oneway (time) effect Random Effect Model
## (Wallace-Hussain's transformation)
##
## Call:
## plm(formula = model, data = datapanel, effect = "time", model = "random",
## random.method = "walhus", index = c("Provinsi", "Tahun"))
##
## Balanced Panel: n = 34, T = 3, N = 102
##
## Effects:
## var std.dev share
## idiosyncratic 954417.2 976.9 1
## time 0.0 0.0 0
## theta: 0
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -1262.45 -587.34 -216.35 168.41 3112.41
##
## Coefficients:
## Estimate Std. Error z-value Pr(>|z|)
## (Intercept) -2896.9852 1966.8288 -1.4729 0.1407721
## IPM 150.9325 40.5246 3.7245 0.0001957 ***
## PDRB -2.2112 23.1730 -0.0954 0.9239816
## Rata.rata_lama_sekolah -817.8994 170.6326 -4.7933 1.64e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 117560000
## Residual Sum of Squares: 94936000
## R-Squared: 0.19241
## Adj. R-Squared: 0.16769
## Chisq: 23.3486 on 3 DF, p-value: 3.4159e-05kable(tidy(rem_gls_time), digits=3,
caption="GLS_Random Effect Model_Waktu")| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | -2896.985 | 1966.829 | -1.473 | 0.141 |
| IPM | 150.933 | 40.525 | 3.724 | 0.000 |
| PDRB | -2.211 | 23.173 | -0.095 | 0.924 |
| Rata.rata_lama_sekolah | -817.899 | 170.633 | -4.793 | 0.000 |
tidy_ranef_time <- tidy(ranef(rem_gls_time, effect="time"))## Warning: 'tidy.numeric' is deprecated.
## See help("Deprecated")kable(tidy_ranef_time, caption = "Pengaruh Acak Waktu", col.names = c("Tahun", "Pengaruh Acak"))| Tahun | Pengaruh Acak |
|---|---|
| 2019 | 0 |
| 2020 | 0 |
| 2021 | 0 |
Berdasarkan hasil tersebut menunjukkan pengaruh efek waktu bernilai 0 pada seluruh periode amatan. Hal ini berarti efek waktu bersifat invariant atau tidak ada pengaruh spesifik waktu.
MLE (One Way REM)
Pendugaan model random efect model dengan menggunakan MLE
dengan menggunakan package lme4 dan
lmerTest.
One Way Individual
library(lme4)## Loading required package: Matrix##
## Attaching package: 'Matrix'## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpacklibrary(lmerTest)##
## Attaching package: 'lmerTest'## The following object is masked from 'package:lme4':
##
## lmer## The following object is masked from 'package:stats':
##
## steprem_MLE_ind = lmer(Jml_Penduduk_Miskin ~ IPM+PDRB+Rata.rata_lama_sekolah + (1|Provinsi),
data=datapanel,
REML=FALSE)
summary(rem_MLE_ind)## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: Jml_Penduduk_Miskin ~ IPM + PDRB + Rata.rata_lama_sekolah + (1 |
## Provinsi)
## Data: datapanel
##
## AIC BIC logLik deviance df.resid
## 1412.1 1427.8 -700.0 1400.1 96
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.1813 -0.2042 -0.0708 0.3238 3.0098
##
## Random effects:
## Groups Name Variance Std.Dev.
## Provinsi (Intercept) 1386914 1177.67
## Residual 6069 77.91
## Number of obs: 102, groups: Provinsi, 34
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -460.798 2845.861 89.331 -0.162 0.8717
## IPM -23.307 51.361 94.963 -0.454 0.6510
## PDRB -2.818 2.129 67.607 -1.323 0.1902
## Rata.rata_lama_sekolah 335.885 144.528 95.535 2.324 0.0222 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) IPM PDRB
## IPM -0.954
## PDRB 0.165 -0.231
## Rt.rt_lm_sk 0.525 -0.751 0.296Kolom Groups menerangkan bahwa terdapat intersep unit individu yang bersifat random. Untuk dapat melihat nilai dari efek acak unit individu dengan menggunakan fungsi berikut :
coef(rem_MLE_ind)$Provinsi## (Intercept) IPM PDRB Rata.rata_lama_sekolah
## ACEH -609.354879 -23.30738 -2.817622 335.8854
## BALI -1068.140130 -23.30738 -2.817622 335.8854
## BANTEN -502.238590 -23.30738 -2.817622 335.8854
## BENGKULU -986.190625 -23.30738 -2.817622 335.8854
## DI YOGYAKARTA -848.467829 -23.30738 -2.817622 335.8854
## DKI JAKARTA -1391.058967 -23.30738 -2.817622 335.8854
## GORONTALO -825.147506 -23.30738 -2.817622 335.8854
## JAMBI -913.240444 -23.30738 -2.817622 335.8854
## JAWA BARAT 2672.801604 -23.30738 -2.817622 335.8854
## JAWA TENGAH 3032.944868 -23.30738 -2.817622 335.8854
## JAWA TIMUR 3404.969230 -23.30738 -2.817622 335.8854
## KALIMANTAN BARAT -523.255090 -23.30738 -2.817622 335.8854
## KALIMANTAN SELATAN -922.569076 -23.30738 -2.817622 335.8854
## KALIMANTAN TENGAH -1079.980748 -23.30738 -2.817622 335.8854
## KALIMANTAN TIMUR -1259.918197 -23.30738 -2.817622 335.8854
## KALIMANTAN UTARA -1312.837130 -23.30738 -2.817622 335.8854
## KEP. BANGKA BELITUNG -958.154941 -23.30738 -2.817622 335.8854
## KEP. RIAU -1487.628427 -23.30738 -2.817622 335.8854
## LAMPUNG -6.268028 -23.30738 -2.817622 335.8854
## MALUKU -1387.889546 -23.30738 -2.817622 335.8854
## MALUKU UTARA -1324.244343 -23.30738 -2.817622 335.8854
## NUSA TENGGARA BARAT -130.350396 -23.30738 -2.817622 335.8854
## NUSA TENGGARA TIMUR 118.619435 -23.30738 -2.817622 335.8854
## PAPUA 82.942037 -23.30738 -2.817622 335.8854
## PAPUA BARAT -813.826158 -23.30738 -2.817622 335.8854
## RIAU -868.103893 -23.30738 -2.817622 335.8854
## SULAWESI BARAT -937.893182 -23.30738 -2.817622 335.8854
## SULAWESI SELATAN -347.588352 -23.30738 -2.817622 335.8854
## SULAWESI TENGAH -916.190023 -23.30738 -2.817622 335.8854
## SULAWESI TENGGARA -1046.114760 -23.30738 -2.817622 335.8854
## SULAWESI UTARA -1290.843538 -23.30738 -2.817622 335.8854
## SUMATERA BARAT -972.721533 -23.30738 -2.817622 335.8854
## SUMATERA SELATAN -30.685413 -23.30738 -2.817622 335.8854
## SUMATERA UTARA -218.511509 -23.30738 -2.817622 335.8854ranef(rem_MLE_ind) #nilai yang disajikan merupakan perbedaan dengan base intersepnya## $Provinsi
## (Intercept)
## ACEH -148.55676
## BALI -607.34201
## BANTEN -41.44047
## BENGKULU -525.39251
## DI YOGYAKARTA -387.66971
## DKI JAKARTA -930.26085
## GORONTALO -364.34939
## JAMBI -452.44232
## JAWA BARAT 3133.59972
## JAWA TENGAH 3493.74299
## JAWA TIMUR 3865.76735
## KALIMANTAN BARAT -62.45697
## KALIMANTAN SELATAN -461.77096
## KALIMANTAN TENGAH -619.18263
## KALIMANTAN TIMUR -799.12008
## KALIMANTAN UTARA -852.03901
## KEP. BANGKA BELITUNG -497.35682
## KEP. RIAU -1026.83031
## LAMPUNG 454.53009
## MALUKU -927.09143
## MALUKU UTARA -863.44622
## NUSA TENGGARA BARAT 330.44772
## NUSA TENGGARA TIMUR 579.41756
## PAPUA 543.74016
## PAPUA BARAT -353.02804
## RIAU -407.30577
## SULAWESI BARAT -477.09506
## SULAWESI SELATAN 113.20977
## SULAWESI TENGAH -455.39190
## SULAWESI TENGGARA -585.31664
## SULAWESI UTARA -830.04542
## SUMATERA BARAT -511.92341
## SUMATERA SELATAN 430.11271
## SUMATERA UTARA 242.28661
##
## with conditional variances for "Provinsi"Perbandingan Model
# Dibagian Tambahanberdasarkan hasil pendugaan dengan metode GLS dan MLE memberikan nilai dugaan koefisien dan perolehan kebaikan model dengan kriteria AIC yang jauh berbeda.
One Way Time
library(lme4)
library(lmerTest)
re_MLE_time = lmer(Jml_Penduduk_Miskin ~ IPM+PDRB+Rata.rata_lama_sekolah + (1|Tahun),
data=datapanel,
REML=FALSE)## boundary (singular) fit: see help('isSingular')summary(re_MLE_time)## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: Jml_Penduduk_Miskin ~ IPM + PDRB + Rata.rata_lama_sekolah + (1 |
## Tahun)
## Data: datapanel
##
## AIC BIC logLik deviance df.resid
## 1703.3 1719.1 -845.7 1691.3 96
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.3086 -0.6088 -0.2243 0.1746 3.2261
##
## Random effects:
## Groups Name Variance Std.Dev.
## Tahun (Intercept) 0 0.0
## Residual 930749 964.8
## Number of obs: 102, groups: Tahun, 3
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -2896.985 1927.878 102.000 -1.503 0.136011
## IPM 150.933 39.722 102.000 3.800 0.000247 ***
## PDRB -2.211 22.714 102.000 -0.097 0.922642
## Rata.rata_lama_sekolah -817.899 167.253 102.000 -4.890 3.77e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) IPM PDRB
## IPM -0.883
## PDRB -0.100 0.102
## Rt.rt_lm_sk 0.401 -0.782 -0.107
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')Kolom Groups menerangkan bahwa terdapat intersep unit individu yang bersifat random. Untuk dapat melihat nilai dari efek acak unit individu dengan menggunakan fungsi berikut :
ranef(re_MLE_time)## $Tahun
## (Intercept)
## 2019 0
## 2020 0
## 2021 0
##
## with conditional variances for "Tahun"Tambahan
Selain dengan package lmer, pemodelan data
panel dengan pendugaan MLE dapat dilakukan dengan menggunakan
package pglm dan nlme.
Pemodelan MLE Lainnya
library(pglm)## Loading required package: maxLik## Loading required package: miscTools##
## Please cite the 'maxLik' package as:
## Henningsen, Arne and Toomet, Ott (2011). maxLik: A package for maximum likelihood estimation in R. Computational Statistics 26(3), 443-458. DOI 10.1007/s00180-010-0217-1.
##
## If you have questions, suggestions, or comments regarding the 'maxLik' package, please use a forum or 'tracker' at maxLik's R-Forge site:
## https://r-forge.r-project.org/projects/maxlik/rem.mle <- pglm(Jml_Penduduk_Miskin ~ IPM+PDRB+Rata.rata_lama_sekolah, data = datapanel, model="random",
effect="individual",
family = gaussian)
summary(rem.mle)## --------------------------------------------
## Maximum Likelihood estimation
## Newton-Raphson maximisation, 150 iterations
## Return code 4: Iteration limit exceeded (iterlim)
## Log-Likelihood: -702.2909
## 6 free parameters
## Estimates:
## Estimate Std. error t value Pr(> t)
## (Intercept) -140.679 3207.334 -0.044 0.96501
## IPM -38.008 57.299 -0.663 0.50712
## PDRB -2.444 2.124 -1.151 0.24983
## Rata.rata_lama_sekolah 419.843 152.543 2.752 0.00592 **
## sd.id 1580.732 5.438 290.680 < 2e-16 ***
## sd.idios 77.172 6.676 11.559 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## --------------------------------------------Dengan menggunakan library nlme
library(nlme)##
## Attaching package: 'nlme'## The following object is masked from 'package:lme4':
##
## lmList## The following object is masked from 'package:dplyr':
##
## collapsereML <- lme(Jml_Penduduk_Miskin ~ IPM+PDRB+Rata.rata_lama_sekolah, data = datapanel, random = ~1 | Provinsi)
summary(reML)## Linear mixed-effects model fit by REML
## Data: datapanel
## AIC BIC logLik
## 1375.597 1391.107 -681.7986
##
## Random effects:
## Formula: ~1 | Provinsi
## (Intercept) Residual
## StdDev: 1207.967 79.26773
##
## Fixed effects: Jml_Penduduk_Miskin ~ IPM + PDRB + Rata.rata_lama_sekolah
## Value Std.Error DF t-value p-value
## (Intercept) -451.4154 2907.1917 65 -0.1552754 0.8771
## IPM -23.7718 52.4486 65 -0.4532405 0.6519
## PDRB -2.8055 2.1669 65 -1.2946669 0.2000
## Rata.rata_lama_sekolah 338.6225 147.3143 65 2.2986397 0.0248
## Correlation:
## (Intr) IPM PDRB
## IPM -0.955
## PDRB 0.165 -0.231
## Rata.rata_lama_sekolah 0.526 -0.751 0.297
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -5.08978677 -0.20178682 -0.07010777 0.32160181 2.95502303
##
## Number of Observations: 102
## Number of Groups: 34#random komponen
ranef(reML)## (Intercept)
## ACEH -149.95173
## BALI -606.13430
## BANTEN -41.39125
## BENGKULU -525.77219
## DI YOGYAKARTA -385.98312
## DKI JAKARTA -932.56065
## GORONTALO -363.21654
## JAMBI -452.06345
## JAWA BARAT 3134.50087
## JAWA TENGAH 3496.86047
## JAWA TIMUR 3868.57202
## KALIMANTAN BARAT -60.60555
## KALIMANTAN SELATAN -460.87976
## KALIMANTAN TENGAH -619.08734
## KALIMANTAN TIMUR -799.71207
## KALIMANTAN UTARA -853.18651
## KEP. BANGKA BELITUNG -495.58127
## KEP. RIAU -1028.76476
## LAMPUNG 455.57098
## MALUKU -931.37752
## MALUKU UTARA -865.82597
## NUSA TENGGARA BARAT 332.75611
## NUSA TENGGARA TIMUR 579.45329
## PAPUA 544.18610
## PAPUA BARAT -352.96459
## RIAU -407.82800
## SULAWESI BARAT -477.33946
## SULAWESI SELATAN 114.29917
## SULAWESI TENGAH -456.70347
## SULAWESI TENGGARA -586.27878
## SULAWESI UTARA -831.57974
## SUMATERA BARAT -512.29868
## SUMATERA SELATAN 430.70296
## SUMATERA UTARA 240.18473Perbandingan Model
#AIC
AIC1 <- rbind(AIC(rem_gls_ind), AIC(reML))
#koefisien
koef1 <- rbind.data.frame(coef(rem_gls_ind), fixef(reML))
rownames(koef1) <- c("REM_GLS", "REM_MLE")
colnames(koef1) <- c("intersep", "IPM","PDRB","Rata-Rata Lama Sekolah")
koef1$AIC <- AIC1
koef1## intersep IPM PDRB Rata-Rata Lama Sekolah AIC
## REM_GLS -646.3765 -13.69776 -3.073470 278.2680 1199.255
## REM_MLE -451.4154 -23.77185 -2.805453 338.6225 1375.597Pemilihan FEM dan REM
Berdasarkan pengecekan model REM, model yang sesuai adalah REM dengan
efek individu, maka yang akan dilakukan pengecekan dengan uji Hausman
adalah model REM efek individu dengan model FEM individu. Uji Hausman
dengan menggunakan
phtest(fixed effect model, random effect model).
FEM dengan REM GLS
kable(tidy(phtest(fem, rem_gls_ind)), caption="Uji Hausman FEM_ind VS REM_ind GLS")| statistic | p.value | parameter | method | alternative |
|---|---|---|---|---|
| 46.62632 | 0 | 3 | Hausman Test | one model is inconsistent |
FEM dengan REM MLE package pglm
kable(tidy(phtest(fem, rem.mle)),caption="Uji Hausman FEM_ind VS REM_ind MLE")| statistic | p.value | parameter | method | alternative |
|---|---|---|---|---|
| 6.435471 | 0.0922424 | 3 | Hausman Test | one model is inconsistent |
Uji Hausmann Manual
# fem dengan rem_MLE_ind
# matrix covariance dari masing-masing model
V1 <- fem$vcov
V2 <- rem_MLE_ind@vcov_beta[2,2]
# Hitung perbedaan slope koefisien masing-masing model
bdiff <- coef(fem) - coef(rem_gls_ind)[2]
# Hitung perbedaan dari matrix covariance
Vdiff <- V1 - V2
# Hitung Hausman statistic
H <- t(bdiff) %*% solve(Vdiff) %*% bdiff
# Hitung nilai pvalue
pval <- 1 - pchisq(H, length(bdiff))
# Print hasil Hausman test
cat("Hausman test:\n")## Hausman test:cat("H = ", H, "\n")## H = 17.56217cat("p-value = ", pval, "\n")## p-value = 0.0005414392Karena nilai p-value yang diperoleh lebih dari alpha (5%), maka kita Gagal Tolak \(H_{0}\) yang menandakan bahwa model yang sesuai adalah model REM dengan pengaruh individu.
Apabila terindikasi terdapat gejala heteroskedastisitas, maka pengujian Hausman juga dapat dilakukan dengan menggunakan uji Hausman robust dengan menyertakan matriks robust covariance untu mendapatkan hasil pengujian yang konsisten.
kable(tidy(phtest(fem, rem_gls_ind, method = "aux", vcov = vcovHC)),caption="Uji Hausman Robust FEM_ind VS REM_ind GLS")| statistic | p.value | parameter | method | alternative |
|---|---|---|---|---|
| 46.62632 | 0 | 3 | Hausman Test | one model is inconsistent |
Berdasarkan pengujian robust Hausman, keputusan yang diambil tetap gagal Tolak \(H_{0}\) yang berimplikasi model yang lebih baik adalah model pengaruh acak efek individu.
Kemudian, tahapan selanjutnya adalah perlu memeriksa antara model pengaruh acak dengan common effect model dengan uji pengganda lagrang.
cem <- plm(model, datapanel, model="pooling")
c.tw <- plmtest(cem, effect = "twoways", type = c("bp"))
c.ind <- plmtest(cem, effect = "individual", type = c("bp"))
c.time <-plmtest(cem, effect = "time", type = c("bp"))
tidy_tests <- bind_rows(tidy(c.tw), tidy(c.ind), tidy(c.time))
colnames(tidy_tests) <- c("Chisq", "p-value", "df","Jenis-Uji", "H1")
kable(tidy_tests, digits=2,caption = "Hasil Uji LM Breusch-Pagan")| Chisq | p-value | df | Jenis-Uji | H1 |
|---|---|---|---|---|
| 99.86 | 0.0 | 2 | Lagrange Multiplier Test - two-ways effects (Breusch-Pagan) | significant effects |
| 98.78 | 0.0 | 1 | Lagrange Multiplier Test - (Breusch-Pagan) | significant effects |
| 1.09 | 0.3 | 1 | Lagrange Multiplier Test - time effects (Breusch-Pagan) | significant effects |
Berdasarkan pengecekan dengan uji pengganda lagrange, dapat ditunjukkan bahwa terdapat efek panel (Efek two ways, dan efek individu signifikan pada taraf 5%) sehingga pemodelan dengan pooling saja tidak lah cukup.
Dengan demikian model yang terpilih adalah model pengaruh acak dengan efek satu arah.
Final Model
final_model <- rem_gls_ind
kable(tidy(final_model), digits=3,caption="Model Pengaruh Acak Efek Individu")| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | -646.376 | 2828.078 | -0.229 | 0.819 |
| IPM | -13.698 | 51.427 | -0.266 | 0.790 |
| PDRB | -3.073 | 2.287 | -1.344 | 0.179 |
| Rata.rata_lama_sekolah | 278.268 | 150.221 | 1.852 | 0.064 |
kable(tidy_ranef_ind,digits=3, caption = "Pengaruh Acak Individu", col.names = c("Provinsi", "Pengaruh Acak Individu"))| Provinsi | Pengaruh Acak Individu |
|---|---|
| ACEH | -119.059 |
| BALI | -632.002 |
| BANTEN | -42.244 |
| BENGKULU | -517.346 |
| DI YOGYAKARTA | -421.653 |
| DKI JAKARTA | -880.222 |
| GORONTALO | -388.581 |
| JAMBI | -460.361 |
| JAWA BARAT | 3114.722 |
| JAWA TENGAH | 3428.180 |
| JAWA TIMUR | 3806.751 |
| KALIMANTAN BARAT | -101.980 |
| KALIMANTAN SELATAN | -480.539 |
| KALIMANTAN TENGAH | -621.185 |
| KALIMANTAN TIMUR | -785.735 |
| KALIMANTAN UTARA | -827.901 |
| KEP. BANGKA BELITUNG | -534.644 |
| KEP. RIAU | -985.355 |
| LAMPUNG | 432.377 |
| MALUKU | -837.159 |
| MALUKU UTARA | -813.776 |
| NUSA TENGGARA BARAT | 281.401 |
| NUSA TENGGARA TIMUR | 577.657 |
| PAPUA | 532.581 |
| PAPUA BARAT | -355.380 |
| RIAU | -396.018 |
| SULAWESI BARAT | -472.792 |
| SULAWESI SELATAN | 90.417 |
| SULAWESI TENGAH | -428.045 |
| SULAWESI TENGGARA | -565.011 |
| SULAWESI UTARA | -797.421 |
| SUMATERA BARAT | -503.794 |
| SUMATERA SELATAN | 417.502 |
| SUMATERA UTARA | 286.614 |
Model Akhir
\[ \begin{aligned} I\hat nv_{it} &= -646.376-13.698IPM_{it}+f_{i}\\ \end{aligned} \]
Setiap penurunan satu satuan IPM maka akan menyebabkan Jumlah Penduduk Miskin turun sebesar -13.698 dengan mengganggap peubah lain konstan (apabila memasukkan peubah lain).
Namun, karena nilai intersep pada model pengaruh acak berbeda-beda pada setiap unit individu, maka pada setiap unit individu akan memiliki persamaan masing-masing sesuai dengan pengaruh acak yang dimilikinya, misalnya pada Provinsi 1 maka :
\[ \begin{aligned} I\hat nv_{it} &= -646.376-13.698IPM_{it}-119.059\\ \end{aligned} \]
Demikian pula untuk persamaan model pengaruh acak individu pada Provinsi lainnya ( 2,3,…10)