1 Pendahuluan

Praktikum ini kita akan mempelajari tentang menyusun matriks bobot spasial sampai menduga parameter dengan tiga cara: OLS, Instrumental Variables (IV), dan Maximum Likelihood (ML). Datanya hanya lima wilayah, jadi setiap angka bisa dicocokkan langsung dengan slide kuliah. Tujuannya memahami prosedurnya, bukan menarik kesimpulan substantif dari lima titik data.

Setelah praktikum ini, diharapkan bisa membuat matriks bobot spasial di R, menghitung spatial lag Wy, menduga model SAR dengan ketiga metode tersebut, dan membaca output lagsarlm().

Paket yang dipakai:

library(ggplot2)     # grafik
library(spdep)       # objek ketetanggaan (nb, listw)
library(spatialreg)  # lagsarlm(), stsls()
library(knitr)       # tabel

biru   <- "#1b6ca8"
oranye <- "#e07a1f"
theme_set(theme_minimal(base_size = 13) +
            theme(plot.title = element_text(face = "bold"),
                  plot.subtitle = element_text(color = "grey35"),
                  panel.grid.minor = element_blank()))

Jika ada paket yang belum terpasang, jalankan install.packages(c("ggplot2", "spdep", "spatialreg", "knitr")) terlebih dahulu.

2 Data dan Matriks Bobot Spasial

Kasusnya begini: kita ingin menjelaskan pertumbuhan produktivitas (Y) di lima wilayah. Penjelas yang paling natural adalah pertumbuhan output wilayah itu sendiri (X1). Namun produktivitas suatu wilayah mungkin juga dipengaruhi produktivitas wilayah di sekitarnya, misalnya lewat perpindahan tenaga kerja, teknologi, atau perdagangan. Pengaruh seperti ini ditangkap model SAR dengan menambahkan spatial lag pada variabel dependen:

\[\mathbf{y} = \rho \mathbf{W}\mathbf{y} + \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon}, \qquad \boldsymbol{\varepsilon} \sim N(\mathbf{0}, \sigma^2 \mathbf{I}_n)\]

Untuk wilayah ke-\(i\), persamaannya menjadi \(y_i = \beta_1 + \beta_2 x_i + \rho \sum_{j} w_{ij} y_j + \varepsilon_i\). Matriks \(\mathbf{W}\) memberi tahu R siapa bertetangga dengan siapa, dan \(\rho\) mengukur seberapa kuat nilai \(y\) di tetangga terbawa ke wilayah \(i\). Bila \(\mathbf{W}\) distandarisasi per baris, \((\mathbf{Wy})_i\) tidak lain adalah rata-rata \(y\) di wilayah-wilayah tetangga.

Data dan koordinat kelima wilayah kita buat langsung di R.

Data <- data.frame(
  Lokasi = LETTERS[1:5],
  X1     = c(0.6, 1.0, 1.6, 2.6, 2.2),   # pertumbuhan output
  Y      = c(0.4, 0.6, 0.9, 1.1, 1.2)    # pertumbuhan produktivitas
)

Koordinat <- data.frame(
  Lokasi   = LETTERS[1:5],
  Datar    = c(1, 3, 3, 5, 7),
  Vertikal = c(3, 1, 3, 3, 3)
)
Data
##   Lokasi  X1   Y
## 1      A 0.6 0.4
## 2      B 1.0 0.6
## 3      C 1.6 0.9
## 4      D 2.6 1.1
## 5      E 2.2 1.2

(Di kelas data dibaca dari Excel dengan readxl::read_excel(). Di sini data dibuat langsung supaya file ini bisa dijalankan tanpa mengatur path folder.)

Matriks ketetanggaan kita susun dari matriks 0/1, angka 1 berarti dua wilayah bertetangga. Setelah itu tiap baris dibagi jumlah barisnya (standarisasi baris) sehingga setiap baris berjumlah 1.

A <- matrix(c(0, 1, 1, 0, 0,     # A bertetangga dengan B, C
              1, 0, 1, 1, 0,     # B bertetangga dengan A, C, D
              1, 1, 0, 1, 0,     # C bertetangga dengan A, B, D
              0, 1, 1, 0, 1,     # D bertetangga dengan B, C, E
              0, 0, 0, 1, 0),    # E hanya bertetangga dengan D
            nrow = 5, byrow = TRUE,
            dimnames = list(Data$Lokasi, Data$Lokasi))

W <- A / rowSums(A)
round(W, 3)
##       A     B     C     D     E
## A 0.000 0.500 0.500 0.000 0.000
## B 0.333 0.000 0.333 0.333 0.000
## C 0.333 0.333 0.000 0.333 0.000
## D 0.000 0.333 0.333 0.000 0.333
## E 0.000 0.000 0.000 1.000 0.000
rowSums(W)
## A B C D E 
## 1 1 1 1 1

Hasilnya sama dengan matriks \(\mathbf{W}\) di slide: baris A berisi 1/2 untuk B dan C, baris E berisi 1 untuk D, dan seterusnya. Pembobotan lain yang umum dipakai adalah invers jarak, yaitu \(w_{ij} = 1/d_{ij}\), sehingga wilayah yang lebih dekat mendapat bobot lebih besar.

jarak <- as.matrix(dist(Koordinat[, c("Datar", "Vertikal")]))
dimnames(jarak) <- list(Data$Lokasi, Data$Lokasi)

