統計ER

R, EZR, SPSS, KH Coder を使ったデータ分析方法を紹介するブログ。ニッチな内容が多め

トップ > 一致度 > Bland-Altman Plotって何？

Bland-Altman Plotって何？

一致度

Bland-Altman Plotは二つの測定系の結果が一致しているかどうかを確認する図。

回帰直線を合わせると不一致に傾向がないかどうか確認できる。

統計ソフトRでBland-Altman Plotを描いてみた。

Bland-Altman Plotはどうやって描くか？

2つの測定値の平均値をX軸に、差をY軸にしてプロットする。

パッケージを使って描く場合

blandrパッケージを使うと簡単に描ける。

まずblandrパッケージをインストールする。

install.packages("blandr")

パッケージのインストールはこちらも参照のこと。 toukeier.hatenablog.com

インストールが終了したら呼び出す。これで準備完了。

library(blandr)

Bland and Altman 1986 論文からPEFR (Peak Expiratory Flow Rate)のデータを使ってグラフを描く。bland.altman.PEFR.1986データセットの一列目（Wright Meter の測定一回目）と三列目（Mini Wright Meter の測定一回目）を使う。

dat0 <- bland.altman.PEFR.1986
dat <- blandr.data.preparation(dat0[,1],dat0[,3],sig.level=0.95)

datの内容をsummary(dat)で見てみる。二つの測定方法はmethod1とmethod2という名前になる。

> summary(dat)
    method1         method2     
 Min.   :178.0   Min.   :259.0  
 1st Qu.:417.0   1st Qu.:380.0  
 Median :434.0   Median :445.0  
 Mean   :450.4   Mean   :452.5  
 3rd Qu.:494.0   3rd Qu.:512.0  
 Max.   :656.0   Max.   :658.0

Bland-Altman Plotは以下のように記述すると描ける。キモはblandr.draw(method1, method2)。ciShading=TRUEでBiasの範囲（青）とLower limit of agreementの範囲（赤）、Upper limit of agreementの範囲（緑）が塗られる。Limit of agreementよりも離れた点がないかどうか、離れた点が多いかどうかで、一致している具合を判断する。その際、Limit of agreementにも幅があることも考慮に入れる。

with(dat, blandr.draw(method1, method2, 
 ciShading=TRUE, 
 plotTitle="Difference against mean for PEFR data"))

f:id:toukeier:20190505203022p:plain

blandr.display.and.draw()を使うと、グラフの描画とともに分析結果を表示してくれる。また、blandr.output.text()はテキスト出力だけ行う。

with(dat, blandr.display.and.draw(method1, method2, 
 ciShading=TRUE, 
 plotTitle="Difference against mean for PEFR data"))

with(dat, blandr.output.text(method1, method2))

Biasが2つの測定方法の差の平均値。その95%の分布範囲（±1.96SD）がLimits of agreement（LOA）。上限値がUpper LOA（ULoA）、下限値がLower LOA（LLoA）。Bias、Upper LOA、Lower LOAそれぞれの95%信頼区間が計算されている。

Number of comparisons:  17 
Maximum value for average measures:  654 
Minimum value for average measures:  218.5 
Maximum value for difference in measures:  73 
Minimum value for difference in measures:  -81 

Bias:  -2.117647 
Standard deviation of bias:  38.76513 

Standard error of bias:  9.401925 
Standard error for limits of agreement:  16.39491 

Bias:  -2.117647 
Bias- upper 95% CI:  17.81354 
Bias- lower 95% CI:  -22.04884 

Upper limit of agreement:  73.86201 
Upper LOA- upper 95% CI:  108.6177 
Upper LOA- lower 95% CI:  39.10636 

Lower limit of agreement:  -78.0973 
Lower LOA- upper 95% CI:  -43.34165 
Lower LOA- lower 95% CI:  -112.8529 

Derived measures:  
Mean of differences/means:  -1.158314 
Point estimate of bias as proportion of lowest average:  -0.9691749 
Point estimate of bias as proportion of highest average -0.3237992 
Spread of data between lower and upper LoAs:  151.9593 
Bias as proportion of LoA spread:  -1.393562 

