베이지안 공간 모형과 INLA로 대규모 공간 자료 추정하기

공간 상관을 잠재장으로 분리하는 방법

대규모 자료, 불균일한 표본 위치, 이상치가 함께 있는 공간 데이터에서는 관측값만으로 패턴을 설명하기 어렵다. 베이지안 공간 모형은 관측층의 우도, 공간적 상관을 담는 잠재층, 하이퍼파라미터 사전분포를 위계적으로 묶어 이 문제를 다룬다.

잠재층은 가우시안 과정(GP)이나 가우시안 마르코프 무작위장(GMRF)으로 표현한다. 격자형 데이터에는 ICAR 또는 BYM2를, 연속공간에는 SPDE 기반 표현을 적용할 수 있다.

INLA(Integrated Nested Laplace Approximation)는 잠재 가우시안 모형(LGM)의 주변사후분포를 라플라스 근사로 구하는 결정론적 추론 기법이다. 희소 정밀도 행렬을 이용한 수치 선형대수 최적화와 저차원 하이퍼파라미터 적분, 잠재장 주변화를 결합한다. MCMC와 비교해 시간과 메모리를 크게 줄일 수 있다.

SPDE(확률편미분방정식) 접근은 Matérn 공분산을 가진 GP를 GMRF로 근사한다. 삼각형 메쉬의 기저함수 표현을 사용하므로 희소 구조를 확보할 수 있다. 이때 실용 범위는 ρ = sqrt(8ν)/κ, 마진 표준편차는 σ로 파라미터화하며, PC-사전(Penalized Complexity prior)으로 ρ와 σ에 해석적 제약을 둘 수 있다.

모형 설계에서 함께 결정할 항목

선형예측자 η = Xβ + Zu는 고정효과와 공간·시공간 임의효과를 결합한다. 연결함수를 통해 포아송, 베타-바이노미얼, 제로-인플레이트 같은 비가우시안 우도도 사용할 수 있으며, 누락값·검출한계·오프셋을 모형 내부에서 일관되게 처리한다.

계산 흐름은 메쉬 설계, SPDE 또는 ICAR 구성, 설계행렬과 A-행렬을 포함한 스택(Stack) 생성, 사전 설정, 근사추론, 후처리와 예측으로 이어진다. 희소 정밀도 행렬의 촐레스키 분해와 라플라스 근사를 활용하며, 메쉬와 연결성에 따라 대규모 자료에서 준선형 O(n^1.5)~O(n^2) 스케일링을 달성한다.

격자형 문제에서는 인접 행렬을 바탕으로 ICAR/Besag 또는 BYM2의 공간 평활화를 구성한다. BYM2는 스케일 정규화와 가중 혼합을 사용한다. 연속공간에서는 SPDE(Matŕn)가 불규칙 관측과 경계조건에 적합하며, 경계·홀 처리와 다중 해상도 메쉬 구성도 지원한다.

사전은 추론 결과를 안정시키는 설계 요소다. 예를 들어 P(ρ < 50km)=0.05, P(σ > 1)=0.01 등의 해석적 하이퍼파라미터 사전으로 도메인 지식을 반영할 수 있다. 과적합 억제와 사후 안정성 향상을 위해 민감도 분석도 함께 수행한다.

예측력은 CPO/LPML, WAIC/DIC, 교차검증으로 평가하고 PIT 히스토그램으로 보정 상태를 살핀다. KLD, 사후 경계효과, 메쉬 민감도 분석은 근사오차와 구조적 가정 검증에 사용한다.

질병 지도부터 수요 예측까지

공중보건의 질병 지도화에서는 지역별 발생 건수, 인구 오프셋, 인접 행렬을 입력으로 사용한다. BYM2 공간 임의효과에 사회경제적 공변량을 결합하고, PC-사전으로 과평활을 제어한다. 지역별 상대위험(RR)의 후방 평균과 95% 신뢰구간을 출력해 과다·과소 위험 지역을 찾는다.

환경·대기질 예측에서는 센서 위치와 측정값, 지형·기상 공변량을 사용한다. 불규칙 메쉬 위에 Matérn-SPDE와 비선형 효과(GAM형 spline)를 결합해 격자 지도 예측값, 불확실성 지도, exceedance 확률지도를 만든다.