W_inv <- 1 / jarak
diag(W_inv) <- 0                      # wilayah tidak bertetangga dengan dirinya sendiri
W_inv_std <- W_inv / rowSums(W_inv)   # versi terstandarisasi baris

round(W_inv, 4)
##        A      B    C      D      E
## A 0.0000 0.3536 0.50 0.2500 0.1667
## B 0.3536 0.0000 0.50 0.3536 0.2236
## C 0.5000 0.5000 0.00 0.5000 0.2500
## D 0.2500 0.3536 0.50 0.0000 0.5000
## E 0.1667 0.2236 0.25 0.5000 0.0000
round(W_inv_std, 4)
##        A      B      C      D      E
## A 0.0000 0.2783 0.3936 0.1968 0.1312
## B 0.2471 0.0000 0.3495 0.2471 0.1563
## C 0.2857 0.2857 0.0000 0.2857 0.1429
## D 0.1559 0.2205 0.3118 0.0000 0.3118
## E 0.1462 0.1961 0.2192 0.4385 0.0000

Dengan \(\mathbf{W}\) di tangan, spatial lag diperoleh dengan mengalikannya pada vektor Y. Kita simpan juga \(\mathbf{WX}\), yang nanti dibutuhkan sebagai instrumen.

Data$Wy_tetangga <- as.numeric(W %*% Data$Y)
Data$Wy_inv      <- as.numeric(W_inv_std %*% Data$Y)
Data$WX          <- as.numeric(W %*% Data$X1)
Data
##   Lokasi  X1   Y Wy_tetangga    Wy_inv  WX
## 1      A 0.6 0.4        0.75 0.8952244 1.3
## 2      B 1.0 0.6        0.80 0.8727521 1.6
## 3      C 1.6 0.9        0.70 0.7714286 1.4
## 4      D 2.6 1.1        0.90 0.8494460 1.6
## 5      E 2.2 1.2        1.10 0.8557866 2.6

Kolom Wy_tetangga harus berisi 0,75; 0,80; 0,70; 0,90; dan 1,10, sama seperti di slide. Sebagai contoh, wilayah A punya dua tetangga, B (0,6) dan C (0,9), sehingga rata-ratanya 0,75.

Supaya lebih terbayang, kita gambar dulu strukturnya. Garis menghubungkan wilayah yang bertetangga dan warna lingkaran menunjukkan nilai Y.

garis <- which(A == 1 & upper.tri(A), arr.ind = TRUE)
garis <- data.frame(x    = Koordinat$Datar[garis[, 1]],    y    = Koordinat$Vertikal[garis[, 1]],
                    xend = Koordinat$Datar[garis[, 2]],    yend = Koordinat$Vertikal[garis[, 2]])

peta <- merge(Koordinat, Data, by = "Lokasi")
peta$teks <- ifelse(peta$Y > 0.8, "white", "grey15")

ggplot() +
  geom_segment(data = garis, aes(x = x, y = y, xend = xend, yend = yend),
               color = "grey60", linewidth = 0.9) +
  geom_point(data = peta, aes(Datar, Vertikal, fill = Y),
             shape = 21, size = 14, color = "white", stroke = 1.5) +
  geom_text(data = peta, aes(Datar, Vertikal, label = Lokasi, color = teks),
            fontface = "bold", size = 5) +
  scale_color_identity() +
  scale_fill_gradient(low = "#c6dbef", high = "#08306b", name = "Produktivitas (Y)") +
  coord_equal(xlim = c(0, 8), ylim = c(0, 4)) +
  labs(title = "Peta Ketetanggaan Lima Wilayah",
       subtitle = "Garis menunjukkan wilayah yang bertetangga",
       x = NULL, y = NULL) +
  theme(axis.text = element_blank(), panel.grid = element_blank())

Wilayah D dan E, yang produktivitasnya paling tinggi, saling bertetangga, sementara A dan B yang rendah juga berdekatan. Kesan ini bisa dilihat lebih jelas dengan memplotkan Y terhadap rata-rata Y tetangganya.

ggplot(Data, aes(Y, Wy_tetangga)) +
  geom_vline(xintercept = mean(Data$Y), linetype = "dashed", color = "grey60") +
  geom_hline(yintercept = mean(Data$Wy_tetangga), linetype = "dashed", color = "grey60") +
  geom_smooth(method = "lm", se = FALSE, color = oranye, linewidth = 1) +
  geom_point(size = 4.5, color = biru) +
  geom_text(aes(label = Lokasi), nudge_y = 0.035, fontface = "bold") +
  labs(title = "Y terhadap Spatial Lag Wy (W ketetanggaan)",
       subtitle = "Garis putus-putus menunjukkan rata-rata masing-masing sumbu",
       x = "Produktivitas wilayah (Y)", y = "Rata-rata produktivitas tetangga (Wy)")

Garisnya naik: wilayah dengan Y tinggi cenderung dikelilingi tetangga dengan Y tinggi, dan sebaliknya. Pola mengelompok seperti ini disebut autokorelasi spasial positif, dan inilah alasan kita memasukkan \(\rho\mathbf{Wy}\) ke dalam model. Perlu diingat bahwa dengan \(n = 5\) plot ini hanya ilustrasi. Pada data sungguhan, keberadaan autokorelasi spasial diuji dulu dengan Moran’s I dan uji Lagrange Multiplier sebelum model spasial dipilih.

3 SAR dengan OLS