Bias: 
 -2.117647  ( -22.04884  to  17.81354 ) 
ULoA: 
 73.86201  ( 39.10636  to  108.6177 ) 
LLoA: 
 -78.0973  ( -112.8529  to  -43.34165 )

blandr.statistics()で全統計値が出力される。

> (ba.stats <- with(dat, blandr.statistics(method1, method2)))
$means
 [1] 503.0 412.5 518.0 431.0 488.0 578.5 388.5 411.0 654.0 439.0 424.5 641.0 263.5
[14] 477.5 218.5 386.5 439.0

$differences
 [1] -18 -35  -4   6 -24 -43  49  62  -8 -12 -15  30   7   1 -81  73 -24

$method1
 [1] 494 395 516 434 476 557 413 442 650 433 417 656 267 478 178 423 427

$method2
 [1] 512 430 520 428 500 600 364 380 658 445 432 626 260 477 259 350 451

$sig.level
[1] 0.95

$sig.level.convert.to.z
[1] 1.959964

$bias
[1] -2.117647

$biasUpperCI
[1] 17.81354

$biasLowerCI
[1] -22.04884

$biasStdDev
[1] 38.76513

$biasSEM
[1] 9.401925

$LOA_SEM
[1] 16.39491

$upperLOA
[1] 73.86201

$upperLOA_upperCI
[1] 108.6177

$upperLOA_lowerCI
[1] 39.10636

$lowerLOA
[1] -78.0973

$lowerLOA_upperCI
[1] -43.34165

$lowerLOA_lowerCI
[1] -112.8529

$proportion
 [1]  -3.5785288  -8.4848485  -0.7722008   1.3921114  -4.9180328  -7.4330164
 [7]  12.6126126  15.0851582  -1.2232416  -2.7334852  -3.5335689   4.6801872
[13]   2.6565465   0.2094241 -37.0709382  18.8874515  -5.4669704

$no.of.observations
[1] 17

$regression.equation
[1] "y(differences) = 0.029 x(means) + -15"

$regression.fixed.slope
means 
0.029 

$regression.fixed.intercept
(Intercept) 
        -15

$regression.equationは、推定回帰直線で、切片（intercept: $regression.fixed.intercept）と傾き（slope: $regression.fixed.slope）が計算されている。これを使ってmethod1とmethod2の平均と差の回帰直線を先ほどのグラフに乗せると以下のようになる。

abline(a=ba.stats$regression.fixed.intercept,b=ba.stats$regression.fixed.slope,lwd=2)
text(400,20,ba.stats$regression.equation)

f:id:toukeier:20190505204049p:plain

回帰分析の結果からもわかる通り、傾きは限りなくゼロに近く、平均の大小で、差の大小が決まるというような傾向はみられない。したがって、系統誤差（一定の傾向を持ったズレ）はないと言える。回帰分析の前に平均meanと差differenceを作成しておく。

dat$difference <- dat$method1 - dat$method2
dat$mean <- (dat$method1+dat$method2)/2

meanのEstimateは、p=0.749で傾きがゼロの帰無仮説が棄却できず、傾きゼロでないとは言えない。

> dat$difference <- dat$method1 - dat$method2
> dat$mean <- (dat$method1+dat$method2)/2
> lm1 <- lm(difference~mean,data=dat)
> summary(lm1)

Call:
lm(formula = difference ~ mean, data = dat)

Residuals:
    Min      1Q  Median      3Q     Max 
-72.201 -21.526  -9.526  14.508  76.980 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)
(Intercept) -15.06750   40.97628  -0.368    0.718
mean          0.02869    0.08821   0.325    0.749

Residual standard error: 39.9 on 15 degrees of freedom
Multiple R-squared:  0.007002,  Adjusted R-squared:  -0.0592 
F-statistic: 0.1058 on 1 and 15 DF,  p-value: 0.7495

blandrパッケージを使わずstep by stepでやってみた

