Econometric Models in Practice
OLS, IV, Panel Data Model (Fixed Effects), DiD, Synthetic Control, and VAR
I Gede Nyoman Mindra Jaya, Ph.D
03 November 2025
Kata Pengantar
Puji syukur saya panjatkan ke hadirat Tuhan Yang Maha Esa, karena atas rahmat dan karunia-Nya buku Econometric Models in Practice ini dapat tersusun. Buku ini disusun sebagai panduan praktis bagi mahasiswa, peneliti, dan praktisi yang ingin memahami serta menerapkan model ekonometrika modern dalam analisis data nyata. Pembahasan dimulai dari model dasar seperti Ordinary Least Squares (OLS), kemudian berlanjut ke model yang lebih kompleks seperti Instrumental Variables (IV), Fixed Effects Panel Model, Difference-in-Differences (DiD), Synthetic Control, hingga Vector Autoregression (VAR).
Pendekatan yang digunakan dalam buku ini menekankan pada keseimbangan antara konsep teoritis dan penerapan empiris. Setiap model tidak hanya dijelaskan secara matematis, tetapi juga disertai contoh implementasi menggunakan perangkat lunak statistik yang umum digunakan. Harapannya, pembaca tidak hanya memahami rumus, tetapi juga mampu menghubungkan model dengan konteks ekonomi dan kebijakan publik yang relevan.
Saya menyadari bahwa kesempurnaan adalah proses yang berkelanjutan. Oleh karena itu, kritik dan saran yang membangun sangat saya harapkan untuk penyempurnaan edisi berikutnya. Semoga buku ini dapat menjadi referensi yang bermanfaat dalam memperkuat literasi ekonometrika di kalangan akademisi dan praktisi di Indonesia.
Bandung, I Gede Nyoman Mindra Jaya, Ph.D
1 Ordinary Least Squares (OLS)
Model dasar untuk relasi linier:
\[ y_i = \beta_0 + \beta_1 x_{1i} + \beta_2 x_{2i} + \varepsilon_i. \]
Konsep Model Regresi Linier dan Nonlinier
Regresi memetakan hubungan variabel dependen \(y\) dengan satu/lebih prediktor \(x\). Intuisi sederhananya: kita mencari garis/hiperbidang yang meringkas pola rerata \(y\) terhadap \(x\).
1.1 Regresi Linier (linier terhadap parameter)
Model klasik: \[ y_i = \beta_0 + \beta_1 x_{1i} + \cdots + \beta_k x_{ki} + \varepsilon_i,\quad E(\varepsilon_i\mid X)=0,\ \text{Var}(\varepsilon_i\mid X)=\sigma^2. \] “Linier” terhadap parameter berarti \(\beta\) hanya muncul pangkat satu. Bentuk fungsi terhadap \(x\) boleh fleksibel asalkan tetap linier pada \(\beta\): misalnya polinomial \(y=\beta_0+\beta_1 x+\beta_2 x^2+e\) tetap linier pada parameter.
Contoh fitur nonlinier tapi linier-dalam-parameter (basis):
- Polinomial: \(x, x^2, x^3\)
- Transformasi: \(\log x\), splines, dummy variabel (kategori)
- Interaksi: \(x_1\times x_2\)
1.2 Regresi Nonlinier (nonlinier terhadap parameter)
Bentuk umum: \(y_i = f(x_i,\,\beta) + \varepsilon_i\) dengan \(f\) nonlinier pada \(\beta\) (mis. eksponensial, Michaelis–Menten). Estimasi biasanya via Nonlinear Least Squares atau MLE khusus.
set.seed(1)
x <- seq(0, 4, length.out = 200)
# fungsi "benar" (tanpa noise)
f_lin <- 2 + 1.5*x
f_non <- 2 + exp(0.8*x)
# data observasi (dengan noise)
y_lin <- f_lin + rnorm(length(x), 0, 0.5)
y_non <- f_non + rnorm(length(x), 0, 0.5)
plot(x, y_lin, pch = 16, cex = 0.5, col = "grey50",
ylab = "y", main = "Relasi Linier vs Nonlinier")
points(x, y_non, pch = 1, cex = 0.5, col = "grey50")
lines(x, f_lin, lwd = 2) # garis mulus linier
lines(x, f_non, lwd = 2, lty = 2) # garis mulus eksponensial
legend("topleft",
c("Data linier (noisy)", "Data nonlinier (noisy)", "Linier (teoretis)", "Eksponensial (teoretis)"),
pch = c(16, 1, NA, NA), lty = c(NA, NA, 1, 2), lwd = c(NA, NA, 2, 2),
pt.cex = 0.8, bty = "n")
Prinsip kunci: jika hubungan tidak linier terhadap \(x\), sering kali rekayasa fitur (polinomial, log, spline) cukup tanpa pindah ke regresi nonlinier murni.
1.3 Metode Estimasi: OLS dan Maximum Likelihood
1.3.1 Ordinary Least Squares (OLS)
OLS memilih \(\hat\beta\) yang meminimalkan jumlah kuadrat residu: \[ \hat\beta = \arg\min_\beta (y-X\beta)’(y-X\beta) \Rightarrow \hat\beta=(X’X)^{-1}X’y. \] Di bawah asumsi klasik (Gauss–Markov: linearitas, eksogenitas, ragam konstan, tak ada multikolinearitas sempurna), \(\hat\beta\) adalah BLUE (Best Linear Unbiased Estimator).
1.3.2 Maximum Likelihood Estimation (MLE)
Asumsikan bahwa kesalahan berdistribusi normal:
\[
\varepsilon \sim N(0, \sigma^2 I).
\]
Likelihood function:
\[ L(\beta, \sigma^2) \propto (\sigma^2)^{-n/2} \exp\Big\{ -\tfrac{1}{2\sigma^2}(y - X\beta)'(y - X\beta) \Big\}. \]
Memaksimalkan log-likelihood menghasilkan
\[
\hat{\beta}_{\text{MLE}} = \hat{\beta}_{\text{OLS}}.
\]
Artinya, di bawah asumsi error normal, estimasi OLS juga merupakan
estimasi MLE;
perbedaan hanya muncul pada estimasi varians error
\(\sigma^2\), yaitu pembagi \(n\) vs \(n - k -
1\).
Catatan praktik:
Saat terjadi heteroskedastisitas atau autokorelasi, \(\hat{\beta}_{\text{OLS}}\) tetap tak
bias,
namun standar error tidak valid. Untuk mengoreksinya,
gunakan robust standard error
(seperti Heteroskedasticity-Consistent / HC, atau Newey–West)
atau gunakan pendekatan GLS/WLS.
1.3.3 Uji Hipotesis
1.3.3.1 Uji t untuk koefisien individual
Hipotesis: \(H_0: \beta_j = b_0\). Statistik: \[ t = \frac{\hat\beta_j - b_0}{SE(\hat\beta_j)} \sim t_{(n-k-1)}. \] P-value kecil \(\Rightarrow\) tolak \(H_0\). Interval kepercayaan \((1-\alpha)\%\): \(\hat\beta_j \pm t_{\alpha/2}\cdot SE\).
1.3.3.2 Uji F untuk restriksi simultan
Hipotesis linear: \(H_0: R\beta=r\) dengan \(q\) restriksi. Statistik F dari model terbatas vs penuh menilai signifikansi gabungan.
1.3.3.3 Uji kecocokan umum (overall)
\(H_0: \beta_1=\cdots=\beta_k=0\). Statistik F pada tabel
summary(lm) — menilai apakah prediktor bersama-sama
menjelaskan variasi \(y\).
1.4 Validitas Model: Uji Asumsi Klasik
Asumsi klasik memastikan inferensi sah: standar error, p-value, dan CI akurat.
- Linearitas: spesifikasi benar (dalam parameter), tidak ada mis‑spesifikasi sistematis.
- Normalitas error: terutama penting untuk inference pada sampel kecil.
- Homoskedastisitas: ragam residual konstan pada seluruh rentang fitted.
- Nonautokorelasi: residu tak berkorelasi (penting pada deret waktu/panel).
- Nonmultikolinearitas: tidak ada kolinearitas (atau minimal) antar prediktor.
Remedi umum bila asumsi gagal: transformasi variabel (log/Box‑Cox), tambah fitur (polinomial/interaksi), robust SE (HC0–HC3), Newey–West untuk autokorelasi, atau model alternatif (GLS, ARIMA errors).
1.5 Contoh Lengkap: Estimasi, Uji Hipotesis, Uji Asumsi
1.5.1 Bangun Data Simulasi
Kita buat data dengan dua prediktor, salah satunya sedikit kolinear dengan yang lain, plus heteroskedastisitas ringan agar diagnosa bermakna.
n <- 250
X1 <- rnorm(n, 50, 10)
X2_raw <- rnorm(n, 30, 5)
X2 <- 0.6*scale(X1) + 0.8*scale(X2_raw) # sengaja kolinear sedang
sigma <- 2 + 0.02*(X1 - mean(X1))^2 # heteroskedastisitas lemah
eps <- rnorm(n, 0, sigma)
Y <- 10 + 0.5*X1 - 0.3*X2 + eps
dat <- tibble(Y, X1, X2, sigma)
pairs(~ Y + X1 + X2, data = dat, main = "Scatterplot Matrix")
1.5.2 Estimasi Model OLS
fit <- lm(Y ~ X1 + X2, data = dat)
summary(fit)##
## Call:
## lm(formula = Y ~ X1 + X2, data = dat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -21.4946 -2.3506 0.0966 2.2356 17.7011
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 13.81273 1.80905 7.635 4.92e-13 ***
## X1 0.42005 0.03614 11.623 < 2e-16 ***
## X2 -0.32006 0.37585 -0.852 0.395
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.735 on 247 degrees of freedom
## Multiple R-squared: 0.4538, Adjusted R-squared: 0.4494
## F-statistic: 102.6 on 2 and 247 DF, p-value: < 2.2e-16confint(fit) # interval kepercayaan koefisien## 2.5 % 97.5 %
## (Intercept) 10.2495975 17.3758554
## X1 0.3488746 0.4912338
## X2 -1.0603371 0.4202134Interpretasi cepat.
- Tanda koefisien selaras dengan data generatif (\(+\) untuk \(X1\),
\(-\) untuk \(X2\)). - Perhatikan
Std. ErrordanPr(>|t|)untuk signifikansi. -R-squaredmengukur proporsi variasi \(Y\) yang dijelaskan.
1.5.3 Prediksi & Interval
newdf <- tibble(X1 = c(40, 60))
newdf$X2 <- matrix(c(25, 35), ncol = 1) # bikin nmatrix.1 juga
predict(fit, newdf, interval = "confidence")## fit lwr upr
## 1 22.61335 3.669021 41.55767
## 2 27.81381 2.368647 53.25898predict(fit, newdf, interval = "prediction")## fit lwr upr
## 1 22.61335 1.497742 43.72895
## 2 27.81381 0.713300 54.91432Catatan. Interval confidence untuk rerata bersifat lebih sempit daripada prediction untuk observasi baru.
1.5.4 Uji t & Uji F
# Uji t sudah ada di summary; contoh uji F gabungan:
# H0: beta_X1 = beta_X2 = 0
fit0 <- lm(Y ~ 1, data = dat)
anova(fit0, fit) # perbandingan model terbatas vs penuh (uji F)## Analysis of Variance Table
##
## Model 1: Y ~ 1
## Model 2: Y ~ X1 + X2
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 249 10139.1
## 2 247 5538.1 2 4601.1 102.6 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 11.5.5 Diagnostik Asumsi Klasik
1.5.5.1 Linearitas & Spesifikasi
Plot residual vs fitted semestinya acak; residualPlots
memeriksa nonlinieritas.
par(mfrow=c(1,2))
plot(fit, which = 1) # Residuals vs Fitted
car::residualPlots(fit) # Component + Residual (CERES) plots
## Test stat Pr(>|Test stat|)
## X1 1.0486 0.2954
## X2 -0.8251 0.4101
## Tukey test 1.1241 0.2610par(mfrow=c(1,1))1.5.5.2 Normalitas Error
Gunakan QQ-plot dan uji formal (Shapiro & Anderson–Darling).
plot(fit, which = 2) # QQ-plot
shapiro.test(residuals(fit)) # sensitif untuk n besar##
## Shapiro-Wilk normality test
##
## data: residuals(fit)
## W = 0.94741, p-value = 7.73e-08nortest::ad.test(residuals(fit))##
## Anderson-Darling normality test
##
## data: residuals(fit)
## A = 2.9442, p-value = 2.043e-071.5.5.3 Homoskedastisitas
Breusch–Pagan (bentuk dasar & White-style).
lmtest::bptest(fit) # BP test (fitted)##
## studentized Breusch-Pagan test
##
## data: fit
## BP = 10.255, df = 2, p-value = 0.005932lmtest::bptest(fit, ~ fitted(fit) + I(fitted(fit)^2)) # White-type##
## studentized Breusch-Pagan test
##
## data: fit
## BP = 58.132, df = 2, p-value = 2.381e-13plot(fit, which = 3) # Scale-Location plot
Robust SE (HC) bila heteroskedastik:
coeftest(fit, vcov = vcovHC(fit, type = "HC3"))##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 13.812726 2.370668 5.8265 1.758e-08 ***
## X1 0.420054 0.049091 8.5567 1.246e-15 ***
## X2 -0.320062 0.406506 -0.7873 0.4318
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 11.5.5.4 Nonautokorelasi (opsional di cross-section)
Durbin–Watson dan Breusch–Godfrey.
lmtest::dwtest(fit)##
## Durbin-Watson test
##
## data: fit
## DW = 1.953, p-value = 0.3571
## alternative hypothesis: true autocorrelation is greater than 0lmtest::bgtest(fit, order = 2) # uji autokorelasi orde 2##
## Breusch-Godfrey test for serial correlation of order up to 2
##
## data: fit
## LM test = 0.77854, df = 2, p-value = 0.6776Untuk deret waktu, pertimbangkan menambah lag residual atau gunakan Newey–West SE:
coeftest(fit, vcov = NeweyWest(fit, lag = 2, prewhite = FALSE))1.5.5.5 Nonmultikolinearitas
VIF dan condition number.
car::vif(fit)## X1 X2
## 1.663999 1.663999X <- model.matrix(fit)
kappa(qr.R(qr(X))) # condition number (>>30 indikasi masalah)## [1] 273.1844Remedi: buang/rekayasa fitur, regularisasi (Ridge/Lasso), atau kumpulkan lebih banyak data.
1.5.5.6 Outlier & Leverage
Cook’s distance, leverage, dan pengaruh.
par(mfrow=c(1,2))
plot(fit, which = 4) # Cook's distance
plot(fit, which = 5) # Residuals vs Leverage
par(mfrow=c(1,1))
influence.measures(fit)$infmat[1:5,]## dfb.1_ dfb.X1 dfb.X2 dffit cov.r cook.d
## 1 0.051030617 -0.065363990 -0.002496740 -0.11851524 1.000909 4.670010e-03
## 2 -0.085733056 0.095141837 -0.040630593 0.10988430 1.025760 4.031613e-03
## 3 0.009082744 -0.008283275 0.009692429 0.01149691 1.030027 4.423737e-05
## 4 0.013489488 -0.019628988 0.044407632 -0.05785212 1.019275 1.118744e-03
## 5 0.010630488 -0.007819191 0.005232231 0.01929349 1.016156 1.245448e-04
## hat
## 1 0.00853456
## 2 0.02019325
## 3 0.01735334
## 4 0.01065226
## 5 0.004788311.6 Estimasi dengan MLE (verifikasi ringkas)
Dengan asumsi normalitas error, OLS = MLE untuk \(\beta\). Kita verifikasi dengan optimisasi numerik sederhana.
negloglik <- function(par, y, X){
beta <- par[1:ncol(X)]
logs2 <- par[length(par)]
s2 <- exp(logs2) # pastikan positif
e <- y - X %*% beta
n <- length(y)
0.5*n*log(2*pi) + 0.5*n*log(s2) + 0.5*sum(e^2)/s2
}
X <- model.matrix(fit)
y <- dat$Y
init <- c(coef(fit), log(var(residuals(fit))))
opt <- optim(init, negloglik, y=y, X=X, method="BFGS", hessian=TRUE)
beta_mle <- opt$par[1:ncol(X)]
cbind(OLS=coef(fit), MLE=beta_mle)## OLS MLE
## (Intercept) 13.8127264 13.8127264
## X1 0.4200542 0.4200542
## X2 -0.3200619 -0.3200619Hasil \(\beta\) dari MLE praktis sama dengan OLS (perbedaan numerik kecil).
1.7 Rangkuman & Praktik Baik
- Gunakan grafik diagnostik terlebih dahulu; uji formal melengkapi, bukan menggantikan.
- Saat asumsi gagal: gunakan robust SE (HC3) untuk heteroskedastisitas, Newey–West untuk autokorelasi, atau pertimbangkan spesifikasi ulang.
- Dokumentasikan keputusan model (transformasi, outlier treatment, pemilihan fitur) agar reproducible.
1.8 Lampiran: Kode Ringkas One‑Shot
fit <- lm(Y ~ X1 + X2, data = dat)
list(
summary = summary(fit),
confint = confint(fit),
t_overallF = anova(lm(Y ~ 1, dat), fit),
bp = bptest(fit),
bp_white = bptest(fit, ~ fitted(fit) + I(fitted(fit)^2)),
dw = dwtest(fit),
bg = bgtest(fit, order = 2),
vif = car::vif(fit),
robust = coeftest(fit, vcov = vcovHC(fit, type = "HC3"))
)Penutup. OLS tetap rajanya untuk model linier—cepat, transparan, dan kuat—asal rumah tangga asumsi rapi. Saat dunia nyata berisik, naikkan perisai robust SE dan akal sehat spesifikasi.
n <- 200
x1 <- rnorm(n)
x2 <- runif(n, -1, 1)
y <- 3 + 2*x1 - 1.5*x2 + rnorm(n, sd = 0.6)
ols_model <- lm(y ~ x1 + x2)
summary(ols_model)##
## Call:
## lm(formula = y ~ x1 + x2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.78896 -0.34550 0.04367 0.38220 1.33299
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.99707 0.04253 70.47 <2e-16 ***
## x1 1.99354 0.04074 48.94 <2e-16 ***
## x2 -1.53342 0.07308 -20.98 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6001 on 197 degrees of freedom
## Multiple R-squared: 0.9344, Adjusted R-squared: 0.9338
## F-statistic: 1404 on 2 and 197 DF, p-value: < 2.2e-16Catatan interpretasi. Koefisien menyatakan perubahan rata-rata pada \(y\) akibat perubahan satu unit pada kovariat terkait, ceteris paribus.
2 Endogeneity dalam Ekonometrika
2.1 Model Regresi Dasar dan Asumsi Eksogenitas
Kita mulai dari model linier sederhana:
\[ Y_i = \beta_0 + \beta_1 X_i + \epsilon_i \]
dengan
- \(Y_i\): variabel dependen (outcome),
- \(X_i\): variabel penjelas / regressor,
- \(\epsilon_i\): error term (faktor tak teramati yang ikut memengaruhi \(Y_i\)).
Estimasi biasa dilakukan dengan Ordinary Least Squares (OLS).
Agar OLS menghasilkan estimator tak-bias (unbiased) dan
konsisten, kita butuh asumsi kunci yang disebut
eksogenitas:
\[ \mathbb{E}[\epsilon_i \mid X_i] = 0 \quad \text{atau ekuivalen} \quad \text{Cov}(X_i, \epsilon_i)=0. \]
Artinya:
- tidak ada korelasi antara variabel penjelas \(X_i\) dan komponen gangguan \(\epsilon_i\),
- semua variasi \(X_i\) dianggap “bersih”, tidak terkontaminasi faktor tak teramati di dalam \(\epsilon_i\).
Inilah kondisi ideal yang diajarkan di buku teks ekonometrika intro: dunia patuh, data jinak, dosen bahagia.
2.2 Definisi Endogeneity
Endogeneity terjadi ketika asumsi eksogenitas dilanggar, yaitu:
\[ \text{Cov}(X_i, \epsilon_i) \neq 0. \]
Secara verbal:
Variabel penjelas \(X_i\) “tercemar” oleh sesuatu yang juga ada di error \(\epsilon_i\), sehingga \(X_i\) sebagian mengandung sinyal yang bukan murni sebab \(Y_i\), tapi justru “reaksi bersama” dari mekanisme lain.
Akibatnya, OLS tidak lagi bisa diperlakukan sebagai estimator hubungan kausal \(X \rightarrow Y\).
Kita menyebut \(X_i\) sebagai endogenous regressor (regresor endogen).
2.3 Sumber-Sumber Umum Endogeneity
(a) Omitted Variable Bias (Variabel Terlewat) Ada variabel penting yang memengaruhi \(Y_i\) tetapi tidak masuk model, dan variabel itu juga berkorelasi dengan \(X_i\).
Contoh:
- Model: gaji (Y) ~ pendidikan (X).
- Variabel kemampuan kognitif bawaan / skill sosial masuk ke \(\epsilon_i\), tidak diukur.
- Orang dengan kemampuan tinggi: cenderung sekolah lebih lama dan bergaji lebih tinggi.
- Jadi pendidikan \(X\) berkorelasi dengan \(\epsilon_i\) → endogeneity.
Secara skematik:
\[ \text{Kemampuan} \rightarrow \text{Pendidikan} \rightarrow \text{Gaji} \] dan \[ \text{Kemampuan} \rightarrow \text{Gaji (langsung)}. \]
Kemampuan ada di error. Masalah.
(b) Simultanitas / Reverse Causality \(X_i\) memengaruhi \(Y_i\), tetapi \(Y_i\) juga memengaruhi \(X_i\). Hubungan dua arah.
Contoh klasik penawaran–permintaan:
- Harga (P) dan Kuantitas (Q) saling ditentukan secara simultan.
- Kalau kita regresikan Q pada P secara naif, error dalam persamaan permintaan akan berkorelasi dengan P karena P sendiri hasil ekuilibrium pasar, bukan variabel bebas murni.
Ini juga muncul di:
- produktivitas perusahaan ↔︎ upah pekerja,
- kesehatan ↔︎ pendapatan rumah tangga.
Measurement Error (Kesalahan Pengukuran) Variabel \(X_i\) dicatat dengan noise.
Misal kita ingin “jam belajar”, tapi data berasal dari self-report siswa yang cenderung melebihkan demi terlihat rajin. Distorsi ini masuk ke error dan membuat OLS bias.
Dynamic Panel Bias (Nickell Bias) Dalam model panel dinamis:
\[ Y_{it} = \alpha Y_{i,t-1} + \beta X_{it} + \mu_i + \varepsilon_{it}, \]
lag dependen \(Y_{i,t-1}\) berkorelasi dengan efek tetap individu \(\mu_i\), sehingga endogen.
2.4 Konsekuensi Endogeneity
2.4.1 Bias
Estimator OLS untuk \(\beta_1\) jadi bias:
\[ \hat{\beta}_1^{OLS} \xrightarrow{E} \beta_1 + \frac{\text{Cov}(X_i,\epsilon_i)}{\text{Var}(X_i)}. \]
2.4.2 Inkonsistensi
Bahkan dengan n → ∞, estimator tetap tidak mendekati nilai benar.
2.4.3 Inferensi Salah
Uji-t dan p-value tetap bisa terlihat “signifikan” padahal model salah spesifikasi. Ini ilusi presisi.
2.5 Strategi Solusi
2.5.1 Tambahkan Variabel Kontrol yang Hilang
Masukkan variabel yang relevan agar sumber korelasi berkurang.
2.6 4.4 Instrumental Variables (IV) / Two-Stage Least Squares (2SLS)
Gunakan variabel instrumen \(Z\) yang relevan dan eksogen:
Tahap 1: \[ X_i = \pi_0 + \pi_1 Z_i + v_i \] Tahap 2: \[ Y_i = \beta_0 + \beta_1 \hat{X}_i + e_i \]
2.6.1 Fixed Effects (FE) pada Data Panel
Hilangkan efek individu konstan \(\mu_i\) dengan transformasi within.
2.6.2 Difference-in-Differences (DiD)
Gunakan perubahan sebelum–sesudah antar kelompok untuk mengidentifikasi efek kausal.
2.6.3 GMM Arellano–Bond
Gunakan lag terdalam sebagai instrumen internal untuk model panel dinamis.
2.7 Ringkasan Konseptual
- Asumsi ideal OLS: \(\text{Cov}(X,\epsilon)=0\).
- Endogeneity = pelanggaran asumsi itu.
- Konsekuensi: bias, inkonsistensi, inferensi menipu.
- Solusi: kontrol, fixed effects, DiD, IV, GMM.
Endogeneity bukan hanya masalah teknis, tapi juga persoalan filsafat kausalitas: setiap solusi menuntut asumsi baru yang harus dijustifikasi secara ilmiah.
3 Instrumental Variable (IV) Regression
3.1 Definisi & Aplikasi
Instrumental Variable (IV) Regression digunakan ketika variabel penjelas utama endogen, yaitu berkorelasi dengan error (\(\text{Cov}(x,u)\neq 0\)). Endogenitas muncul karena omitted variables, measurement error, atau simultaneity (keterbalikan kausal).
3.2 Mengapa Instrumental Variable (IV) Diperlukan?
Bayangkan kita ingin mengestimasi pengaruh pendidikan \(X\) terhadap pendapatan \(Y\). Masalahnya, pendidikan sering berkorelasi dengan hal-hal yang tak teramati, seperti kemampuan, motivasi, atau lingkungan keluarga \(u\). Model sederhananya:
\[ Y_i = \beta_0 + \beta_1 X_i + \epsilon_i. \]
Idealnya, \(\mathrm{Cov}(X_i, \epsilon_i) = 0\). Namun jika ternyata \(\mathrm{Cov}(X_i, \epsilon_i) \neq 0\), maka asumsi OLS rusak. Akibatnya:
- Koefisien \(\hat{\beta}_1\) menjadi bias (menyimpang dari nilai sebenarnya).
- Dan inkonsisten (tidak akan mendekati nilai benar meski ukuran sampel diperbesar).
Sumber endogenitas dapat berasal dari:
- Omitted variable (ada variabel penting tak dimasukkan ke model),
- Measurement error (kesalahan pengukuran pada \(X\)),
- Simultaneity (hubungan dua arah: \(Y\) memengaruhi \(X\) juga).
Untuk itulah muncul Instrumental Variable (IV) — alat bantu untuk “menyaring” pengaruh murni \(X\) yang benar-benar eksogen.
Apa itu variabel instrumen?
Variabel instrumen adalah “solusi kreatif” untuk masalah ini: kita cari variabel \(Z\) yang:
- Terkorelasi dengan \(X\) (artinya:
\(Z\) mempengaruhi atau berkaitan
dengan \(X\)) → disebut
relevansi.
- Tidak berkorelasi dengan kesalahan \(u\) dan tidak mempengaruhi \(Y\) kecuali melalui pengaruhnya ke \(X\) → disebut eksogenitas atau exclusion restriction.
Jadi inti: \(Z\) membantu kita mendapatkan “variasi murni” dalam \(X\) yang bebas dari bias karena \(u\). Dengan variasi ini kita bisa mengestimasi \(\beta\) dengan lebih kredibel
Contoh sederhana
Misalnya kita ingin tahu: apakah merokok \((X)\) menyebabkan kesehatan menurun \((Y)\)?
Tapi kita tahu: banyak variabel tersembunyi seperti “stress”,
“kepribadian”, “kondisi sosio-ekonomi” yang bisa mempengaruhi baik
merokok maupun kesehatan → masalah endogenitas.
Kita bisa memilih instrumen misalnya: tarif cukai rokok di suatu
wilayah \((Z)\). Alasan:
- Tarif cukai mempengaruhi kemungkinan seseorang merokok (relevan).
:contentReferenceoaicite:1
- Kita anggap tarif cukai tidak langsung mempengaruhi kesehatan
kecuali melalui merokok (validitas eksklusinya) — ini asumsi
besar!
- Jika asumsi ini kira-kira benar, maka kita bisa menggunakan tarif
cukai sebagai instrumen untuk merokok dan mengestimasi efek merokok
terhadap kesehatan. :contentReferenceoaicite:2
(Disclaimer: dalam praktik, seringkali asumsi eksklusinya diperdebatkan — karena tarif bisa pula mempengaruhi kesehatan melalui jalur lain, misalnya melalui pengeluaran rumah tangga, stress fiskal, dsb.)
Catatan
- Jika instrumen \(Z\) tidak
cukup kuat dalam mempengaruhi \(X\) (instrumen lemah), maka estimasi IV
bisa punya varians besar atau bahkan bias besar.
- Asumsi bahwa \(Z\) tidak
memiliki jalur langsung ke \(Y\)
selain lewat \(X\) (eksclusion)
tidak bisa sepenuhnya diuji secara statistik — banyak
bergantung pada argumen substantif dari konteks.
- Estimasi IV biasanya memberikan efek untuk sub-populasi yang memang terpengaruh oleh instrumen (termasuk konsep LATE — Local Average Treatment Effect) bukan otomatis efek rata-rata seluruh populasi.
Mengapa ini penting? Karena dalam banyak
penelitian, kita ingin tahu “apa yang terjadi jika kita
mengubah \(X\)?” — itu
pertanyaan sebab-akibat.
Tapi banyak gap: orang memilih sendiri, ada faktor tersembunyi, penyakit
bisa mempengaruhi pilihan, dsb. Instrumen adalah seperti “eksternal
shove” ke \(X\) yang kita harap cukup
acak untuk memberi kredibilitas sebab-akibatnya.
Bayangkan: kamu ingin tahu apakah tidur lebih lama bikin nilai kuliah
kamu naik, tapi tentu banyak faktor (motivasi, kesehatan,
murid-lebih-bling). Jika suatu universitas mengubah kebijakan jam tidur
mahasiswa secara acak karena pembangunan kampus — itu bisa jadi
instrumen. (Ya, agak absurd tapi you get the idea.)
Catatan:
- Variabel instrumen = \(Z\).
- Tujuannya: membantu kita mengestimasi efek sebab-akibat dari \(X\) ke \(Y\) ketika \(X\) bermasalah (endogen).
- Syarat utama:
- \(Z\) berkorelasi dengan \(X\) (relevansi).
- \(Z\) tidak berkorelasi dengan
kesalahan \(u\) dan pengaruh \(Z\) ke \(Y\) hanya lewat \(X\) (eksclusion).
- \(Z\) berkorelasi dengan \(X\) (relevansi).
- Digunakan via 2-step (2SLS) dalam banyak aplikasi.
- Harus diingat: asumsi kuat, instrumen sulit ditemukan, dan interpretasi harus hati-hati.
Ringkasnya: OLS menjawab korelasi, IV berupaya mendekati kausalitas ketika ada endogenitas.
Catatan Instrumen \(z\) diperlukan untuk mengisolasi variasi \(x\) yang eksogen terhadap error:
- Relevansi: \(\text{Cov}(z, x) \neq 0\) — instrumen harus kuat memprediksi \(x\).
- Validitas (Eksogenitas): \(\text{Cov}(z, u) = 0\) — instrumen tidak memengaruhi \(y\) kecuali melalui \(x\).
Aplikasi umum: pengaruh pendidikan terhadap upah (instrumen: jarak ke kampus), permintaan–penawaran (harga endogen; instrumen: shock biaya), kebijakan publik (intensitas program, kuota, lotere), kesehatan (ketersediaan fasilitas).
Intuisi: IV “meminjam” variasi \(x\) yang datang dari \(z\) (yang bersih dari noise endogen), sehingga estimasi efek kausal \(x \to y\) menjadi konsisten.
3.3 Detail: apa Itu Instrumen?
Sebuah variabel \(Z\) disebut instrumen yang valid jika memenuhi dua syarat pokok:
- Relevansi: \(\mathrm{Cov}(Z, X) \neq 0\) (instrumen mendorong variasi pada \(X\)).
- Eksogenitas: \(\mathrm{Cov}(Z, u) = 0\) (instrumen tidak berkorelasi dengan error).
Dengan kata lain, \(Z\) memengaruhi \(Y\) hanya melalui \(X\), bukan lewat jalur lain. Secara intuitif: instrumen harus “mendorong” \(X\), tetapi tidak “menyentuh” \(Y\) secara langsung.
3.4 Struktur Model (Structural Equation) dan 2SLS
Model dengan variabel endogen disebut structural model. Misalnya:
\[ Y = \beta_0 + \beta_1 X + u, \]
dengan \(X\) endogen. Karena \(X\) berkorelasi dengan \(u\), OLS tidak dapat dipakai secara valid. Solusi umum adalah Two-Stage Least Squares (2SLS) menggunakan instrumen \(Z\).
Tahap 1 (First Stage): bangun reduced form untuk \(X\) \[ X = \pi_0 + \pi_1 Z + v \quad \Rightarrow \quad \hat{X} = \hat{\pi}_0 + \hat{\pi}_1 Z. \]
Tahap 2 (Second Stage): gunakan \(\hat{X}\) pada persamaan struktural \[ Y = \beta_0 + \beta_1 \hat{X} + e. \]
Secara aljabar, 2SLS/IV dapat ditulis sebagai \[ \hat{\beta}_{IV} = (X'P_Z X)^{-1} X' P_Z Y, \quad \text{dengan } P_Z = Z (Z'Z)^{-1} Z'. \]
Interpretasi: \(P_Z\) mem-proyeksi \(X\) ke ruang yang dijelaskan oleh \(Z\) (bagian yang eksogen), lalu efeknya diukur pada \(Y\).
Catatan praktik: kekuatan instrumen penting untuk diperiksa (mis. first-stage F-statistic > 10 sebagai aturan praktis).
3.5 Inti Filosofinya
Instrumen bekerja seperti filter kausalitas: memisahkan sinyal murni \(X\) dari gangguan \(u\). Tanpa instrumen, kita mudah terjebak pada interpretasi korelasional belaka. Dengan instrumen, kita punya pegangan untuk mengatakan “bagian dari \(X\) yang didorong oleh \(Z\) benar-benar memengaruhi \(Y\).”
IV bukan tongkat sihir: ia membantu ketika instrumen valid dan kuat. Diagnosis (relevansi dan eksogenitas) adalah bagian dari sains, bukan sekadar ritual.
3.6 Detail Model
Model struktural sederhana (satu regressor endogen):
\[ \begin{aligned} \text{(Struktural)}\quad & y = \beta_0 + \beta_1 x + w'\gamma + u, \qquad \text{dengan } \mathrm{Cov}(x,u)\neq 0, \\\\ \text{(First stage)}\quad & x = \pi_0 + \pi_1 z + w'\delta + v, \qquad \text{dengan } \mathrm{Cov}(z,u)=0,~\mathrm{Cov}(z,v)\neq 0, \end{aligned} \]
di mana \(w\) adalah kovariat eksogen (kontrol).
Identifikasi:
- Just-identified: jumlah instrumen eksklusif (di luar
\(w\)) = jumlah regressor
endogen.
- Overidentified: instrumen lebih banyak dari endogen
(memungkinkan uji validitas — Sargan/Hansen).
Target estimand: pada model linear dengan kesetaraan struktural, IV menaksir ATE (average treatment effect) jika seluruh asumsi terpenuhi. Pada kasus biner dengan kepatuhan heterogen dan monotonicity, interpretasi sering kali LATE (local average treatment effect).
3.7 Metode Estimasi: 1SLS & 2SLS
3.7.1 One-Stage Least Squares (1SLS) — Estimator IV langsung
Secara aljabar: \[ \hat\beta^{IV} = (X'P_Z X)^{-1} X'P_Z y, \qquad P_Z = Z\,(Z'Z)^{-1}Z'. \] Di sini \(X\) memuat kolom endogen dan eksogen, \(Z\) memuat semua instrumen (instrumen eksklusif + eksogen). 1SLS menghitung proyeksi \(P_Z y\) dan \(P_Z X\) lalu OLS pada data terproyeksi. Pada kasus just-identified, 1SLS identik dengan 2SLS.
3.7.2 Two-Stage Least Squares (2SLS)
Tahap-1 (first stage): regresikan setiap regressor
endogen pada seluruh instrumen \(Z\)
dan eksogen \(w\); simpan nilai
terprediksi \(\hat x\).
Tahap-2: regresikan \(y\) pada \(\hat
x\) dan \(w\) dengan OLS.
Secara numerik, 2SLS \(=\) 1SLS;
perbedaan hanya pada cara implementasi.
Diagnostik penting:
- Uji kelemahan instrumen: F-statistik tahap-1
(rule-of-thumb \(F>10\)).
- Uji overidentifikasi: Sargan atau Hansen J (jika
overidentified).
- Uji endogenitas: Wu–Hausman (apakah OLS vs IV berbeda
signifikan).
- SE robust: gunakan HC
(heteroskedasticity-consistent) karena IV sensitif terhadap
heteroskedastisitas.
3.7.3 Ilustrasi Detail di R
Kita simulasikan endogenitas: \(x\) mengandung komponen error struktural \(u\) (melalui \(v\)). Instrumen \(z\) memengaruhi \(x\) tapi tidak \(y\) secara langsung.
n <- 1000
# Kovariat eksogen (kontrol)
w1 <- rnorm(n)
w2 <- rbinom(n, 1, 0.4)
# Instrumen eksklusif yang relevan dan valid
z <- 0.5*rnorm(n) + 0.5*w1 + rnorm(n, sd = 0.5)
# Error struktural dan korelasi endogenitas
u <- rnorm(n)
v <- 0.6*u + rnorm(n) # membuat x berkorelasi dg u via v
# First stage: x dipengaruhi instrumen z, kontrol w, dan v yang kor. dg u
x <- 0.8 + 0.9*z + 0.5*w1 - 0.4*w2 + v
# Structural: y dipengaruhi x, kontrol w, dan u
y <- 1.5 + 1.2*x + 0.7*w1 - 0.5*w2 + u
dat <- tibble(y, x, z, w1, w2 = factor(w2))3.7.4 OLS (bias karena endogenitas)
ols_fit <- lm(y ~ x + w1 + w2, data = dat)
summary(ols_fit)##
## Call:
## lm(formula = y ~ x + w1 + w2, data = dat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.0455 -0.6060 -0.0229 0.6353 3.2190
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.24126 0.04258 29.150 < 2e-16 ***
## x 1.49459 0.02220 67.313 < 2e-16 ***
## w1 0.43364 0.03477 12.470 < 2e-16 ***
## w21 -0.29325 0.05980 -4.904 1.1e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.92 on 996 degrees of freedom
## Multiple R-squared: 0.9005, Adjusted R-squared: 0.9002
## F-statistic: 3005 on 3 and 996 DF, p-value: < 2.2e-16coeftest(ols_fit, vcov = vcovHC(ols_fit, type = "HC3")) # SE robust##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.241256 0.040811 30.4151 < 2.2e-16 ***
## x 1.494595 0.021420 69.7742 < 2.2e-16 ***
## w1 0.433642 0.034569 12.5443 < 2.2e-16 ***
## w21 -0.293248 0.059351 -4.9409 9.118e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Karena \(x\) endogen, koefisien OLS untuk \(x\) cenderung bias. Bandingkan nanti dengan IV.
3.7.5 2SLS dengan
AER::ivreg
Spesifikasi formula: y ~ x + w1 + w2 | z + w1 + w2. Di
kanan tanda | taruh semua instrumen
(instrumen eksklusif + eksogen yang selalu valid).
iv_2sls <- AER::ivreg(y ~ x + w1 + w2 | z + w1 + w2, data = dat)
summary(iv_2sls, diagnostics = TRUE) # tampilkan uji endogenitas, weak IV, overid (jika ada)##
## Call:
## AER::ivreg(formula = y ~ x + w1 + w2 | z + w1 + w2, data = dat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.42894 -0.71128 -0.02568 0.76344 3.88219
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.60376 0.06523 24.587 < 2e-16 ***
## x 1.04099 0.05773 18.032 < 2e-16 ***
## w1 0.85163 0.06286 13.547 < 2e-16 ***
## w21 -0.52806 0.07603 -6.946 6.79e-12 ***
##
## Diagnostic tests:
## df1 df2 statistic p-value
## Weak instruments 1 996 264.6 <2e-16 ***
## Wu-Hausman 1 995 124.6 <2e-16 ***
## Sargan 0 NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.096 on 996 degrees of freedom
## Multiple R-Squared: 0.8588, Adjusted R-squared: 0.8584
## Wald test: 1162 on 3 and 996 DF, p-value: < 2.2e-16coeftest(iv_2sls, vcov = vcovHC(iv_2sls, type = "HC3"))##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.603762 0.065725 24.4012 < 2.2e-16 ***
## x 1.040992 0.060797 17.1225 < 2.2e-16 ***
## w1 0.851628 0.064668 13.1693 < 2.2e-16 ***
## w21 -0.528062 0.075412 -7.0024 4.622e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Membaca diagnostics
(summary(..., diagnostics=TRUE)):
- Weak instruments: lihat “Weak instruments” /
“First-stage F”.
- Wu–Hausman: menguji endogenitas \(x\).
- Sargan (jika ada >1 instrumen eksklusif): menguji
validitas gabungan instrumen.
3.7.6 1SLS manual (proyeksi dengan matriks)
Konstruksi \(Z\) dan \(X\) lalu hitung \(\hat\beta^{IV}=(X'P_ZX)^{-1}X'P_Z y\).
# Desain matriks: X (konstanta, x endogen, kontrol), Z (konstanta, instrumen, kontrol)
X <- model.matrix(~ x + w1 + w2, data = dat) # termasuk intercept
Z <- model.matrix(~ z + w1 + w2, data = dat)
PZ <- Z %*% solve(t(Z) %*% Z) %*% t(Z)
beta_iv <- solve(t(X) %*% PZ %*% X) %*% (t(X) %*% PZ %*% dat$y)
colnames(beta_iv) <- "coef_IV_1SLS"
beta_iv## coef_IV_1SLS
## (Intercept) 1.6037624
## x 1.0409924
## w1 0.8516277
## w21 -0.5280620Untuk SE pada 1SLS manual, gunakan sandwich formula IV (tidak seremeh OLS). Praktisnya, gunakan
ivreg/ivfixedyang sudah mengimplementasikan varian robust.
3.7.7 2SLS manual (dua tahap)
Tahap-1: prediksi \(\hat x\) dari \((z, w)\). Tahap-2: regresikan \(y\) pada \(\hat x\) dan \(w\) dengan OLS.
# Tahap-1
fs_fit <- lm(x ~ z + w1 + w2, data = dat)
summary(fs_fit) # cek F-stat (kekuatan instrumen)##
## Call:
## lm(formula = x ~ z + w1 + w2, data = dat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.5840 -0.8050 0.0564 0.7821 3.2299
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.79908 0.04910 16.274 < 2e-16 ***
## z 0.85511 0.05257 16.267 < 2e-16 ***
## w1 0.50354 0.04396 11.455 < 2e-16 ***
## w21 -0.47685 0.07448 -6.402 2.36e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.167 on 996 degrees of freedom
## Multiple R-squared: 0.4963, Adjusted R-squared: 0.4948
## F-statistic: 327.2 on 3 and 996 DF, p-value: < 2.2e-16xhat <- fs_fit$fitted.values
# Tahap-2
tsls_fit <- lm(y ~ xhat + w1 + w2, data = dat)
summary(tsls_fit)##
## Call:
## lm(formula = y ~ xhat + w1 + w2, data = dat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.8724 -1.4524 0.0327 1.4399 6.5000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.6038 0.1235 12.987 < 2e-16 ***
## xhat 1.0410 0.1093 9.525 < 2e-16 ***
## w1 0.8516 0.1190 7.156 1.61e-12 ***
## w21 -0.5281 0.1439 -3.669 0.000257 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.075 on 996 degrees of freedom
## Multiple R-squared: 0.4941, Adjusted R-squared: 0.4925
## F-statistic: 324.2 on 3 and 996 DF, p-value: < 2.2e-16# Catatan: SE di tahap-2 OLS biasa tidak tepat untuk 2SLS.
# Solusi praktis: pakai ivreg (sudah menghitung SE 2SLS yang benar/robust).3.7.8 Prediksi & Efek Kausal
IV fokus pada efek kausal rata-rata dari \(x\) terhadap \(y\) (di bawah asumsi). Prediksi out-of-sample bisa dilakukan, tetapi ingat bahwa IV menekankan identifikasi kausal, bukan akurasi prediksi semata.
cbind(
OLS = coef(ols_fit)["x"],
IV_2SLS = coef(iv_2sls)["x"]
)## OLS IV_2SLS
## x 1.494595 1.040992Bandingkan nilai koefisien
x: OLS vs IV. OLS bias menuju arah korelasi \(x\)–\(u\). IV mendekati nilai benar (1.2) karena variasi dari instrumen bebas dari \(u\).
3.8 Uji-Uji Tambahan
- Kekuatan instrumen: Targetkan first-stage \(F>10\) (rule-of-thumb). Instrumen lemah
→ estimasi IV inflasi varians & bias finite-sample.
- Validitas instrumen: Uji overidentifikasi
(Sargan/Hansen) bila overidentified. Ingat, uji ini tidak memvalidasi
instrumen secara absolut — hanya konsistensi internal.
- Endogenitas: Wu–Hausman membandingkan OLS vs IV.
Jika H0 ditolak, gunakan IV.
- Robust SE: Selalu laporkan SE robust (HC).
- Interpretasi: Pada instrumen biner dengan kepatuhan heterogen, IV menaksir LATE pada compliers (di bawah monotonicity).
3.9 Ringkasan
- IV menyelesaikan endogenitas dengan instrumen yang
relevan dan valid.
- 1SLS (rumus proyeksi) dan 2SLS
(dua tahap) setara secara numerik;
AER::ivregmemudahkan implementasi dan menyediakan diagnostics.
- Selalu cek first-stage F, uji overid, dan SE robust. Dokumentasikan asumsi dan justifikasi instrumen secara substantif.
Diperlukan saat variabel penjelas endogen (berkorelasi dengan error). Gunakan instrumen \(z\) yang relevan dan valid.
Model struktural:
\[ y_i = \beta_0 + \beta_1 x_i + \epsilon_i, \quad x_i \text{ endogen}. \]
n <- 500
u <- rnorm(n)
z <- rnorm(n) # instrumen
x <- 0.9*z + 0.5*u + rnorm(n) # endogen (terkait u)
y <- 2 + 1.5*x + u
iv_model <- AER::ivreg(y ~ x | z)
summary(iv_model, diagnostics = TRUE) # tampilkan diagnosa instrumen##
## Call:
## AER::ivreg(formula = y ~ x | z)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.30779 -0.67530 0.04916 0.73319 2.60398
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.00524 0.04580 43.78 <2e-16 ***
## x 1.45045 0.05816 24.94 <2e-16 ***
##
## Diagnostic tests:
## df1 df2 statistic p-value
## Weak instruments 1 498 223.11 < 2e-16 ***
## Wu-Hausman 1 497 50.77 3.67e-12 ***
## Sargan 0 NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.022 on 498 degrees of freedom
## Multiple R-Squared: 0.8508, Adjusted R-squared: 0.8505
## Wald test: 622 on 1 and 498 DF, p-value: < 2.2e-16Ringkas. IV mereduksi bias kausal saat ada endogenitas; cek kekuatan instrumen (uji F first-stage) dan validitas (overidentification bila (>)1 instrumen).
4 Regresi dengan Variabel Dummy
Variabel dummy (juga disebut indicator atau binary variables) adalah variabel bernilai 0/1 yang kita gunakan untuk memasukkan informasi kategorikal ke dalam model regresi linier. Dengan dummy, kita bisa memodelkan:
- Varying intercept (intersep berbeda): garis regresi paralel untuk setiap kelompok.
- Varying slope (kemiringan berbeda): respons terhadap perubahan kovariat kontinu berbeda antar kelompok.
- Interaksi: pengaruh suatu kovariat bergantung pada tingkat kovariat lain, baik keduanya dummy, ataupun kombinasi dummy–kontinu, bahkan antar banyak kategori.
Tujuan dokumen ini: memberi panduan praktis—lengkap dengan
data simulasi, kode R, dan
visualisasi—untuk semua situasi umum pemodelan dengan
variabel dummy. Kita akan memakai sintaks lm() agar
interpretasi tetap jernih. Untuk visualisasi kita andalkan
ggplot2.
Kenapa perlu dummy? Karena dunia nyata jarang 100% kontinu. Kebijakan, jenis kelamin, perlakuan eksperimen, cabang toko, semester, provinsi—semuanya kategori. Tanpa dummy, model tidak bisa membedakan kelompok dan akan “meratakan” efeknya jadi satu garis saja. Dengan dummy, kita memberi hak suara kepada tiap kelompok.
Catatan gaya: kita menjaga nada ringkas dan lugas, sesekali menyelipkan humor hemat—secukupnya untuk membuat konsep lengket di kepala tanpa membuatnya lengket di keyboard.
4.1 Paket dan Data Simulasi
Kita mulai dengan memuat paket dan membuat satu taman bermain
data yang kaya: ada satu kovariat kontinu (x),
beberapa faktor (misal kelompok 3 level), perlakuan biner
(treat), dan satu faktor tambahan (region 4
level). Sinyal dan noise diatur agar interpretasi koefisien bisa diuji
lewat grafik.
library(tidyverse)
library(broom)
library(modelr)
library(scales)
library(kableExtra)
set.seed(123)4.2 Membuat Data Simulasi Dasar
npengamatan per kelompok agar sampling cukup.xsebagai kovariat kontinu (misal suhu, usia, pengeluaran iklan).kelompokberlevel"A","B","C"untuk contoh multi-kategori.treatsebagai dummy 0/1 (misal sebelum/sesudah kebijakan, kontrol/perlakuan).regionempat level untuk contoh banyak dummy (dan jebakan dummy variable trap).
4.3 Skenario 1: Intersep Berbeda (garis paralel)
Kita membuat true model yang intersepnya berbeda per
kelompok, tetapi kemiringannya sama untuk
x.
n <- 300
dat_vi <- tibble(
kelompok = factor(rep(c("A","B","C"), each = n/3)),
x = rnorm(n, mean = 10, sd = 3)
) %>%
mutate(
# Intersep dasar
beta0 = 5,
# Tambahan intersep per kelompok
delta = case_when(
kelompok == "A" ~ 0,
kelompok == "B" ~ 3,
kelompok == "C" ~ -2
),
beta1 = 1.2, # kemiringan sama untuk semua kelompok
eps = rnorm(n, 0, 2),
y = beta0 + delta + beta1 * x + eps
)
dat_vi %>% head() %>% knitr::kable() %>% kable_styling(full_width = FALSE)| kelompok | x | beta0 | delta | beta1 | eps | y |
|---|---|---|---|---|---|---|
| A | 8.318573 | 5 | 0 | 1.2 | -1.4304844 | 13.55180 |
| A | 9.309468 | 5 | 0 | 1.2 | -1.5053779 | 14.66598 |
| A | 14.676125 | 5 | 0 | 1.2 | -1.8770774 | 20.73427 |
| A | 10.211525 | 5 | 0 | 1.2 | -2.1050266 | 15.14880 |
| A | 10.387863 | 5 | 0 | 1.2 | -0.8743191 | 16.59112 |
| A | 15.145195 | 5 | 0 | 1.2 | 0.6623583 | 23.83659 |
4.4 Skenario 2: Kemiringan Berbeda (interaksi dummy–kontinu)
Sekarang kemiringan terhadap x bergantung pada
kelompok. Intersep sama, tetapi slope berbeda.
n <- 300
dat_vs <- tibble(
kelompok = factor(rep(c("A","B","C"), each = n/3)),
x = rnorm(n, mean = 10, sd = 3)
) %>%
mutate(
beta0 = 8,
# Slope berbeda per kelompok
beta1 = case_when(
kelompok == "A" ~ 0.8,
kelompok == "B" ~ 1.5,
kelompok == "C" ~ 0.3
),
eps = rnorm(n, 0, 2),
y = beta0 + beta1 * x + eps
)
dat_vs %>% head() %>% knitr::kable() %>% kable_styling(full_width = FALSE)| kelompok | x | beta0 | beta1 | eps | y |
|---|---|---|---|---|---|
| A | 13.222037 | 8 | 0.8 | -2.0282283 | 16.54940 |
| A | 9.917959 | 8 | 0.8 | -1.5826278 | 14.35174 |
| A | 9.900009 | 8 | 0.8 | 0.5991874 | 16.51919 |
| A | 5.451797 | 8 | 0.8 | 3.2781038 | 15.63954 |
| A | 12.371156 | 8 | 0.8 | 2.1692340 | 20.06616 |
| A | 9.367797 | 8 | 0.8 | -1.2491349 | 14.24510 |
4.5 Skenario 3: Interaksi Dummy–Dummy (dua kategori saling memengaruhi)
Di sini ada treat (0/1) dan kelompok
(A/B/C). Efek treat tidak sama di setiap
kelompok—inilah interaksi dummy–dummy.
n <- 600
dat_dd <- tibble(
kelompok = factor(rep(c("A","B","C"), each = n/3)),
treat = rbinom(n, 1, 0.5),
x = rnorm(n, 10, 3) # ikutkan 1 kovariat kontinu agar realistis
) %>%
mutate(
beta0 = 6,
beta1 = 1.0, # slope x (sama)
# Efek treat baseline (pada kelompok A)
tau_A = 2.0,
# Modifikasi efek treat pada B & C
tau_mod = case_when(
kelompok == "A" ~ 0.0,
kelompok == "B" ~ 1.5, # treat lebih kuat di B
kelompok == "C" ~ -1.0 # treat lebih lemah di C
),
eps = rnorm(n, 0, 2),
y = beta0 + beta1 * x + (tau_A + tau_mod) * treat + eps
)
dat_dd %>% head() %>% knitr::kable() %>% kable_styling(full_width = FALSE)| kelompok | treat | x | beta0 | beta1 | tau_A | tau_mod | eps | y |
|---|---|---|---|---|---|---|---|---|
| A | 1 | 7.537040 | 6 | 1 | 2 | 0 | -0.4002940 | 15.13675 |
| A | 1 | 9.078228 | 6 | 1 | 2 | 0 | 0.7756405 | 17.85387 |
| A | 0 | 7.293706 | 6 | 1 | 2 | 0 | 1.5878367 | 14.88154 |
| A | 1 | 11.881206 | 6 | 1 | 2 | 0 | -0.2810279 | 19.60018 |
| A | 1 | 13.361065 | 6 | 1 | 2 | 0 | 0.9116104 | 22.27268 |
| A | 0 | 16.381641 | 6 | 1 | 2 | 0 | -2.2911458 | 20.09049 |
4.6 Skenario 4: Banyak Dummy (misal region 4 level) + Dummy Variable Trap
Kita buat contoh dengan faktor region berlevel empat.
Kita akan lihat bagaimana R otomatis membuat k-1 dummy
dan kenapa jangan memasukkan semua dummy
plus intercept (itu menciptakan kolinearitas
sempurna).
n <- 400
dat_region <- tibble(
region = factor(sample(paste0("R",1:4), n, replace = TRUE)),
x = rnorm(n, 9, 2.5)
) %>%
mutate(
beta0 = 7,
# Intersep berbeda per region, baseline R1 = 0
delta = case_when(
region == "R1" ~ 0,
region == "R2" ~ 2.0,
region == "R3" ~ -1.5,
region == "R4" ~ 0.5
),
beta1 = 1.1,
eps = rnorm(n, 0, 1.8),
y = beta0 + delta + beta1 * x + eps
)
dat_region %>%
count(region) %>%
knitr::kable(col.names = c("Region","n")) %>%
kable_styling(full_width = FALSE)| Region | n |
|---|---|
| R1 | 106 |
| R2 | 93 |
| R3 | 98 |
| R4 | 103 |
4.7 Regresi: Intersep Berbeda (Garis Paralel)
Model umum: \[
y_i = \beta_0 + \beta_1 x_i + \gamma_1 D_{B,i} + \gamma_2 D_{C,i} +
\varepsilon_i,
\] dengan baseline A (jadi tidak ada \(D_A\)). Koefisien \(\gamma_1\) dan \(\gamma_2\) adalah pergeseran
intersep relatif terhadap A. Slope \(\beta_1\) sama untuk semua kelompok.
m_vi <- lm(y ~ x + kelompok, data = dat_vi)
summary(m_vi)##
## Call:
## lm(formula = y ~ x + kelompok, data = dat_vi)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.5530 -1.2101 0.0667 1.3058 5.2295
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.31350 0.46156 11.512 < 2e-16 ***
## x 1.16242 0.04059 28.637 < 2e-16 ***
## kelompokB 3.26183 0.28106 11.605 < 2e-16 ***
## kelompokC -2.00877 0.28005 -7.173 5.91e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.98 on 296 degrees of freedom
## Multiple R-squared: 0.7845, Adjusted R-squared: 0.7823
## F-statistic: 359.1 on 3 and 296 DF, p-value: < 2.2e-164.7.1 Interpretasi Koefisien
(Intercept): intersep untuk baseline (A).kelompokB: selisih intersep antaraBdanA.kelompokC: selisih intersep antaraCdanA.x: kemiringan yang sama untuk semua kelompok.
Bandingkan dengan nilai simulasi: kita set delta_B=+3
dan delta_C=-2—lihat apakah koefisien mendekati ini.
4.7.2 Visualisasi: Garis Paralel
dat_vi %>%
ggplot(aes(x, y, color = kelompok)) +
geom_point(alpha = 0.35) +
geom_smooth(method = "lm", se = FALSE) +
labs(title = "Varying Intercept: Garis Regresi Paralel",
subtitle = "Intersep berbeda per kelompok, slope sama",
x = "x (kovariat kontinu)", y = "y (respon)") +
theme_minimal()
4.7.3 Prediksi Per Kelompok
Gunakan modelr::add_predictions untuk menggambar garis
prediksi terpisah—tetap sejajar.
grid_vi <- dat_vi %>%
group_by(kelompok) %>%
summarize(x = seq(min(x), max(x), length.out = 50), .groups = "drop")
pred_vi <- grid_vi %>%
add_predictions(m_vi)
pred_vi %>% head() %>% knitr::kable() %>% kable_styling(full_width = FALSE)| kelompok | x | pred |
|---|---|---|
| A | 3.072493 | 8.885044 |
| A | 3.347789 | 9.205055 |
| A | 3.623085 | 9.525066 |
| A | 3.898382 | 9.845076 |
| A | 4.173678 | 10.165087 |
| A | 4.448974 | 10.485098 |
pred_vi %>%
ggplot(aes(x, pred, color = kelompok)) +
geom_line(size = 1.2) +
labs(title = "Garis Prediksi per Kelompok (Intersep Berbeda)",
x = "x", y = "Prediksi y") +
theme_minimal()
4.8 Regresi: Kemiringan Berbeda (Interaksi Dummy–Kontinu)
Model umum: \[ y_i = \beta_0 + \beta_1 x_i + \gamma_1 D_{B,i} + \gamma_2 D_{C,i} + \delta_1 (x_i \cdot D_{B,i}) + \delta_2 (x_i \cdot D_{C,i}) + \varepsilon_i. \]
x:kelompokBmengukur selisih slope B vs A.x:kelompokCmengukur selisih slope C vs A.- Intersep untuk B/C juga bisa berbeda jika kita sertakan
kelompoktanpa interaksi—dan memang kita sertakan (sehingga model lengkap).
m_vs <- lm(y ~ x * kelompok, data = dat_vs)
summary(m_vs)##
## Call:
## lm(formula = y ~ x * kelompok, data = dat_vs)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.8955 -1.4285 -0.0354 1.0557 6.3866
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.09095 0.66795 12.113 < 2e-16 ***
## x 0.78646 0.06658 11.812 < 2e-16 ***
## kelompokB -0.88729 0.98845 -0.898 0.370
## kelompokC -0.03020 0.96732 -0.031 0.975
## x:kelompokB 0.81397 0.09501 8.567 6.09e-16 ***
## x:kelompokC -0.49610 0.09311 -5.328 1.98e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.044 on 294 degrees of freedom
## Multiple R-squared: 0.9015, Adjusted R-squared: 0.8998
## F-statistic: 537.9 on 5 and 294 DF, p-value: < 2.2e-164.8.1 Interpretasi
(Intercept): intersep baseline (A) saatx=0. Jikaxjauh dari nol, pertimbangkan mean-centering padaxagar intersep lebih bermakna.x: slope untuk baseline (A).kelompokB,kelompokC: perbedaan intersep padax=0antara B/C dan A.x:kelompokB,x:kelompokC: perbedaan slope B/C vs A.
4.8.2 Visualisasi: Garis Berpotongan
dat_vs %>%
ggplot(aes(x, y, color = kelompok)) +
geom_point(alpha = 0.35) +
geom_smooth(method = "lm", se = FALSE) +
labs(title = "Varying Slope: Garis dengan Kemiringan Berbeda",
subtitle = "Interaksi x * kelompok mengizinkan slope berbeda",
x = "x", y = "y") +
theme_minimal()
4.8.3 Centring: Membuat Intersep Lebih Bermakna
dat_vs_c <- dat_vs %>% mutate(x_c = x - mean(x))
m_vs_c <- lm(y ~ x_c * kelompok, data = dat_vs_c)
broom::tidy(m_vs_c) %>% knitr::kable(digits = 3) %>% kable_styling(full_width = FALSE)| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | 15.995 | 0.207 | 77.256 | 0 |
| x_c | 0.786 | 0.067 | 11.812 | 0 |
| kelompokB | 7.293 | 0.291 | 25.021 | 0 |
| kelompokC | -5.016 | 0.291 | -17.219 | 0 |
| x_c:kelompokB | 0.814 | 0.095 | 8.567 | 0 |
| x_c:kelompokC | -0.496 | 0.093 | -5.328 | 0 |
Dengan x_c, (Intercept) sekarang
merepresentasikan rata-rata y pada x
rata-rata untuk baseline. Koefisien interaksi tidak berubah
interpretasinya (tetap selisih slope).
4.9 Interaksi Dummy–Dummy (Efek Perlakuan Berbeda per Kelompok)
Model: \[ y = \beta_0 + \beta_1 x + \tau \cdot \text{treat} + \gamma_B D_B + \gamma_C D_C + \eta_B (\text{treat}\cdot D_B) + \eta_C (\text{treat}\cdot D_C) + \varepsilon. \]
treatadalah dummy 0/1, baseline efeknya berlaku pada kelompok A.treat:kelompokBadalah tambahan efek treat khusus di B (di atas baseline A).treat:kelompokCtambahan efek treat khusus di C.
m_dd <- lm(y ~ x + treat * kelompok, data = dat_dd)
summary(m_dd)##
## Call:
## lm(formula = y ~ x + treat * kelompok, data = dat_dd)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.5744 -1.2409 0.0103 1.3311 5.8518
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.54592 0.33554 16.528 < 2e-16 ***
## x 1.00886 0.02634 38.301 < 2e-16 ***
## treat 2.25936 0.27511 8.213 1.36e-15 ***
## kelompokB 0.55336 0.27145 2.039 0.04194 *
## kelompokC 0.36787 0.27502 1.338 0.18153
## treat:kelompokB 0.92275 0.38884 2.373 0.01796 *
## treat:kelompokC -1.28031 0.38872 -3.294 0.00105 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.943 on 593 degrees of freedom
## Multiple R-squared: 0.7419, Adjusted R-squared: 0.7392
## F-statistic: 284 on 6 and 593 DF, p-value: < 2.2e-164.9.1 Menghitung Efek Treat per Kelompok
library(dplyr)
library(broom)
library(knitr)
library(kableExtra)
coefs <- broom::tidy(m_dd) %>%
dplyr::select(term, estimate)
tau_A <- coefs %>% filter(term == "treat") %>% pull(estimate)
tau_B <- tau_A + (coefs %>% filter(term == "treat:kelompokB") %>% pull(estimate))
tau_C <- tau_A + (coefs %>% filter(term == "treat:kelompokC") %>% pull(estimate))
tibble(
Kelompok = c("A", "B", "C"),
Efek_Treat = c(tau_A, tau_B, tau_C)
) %>%
knitr::kable(digits = 3) %>%
kableExtra::kable_styling(full_width = FALSE)| Kelompok | Efek_Treat |
|---|---|
| A | 2.259 |
| B | 3.182 |
| C | 0.979 |
4.9.2 Visualisasi: Difference-in-Differences ala Sederhana
dat_dd %>%
group_by(kelompok, treat) %>%
summarize(mean_y = mean(y), .groups = "drop") %>%
mutate(treat = factor(treat, levels = c(0,1), labels = c("Kontrol","Treat"))) %>%
ggplot(aes(treat, mean_y, group = kelompok, linetype = kelompok)) +
geom_point(size = 3) +
geom_line() +
labs(title = "Efek Treat Berbeda per Kelompok",
subtitle = "Garis menunjukkan pergeseran mean saat treat=1 vs 0",
x = "Status", y = "Rata-rata y") +
theme_minimal()
4.10 Banyak Dummy: Kategori dengan Banyak Level
Jika region punya 4 level, R akan membangun otomatis 3
dummy (k-1) dengan baseline R1. Lihat
model.matrix() untuk memahami konstruksinya.
MM <- model.matrix(~ region, data = dat_region)
head(MM) %>% knitr::kable() %>% kable_styling(full_width = FALSE)| (Intercept) | regionR2 | regionR3 | regionR4 |
|---|---|---|---|
| 1 | 1 | 0 | 0 |
| 1 | 0 | 0 | 1 |
| 1 | 0 | 0 | 0 |
| 1 | 0 | 0 | 0 |
| 1 | 0 | 0 | 0 |
| 1 | 0 | 0 | 0 |
4.10.1 Pasang Model dan Interpretasi
m_region <- lm(y ~ x + region, data = dat_region)
summary(m_region)##
## Call:
## lm(formula = y ~ x + region, data = dat_region)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.6656 -1.0873 0.1228 1.2255 4.4737
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.65097 0.35811 21.365 < 2e-16 ***
## x 1.06411 0.03446 30.879 < 2e-16 ***
## regionR2 1.46181 0.25310 5.776 1.55e-08 ***
## regionR3 -1.78943 0.24857 -7.199 3.08e-12 ***
## regionR4 0.16767 0.24540 0.683 0.495
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.773 on 395 degrees of freedom
## Multiple R-squared: 0.7264, Adjusted R-squared: 0.7236
## F-statistic: 262.2 on 4 and 395 DF, p-value: < 2.2e-16Interpretasi: - (Intercept): intersep untuk
R1 (baseline). -
regionR2,regionR3,regionR4:
perbedaan intersep relatif R1. - x: slope yang
sama untuk semua region (karena belum ada interaksi dengan
region).
4.10.2 Visualisasi
dat_region %>%
ggplot(aes(x, y, color = region)) +
geom_point(alpha = 0.35) +
geom_smooth(method = "lm", se = FALSE) +
labs(title = "Banyak Dummy (Region): Intersep Berbeda",
x = "x", y = "y") +
theme_minimal()
4.11 Interaksi Lebih Lanjut
4.11.1 Interaksi Faktor × Faktor (k × m level)
Jika dua faktor masing-masing memiliki \(k\) dan \(m\) level, interaksi akan menghasilkan \((k-1)\times(m-1)\) koefisien tambahan (di atas main effects), karena semua koefisien didefinisikan relatif terhadap baseline ganda.
Contoh: kelompok (A/B/C) × region (R1–R4)
memerlukan banyak koefisien. Gunakan ini hanya bila
perlu, atau ketika ada hipotesis teoritis yang meyakinkan
mengenai heterogenitas efek.
dat_big <- expand_grid(
kelompok = factor(c("A","B","C")),
region = factor(paste0("R",1:4))
) %>%
mutate(
# simulasi intersep unik per kombinasi (sekadar contoh)
mu = c(5,6,7,8, 4,6,6.5,7.5, 5.5,5,6.7,7.2)
)
# Gabungkan dengan x acak dan noise
set.seed(42)
dat_big <- dat_big %>%
group_by(kelompok, region) %>%
do({
tibble(
x = rnorm(80, 10, 2.5),
y = .$mu + 1.1 * x + rnorm(80,0,2)
)
}) %>% ungroup()
m_big <- lm(y ~ x + kelompok * region, data = dat_big)
broom::glance(m_big) %>% knitr::kable() %>% kable_styling(full_width = FALSE)| r.squared | adj.r.squared | sigma | statistic | p.value | df | logLik | AIC | BIC | deviance | df.residual | nobs |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 0.7223359 | 0.7188174 | 1.965931 | 205.2996 | 0 | 12 | -2004.564 | 4037.128 | 4105.265 | 3660.047 | 947 | 960 |
4.11.2 Plot Interaksi
dat_big %>%
group_by(kelompok, region) %>%
summarize(mean_y = mean(y), .groups = "drop") %>%
ggplot(aes(region, mean_y, fill = kelompok)) +
geom_col(position = position_dodge(width = 0.7)) +
labs(title = "Rata-rata y per Kombinasi Kelompok × Region",
x = "Region", y = "Mean y") +
theme_minimal()
4.12 Mengapa Menggunakan Dummy?
- Merepresentasikan kategori: Tanpa dummy, faktor kategorikal tak akan masuk ke regresi linier klasik.
- Kontras: Dummy memungkinkan kita mengestimasi perbedaan rata-rata antar kelompok (intersep berbeda).
- Heterogenitas efek: Interaksi dummy–kontinu dan dummy–dummy membuat kita bisa menangkap efek yang tidak homogen (misalnya kebijakan lebih efektif di satu provinsi daripada yang lain).
- Keterbacaan kebijakan: Koefisien dummy adalah selisih yang mudah dijelaskan ke non-teknis.
- Kontrol confounding: Menambahkan dummy untuk kategori yang relevan dapat mengurangi bias karena variabel konfonder terklasifikasi (misal tahun, musim, batch eksperimen).
4.13 Praktik Terbaik dan Jebakan
- Baseline yang jelas: Ingat, interpretasi dummy selalu relatif terhadap baseline. Pilih baseline yang paling umum atau paling bermakna substantif.
- Dummy Variable Trap: Jangan masukkan semua dummy dan intersep sekaligus. R default aman, tetapi hati-hati saat membuat dummy manual.
- Unbalanced data: Jika satu level jarang, estimasi interaksi bisa tak stabil. Pertimbangkan penggabungan level atau regularisasi.
- Centering: Untuk interaksi dummy–kontinu,
centering
xmemudahkan interpretasi intersep. - P-value vs ukuran efek: Jangan terpaku pada signifikansi. Plot selalu membantu. Ukuran efek (perbedaan rata-rata, slope) lebih komunikatif.
- Robust SE: Pada data heteroskedastik, gunakan
sandwich::vcovHC+lmtest::coeftestuntuk SE robust. - Kausalitas: Dummy treat tidak otomatis kausal. Butuh desain (randomisasi, DiD, IV, RDD). Di sini kita fokus deskriptif/regresif.
library(lmtest)
library(sandwich)
coeftest(m_vi, vcov = vcovHC(m_vi, type = "HC3"))##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.313504 0.475758 11.168 < 2.2e-16 ***
## x 1.162424 0.039543 29.396 < 2.2e-16 ***
## kelompokB 3.261833 0.293718 11.105 < 2.2e-16 ***
## kelompokC -2.008766 0.282171 -7.119 8.263e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 14.14 Ringkasan Interpretasi Koefisien
- Model tanpa interaksi:
\(\hat{y} = \hat{\beta}_0 + \hat{\beta}_1 x + \hat{\gamma}_B D_B + \hat{\gamma}_C D_C\)Intercept: mean untuk baseline saat \(x=0\).
x: slope umum.
D_B: selisih mean (intersep) B vs A.
D_C: selisih mean (intersep) C vs A.
- Model dengan interaksi dummy–kontinu:
\(\hat{y} = \hat{\beta}_0 + \hat{\beta}_1 x + \hat{\gamma}_B D_B + \hat{\gamma}_C D_C + \hat{\delta}_B xD_B + \hat{\delta}_C xD_C\)- Slope A = \(\hat{\beta}_1\)
- Slope B = \(\hat{\beta}_1 +
\hat{\delta}_B\)
- Slope C = \(\hat{\beta}_1 + \hat{\delta}_C\)
- Slope A = \(\hat{\beta}_1\)
- Model treat × kelompok:
Efek treat pada A =treat
Efek treat pada B =treat+treat:kelompokB
Efek treat pada C =treat+treat:kelompokC
4.15 Contoh Aplikasi Mini
4.15.1 Evaluasi Kebijakan (Sebelum–Sesudah vs Kelompok)
Misalkan treat=1 sesudah kebijakan. Interaksi
treat * kelompok memeriksa apakah dampak kebijakan berbeda
antar kelompok. Cocok untuk heterogeneous treatment effects
(HTE).
4.15.2 Bisnis & Pemasaran
kelompok = kanal penjualan; x = anggaran
iklan. Interaksi x * kanal memberitahu kanal mana yang
lebih responsif terhadap iklan (slope lebih curam).
4.16 Kesehatan Lingkungan
region sebagai provinsi; x = polutan;
kelompok = kebijakan emisi (ya/tidak). Interaksi menangkap
apakah hubungan polutan–outcome tergantung status kebijakan.
4.17 Diagnostik
Regresi linier standar mensyaratkan asumsi residual “sehat”. Kita cek cepat pada salah satu model.
par(mfrow = c(1,2))
plot(m_vi, which = 1) # residual vs fitted
plot(m_vi, which = 2) # QQ plot
par(mfrow = c(1,1))4.18 Epilog
Variabel dummy itu seperti saklar: kecil, tapi mengubah aliran arus
interpretasi. Dengan intersep berbeda, kita akui tiap kelompok punya
“titik awal” berbeda. Dengan slope berbeda, kita izinkan kelompok
merespons perubahan x secara unik. Dengan interaksi, kita
menangkap dialektika antar variabel—kapan efek ini tergantung pada
yang itu.
Jika Anda mengajar, plot–plot di atas adalah alat retoris yang sangat kuat: tidak sekadar berkata “koefisien beda”, tetapi menunjukkan garis yang benar‑benar beda. Dan bila suatu saat ada yang bilang “dummy itu cuma 0/1”, Anda sudah tahu balasannya: betul—namun 0/1 yang bisa menceritakan 1001 cerita.
5 Panel Data
Data panel adalah dua dimensi kebahagiaan statistik: kita mengamati banyak unit (individu, perusahaan, provinsi) berulang kali pada banyak periode waktu. Keuntungan utama:
- Mengendalikan heterogenitas tak-teramati yang konstan terhadap waktu (mis. bakat manajer, kultur daerah).
- Mengurangi kolinearitas dan meningkatkan degrees of freedom.
- Mengizinkan dinamika temporal (lag, tren) dan struktur korelasi yang realistis.
Dokumen ini menyajikan panduan komprehensif untuk tiga model klasik regresi panel:
- Pooled OLS (semua data “digabung” seakan-akan cross-section besar).
- Fixed Effects (FE) — efek tetap per unit dan/atau waktu.
- Random Effects (RE) — efek acak per unit dan/atau waktu.
Kita juga membahas perbedaan format data wide vs long, proses transformasi, pilihan model, uji-uji penting (LM Breusch–Pagan, F-test FE, Hausman), serta praktik terbaik (SE robust, klaster, diagnostik dependensi serial dan ketergantungan silang). Semua disertai simulasi data, kode R, dan visualisasi.
5.1 Paket yang Digunakan
library(tidyverse)
library(plm) # jagoan panel
library(lmtest) # uji diagnostik
library(sandwich) # robust vcov
library(broom) # ringkasan tabel
library(knitr)
library(kableExtra)
library(scales)
set.seed(123)5.2 Struktur Data Panel: Wide vs Long
Sebelum membicarakan model, pastikan format data jelas.
5.2.1 Format Long (Panel “Ramah-Model”)
Format long memiliki kolom pengenal unit (mis.
id), pengenal waktu (mis. year), lalu
kolom-kolom variabel:
id | year | y | x1 | x2 | z ...
1 | 2018 | . | . | . | .
1 | 2019 | . | . | . | .
...
2 | 2018 | . | . | . | .Sebagian besar fungsi panel (plm, fixest,
lfe) menyukai format ini.
5.2.2 Format Wide (Melebar ke Samping)
Format wide menaruh setiap periode sebagai kolom berbeda:
id | y_2018 | y_2019 | y_2020 | x1_2018 | x1_2019 | ...
1 | . | . | . | . | . | ...
2 | . | . | . | . | . | ...Format ini enak untuk melihat snapshot antar-periode sekaligus, tetapi perlu diubah ke long untuk dianalisis dengan mudah.
5.2.3 Transformasi Wide ↔︎ Long
# Contoh mini untuk transformasi
wide_df <- tibble(
id = 1:3,
y_2018 = c(10, 12, 9),
y_2019 = c(11, 15, 10),
x_2018 = c(3, 4, 3.5),
x_2019 = c(3.3, 4.2, 3.7)
)
# Wide -> Long
long_df <- wide_df %>%
pivot_longer(cols = -id, names_to = "var", values_to = "value") %>%
separate(var, into = c("name", "year"), sep = "_") %>%
pivot_wider(names_from = name, values_from = value) %>%
mutate(year = as.integer(year)) %>%
arrange(id, year)
long_df %>% kable() %>% kable_styling(full_width = FALSE)| id | year | y | x |
|---|---|---|---|
| 1 | 2018 | 10 | 3.0 |
| 1 | 2019 | 11 | 3.3 |
| 2 | 2018 | 12 | 4.0 |
| 2 | 2019 | 15 | 4.2 |
| 3 | 2018 | 9 | 3.5 |
| 3 | 2019 | 10 | 3.7 |
Ingat: long adalah format kerja utama untuk model panel.
5.3 Simulasi Data Panel
Kita bangun data yang “masuk akal” untuk mendemonstrasikan pooled, FE, dan RE.
- Unit: 60 perusahaan (id = 1..60).
- Periode: 6 tahun (t = 2018..2023).
- Variabel:
y(kinerja),x1(belanja R&D),x2(leverage),z(shock makro tahunan).
- Ada heterogenitas unit via
alpha_i(waktu-tetap) dan tren waktu vialambda_t.
library(dplyr)
library(tidyr)
library(kableExtra)
set.seed(123)
N <- 60 # unit
Tt <- 6 # time periods
years <- 2018:(2018 + Tt - 1)
panel <- expand.grid(id = 1:N, year = years) %>%
arrange(id, year)
# Komponen tak teramati: efek tetap unit dan efek tetap waktu
alpha_i <- rnorm(N, sd = 1.2) # heterogenitas per unit (tetap)
lambda_t <- seq(-0.5, 0.5, length.out = Tt) # tren waktu linear halus
# (Opsional) buat kamus tahun -> lambda
year_map <- tibble(
year = years,
lambda = lambda_t
)
# Kovariat observabel + gabungkan komponen
# Gunakan match / left_join agar setiap baris mendapatkan scalar 'alpha' dan 'lambda'
panel <- panel %>%
mutate(
x1 = rlnorm(n(), meanlog = log(5), sdlog = 0.3), # eksponensial-positif
x2 = rnorm(n(), mean = 0, sd = 1), # leverage ~ N(0,1)
z = rep(rnorm(Tt, 0, 0.5), each = N) # shock makro tahunan (sama untuk semua unit pada t)
) %>%
# alpha_i per unit (scalar per baris)
mutate(alpha = alpha_i[id]) %>%
# lambda_t per tahun (scalar per baris) — pilih salah satu:
# 1) via match:
mutate(lambda = lambda_t[match(year, years)])
# 2) atau pakai join (komentari baris match di atas, lalu pakai ini):
# left_join(year_map, by = "year")
# (setelah join, 'lambda' sudah ada sebagai kolom)
# True DGP (data generating process)
beta0 <- 1.5
beta1 <- 0.6 # x1
beta2 <- -0.4 # x2
beta3 <- 0.8 # z
# Error & outcome
panel <- panel %>%
mutate(
u = rnorm(n(), 0, 0.8),
y = beta0 + beta1*x1 + beta2*x2 + beta3*z + alpha + lambda + u
)
# Tampilkan cuplikan
panel %>%
head() %>%
knitr::kable() %>%
kable_styling(full_width = FALSE)| id | year | x1 | x2 | z | alpha | lambda | u | y |
|---|---|---|---|---|---|---|---|---|
| 1 | 2018 | 5.603155 | 1.7795029 | 0.4067002 | -0.6725708 | -0.5 | 0.6379044 | 3.940785 |
| 1 | 2019 | 4.300541 | 0.2864244 | 0.4067002 | -0.6725708 | -0.3 | -1.5665641 | 1.751980 |
| 1 | 2020 | 4.524358 | 0.1263159 | 0.4067002 | -0.6725708 | -0.1 | -1.5090601 | 2.207818 |
| 1 | 2021 | 3.683507 | 1.2722668 | 0.4067002 | -0.6725708 | 0.1 | -0.5230239 | 2.430963 |
| 1 | 2022 | 3.625168 | -0.7184662 | 0.4067002 | -0.6725708 | 0.3 | 0.3155159 | 4.230792 |
| 1 | 2023 | 5.476666 | -0.4503386 | 0.4067002 | -0.6725708 | 0.5 | -0.7308528 | 4.388071 |
Secara desain, alpha_i berkorelasi
positif dengan x1 (karena firm yang
“bagus” cenderung belanja R&D lebih besar). Untuk memperjelas
endogenitas halus ini, kita tambahkan korelasi ringan:
# Memperkenalkan korelasi ringan antara alpha_i dan x1:
panel <- panel %>% group_by(id) %>%
mutate(x1 = x1 + 0.2*alpha[1]) %>%
ungroup()
# Rehitung y agar konsisten
panel <- panel %>% mutate(y = beta0 + beta1*x1 + beta2*x2 + beta3*z + alpha + lambda + u)5.4 Visualisasi
ggplot(panel, aes(year, y, group = id)) +
geom_line(alpha = 0.2) +
stat_summary(aes(group = 1), fun = mean, geom = "line", size = 1.2) +
labs(title = "Rata-Rata dan Trajektori Unit (Tipis)",
subtitle = "Garis tebal = rata-rata agregat tahunan",
x = "Tahun", y = "y") +
theme_minimal()
5.5 Mempersiapkan Objek Panel
plm menyukai objek pdata.frame dengan
indeks (id, year).
pdata <- pdata.frame(panel, index = c("id","year"))5.6 Pooled OLS
Anggap data panel cuma data cross-section besar. Abaikan efek tetap/acak.
Model: \[ y_{it} = \beta_0 + \beta_1 x1_{it} + \beta_2 x2_{it} + \beta_3 z_t + \varepsilon_{it} \]
pooled <- plm(y ~ x1 + x2 + z, data = pdata, model = "pooling")
summary(pooled)## Pooling Model
##
## Call:
## plm(formula = y ~ x1 + x2 + z, data = pdata, model = "pooling")
##
## Balanced Panel: n = 60, T = 6, N = 360
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -3.312435 -0.958148 0.012376 0.901359 4.337458
##
## Coefficients:
## Estimate Std. Error t-value Pr(>|t|)
## (Intercept) 1.076380 0.244307 4.4058 1.396e-05 ***
## x1 0.702449 0.044590 15.7535 < 2.2e-16 ***
## x2 -0.364761 0.072186 -5.0531 6.963e-07 ***
## z 0.769942 0.164609 4.6774 4.132e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 1223.3
## Residual Sum of Squares: 666.95
## R-Squared: 0.45481
## Adj. R-Squared: 0.45021
## F-statistic: 98.9934 on 3 and 356 DF, p-value: < 2.22e-16Masalah potensial:
- Omitted Heterogeneity: jika \(x\) berkorelasi dengan karakteristik unit yang tak-teramati (konsisten waktu), OLS menjadi bias.
- Serial Correlation: error dalam unit yang sama lintas waktu sering berkorelasi.
- Heteroskedastisitas: varians tak konstan antar unit/waktu.
Solusi minimal: gunakan SE robust berklaster
(klaster pada id).
coeftest(pooled, vcov = vcovHC(pooled, type = "HC1", cluster = "group"))##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.076380 0.252806 4.2577 2.645e-05 ***
## x1 0.702449 0.043189 16.2647 < 2.2e-16 ***
## x2 -0.364761 0.062895 -5.7995 1.469e-08 ***
## z 0.769942 0.322614 2.3866 0.01753 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 15.7 Fixed Effects (FE)
5.7.1 Fixed Effect Model vs Endogeneity
Dokumen ini menjelaskan hubungan panel data dengan endogeneity, serta mengapa model efek tetap (fixed effects, FE) kerap dianggap solusi—walau hanya parsial—untuk isu endogeneity. Pembahasan mencakup intuisi, formulasi matematis, equivalence dua pendekatan FE (within vs LSDV), keterbatasan FE, dan demonstrasi simulasi kecil di R.
Ringkasan Eksekutif
- Endogeneity terjadi ketika variabel penjelas berkorelasi dengan galat, lazim karena omitted variable bias, simultaneity, atau measurement error.
- Panel data memungkinkan kita mengamati entitas yang sama berulang kali, sehingga karakteristik tak teramati yang tetap dalam waktu (time-invariant) dapat dikendalikan.
- Fixed effects (FE) secara konseptual mengakui adanya efek individu tetap \((\alpha_i)\), tetapi menghilangkannya secara matematis melalui transformasi within/demeaning, atau memasukkannya secara eksplisit melalui dummy (LSDV). Kedua pendekatan memberikan estimasi koefisien \(\beta\) yang ekuivalen.
- FE menangani endogeneity yang bersumber dari omitted variables yang time-invariant, namun tidak menyelesaikan endogeneity karena simultaneity atau measurement error. Untuk itu, gunakan IV, DiD, atau dynamic panel GMM.
- Simulasi kecil menunjukkan: ketika efek individu berkorelasi dengan \(X\), OLS pooled menjadi bias; FE mengoreksi bias tersebut.
Apa itu Endogeneity? Secara umum, endogeneity adalah kondisi ketika variabel penjelas \(X\) berkorelasi dengan error term \(\varepsilon\) dalam model regresi. Sumber utama:
- Omitted Variable Bias (OVB): ada faktor relevan yang tidak dimasukkan ke model.
- Simultaneity: sebab–akibat terjadi dua arah (\(X \leftrightarrow Y\)).
- Measurement Error: kesalahan pengukuran pada variabel kunci (sering di \(X\)).
Konsekuensinya: estimasi OLS menjadi bias dan tidak konsisten.
Di mana Panel Data Membantu? Panel data menggabungkan dimensi individu (atau entitas) dan waktu. Dengan mengamati individu yang sama berulang kali, kita bisa mengontrol faktor tak teramati yang tidak berubah dalam waktu (misal: bakat manajerial, budaya organisasi, lokasi tetap). Inilah sumber endogeneity yang paling sering mengganggu analisis cross-section murni.
3. Model Fixed Effects: Formulasi dan Intuisi Misalkan model dasar: \[ Y_{it} = \alpha_i + \beta X_{it} + \varepsilon_{it}, \quad i=1,\dots,N;\; t=1,\dots,T. \]
- \(\alpha_i\): efek tetap individu (time-invariant, tidak teramati).
- \(X_{it}\): variabel penjelas yang bervariasi lintas waktu dan/atau individu.
- \(\varepsilon_{it}\): galat idiosinkratik yang tidak berkorelasi dengan \(X_{it}\) setelah mengontrol \(\alpha_i\).
Asumsi kunci FE: setelah mengondisikan pada \(\alpha_i\), tidak ada korelasi antara \(X_{it}\) dan \(\varepsilon_{it}\). Korelasi antara \(\alpha_i\) dan \(X_{it}\) diperbolehkan.
Within (Demeaning) Transformation Ambil rata-rata per individu: \[ \bar{Y}_i = \alpha_i + \beta \bar{X}_i + \bar{\varepsilon}_i. \]
Kurangkan dari persamaan asli: \[ (Y_{it} - \bar{Y}_{i}) = \beta (X_{it} - \bar{X}_{i}) + (\varepsilon_{it} - \bar{\varepsilon}_{i}) \]
Karena \(\alpha_i\) bersifat konstan untuk setiap individu \(i\) dan tidak berubah sepanjang waktu, maka ketika setiap variabel dikurangi dengan rata-ratanya (\(\bar{Y}_{i}, \bar{X}_{i}\)), komponen \(\alpha_i\) akan tereliminasi secara sempurna dari persamaan. Dengan demikian, variasi yang tersisa hanya mencerminkan perubahan dalam individu (within variation), bukan perbedaan antarindividu.
Inilah maksud frasa “FE menghilangkan \(\alpha_i\)”.
LSDV (Least-Squares Dummy Variables) Alternatif ekuivalen: masukkan dummy untuk tiap individu (kecuali satu untuk menghindari dummy trap): \[ Y_{it} = \beta X_{it} + \sum_{i=1}^{N-1} \alpha_i D_i + \varepsilon_{it}. \] Di sini, \(\alpha_i\) dimasukkan secara eksplisit. Secara teori, estimasi \(\beta\) dari LSDV ekuivalen dengan hasil within transformation.
Inti klarifikasi: FE mengakui \(\alpha_i\) (efek individu tetap), namun menghapusnya secara aljabar (via demeaning) saat menaksir \(\beta\).
Kapan FE Menyelesaikan Endogeneity? FE efektif untuk OVB yang time-invariant. Contoh: - “Budaya perusahaan” yang tak terukur tetapi relatif tetap selama periode studi. - Perbedaan lokasi geografis tetap. - Karakteristik bawaan individu (mis. kemampuan kognitif dasar).
Dengan FE, sumber bias ini tereliminasi dari persamaan sehingga \(\beta\) tidak lagi tercemar oleh korelasi \(\alpha_i\)–\(X_{it}\).
Kapan FE Tidak Cukup? FE tidak menyelesaikan: - Simultaneity (reverse causality): jika \(Y\) juga memengaruhi \(X\) pada waktu yang sama, maka \(X_{it}\) tetap berkorelasi dengan galat (meski \(\alpha_i\) sudah “dibersihkan”). - Time-varying omitted variables: ada faktor tak teramati berubah sepanjang waktu yang memengaruhi \(X\) dan \(Y\) sekaligus (mis. kebijakan baru). - Measurement error di \(X\): kesalahan ukur tetap menimbulkan bias.
Solusi pelengkap: Instrumental Variables (IV), Difference-in-Differences (DiD), dan Dynamic Panel GMM (mis. Arellano–Bond).
Demonstrasi Simulasi (R) Simulasi ini menunjukkan:
1) ketika \(\alpha_i\) berkorelasi dengan \(X_{it}\), OLS
pooled menjadi bias;
2) FE memulihkan estimasi \(\beta\) mendekati nilai sebenarnya.
set.seed(123)
library(dplyr)
library(tidyr)
library(ggplot2)
library(fixest) # untuk feols (FE) dan ols
library(broom)
N <- 200 # jumlah individu
Tt <- 6 # jumlah periode
beta_true <- 0.8
# Efek individu yang tak teramati (alpha_i), mis. "bakat" perusahaan
alpha_i <- rnorm(N, mean = 0, sd = 1)
# Variabel pengganggu yang mendorong korelasi antara X dan alpha (time-invariant)
z_i <- scale(alpha_i + rnorm(N, sd = 0.3))[,1]
# Bangun panel
dat <- expand.grid(i = 1:N, t = 1:Tt) %>% arrange(i, t)
# X dipengaruhi oleh z_i (yang berkorelasi dengan alpha_i), plus shock waktu
dat <- dat %>%
mutate(alpha = alpha_i[i],
z = z_i[i],
nu = rnorm(n(), sd = 0.6),
X = 0.7*z + 0.4*runif(n(), -1, 1) + nu)
# Hasil Y
dat <- dat %>%
mutate(epsilon = rnorm(n(), sd = 1),
Y = alpha + beta_true*X + epsilon)
head(dat)## i t alpha z nu X epsilon Y
## 1 1 1 -0.5604756 0.09690107 -0.04413361 -0.06727743 2.2819670 1.6676694
## 2 1 2 -0.5604756 0.09690107 -0.70119085 -0.59945388 -0.4636830 -1.5037218
## 3 1 3 -0.5604756 0.09690107 -0.38084896 -0.11761574 -0.3263536 -0.9809218
## 4 1 4 -0.5604756 0.09690107 -0.01730493 -0.08280823 0.8824932 0.2557710
## 5 1 5 -0.5604756 0.09690107 0.40241758 0.33018052 1.2812861 0.9849549
## 6 1 6 -0.5604756 0.09690107 -0.99032793 -0.98350274 -0.6586819 -2.0059597Estimasi OLS Pooled vs FE
# OLS pooled (abaikan struktur panel)
m_ols <- feols(Y ~ X, data = dat, notes = FALSE)
# FE: kontrol efek tetap individu i
m_fe <- feols(Y ~ X | i, data = dat, notes = FALSE)
# Ringkasan koefisien
tbl <- bind_rows(tidy(m_ols) %>% mutate(model = "Pooled OLS"),
tidy(m_fe) %>% mutate(model = "Fixed Effects (i)")) %>%
filter(term == "X") %>%
mutate(beta_true = beta_true)
tbl## # A tibble: 2 × 7
## term estimate std.error statistic p.value model beta_true
## <chr> <dbl> <dbl> <dbl> <dbl> <chr> <dbl>
## 1 X 1.48 0.0374 39.7 1.39e-220 Pooled OLS 0.8
## 2 X 0.719 0.0488 14.7 1.62e- 44 Fixed Effects (i) 0.8Visualisasi Bias
ggplot(tbl, aes(x = model, y = estimate)) +
geom_point(size = 3) +
geom_errorbar(aes(ymin = estimate - 1.96*std.error,
ymax = estimate + 1.96*std.error),
width = 0.1) +
geom_hline(yintercept = beta_true, linetype = "dashed") +
labs(title = "Perbandingan Estimasi Koefisien X",
subtitle = "Garis putus-putus = nilai sebenarnya (beta_true = 0.8)",
x = NULL, y = "Estimasi koefisien") +
theme_minimal(base_size = 12)
Dalam simulasi ini, Pooled OLS cenderung menyimpang dari \(\beta\) yang benar karena \(X\) berkorelasi dengan \(\alpha_i\). FE mengurangi bias tersebut dengan “menghapus” \(\alpha_i\) melalui transformasi within.
Within vs LSDV: Bukti Ekuivalensi Secara teori,
estimasi \(\hat\beta\) dari within estimator dan LSDV
identik. Perbedaan ada pada implementasi komputasional (efisiensi ketika
\(N\) besar). Paket modern (mis. fixest, lfe,
plm) mengestimasi FE tanpa benar-benar membuat ribuan
dummy, melainkan menggunakan algoritme matriks yang menerapkan
transformasi within secara efisien.
Keterbatasan Praktis dan Saran
- Variasi dalam waktu diperlukan: jika \(X\) tidak berubah dalam waktu untuk tiap \(i\), FE tidak dapat mengidentifikasi koefisiennya.
- Efek waktu bersama (time fixed effects): sering
disarankan untuk mengontrol guncangan makro (mis. kebijakan nasional)
yang memengaruhi semua \(i\) pada periode \(t\). Contoh sintaks:
Y ~ X | i + t. - Standard error: pertimbangkan klasterisasi SE pada tingkat individu (atau two-way clustering) ketika relevan.
- Endogeneity residual: bila simultanitas atau
measurement error masih mungkin, pertimbangkan IV-FE
(
feols(..., iv=...)) atau dynamic panel GMM.
Kalimat Kunci
- “FE keep \(\alpha_i\) conceptually, but wipe it out
algebraically.”
- “FE menyelesaikan OVB time-invariant, bukan
semua bentuk endogeneity.”
- “Kombinasikan FE dengan IV/DiD/GMM untuk sumber endogeneity yang berubah sepanjang waktu atau simultanitas.”
Lampiran: Contoh Sintaks Alternatif Beberapa
alternatif paket/sintaks yang umum: -
plm::plm(Y ~ X, data = dat, index = c("i","t"), model = "within")
- lfe::felm(Y ~ X | i) -
fixest::feols(Y ~ X | i + t) (tambahkan efek waktu
bersama)
Catatan: untuk IV-FE di
fixest, gunakan formula instrumen, contoh:
feols(Y ~ 1 | i, iv = ~ X ~ Z, data = dat)di manaZadalah instrumen yang valid.
5.8 Konsep Fixed Effect Model
5.8.1 FE Satu Arah
Model Fixed Effect (FE) digunakan dalam data panel untuk mengontrol heterogenitas tak-teramati (unobserved heterogeneity) yang bersifat konstan dalam waktu tetapi berbeda antarunit (misalnya antar individu, perusahaan, atau daerah).
Model umum one-way fixed effect (unit-specific) adalah:
\[ (Y_{it} - \bar{Y}_{i}) = \beta (X_{it} - \bar{X}_{i}) + (\varepsilon_{it} - \bar{\varepsilon}_{i}) \]
dengan:
- \(y_{it}\): variabel dependen untuk unit \(i\) pada waktu \(t\).
- \(x1_{it}\): variabel independen yang berubah antar unit dan antar waktu (misalnya, tingkat investasi, output produksi, atau konsumsi energi).
- \(x2_{it}\): variabel independen lain yang juga bervariasi dalam waktu dan antarunit (misalnya upah rata-rata, pendapatan per kapita, atau indeks biaya hidup).
- \(z_t\): variabel yang berubah hanya terhadap waktu (misalnya inflasi nasional, kebijakan makro, atau efek waktu bersama).
- \(\alpha_i\): efek individu tetap yang konstan sepanjang waktu untuk setiap unit, tetapi berbeda antarunit.
- \(\varepsilon_{it}\): komponen error idiosinkratik yang bersifat acak.
Tujuan utama model FE adalah menghilangkan korelasi antara \(\alpha_i\) dan variabel penjelas \(x1_{it}, x2_{it}\), sehingga estimasi koefisien \(\beta\) menjadi tidak bias.
5.8.2 Estimasi dengan Least Squares Dummy Variable (LSDV)
Pendekatan LSDV (Least Squares Dummy Variable) memasukkan efek tetap \(\alpha_i\) secara eksplisit ke dalam model dengan menambahkan dummy variable untuk setiap unit individu.
Modelnya menjadi:
\[ y_{it} = \beta_1 x1_{it} + \beta_2 x2_{it} + \beta_3 z_t + \sum_{i=1}^{N-1} \alpha_i D_i + \varepsilon_{it} \]
dengan:
- \(D_i\): dummy untuk unit ke-\(i\) (1 jika observasi berasal dari unit \(i\), 0 jika tidak).
- \(N-1\): karena satu dummy harus dihapus agar tidak terjadi dummy variable trap (multikolinearitas sempurna).
Secara praktis, model ini diestimasi menggunakan OLS, namun dengan sejumlah besar dummy.
Estimator OLS dapat ditulis sebagai:
\[ \hat{\beta}_{LSDV} = (X'D_M X)^{-1} X'D_M y \]
di mana \(D_M\) adalah matriks yang memasukkan dummy untuk setiap unit individu.
Kelebihan LSDV: - Interpretasi mudah (efek tetap tampak eksplisit dalam hasil regresi). - Cocok untuk data panel dengan \(N\) kecil dan \(T\) besar.
Kelemahan LSDV: - Tidak efisien bila \(N\) besar (karena menambah banyak dummy). - Tidak cocok untuk data panel besar (misalnya ribuan individu).
Contoh sintaks di R menggunakan fixest:
library(fixest)
model_lsdv <- feols(y ~ x1 + x2 + z + factor(id), data = data_panel)
summary(model_lsdv)5.8.3 Estimasi dengan Within (Demeaning / Differencing Mean) Transformation
Pendekatan within transformation atau demeaning menghilangkan efek individu \(\alpha_i\) secara matematis tanpa menambahkan dummy.
Langkah-langkahnya:
(a) Ambil rata-rata per individu \[ \bar{y_i} = \frac{1}{T_i}\sum_{t=1}^{T_i} y_{it}, \quad \bar{x1_i} = \frac{1}{T_i}\sum_{t=1}^{T_i} x1_{it}, \quad \bar{x2_i} = \frac{1}{T_i}\sum_{t=1}^{T_i} x2_{it}, \quad \bar{z_i} = \frac{1}{T_i}\sum_{t=1}^{T_i} z_t \]
(b) Kurangkan setiap nilai aktual dengan rata-ratanya \[ (y_{it} - \bar{y_i}) = \beta_1 (x1_{it} - \bar{x1_i}) + \beta_2 (x2_{it} - \bar{x2_i}) + \beta_3 (z_t - \bar{z_i}) + (\varepsilon_{it} - \bar{\varepsilon_i}) \]
Karena \(\alpha_i\) konstan untuk setiap \(i\): \[ \alpha_i - \bar{\alpha_i} = 0 \]
maka efek individu terhapus dari persamaan.
Regresi dilakukan terhadap variabel yang sudah “dimean-kan”
(demeaned).
Estimator FE (within estimator) diberikan oleh:
\[ \hat{\beta}_{FE} = (X^{*\top} X^*)^{-1} X^{*\top} y^* \]
dengan:
- \(X^* = X_{it} - \bar{X_i}\)
- \(y^* = y_{it} - \bar{y_i}\)
Interpretasi:
Koefisien \(\hat{\beta}_{FE}\)
menggambarkan efek rata-rata dari perubahan \(x1_{it}\) atau \(x2_{it}\) dalam individu yang sama
dari waktu ke waktu, bukan antarindividu.
Perbandingan LSDV dan Within Transformation
| Aspek | LSDV | Within (Demeaning) |
|---|---|---|
| Efek individu \(\alpha_i\) | Dimodelkan eksplisit (dummy) | Dihilangkan secara matematis |
| Implementasi | Tambah dummy untuk setiap individu | Kurangkan rata-rata per individu |
| Efisiensi | Buruk jika N besar | Efisien untuk N besar |
| Interpretasi | Intersep spesifik terlihat | Intersep tidak terlihat langsung |
| Hasil \(\hat{\beta}\) | Identik | Identik |
| Variabel time-invariant | Bisa diestimasi | Tidak bisa diestimasi |
Contoh Implementasi di R
Contoh 1
Berikut contoh data simulasi sederhana dengan dua variabel independen \(x1, x2\), satu efek waktu \(z_t\), dan efek individu \(\alpha_i\).
set.seed(123)
library(dplyr)
library(plm)
# Simulasi data panel
N <- 100 # jumlah unit
Tt <- 5 # jumlah waktu
alpha_i <- rnorm(N, 0, 2) # efek tetap individu
beta <- c(0.5, -0.3, 0.2)
dat <- expand.grid(i = 1:N, t = 1:Tt) %>%
mutate(x1 = runif(N*Tt, 10, 20),
x2 = rnorm(N*Tt, 5, 1.5),
z = rnorm(Tt)[t],
u = rnorm(N*Tt, 0, 1),
alpha = alpha_i[i],
y = beta[1]*x1 + beta[2]*x2 + beta[3]*z + alpha + u)
# Konversi ke data panel
pdata <- pdata.frame(dat, index = c("i", "t"))
# Estimasi dengan LSDV
lsdv_model <- lm(y ~ x1 + x2 + z + factor(i), data = dat)
# Estimasi dengan within transformation (plm)
fe_model <- plm(y ~ x1 + x2 + z, data = pdata, model = "within")
summary(lsdv_model)$coefficients[1:5,]## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.5854085 0.52978028 -2.9925775 2.939124e-03
## x1 0.4984952 0.01652002 30.1752188 8.389282e-105
## x2 -0.3024865 0.03204598 -9.4391402 3.245824e-19
## z 0.1939942 0.03737305 5.1907505 3.350568e-07
## factor(i)2 0.4646629 0.61199709 0.7592568 4.481496e-01summary(fe_model)## Oneway (individual) effect Within Model
##
## Call:
## plm(formula = y ~ x1 + x2 + z, data = pdata, model = "within")
##
## Balanced Panel: n = 100, T = 5, N = 500
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -2.6693563 -0.5574777 0.0071086 0.5214703 2.6382967
##
## Coefficients:
## Estimate Std. Error t-value Pr(>|t|)
## x1 0.498495 0.016520 30.1752 < 2.2e-16 ***
## x2 -0.302486 0.032046 -9.4391 < 2.2e-16 ***
## z 0.193994 0.037373 5.1908 3.351e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 1298.9
## Residual Sum of Squares: 369.63
## R-Squared: 0.71542
## Adj. R-Squared: 0.64231
## F-statistic: 332.686 on 3 and 397 DF, p-value: < 2.22e-16Kedua metode menghasilkan estimasi \(\hat{\beta}\) yang identik, meskipun pendekatannya berbeda.
Interpretasi Hasil Jika \(\hat{\beta}_1 = 0.5\), artinya setiap kenaikan satu unit pada \(x1_{it}\) dalam individu yang sama (misalnya peningkatan investasi perusahaan dari tahun ke tahun) dikaitkan dengan kenaikan 0.5 unit pada \(y_{it}\), setelah mengendalikan efek tetap individu dan variabel lain.
Kesimpulan
- LSDV dan within estimator identik secara matematis
dalam hasil estimasi koefisien \(\beta\).
- LSDV cocok untuk \(N\) kecil, sementara within
estimator efisien untuk \(N\)
besar.
- Keduanya mengatasi endogeneity akibat omitted variable yang
time-invariant, namun tidak menyelesaikan
simultaneity atau measurement error.
- Dalam praktik empiris, penggunaan plm() atau fixest::feols() adalah cara paling efisien untuk mengestimasi model FE modern.
Contoh 2
fe <- plm(y ~ x1 + x2 + z, data = pdata, model = "within", effect = "individual")
summary(fe)## Oneway (individual) effect Within Model
##
## Call:
## plm(formula = y ~ x1 + x2 + z, data = pdata, effect = "individual",
## model = "within")
##
## Balanced Panel: n = 100, T = 5, N = 500
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -2.6693563 -0.5574777 0.0071086 0.5214703 2.6382967
##
## Coefficients:
## Estimate Std. Error t-value Pr(>|t|)
## x1 0.498495 0.016520 30.1752 < 2.2e-16 ***
## x2 -0.302486 0.032046 -9.4391 < 2.2e-16 ***
## z 0.193994 0.037373 5.1908 3.351e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 1298.9
## Residual Sum of Squares: 369.63
## R-Squared: 0.71542
## Adj. R-Squared: 0.64231
## F-statistic: 332.686 on 3 and 397 DF, p-value: < 2.22e-16coeftest(fe, vcov = vcovHC(fe, type = "HC1", cluster = "group"))##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## x1 0.498495 0.015928 31.2968 < 2.2e-16 ***
## x2 -0.302486 0.028804 -10.5016 < 2.2e-16 ***
## z 0.193994 0.035221 5.5079 6.539e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 15.8.4 FE Dua Arah (Unit + Waktu)
Dokumen ini merangkum spesifikasi Fixed Effects (FE) dua arah—unit dan waktu—beserta langkah estimasi dan implementasi R yang ringkas. Notasi galat menggunakan \(\varepsilon\) (bukan \(u\)).
Intuisi: selain heterogenitas tak-teramati antar-unit (\(\alpha_i\)), sering ada shock agregat setiap periode (\(\lambda_t\)). FE dua arah menyerap keduanya sehingga identifikasi koefisien slope \(\beta\) lebih bersih.
5.8.4.1 Spesifikasi Model
Model panel FE dua arah (unit + waktu): \[ Y_{it} \;=\; \alpha_i \;+\; \lambda_t \;+\; X_{it}^{\top}\beta \;+\; \varepsilon_{it}, \qquad i=1,\dots,N,\; t=1,\dots,T. \]
Komponen:
- \(Y_{it}\) : variabel dependen,
- \(X_{it}\) : vektor kovariat (ukuran \(K\times 1\)),
- \(\beta\) : parameter slope (ukuran \(K\times 1\)),
- \(\alpha_i\) : efek tetap unit (konstan sepanjang waktu),
- \(\lambda_t\) : efek tetap waktu (sama untuk semua unit pada periode \(t\)),
- \(\varepsilon_{it}\) : galat idiosinkratik.
Normalisasi identifikasi (pilih salah satu) : \[ \sum_{i=1}^N \alpha_i = 0, \quad \sum_{t=1}^T \lambda_t = 0 \quad \text{atau} \quad \alpha_1 = 0, \; \lambda_1 = 0 . \]
5.8.4.2 Estimasi
Tiga pendekatan umum yang saling ekivalen untuk memperoleh \(\hat{\beta}\).
5.8.4.2.1 Within (Double-Demeaning)
Definisikan rata-rata menurut unit, waktu, dan keseluruhan: \[ \bar{Y}_{i\cdot}=\frac{1}{T}\sum_t Y_{it},\quad \bar{Y}_{\cdot t}=\frac{1}{N}\sum_i Y_{it},\quad \bar{Y}=\frac{1}{NT}\sum_{i,t} Y_{it}. \]
Transformasi double-demeaning: \[ \tilde{Y}_{it} = Y_{it} - \bar{Y}_{i\cdot} - \bar{Y}_{\cdot t} + \bar{Y}, \qquad \tilde{X}_{it} = X_{it} - \bar{X}_{i\cdot} - \bar{X}_{\cdot t} + \bar{X}. \]
Setelah transformasi, efek tetap lenyap: \[ \tilde{Y}_{it} = \tilde{X}_{it}^{\top} \beta + \tilde{\varepsilon}_{it}. \]
Estimator OLS pada data bertanda “tilde”: \[ \hat{\beta}_{\text{FE,2W}} = (\tilde{X}^{\top}\tilde{X})^{-1} \tilde{X}^{\top}\tilde{y}. \]
Catatan proyeksi (opsional) : menggunakan operator pusat \(M_1\) (unit) dan \(M_2\) (waktu), double-demeaning setara dengan proyeksi Kronecker \(M_2 \otimes M_1\).
5.8.4.2.2 LSDV (Least Squares Dummy Variables)
Regresi OLS dengan dummy untuk semua unit dan semua waktu (kecuali satu per kelompok untuk menghindari dummy trap) : \[ Y_{it} = \gamma + \sum_{i=2}^{N} \delta_i \, \mathbb{1}\{i\} + \sum_{t=2}^{T} \tau_t \, \mathbb{1}\{t\} + X_{it}^{\top}\beta + \varepsilon_{it}. \]
Koefisien \(\beta\) identik dengan within; \(\delta_i\) dan \(\tau_t\) merepresentasikan \(\alpha_i\) dan \(\lambda_t\) relatif kategori acuan.
5.8.4.2.3 Frisch–Waugh–Lovell (FWL)
Regress-out \(\alpha_i\) dan \(\lambda_t\) dari \(Y\) dan setiap kolom \(X\), lalu OLS pada residualnya. Hasil \(\hat{\beta}\) sama dengan (A) dan (B).
5.8.4.3 Varian Galat & Standar Error
Panel jarang bersifat i.i.d. Gunakan SE yang tepat: - Cluster by unit: robust terhadap heteroskedastisitas dan autokorelasi dalam unit. - Driscoll–Kraay: robust terhadap dependensi lintas-seksi (cocok saat \(T\) sedang-besar). - Two-way cluster: unit dan waktu sekaligus bila tersedia.
Residual terestimasi: \[ \hat{\varepsilon}_{it} = Y_{it} - \hat{\alpha}_i - \hat{\lambda}_t - X_{it}^{\top}\hat{\beta}. \]
5.8.4.4 Uji dan Diagnosis
- Signifikansi efek tetap waktu: uji F untuk blok dummy waktu (LSDV) atau pembandingan FE satu arah vs dua arah.
- FE vs RE: Hausman test (menguji konsistensi asumsi RE).
- Signifikansi bersama kovariat: uji F pada \(\beta\).
- Multikolinieritas: hati-hati pada kovariat yang hampir konstan pada dimensi unit/waktu.
5.8.4.5 Implementasi di R
Berikut dua jalur populer: fixest (efisien, fleksibel)
dan plm (klasik untuk panel).
Dengan fixest
# install.packages("fixest") # jika belum terpasang
library(fixest)
# Contoh data: gunakan data 'Grunfeld' dari paket 'plm'
# install.packages("plm")
library(plm)
data("Grunfeld", package = "plm")
# Model FE dua arah: efek tetap perusahaan (firm) dan tahun (year)
# Cluster-robust SE by firm (bisa juga 'twoway' dengan fixest >= 0.11)
m_fe2_fixest <- feols(inv ~ value + capital | firm + year, data = Grunfeld, cluster = "firm")
summary(m_fe2_fixest)## OLS estimation, Dep. Var.: inv
## Observations: 200
## Fixed-effects: firm: 10, year: 20
## Standard-errors: Clustered (firm)
## Estimate Std. Error t value Pr(>|t|)
## value 0.117716 0.010824 10.87502 1.7727e-06 ***
## capital 0.357916 0.047848 7.48021 3.7703e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## RMSE: 47.5 Adj. R2: 0.943118
## Within R2: 0.720145Beberapa opsi SE lain (jika tersedia versinya) :
# Driscoll-Kraay (melalui vcovDC dari fixest)
summary(m_fe2_fixest, vcov = "DK")
# Two-way cluster (unit & waktu) jika diperlukan/tersedia
summary(m_fe2_fixest, cluster = c("firm","year"))Dengan plm
# install.packages("plm") # jika belum
library(plm)
# Spesifikasi FE dua arah (within transform)
m_fe2_plm <- plm(inv ~ value + capital, data = Grunfeld,
index = c("firm","year"),
model = "within", effect = "twoways")
summary(m_fe2_plm)## Twoways effects Within Model
##
## Call:
## plm(formula = inv ~ value + capital, data = Grunfeld, effect = "twoways",
## model = "within", index = c("firm", "year"))
##
## 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-16Standar error robust dengan plm:
# install.packages("sandwich"); install.packages("lmtest")
library(sandwich); library(lmtest)
# Cluster by firm
coeftest(m_fe2_plm, vcov = vcovHC(m_fe2_plm, method = "arellano", type = "HC1", cluster = "group"))##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## value 0.117716 0.009761 12.0599 < 2.2e-16 ***
## capital 0.357916 0.043147 8.2952 3.277e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Driscoll-Kraay
coeftest(m_fe2_plm, vcov = vcovSCC(m_fe2_plm, type = "HC1", maxlag = 2)) # maxlag pilih sesuai kebutuhan##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## value 0.117716 0.020539 5.7313 4.481e-08 ***
## capital 0.357916 0.056092 6.3809 1.629e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 15.8.4.5.1 Ringkasan Praktis
Spesifikasi model: \[ Y_{it} = \alpha_i + \lambda_t + X_{it}^{\top}\beta + \varepsilon_{it}. \]
Transformasi double-demeaning: \[ \tilde{Y}_{it} = Y_{it}-\bar{Y}_{i\cdot}-\bar{Y}_{\cdot t}+\bar{Y},\quad \tilde{X}_{it} = X_{it}-\bar{X}_{i\cdot}-\bar{X}_{\cdot t}+\bar{X}. \]
Estimator within dua arah: \[ \hat{\beta}_{\text{FE,2W}} = (\tilde{X}^{\top}\tilde{X})^{-1}\tilde{X}^{\top}\tilde{y}. \]
Pulihkan \(\hat{\alpha}_i\) dan \(\hat{\lambda}_t\) sesuai normalisasi, lalu evaluasi residual \(\hat{\varepsilon}_{it}\).
Laporkan SE yang sesuai (clustered/Driscoll–Kraay/two-way cluster) dan uji relevan.
5.8.5 Lampiran: FWL Ringkas (Opsional)
FWL menyatakan: jika Anda meregresikan \(Y\) dan setiap kolom \(X\) pada semua dummy unit dan waktu, lalu mengambil residualnya (\(\tilde{y}, \tilde{X}\)), maka OLS pada \(\tilde{y} \sim \tilde{X}\) menghasilkan \(\hat{\beta}\) yang sama dengan LSDV dan within.
Secara aljabar, dengan matriks proyeksi \(M\) yang menyingkirkan ruang yang direntang dummy unit dan waktu: \[ \hat{\beta} = (X^{\top} M X)^{-1} X^{\top} M y. \]
fe_tw <- plm(y ~ x1 + x2 + z, data = pdata, model = "within", effect = "twoways")
summary(fe_tw)## Twoways effects Within Model
##
## Call:
## plm(formula = y ~ x1 + x2 + z, data = pdata, effect = "twoways",
## model = "within")
##
## Balanced Panel: n = 100, T = 5, N = 500
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -2.6024620 -0.5481059 -0.0048051 0.5006134 2.7048882
##
## Coefficients:
## Estimate Std. Error t-value Pr(>|t|)
## x1 0.499424 0.016551 30.1748 < 2.2e-16 ***
## x2 -0.300839 0.032236 -9.3325 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 1274
## Residual Sum of Squares: 367.55
## R-Squared: 0.71149
## Adj. R-Squared: 0.63461
## F-statistic: 485.828 on 2 and 394 DF, p-value: < 2.22e-16coeftest(fe_tw, vcov = vcovHC(fe_tw, type = "HC1", cluster = "group"))##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## x1 0.499424 0.015766 31.677 < 2.2e-16 ***
## x2 -0.300839 0.029185 -10.308 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Alternatif manual:
lm(y ~ x1 + x2 + z + factor(id) + factor(year), data = panel), tapi FE viaplmlebih efisien.
5.8.5.0.1 Uji Signifikansi FE vs Pooled
- F-test FE: Apakah efek tetap unit signifikan?
- LM (Breusch–Pagan) nanti untuk RE vs pooled.
# One-way FE vs pooled
pFtest(fe, pooled)##
## F test for individual effects
##
## data: y ~ x1 + x2 + z
## F = -7.7887, df1 = -41, df2 = 397, p-value = NA
## alternative hypothesis: significant effects# Two-ways FE vs pooled (tak langsung, bisa bandingkan dengan model factor manual)
lm_pooled <- lm(y ~ x1 + x2 + z, data = panel)
lm_twoway <- lm(y ~ x1 + x2 + z + factor(id) + factor(year), data = panel)
anova(lm_pooled, lm_twoway) # indikasi kasar## Analysis of Variance Table
##
## Model 1: y ~ x1 + x2 + z
## Model 2: y ~ x1 + x2 + z + factor(id) + factor(year)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 356 666.95
## 2 293 211.05 63 455.9 10.046 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 15.9 Random Effects (RE)
RE memandang heterogenitas unit sebagai komponen acak: \(\alpha_i \sim \mathcal{N}(0,\sigma_{\alpha}^2)\) dan tidak berkorelasi dengan kovariat \(x_{it}\).
Model: \[ y_{it} = \beta_0 + \beta_1 x1_{it} + \beta_2 x2_{it} + \beta_3 z_t + \alpha_i + u_{it} \] dengan asumsi \(\text{Cov}(\alpha_i, x_{it}) = 0\).
re <- plm(y ~ x1 + x2 + z, data = pdata, model = "random", random.method = "swar")
summary(re)## Oneway (individual) effect Random Effect Model
## (Swamy-Arora's transformation)
##
## Call:
## plm(formula = y ~ x1 + x2 + z, data = pdata, model = "random",
## random.method = "swar")
##
## Balanced Panel: n = 100, T = 5, N = 500
##
## Effects:
## var std.dev share
## idiosyncratic 0.9311 0.9649 0.214
## individual 3.4204 1.8494 0.786
## theta: 0.7728
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -2.593311 -0.615791 -0.027389 0.660892 3.189176
##
## Coefficients:
## Estimate Std. Error z-value Pr(>|z|)
## (Intercept) 0.253760 0.343940 0.7378 0.4606
## x1 0.497239 0.016392 30.3343 < 2.2e-16 ***
## x2 -0.303237 0.031803 -9.5350 < 2.2e-16 ***
## z 0.193897 0.037318 5.1958 2.038e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 1397.6
## Residual Sum of Squares: 460.45
## R-Squared: 0.67055
## Adj. R-Squared: 0.66855
## Chisq: 1009.52 on 3 DF, p-value: < 2.22e-16coeftest(re, vcov = vcovHC(re, type = "HC1", cluster = "group"))##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.253760 0.310321 0.8177 0.4139
## x1 0.497239 0.015743 31.5856 < 2.2e-16 ***
## x2 -0.303237 0.029119 -10.4137 < 2.2e-16 ***
## z 0.193897 0.035276 5.4966 6.203e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 15.9.1 RE vs Pooled: LM Breusch–Pagan
Uji apakah varians komponen acak \(\sigma_\alpha^2\) non-zero (artinya RE beralasan).
plmtest(pooled, type = "bp")##
## Lagrange Multiplier Test - (Breusch-Pagan)
##
## data: y ~ x1 + x2 + z
## chisq = 249.51, df = 1, p-value < 2.2e-16
## alternative hypothesis: significant effects5.9.2 Hausman Test: FE vs RE
Kunci keputusan: bila \(\alpha_i\) berkorelasi dengan \(x_{it}\), FE konsisten tetapi RE tidak. Gunakan uji Hausman.
phtest(fe, re)##
## Hausman Test
##
## data: y ~ x1 + x2 + z
## chisq = 0.43892, df = 3, p-value = 0.9321
## alternative hypothesis: one model is inconsistent- Signifikan → tolak H0 (tidak ada korelasi); pilih
FE.
- Tidak signifikan → RE boleh dipakai (lebih efisien jika asumsi terpenuhi).
Dalam simulasi kita,
x1sengaja dibuat agak berkorelasi denganalpha_i; ekspektasi: Hausman signifikan → pilih FE.
5.10 Interpretasi Koefisien
- Pada FE, koefisien \(\beta\) diinterpretasi sebagai hubungan dalam-unit (how changes within the same unit relate to \(y\)).
- Pada RE dan Pooled, interpretasi mencampur variasi antara unit dan dalam unit — berpotensi bias jika ada korelasi \(x\) dengan \(\alpha_i\).
5.11 Interaksi dan Heterogenitas Slope
Model FE standar mengasumsikan slope sama untuk semua unit. Jika Anda
ingin slope berbeda (mis. x1 berefek lebih kuat di sektor
A), gunakan interaksi x1 * sector dan tetap sertakan efek
tetap unit/waktu bila relevan.
# Misal kita definisikan 3 sektor berdasarkan id
panel <- panel %>% mutate(sector = factor(ifelse(id <= 20, "A", ifelse(id <= 40, "B", "C"))))
pdata <- pdata.frame(panel, index = c("id","year"))
fe_int <- plm(y ~ x1*sector + x2 + z, data = pdata, model = "within", effect = "twoways")
summary(fe_int)## Twoways effects Within Model
##
## Call:
## plm(formula = y ~ x1 * sector + x2 + z, data = pdata, effect = "twoways",
## model = "within")
##
## Balanced Panel: n = 60, T = 6, N = 360
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -2.244464 -0.527023 -0.027684 0.461679 2.103522
##
## Coefficients:
## Estimate Std. Error t-value Pr(>|t|)
## x1 0.619239 0.050011 12.3821 < 2.2e-16 ***
## x2 -0.368322 0.049216 -7.4837 8.622e-13 ***
## x1:sectorB -0.032285 0.072973 -0.4424 0.6585
## x1:sectorC -0.034998 0.074177 -0.4718 0.6374
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 553.52
## Residual Sum of Squares: 210.84
## R-Squared: 0.61909
## Adj. R-Squared: 0.53008
## F-statistic: 118.239 on 4 and 291 DF, p-value: < 2.22e-16coeftest(fe_int, vcov = vcovHC(fe_int, type = "HC1", cluster = "group"))##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## x1 0.619239 0.049172 12.5934 < 2.2e-16 ***
## x2 -0.368322 0.044998 -8.1852 8.592e-15 ***
## x1:sectorB -0.032285 0.065465 -0.4932 0.6223
## x1:sectorC -0.034998 0.063755 -0.5489 0.5835
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Hati-hati ledakan parameter kalau level faktor banyak. Gunakan interaksi jika ada hipotesis teoritis yang jelas.
5.12 Visualisasi Efek Tetap: Within vs Between
Pisahkan variasi within (dalam unit) dan between
(antar unit) untuk x1.
wb <- panel %>%
group_by(id) %>%
mutate(x1_between = mean(x1), x1_within = x1 - x1_between) %>%
ungroup()
# Hubungan antar-unit (between): rata-rata y vs rata-rata x1 per unit
between_df <- wb %>% group_by(id) %>% summarize(y_m = mean(y), x1_m = mean(x1))
p1 <- ggplot(between_df, aes(x1_m, y_m)) +
geom_point() + geom_smooth(method = "lm", se = FALSE) +
labs(title = "Between Variation", x = "mean(x1_i)", y = "mean(y_i)") +
theme_minimal()
# Hubungan dalam-unit (within): deviasi dari mean unit
p2 <- ggplot(wb, aes(x1_within, y - ave(y, id, FUN=mean), color = factor(id))) +
geom_point(alpha = 0.2) +
geom_smooth(method = "lm", se = FALSE) +
guides(color = "none") +
labs(title = "Within Variation",
x = "x1 - mean(x1_i)", y = "y - mean(y_i)") +
theme_minimal()
p1; p2
5.13 Diagnostik Lanjutan
Panel sering menghadapi heteroskedastisitas, autokorelasi, dan ketergantungan silang.
5.13.1 Heteroskedastisitas & Autokorelasi
Gunakan SE robust berklaster pada id
(sudah dicontohkan). Bisa juga Driscoll–Kraay (ketahanan
terhadap dependensi silang) pada T moderat-besar.
# SE Driscoll-Kraay (butuh T cukup besar biasanya)
coeftest(fe_tw, vcov = vcovSCC(fe_tw, type = "HC1", maxlag = 2))##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## x1 0.4994243 0.0128222 38.950 < 2.2e-16 ***
## x2 -0.3008392 0.0087336 -34.446 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 15.13.2 Uji Serial Correlation (Breusch–Godfrey) dalam Panel
pbgtest(fe) # uji BG untuk model FE one-way##
## Breusch-Godfrey/Wooldridge test for serial correlation in panel models
##
## data: y ~ x1 + x2 + z
## chisq = 91.826, df = 5, p-value < 2.2e-16
## alternative hypothesis: serial correlation in idiosyncratic errors5.13.3 Uji Cross-Sectional Dependence
pcdtest(fe_tw, test = "cd") # Pesaran CD test##
## Pesaran CD test for cross-sectional dependence in panels
##
## data: y ~ x1 + x2 + z
## z = -1.1274, p-value = 0.2596
## alternative hypothesis: cross-sectional dependenceJika ada ketergantungan silang kuat (mis. shock nasional), pertimbangkan model dua-arah FE (sudah), faktor tak-teramati bersama (CMG), atau standar error Driscoll–Kraay.
5.14 Ringkasan Keputusan Model
- Mulai dengan Pooled sebagai baseline (dengan SE
robust).
- Uji LM (BP): jika signifikan → panel structure
penting (pertimbangkan RE/FE).
- Estimasi RE dan FE.
- Hausman: signifikan → pilih FE;
jika tidak → RE boleh.
- Tambahkan time FE bila ada tren/shock
agregat.
- Gunakan SE klaster atau
Driscoll–Kraay untuk robust inference.
- Plot dan eksplorasi variasi within vs between agar interpretasi tepat sasaran.
5.15 Praktik
- Skalakan/transform variabel yang sangat
skewed (contoh
x1log-normal) atau gunakan log jika interpretasi elastisitas relevan. - Centering membantu interpretasi intersep pada FE dengan interaksi kontinu–kategori.
- Hati-hati endogenitas: FE tidak menyelesaikan endogenitas time-varying (mis. reverse causality). Pertimbangkan IV panel atau desain kausal (DiD, event study).
- Panel tidak seimbang (unbalanced)
boleh dipakai
plm; dokumentasikan alasan kehilangan observasi. - Dokumentasikan indeks dengan jelas
(
id,year) dan pastikan merge/join variabel waktu-spesifik benar (hindari look-ahead bias).
5.16 Lampiran: Model Manual dengan Dummy Ekspisit
Untuk memperjelas ekivalensi FE ⇔ dummy penuh:
lm_manual <- lm(y ~ x1 + x2 + z + factor(id), data = panel) # one-way FE manual
summary(lm_manual)$coef[1:8,] %>% kable(digits=3) %>% kable_styling(full_width = FALSE)| Estimate | Std. Error | t value | Pr(>|t|) | |
|---|---|---|---|---|
| (Intercept) | 1.069 | 0.308 | 3.468 | 0.001 |
| x1 | 0.596 | 0.032 | 18.355 | 0.000 |
| x2 | -0.391 | 0.053 | -7.418 | 0.000 |
| z | -1.141 | 0.758 | -1.506 | 0.133 |
| factor(id)2 | 1.204 | 0.534 | 2.254 | 0.025 |
| factor(id)3 | 3.457 | 0.534 | 6.473 | 0.000 |
| factor(id)4 | 1.678 | 0.534 | 3.140 | 0.002 |
| factor(id)5 | 1.024 | 0.536 | 1.908 | 0.057 |
Koefisien
factor(id)adalah perbedaan intersep unit relatif baseline. Within estimator akan memberikan \(\hat{\beta}\) yang sama untuk kovariat.
5.17 Ekstensi Singkat
- Two-way FE: tambah
factor(year)ataueffect="twoways"(sudah dicontohkan). - Random slopes: izinkan slope berbeda per unit (model hierarkis/lme4) — beda filosofi dari FE/RE klasik.
- Dynamic panels: tambahkan lag
y_{it-1}→ butuh Arellano–Bond/GMM (plm::pgmm).
library(plm)
# pastikan pdata sudah: pdata <- pdata.frame(panel, index = c("id","year"))
# 1) Hanya y_lag sebagai regressor dinamis; instrumen GMM = lag(y, 2) saja
dyn_min <- pgmm(
formula = y ~ lag(y, 1) | lag(y, 2),
data = pdata,
effect = "individual",
model = "onestep", # 1-step lebih stabil saat matriks hampir singular
transformation = "d", # first-difference Arellano-Bond
collapse = TRUE # kurangi proliferasi instrumen
)
summary(dyn_min, robust = TRUE) # ringkasan + vcov robust## Oneway (individual) effect One-step model Difference GMM
##
## Call:
## pgmm(formula = y ~ lag(y, 1) | lag(y, 2), data = pdata, effect = "individual",
## model = "onestep", collapse = TRUE, transformation = "d")
##
## Balanced Panel: n = 60, T = 6, N = 360
##
## Number of Observations Used: 240
## Residuals:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -5.3838 -0.7915 0.4694 0.4386 1.8238 5.7426
##
## Coefficients:
## Estimate Std. Error z-value Pr(>|z|)
## lag(y, 1) -0.99908 0.53013 -1.8846 0.05949 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Sargan test: chisq(0) = 2.179288e-31 (p-value = < 2.22e-16)
## Autocorrelation test (1): normal = 0.1212099 (p-value = 0.90352)
## Autocorrelation test (2): normal = -2.030757 (p-value = 0.04228)
## Wald test for coefficients: chisq(1) = 3.551666 (p-value = 0.059486)mtest(dyn_min, order = 1) # AR(1) harus signifikan##
## Arellano-Bond autocorrelation test of degree 1
##
## data: y ~ lag(y, 1) | lag(y, 2)
## normal = 0.3248, p-value = 0.7453
## alternative hypothesis: autocorrelation presentmtest(dyn_min, order = 2) # AR(2) harus TIDAK signifikan##
## Arellano-Bond autocorrelation test of degree 2
##
## data: y ~ lag(y, 1) | lag(y, 2)
## normal = NaN, p-value = NA
## alternative hypothesis: autocorrelation presentsargan(dyn_min) # Hansen/Sargan idealnya tidak signifikan##
## Sargan test
##
## data: y ~ lag(y, 1) | lag(y, 2)
## chisq = 2.1793e-31, df = 0, p-value < 2.2e-16
## alternative hypothesis: overidentifying restrictions not valid5.18 Kesimpulan
- Pooled OLS: sederhana, tapi rawan bias jika ada heterogenitas unit yang berkorelasi dengan kovariat.
- Fixed Effects: mengontrol semua faktor tak-teramati konstan terhadap waktu pada level unit (dan/atau waktu). Andal saat ada korelasi \(\alpha_i\)–\(x_{it}\).
- Random Effects: efisien jika asumsi kunci (tidak ada korelasi \(\alpha_i\)–\(x_{it}\)) benar; uji Hausman adalah kompasnya.
- Struktur wide vs long memengaruhi workflow teknis: analisis panel modern hampir selalu di long.
- Gunakan robust/clustered SE, periksa serial correlation dan dependensi silang, serta pikirkan desain kausal untuk klaim sebab–akibat.
Dokumen ini mestinya menjadi toolbox ajar/riset. Anda bisa langsung knit untuk kelas, menyalin tabel ke slide, dan menyesuaikan simulasi ke konteks lokal (misalnya kabupaten–kota dan tahun kebijakan). Jika analisis makin “spicy”, jangan ragu menambahkan two-way FE dan tes Driscoll–Kraay—bumbu pedas yang aman untuk inferensi.
6 KUIS
6.1 Petunjuk
Jawablah dua soal aplikasi berikut secara singkat dan jelas. Gunakan notasi matematis atau interpretasi ekonometrika bila diperlukan.
6.2 Soal 1 — Aplikasi OLS (50 poin)
Seorang peneliti ingin menganalisis pengaruh jam belajar terhadap nilai ujian mahasiswa. Ia mengumpulkan data dari 100 mahasiswa dan memperkirakan model berikut:
\[ \text{Score}_i = \beta_0 + \beta_1 \text{StudyHours}_i + \varepsilon_i \]
Hasil estimasi menunjukkan bahwa \(\hat{\beta}_1 = 2.5\) dan signifikan pada tingkat 5%.
Jelaskan interpretasi dari koefisien \(\hat{\beta}_1 = 2.5\) dalam konteks penelitian ini. (25 poin)
Jika variabel kemampuan awal mahasiswa tidak dimasukkan ke dalam model, apa potensi masalah yang muncul? Jelaskan secara singkat. (25 poin)
6.3 Soal 2 — Aplikasi Fixed Effects (50 poin)
Seorang ekonom ingin meneliti pengaruh belanja pemerintah daerah terhadap tingkat kemiskinan di 10 provinsi selama 5 tahun. Model yang digunakan:
\[ \text{Poverty}_{it} = \alpha_i + \beta_1 \text{GovExp}_{it} + \varepsilon_{it} \]
Hasil model Fixed Effects menunjukkan bahwa \(\hat{\beta}_1 = -0.8\) dan signifikan pada tingkat 1%.
Jelaskan makna dari nilai \(\hat{\beta}_1 = -0.8\) dalam konteks penelitian ini. (25 poin)
Mengapa model Fixed Effects lebih sesuai dibanding pooled OLS dalam kasus ini? Berikan alasan singkat. (25 poin)
Kuis dikumpulan pada link gform berikut (10 November 2025 Pukul 16:00 WIB): https://forms.gle/QEwdZz6smK9eW3Jb8
7 Difference-in-Differences (DiD)
Model Difference-in-Differences (DiD) adalah teknik ekonometrika untuk mengukur efek kausal dari suatu kebijakan, intervensi, atau peristiwa terhadap variabel hasil (outcome), dengan memanfaatkan data sebelum vs sesudah pada kelompok perlakuan (treatment) dan kelompok kontrol (control).
Di dokumen ini, kita membedah DiD dari intuisi → formulasi → asumsi → implementasi R (termasuk event study sederhana), plus diagnostics dan uji placebo ringkas.
7.1 Intuisi Dasar
Bayangkan pemerintah memberi pelatihan keterampilan kepada sebagian provinsi (treatment), sementara provinsi lain tidak (control). Kita ingin tahu: apakah program ini meningkatkan pendapatan rata-rata?
Masalah: pendapatan bisa naik karena banyak hal (pertumbuhan ekonomi umum, inflasi, tren waktu). Jika kita hanya sebelum–sesudah pada treatment, efek bisa terkontaminasi tren umum.
Trik DiD: bandingkan perubahan di treatment dengan perubahan di kontrol:
\[ \text{Efek kebijakan} = \big(\Delta Y_{\text{treat}}\big) - \big(\Delta Y_{\text{control}}\big) \]
Kelompok kontrol bertindak sebagai counterfactual—dunia alternatif tanpa kebijakan.
7.2 Formulasi Matematis
Misalkan untuk unit \(i\) dan waktu \(t\):
- \(Y_{it}\): outcome (mis.
pendapatan)
- \(D_i \in \{0,1\}\): indikator
treatment (1 jika unit termasuk kelompok treatment)
- \(Post_t \in \{0,1\}\): indikator periode sesudah kebijakan
7.3 Model DiD Klasik
Model two-by-two dapat dituliskan sebagai:
\[ Y_{it} = \alpha + \beta D_i + \gamma Post_t + \delta (D_i \times Post_t) + \varepsilon_{it}. \]
Koefisien \(\delta\) merepresentasikan efek kausal kebijakan, yaitu selisih perubahan rata-rata outcome antara kelompok treatment dan control:
\[ \text{Efek kebijakan} = \big(\Delta Y_{\text{treat}}\big) - \big(\Delta Y_{\text{control}}\big) \]
Interpretasi
- \(\alpha\): rata-rata
outcome pada kelompok kontrol sebelum kebijakan
- \(\beta\): perbedaan tetap antara
treatment dan control sebelum kebijakan
- \(\gamma\): efek waktu umum
(terjadi di semua kelompok)
- \(\delta\): parameter utama (DiD effect) — estimasi dampak kebijakan hanya pada kelompok treatment setelah intervensi
7.4 Tabel Empat Sel
| Sebelum (0) | Sesudah (1) | Perubahan | |
|---|---|---|---|
| Kontrol (0) | \(\bar{Y}_{C,0}\) | \(\bar{Y}_{C,1}\) | \(\bar{Y}_{C,1} - \bar{Y}_{C,0}\) |
| Treatment (1) | \(\bar{Y}_{T,0}\) | \(\bar{Y}_{T,1}\) | \(\bar{Y}_{T,1} - \bar{Y}_{T,0}\) |
Efek DiD:
\[ \text{DiD} = (\bar{Y}_{T,1} - \bar{Y}_{T,0}) - (\bar{Y}_{C,1} - \bar{Y}_{C,0}) \]
7.5 Asumsi Kunci: Parallel Trends
Tanpa kebijakan, tren outcome di treatment dan kontrol sejajar (paralel):
\[ E[Y_{T,1} - Y_{T,0} \mid \text{no treatment}] = E[Y_{C,1} - Y_{C,0}] \]
Bila tren dasar tidak paralel, \(\hat{\delta}\) berpotensi
bias.
Karenanya, visualisasi dan uji lead (placebo) berguna sebagai
diagnostics.
7.6 Versi dengan Fixed Effects (TWFE)
Untuk data panel (banyak unit × banyak waktu), lazim digunakan Two-Way Fixed Effects (TWFE):
\[ Y_{it} = \alpha_i + \lambda_t + \delta (D_i \times Post_t) + \varepsilon_{it} \]
dengan:
- \(\alpha_i\): unit fixed
effect (mengontrol heterogenitas tetap antar unit),
- \(\lambda_t\): time fixed
effect (mengontrol guncangan waktu bersama),
- \(\delta\): tetap menjadi efek perlakuan (treatment effect).
7.7 Implementasi di R: Contoh Simulasi 2×2
Kita mulai dari skenario klasik 2×2: satu gelombang kebijakan, dua periode (0=sebelum, 1=sesudah), dua grup (kontrol vs treatment).
pkgs <- c("tidyverse", "fixest", "sandwich", "lmtest", "broom")
need <- pkgs[!pkgs %in% rownames(installed.packages())]
if (length(need)) install.packages(need, repos = "https://cloud.r-project.org")
invisible(lapply(pkgs, library, character.only = TRUE))set.seed(123)
N <- 200 # total unit
Tt <- 2 # dua periode: 0=pre, 1=post
id <- rep(1:N, each = Tt)
t <- rep(0:1, times = N)
# separuh jadi treatment:
D_i <- rep(rbinom(N, 1, 0.5), each = Tt)
Post_t <- t
# Komponen efek tetap & tren umum
alpha_i <- rnorm(N, 0, 1) # unit heterogeneity (random intercept) - akan di-absorb oleh FE
alpha_i <- rep(alpha_i, each = Tt)
time_shock <- ifelse(t==1, 1.0, 0.0) # ada kenaikan umum di periode post
# Efek perlakuan sebenarnya (ground truth)
delta_true <- 2.0
# Outcome dihasilkan dengan struktur DiD + noise
y <- 1.0 + 0.3*D_i + 0.8*Post_t + delta_true*(D_i*Post_t) + alpha_i + time_shock + rnorm(N*Tt, 0, 1)
did_2x2 <- tibble(id, t, D_i, Post_t, y) %>%
mutate(group = ifelse(D_i==1, "Treatment", "Control"),
period = ifelse(Post_t==1, "Post", "Pre"))
did_2x2 %>% head()## # A tibble: 6 × 7
## id t D_i Post_t y group period
## <int> <int> <int> <int> <dbl> <chr> <chr>
## 1 1 0 0 0 -0.426 Control Pre
## 2 1 1 0 1 1.34 Control Post
## 3 2 0 1 0 0.618 Treatment Pre
## 4 2 1 1 1 4.30 Treatment Post
## 5 3 0 0 0 0.316 Control Pre
## 6 3 1 0 1 2.88 Control Post7.7.1 Visualisasi Tren Rata-rata
avg <- did_2x2 %>%
group_by(group, period) %>%
summarise(mean_y = mean(y), .groups = "drop")
ggplot(avg, aes(x = period, y = mean_y, group = group)) +
geom_line(size = 1.1) + geom_point(size = 2) +
labs(x = NULL, y = "Rata-rata Y", color = NULL) +
theme_minimal(base_size = 12)
Rata-rata outcome per grup dan periode.
Intepretasi
Sumbu dan Makna Garis
- Sumbu X (horizontal): waktu — sebelum
(Pre) dan sesudah (Post) kebijakan.
- Sumbu Y (vertikal): nilai rata-rata variabel hasil
(outcome), misalnya pendapatan, skor ujian, atau tingkat
kesehatan.
- Dua garis:
- Garis atas biasanya mewakili kelompok
perlakuan (treatment) — unit yang terkena
kebijakan/intervensi.
- Garis bawah mewakili kelompok kontrol — unit yang tidak terkena kebijakan.
- Garis atas biasanya mewakili kelompok
perlakuan (treatment) — unit yang terkena
kebijakan/intervensi.
Apa yang Terlihat dari Grafik
Kedua garis naik dari Pre ke Post, tetapi garis treatment naik lebih curam. Artinya:
- Baik treatment maupun control mengalami peningkatan outcome (ada
time effect).
- Peningkatan pada kelompok treatment lebih besar — itulah efek kausal kebijakan (ΔDiD).
Secara matematis:
\[ \text{DiD} = (\bar{Y}_{T,\,post} - \bar{Y}_{T,\,pre}) - (\bar{Y}_{C,\,post} - \bar{Y}_{C,\,pre}). \]
Garis kontrol berfungsi sebagai counterfactual — apa yang akan terjadi pada kelompok treatment jika tidak ada kebijakan.
Cara Membaca Secara Intuitif
- Jika dua garis paralel sebelum intervensi, asumsi
parallel trends terpenuhi.
- Jika setelah intervensi garis treatment “meninggalkan” garis kontrol, jarak vertikal di titik Post menggambarkan efek kebijakan murni (setelah menghilangkan pengaruh waktu umum).
Ilustrasi Konseptual
Bayangkan:
- Semua provinsi naik pendapatannya sekitar +2 karena
pertumbuhan ekonomi nasional.
- Provinsi yang mendapat pelatihan (treatment) naik
+7.
Perbedaan tambahan +5 itulah efek DiD:
kenaikan ekstra karena kebijakan.
Contoh Reproducible (R)
Contoh kecil berikut mensimulasikan data 2×2 lalu menampilkan grafik mean Pre–Post untuk treatment vs control, lengkap dengan garis counterfactual (putus‑putus) yang menunjukkan “seandainya tidak ada kebijakan”.
set.seed(123)
library(dplyr)
library(ggplot2)
# N unit, 2 periode
n <- 200
dat <- expand.grid(id = 1:n, time = 0:1) %>% as_tibble()
# 50% unit sebagai treatment
dat <- dat %>%
mutate(treat = as.integer(id <= n/2))
# DGP (gunakan nrow(dat) agar aman)
dat <- dat %>%
mutate(
Y = 1.0 + 0.2*treat + 1.0*time + 0.8*(treat*time) + rnorm(nrow(dat), sd = 0.25)
)
# Plot
mean_trend <- dat %>%
group_by(treat, time) %>%
summarise(meanY = mean(Y), .groups = "drop")
ggplot(mean_trend, aes(x = factor(time, labels = c("Pre","Post")),
y = meanY, color = factor(treat))) +
geom_line(size = 1.2) +
geom_point(size = 3) +
labs(x = "Waktu", y = "Rata-rata Y", color = "Kelompok",
title = "Visualisasi Intuisi Difference-in-Differences") +
scale_color_discrete(labels = c("Kontrol","Treatment")) +
theme_minimal()
Kesimpulan Visual
Grafik menunjukkan bahwa kebijakan meningkatkan outcome kelompok treatment lebih besar daripada kontrol. Setelah mengeluarkan efek waktu umum, selisih tambahan itulah efek DiD (\(\Delta \text{DiD}\)).
7.7.2 Estimasi Model DiD Klasik (OLS)
m_ols <- lm(y ~ D_i + Post_t + D_i:Post_t, data = did_2x2)
summary(m_ols)$coefficients## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.1081474 0.1376832 8.0485281 9.875243e-15
## D_i 0.1920121 0.1977016 0.9712217 3.320307e-01
## Post_t 1.6339070 0.1947135 8.3913387 8.594261e-16
## D_i:Post_t 2.0130034 0.2795923 7.1997812 3.054983e-12Koefisien interaksi D_i:Post_t adalah \(\hat{\delta}\) (estimasi efek DiD). Untuk
standard errors yang robust & clustered (per
id), gunakan sandwich + lmtest:
# Cluster-robust SE by id
vcov_cl <- sandwich::vcovCL(m_ols, cluster = ~ id, type = "HC1")
coeftest(m_ols, vcov = vcov_cl)##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.10815 0.14832 7.4713 5.127e-13 ***
## D_i 0.19201 0.19574 0.9809 0.3272
## Post_t 1.63391 0.15350 10.6447 < 2.2e-16 ***
## D_i:Post_t 2.01300 0.21058 9.5591 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 17.7.3 Estimasi TWFE (Unit & Time FE)
# TWFE via fixest::feols
m_twfe <- fixest::feols(y ~ D_i:Post_t | id + t, data = did_2x2, cluster = ~ id)
summary(m_twfe)## OLS estimation, Dep. Var.: y
## Observations: 400
## Fixed-effects: id: 200, t: 2
## Standard-errors: Clustered (id)
## Estimate Std. Error t value Pr(>|t|)
## D_i:Post_t 2.013 0.210319 9.57118 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## RMSE: 0.741589 Adj. R2: 0.739123
## Within R2: 0.315115\(\hat{\delta}\) (koefisien pada
D_i:Post_t) adalah efek perlakuan setelah mengontrol
unit FE dan time FE, dengan
clustered SE pada level unit.
7.7.4 Diagnostik: Parallel Trends via Event-Study Sederhana
Pada 2×2, event study terbatas. Namun kita bisa mensimulasikan multi-periode dan memeriksa lead (efek sebelum kebijakan) sebagai placebo. Jika lead signifikan, indikasi pelanggaran tren paralel.
set.seed(42)
N <- 300
Tt <- 6 # lebih banyak periode
id <- rep(1:N, each = Tt)
t <- rep(1:Tt, times = N)
# separuh unit menerima treatment mulai periode 4 (statis, tidak bertahap)
D_i <- rep(rbinom(N, 1, 0.5), each = Tt)
Post_t <- as.integer(t >= 4)
# FE + tren umum
alpha_i <- rep(rnorm(N, 0, 1), each = Tt)
lambda_t <- rep(seq(0, 1.0, length.out = Tt), times = N) # tren umum linear ringan
delta_true <- 1.5
y <- 2 + 0.2*D_i + 0.6*Post_t + delta_true*(D_i*Post_t) + alpha_i + lambda_t + rnorm(N*Tt, 0, 1)
panel <- tibble(id, t, D_i, Post_t, y)7.7.5 Event-time & Leads/Lags
Definisikan event time untuk unit treatment: \(\tau = t - t\_0\) dengan \(t\_0=4\). Untuk kontrol, kita bisa set
NA (tidak pernah treat). Lalu kita bentuk variabel
lead/lag sebagai faktor dan estimasi model dengan
time FE & unit FE, mengabaikan
periode referensi (misal \(\tau=-1\)).
# --- PARAMETER EVENT TIME ---
t0 <- 4
# 1) Event time utk SEMUA unit (kontrol juga dihitung relatif t0),
# agar tidak ada NA yang bikin level hilang.
panel <- panel %>%
dplyr::mutate(
event_time_all = t - t0
)
# 2) Bucket lead/lag yang sama utk semua unit
panel <- panel %>%
dplyr::mutate(
event_bucket = dplyr::case_when(
event_time_all <= -3L ~ "lead_3ormore",
event_time_all == -2L ~ "lead_2",
event_time_all == -1L ~ "lead_1", # baseline yang DI-DROP sbg ref
event_time_all == 0L ~ "lag_0",
event_time_all == 1L ~ "lag_1",
event_time_all == 2L ~ "lag_2",
event_time_all >= 3L ~ "lag_3ormore"
),
event_bucket = factor(
event_bucket,
levels = c("lead_3ormore","lead_2","lead_1",
"lag_0","lag_1","lag_2","lag_3ormore")
)
)
# 3) Pastikan level referensi ADA pada treated; jika tidak, fallback otomatis
avail_ref <- levels(droplevels(panel$event_bucket[panel$D_i == 1L]))
ref_level <- if ("lead_1" %in% avail_ref) "lead_1" else
if ("lead_2" %in% avail_ref) "lead_2" else
if ("lead_3ormore" %in% avail_ref) "lead_3ormore" else
stop("Tidak ada level 'lead_x' pada unit treated. Periksa t0 & konstruksi data.")
# 4) ESTIMASI EVENT-STUDY (pakai SEMUA unit; FE: id + t; cluster: id)
m_es <- fixest::feols(
y ~ i(event_bucket, D_i, ref = ref_level) | id + t,
data = panel,
cluster = ~ id
)
# --- EKSTRAK KOEFISIEN: cari term yang dimulai "event_bucket::" ---
est <- broom::tidy(m_es) %>%
dplyr::filter(stringr::str_detect(term, "^event_bucket::")) %>%
dplyr::select(term, estimate, std.error, p.value)
est## # A tibble: 5 × 4
## term estimate std.error p.value
## <chr> <dbl> <dbl> <dbl>
## 1 event_bucket::lead_3ormore:D_i -0.208 0.165 2.08e- 1
## 2 event_bucket::lead_2:D_i -0.0311 0.164 8.50e- 1
## 3 event_bucket::lag_0:D_i 1.34 0.163 6.77e-15
## 4 event_bucket::lag_1:D_i 1.56 0.169 4.22e-18
## 5 event_bucket::lag_2:D_i 1.38 0.162 7.95e-16Interpretasi: koefisien pada leads (mis.
lead_3ormore, lead_2) seharusnya tidak berbeda
signifikan dari nol (tidak ada efek sebelum treatment).
Koefisien pada lags (lag_0, lag_1,
dst.) menunjukkan dinamika efek setelah kebijakan.
# Map nama term -> indeks waktu (-3, -2, 0, 1, 2, 3)
map_time <- c(
"lead_3ormore" = -3,
"lead_2" = -2,
# "lead_1" adalah REFERENCE -> tidak muncul sebagai koefisien
"lag_0" = 0,
"lag_1" = 1,
"lag_2" = 2,
"lag_3ormore" = 3
)
# Ambil nama bucket dari "event_bucket::<bucket>:D_i"
est_plot <- est %>%
dplyr::mutate(
bucket = term |>
stringr::str_remove("^event_bucket::") |>
stringr::str_remove(":D_i$"),
time = dplyr::case_when(
bucket == "lead_3ormore" ~ -3,
bucket == "lead_2" ~ -2,
# "lead_1" adalah baseline (ref), jadi TIDAK muncul di koefisien
bucket == "lag_0" ~ 0,
bucket == "lag_1" ~ 1,
bucket == "lag_2" ~ 2,
bucket == "lag_3ormore" ~ 3,
TRUE ~ NA_real_
),
ci_lo = estimate - 1.96*std.error,
ci_hi = estimate + 1.96*std.error
) %>%
dplyr::arrange(time)
# Cek cepat: kalau semua NA, berarti label-nya beda—lihat names(coef(m_es))
if (all(is.na(est_plot$time))) {
message("Label tidak terpetakan. Cek names(coef(m_es)) untuk pola aktual:")
print(names(coef(m_es)))
}
ggplot2::ggplot(est_plot, ggplot2::aes(x = time, y = estimate)) +
ggplot2::geom_hline(yintercept = 0, linetype = 2) +
ggplot2::geom_point(size = 2) +
ggplot2::geom_errorbar(ggplot2::aes(ymin = ci_lo, ymax = ci_hi), width = 0.1) +
ggplot2::scale_x_continuous(breaks = -3:3) +
ggplot2::labs(
x = "Event time (t - t0)",
y = paste0("ATT (ref = ", ref_level, ")"),
title = "Event-study (TWFE; SE dikelompokkan @ id)"
) +
ggplot2::theme_minimal(base_size = 12)
Event-study: koefisien lead/lag relatif terhadap periode referensi (lead_1).
7.8 Ringkasan Praktis
- Gunakan OLS interaksi untuk 2×2, tambahkan clustered SE (umumnya per unit) untuk inferensi yang valid.
- Untuk panel banyak periode, gunakan TWFE dengan unit dan time FE; cek parallel trends via event study (uji lead sebagai placebo).
- Waspadai staggered adoption (kebijakan diimplementasi pada waktu berbeda). Untuk itu, pertimbangkan estimator modern (Callaway–Sant’Anna; Sun–Abraham).
7.9 Lampiran: Varian Lanjutan (Opsional)
Di sini contoh staggered adoption minimalis memakai
fixest::sunab() (Sun–Abraham).
set.seed(99)
N <- 400
Tt <- 8
id <- rep(1:N, each = Tt)
t <- rep(1:Tt, times = N)
# Unit dibagi ke beberapa gelombang treatment (g): 3,5,7; sebagian tak pernah treated (g=Inf)
groups <- sample(c(3,5,7,Inf), N, replace = TRUE, prob = c(.3,.3,.2,.2))
g <- rep(groups, each = Tt)
treated <- as.integer(t >= g) # akan 0 bila g=Inf
alpha_i <- rep(rnorm(N, 0, 1), each = Tt)
lambda_t <- rep(sin(2*pi*(1:Tt)/Tt), times = N) # shock waktu non-linear
tau <- 1.2 # ATT rata-rata
y <- 1 + 0.5*treated + tau*treated + alpha_i + lambda_t + rnorm(N*Tt, 0, 1)
dat <- tibble(id, t, y, g, treated)
# Sun-Abraham 'sunab' with fixest
m_sa <- feols(y ~ sunab(g, t) | id + t, data = dat, cluster = ~ id)
iplot(m_sa, ref.line = 0) # plot event-time ATT7.10 Intuisi dan Ringkasan Konsep (Tambahan)
Difference-in-Differences (DiD) mengukur efek kausal kebijakan dengan membandingkan perubahan outcome pada kelompok treatment terhadap perubahan pada kelompok control dari sebelum ke sesudah kebijakan. Intuisi kasarnya:
\[ \text{Efek DiD} = (\Delta Y_{\text{treat}}) - (\Delta Y_{\text{control}}). \]
Asumsi kunci: parallel trends — tanpa kebijakan, tren outcome kedua kelompok akan bergerak sejajar.
7.11 Simulasi Data dan Visualisasi Tren
set.seed(123)
library(dplyr)
library(ggplot2)
# N unit, 2 periode (0 = sebelum, 1 = sesudah)
n <- 200
dat <- expand.grid(id = 1:n, time = 0:1) %>% as_tibble()
# 50% unit sebagai treatment
dat <- dat %>%
mutate(treat = as.integer(id <= n/2))
# DGP:
# baseline 10, beda baseline treatment 2, efek waktu umum 3,
# efek kebijakan (DiD) = 5 (hanya untuk treatment pada periode sesudah)
dat <- dat %>%
mutate(
Y = 10 + 2*treat + 3*time + 5*(treat*time) + rnorm(nrow(.), sd = 3)
# ^^^^^^^^^^ panjang noise tepat = jumlah baris data
)
# Rata-rata per grup-waktu untuk plot
mean_trend <- dat %>%
group_by(treat, time) %>%
summarise(meanY = mean(Y), .groups = "drop")
ggplot(mean_trend, aes(x = time, y = meanY, color = factor(treat))) +
geom_line(size = 1.2) +
geom_point(size = 3) +
scale_x_continuous(breaks = c(0,1), labels = c("Sebelum","Sesudah")) +
scale_color_discrete(labels = c("Kontrol","Treatment")) +
labs(x = "Waktu", y = "Rata-rata Y", color = "Kelompok",
title = "Visualisasi Intuisi Difference-in-Differences") +
theme_minimal()
# (Opsional) Sanity check: regresi DiD klasik
did_lm <- lm(Y ~ treat + time + treat:time, data = dat)
summary(did_lm)$coefficients## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.677360 0.2904595 33.317413 6.311389e-117
## treat 2.593858 0.4107718 6.314597 7.278946e-10
## time 3.213972 0.4107718 7.824227 4.696510e-14
## treat:time 4.876206 0.5809190 8.393951 8.433604e-16# Koefisien interaksi (treat:time) ≈ 5 (sesuai DGP)7.12 Estimasi Model DiD (Regresi Interaksi)
Model dasar DiD:
\[ Y_{it} = \alpha + \beta\,\text{treat}_i + \gamma\,\text{post}_t + \delta(\text{treat}_i \times \text{post}_t) + \varepsilon_{it}. \]
Di simulasi ini, \(\text{post}_t = \text{time}\). Koefisien interaksi (\(\delta\)) adalah efek DiD.
m_did <- lm(Y ~ treat + time + treat:time, data = dat)
broom::tidy(m_did)## # A tibble: 4 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 9.68 0.290 33.3 6.31e-117
## 2 treat 2.59 0.411 6.31 7.28e- 10
## 3 time 3.21 0.411 7.82 4.70e- 14
## 4 treat:time 4.88 0.581 8.39 8.43e- 16Interpretasi cepat: koefisien treat:time adalah estimasi
efek kebijakan (DiD).
7.13 Standard Error Berkelompok
Pada praktiknya, error sering berkorelasi dalam unit/kelompok.
Gunakan clustered robust SE (misalnya dikelompokkan
pada id). Jika paket sandwich dan
lmtest tersedia, hasil dengan SE berkelompok akan
ditampilkan.
if (requireNamespace("sandwich", quietly = TRUE) && requireNamespace("lmtest", quietly = TRUE)) {
vcov_cl <- sandwich::vcovCL(m_did, cluster = dat$id, type = "HC1")
lmtest::coeftest(m_did, vcov = vcov_cl)
} else {
message("Lewati SE berkelompok: paket 'sandwich' dan/atau 'lmtest' tidak tersedia.")
}##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.67736 0.29046 33.3173 < 2.2e-16 ***
## treat 2.59386 0.39943 6.4939 2.511e-10 ***
## time 3.21397 0.41687 7.7097 1.029e-13 ***
## treat:time 4.87621 0.59210 8.2355 2.631e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 17.14 Versi Panel dengan Fixed Effects (TWFE)
Untuk data panel (unit berulang antar waktu), model two-way fixed effects (TWFE) menambahkan efek tetap unit dan waktu:
\[ Y_{it} = \alpha_i + \lambda_t + \delta(\text{treat}_i \times \text{post}_t) + \varepsilon_{it}. \]
Implementasi sederhana (demeaning manual) untuk memberi gambaran ide dasarnya:
# Buat faktor unit dan waktu
dat <- dat %>% mutate(id_f = factor(id), t_f = factor(time))
# Estimasi TWFE via regresi dengan dummy unit & waktu (tanpa intercept duplikat)
m_twfe <- lm(Y ~ treat:time + id_f + t_f, data = dat)
broom::tidy(m_twfe) %>%
filter(term == "treat:time") # sorot koefisien interaksi utama## # A tibble: 1 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 treat:time 4.88 0.591 8.25 2.24e-14Catatan: Untuk aplikasi serius, pertimbangkan paket khusus seperti fixest atau lfe untuk efisiennya estimasi FE skala besar serta ragam opsi SE.
7.15 Pemeriksaan (Kualitatif) Parallel Trends
Dengan hanya dua periode, pemeriksaan bersifat kualitatif melalui plot rata-rata per grup-waktu. Untuk banyak periode, gunakan event study agar bisa melihat dinamika pra- dan pasca-kebijakan.
ggplot(mean_trend, aes(x = time, y = meanY, group = factor(treat), linetype = factor(treat))) +
geom_line(size = 1.1) + geom_point(size = 3) +
scale_x_continuous(breaks = c(0,1), labels = c("Sebelum","Sesudah")) +
labs(title = "Pemeriksaan Kualitatif Parallel Trends",
x = "Waktu", y = "Rata-rata Outcome", linetype = "Kelompok") +
theme_bw()
7.16 Ringkasan dan Praktik Baik
- Kunci DiD: koefisien interaksi adalah estimasi efek kebijakan.
- Asumsi inti: parallel trends. Uji/pertimbangkan via pra-tren (bila banyak periode).
- SE berkelompok: hampir wajib saat panel/clustered sampling.
- Hati-hati saat ada staggered treatment (waktu treatment berbeda): gunakan metode modern (Callaway–Sant’Anna, Sun–Abraham, Synthetic DiD) untuk menghindari bias agregasi.
7.17 Lampiran: Replikasi Rumus Empat Sel
Rumus DiD empat sel:
\[ (\bar{Y}_{T,1}-\bar{Y}_{T,0})-(\bar{Y}_{C,1}-\bar{Y}_{C,0}). \]
cells <- dat %>%
group_by(treat, time) %>%
summarise(meanY = mean(Y), .groups = "drop") %>%
pivot_wider(names_from = time, values_from = meanY, names_prefix = "t")
cells## # A tibble: 2 × 3
## treat t0 t1
## <int> <dbl> <dbl>
## 1 0 9.68 12.9
## 2 1 12.3 20.4did_four_cells <- (cells %>% filter(treat == 1) %>% pull(t1) -
cells %>% filter(treat == 1) %>% pull(t0)) -
(cells %>% filter(treat == 0) %>% pull(t1) -
cells %>% filter(treat == 0) %>% pull(t0))
did_four_cells## [1] 4.876206Mengevaluasi efek kebijakan/perlakuan sebelum–sesudah pada kelompok perlakuan vs kontrol. Asumsi kunci: parallel trends.
\[ Y_{it} = \alpha + \delta D_i + \gamma T_t + \beta(D_i \times T_t) + \varepsilon_{it}. \]
set.seed(42)
n_g <- 100
group <- rep(c(0,1), each = n_g) # 0: kontrol, 1: perlakuan
time <- rep(rep(c(0,1), each = n_g/2), 2) # 0: sebelum, 1: sesudah
tau <- 2 # efek perlakuan sejati
y <- 5 + 1*group + 1*time + tau*(group*time) + rnorm(2*n_g, 0, 1)
did_df <- tibble(y, group, time)
did_model <- fixest::feols(y ~ group*time, data = did_df)
summary(did_model)## OLS estimation, Dep. Var.: y
## Observations: 200
## Standard-errors: IID
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.964328 0.138292 35.89743 < 2.2e-16 ***
## group 0.884421 0.195574 4.52217 1.0585e-05 ***
## time 1.136373 0.195574 5.81044 2.4874e-08 ***
## group:time 1.991162 0.276584 7.19912 1.2635e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## RMSE: 0.968044 Adj. R2: 0.703117Interpretasi. Koefisien interaksi
group:time adalah taksiran efek perlakuan rata-rata.
8 Synthetic Control Model (SCM)
Synthetic Control Method (SCM) adalah metode kuasi-eksperimen untuk mengevaluasi dampak kausal kebijakan/peristiwa pada satu unit yang terpapar (treated unit) ketika hanya tersedia beberapa unit pembanding. Inti gagasannya: membangun “kontrol sintetis” sebagai kombinasi cembung dari unit-unit donor (kontrol) sehingga jejak (trajectory) pra-intervensi dari treated tercocokkan sebaik mungkin.
Secara informal, Synthetic Control Method (SCM) bertujuan memperkirakan counterfactual \(Y_{1t}(0)\)—yakni hasil yang mungkin terjadi pada unit treated jika tidak ada intervensi atau kebijakan. Estimatornya didefinisikan sebagai:
\[ \hat{Y}^{\text{SC}}_{1t}(0) = \sum_{j=2}^{J+1} w_j \, Y_{jt}, \]
dengan syarat \(w_j \ge 0\) dan
\(\sum_{j=2}^{J+1} w_j = 1\).
Bobot \(w_j\) dipilih sedemikian rupa
sehingga kombinasi linier dari unit kontrol (\(j = 2, \dots, J+1\)) menghasilkan
hasil dan predictors pra-intervensi yang sedekat
mungkin dengan milik unit treated.
8.1 Notasi & Estimand
- Unit \(i=1\) treated pada \(t \ge T_0\) (setelah kebijakan); unit \(i=2,\dots,J+1\) adalah donor pool (kontrol).
- Outcome teramati: \(Y_{it}\).
- Efek kausal (ATT) pada waktu \(t \ge T_0\): \[ \tau_t = Y_{1t}(1) - Y_{1t}(0) \approx Y_{1t} - \hat Y^{\text{SC}}_{1t}(0). \]
SCM meminimalkan jarak pra-intervensi antara treated dan kombinasi sintetis terhadap sekumpulan predictors (rata-rata pra, lag outcome, kovariat). Matriks \(V\) (bobot untuk predictors) dan vektor \(w\) (bobot donor) dipilih melalui prosedur dua tingkat (nested optimization).
8.2 Estimasi (Garis Besar)
- Pilih donor pool (unit kontrol yang relevan, tidak
ikut terpapar).
- Tentukan periode pra-intervensi \([1, \dots,
T_0-1]\) dan pasca-intervensi \([T_0, \dots,
T]\).
- Susun predictors: ringkasan outcome pra
(rata-rata, lag tertentu) dan kovariat exogenous (opsional).
- Optimasi bobot \(w\) untuk meminimalkan jarak
pra-intervensi sesuai bobot predictors \(V\).
- Hitung kontrol sintetis dan efek
\(\tau_t\) untuk \(t \ge T_0\).
- Diagnostik & inferensi: path plot, gaps plot, placebo (permutation) test, MSPE (Mean Squared Prediction Error) pra vs pasca, leave-one-out.
8.3 Analisis Dengan R
This R Markdown is adapted from Danilo Freire’s “A Simple SCM Tutorial” (2018). Tujuannya: file yang siap knit tanpa error umum seperti
unit.variable not found as numeric variable in foo.
8.3.1 Instal paket (opsional, jalankan sekali saja)
# Jalankan baris-baris ini jika paket belum terpasang (hapus komentar #)
# install.packages("Synth", repos = "https://cloud.r-project.org")
# install.packages("gsynth", repos = "https://cloud.r-project.org")8.3.2 Muat paket dan set opsi
library(Synth)
library(gsynth)
# Reproducibility
set.seed(1)
# Knitr options (opsional)
knitr::opts_chunk$set(message = FALSE, warning = FALSE, dpi = 120)8.3.3 Simulasi Data
Kita buat data untuk 10 state dan 30 tahun. State
A menerima perlakuan mulai tahun 15. Kita
pastikan kolom pengenal unit numeric agar tidak muncul
error
unit.variable not found as numeric variable in foo.
year <- rep(1:30, 10)
state <- rep(LETTERS[1:10], each = 30)
X1 <- round(rnorm(300, mean = 2, sd = 1), 2)
X2 <- round(rbinom(300, 1, 0.5) + rnorm(300), 2)
Y <- round(1 + 2*X1 + rnorm(300), 2)
df <- data.frame(Y = Y, X1 = X1, X2 = X2, state = state, year = year,
stringsAsFactors = FALSE)
# Pastikan tipe data yang sesuai
df$Y <- as.numeric(df$Y)
df$X1 <- as.numeric(df$X1)
df$X2 <- as.numeric(df$X2)
df$year <- as.numeric(df$year)
df$state <- as.character(df$state)
# Kode numerik untuk unit (WAJIB numeric untuk Synth::dataprep)
df$state.num <- as.numeric(factor(df$state, levels = LETTERS[1:10]))
# Perlakuan untuk State A mulai tahun 15
df$T <- ifelse(df$state == "A" & df$year >= 15, 1, 0)
# Tambahkan efek treatment sebesar 20 pada outcome setelah tahun 15 untuk A
df$Y <- ifelse(df$state == "A" & df$year >= 15, df$Y + 20, df$Y)
# Cek struktur
str(df)## 'data.frame': 300 obs. of 7 variables:
## $ Y : num 2.29 4.51 2.07 8.87 4.37 1.32 8 7.49 6.98 3.72 ...
## $ X1 : num 1.37 2.18 1.16 3.6 2.33 1.18 2.49 2.74 2.58 1.69 ...
## $ X2 : num 1.96 0.4 -0.75 -0.56 -0.45 1.06 0.51 -2.1 0 0.54 ...
## $ state : chr "A" "A" "A" "A" ...
## $ year : num 1 2 3 4 5 6 7 8 9 10 ...
## $ state.num: num 1 1 1 1 1 1 1 1 1 1 ...
## $ T : num 0 0 0 0 0 0 0 0 0 0 ...head(df)## Y X1 X2 state year state.num T
## 1 2.29 1.37 1.96 A 1 1 0
## 2 4.51 2.18 0.40 A 2 1 0
## 3 2.07 1.16 -0.75 A 3 1 0
## 4 8.87 3.60 -0.56 A 4 1 0
## 5 4.37 2.33 -0.45 A 5 1 0
## 6 1.32 1.18 1.06 A 6 1 08.3.4 Estimasi dengan Synth
Kita siapkan objek dataprep dan jalankan
synth. Perhatikan bahwa: - unit.variable
harus numeric. - treatment.identifier
di-set ke 1 (karena A dikodekan
menjadi 1). - controls.identifier adalah 2 hingga 10.
dataprep.out <- Synth::dataprep(
foo = df,
predictors = c("X1", "X2"),
dependent = "Y",
unit.variable = "state.num",
time.variable = "year",
unit.names.variable = "state",
treatment.identifier = 1,
controls.identifier = 2:10,
time.predictors.prior = 1:14,
time.optimize.ssr = 1:14,
time.plot = 1:30
)
synth.out <- Synth::synth(dataprep.out)##
## X1, X0, Z1, Z0 all come directly from dataprep object.
##
##
## ****************
## searching for synthetic control unit
##
##
## ****************
## ****************
## ****************
##
## MSPE (LOSS V): 9.831789
##
## solution.v:
## 0.3888387 0.6111613
##
## solution.w:
## 0.1115941 0.1832781 0.1027237 0.312091 0.06096758 0.03509706 0.05893735 0.05746256 0.07784853# Tabel ringkas
synth.tables <- Synth::synth.tab(
dataprep.res = dataprep.out,
synth.res = synth.out
)
synth.tables## $tab.pred
## Treated Synthetic Sample Mean
## X1 2.028 2.028 2.017
## X2 0.513 0.513 0.394
##
## $tab.v
## v.weights
## X1 0.389
## X2 0.611
##
## $tab.w
## w.weights unit.names unit.numbers
## 2 0.112 B 2
## 3 0.183 C 3
## 4 0.103 D 4
## 5 0.312 E 5
## 6 0.061 F 6
## 7 0.035 G 7
## 8 0.059 H 8
## 9 0.057 I 9
## 10 0.078 J 10
##
## $tab.loss
## Loss W Loss V
## [1,] 9.761706e-12 9.8317898.3.5 Plot Jalur (Path Plot)
path.plot(
synth.res = synth.out,
dataprep.res = dataprep.out,
Ylab = "Y",
Xlab = "Year",
Legend = c("State A", "Synthetic State A"),
Legend.position = "topleft"
)
abline(v = 15, lty = 2)
8.3.6 Plot Gap
gaps.plot(
synth.res = synth.out,
dataprep.res = dataprep.out,
Ylab = "Gap",
Xlab = "Year",
Ylim = c(-30, 30),
Main = ""
)
abline(v = 15, lty = 2)
8.3.7 Estimasi dengan gsynth
gsynth memungkinkan multiple treated units dan
komponen faktor laten (iterative fixed effects). Kita
gunakan two-way fixed effects dengan cross-validation
untuk memilih jumlah faktor r. Agar knitting mulus di
berbagai lingkungan, kita set parallel = FALSE dan
mengurangi nboots.
gs.out <- gsynth(
formula = Y ~ T + X1 + X2,
data = df,
index = c("state", "year"),
force = "two-way",
CV = TRUE,
r = c(0, 5),
se = TRUE,
inference = "parametric",
nboots = 200, # kurangi untuk mempercepat knit
parallel = FALSE # TRUE di mesin lokal jika ingin lebih cepat
)## Cross-validating ...
## r = 0; sigma2 = 1.13533; IC = 0.95632; PC = 0.96713; MSPE = 1.65502
## r = 1; sigma2 = 0.96885; IC = 1.54420; PC = 4.30644; MSPE = 1.33375
## r = 2; sigma2 = 0.81855; IC = 2.08062; PC = 6.58556; MSPE = 1.27341*
## r = 3; sigma2 = 0.71670; IC = 2.61125; PC = 8.35187; MSPE = 1.79319
## r = 4; sigma2 = 0.62823; IC = 3.10156; PC = 9.59221; MSPE = 2.02301
## r = 5; sigma2 = 0.55497; IC = 3.55814; PC = 10.48406; MSPE = 2.79596
##
## r* = 2
##
## Simulating errors .....Bootstrapping ...
## ..gs.out## Call:
## gsynth.formula(formula = Y ~ T + X1 + X2, data = df, index = c("state",
## "year"), force = "two-way", r = c(0, 5), CV = TRUE, se = TRUE,
## nboots = 200, inference = "parametric", parallel = FALSE)
##
## Average Treatment Effect on the Treated:
## Estimate S.E. CI.lower CI.upper p.value
## ATT.avg 20.45 0.7141 19.06 21.85 0
##
## ~ by Period (including Pre-treatment Periods):
## ATT S.E. CI.lower CI.upper p.value n.Treated
## -13 -0.43323 1.0923 -2.5741 1.7076 6.916e-01 0
## -12 -0.78276 0.9471 -2.6390 1.0735 4.085e-01 0
## -11 -0.19199 0.9309 -2.0165 1.6325 8.366e-01 0
## -10 -0.98289 1.0057 -2.9540 0.9882 3.284e-01 0
## -9 -0.94504 0.8148 -2.5419 0.6519 2.461e-01 0
## -8 -0.23970 0.9051 -2.0136 1.5342 7.911e-01 0
## -7 1.74468 1.1957 -0.5988 4.0882 1.445e-01 0
## -6 1.73998 1.2365 -0.6836 4.1636 1.594e-01 0
## -5 0.63689 0.9037 -1.1343 2.4081 4.810e-01 0
## -4 -0.03246 0.8872 -1.7713 1.7063 9.708e-01 0
## -3 -0.25405 0.8082 -1.8380 1.3299 7.533e-01 0
## -2 0.10273 0.8376 -1.5390 1.7445 9.024e-01 0
## -1 -0.23391 0.5835 -1.3776 0.9097 6.885e-01 0
## 0 -0.12827 0.9169 -1.9253 1.6688 8.887e-01 0
## 1 21.19473 1.7697 17.7263 24.6632 0.000e+00 1
## 2 22.84214 1.2183 20.4544 25.2299 0.000e+00 1
## 3 20.36902 1.2826 17.8552 22.8829 0.000e+00 1
## 4 22.33201 1.9629 18.4848 26.1793 0.000e+00 1
## 5 19.68346 2.4703 14.8417 24.5252 1.554e-15 1
## 6 18.08863 0.8339 16.4541 19.7231 0.000e+00 1
## 7 19.75869 1.3372 17.1378 22.3796 0.000e+00 1
## 8 18.92529 1.9750 15.0543 22.7963 0.000e+00 1
## 9 21.74294 0.9219 19.9360 23.5499 0.000e+00 1
## 10 21.64782 0.9279 19.8291 23.4665 0.000e+00 1
## 11 19.39832 1.3055 16.8396 21.9570 0.000e+00 1
## 12 19.31785 0.8282 17.6946 20.9411 0.000e+00 1
## 13 22.28342 1.1193 20.0895 24.4773 0.000e+00 1
## 14 20.38876 0.9184 18.5888 22.1888 0.000e+00 1
## 15 18.92550 1.5937 15.8019 22.0491 0.000e+00 1
## 16 20.37687 1.6095 17.2223 23.5314 0.000e+00 1
##
## Coefficients for the Covariates:
## beta S.E. CI.lower CI.upper p.value
## X1 1.93082 0.06295 1.80745 2.0542 0.0000
## X2 0.03809 0.04610 -0.05227 0.1284 0.4087# Plot utama
plot(gs.out)
# Counterfactual (trajektori aktual vs counterfactual untuk unit treated)
plot(gs.out, type = "counterfactual")
# Estimasi untuk semua kontrol (raw = "all")
plot(gs.out, type = "counterfactual", raw = "all")
8.3.8 Menghindari Error Umum
unit.variable not found as numeric variable in foo
Pastikan kolom yang dipakai sebagaiunit.variablenumeric (contoh:state.num). Jika masih gagal, cek apakah nama kolom sesuai (ejaan tepat).X1,X0,Z1,Z0 were individually input(pesan informasi)
Itu normal saatsynth()dipanggil menggunakan objekdataprep.Plot tidak muncul saat knit
Pastikan chunk tidakecho=FALSEdan tidakfig.show='hide'.
9 Vector Autoregressive Model (VAR)
9.1 Apa itu VAR?
Vector Autoregression (VAR) adalah model untuk beberapa variabel time‑series yang saling mempengaruhi secara dinamis. Setiap variabel dijelaskan oleh lag dirinya dan lag variabel lain. Intuisinya: ekonomi itu ngobrol—bukan monolog—jadi modelnya juga multi-way conversations.
Secara ringkas, VAR cocok ketika:
- kita punya beberapa seri yang mungkin saling memengaruhi,
- tidak ingin memaksakan struktur kausal yang kaku di awal,
- fokus pada peramalan, impulse response, dan variance decomposition.
9.2 Spesifikasi Model
Untuk vektor \(y_t\in\mathbb{R}^K\) dan lag orde \(p\), model VAR(p): \[ y_t = c + A_1 y_{t-1} + A_2 y_{t-2} + \cdots + A_p y_{t-p} + u_t,\quad u_t\sim \text{i.i.d. }(0,\,\Sigma_u), \] dengan: - \(c\) vektor intersep \((K\times 1)\), - \(A_j\) matriks koefisien \((K\times K)\), - \(u_t\) inovasi white noise kovarian \(\Sigma_u\).
Bentuk persamaan per variabel (misal \(K=2\)): \[ \begin{aligned} y_{1t} &= c_1 + \sum_{j=1}^p a^{(11)}j\,y{1,t-j} + \sum_{j=1}^p a^{(12)}j\,y{2,t-j} + u_{1t},\\ y_{2t} &= c_2 + \sum_{j=1}^p a^{(21)}j\,y{1,t-j} + \sum_{j=1}^p a^{(22)}j\,y{2,t-j} + u_{2t}. \end{aligned} \]
9.3 Ide yang Melandasi
- Endogenitas bersama: semua variabel dianggap endogen dalam sistem.
- Parsimonious but rich: cukup fleksibel menangkap dinamika silang, tanpa mengasumsikan arah kausal a priori.
- Analisis dinamika kebijakan/shock: lewat Impulse Response Function (IRF) dan Forecast Error Variance Decomposition (FEVD).
9.4 Asumsi Penting
- Stasioneritas (weak/kovarian): nilai harapan
konstan, varians konstan, kovarians hanya fungsi dari lag. Secara
praktis: akar polinomial karakteristik (akar dari companion
matrix) berada di luar lingkaran satuan (|root|
> 1).
- Inovasi white noise: \(u_t\) tak berkorelasi
serial, \(\mathbb{E}(u_t)=0\), kovarian konstan \(\Sigma_u\).
- Tidak ada perfect multicollinearity di
regresor (lag-lag).
- Ukuran sampel memadai relatif terhadap parameter (VAR bisa cepat “haus data”).
Bila data punya unit root (non‑stationary) dan tidak kointegrasi, gunakan differences. Bila terkointegrasi, gunakan VECM (VAR terbatas) dan, bila perlu, transform ke SVAR/IRF struktural.
9.5 Estimasi Parameter
- VAR diestimasi persamaan-per-persamaan memakai OLS karena regresor sama di setiap persamaan (blok lag vektor \(y\)).
- OLS konsisten & efisien (di kelas seemingly unrelated regressions dengan identik regresor).
9.6 Pemilihan Orde Lag
Pilih \(p\) via kriteria informasi: AIC, BIC (SC), HQ. Secara umum: - AIC cenderung memilih model lebih besar (baik untuk forecast), - BIC/SC lebih hemat parameter (baik untuk penjelasan).
9.7 Diagnostik & Stabilitas
- Serial correlation inovasi: uji Portmanteau /
LM.
- Normalitas residu (opsional untuk inferensi
tertentu).
- ARCH effects (heteroskedastisitas bersyarat).
- Stabilitas: periksa akar (roots) VAR (harus di luar unit circle).
9.8 Granger Causality
Variabel \(x\) Granger-causes \(y\) bila lag dari \(x\) meningkatkan kemampuan memprediksi \(y\) di luar lag \(y\) sendiri (uji Wald pada blok koefisien lag \(x\) pada persamaan \(y\)).
9.9 IRF & FEVD
- IRF: respons dinamis variabel terhadap
kejutan satu-satuan pada salah satu inovasi. Karena
inovasi saling berkorelasi, perlu ortogonalization
(umum: Cholesky — urutan variabel penting!).
- FEVD: proporsi varians kesalahan ramalan pada horizon tertentu yang dijelaskan oleh setiap inovasi.
9.10 Contoh Kasus di R
(paket vars)
Kita gunakan data Canada bawaan paket
vars (indikator makro Kanada: e (employment),
prod (productivity), rw (real wages),
U (unemployment)). Agar ringkas dan stasioner, kita pakai
perubahan log (growth rates) untuk tiga seri pertama,
dan diferensiasi U.
library(vars)
library(tidyverse)
library(urca) # unit root tests9.10.1 Siapkan Data
data("Canada", package = "vars") # ts quarterly
# Transformasi: growth rate (diff log) untuk e, prod, rw; dan diff untuk U
can_df <- as_tibble(Canada) |>
mutate(
le = log(e),
lprod= log(prod),
lrw = log(rw)
) |>
mutate(
d_le = c(NA, diff(le)),
d_lprod = c(NA, diff(lprod)),
d_lrw = c(NA, diff(lrw)),
d_U = c(NA, diff(U))
) |>
drop_na(d_le, d_lprod, d_lrw, d_U)
glimpse(can_df)## Rows: 83
## Columns: 11
## $ e <dbl> 929.8040, 930.3184, 931.4277, 932.6620, 933.5509, 933.5315, 93…
## $ prod <dbl> 404.6398, 403.8149, 404.2158, 405.0467, 404.4167, 402.8191, 40…
## $ rw <dbl> 388.1358, 390.5401, 393.9638, 396.7647, 400.0217, 400.7515, 40…
## $ U <dbl> 7.70, 7.47, 7.27, 7.37, 7.13, 7.40, 8.33, 8.83, 10.43, 12.20, …
## $ le <dbl> 6.834974, 6.835527, 6.836719, 6.838043, 6.838996, 6.838975, 6.…
## $ lprod <dbl> 6.002997, 6.000957, 6.001949, 6.004002, 6.002446, 5.998488, 5.…
## $ lrw <dbl> 5.961355, 5.967531, 5.976259, 5.983343, 5.991519, 5.993342, 6.…
## $ d_le <dbl> 2.080985e-04, 5.530851e-04, 1.191677e-03, 1.324312e-03, 9.5266…
## $ d_lprod <dbl> -1.794141e-03, -2.040810e-03, 9.922648e-04, 2.053575e-03, -1.5…
## $ d_lrw <dbl> 0.0051652516, 0.0061755129, 0.0087283848, 0.0070843162, 0.0081…
## $ d_U <dbl> 0.17, -0.23, -0.20, 0.10, -0.24, 0.27, 0.93, 0.50, 1.60, 1.77,…9.10.2 Uji Akar Unit (sekilas)
adf_le <- ur.df(can_df$d_le, type = "drift", selectlags = "AIC")
adf_lprod <- ur.df(can_df$d_lprod, type = "drift", selectlags = "AIC")
adf_lrw <- ur.df(can_df$d_lrw, type = "drift", selectlags = "AIC")
adf_U <- ur.df(can_df$d_U, type = "drift", selectlags = "AIC")
summary(adf_le)##
## ###############################################
## # Augmented Dickey-Fuller Test Unit Root Test #
## ###############################################
##
## Test regression drift
##
##
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.716e-04 -2.668e-04 -2.083e-05 2.900e-04 7.841e-04
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.491e-04 5.622e-05 2.651 0.00971 **
## z.lag.1 -3.595e-01 7.927e-02 -4.536 2.05e-05 ***
## z.diff.lag 3.241e-01 1.076e-01 3.011 0.00350 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.0004141 on 78 degrees of freedom
## Multiple R-squared: 0.2251, Adjusted R-squared: 0.2052
## F-statistic: 11.33 on 2 and 78 DF, p-value: 4.789e-05
##
##
## Value of test-statistic is: -4.536 10.2891
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau2 -3.51 -2.89 -2.58
## phi1 6.70 4.71 3.86Dari praktik umum, growth/difference biasanya stasioner. Kita lanjutkan dengan VAR pada \((d\_le, d\_lprod, d\_lrw, d\_U)\).
9.10.3 Pilih Lag Optimum
library(dplyr)
Y <- as.matrix(dplyr::select(can_df, d_le, d_lprod, d_lrw, d_U))
lag_sel <- VARselect(Y, lag.max = 8, type = "const")
lag_sel$selection # Rekomendasi AIC/SC/HQ## AIC(n) HQ(n) SC(n) FPE(n)
## 2 1 1 2lag_sel$criteria # Semua nilai kriteria## 1 2 3 4 5
## AIC(n) -4.416680e+01 -4.419681e+01 -4.395284e+01 -4.397229e+01 -4.373554e+01
## HQ(n) -4.392005e+01 -4.375264e+01 -4.331127e+01 -4.313330e+01 -4.269915e+01
## SC(n) -4.354881e+01 -4.308441e+01 -4.234605e+01 -4.187110e+01 -4.113995e+01
## FPE(n) 6.590898e-20 6.420800e-20 8.273024e-20 8.260333e-20 1.077925e-19
## 6 7 8
## AIC(n) -4.349295e+01 -4.328747e+01 -4.326637e+01
## HQ(n) -4.225915e+01 -4.185626e+01 -4.163775e+01
## SC(n) -4.040297e+01 -3.970309e+01 -3.918759e+01
## FPE(n) 1.436324e-19 1.880081e-19 2.097188e-199.10.4 Estimasi VAR
p_opt <- lag_sel$selection[1] # ambil satu (misal AIC)
var_mod <- VAR(Y, p = p_opt, type = "const")
summary(var_mod)##
## VAR Estimation Results:
## =========================
## Endogenous variables: d_le, d_lprod, d_lrw, d_U
## Deterministic variables: const
## Sample size: 81
## Log Likelihood: 1346.623
## Roots of the characteristic polynomial:
## 0.667 0.667 0.5478 0.5478 0.4377 0.4377 0.3887 0.3448
## Call:
## VAR(y = Y, p = p_opt, type = "const")
##
##
## Estimation results for equation d_le:
## =====================================
## d_le = d_le.l1 + d_lprod.l1 + d_lrw.l1 + d_U.l1 + d_le.l2 + d_lprod.l2 + d_lrw.l2 + d_U.l2 + const
##
## Estimate Std. Error t value Pr(>|t|)
## d_le.l1 9.357e-01 1.522e-01 6.149 3.91e-08 ***
## d_lprod.l1 7.723e-02 2.740e-02 2.818 0.00623 **
## d_lrw.l1 -1.305e-02 2.136e-02 -0.611 0.54327
## d_U.l1 9.867e-05 2.065e-04 0.478 0.63419
## d_le.l2 -3.723e-01 1.641e-01 -2.269 0.02629 *
## d_lprod.l2 9.035e-03 2.809e-02 0.322 0.74864
## d_lrw.l2 -2.026e-02 2.083e-02 -0.973 0.33393
## d_U.l2 -6.821e-05 2.156e-04 -0.316 0.75265
## const 2.260e-04 9.671e-05 2.337 0.02223 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Residual standard error: 0.0003926 on 72 degrees of freedom
## Multiple R-Squared: 0.6497, Adjusted R-squared: 0.6108
## F-statistic: 16.69 on 8 and 72 DF, p-value: 1.041e-13
##
##
## Estimation results for equation d_lprod:
## ========================================
## d_lprod = d_le.l1 + d_lprod.l1 + d_lrw.l1 + d_U.l1 + d_le.l2 + d_lprod.l2 + d_lrw.l2 + d_U.l2 + const
##
## Estimate Std. Error t value Pr(>|t|)
## d_le.l1 -0.4012469 0.6212914 -0.646 0.52044
## d_lprod.l1 0.2124687 0.1118881 1.899 0.06158 .
## d_lrw.l1 0.0544832 0.0872323 0.625 0.53422
## d_U.l1 -0.0022976 0.0008431 -2.725 0.00806 **
## d_le.l2 -1.2057650 0.6700213 -1.800 0.07611 .
## d_lprod.l2 -0.0466083 0.1146841 -0.406 0.68565
## d_lrw.l2 -0.1537768 0.0850367 -1.808 0.07473 .
## d_U.l2 -0.0002858 0.0008803 -0.325 0.74640
## const 0.0012080 0.0003949 3.059 0.00312 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Residual standard error: 0.001603 on 72 degrees of freedom
## Multiple R-Squared: 0.2571, Adjusted R-squared: 0.1746
## F-statistic: 3.115 on 8 and 72 DF, p-value: 0.004458
##
##
## Estimation results for equation d_lrw:
## ======================================
## d_lrw = d_le.l1 + d_lprod.l1 + d_lrw.l1 + d_U.l1 + d_le.l2 + d_lprod.l2 + d_lrw.l2 + d_U.l2 + const
##
## Estimate Std. Error t value Pr(>|t|)
## d_le.l1 -0.0696036 0.7978308 -0.087 0.9307
## d_lprod.l1 -0.1914160 0.1436810 -1.332 0.1870
## d_lrw.l1 0.2577981 0.1120192 2.301 0.0243 *
## d_U.l1 0.0011986 0.0010827 1.107 0.2719
## d_le.l2 1.1734908 0.8604072 1.364 0.1769
## d_lprod.l2 -0.3747133 0.1472714 -2.544 0.0131 *
## d_lrw.l2 0.1619206 0.1091998 1.483 0.1425
## d_U.l2 -0.0002132 0.0011305 -0.189 0.8510
## const 0.0010521 0.0005071 2.075 0.0416 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Residual standard error: 0.002058 on 72 degrees of freedom
## Multiple R-Squared: 0.3933, Adjusted R-squared: 0.3259
## F-statistic: 5.835 on 8 and 72 DF, p-value: 9.696e-06
##
##
## Estimation results for equation d_U:
## ====================================
## d_U = d_le.l1 + d_lprod.l1 + d_lrw.l1 + d_U.l1 + d_le.l2 + d_lprod.l2 + d_lrw.l2 + d_U.l2 + const
##
## Estimate Std. Error t value Pr(>|t|)
## d_le.l1 -562.56339 112.51059 -5.000 3.9e-06 ***
## d_lprod.l1 -61.92467 20.26199 -3.056 0.00314 **
## d_lrw.l1 17.16277 15.79702 1.086 0.28090
## d_U.l1 -0.16069 0.15268 -1.053 0.29608
## d_le.l2 14.94288 121.33516 0.123 0.90233
## d_lprod.l2 -6.39511 20.76831 -0.308 0.75903
## d_lrw.l2 47.22018 15.39942 3.066 0.00305 **
## d_U.l2 -0.25913 0.15942 -1.625 0.10844
## const 0.08461 0.07151 1.183 0.24063
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Residual standard error: 0.2903 on 72 degrees of freedom
## Multiple R-Squared: 0.6063, Adjusted R-squared: 0.5626
## F-statistic: 13.86 on 8 and 72 DF, p-value: 5.711e-12
##
##
##
## Covariance matrix of residuals:
## d_le d_lprod d_lrw d_U
## d_le 1.541e-07 -4.322e-08 -3.850e-09 -7.773e-05
## d_lprod -4.322e-08 2.569e-06 3.696e-07 -1.687e-06
## d_lrw -3.850e-09 3.696e-07 4.237e-06 1.179e-04
## d_U -7.773e-05 -1.687e-06 1.179e-04 8.426e-02
##
## Correlation matrix of residuals:
## d_le d_lprod d_lrw d_U
## d_le 1.000000 -0.068693 -0.004765 -0.682158
## d_lprod -0.068693 1.000000 0.112010 -0.003627
## d_lrw -0.004765 0.112010 1.000000 0.197328
## d_U -0.682158 -0.003627 0.197328 1.0000009.10.5 Diagnostik
# Serial correlation (Portmanteau)
serial.test(var_mod, lags.pt = 12, type = "PT.asymptotic")##
## Portmanteau Test (asymptotic)
##
## data: Residuals of VAR object var_mod
## Chi-squared = 125.43, df = 160, p-value = 0.98# Normalitas residual (Jarque–Bera)
normality.test(var_mod)## $JB
##
## JB-Test (multivariate)
##
## data: Residuals of VAR object var_mod
## Chi-squared = 11.508, df = 8, p-value = 0.1745
##
##
## $Skewness
##
## Skewness only (multivariate)
##
## data: Residuals of VAR object var_mod
## Chi-squared = 8.4598, df = 4, p-value = 0.07611
##
##
## $Kurtosis
##
## Kurtosis only (multivariate)
##
## data: Residuals of VAR object var_mod
## Chi-squared = 3.048, df = 4, p-value = 0.5498# ARCH LM test
arch.test(var_mod, lags.multi = 6)##
## ARCH (multivariate)
##
## data: Residuals of VAR object var_mod
## Chi-squared = 611.88, df = 600, p-value = 0.3595# Stabilitas (akar di luar unit circle)
mod_roots <- roots(var_mod); mod_roots## [1] 0.6670496 0.6670496 0.5478293 0.5478293 0.4376757 0.4376757 0.3886854
## [8] 0.3448158plot(mod_roots, xlab="Real", ylab="Imaginary", main="Roots of the Companion Matrix (Unit Circle)")
abline(h=0, v=0, col="gray70")
9.10.6 Granger Causality
Contoh: apakah d_lprod “Granger‑menyebabkan”
d_le?
causality(var_mod, cause = "d_lprod")## $Granger
##
## Granger causality H0: d_lprod do not Granger-cause d_le d_lrw d_U
##
## data: VAR object var_mod
## F-Test = 3.4444, df1 = 6, df2 = 288, p-value = 0.002662
##
##
## $Instant
##
## H0: No instantaneous causality between: d_lprod and d_le d_lrw d_U
##
## data: VAR object var_mod
## Chi-squared = 2.1809, df = 3, p-value = 0.53579.10.7 IRF (Impulse Responses)
Orde variabel memengaruhi hasil (Cholesky). Kita gunakan urutan (d_le, d_lprod, d_lrw, d_U).
irf_obj <- irf(var_mod, impulse = "d_lprod", response = c("d_le","d_lrw","d_U"),
n.ahead = 12, boot = TRUE, ci = 0.95, runs = 500)
plot(irf_obj)
9.10.8 FEVD
fevd_obj <- fevd(var_mod, n.ahead = 12)
plot(fevd_obj)
9.10.9 Peramalan
fc <- predict(var_mod, n.ahead = 8, ci = 0.95)
# Plot untuk satu seri contoh
plot(fc, names = "d_le")
9.11 Catatan Praktik & Tips
- Skala & transformasi penting: stabilkan varians (log), buat stasioner (diff atau gunakan VECM bila kointegrasi).
- Periksa sensitivitas terhadap pemilihan lag dan urutan variabel (untuk IRF berbasis Cholesky).
- Untuk SVAR (restriksi struktural berbasis teori),
lihat paket
svars/vars::SVAR. - Jangan lupa cek overfitting: VAR mudah “lapar parameter” (\(K^2p\) koefisien + konstanta).
9.12 Ringkasan Satu
VAR memodelkan interaksi dinamis multi‑variabel secara fleksibel. Estimasi mudah (OLS per persamaan), evaluasi pakai uji serial correlation/normalitas/ARCH, stabilitas via akar. Kegunaan utama: IRF, FEVD, peramalan. Saat ada kointegrasi, beralih ke VECM dan/atau SVAR untuk interpretasi struktural.
10 Ringkasan Model
| Model | Fokus Analisis | Asumsi Kunci | Kegunaan Utama |
|---|---|---|---|
| OLS | Hubungan linier rata-rata | Eksogenitas, homoskedastisitas | Cross‑section & baseline |
| IV | Kausalitas dgn endogenitas | Instrumen relevan & valid | Causal inference |
| FE (Panel) | Kontrol heterogenitas tetap | Efek tetap konstan per unit | Panel data |
| DiD | Evaluasi kebijakan | Parallel trends | Eksperimen alami |
| Synthetic Control | Unit sintetik pembanding | Good pre‑treatment fit | Studi kebijakan per-unit |
| VAR | Dinamika multivariat | (Kira‑)stasioneritas | Makro/keuangan/DERET waktu |
Untuk komputasi dapat dicek di link: https://project25.shinyapps.io/ekonometrika/
11 Referensi Singkat
11.1 Dasar Ekonometrika dan Panel Data
Wooldridge, J. M. (2010). Econometric Analysis of Cross Section and Panel Data (2nd ed.). MIT Press.
Angrist, J. D., & Pischke, J.-S. (2009). Mostly Harmless Econometrics: An Empiricist’s Companion. Princeton University Press.
11.2 Metode Synthetic Control (SCM)
Abadie, A., Diamond, A., & Hainmueller, J. (2010). Synthetic control methods for comparative case studies: Estimating the effect of California’s tobacco control program. Journal of the American Statistical Association, 105(490), 493–505. https://doi.org/10.1198/jasa.2009.ap08746
Abadie, A. (2021). Using Synthetic Controls: Methods and Practice. Cambridge University Press.
Ben-Michael, E., Feller, A., & Rothstein, J. (2021). The augmented synthetic control method. Journal of the American Statistical Association, 116(536), 1789–1803. https://doi.org/10.1080/01621459.2021.1929245
11.3 Perluasan SCM dan Panel Dinamis (DiD & Event Study)
Imai, K., & Kim, I. S. (2021). On the use of two-way fixed effects models for causal inference with heterogeneous treatment effects. American Journal of Political Science, 65(3), 849–866. https://doi.org/10.1111/ajps.12674
Callaway, B., & Sant’Anna, P. H. C. (2021). Difference-in-differences with multiple time periods. Journal of Econometrics, 225(2), 200–230. https://doi.org/10.1016/j.jeconom.2020.12.001
Sun, L., & Abraham, S. (2021). Estimating dynamic treatment effects in event studies with heterogeneous treatment effects. Journal of Econometrics, 225(2), 175–199. https://doi.org/10.1016/j.jeconom.2020.09.006
Lütkepohl, H. (2005). New Introduction to Multiple Time Series Analysis. Springer.
11.4 VAR
- Lütkepohl, H. (2005). New Introduction to Multiple Time Series
Analysis. Springer.
- Kilian, L., & Lütkepohl, H. (2017). Structural Vector
Autoregressive Analysis. Cambridge University Press.
- Pfaff, B. (2008). VAR, SVAR and SVEC Models: Implementation
Within R Package vars. Journal of Statistical Software,
27(4).
- Stock, J. H., & Watson, M. W. (2001). Vector Autoregressions. Journal of Economic Perspectives, 15(4), 101–115.