Bland-Altman plot 이해하고 그리기(Python) - 두 검사 방법의 일치도 평가 방법(Method comparison)

두 방법의 일치도 평가(Method comparison)을 위한 통계 기법의 필요성

임상 연구를 할 때 두 측정 방법의 일치도를 평가해야하는 경우가 있다. 예를 들어 A 기기로 측정한 결과와 B 기기로 측정한 결과가 있을 때, 두 결과가 일치함/비슷함/동등함을 보여주어야 하는 것이다. 일반적으로는 A 기기가 일종의 'gold standard'이자 전통적인 기기, B 기기가 우리가 도입하고픈 새로운 기기일 경우가 많아서, 이 새로운 기기가 전통적인 기기와 비슷한 결과를 낸다는 것을 주장하기 위해 사용되기도 한다. 

 

이를 위해 일반적인 t-test 나 상관관계 분석을 해서는 안된다. 즉 t-test를 통해 두 측정값이 차이가 없다는 것을 보여주는 방법도 안되고, 높고 유의한 상관계수를 보여주어서도 안된다는 것이다. t-test는 우선 영가설을 기각하는 방식의 가설 검증 기법이므로, 차이가 있다는 가설이 기각되었다고 해서 '차이가 없다는 사실이 통계적으로 유의하다' 주장할 수 없기 때문이다. 한편, correlation의 경우 같은 방향으로 비슷한 정도로 움직인다 정도는 보여줄 수 있지만, 두 측정값이 1:1 관계가 아니더라도 상관관계는 높게 나올 수 있다. 그러나 실제 두 검사 방법의 '일치도'를 보여주어야 하는 method comparison 기법으로는 적절하지 않다. 

 

Bland-Altman Plot 

구성요소 이해하기

Bland-Altman plot은 두 방법으로 특정된 측정값들의 일치도를 평가하기 위해 사용되는 통계적 방법이다.

일단 데이터로는, 같은 대상(subject)에 대해 A, B 방법으로 측정한 값이 있을 것이다.

우선 그림을 보며 이해해보자.

<X, Y axis>

1. X축: A와 B로 측정한 측정값들의 평균값

2. Y축: 두 측정값의 차이

즉 한 샘플에 대해서 두 개의 측정값이 있으면, 둘의 평균과 차이 값을 구해서 각각 x, y축 삼아 그리는 것이다. 그 데이터들을 그리면 위의 dot 들이 그려지게 된다. 

<Mean difference>

3. 진한 검은색 선: 두 측정 값의 차이의 평균 (mean difference). 추정된 편향(estimated bias를 의미함)

4. 진한 검은색 선 바로 위의 점선: mean difference =0 선. 두 측정 값 차이의 평균이 0일 때의 선 (일종의 기준 선)

진한 검은색 선인 mean difference 선은 4번 mean difference = 0 선과 비교해서 보아야 한다. mean difference = 0 선은 bias 가 없다는 것을 의미한다. 즉 검은 실선이 0에 가까울 수록 두 방법 간에 bias 가 작다는 것을 의미하고, mean difference가 양수나 음수이면 양의 bias 혹은 음의 방향으로 bias가 있다는 것을 의미한다. mean difference가 0보다 유의하게 다르면 fixed bias, 즉 systematic difference가 있다고 봐야한다. 보통 mean difference의 confidence interval 이 0을 포함하고 있으면 bias가 유의하지 않다고 간주할 수 있다. 

<Limit of agreements>

5. 가장 위쪽, 아래쪽 점선: Limit of agreements (일치도의 한계?). 보동 95% limit of agreement를 사용한다. 두 측정값의 차이의 standard deviation을 구한 뒤, 1.96을 곱한 후 mean에서 더하고 빼주어 구한다(mean+/-1.96SD). 보통 두 가지 방법에 의한 측정이 얼마나 떨어지는져 있는지 보여준다. 이 상한과 하한 값을 바탕으로 임상적으로 받아들여질 수 있을만한 차이인지 아닌지 판단하게 된다. 어떻게 활용하는지는 아래에서 보다 구체적으로 다루겠다.