差の平均値bias、差の標準偏差sd.diff、差の標準誤差se.diff、Limits of agreementの標準誤差se.LOA、biasの95%信頼区間bias.ci、biasのlimits of agreement bias.LOA、Upper LOAの95%信頼区間 Upper LOA.ci、Lower LOAの95%信頼区間 Lower LOA.ciを計算する。

(bias <- mean(dat$difference))
(sd.diff <- sd(dat$difference))

(n <- nrow(dat))
(se.diff <- sd.diff/sqrt(n))
(se.LOA <- sqrt((1/n+qnorm(0.05/2,lower.tail=F)^2/(2*(n-1)))*sd.diff^2))

(bias <- mean(dat$difference))
(bias.ci <- bias + c(-1,1)*qt(0.05/2,lower.tail=F,df=n-1)*se.diff)

(bias.LOA <- bias + c(-1,1)*1.96*sd.diff)
(UpperLOA.ci <- bias.LOA[2] + c(-1,1)*qt(0.05/2,lower.tail=F,df=n-1)*se.LOA)
(LowerLOA.ci <- bias.LOA[1] + c(-1,1)*qt(0.05/2,lower.tail=F,df=n-1)*se.LOA)

上記のwith(dat, blandr.output.text(method1, method2))の結果と一致していることが確認できる。

> (bias <- mean(dat$difference))
[1] -2.117647
> (sd.diff <- sd(dat$difference))
[1] 38.76513
> 
> (n <- nrow(dat))
[1] 17
> (se.diff <- sd.diff/sqrt(n))
[1] 9.401925
> (se.LOA <- sqrt((1/n+qnorm(0.05/2,lower.tail=F)^2/(2*(n-1)))*sd.diff^2))
[1] 16.39491
> 
> (bias <- mean(dat$difference))
[1] -2.117647
> (bias.ci <- bias + c(-1,1)*qt(0.05/2,lower.tail=F,df=n-1)*se.diff)
[1] -22.04884  17.81354
> 
> (bias.LOA <- bias + c(-1,1)*1.96*sd.diff)
[1] -78.09730  73.86201
> (UpperLOA.ci <- bias.LOA[2] + c(-1,1)*qt(0.05/2,lower.tail=F,df=n-1)*se.LOA)
[1]  39.10636 108.61766
> (LowerLOA.ci <- bias.LOA[1] + c(-1,1)*qt(0.05/2,lower.tail=F,df=n-1)*se.LOA)
[1] -112.85295  -43.34165

Bland-Altman Plotを以下のスクリプトで描くと、blandr.draw()に近いプロットが描ける。95%信頼区間の範囲を塗りつぶす代わりに、95%信頼区間の上限下限を色付きの点線で描いた。

plot(difference~mean, data=dat, ylim=c(-120,120), xlim=c(200,650),las=1,
ylab="Differences",xlab="Means",main="Difference against mean for PEFR data",pch=20)
abline(h=0)
abline(h=c(bias,bias.ci),lty=c(2,3,3),col=c(1,4,4))
abline(h=c(bias.LOA[2],UpperLOA.ci),lty=c(2,3,3),col=c(1,3,3))
abline(h=c(bias.LOA[1],LowerLOA.ci),lty=c(2,3,3),col=c(1,2,2))

f:id:toukeier:20190505205807p:plain

さらに回帰直線を重ねると以下のようになる。回帰直線はグレーで描きいれた。

lm1 <- lm(difference~mean,data=dat)
abline(a=lm1$coeff[1],b=lm1$coeff[2],col="gray",lwd=2)
text(350,-15,"y = 0.29 x - 15",col="gray")

f:id:toukeier:20190505205951p:plain

まとめ

二つの測定方法の一致を確認する分析方法のBland-Altman Plotを統計ソフトRで描いてみた。blandrパッケージはggplot2を使ったきれいなグラフが描けるのでいいね！また、パッケージを使わずにstep by stepで確認してみたらと理解が深まった。

参考

Bland-Altman Plotを発明したJ. Martin Bland教授のウェブサイト

https://www-users.york.ac.uk/~mb55/meas//ba.htm#top

blandrパッケージのマニュアル

https://cran.r-project.org/web/packages/blandr/blandr.pdf

↑