물류·수요 예측은 시간스탬프 좌표 이벤트와 휴일·날씨 피처를 입력으로 삼는다. 시간 AR(1)과 공간 SPDE를 결합하고 주·일 주기성 성분을 더해 단기 수요 히트맵과 고신뢰 구역 기반 차량 재배치 정책을 산출한다.

산불·재해 위험도에서는 연료지수, 풍속, 경사, 역사 화재 위치를 활용한다. 베타-바이노미얼 또는 제로-인플레이트 모형에 SPDE 공간 상관을 반영해 위험 확률 지도와 자원 배치 우선순위를 만든다.

동등 복잡도 MCMC와 비교하면 추론 시간을 10~~100배 단축할 수 있으며, 데이터·모형·메쉬 조건에 따라 달라진다. 수십만 관측 규모를 처리할 수 있고, 희소 구조 활용으로 메모리 사용을 3~~10배 줄일 수 있다. PC-사전은 과적합과 발산 위험을 낮추며, 결정론적 결과의 재현성과 하이퍼파라미터 탐색 자동화는 운영비용 절감에 기여한다. 예측 평균과 불확실성을 함께 제공하므로 위험 기반 정책 최적화에도 활용할 수 있다.

메쉬와 진단을 반복하는 추론 흐름

입력은 공간 좌표 또는 격자, 관측값, 공변량, 경계 폴리곤, 인접행렬 또는 메쉬 설계 규칙이다. 메쉬 생성 시 최대 요소 크기와 경계 버퍼, 홀 처리를 정하고, 우도·연결함수·PC 사전을 선택한다. 관측과 예측의 설계행렬(A), 오프셋, 가중치를 스택으로 구성한 뒤 INLA의 근사 전략(laplace/gaussian)과 정밀도 분해를 적용한다.

사후 요약, Raster/Vector 지도화 레이어, exceedance 확률, 의사결정 지표를 출력한다. CPO, PIT, WAIC, KLD, 메쉬 민감도를 진단에 포함하고, CPO 실패나 PIT 왜곡이 나타나면 메쉬·사전·근사 전략을 조정한다.

⤢✕

INLA, MCMC, 변분추론의 선택 기준

지표	INLA	MCMC(NUTS)	변분추론
성능(시간)	빠름, 대규모에 우수	정확하나 느림	매우 빠름
확장성	희소구조 활용로 높은 확장성	고차원 잠재장에 취약	매우 높음
일관성/정확성	주변사후 근사 정확도 높음(LGM)	기준선, 매우 정확	하한 편향 위험
안정성	결정론적, 재현성 높음	튜닝·수렴 민감	초기화·설계 민감
운영 편의	메쉬·사전 설계 중심	시뮬 길이·진단 복잡	설정 간단하나 진단 난해

INLA는 Latent Gaussian Model 범주에서 강점이 있다. 비가우시안 잠재장이나 복잡 상호작용에는 추가 근사나 대안을 검토해야 한다.

R-INLA로 연속공간 예측 구성하기

환경은 R >= 4.2, INLA 패키지(안정판)이며 sf/sp 패키지는 선택적으로 사용한다. 설치 명령은 다음과 같다.

install.packages("INLA", repos="https://inla.r-inla-download.org/R/stable")

아래 예제는 연속공간 SPDE와 가우시안 우도를 결합해 임의 지점 예측을 수행한다.

# 패키지
library(INLA)

set.seed(2025)

# 1) 모의 데이터 생성
n <- 500
coords <- cbind(runif(n, 0, 100), runif(n, 0, 100))
beta0 <- 1.0
true_range <- 20  # km 가정
true_sigma  <- 1.0
noise_sd    <- 0.3

# 예측 그리드(출력용)
grid_xy <- expand.grid(x = seq(0, 100, by=5), y = seq(0, 100, by=5))
pred_coords <- as.matrix(grid_xy)

# 2) 메쉬 구성
mesh <- inla.mesh.2d(loc = coords, max.edge = c(5, 15), cutoff = 1, offset = c(5, 20))

# 3) SPDE(Matern)와 PC-사전 설정
spde <- inla.spde2.pcmatern(
  mesh = mesh,
  alpha = 2,                    # nu ~ 1
  prior.range = c(30, 0.05),    # P(range < 30) = 0.05
  prior.sigma = c(1.5, 0.05)    # P(sigma > 1.5) = 0.05
)