해석하기(+유의점)

1. 각 차이값들이 정규분포를 따라야 한다(normality 가정).

본 분석은 두 측정치 간의 차이값들이 정규분포를 따른다는 가정을 하고 있다. 이에 따르면 실제로 95%의 데이터가 95% limit of agreement 내에 들어와야 한다. 이것부터 먼저 확인하자. 가정에 어긋나면 로그화 등의 조치를 취해서 정규분포를 따를 수 있도록 해야한다. 정규성을 검증할 수 있는 통계적 방법까지 제대로 사용한다면 best.

2. Bland-Altman Plot은 그 자체로 두 방법 간의 일치도가 충분한지 알려주지 않는다.  사전 기준이 필요하다.

예를 들어 t-test의 경우에는 일반적으로 p-value가 0.05 미만이면 유의하다고 볼 수 있는데, 여기엔 그런게 없다. mean difference의 confidence interval이 0을 포함하지 않아 systematic bias가 있어도, 이것이 A 방법 대신 B 방법을 사용하지 못한다는 것을 시사하지는 않는다. 만약 bias가 일관적이기만 하면, B방법으로 측정한 후 그 bias값 만큼 숫자를 빼서 사용하면 되기 때문이다. Limit of agreements도 마찬가지이다. 그 자체로 Limit of agreement 범위가 너무 넓다거나 좁다거나의 사실을 알려주지 않는다. 

 

즉, Bland-Altman plot은 그저 bias와 range를 알려줄 뿐이다. 그것들이 충분한지는 말해주지 않는다.

 

중요한 것은 사전에(a priori) 임상적인 기준을 만드는 것이다. 기대되는 limit of agreement, 최대로 수용가능한 limit of agreement 등을 정해야 한다. 이것은 보통 다른 임상적인 기준에 의해 만들어질 수 있다. 기존의 논문을 참조하여 설정할 수 도 있고, 분석의 목적에 따라 달라진다. 즉, 사전에 마련된 임상적 기준(최대로 수용가능한 limit of agreement의 범위)를 준비한 뒤, Bland-Altman plot을 그린다 -> 실제 limit of agreement가 사전 기준 내에 있으면 두 방법은 충분한 일치도를 보인다고 할 수 있고, 그렇지 않으면 일치도를 보이지 않는다(서로 바꿔 사용할 수 없다)고 결론내릴 수 있다. Bias를 허용할 수 있는지 없는지, 허용한다면 어느 정도를 허용할 수 있는지 또한 사전에 임상적 기준으로 마련되어야 한다. 

3. mean difference, limit of agreements 의 confidence interval (CI)를 구하는 것이 좋다.

결국 모두 표본을 기반으로 해석하고 추론해야한다. 측정값의 precision을 구해서 같이 plotting하면 보다 정확한 정보를 전달할 수 있다. 특히 해석을 할 때도 단순히 mean difference가 얼마인지, 0과 얼마나 떨어져 있는지 보기보다는 mean difference의 범위(CI)가 0을 포함하는지 등을 말하는 것이 좋다. limit of agreement의 CI 또한 그것이 사전에 마련한 기준을 포함하는지, 기준 내에 있는지 등을 판단하는데 효과적으로 활용될 수 있다. 또한 CI가 너무 넓으면 sample이 부족하다는 걸 의미할 수도 있기에.. (특히 통계에 빡센 리뷰어를 만나면 이걸 요청받을 수 있다 ^^) 

 

종합하면,

분석의 순서는 다음과 같다.

1) 분석 목적에 따른 임상적 기준 마련 

2) 데이터 수집

3) 정규성 확인

4) Bland-Altman plot 그리기

5) 1)에서 마련한 기준에 부합하는지 해석

 

 

Python 으로 계산하고 그리기

기본적으로 statsmodels 에서 mean_diff_plot function을 제공하고 있다 (링크).

import statsmodels.api as sm
import numpy as np
import matplotlib.pyplot as plt

