Khi phân tích dữ liệu không gian, một câu hỏi thường gặp là: làm thế nào để đo lường và trực quan hóa tác động của một thay đổi tại một vị trí lan tỏa đến các vùng lân cận? Bài viết này trình bày quy trình xây dựng mô hình tự hồi quy không gian (SAR), giải thích cơ chế sinh ra hiệu ứng lan truyền, và hướng dẫn chi tiết cách tạo một ảnh động (animated GIF) minh họa quá trình này bằng Stata. Kết quả cuối cùng là một biểu đồ động cho thấy tác động của việc tăng tỷ lệ thất nghiệp tại Dallas lan rộng ra các hạt lân cận ở Texas.
Mô hình tự hồi quy không gian (SAR)
Giả sử chúng ta muốn nghiên cứu tỷ lệ giết người (hrate) ở các hạt Texas như một hàm của tỷ lệ thất nghiệp (unemployment). Mô hình hồi quy tuyến tính chuẩn chỉ xem xét mối quan hệ trong cùng một đơn vị quan sát:
hrateᵢ = β₀ + β₁ unemploymentᵢ + εᵢ
Tuy nhiên, thực tế tỷ lệ tội phạm tại một hạt thường chịu ảnh hưởng từ các hạt xung quanh. Để nắm bắt điều này, mô hình SAR đưa vào một thành phần không gian:
hrateᵢ = γ₁ Σ Wᵢⱼ hrateⱼ + β₁ unemploymentᵢ + β₀ + εᵢ
Trong đó:
- Wᵢⱼ là phần tử của ma trận trọng số không gian W. Giá trị dương nếu hạt j giáp hạt i, bằng 0 nếu không giáp hoặc j = i.
- Σ Wᵢⱼ hrateⱼ là tổng có trọng số của tỷ lệ giết người ở các hạt lân cận, đại diện cho hiệu ứng không gian.
Ma trận W thường được chuẩn hóa (normalized contiguity matrix) để tổng mỗi hàng bằng 1, giúp tham số λ (gamma trong ký hiệu bài gốc) dễ hiểu hơn như một hệ số lan truyền.
Ước lượng mô hình bằng Stata
Stata cung cấp bộ lệnh sp chuyên dụng cho dữ liệu không gian. Quy trình bao gồm: tải dữ liệu, liên kết file hình dạng (shapefile), tạo ma trận trọng số, và ước lượng.
1/* Tải dữ liệu tỷ lệ giết người năm 1990 và lọc Texas */
2copy http://www.stata-press.com/data/r15/homicide1990.dta ., replace
3use homicide1990
4keep if sname == "Texas"
5save texas, replace
6
7/* Tải shapefile liên kết và kiểm tra cài đặt không gian */
8copy http://www.stata-press.com/data/r15/homicide1990_shp.dta ., replace
9spset
10
11/* Tạo ma trận kề chuẩn hóa */
12spmatrix create contiguity W
13
14/* Ước lượng SAR bằng GS2SLS */
15spregress hrate unemployment, dvarlag(W) gs2slsKết quả ước lượng cho thấy tham số không gian λ (hiển thị là W.hrate) bằng 0.341 (p = 0.075), cho thấy bằng chứng về sự tự tương quan không gian. Hệ số thất nghiệp là 0.458 (p = 0.003).
Cơ chế sinh ra hiệu ứng lan truyền
Tại sao mô hình SAR lại tạo ra hiệu ứng lan truyền? Câu trả lời nằm ở dạng ma trận của mô hình:
y = λ W y + Xβ + ε
Giải cho vectơ kết quả y ta được:
y = (I - λW)⁻¹ Xβ + ε
Kỳ vọng có điều kiện E(y|X) = (I - λW)⁻¹ Xβ. Khi thay đổi X từ X₀ (dữ liệu gốc) sang X₁ (gán thất nghiệp Dallas = 10%), sự thay đổi kỳ vọng là:
ΔE(y) = (I - λW)⁻¹ ΔX β (1)
Điều kiện ổn định của mô hình SAR cho phép khai triển ma trận nghịch đảo thành chuỗi lũy thừa (Neumann series):
(I - λW)⁻¹ = I + λW + λ²W² + λ³W³ + ... (2)
Thay (2) vào (1) ta được biểu thức cốt lõi cho ảnh động:
ΔE(y) = ΔXβ + λW ΔXβ + λ²W² ΔXβ + λ³W³ ΔXβ + ... (3)
Ý nghĩa từng hạng tử:
- Hạng 1 (ΔXβ): Tác động trực tiếp, chỉ ảnh hưởng đến Dallas.
- Hạng 2 (λW ΔXβ): Tác động lan truyền đến các hạt giáp Dallas.
- Hạng 3 (λ²W² ΔXβ): Tác động lan truyền đến các hạt giáp các hạt giáp Dallas (hàng xóm cấp 2).
- Tiếp tục vô hạn, độ lớn giảm dần theo lũy thừa λ.
Tạo ảnh động minh họa lan truyền
Ý tưởng: Vẽ lần lượt các biểu đồ tích lũy các hạng tử trong (3). Biểu đồ 1 chỉ có hạng 1, biểu đồ 2 cộng thêm hạng 2, biểu đồ 3 cộng thêm hạng 3... Cuối cùng ghép chuỗi thành GIF.
Mã Stata dưới đây thực hiện chính xác logic toán học trên, sử dụng spgenerate để nhân ma trận W lặp lại.
1/* 1. Lấy tham số lambda đã ước lượng */
2local lambda = _b[W:hrate]
3
4/* 2. Tính Xβ cho dữ liệu gốc (X0) và dữ liệu sửa (X1) */
5predict xb0, xb
6replace unemployment = 10 if cname == "Dallas"
7predict xb1, xb
8
9/* 3. Hạng tử đầu tiên: ΔXβ */
10generate dy = xb1 - xb0
11format dy %9.2f
12
13/* 4. Khởi tạo cho vòng lặp: W^0*y = dy, lambda^0*W^0*y = dy */
14generate Wy = dy
15generate lamWy = dy
16
17/* 5. Vẽ khung hình đầu tiên (chỉ tác động trực tiếp) */
18grmap dy
19graph export dy_0.png, replace
20local input dy_0.png
21
22/* 6. Vòng lặp tính và vẽ 20 khung hình tiếp theo */
23forvalues p = 1/20 {
24 /* Tính W * Wy_trước = W^p * dy */
25 spgenerate tmp = W*Wy
26
27 /* Tính λ^p * W^p * dy */
28 replace lamWy = `lambda'^`p' * tmp
29
30 /* Cập nhật Wy cho vòng sau */
31 replace Wy = tmp
32
33 /* Cộng dồn vào hiệu ứng tổng */
34 replace dy = dy + lamWy
35
36 /* Vẽ và xuất file */
37 grmap dy
38 graph export dy_`p'.png, replace
39 local input `input' dy_`p'.png
40 drop tmp
41}
42
43/* 7. Ghép thành GIF (cần cài đặt ImageMagick hoặc FFmpeg) */
44/* Linux/macOS với ImageMagick: */
45shell convert -delay 150 -loop 0 `input' glsp.gif
46/* Windows với FFmpeg (ví dụ): */
47/* shell ffmpeg -framerate 5 -i dy_%d.png -loop 0 glsp.gif */
48
49/* 8. Dọn dẹp file trung gian */
50shell rm -fR *.pngGiải thích các bước quan trọng:
- Dòng 26:
spgenerate tmp = W*Wythực hiện phép nhân ma trận thưa W với vectơ Wy hiện tại, hiệu quả tính W⁽ᵖ⁾ ΔXβ mà không cần materialize ma trận đầy đủ. - Dòng 27: Nhân với
lambda^pđể có hạng tử λ⁽ᵖ⁾ W⁽ᵖ⁾ ΔXβ. - Dòng 29: Tích lũy vào
dyđể vẽ hiệu ứng cộng dồn. - Dòng 37: Lệnh
shellgọi công cụ bên ngoài (ImageMagickconvert) để ghép các frame.pngthành.gif.
Giá trị đắt giá ✨
Hiệu ứng lan truyền trong mô hình SAR không chỉ là một khái niệm lý thuyết trừu tượng. Bằng cách phân rã ma trận nghịch đảo (I - λW)⁻¹ thành chuỗi lũy thừa, chúng ta biến một phép toán đại số tuyến tính thành một câu chuyện trực quan: sự thay đổi tại một điểm xuất phát, lan qua hàng xóm cấp 1, sau đó cấp 2, cấp 3... và dần tàn đi theo lũy thừa của λ. Kỹ thuật vẽ từng khung hình tích lũy này áp dụng rộng rãi cho mọi mô hình không gian có dạng giảm dần tương tự (SDM, SLX), giúp người xem "thấy" được không gian lan tỏa như một quá trình động học thay vì một kết quả tĩnh cuối cùng.
Câu hỏi tư duy
Nếu tham số không gian λ được ước lượng âm (ví dụ -0.3) thay vì dương, ảnh động sẽ thay đổi như thế nào về mặt mẫu hình màu sắc (hot/cold spots) trên bản đồ? Hãy thử sửa local lambda = -0.3 trong mã và đoán hình dung kết quả trước khi chạy.
