Review: Multivariate meta-analysis for non-linear and other multi-parameter associations

Author: A. Gasparrini, B. Armstrong, and M. G. Kenward

Journal: Statistics in Medicine (2012)

요약

Multivariate regression을 응용해서, outcome이 다변수 벡터인 경우도 메타분석이 가능하다.
첫번째 저자가 내용을 정리해 R package "mvmeta"를 만들었으니, 활용하면 된다.
이후에 나온 패키지 "mixmeta"는 이 모델에서 좀 더 일반화된 모델을 다룬다.

소개

지난 글에서 meta-regression의 효율적인 계산을 위해 DLNM의 coefficient를 적절하게 reducing하는 방법을 알아봤다. 이 글에서는 그 coefficient를 가지고 meta-regression을 하는 방법을 알아본다. DLNM의 coefficient를 추정하는 과정은 first-stage, 그 coefficient를 meta-predictor를 통해 회귀분석하는 것을 second-stage라고 했었다. First-stage는 이 글에서 다루지 않을 것이므로, first-stage에서 얻어진 coefficient는 여기서 outcome이라고 부른다.

메타 분석을 하는 이유는, study(e.g., 지역)간 effect size에 차이가 있을 수 있기 때문에 그것을 meta-predictor로 설명하기 위함이다. 이 논문에서는 주로 temperature-mortality와 같은 exposure-response relationship에 중점을 두지만, 계량경제학(실증) 같은 분야에서도 사용될 만 하다고 생각한다. 유명한 책 "총, 균, 쇠"에서 문명의 발달이 인종의 차이에서 유래한 것이 아니고 지역적인 차이에서 유래했다고 말했듯이, 개발도상국과 선진국 사이의 차이를 지역 변수로 설명할 수 있지 않을까? 라고 나이브하게 상상해본다.

모델 설명

거두절미하고, 메타-regression이 아닌, 메타 분석의 모델부터 보자. ${\hat{θ}}_{i} \sim N_{k} (θ, S_{i} + Ψ)$ , $i = 1, \dots, m$ . 여기서 $i$ 는 study에 대한 인덱스, 그리고 $k$ 는 ${\hat{θ}}_{i}$ 의 길이다. $S_{i}$ 는 first-stage에서 추정된 ${\hat{θ}}_{i}$ 의 분산이다. 그럼 자연스럽게 $Ψ$ 는 study 사이의 variability가 되고, $S_{i}$ 는 study $i$ 안에서 발생한 variability를 의미하게 된다. 이 모델을 살짝 확장하면 meta-regression 모델 ${\hat{θ}}_{i} \sim N_{k} (X_{i} β, S_{i} + Ψ)$ 가 되는데, $X_{i}$ 는 meta-predictor로 이루어진 행렬이고 $β$ 는 meta-regression의 coefficient가 된다. Multivariate regression을 하기 위해서 $X_{i} = I_{k} \otimes x_{i}^{t}$ 이고 kronecker product의 정의는 지난 글과는 달리 wikipedia의 정의를 따른다. $x_{i}$ 는 당연하게도 meta-predictor로 이루어진 $p$ 길이의 열벡터이다. 보통 intercept를 포함하기 때문에 맨 첫번째 성분은 1인 경우가 많다.

이번엔 위 모델을 뜯어보면서 왜 저렇게 생겼는지 정당성을 부여할 차례이다. 가장 critical한 질문으로, 왜 저런 정규분포 가정을 하게 되었을까? 주로 first-stage에서 사용되는 DLNM coefficient는 quasi-likelihood를 통해 추정되기 때문에, sample이 적당히 많다면 large sample theory에 의해서 coefficient의 joint distribution을 정규분포로 근사할 수 있다. (지도교수님이 지적해주신 부분, 정확한 process는 내가 정확하게 알고 있지는 않으니 시간 날 때 다시 공부하자.) 분산이 저렇게 생긴 이유는 between-study heterogeneity $Ψ$ 와 within-study variability $S_{i}$ 는 서로 독립이라는 가정 때문이다. 평균이 study마다 달라지는게 meta-regression 모델이고, 그 평균은 meta-predictor에 의해서 설명되는 모양새다. 저자들은 linear-mixed model에서 영감을 받았다고 한다. Meta-regression 모델을 ${\hat{θ}}_{i} | u_{i} \sim N_{k} (X_{i} β + u_{i}, S_{i})$ , $u_{i} \sim N_{k} (0, Ψ)$ 라고 다시 쓸 수 있는데, 여기서 $X_{i} β$ 는 fixed-effect이고 $u_{i}$ 가 random-effect다. 고정 효과와 랜덤 효과가 더해져 "mixed-effect"가 되어 mixed model이라는 이름이 붙었다.

추정 방법

