=====================================================================
ANALISIS DATA PANEL MODEL LINIER - DATASET GRUNFELD (paket plm)
=====================================================================
0. PAKET
paket <- c("plm", "lmtest", "sandwich", "car", "stargazer")
for (p in paket) {
if (!requireNamespace(p, quietly = TRUE)) install.packages(p)
library(p, character.only = TRUE)
}
## Warning: package 'plm' was built under R version 4.5.3
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
## Warning: package 'sandwich' was built under R version 4.5.2
## Warning: package 'car' was built under R version 4.5.2
## Loading required package: carData
## Warning: package 'carData' was built under R version 4.5.2
## Warning: package 'stargazer' was built under R version 4.5.2
##
## 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=stargazer
———————————————————————
1. DATA DAN STRUKTUR PANEL
———————————————————————
data("Grunfeld", package = "plm")
# inv = investasi, value = nilai pasar perusahaan, capital = stok modal
head(Grunfeld)
## firm year inv value capital
## 1 1 1935 317.6 3078.5 2.8
## 2 1 1936 391.8 4661.7 52.6
## 3 1 1937 410.6 5387.1 156.9
## 4 1 1938 257.7 2792.2 209.2
## 5 1 1939 330.8 4313.2 203.4
## 6 1 1940 461.2 4643.9 207.2
str(Grunfeld)
## 'data.frame': 200 obs. of 5 variables:
## $ firm : int 1 1 1 1 1 1 1 1 1 1 ...
## $ year : int 1935 1936 1937 1938 1939 1940 1941 1942 1943 1944 ...
## $ inv : num 318 392 411 258 331 ...
## $ value : num 3078 4662 5387 2792 4313 ...
## $ capital: num 2.8 52.6 156.9 209.2 203.4 ...
pdata <- pdata.frame(Grunfeld, index = c("firm", "year"))
pdim(pdata) # cek balanced / unbalanced, N, T, n
## Balanced Panel: n = 10, T = 20, N = 200
is.pbalanced(pdata)
## [1] TRUE
summary(pdata)
## firm year inv value
## 1 :20 1935 : 10 Min. : 0.93 Min. : 58.12
## 2 :20 1936 : 10 1st Qu.: 33.56 1st Qu.: 199.97
## 3 :20 1937 : 10 Median : 57.48 Median : 517.95
## 4 :20 1938 : 10 Mean : 145.96 Mean :1081.68
## 5 :20 1939 : 10 3rd Qu.: 138.04 3rd Qu.:1679.85
## 6 :20 1940 : 10 Max. :1486.70 Max. :6241.70
## (Other):80 (Other):140
## capital
## Min. : 0.80
## 1st Qu.: 79.17
## Median : 205.60
## Mean : 276.02
## 3rd Qu.: 358.10
## Max. :2226.30
##
———————————————————————
2. EKSPLORASI DATA
———————————————————————
# Statistik deskriptif
stargazer(Grunfeld[, c("inv", "value", "capital")], type = "text",
title = "Statistik Deskriptif", digits = 2)
##
## Statistik Deskriptif
## ==============================================
## Statistic N Mean St. Dev. Min Max
## ----------------------------------------------
## inv 200 145.96 216.88 0.93 1,486.70
## value 200 1,081.68 1,314.47 58.12 6,241.70
## capital 200 276.02 301.10 0.80 2,226.30
## ----------------------------------------------
# Korelasi antar variabel
cor(Grunfeld[, c("inv", "value", "capital")])
## inv value capital
## inv 1.0000000 0.8569329 0.6625653
## value 0.8569329 1.0000000 0.4887058
## capital 0.6625653 0.4887058 1.0000000
# Rata-rata per perusahaan dan per tahun
aggregate(inv ~ firm, data = Grunfeld, FUN = mean)
## firm inv
## 1 1 608.0200
## 2 2 410.4750
## 3 3 102.2900
## 4 4 86.1235
## 5 5 61.8025
## 6 6 55.4110
## 7 7 47.5955
## 8 8 42.8915
## 9 9 41.8890
## 10 10 3.0845
aggregate(inv ~ year, data = Grunfeld, FUN = mean)
## year inv
## 1 1935 72.746
## 2 1936 101.607
## 3 1937 122.481
## 4 1938 77.555
## 5 1939 80.526
## 6 1940 113.265
## 7 1941 139.719
## 8 1942 122.665
## 9 1943 117.790
## 10 1944 120.925
## 11 1945 124.159
## 12 1946 161.359
## 13 1947 147.135
## 14 1948 153.948
## 15 1949 139.244
## 16 1950 151.061
## 17 1951 199.583
## 18 1952 224.033
## 19 1953 275.583
## 20 1954 273.781
# Visualisasi
par(mfrow = c(1, 2))
boxplot(inv ~ firm, data = Grunfeld, main = "Investasi per Perusahaan",
xlab = "Firm", ylab = "inv")
plot(tapply(Grunfeld$inv, Grunfeld$year, mean), type = "b",
xaxt = "n", main = "Rata-rata Investasi per Tahun",
xlab = "Tahun", ylab = "inv")
axis(1, at = 1:20, labels = sort(unique(Grunfeld$year)))
par(mfrow = c(1, 1))
———————————————————————
3. ESTIMASI TIGA MODEL
———————————————————————
form <- inv ~ value + capital
cem <- plm(form, data = pdata, model = "pooling") # Common Effect
fem <- plm(form, data = pdata, model = "within") # Fixed Effect
rem <- plm(form, data = pdata, model = "random") # Random Effect
summary(cem)
## Pooling Model
##
## Call:
## plm(formula = form, data = pdata, model = "pooling")
##
## Balanced Panel: n = 10, T = 20, N = 200
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -291.6757 -30.0137 5.3033 34.8293 369.4464
##
## Coefficients:
## Estimate Std. Error t-value Pr(>|t|)
## (Intercept) -42.7143694 9.5116760 -4.4907 1.207e-05 ***
## value 0.1155622 0.0058357 19.8026 < 2.2e-16 ***
## capital 0.2306785 0.0254758 9.0548 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 9359900
## Residual Sum of Squares: 1755900
## R-Squared: 0.81241
## Adj. R-Squared: 0.8105
## F-statistic: 426.576 on 2 and 197 DF, p-value: < 2.22e-16
summary(fem)
## Oneway (individual) effect Within Model
##
## Call:
## plm(formula = form, data = pdata, model = "within")
##
## Balanced Panel: n = 10, T = 20, N = 200
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -184.00857 -17.64316 0.56337 19.19222 250.70974
##
## Coefficients:
## Estimate Std. Error t-value Pr(>|t|)
## value 0.110124 0.011857 9.2879 < 2.2e-16 ***
## capital 0.310065 0.017355 17.8666 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 2244400
## Residual Sum of Squares: 523480
## R-Squared: 0.76676
## Adj. R-Squared: 0.75311
## F-statistic: 309.014 on 2 and 188 DF, p-value: < 2.22e-16
summary(rem)
## Oneway (individual) effect Random Effect Model
## (Swamy-Arora's transformation)
##
## Call:
## plm(formula = form, data = pdata, model = "random")
##
## Balanced Panel: n = 10, T = 20, N = 200
##
## Effects:
## var std.dev share
## idiosyncratic 2784.46 52.77 0.282
## individual 7089.80 84.20 0.718
## theta: 0.8612
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -177.6063 -19.7350 4.6851 19.5105 252.8743
##
## Coefficients:
## Estimate Std. Error z-value Pr(>|z|)
## (Intercept) -57.834415 28.898935 -2.0013 0.04536 *
## value 0.109781 0.010493 10.4627 < 2e-16 ***
## capital 0.308113 0.017180 17.9339 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 2381400
## Residual Sum of Squares: 548900
## R-Squared: 0.7695
## Adj. R-Squared: 0.76716
## Chisq: 657.674 on 2 DF, p-value: < 2.22e-16
# Efek individu (intersep tiap perusahaan) pada FEM
fixef(fem)
## 1 2 3 4 5 6 7 8
## -70.2967 101.9058 -235.5718 -27.8093 -114.6168 -23.1613 -66.5535 -57.5457
## 9 10
## -87.2223 -6.5678
summary(fixef(fem))
## Estimate Std. Error t-value Pr(>|t|)
## 1 -70.2967 49.7080 -1.4142 0.15896
## 2 101.9058 24.9383 4.0863 6.485e-05 ***
## 3 -235.5718 24.4316 -9.6421 < 2.2e-16 ***
## 4 -27.8093 14.0778 -1.9754 0.04969 *
## 5 -114.6168 14.1654 -8.0913 7.141e-14 ***
## 6 -23.1613 12.6687 -1.8282 0.06910 .
## 7 -66.5535 12.8430 -5.1821 5.629e-07 ***
## 8 -57.5457 13.9931 -4.1124 5.848e-05 ***
## 9 -87.2223 12.8919 -6.7657 1.635e-10 ***
## 10 -6.5678 11.8269 -0.5553 0.57933
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Komponen varians pada REM
ercomp(rem)
## var std.dev share
## idiosyncratic 2784.46 52.77 0.282
## individual 7089.80 84.20 0.718
## theta: 0.8612
# Perbandingan tabel
stargazer(cem, fem, rem, type = "text",
column.labels = c("CEM", "FEM", "REM"),
model.names = FALSE)
##
## =========================================================================
## Dependent variable:
## ------------------------------------------------------------
## inv
## CEM FEM REM
## (1) (2) (3)
## -------------------------------------------------------------------------
## value 0.116*** 0.110*** 0.110***
## (0.006) (0.012) (0.010)
##
## capital 0.231*** 0.310*** 0.308***
## (0.025) (0.017) (0.017)
##
## Constant -42.714*** -57.834**
## (9.512) (28.899)
##
## -------------------------------------------------------------------------
## Observations 200 200 200
## R2 0.812 0.767 0.770
## Adjusted R2 0.811 0.753 0.767
## F Statistic 426.576*** (df = 2; 197) 309.014*** (df = 2; 188) 657.674***
## =========================================================================
## Note: *p<0.1; **p<0.05; ***p<0.01
———————————————————————
4. UJI PEMILIHAN MODEL
———————————————————————
chow <- pFtest(fem, cem) # CEM vs FEM
bp_lm <- plmtest(cem, type = "bp") # CEM vs REM
hausman <- phtest(fem, rem) # FEM vs REM
print(chow)
##
## F test for individual effects
##
## data: form
## F = 49.177, df1 = 9, df2 = 188, p-value < 2.2e-16
## alternative hypothesis: significant effects
print(bp_lm)
##
## Lagrange Multiplier Test - (Breusch-Pagan)
##
## data: form
## chisq = 798.16, df = 1, p-value < 2.2e-16
## alternative hypothesis: significant effects
print(hausman)
##
## Hausman Test
##
## data: form
## chisq = 2.3304, df = 2, p-value = 0.3119
## alternative hypothesis: one model is inconsistent
# Kesimpulan otomatis (alpha = 5%)
alpha <- 0.05
cat("\n===== KESIMPULAN UJI PEMILIHAN MODEL =====\n")
##
## ===== KESIMPULAN UJI PEMILIHAN MODEL =====
cat("Uji Chow : p =", signif(chow$p.value, 4), "->",
ifelse(chow$p.value < alpha, "FEM lebih baik dari CEM", "CEM cukup"), "\n")
## Uji Chow : p = 8.7e-45 -> FEM lebih baik dari CEM
cat("Uji BP-LM : p =", signif(bp_lm$p.value, 4), "->",
ifelse(bp_lm$p.value < alpha, "REM lebih baik dari CEM", "CEM cukup"), "\n")
## Uji BP-LM : p = 1.354e-175 -> REM lebih baik dari CEM
cat("Uji Hausman : p =", signif(hausman$p.value, 4), "->",
ifelse(hausman$p.value < alpha, "FEM lebih baik dari REM",
"REM lebih baik dari FEM"), "\n")
## Uji Hausman : p = 0.3119 -> REM lebih baik dari FEM
model_terpilih <- if (chow$p.value >= alpha && bp_lm$p.value >= alpha) {
"CEM"
} else if (chow$p.value < alpha && bp_lm$p.value >= alpha) {
"FEM"
} else if (chow$p.value >= alpha && bp_lm$p.value < alpha) {
"REM"
} else {
ifelse(hausman$p.value < alpha, "FEM", "REM")
}
cat("Model terpilih:", model_terpilih, "\n\n")
## Model terpilih: REM
———————————————————————
5. TWO-WAY EFFECTS (INDIVIDU + WAKTU)
———————————————————————
fem_tw <- plm(form, data = pdata, model = "within", effect = "twoways")
summary(fem_tw)
## Twoways effects Within Model
##
## Call:
## plm(formula = form, data = pdata, effect = "twoways", model = "within")
##
## Balanced Panel: n = 10, T = 20, N = 200
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -162.6094 -19.4710 -1.2669 19.1277 211.8420
##
## Coefficients:
## Estimate Std. Error t-value Pr(>|t|)
## value 0.117716 0.013751 8.5604 6.653e-15 ***
## capital 0.357916 0.022719 15.7540 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 1615600
## Residual Sum of Squares: 452150
## R-Squared: 0.72015
## Adj. R-Squared: 0.67047
## F-statistic: 217.442 on 2 and 169 DF, p-value: < 2.22e-16
pFtest(fem_tw, fem) # perlu efek waktu? p < 0.05 -> ya
##
## F test for twoways effects
##
## data: form
## F = 1.4032, df1 = 19, df2 = 169, p-value = 0.1309
## alternative hypothesis: significant effects
# Uji LM untuk efek individu, waktu, dan keduanya
plmtest(cem, effect = "individual", type = "bp")
##
## Lagrange Multiplier Test - (Breusch-Pagan)
##
## data: form
## chisq = 798.16, df = 1, p-value < 2.2e-16
## alternative hypothesis: significant effects
plmtest(cem, effect = "time", type = "bp")
##
## Lagrange Multiplier Test - time effects (Breusch-Pagan)
##
## data: form
## chisq = 6.4539, df = 1, p-value = 0.01107
## alternative hypothesis: significant effects
plmtest(cem, effect = "twoways", type = "bp")
##
## Lagrange Multiplier Test - two-ways effects (Breusch-Pagan)
##
## data: form
## chisq = 804.62, df = 2, p-value < 2.2e-16
## alternative hypothesis: significant effects
# REM two-way (opsional)
rem_tw <- plm(form, data = pdata, model = "random", effect = "twoways")
phtest(fem_tw, rem_tw)
##
## Hausman Test
##
## data: form
## chisq = 13.46, df = 2, p-value = 0.001194
## alternative hypothesis: one model is inconsistent
———————————————————————
6. UJI ASUMSI
———————————————————————
# 6a. Heteroskedastisitas (Breusch-Pagan, dengan dummy unit)
bptest(inv ~ value + capital + factor(firm), data = Grunfeld)
##
## studentized Breusch-Pagan test
##
## data: inv ~ value + capital + factor(firm)
## BP = 85.836, df = 11, p-value = 1.086e-13
# 6b. Autokorelasi serial
pbgtest(fem) # Breusch-Godfrey panel
##
## Breusch-Godfrey/Wooldridge test for serial correlation in panel models
##
## data: form
## chisq = 65.063, df = 20, p-value = 1.14e-06
## alternative hypothesis: serial correlation in idiosyncratic errors
pwartest(fem) # Wooldridge
##
## Wooldridge's test for serial correlation in FE panels
##
## data: fem
## F = 76.929, df1 = 1, df2 = 188, p-value = 1.057e-15
## alternative hypothesis: serial correlation
pdwtest(fem) # Durbin-Watson panel
##
## Durbin-Watson test for serial correlation in panel models
##
## data: form
## DW = 1.0789, p-value = 3.184e-11
## alternative hypothesis: serial correlation in idiosyncratic errors
# 6c. Dependensi antarunit (cross-sectional dependence)
pcdtest(fem, test = "cd") # Pesaran CD
##
## Pesaran CD test for cross-sectional dependence in panels
##
## data: inv ~ value + capital
## z = 4.6612, p-value = 3.144e-06
## alternative hypothesis: cross-sectional dependence
pcdtest(fem, test = "lm") # Breusch-Pagan LM
##
## Breusch-Pagan LM test for cross-sectional dependence in panels
##
## data: inv ~ value + capital
## chisq = 246.33, df = 45, p-value < 2.2e-16
## alternative hypothesis: cross-sectional dependence
# 6d. Multikolinearitas
vif(lm(inv ~ value + capital, data = Grunfeld))
## value capital
## 1.313773 1.313773
cor(Grunfeld[, c("value", "capital")])
## value capital
## value 1.0000000 0.4887058
## capital 0.4887058 1.0000000
# 6e. Normalitas residual (informatif)
res_fem <- residuals(fem)
shapiro.test(as.numeric(res_fem))
##
## Shapiro-Wilk normality test
##
## data: as.numeric(res_fem)
## W = 0.90975, p-value = 1.1e-09
qqnorm(as.numeric(res_fem)); qqline(as.numeric(res_fem))
# 6f. Akar unit / stasioneritas
purtest(pdata$inv, test = "levinlin", exo = "intercept",
lags = "AIC", pmax = 4)
## Warning in selectT(l, theTs): the time series is short
##
## Levin-Lin-Chu Unit-Root Test (ex. var.: Individual Intercepts)
##
## data: pdata$inv
## z = 3.1416, p-value = 0.9992
## alternative hypothesis: stationarity
purtest(pdata$value, test = "levinlin", exo = "intercept",
lags = "AIC", pmax = 4)
## Warning in selectT(l, theTs): the time series is short
##
## Levin-Lin-Chu Unit-Root Test (ex. var.: Individual Intercepts)
##
## data: pdata$value
## z = 1.7905, p-value = 0.9633
## alternative hypothesis: stationarity
purtest(pdata$capital, test = "levinlin", exo = "intercept",
lags = "AIC", pmax = 4)
## Warning in selectT(l, theTs): the time series is short
##
## Levin-Lin-Chu Unit-Root Test (ex. var.: Individual Intercepts)
##
## data: pdata$capital
## z = 3.5537, p-value = 0.9998
## alternative hypothesis: stationarity
# Alternatif (lebih toleran terhadap heterogenitas antar unit)
purtest(pdata$inv, test = "ips", exo = "intercept", lags = "AIC", pmax = 4)
##
## Im-Pesaran-Shin Unit-Root Test (ex. var.: Individual Intercepts)
##
## data: pdata$inv
## Wtbar = 3.0768, p-value = 0.999
## alternative hypothesis: stationarity
———————————————————————
7. PENANGANAN PELANGGARAN: STANDARD ERROR ROBUST
———————————————————————
# Cluster per unit (robust terhadap heteroskedastisitas & autokorelasi)
coeftest(fem, vcov = vcovHC(fem, type = "HC1", cluster = "group"))
##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## value 0.110124 0.014414 7.6398 1.076e-12 ***
## capital 0.310065 0.050043 6.1959 3.576e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Driscoll-Kraay (jika ada cross-sectional dependence)
coeftest(fem, vcov = vcovSCC(fem, type = "HC1"))
##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## value 0.110124 0.017775 6.1954 3.586e-09 ***
## capital 0.310065 0.034996 8.8601 6.017e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Perbandingan SE biasa vs robust
se_biasa <- sqrt(diag(vcov(fem)))
se_cluster <- sqrt(diag(vcovHC(fem, type = "HC1", cluster = "group")))
se_dk <- sqrt(diag(vcovSCC(fem, type = "HC1")))
print(round(cbind(SE_biasa = se_biasa, SE_cluster = se_cluster,
SE_DriscollKraay = se_dk), 4))
## SE_biasa SE_cluster SE_DriscollKraay
## value 0.0119 0.0144 0.0178
## capital 0.0174 0.0500 0.0350
# Tabel akhir dengan SE robust
stargazer(fem, fem, type = "text",
se = list(se_biasa, se_cluster),
column.labels = c("FEM (SE biasa)", "FEM (SE cluster)"),
model.names = FALSE)
##
## =========================================================
## Dependent variable:
## -------------------------------
## inv
## FEM (SE biasa) FEM (SE cluster)
## (1) (2)
## ---------------------------------------------------------
## value 0.110*** 0.110***
## (0.012) (0.014)
##
## capital 0.310*** 0.310***
## (0.017) (0.050)
##
## ---------------------------------------------------------
## Observations 200 200
## R2 0.767 0.767
## Adjusted R2 0.753 0.753
## F Statistic (df = 2; 188) 309.014*** 309.014***
## =========================================================
## Note: *p<0.1; **p<0.05; ***p<0.01
———————————————————————
8. MODEL AKHIR DAN INTERPRETASI
———————————————————————
Ganti ‘fem’ dengan model_terpilih bila hasilnya REM / CEM
model_akhir <- switch(model_terpilih, CEM = cem, FEM = fem, REM = rem)
cat("\n===== MODEL AKHIR:", model_terpilih, "(SE cluster per unit) =====\n")
##
## ===== MODEL AKHIR: REM (SE cluster per unit) =====
print(coeftest(model_akhir,
vcov = vcovHC(model_akhir, type = "HC1", cluster = "group")))
##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -57.834415 23.627502 -2.4478 0.01525 *
## value 0.109781 0.013083 8.3914 9.224e-15 ***
## capital 0.308113 0.052283 5.8932 1.618e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# R-squared model akhir
r.squared(model_akhir)
## [1] 0.7695027
# Nilai prediksi vs aktual (FEM)
aktual <- as.numeric(pdata$inv)
fitted_fem <- aktual - as.numeric(residuals(fem))
plot(fitted_fem, aktual, pch = 16, xlab = "Fitted", ylab = "Aktual",
main = "Aktual vs Fitted (FEM)")
abline(0, 1, col = "red")
