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{\theta}_i \sim N_k (\theta, S_i +\Psi)\), \(i = 1, \cdots, m\). 여기서 \(i\)는 study에 대한 인덱스, 그리고 \(k\)는 \(\hat{\theta}_i\)의 길이다. \(S_i\)는 first-stage에서 추정된 \(\hat{\theta}_i\)의 분산이다. 그럼 자연스럽게 \(\Psi\)는 study 사이의 variability가 되고, \(S_i\)는 study \(i\)안에서 발생한 variability를 의미하게 된다. 이 모델을 살짝 확장하면 meta-regression 모델 \(\hat{\theta}_i \sim N_k (X_i \beta, S_i +\Psi)\)가 되는데, \(X_i\)는 meta-predictor로 이루어진 행렬이고 \(\beta\)는 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 \(\Psi\)와 within-study variability \(S_i\)는 서로 독립이라는 가정 때문이다. 평균이 study마다 달라지는게 meta-regression 모델이고, 그 평균은 meta-predictor에 의해서 설명되는 모양새다. 저자들은 linear-mixed model에서 영감을 받았다고 한다. Meta-regression 모델을 \(\hat{\theta}_i \lvert u_i \sim N_k(X_i \beta + u_i, S_i)\), \(u_i \sim N_k(0, \Psi)\)라고 다시 쓸 수 있는데, 여기서 \(X_i \beta\)는 fixed-effect이고 \(u_i\)가 random-effect다. 고정 효과와 랜덤 효과가 더해져 "mixed-effect"가 되어 mixed model이라는 이름이 붙었다.

추정 방법

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

가설검정 및 관련된 통계량

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

Study 사이에도 결과에 영향을 주는 variability가 있을 수 있고, 그것을 \(\Psi\)로 썼다. 따라서 귀무가설 \(\Psi = 0\)에 자연스럽게 관심이 생기는데, 이는 Cochran Q test를 통해 검정할 수 있다. 통계량은 \(Q = \sum_i ( \hat{\theta}_i - X_i \hat{\beta} )^t S_i^{-1} ( \hat{\theta}_i - X_i \hat{\beta} )\)로 정의되며 이는 점근적으로 \(\chi^2\), \(df = n-q\)를 따른다. Linear regression의 residual sum of squares (RSS)에서 파생된 것 처럼 생겼다. \((\hat{\theta}_i - X_i \hat{\beta})\)는 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{\beta}\)를 계산하는 방법이고, 후자는 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.