library(plm)
## Warning: package 'plm' was built under R version 4.3.3
library(lmtest)
## Warning: package 'lmtest' was built under R version 4.3.2
## Loading required package: zoo
## Warning: package 'zoo' was built under R version 4.3.2
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
library(car)
## Warning: package 'car' was built under R version 4.3.3
## Loading required package: carData
library(readxl)
## Warning: package 'readxl' was built under R version 4.3.3
data <- readxl::read_excel("D:\\NSC\\2023 Data Panel NTT Ketahanan Pangan.xlsx")
head(data)
## # A tibble: 6 × 9
## No Kab Tahun AHH IKP IPM PDRB Prod.Pangan IKP_cat
## <dbl> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 293 Sumba Barat 2020 67.1 63.8 63.5 2351536. 5242 4
## 2 293 Sumba Barat 2021 67.1 64.6 63.8 2416824. 21216 4
## 3 293 Sumba Barat 2022 67.4 68.1 64.4 2583328. 18571 5
## 4 294 Sumba Timur 2020 65.1 66.7 65.5 6254369. 31354 5
## 5 294 Sumba Timur 2021 65.2 66.2 65.7 6403410. 33357 5
## 6 294 Sumba Timur 2022 65.4 69.3 66.2 6825709. 31040 5
pdata <- pdata.frame(data, index = c("Kab","Tahun"))
head(pdata)
## No Kab Tahun AHH IKP IPM PDRB Prod.Pangan IKP_cat
## Alor-2020 299 Alor 2020 61.48 63.30 61.33 3052785 55477 4
## Alor-2021 299 Alor 2021 61.64 60.82 61.37 3166850 45011 4
## Alor-2022 299 Alor 2022 61.99 58.97 62.26 3362313 56341 4
## Belu-2020 298 Belu 2020 64.61 71.06 62.68 4549707 59687 5
## Belu-2021 298 Belu 2021 64.89 71.89 62.77 4699017 49834 5
## Belu-2022 298 Belu 2022 65.28 70.48 63.22 5037063 65998 5
# Model FEM (Fixed Effect)
fem <- plm(IKP ~ AHH + PDRB + Prod.Pangan, data=pdata, model="within")
summary(fem)
## Oneway (individual) effect Within Model
##
## Call:
## plm(formula = IKP ~ AHH + PDRB + Prod.Pangan, data = pdata, model = "within")
##
## Balanced Panel: n = 22, T = 3, N = 66
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -6.700012 -0.811781 0.016126 0.864163 5.164580
##
## Coefficients:
## Estimate Std. Error t-value Pr(>|t|)
## AHH 2.2827e+00 2.1199e+00 1.0768 0.2879
## PDRB -3.1641e-07 1.3913e-06 -0.2274 0.8212
## Prod.Pangan 1.1003e-05 8.5189e-06 1.2916 0.2037
##
## Total Sum of Squares: 183.91
## Residual Sum of Squares: 168.75
## R-Squared: 0.082422
## Adj. R-Squared: -0.4547
## F-statistic: 1.22762 on 3 and 41 DF, p-value: 0.31192
# Model REM (Random Effect)
rem <- plm(IKP ~ AHH + PDRB + Prod.Pangan, data=pdata, model="random")
summary(rem)
## Oneway (individual) effect Random Effect Model
## (Swamy-Arora's transformation)
##
## Call:
## plm(formula = IKP ~ AHH + PDRB + Prod.Pangan, data = pdata, model = "random")
##
## Balanced Panel: n = 22, T = 3, N = 66
##
## Effects:
## var std.dev share
## idiosyncratic 4.116 2.029 0.103
## individual 35.722 5.977 0.897
## theta: 0.8077
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -8.387162 -0.949613 0.058879 1.060229 3.591477
##
## Coefficients:
## Estimate Std. Error z-value Pr(>|z|)
## (Intercept) -3.2179e+01 3.8659e+01 -0.8324 0.40519
## AHH 1.4822e+00 5.8993e-01 2.5125 0.01199 *
## PDRB 1.7881e-07 2.7829e-07 0.6425 0.52053
## Prod.Pangan 9.6076e-06 7.6536e-06 1.2553 0.20937
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 286.19
## Residual Sum of Squares: 244.36
## R-Squared: 0.14614
## Adj. R-Squared: 0.10483
## Chisq: 10.6118 on 3 DF, p-value: 0.014022
# Uji Chow (Pooled vs FEM)
pooling <- plm(IKP ~ AHH + PDRB + Prod.Pangan, data=pdata, model="pooling")
pFtest(fem, pooling)
##
## F test for individual effects
##
## data: IKP ~ AHH + PDRB + Prod.Pangan
## F = 23.237, df1 = 21, df2 = 41, p-value < 2.2e-16
## alternative hypothesis: significant effects
# Uji Hausman (FEM vs REM)
phtest(fem, rem)
##
## Hausman Test
##
## data: IKP ~ AHH + PDRB + Prod.Pangan
## chisq = 0.33242, df = 3, p-value = 0.9538
## alternative hypothesis: one model is inconsistent
# Lagrange Multiplier Test (LM Test)
plmtest(pooling, type = "bp")
##
## Lagrange Multiplier Test - (Breusch-Pagan)
##
## data: IKP ~ AHH + PDRB + Prod.Pangan
## chisq = 51.178, df = 1, p-value = 8.436e-13
## alternative hypothesis: significant effects
# --- Uji Asumsi Klasik ---
# Multikolinearitas (VIF)
vif(lm(IKP ~ AHH + PDRB + Prod.Pangan, data=data))
## AHH PDRB Prod.Pangan
## 1.119175 1.116783 1.003461
# Normalitas residual
resid_rem <- residuals(rem)
shapiro.test(resid_rem)
##
## Shapiro-Wilk normality test
##
## data: resid_rem
## W = 0.93103, p-value = 0.001222
tseries::jarque.bera.test(resid_rem)
## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo
##
## Jarque Bera Test
##
## data: resid_rem
## X-squared = 53.212, df = 2, p-value = 2.787e-12
# Heteroskedastisitas (Breusch-Pagan)
bptest(rem)
##
## studentized Breusch-Pagan test
##
## data: rem
## BP = 4.3696, df = 3, p-value = 0.2242
# Standardisasi residu ke Z-score
resid_std <- (resid_rem - mean(resid_rem)) / sd(resid_rem)
# Kolmogorov–Smirnov Test terhadap Normal(0,1)
ks.test(resid_std, "pnorm", mean=0, sd=1)
##
## Exact one-sample Kolmogorov-Smirnov test
##
## data: resid_std
## D = 0.088546, p-value = 0.6461
## alternative hypothesis: two-sided