Cara yang paling mudah untuk menduga model SAR adalah memperlakukan Wy sebagai variabel penjelas biasa lalu menjalankan lm(). Cara ini keliru, dan kekeliruannya bisa ditunjukkan secara aljabar. Dugaan OLS untuk \(\rho\) adalah \(\hat{\rho} = [(\mathbf{Wy})'\mathbf{Wy}]^{-1}(\mathbf{Wy})'\mathbf{y}\). Dengan mensubstitusi \(\mathbf{y} = \rho \mathbf{Wy} + \boldsymbol{\varepsilon}\) diperoleh

\[\hat{\rho} = \rho + [(\mathbf{Wy})'\mathbf{Wy}]^{-1}(\mathbf{Wy})'\boldsymbol{\varepsilon}.\]

Suku kedua tidak akan menghilang. Karena \(\mathbf{y} = (\mathbf{I} - \rho\mathbf{W})^{-1}(\mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon})\), kita punya \((\mathbf{Wy})'\boldsymbol{\varepsilon} = \boldsymbol{\varepsilon}'\mathbf{W}(\mathbf{I}-\rho\mathbf{W})^{-1}\boldsymbol{\varepsilon}\), yaitu bentuk kuadratik dalam \(\boldsymbol{\varepsilon}\) yang tidak bernilai nol selama \(\rho \neq 0\). Maka \(\operatorname{plim}\hat{\rho} \neq \rho\). Dengan kata lain \(\mathbf{Wy}\) berkorelasi dengan galat (endogen), sehingga penduga OLS bias dan tidak konsisten (Anselin, 1988). Bias ini tidak hilang walaupun datanya diperbanyak.

Meski begitu, kita tetap menjalankannya untuk melihat seperti apa hasilnya. Sebagai pembanding disertakan juga regresi biasa tanpa unsur spasial.

ols_biasa    <- lm(Y ~ X1, data = Data)
ols_tetangga <- lm(Y ~ X1 + Wy_tetangga, data = Data)
ols_inv      <- lm(Y ~ X1 + Wy_inv, data = Data)

kable(round(summary(ols_biasa)$coefficients, 4),
      caption = "Regresi biasa: Y ~ X1")
Regresi biasa: Y ~ X1
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.2165 0.1217 1.7793 0.1732
X1 0.3897 0.0691 5.6434 0.0110 
kable(round(summary(ols_tetangga)$coefficients, 4),
      caption = "OLS dengan Wy, W ketetanggaan terstandar")
OLS dengan Wy, W ketetanggaan terstandar
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.0176 0.3606 -0.0487 0.9656
X1 0.3464 0.0979 3.5378 0.0714
Wy_tetangga 0.3568 0.5107 0.6987 0.5570 
kable(round(summary(ols_inv)$coefficients, 4),
      caption = "OLS dengan Wy, W invers jarak terstandar")
OLS dengan Wy, W invers jarak terstandar
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.9494 1.3183 0.7201 0.5462
X1 0.3726 0.0844 4.4153 0.0477
Wy_inv -0.8311 1.4867 -0.5590 0.6324

Angka-angka ini sama dengan slide: untuk W ketetanggaan \(\hat\rho = 0{,}3568\), \(\hat\beta_0 = -0{,}0176\), dan \(\hat\beta_1 = 0{,}3464\); untuk W invers jarak \(\hat\rho = -0{,}8311\), \(\hat\beta_0 = 0{,}9494\), dan \(\hat\beta_1 = 0{,}3726\). Perhatikan bahwa tanda \(\hat\rho\) berbalik hanya karena \(\mathbf{W}\) diganti, dan galat bakunya besar (nilai-\(p\) untuk \(\hat\rho\) masing-masing 0,557 dan 0,632). Dugaan seperti ini tidak stabil.

Bias OLS tidak mudah terlihat dari lima data, jadi kita buat dunia buatan yang nilai \(\rho\) sebenarnya kita ketahui, yaitu \(\rho = 0{,}6\), pada 100 wilayah berbentuk kisi 10×10. Data dibangkitkan 300 kali dari model SAR, dan setiap kali \(\rho\) diduga dengan OLS.

set.seed(2026)
W_sim  <- nb2mat(cell2nb(10, 10), style = "W")
n_sim  <- nrow(W_sim)
rho_benar  <- 0.6
invers_IrW <- solve(diag(n_sim) - rho_benar * W_sim)

rho_ols <- replicate(300, {
  x   <- rnorm(n_sim)
  eps <- rnorm(n_sim)
  y   <- as.numeric(invers_IrW %*% (1 + 2 * x + eps))
  Wy  <- as.numeric(W_sim %*% y)
  coef(lm(y ~ x + Wy))["Wy"]
})

ggplot(data.frame(rho = rho_ols), aes(rho)) +
  geom_histogram(bins = 25, fill = biru, color = "white", alpha = 0.9) +
  geom_vline(xintercept = rho_benar, color = "red3", linewidth = 1.1) +
  geom_vline(xintercept = mean(rho_ols), color = "grey15", linewidth = 1.1, linetype = "dashed") +
  annotate("text", x = rho_benar, y = Inf, label = "ρ sebenarnya = 0,6",
           hjust = 1.1, vjust = 2, color = "red3", fontface = "bold") +
  annotate("text", x = mean(rho_ols), y = Inf, label = "rata-rata dugaan OLS",
           hjust = -0.1, vjust = 4, color = "grey15", fontface = "bold") +
  labs(title = "Sebaran Dugaan ρ dengan OLS (300 simulasi)",
       subtitle = "Penduga yang tidak bias akan berpusat tepat di garis merah",
       x = expression(hat(rho)[OLS]), y = "Frekuensi")

c(rho_sebenarnya = rho_benar, rata_rata_dugaan_OLS = mean(rho_ols))
##       rho_sebenarnya rata_rata_dugaan_OLS 
##             0.600000             0.673505

Jika OLS tidak bias, histogram akan berpusat di garis merah. Jarak antara garis putus-putus (rata-rata dugaan) dan garis merah itulah bias yang dibicarakan di atas. Karena itu uji-\(t\) dan uji-\(F\) dari lm() tidak valid untuk model SAR, dan kita membutuhkan penduga yang konsisten. Dua pilihannya adalah IV dan ML.

4 SAR dengan Instrumental Variables

Sumber masalah OLS adalah Wy yang berkorelasi dengan galat. Gagasan IV adalah menggantinya dengan versi yang bersih dari galat, yaitu hasil prediksi Wy dari variabel yang eksogen: \(\mathbf{X}\) dan \(\mathbf{WX}\). Prosesnya dua tahap (sehingga sering disebut Two-Stage Least Squares): pertama memprediksi Wy, kedua meregresikan \(\mathbf{y}\) pada hasil prediksi tersebut.

Secara matriks, definisikan \(\mathbf{Z} = [\mathbf{X}\;\; \mathbf{Wy}]\) dan \(\boldsymbol{\theta} = [\boldsymbol{\beta}'\;\; \rho]'\) sehingga model menjadi \(\mathbf{y} = \mathbf{Z}\boldsymbol{\theta} + \boldsymbol{\varepsilon}\). Dengan matriks instrumen \(\mathbf{Q} = [\mathbf{X}\;\; \mathbf{WX}]\),

\[\hat{\mathbf{Z}} = \mathbf{Q}(\mathbf{Q}'\mathbf{Q})^{-1}\mathbf{Q}'\mathbf{Z}, \qquad \hat{\boldsymbol{\theta}}_{IV} = (\hat{\mathbf{Z}}'\mathbf{Z})^{-1}\hat{\mathbf{Z}}'\mathbf{y}.\]

Matriks kovariansnya adalah \(\widehat{\operatorname{Cov}}(\hat{\boldsymbol\theta}) = \hat\sigma^2(\hat{\mathbf{Z}}'\hat{\mathbf{Z}})^{-1}\) dengan \(\hat\sigma^2 = \mathbf{e}'\mathbf{e}/n\) dan \(\mathbf{e} = \mathbf{y} - \mathbf{Z}\hat{\boldsymbol\theta}\). Rumus-rumus ini kita terjemahkan langsung ke R.

y <- Data$Y
n <- length(y)

Z <- cbind("(Intercept)" = 1, X1 = Data$X1, Wy = Data$Wy_tetangga)
Q <- cbind(1, Data$X1, Data$WX)

# tahap 1: Wy diganti hasil prediksinya
Zhat <- Q %*% solve(t(Q) %*% Q) %*% t(Q) %*% Z
round(Zhat, 4)
##      (Intercept)  X1     Wy
## [1,]           1 0.6 0.7109
## [2,]           1 1.0 0.8025
## [3,]           1 1.6 0.7730
## [4,]           1 2.6 0.8607
## [5,]           1 2.2 1.1029
# tahap 2: dugaan theta
theta <- solve(t(Zhat) %*% Z) %*% t(Zhat) %*% y
round(theta, 4)
##                [,1]
## (Intercept) -0.2057
## X1           0.3116
## Wy           0.6437

Hasilnya \(\hat\beta_0 = -0{,}2057\), \(\hat\beta_1 = 0{,}3116\), dan \(\hat\rho = 0{,}6437\), persis seperti di slide. Kolom ketiga Zhat (0,7109; 0,8025; 0,7730; 0,8607; 1,1029) adalah Wy yang sudah “dibersihkan”. Perubahannya dari nilai aslinya bisa digambar sebagai berikut.

df_iv <- data.frame(Lokasi = Data$Lokasi,
                    asli   = Data$Wy_tetangga,
                    murni  = Zhat[, 3])

ggplot(df_iv, aes(asli, murni)) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "grey50") +
  geom_segment(aes(xend = asli, yend = asli), color = oranye, linewidth = 1) +
  geom_point(size = 4.5, color = biru) +
  geom_text(aes(label = Lokasi), nudge_y = 0.03, fontface = "bold") +
  labs(title = "Wy sebelum dan sesudah Diprediksi dari X dan WX",
       subtitle = "Garis oranye menunjukkan pergeseran dari Wy asli ke hasil prediksinya",
       x = "Wy asli", y = "Wy hasil prediksi")

Titik yang tepat di garis diagonal berarti nilainya tidak berubah; makin panjang garis oranye, makin jauh Wy bergeser. Selanjutnya kita hitung ragam galat, galat baku, dan uji-\(z\) untuk setiap koefisien.

e         <- y - Z %*% theta
sigma2    <- sum(e^2) / n
cov_theta <- sigma2 * solve(t(Zhat) %*% Zhat)

se <- sqrt(diag(cov_theta))
z  <- as.numeric(theta) / se
p  <- 2 * (1 - pnorm(abs(z)))

hasil_iv <- data.frame(
  Koefisien = rownames(theta), Estimasi = as.numeric(theta),
  SE = se, z = z, p_value = p,
  Keputusan = ifelse(abs(z) > qnorm(0.975), "Tolak H0", "Gagal tolak H0")
)
kable(hasil_iv, digits = 4, row.names = FALSE,
      caption = paste0("Hasil IV (sigma^2 = ", round(sigma2, 4),
                       ", sigma = ", round(sqrt(sigma2), 4), ")"))
Hasil IV (sigma^2 = 0.0072, sigma = 0.0851)
Koefisien Estimasi SE z p_value Keputusan
(Intercept) -0.2057 0.2621 -0.7849 0.4325 Gagal tolak H0
X1 0.3116 0.0688 4.5306 0.0000 Tolak H0
Wy 0.6437 0.3748 1.7174 0.0859 Gagal tolak H0

Sebagai pemeriksaan, hasil manual ini bisa dibandingkan dengan fungsi bawaan stsls(). Argumen W2X = FALSE dipakai supaya instrumennya sama persis dengan yang kita susun, yaitu \([\mathbf{X}\;\mathbf{WX}]\). Koefisiennya harus sama, sedangkan galat bakunya boleh sedikit berbeda karena penyebut ragam galat yang dipakai bisa berbeda.

iv_paket <- stsls(Y ~ X1, data = Data, listw = mat2listw(W, style = "W"), W2X = FALSE)
round(coef(iv_paket), 4)    # urutan: Rho, (Intercept), X1
##         Rho (Intercept)          X1 
##      0.6437     -0.2057      0.3116

Dengan nilai kritis \(z_{0,975} = 1{,}96\), intersep tidak signifikan (\(|z| = 0{,}785\)), koefisien X1 signifikan (\(|z| = 4{,}530\)), dan \(\rho\) belum signifikan (\(|z| = 1{,}717\)). Persamaan dugaannya adalah \(\hat y_i = -0{,}206 + 0{,}312\,x_i + 0{,}644\sum_j w_{ij}y_j\) dengan \(\hat\sigma^2 = 0{,}0072\). Kenaikan output sebesar satu satuan berasosiasi dengan kenaikan produktivitas sekitar 0,31 satuan, dengan efek tetangga bertanda positif tetapi belum bisa dibuktikan secara statistik. Nilai \(\hat\rho_{IV} = 0{,}644\) juga jauh berbeda dari \(\hat\rho_{OLS} = 0{,}357\), dan selisih itulah koreksi terhadap bias OLS.

5 SAR dengan Maximum Likelihood

Metode ML bertanya nilai \(\rho\), \(\boldsymbol\beta\), dan \(\sigma^2\) mana yang membuat data yang teramati paling mungkin muncul. Dengan \(\boldsymbol\varepsilon = (\mathbf{I} - \rho\mathbf{W})\mathbf{y} - \mathbf{X}\boldsymbol\beta\) dan Jacobian \(J = |\mathbf{I} - \rho\mathbf{W}|\), fungsi log-likelihood-nya adalah

\[\ln L(\boldsymbol\beta, \rho, \sigma^2 \mid \mathbf{y}, \mathbf{X}) = -\frac{n}{2}\ln(2\pi) - \frac{n}{2}\ln\sigma^2 + \ln|\mathbf{I} - \rho\mathbf{W}| - \frac{1}{2\sigma^2}\boldsymbol\varepsilon'\boldsymbol\varepsilon.\]

Dibandingkan dengan regresi biasa, ada tambahan suku \(\ln|\mathbf{I} - \rho\mathbf{W}|\). Suku inilah yang memperhitungkan bahwa semua wilayah saling memengaruhi secara simultan, dan karenanya ML tidak jatuh ke masalah yang dialami OLS. Memaksimumkan \(\ln L\) sama dengan meminimumkan \(\mathbf{e}'\mathbf{e}\) yang dikoreksi log-Jacobian; hasilnya sama dengan OLS hanya jika \(\rho \to 0\). Sifat konsistensi dan efisiensi asimtotik berlaku bila kondisi keteraturan terpenuhi, khususnya \(J > 0\).

Karena \(\boldsymbol\beta\) dan \(\sigma^2\) dapat dinyatakan sebagai fungsi dari \(\rho\), pemaksimuman dilakukan dalam empat langkah pada concentrated likelihood. Mari kita kerjakan satu per satu. Langkah pertama dan kedua adalah meregresikan \(\mathbf{y}\) dan \(\mathbf{Wy}\) masing-masing pada \(\mathbf{X}\) dan menyimpan residunya.

X    <- cbind(1, Data$X1)
XtXi <- solve(t(X) %*% X)

# langkah 1: regresi y pada X
b_O <- XtXi %*% t(X) %*% y
e_O <- y - X %*% b_O

# langkah 2: regresi Wy pada X
Wy  <- W %*% y
b_L <- XtXi %*% t(X) %*% Wy
e_L <- Wy - X %*% b_L

list(b_O = round(as.numeric(b_O), 4), b_L = round(as.numeric(b_L), 4))
## $b_O
## [1] 0.2165 0.3897
## 
## $b_L
## [1] 0.6559 0.1213

Nilainya harus sama dengan slide: \(\hat\beta_O = (0{,}2165;\ 0{,}3897)\) dan \(\hat\beta_L = (0{,}6562;\ 0{,}1211)\). Langkah ketiga mencari \(\rho\) yang memaksimumkan \(\ln L_c(\rho) = C - \frac{n}{2}\ln\!\big[\tfrac{1}{n}(\mathbf{e}_O - \rho\mathbf{e}_L)'(\mathbf{e}_O - \rho\mathbf{e}_L)\big] + \ln|\mathbf{I} - \rho\mathbf{W}|\), dan langkah keempat menghitung \(\hat{\boldsymbol\beta}_{ML} = \hat{\boldsymbol\beta}_O - \hat\rho_{ML}\hat{\boldsymbol\beta}_L\) serta \(\hat\sigma^2_{ML}\).

# langkah 3: maksimumkan concentrated log-likelihood terhadap rho
lnLc <- function(rho) {
  e <- e_O - rho * e_L
  -n / 2 * log(sum(e^2) / n) +
    as.numeric(determinant(diag(n) - rho * W, logarithm = TRUE)$modulus)
}

opt    <- optimize(lnLc, interval = c(-0.99, 0.99), maximum = TRUE)
rho_ml <- opt$maximum

# langkah 4: dugaan beta dan sigma^2
b_ML      <- b_O - rho_ml * b_L
e_ML      <- e_O - rho_ml * e_L
sigma2_ml <- sum(e_ML^2) / n

data.frame(rho = rho_ml, beta0 = b_ML[1], beta1 = b_ML[2],
           sigma2 = sigma2_ml, lnLc = opt$objective) |>
  kable(digits = 4, caption = "Hasil ML manual, W ketetanggaan")
Hasil ML manual, W ketetanggaan
rho beta0 beta1 sigma2 lnLc
0.2872 0.0281 0.3549 0.0063 12.573

Hasilnya \(\hat\rho_{ML} = 0{,}2872\), \(\hat{\boldsymbol\beta}_{ML} = (0{,}0280;\ 0{,}3549)\), \(\hat\sigma^2_{ML} = 0{,}0063\), dan \(\ln L_c^* = 12{,}5729\) (tanpa konstanta \(C\)), sama dengan slide. Proses pencarian \(\hat\rho\) tadi bisa dilihat langsung dengan menghitung \(\ln L_c(\rho)\) pada banyak nilai \(\rho\). Kita lakukan untuk kedua jenis \(\mathbf{W}\).

profil_lnLc <- function(Wm, label) {
  XtXi <- solve(t(X) %*% X)
  e_O  <- y - X %*% (XtXi %*% t(X) %*% y)
  Wy   <- Wm %*% y
  e_L  <- Wy - X %*% (XtXi %*% t(X) %*% Wy)

  f <- function(rho) {
    e <- e_O - rho * e_L
    -n / 2 * log(sum(e^2) / n) +
      as.numeric(determinant(diag(n) - rho * Wm, logarithm = TRUE)$modulus)
  }
  ev    <- Re(eigen(Wm)$values)
  batas <- c(1 / min(ev), 1 / max(ev))          # rentang rho dengan J > 0
  grid  <- seq(batas[1] + 0.01, batas[2] - 0.01, length.out = 400)
  opt   <- optimize(f, interval = range(grid), maximum = TRUE)

  list(kurva = data.frame(rho = grid, lnLc = sapply(grid, f), W = label),
       puncak = data.frame(rho = opt$maximum, lnLc = opt$objective, W = label))
}

p1 <- profil_lnLc(W,     "W ketetanggaan (terstandar)")
p2 <- profil_lnLc(W_inv, "W invers jarak")

kurva  <- rbind(p1$kurva,  p2$kurva)
puncak <- rbind(p1$puncak, p2$puncak)

ggplot(kurva, aes(rho, lnLc)) +
  geom_line(color = biru, linewidth = 1.1) +
  geom_vline(data = puncak, aes(xintercept = rho), color = oranye, linetype = "dashed") +
  geom_point(data = puncak, color = oranye, size = 4) +
  geom_text(data = puncak, aes(label = paste0("ρ-topi = ", round(rho, 3))),
            color = oranye, fontface = "bold", vjust = -1.1, hjust = 0.5) +
  geom_vline(xintercept = 0, color = "grey60") +
  facet_wrap(~ W, scales = "free") +
  labs(title = "Concentrated Log-Likelihood terhadap ρ",
       subtitle = "Titik oranye adalah nilai ρ yang memaksimumkan likelihood",
       x = expression(rho), y = expression(ln~L[c](rho))) +
  theme(strip.text = element_text(face = "bold"))

Puncak kurva adalah dugaan ML untuk \(\rho\). Garis abu-abu di \(\rho = 0\) menandai kondisi tanpa efek spasial: makin jauh dan makin tajam puncaknya dari garis itu, makin kuat indikasi bahwa \(\rho \neq 0\). Dalam praktik, pencarian ini tidak perlu dikerjakan manual. Fungsi lagsarlm() dari paket spatialreg melakukan semuanya sekaligus.

lw_tetangga <- mat2listw(unname(W), style = "W")

Hasil <- lagsarlm(Y ~ X1, data = Data, listw = lw_tetangga)
summary(Hasil)
## 
## Call:lagsarlm(formula = Y ~ X1, data = Data, listw = lw_tetangga)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.109224 -0.056418 -0.012722  0.075287  0.103076 
## 
## Type: lag 
## Coefficients: (asymptotic standard errors) 
##             Estimate Std. Error z value  Pr(>|z|)
## (Intercept) 0.028118   0.164543  0.1709    0.8643
## X1          0.354865   0.056629  6.2664 3.694e-10
## 
## Rho: 0.28717, LR test value: 0.86648, p-value: 0.35193
## Asymptotic standard error: 0.21686
##     z-value: 1.3242, p-value: 0.18543
## Wald statistic: 1.7536, p-value: 0.18543
## 
## Log likelihood: 5.478291 for lag model
## ML residual variance (sigma squared): 0.0063135, (sigma: 0.079458)
## Number of observations: 5 
## Number of parameters estimated: 4 
## AIC: -2.9566, (AIC for lm: -4.0901)
## LM test for residual autocorrelation
## test value: 3.9562, p-value: 0.046698

Cara membaca outputnya: blok Coefficients berisi intersep dan koefisien X1 beserta uji-\(z\); baris Rho berisi \(\hat\rho\) yang diuji dengan uji LR dan uji Wald; sedangkan AIC dibandingkan dengan AIC for lm, yaitu AIC model tanpa efek spasial (makin kecil makin baik). Angkanya sama dengan hitungan manual: \(\hat\rho = 0{,}2872\), intersep \(0{,}0281\), dan \(\hat\beta_1 = 0{,}3549\).

Pada W ketetanggaan ini, koefisien X1 sangat signifikan (\(z = 6{,}27\); \(p < 0{,}001\)), tetapi \(\rho\) tidak: uji LR memberi \(p = 0{,}352\) dan uji Wald \(p = 0{,}185\), sehingga \(H_0: \rho = 0\) tidak ditolak. AIC model lag (\(-2{,}957\)) bahkan lebih besar daripada AIC model biasa (\(-4{,}090\)), artinya untuk data ini model SAR tidak lebih baik daripada regresi klasik.

Sekarang kita ulangi dengan W invers jarak. Mengikuti slide, matriks ini tidak distandarisasi per baris.

Hasil2 <- lagsarlm(Y ~ X1, data = Data, listw = mat2listw(unname(W_inv)))
summary(Hasil2)
## 
## Call:lagsarlm(formula = Y ~ X1, data = Data, listw = mat2listw(unname(W_inv)))
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -0.0890398 -0.0682359  0.0092734  0.0439466  0.1040557 
## 
## Type: lag 
## Coefficients: (asymptotic standard errors) 
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.600413   0.245876  2.4419  0.01461
## X1          0.397047   0.043479  9.1320  < 2e-16
## 
## Rho: -0.32574, LR test value: 1.8613, p-value: 0.17248
## Asymptotic standard error: 0.19886
##     z-value: -1.638, p-value: 0.10142
## Wald statistic: 2.6831, p-value: 0.10142
## 
## Log likelihood: 5.975688 for lag model
## ML residual variance (sigma squared): 0.0050858, (sigma: 0.071315)
## Number of observations: 5 
## Number of parameters estimated: 4 
## AIC: -3.9514, (AIC for lm: -4.0901)
## LM test for residual autocorrelation
## test value: 1.9722, p-value: 0.16021

Kali ini \(\hat\rho = -0{,}3257\), bertanda negatif, dengan \(\hat\beta_0 = 0{,}6004\), \(\hat\beta_1 = 0{,}3970\), dan \(\hat\sigma^2 = 0{,}0051\). Tanda negatif akan berarti produktivitas tinggi di tetangga justru menekan wilayah ini, tetapi \(\rho\) tidak signifikan (LR \(p = 0{,}172\); Wald \(p = 0{,}101\)) sehingga tidak perlu ditafsirkan lebih jauh. AIC-nya (\(-3{,}951\)) juga masih sedikit di atas AIC model biasa.

Satu catatan interpretasi: pada model SAR, \(\beta_k\) bukan efek marjinal total, karena perubahan \(x_i\) merambat ke wilayah lain lewat pengganda spasial \((\mathbf{I} - \rho\mathbf{W})^{-1}\). Penguraian menjadi efek langsung, tidak langsung, dan total adalah topik lanjutan. Perlu diingat pula bahwa semua uji di sini bersifat asimtotik, jadi dengan \(n = 5\) nilai-\(p\) hanya berfungsi sebagai ilustrasi prosedur.

6 Perbandingan Hasil

Semua hasil tadi kita kumpulkan dalam satu tabel.

perbandingan <- data.frame(
  Metode = c("OLS", "OLS", "OLS terkoreksi (IV)", "ML manual", "ML lagsarlm", "ML lagsarlm"),
  W      = c("Ketetanggaan", "Invers jarak std", "Ketetanggaan",
             "Ketetanggaan", "Ketetanggaan", "Invers jarak"),
  rho = unname(c(coef(ols_tetangga)["Wy_tetangga"], coef(ols_inv)["Wy_inv"],
                 theta[3], rho_ml, Hasil$rho, Hasil2$rho)),
  b0  = unname(c(coef(ols_tetangga)["(Intercept)"], coef(ols_inv)["(Intercept)"],
                 theta[1], b_ML[1], coef(Hasil)["(Intercept)"], coef(Hasil2)["(Intercept)"])),
  b1  = unname(c(coef(ols_tetangga)["X1"], coef(ols_inv)["X1"],
                 theta[2], b_ML[2], coef(Hasil)["X1"], coef(Hasil2)["X1"]))
)
perbandingan$Label    <- paste0(perbandingan$Metode, "\n(W ", tolower(perbandingan$W), ")")
perbandingan$Kelompok <- ifelse(perbandingan$W == "Ketetanggaan", "W ketetanggaan", "W invers jarak")

kable(perbandingan[, c("Metode", "W", "rho", "b0", "b1")], digits = 4,
      col.names = c("Metode", "Matriks W", "ρ", "β0 (intersep)", "β1 (X1)"))
Metode Matriks W ρ β0 (intersep) β1 (X1)
OLS Ketetanggaan 0.3568 -0.0176 0.3464
OLS Invers jarak std -0.8311 0.9494 0.3726
OLS terkoreksi (IV) Ketetanggaan 0.6437 -0.2057 0.3116
ML manual Ketetanggaan 0.2872 0.0281 0.3549
ML lagsarlm Ketetanggaan 0.2872 0.0281 0.3549
ML lagsarlm Invers jarak -0.3257 0.6004 0.3970

Dugaan \(\rho\) dari tiap metode digambarkan berikut.

perbandingan$Label <- factor(perbandingan$Label, levels = rev(perbandingan$Label))

ggplot(perbandingan, aes(rho, Label, color = Kelompok)) +
  geom_vline(xintercept = 0, color = "grey50") +
  geom_segment(aes(x = 0, xend = rho, yend = Label), linewidth = 1) +
  geom_point(size = 5) +
  geom_text(aes(label = round(rho, 3)), nudge_x = 0.07, hjust = 0, size = 3.8,
            show.legend = FALSE) +
  scale_color_manual(values = c("W ketetanggaan" = biru, "W invers jarak" = oranye), name = NULL) +
  scale_x_continuous(expand = expansion(mult = c(0.1, 0.2))) +
  labs(title = "Dugaan ρ menurut Metode Estimasi",
       subtitle = "Garis abu-abu menandai ρ = 0 (tanpa efek spasial)",
       x = expression(hat(rho)), y = NULL) +
  theme(legend.position = "bottom")

Sementara dugaan \(\beta_1\) sebagai berikut.

ggplot(perbandingan, aes(b1, Label, color = Kelompok)) +
  geom_point(size = 5) +
  geom_text(aes(label = round(b1, 3)), nudge_x = 0.004, hjust = 0, size = 3.8,
            show.legend = FALSE) +
  scale_color_manual(values = c("W ketetanggaan" = biru, "W invers jarak" = oranye), name = NULL) +
  scale_x_continuous(expand = expansion(mult = c(0.1, 0.2))) +
  labs(title = "Dugaan β1 menurut Metode Estimasi",
       x = expression(hat(beta)[1]), y = NULL) +
  theme(legend.position = "bottom")

Dua hal langsung terlihat. Pertama, \(\hat\beta_1\) relatif stabil di kisaran 0,31 sampai 0,40 apa pun metode dan \(\mathbf{W}\)-nya, sedangkan \(\hat\rho\) bergerak dari \(-0{,}83\) (OLS, W invers jarak) hingga \(0{,}64\) (IV). Pengaruh tetangga jauh lebih sulit diduga secara presisi daripada pengaruh variabel penjelas, dan ini sejalan dengan galat baku \(\hat\rho\) yang besar di semua metode.

Kedua, metode yang berbeda memberi \(\hat\rho\) yang berbeda walaupun \(\mathbf{W}\)-nya sama: OLS 0,357, IV 0,644, dan ML 0,287. Karena OLS tidak konsisten, hasilnya tidak dijadikan acuan. ML manual dan lagsarlm() memberi hasil yang identik, yang menegaskan bahwa fungsi tersebut memang mengerjakan prosedur empat langkah tadi. Pergantian \(\mathbf{W}\) dari ketetanggaan ke invers jarak juga mengubah tanda \(\hat\rho\), jadi dalam penelitian nyata pilihan \(\mathbf{W}\) harus didasarkan pada teori dan hasilnya sebaiknya diperiksa dengan beberapa spesifikasi. Untuk model SAR, ML melalui lagsarlm() adalah pilihan standar; IV berguna terutama bila asumsi kenormalan galat diragukan atau ada variabel endogen lain di dalam model.

Pada kedua \(\mathbf{W}\), \(H_0: \rho = 0\) tidak ditolak dan AIC model lag tidak lebih baik daripada model klasik. Jadi pada lima wilayah ini belum ada bukti kuat bahwa efek spasial diperlukan, dan memang bukan itu tujuan latihan ini.

7 Latihan Mandiri

  1. Hitung \((\mathbf{Wy})_C\) dengan tangan (W ketetanggaan), lalu cocokkan dengan kolom Wy_tetangga pada baris C.
  2. Ulangi lagsarlm() dengan W invers jarak yang sudah distandarisasi baris (W_inv_std). Apakah \(\hat\rho\) berubah, dan apakah kesimpulan tentang signifikansinya ikut berubah?
  3. Tambahkan hubungan ketetanggaan A–D pada matriks A, lalu ulangi analisis ML. Bagaimana pengaruhnya terhadap \(\hat\rho\)?
  4. Ubah rho_benar pada simulasi menjadi 0,2 dan 0,8. Pada nilai \(\rho\) mana bias OLS terlihat paling besar?
  5. Tuliskan interpretasi summary(Hasil) dalam satu paragraf, seolah-olah untuk laporan.

Referensi. Anselin, L. (1988). Spatial Econometrics: Methods and Models. Kluwer Academic. Slide kuliah Regresi Spasial, Pertemuan 2: Model SAR.