모델을 알아봤으면, 이젠 parameter를 추정할 차례다. Study의 heterogeneity를 설명하는 $Ψ$ 는 $k \times k$ 모양의 행렬인데, upper-triangular 부분 성분의 개수 $k (k + 1) / 2$ 개의 값이 필요하다. $β$ 는 meta-predictor가 각 outcome에 끼치는 영향을 말해주는데, $X_{i}$ 의 column 개수인 $k \times p$ 만큼의 값이 필요하다. 따라서 reducing-dlnm에 나와있는 방법대로 $k$ 를 줄이면 계산 속도가 빨라지는 것이다. 추정은 여러가지 방법으로 할 수 있는데, 이 논문에서는 maximum likelihood (MLE) 그리고 restricted maximum likelihood (REML)를 사용했다. 편미분을 사용해서 closed form을 구하기 힘든 탓에, iterative approach를 통해 수렴을 만들어 추정한다. 먼저 $Ψ$ 를 안다고 가정하고 고정시킨 후 likelihood를 $β$ 에 대해 최대화한다. 구해진 maximizer를 원래 likelihood에 대입해서 $Ψ$ 에 대한 식으로 만든 후 다시 $Ψ$ 에 대한 최적화 과정을 거쳐 maximizer를 구한다. 그리고 이를 수렴할 때 까지 반복하여 REML을 구한다.

가설검정 및 관련된 통계량

가장 궁금한 것은, 실제로 meta-predictor를 포함한 분석이 의미가 있었는지다. 즉 귀무가설 $β = 0$ 을 검정해봐야 한다. 다행히 추정 과정이 Likelihood에 기반했기 때문에 Wald test나 Likelihood Ratio test를 사용할 수 있다. R을 사용해 모델을 학습시키고 summary 함수에 학습한 모델을 입력하면 해당 통계량을 알 수 있다. 정확하게 계산하는 방법은 현재 가물가물하니까, 학기 시작하면 복습하도록 하자.

Study 사이에도 결과에 영향을 주는 variability가 있을 수 있고, 그것을 $Ψ$ 로 썼다. 따라서 귀무가설 $Ψ = 0$ 에 자연스럽게 관심이 생기는데, 이는 Cochran Q test를 통해 검정할 수 있다. 통계량은 $Q = \sum_{i} ({\hat{θ}}_{i} - X_{i} \hat{β})^{t} S_{i}^{- 1} ({\hat{θ}}_{i} - X_{i} \hat{β})$ 로 정의되며 이는 점근적으로 $χ^{2}$ , $d f = n - q$ 를 따른다. Linear regression의 residual sum of squares (RSS)에서 파생된 것 처럼 생겼다. $({\hat{θ}}_{i} - X_{i} \hat{β})$ 는 random effect가 없을 때 편차로 볼 수 있다. 따라서 $Q$ 가 의미하는 것은, 전체 variability의 within-study variability에 대한 비율이 된다. 이 값이 크면 $S_{i}$ 로는 다 설명할 수 없는 variability가 있다, 즉 between-study variability 존재의 증거가 되고, 반대로 작으면 대부분의 variability가 within-study variability로 설명된다는 얘기다.

$1 - R^{2}$ 를 떠올려보자. 이 통계량은 전체 variability에서 linear regression으로 설명되지 않은 variability의 비율을 말해준다. 따라서 이 값이 낮으면 linear regression 모델이 주어진 데이터를 잘 표현한다고 말할 수 있었다. (물론 아주 기본적인 선형모델에서.) 비슷하게, 전체 variability 중 meta-regression으로 설명하지 못한 variability의 비율은 $(Q - n + q) / Q$ 로 잴 수 있다. 이 값을 $I^{2}$ 라고 한다. Meta-predictor를 추가해서 $I^{2}$ 가 줄어든다면 좋은 meta-predictor라고 할 수 있는 근거가 된다.

정리:

Meta-predictor의 significance를 위해 Wald-test, likelihood ratio test
Between-study variability (heterogeneity)의 significance를 위해 Cochran Q-test
전체 variability 중 heterogeneity 때문에 발생한 variability의 비율은 $I^{2}$

예측

Fixed effect만 가지고 prediction을 계산하는 방법과, random effect를 고려하여 계산하는 방법이 있다. 전자는 주어진 meta-predictor $x_{0}$ 를 가지고 design matrix $X_{0}$ 를 만든 후 $X_{0} \hat{β}$ 를 계산하는 방법이고, 후자는 Best linear unbiased predictor (BLUP)를 계산해 사용한다. Bayesian hierarchical model처럼 BLUP는 study 간 정보를 종합하여 좀 더 shrink 된 predicted value를 뱉는다. BLUP의 생김새를 뜯어보면 Gaussian 분포의 conditional mean 처럼 생겼는데, 복습하면서 그 모양을 추측해보도록 하자.

References

Gasparrini A, Armstrong B, Kenward MG. Multivariate meta-analysis for non-linear and other multi-parameter associations. Stat Med. 2012 Dec 20;31(29):3821-39. doi: 10.1002/sim.5471. Epub 2012 Jul 16. PMID: 22807043; PMCID: PMC3546395.