# 임의의 데이터 만들기
np.random.seed(9999)
m1 = np.random.random(20) # 첫번째 측정값
m2 = np.random.random(20) # 두번째 측정값

f, ax = plt.subplots(1, figsize = (8,5))
sm.graphics.mean_diff_plot(m1, m2, ax = ax)
plt.show()

 

https://www.statsmodels.org/dev/generated/statsmodels.graphics.agreement.mean_diff_plot.html

mean_diff_plot function을 이용하면 위의 plot을 얻을 수 있는데, 문제는 confidence interval을 그려주지 않는다는 점 ㅜㅜ 

그래서 소스코드를 변형해보았다.

 

원본 소스코드 링크는 여기. 비교해보면 차이가 보일텐데, CI 계산하는 부분만 추가되었고 주석을 좀 더 추가해보았다.

import numpy as np
import scipy.stats as stats
# source: https://www.statsmodels.org/dev/_modules/statsmodels/graphics/agreement.html#mean_diff_plot

def mean_diff_plot_ci(m1, m2, sd_limit=1.96, ax=None, fig=None, scatter_kwds=None,
                   mean_line_kwds=None, limit_lines_kwds=None):
    
    ## --- check inputs
    if len(m1) != len(m2):
        raise ValueError('m1 does not have the same length as m2.')
    if sd_limit < 0:
        raise ValueError('sd_limit ({}) is less than 0.'.format(sd_limit))
    
    ## --- calculate means of two measures/statistics of difference between two measures
    # get sample size 
    n = len(m1)

    # calculate means of two measures
    means = np.mean([m1, m2], axis=0)

    # calculate statistics of difference between two measures
    diffs = m1 - m2
    mean_diff = np.mean(diffs)
    std_diff = np.std(diffs, axis=0)
    
    # variance of difference between two measures
    var = std_diff**2
    # Standard error of the mean diff
    se_mean_diff = np.sqrt(var/n)
    # standard error of the limits of agreement
    se_loas = np.sqrt(2.92*var/n)
    
    # Endpoints of the range that contains 95% of the Student’s t distribution
    t_interval = stats.t.interval(confidence=0.95, df=n - 1)
    # Confidence intervals of mean of difference between two measures
    ci_bias = mean_diff + np.array(t_interval) * se_mean_diff

    
    # Plot means, means against diffs, etc.
    scatter_kwds = scatter_kwds or {}
    if 's' not in scatter_kwds:
        scatter_kwds['s'] = 20
    mean_line_kwds = mean_line_kwds or {}
    limit_lines_kwds = limit_lines_kwds or {}
    for kwds in [mean_line_kwds, limit_lines_kwds]:
        if 'color' not in kwds:
            kwds['color'] = 'gray'
        if 'linewidth' not in kwds:
            kwds['linewidth'] = 1
    if 'linestyle' not in mean_line_kwds:
        kwds['linestyle'] = '--'
    if 'linestyle' not in limit_lines_kwds:
        kwds['linestyle'] = ':'
        
    ax.axhline(0, color='black', alpha=0.8)  # plot mean difference = 0 line.
    ax.scatter(means, diffs, **scatter_kwds) # Plot the means against the diffs.
    ax.axhline(mean_diff, **mean_line_kwds)  # draw mean line.

    # Annotate mean line with mean difference.
    ax.annotate('mean diff:\n{}'.format(np.round(mean_diff, 2)),
                xy=(0.99, 0.5),
                horizontalalignment='right',
                verticalalignment='center',
                fontsize=14,
                xycoords='axes fraction')
    
    ## --- get limit of agreements and CI of them + plot
    if sd_limit > 0:
        half_ylim = (1.5 * sd_limit) * std_diff
        half_ylim = (1.5 * sd_limit) * std_diff
        
        ax.set_ylim(mean_diff - half_ylim,
                    mean_diff + half_ylim)
        limit_of_agreement = sd_limit * std_diff
        lower = mean_diff - limit_of_agreement; 
        upper = mean_diff + limit_of_agreement
        
        ci_upperloa = upper + np.array(t_interval) * se_loas
        ci_lowerloa = lower + np.array(t_interval) * se_loas

        # plot limit of agreements
        for j, lim in enumerate([lower, upper]):
            ax.axhline(lim, **limit_lines_kwds)
        ax.annotate(f'-{sd_limit} SD: {lower:0.2g}',
                    xy=(0.99, 0.07),
                    horizontalalignment='right',
                    verticalalignment='bottom',
                    fontsize=14,
                    xycoords='axes fraction')
        ax.annotate(f'+{sd_limit} SD: {upper:0.2g}',
                    xy=(0.99, 0.92),
                    horizontalalignment='right',
                    fontsize=14,
                    xycoords='axes fraction')

    elif sd_limit == 0:
        half_ylim = 3 * std_diff
        ax.set_ylim(mean_diff - half_ylim,
                    mean_diff + half_ylim)
        


    # Plot the confidence intervals
    left, right = ax.get_xlim()
    ax.plot([left] * 2, list(ci_upperloa), c='grey', ls='--', alpha=0.5)
    ax.plot([left] * 2, list(ci_bias), c='grey', ls='--', alpha=0.5)
    ax.plot([left] * 2, list(ci_lowerloa), c='grey', ls='--', alpha=0.5)

    # Plot the confidence intervals' caps
    domain=right-left
    x_range = [left - domain * 0.025, left + domain * 0.025]
    ax.plot(x_range, [ci_upperloa[1]] * 2, c='grey', ls='--', alpha=0.5)
    ax.plot(x_range, [ci_upperloa[0]] * 2, c='grey', ls='--', alpha=0.5)
    ax.plot(x_range, [ci_bias[1]] * 2, c='grey', ls='--', alpha=0.5)
    ax.plot(x_range, [ci_bias[0]] * 2, c='grey', ls='--', alpha=0.5)
    ax.plot(x_range, [ci_lowerloa[1]] * 2, c='grey', ls='--', alpha=0.5)
    ax.plot(x_range, [ci_lowerloa[0]] * 2, c='grey', ls='--', alpha=0.5)
    interval = 0.01 if right > 1 else 0.0000001
    ax.fill_between(np.arange(left, right, interval), ci_upperloa[0], ci_upperloa[1], alpha=0.3, color='lightgrey')
    ax.fill_between(np.arange(left, right, interval), ci_lowerloa[0], ci_lowerloa[1], alpha=0.3, color='lightgrey')

    # labels and ticks
    ax.set_ylabel('Difference', fontsize=15)
    ax.set_xlabel('Means', fontsize=15)
    ax.tick_params(labelsize=13)
    fig.tight_layout()
    return fig

