Skip to content

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)

요약

  1. Multivariate regression을 응용해서, outcome이 다변수 벡터인 경우도 메타분석이 가능하다.
  2. 첫번째 저자가 내용을 정리해 R package "mvmeta"를 만들었으니, 활용하면 된다.
  3. 이후에 나온 패키지 "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이 아닌, 메타 분석의 모델부터 보자. θ^iNk(θ,Si+Ψ), i=1,,m. 여기서 i는 study에 대한 인덱스, 그리고 kθ^i의 길이다. Si는 first-stage에서 추정된 θ^i의 분산이다. 그럼 자연스럽게 Ψ는 study 사이의 variability가 되고, Si는 study i안에서 발생한 variability를 의미하게 된다. 이 모델을 살짝 확장하면 meta-regression 모델 θ^iNk(Xiβ,Si+Ψ)가 되는데, Xi는 meta-predictor로 이루어진 행렬이고 β는 meta-regression의 coefficient가 된다. Multivariate regression을 하기 위해서 Xi=Ikxit이고 kronecker product의 정의는 지난 글과는 달리 wikipedia의 정의를 따른다. xi는 당연하게도 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 Si는 서로 독립이라는 가정 때문이다. 평균이 study마다 달라지는게 meta-regression 모델이고, 그 평균은 meta-predictor에 의해서 설명되는 모양새다. 저자들은 linear-mixed model에서 영감을 받았다고 한다. Meta-regression 모델을 θ^i|uiNk(Xiβ+ui,Si), uiNk(0,Ψ)라고 다시 쓸 수 있는데, 여기서 Xiβ는 fixed-effect이고 ui가 random-effect다. 고정 효과와 랜덤 효과가 더해져 "mixed-effect"가 되어 mixed model이라는 이름이 붙었다.

추정 방법

모델을 알아봤으면, 이젠 parameter를 추정할 차례다. Study의 heterogeneity를 설명하는 Ψk×k 모양의 행렬인데, upper-triangular 부분 성분의 개수 k(k+1)/2개의 값이 필요하다. β는 meta-predictor가 각 outcome에 끼치는 영향을 말해주는데, Xi의 column 개수인 k×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=i(θ^iXiβ^)tSi1(θ^iXiβ^)로 정의되며 이는 점근적으로 χ2, df=nq를 따른다. Linear regression의 residual sum of squares (RSS)에서 파생된 것 처럼 생겼다. (θ^iXiβ^)는 random effect가 없을 때 편차로 볼 수 있다. 따라서 Q가 의미하는 것은, 전체 variability의 within-study variability에 대한 비율이 된다. 이 값이 크면 Si로는 다 설명할 수 없는 variability가 있다, 즉 between-study variability 존재의 증거가 되고, 반대로 작으면 대부분의 variability가 within-study variability로 설명된다는 얘기다.

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

정리:

  1. Meta-predictor의 significance를 위해 Wald-test, likelihood ratio test
  2. Between-study variability (heterogeneity)의 significance를 위해 Cochran Q-test
  3. 전체 variability 중 heterogeneity 때문에 발생한 variability의 비율은 I2

예측

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


References

  1. 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.