blog

pythonでTOST法による同等性検証の実装¶

1. 初めに¶
2. シチュエーション¶
3. TOST法¶
4. 解析¶
5. 補足¶
6. まとめ¶

1. 初めに¶

同等性検証は、差がないことを検定する手法の一つです。
ここでは、とある製造業でのシチュエーションを仮定し、その使い方を紹介します。

2. シチュエーション¶

あなたはとある製造業で製造技術を担当しており、現在携わっているコスト削減業務において、新規工法の導入を検討しています。
ただし、目的はあくまでコスト削減業務であり、品質は導入前後で同等を担保する事が前提となっています。
新規工法の導入にめどが立ったため、品質が担保できるかどうかテストを行うことになりました。

3. TOST法¶

独立二群のt検定は「有意差がある」もしくは「有意差があるとは言えない」の二択であり、同等性を検定することはできません。
一方で、「差がある」ことをうまく使えば、同等性を示すことが可能です。それが、TOST法と呼ばれる手法です。
TOST法は「Two One-Side Test」の略で、その名の通り「Two：2回」の「One-Side Test：片側検定」からなります。
以下に、その手順と原理を説明します。

比較元となるデータ（今回のケースだと新規工法導入前）の評価指標について平均値を算出する。
1で求めた平均値に対して、どの程度のブレであれば許容できるか検討する。
片側検定１回目：比較対象の平均値が「比較元ー許容幅」以上であることを検定する。
片側検定２回目：比較対象の平均値が「比較元＋許容幅」以下であることを検定する。

今回は、品質保証部と協議の結果、以下内容で検証を進めることになりました。

測定項目は製品の寸法（幅）の平均値とする。
サンプルサイズは、新規工法導入前後でそれぞれn=50とする。
新規工法導入前後で平均値の差が±5%以内に収まっていれば、同等の品質が担保されていると判断する。
同等性検証を行う上でのp値の閾値は、多重検定問題を考慮し、0.025とする。

[補足]
実際には、平均値の差が何％以内に収まっていれば同等とみなすのかを決めることが一番難しいです。計算自体はデータがそろっていればpythonで簡単に実装できるためです。
例えば、ある機器を用いた測定値で同等性検証を行う場合、その機器の測定誤差を参考にする方法が考えられます。

4. 解析¶

比較元データと最初のテストデータとして、以下のデータを得ました。
※どのようなデータかは最後に種明かしします。

In [1]:

from scipy import stats
import numpy as np
import pingouin as pg
import plotly.graph_objects as go

In [ ]:

array_ref = [
    11.8, 10.4, 11.0, 12.2, 11.9, 9.0, 11.0, 9.8, 9.9, 10.4, 10.1, 11.5, 10.8,
    10.1, 10.4, 10.3, 11.5, 9.8, 10.3, 9.1, 7.4, 10.7, 10.9, 9.3, 12.3, 8.5,
    10.0, 9.8, 11.5, 11.5, 10.2, 10.4, 9.1, 8.0, 9.7, 10.2, 11.2, 11.2, 9.6,
    9.7, 9.0, 8.6, 8.3, 12.0, 9.5, 9.6, 8.7, 10.8, 8.4, 9.8
    ]

array_test_1st = [
    11.8, 11.2, 9.0, 10.5, 11.1, 9.8, 10.8, 10.6, 10.5, 10.3, 10.9, 11.7, 9.5,
    11.5, 10.7, 10.9, 9.4, 10.6, 12.0, 9.4, 8.5, 8.8, 10.8, 12.9, 11.6, 12.2,
    10.6, 11.9, 10.2, 11.1, 10.2, 10.0, 10.6, 10.0, 11.8, 10.7, 10.9, 10.2, 11.8,
    9.8, 11.2, 10.1, 9.6, 10.0, 9.7, 10.3, 10.2, 10.8, 11.1, 10.4
    ]

どのような分布になっているか、まずはヒストグラムで確認してみましょう。
また、それぞれの基本的な統計量も併せて表示させます。

In [6]:

fig = go.Figure()

fig.add_trace(go.Histogram(
    x = array_ref,
    name = "ref",
    opacity = 0.6
))

fig.add_trace(go.Histogram(
    x = array_test_1st,
    name = "test_1st",
    opacity = 0.6
))

fig.update_layout(barmode = "overlay")
fig.update_xaxes(title = "width_mean(mm)", range = [5, 15])
fig.update_yaxes(title = "count(n)")
fig.add_annotation(text = 
                   f"""
                   mean_ref = {np.mean(array_ref):.2f}mm<br>
                   std_ref = {np.std(array_ref, ddof = 1):.2f}<br>
                   mean_test_1st = {np.mean(array_test_1st):.2f}mm<br>
                   std_test_1st = {np.std(array_test_1st, ddof = 1):.2f}<br>
                   """,
                   align = "left", x = 14, y = 8, showarrow = False)
fig.show()

新規工法導入前の平均値は10.14mmであるのに対し、導入後は10.60mmでした。10.14mmの5%上は10.647であり、平均値だけ見れば同等と判断出来そうです。
一方で、両者ともにそれなりに裾野の厚い分布となっており、本当に平均値だけ見て判断していいでしょうか。
そこで、TOST法による同等性検証を実施します。実装は非常にシンプルで、pingouinという統計解析ライブラリを用いれば1行で完結します。