# 잠재장 인덱스
s.index <- inla.spde.make.index("spatial", n.spde = spde$n.spde)

# 4) 설계행렬 A(관측/예측)
A_obs  <- inla.spde.make.A(mesh, loc = coords)
A_pred <- inla.spde.make.A(mesh, loc = pred_coords)

# 5) 스택 구성
stack_obs <- inla.stack(
  data = list(y = NA),   # y는 추후 생성
  A = list(A_obs, 1),
  effects = list(s.index, data.frame(b0 = 1)),
  tag = "obs"
)

stack_pred <- inla.stack(
  data = list(y = NA),
  A = list(A_pred, 1),
  effects = list(s.index, data.frame(b0 = 1)),
  tag = "pred"
)

# 6) 잠재장 샘플로 모의 관측 생성(ground truth 실험 용)
Q <- inla.spde.precision(spde, theta = c(log( sqrt(8) / true_range ), log(true_sigma)))
w <- as.numeric(inla.qsample(1, Q))  # GMRF 샘플
eta <- as.vector(A_obs %*% w + beta0)
y <- eta + rnorm(n, 0, noise_sd)

# 관측 y를 스택에 주입
stack_obs <- inla.stack(
  data = list(y = y),
  A = list(A_obs, 1),
  effects = list(s.index, data.frame(b0 = 1)),
  tag = "obs"
)

# 전체 스택
stack_all <- inla.stack(stack_obs, stack_pred)
A_all <- inla.stack.A(stack_all)
y_all <- inla.stack.data(stack_all)$y

# 7) 모형 식
formula <- y ~ 0 + b0 + f(spatial, model = spde)

# 8) INLA 추론
res <- inla(
  formula,
  data = inla.stack.data(stack_all),
  family = "gaussian",
  control.predictor = list(A = A_all, compute = TRUE),
  control.compute = list(dic = TRUE, waic = TRUE, cpo = TRUE, config = TRUE),
  control.inla = list(strategy = "laplace")  # 필요 시 "gaussian", "simplified.laplace" 비교
)

# 9) 관측/예측 분리
idx_pred <- inla.stack.index(stack_all, tag = "pred")$data
pred_mean <- res$summary.linear.predictor$mean[idx_pred]
pred_sd   <- res$summary.linear.predictor$sd[idx_pred]

# 10) 결과 요약
print(res$summary.fixed)             # b0 추정
print(res$summary.hyperpar)          # range/sigma 관련 하이퍼파라미터
print(summary(res$cpo$failure))      # CPO 실패율 확인(0에 근접 권장)

# 11) 간단한 시각 확인(선택)
# image(matrix(pred_mean, nrow=length(unique(grid_xy$x))), main="예측 평균")

메쉬의 max.edge와 cutoff는 관측 밀도와 경계 복잡도에 맞춰 조정한다. 지나치게 세분화하면 계산 비용이 증가한다. CPO 실패율이 높거나 PIT 편향이 있다면 사전을 다시 설정하고, 메쉬를 개선하거나 control.inla 전략 변경 및 이상치 모델링(robust likelihood)을 검토한다.

정확도와 운영 비용 사이의 조정

관측이 밀집한 구역에는 고해상도 메쉬를, 외곽에는 저해상도 메쉬를 적용할 수 있다. 경계 버퍼는 경계왜곡 완화에 도움이 된다. 해상도↑는 정확도↑와 함께 시간·메모리↑라는 트레이드오프를 만든다.

PC-사전은 도메인 지식에 따른 확률적 제약으로 교정한다. 제약이 지나치면 신호가 약해지고, 너무 완화하면 과적합 위험이 커진다. laplace는 정확도에, gaussian은 속도에 강점이 있으며 simplified.laplace는 그 사이의 절충안이다. 모형 복잡성과 데이터 크기를 기준으로 선택한다.

CPO나 WAIC의 개선이 미미하다면 모형을 복잡하게 만들기 전에 공변량 품질, 정규화, 좌표 시스템 오류를 먼저 점검한다. SPDE 또는 ICAR, PC-사전, 반복 가능한 진단 루프를 함께 운영해야 대규모 공간 자료에서도 신뢰 가능한 예측과 불확실성 지도를 제공할 수 있다.