또 차이가 있다면 기존 statsmodel 함수에는 ax만 넣으면 되는데 변형된 함수에는 fig까지 넣어주어야 한다. 즉 적용 코드는 아래와 같다.

f, ax = plt.subplots(1, figsize = (8,5))
mean_diff_plot_ci(m1, m2, ax = ax, fig=f)
plt.show()

코드를 돌리면 위와 같은 결과를 얻을 수 있다! 

좌측의 점선 error bar가 각각의 confidence interval이다.

 

 

참고사항

• 또 다른 method comparison 방법으로 intraclass correlation coefficient를 활용할 수 있다. 일반적인 상관계수와는 다름. 그러나 종종 이것도 상관계수라고 안 좋아하는 reviewer가 있을 수 있다.
• r code이지만 bland-altman plot을 위한 sample size 계산 코드도 있다(링크)

참고 자료

• 도움이 가장 많이된 논문!!!! (링크) 쉽고 친절하고 체계적이다 ㅜㅜ 
• 위키도 참고했다.
• 좀 더 step-by-step 계산 과정 및 코드를 보려면 이곳을 참고. 나는 confidence interval 파트를 참고했다. 
• R code로 그리고 싶으면 여기를 참고.
• Bland-Altman plot에 대해 다룬 영문 블로그(링크)