In [7]:

mean_ref = np.mean(array_ref)
mergin_percent = 0.05
mergin = mean_ref * mergin_percent
res = pg.tost(x = array_ref, y = array_test_1st, bound = mergin, correction = False)
res

Out[7]:

	bound	dof	pval
TOST	0.5072	98	0.40997

boundは許容マージン幅、dofは自由度、pvalはp値を示し、特にp値は2回行う検定のうち大きいほうのp値を返します。
p値を確認すると、0.41程度となり、閾値の0.025より大幅に大きい値となりました。平均値だけ見ると問題なさそうでしたが、
裾野が厚い分、同等とは言えないという結果になってしまったようです。

そこで、新規工法のパラメータを調整し、再度テストを実施し、以下のデータを得ました。

In [10]:

array_test_2nd = [
    11.0, 10.3, 10.5, 7.8, 9.0, 10.7, 11.0, 9.1, 10.6, 10.5, 10.4, 9.6, 9.3, 
    9.3, 8.8, 10.3, 10.1, 11.7, 11.2, 10.5, 9.9, 6.9, 11.2, 11.4, 9.9, 10.0,
    10.0, 11.7, 10.2, 8.0, 9.2, 9.0, 11.0, 12.2, 9.2, 8.5, 10.6, 8.5, 10.4,
    9.0, 9.8, 10.4, 9.5, 11.6, 10.7, 10.7, 10.5, 11.5, 11.3, 11.5
    ]

前回のテストデータとも合わせ、ヒストグラムおよび基本統計量を表示させます。

In [61]:

fig = go.Figure()

fig.add_trace(go.Histogram(
    x = array_ref,
    name = "ref",
    opacity = 0.6
))

fig.add_trace(go.Histogram(
    x = array_test_1st,
    name = "test_1st",
    opacity = 0.6
))

fig.add_trace(go.Histogram(
    x = array_test_2nd,
    name = "test_2nd",
    opacity = 0.6
))

fig.update_layout(barmode = "overlay")
fig.update_xaxes(title = "width_mean(mm)", range = [5, 15])
fig.update_yaxes(title = "count(n)")
fig.add_annotation(text = 
                   f"""
                   mean_ref = {np.mean(array_ref):.2f}mm<br>
                   std_ref = {np.std(array_ref, ddof = 1):.2f}<br>
                   mean_test_1st = {np.mean(array_test_1st):.2f}mm<br>
                   std_test_1st = {np.std(array_test_1st, ddof = 1):.2f}<br>
                   mean_test_2nd = {np.mean(array_test_2nd):.2f}mm<br>
                   std_test_2nd = {np.std(array_test_2nd, ddof = 1):.2f}<br>
                   """,
                   align = "left", x = 14, y = 8, showarrow = False)
fig.show()

最初のテストよりも平均値が新規工法導入前に近くなっています。平均値としては申し分ないでしょう。
一方で、値が小さい側に少々裾野が厚いようにも見えます。TOST法を適用し、p値の観点からも同等と言えるか検定してみましょう。

In [62]:

res = pg.tost(x = array_ref, y = array_test_2nd, bound = mergin, correction = False)
res

Out[62]:

	bound	dof	pval
TOST	0.5072	98	0.01773

p値としても0.025を下回り、新規工法導入前後で品質は同等であると統計的に判断できる結果となりました。
この結果を持って、新規工法を導入できることとなりました。

5. 補足¶

今回の検証において用いたデータセットは以下コードで作成しました。リファレンスの平均値を10と設定していたので、許容幅を5%とすれば、±0.5のずれまで許容されることになります。
1stデータは許容幅のライン上、2ndデータは許容幅に収まる平均値に設定していました。解析結果はこのデータセットの前提にもそぐう内容となっていることがわかります。

In [ ]:

array_ref = np.round(stats.norm.rvs(loc = 10, scale = 1, size = 50, random_state = 0), 1)
array_test_1st = np.round(stats.norm.rvs(loc = 10.5, scale = 1, size = 50, random_state = 10), 1)
array_test_2nd = np.round(stats.norm.rvs(loc = 10.1, scale = 1, size = 50, random_state = 20), 1)

今回TOST法の実装はpingouinライブラリを使用しましたが、中身は通常のt検定となんら変わりません。以下、おなじみscipy.statsでも同じ結果となるか確認します。

In [70]:

res_stats_greater = stats.ttest_ind(a = array_ref, b = array_test_2nd - mergin, equal_var = True, alternative = "greater")
res_stats_less = stats.ttest_ind(a = array_ref, b = array_test_2nd + mergin, equal_var = True, alternative = "less")
print(f"「比較元-許容幅」以上に対するp値: {res_stats_greater[1]:.4f}")
print(f"「比較元+許容幅」以下に対するp値: {res_stats_less[1]:.4f}")

「比較元-許容幅」以上に対するp値: 0.0105
「比較元+許容幅」以下に対するp値: 0.0177

2回目のt検定、つまり「比較元＋許容幅」以下に対するp値の方が大きく、この値がpingouinのTOST結果と一致しています。

6. まとめ¶

いかがだったでしょうか。同等性検証という名前を聞くと仰々しいですが、中身はいたってシンプルです。
一番難しいのは、補足でも書いたとおり、許容幅の設定です。ここをいかに関係部署と握れるかがキモと言っても過言ではありません。