R-Tableau勉強記３：t検定を実装する

今回は2標本のt検定を取り上げます。

データはSample SuperStoreを使用します。

また今回使用したWorkbookは以下からダウンロードできます。

https://drive.google.com/file/d/11SnJYnJagLXdMLa74V-YxE6KxmBoVAp0/view?usp=sharing

参考

Logics of Blue: t検定の考え方

Logics of Blue; 1章 t検定（数式なし）

Logics of Blue: 2章 t検定（数式あり）

R-Source: 65. 二標本検定

Introductory Statistics with R (古い本ですが統計学のR実装を学ぶにはいい本でした。そのうちもう少し最近の本を参照するようにします）

---

2標本のt検定とは（詳細は割愛）

詳しい話は参考資料を見て頂くとして、ざっくり言えば「2つのデータの平均値の違いが統計的に有意かどうか調べる」ための手法と理解しています。

例えば今回はSuperStoreの「WestとEastのOrder IDごとの合計Salesの平均値」を比較します。その結果は以下のように、約$34の差がありました。

この$34は統計的に有意な差と言えるのかどうか（厳密には両データの平均値の差がゼロであるという仮説が棄却されるかどうか）を知るために2標本t検定を使用します。

実際にはデータが正規分布かどうか、2つのデータの分散が等しいかを気にかけるべきなのですが、まずは実装を考えましょう。

---

t検定の実装

以降、2標本t検定を単にt検定と呼びます。

まずは今回扱うデータを見てみます。

今回の入力には「OrderIDごとの合計売上」を使用したいので、Order IDをViz-LODに含めます。

このあたりの内容についてはこちらの記事を参照ください。

次に、計算フィールドでt検定の結果を返す実装を考えるのですが、まずはR上でt検定の実行結果がどのような形で得られるのかを確認します。

（等分散かどうかの確認は後程行うとして、まずは等分散を仮定した場合のt検定の結果を見てみます）

上記はR Studio上でwest: RegionがWestのOrderIDごとの合計Salesのデータ、east: RegionがEastのOrderIDごとの合計Salesのデータを等分散を仮定したt検定にかけた結果です。

Tableauの実装上で大事なことは2枚目で、t検定は出力をリスト型で出します。

今回Tableauにはp値（データの平均値の差がゼロである確率）が欲しいので、上記を見るとp値は以下のような実装から得られることが分かります。

p_value = t.test(west$Sales,east$Sales,var.equal=T)$p.value

方針は分かったので、Tableau上で入力ベクトルを各計算フィールドで持たせれば良さそうなので、Tableauでの実装を考えていきます。

まず、入力ベクトルのサイズを確かめます。

今回は各RegionのOrder IDごとの合計Salesについて扱いたいので、COUNTD(Order ID)がそのまま想定されるベクトルサイズ（t検定に使うサンプル数）になります。

ということで、以下のような計算式とビューを作成し、また入力ベクトルのサイズを見てみましょう。

West Sales（East Salesも同様)

IF [Region] = 'West'

THEN [Sales]

END

vector size (west)（Eastも同様)

SCRIPT_REAL(

'length(.arg1)'

, SUM([West Sales])

)

ラベルには対応するvector sizeを使用しています。

また表計算はOrder IDを使用して計算させます。

ここでVizLODにRegionを使用していないことに注意してください。

二つのデータはMeasure Namesで分けています（これはデータをLODの意味で分けないので、このビューのViz-LODはOrder IDのみです）。

さて、本来East/Westに対して470/538が想定されていたのに対して、入力ベクトルは二つの値の合計値を返してきました。これは何故でしょうか？

答えは「COUNTDを計算した表と違い、上のビューではRegionがViz-LODにないから」です。つまりSCRIPT_REALでRに入力されるベクトルは、計算式上でRegionごとに分けようとしても、実際にはデータがRegionで分けられていないため、両Regionに含まれるOrder ID全てが入力されるので、上記のようにベクトルサイズが1008となります。

t検定の式（というよりt値の式）には2データのサンプルサイズを使用するため、入力ベクトルのサイズは正しくなければいけません。

したがって、vector sizeの式を修正する必要があります。

vector size (west)（修正版）

SCRIPT_REAL(

'length(subset(.arg1,!is.na(.arg1)))'

, SUM([West Sales])

)

修正版では、各入力から「NA: 欠損値を抜いたベクトルを使用する」ための記述を足しています。

（ところでTableauからRに入る際の欠損値はNULLでなくNAなのですね）

実際に修正版の結果を見てみると、以下のようになります。

たしかに入力ベクトルのサイズが期待しているものと一致しました。

正しい入力ベクトルが得られたので、実際にt検定でのp値を出してみます。

p-value of t_test

SCRIPT_REAL(

't.test(

subset(.arg1,!is.na(.arg1))

, subset(.arg2,!is.na(.arg2))

, var.equal=T)$p.value'

, SUM([East Sales]), SUM([West Sales])

)

R Studioでの実行結果と一致しているので、正しい結果が得られていることが分かります。

---

余談：実は使う関数によっては自動でNAを省いてくれる

subset()を使用して欠損値を省いて正しいベクトルを入力しましょう、という議論を上でしましたが、実は今回のケースでは欠損値処理をしなくても正しいp値が返ってきます。

というのも、t.test()にはna.actionというパラメータがあり、こちら欠損値をどう処理するかを指定するパラメータです。

デフォルトでは欠損値を省いて検定するようになっているのですが、そもそも入力をある程度コントロールした方が良いのでは、と思い上記の議論を書きました。

---

さて、上記までで「2019年度のOrder IDごとのSalesの平均値は、EastとWestで有意な差があるか」をTableau上で確かめることができるようになりました。

ここまでの内容を元にして、月次で有意な差があったかを確かめられるビューを作成できます。

表計算の設定含め、詳細はWorkbookをダウンロード頂ければと思います。

---

等分散性の検定

ところで、ここまで本当にWestとEastでOrder IDごとのSalesにおける分散が等しいのかを確かめていませんでした。

等分散性を確かめるにはF検定を使用します。

Rには簡単に等分散性を確かめる関数var.test()があるので、実際に使用してみましょう。

p-value of F_test

SCRIPT_REAL(

'var.test(

subset(.arg1,!is.na(.arg1))

, subset(.arg2,!is.na(.arg2))

, var.equal=T)$p.value'

, SUM([East Sales]), SUM([West Sales])

)

p_value < 0.05 (F_test)

//H1: true difference in means is not equal to 0

IF [p-value of F_test] < 0.05

THEN 'Diff. Var.'

ELSE 'Same Var.'

END

見てみると、結構等分散性は月によって変わるようです。

したがって等分散性に合わせてt検定の中身を変える実装が必要そうです。

ところで、t検定の等分散性仮定はt.test()関数のvar.equalパラメータで制御します。

したがってF検定の結果をT/Fで渡す計算フィールドを用意してあげれば良いということになります。

F-test Result

//H1: true difference in means is not equal to 0

IF [p-value of F_test] < 0.05

THEN 'F' //H0 is rejected

ELSE 'T' //H0 is not rejected

END

p-value of t_test（修正版）

SCRIPT_REAL(

't.test(

subset(.arg1,!is.na(.arg1))

, subset(.arg2,!is.na(.arg2))

, var.equal=' + [F-test Result] +')$p.value'

, SUM([East Sales]), SUM([West Sales])

)

ところで結論から言えば、等分散性にかかわらずt検定の結果は変わりませんでした。

上記ビューではEastとWestのOrder IDの個別カウントを追加表記しているのですが、見る限り

• 平均値の差が大きく出ていてもOrder数（サンプル数）が少ないので、結果的にt検定で有意な差が出にくい（かつ正規分布の過程がしにくい：気になる方はt検定のロバスト性と中心極限定理をお調べ下さい）。

• サンプル数が大きくても、そもそもの平均値差が小さいので差が有意でないと判断される。

上記2点が出てくるのかなと。結論としては、本質的にはEastとWestではOrderごとの平均Salesに有意な差は基本的には無さそうです。

じゃあRegion全部の差の有意性を見たいときはどうするんだ、という話は別の日に。

気になる方は分散分析（ANOVA）でお調べ下さい。近々ブログを書きます。

---

2020/7/1追記：異常値の処理について

今回の本題であったt検定の実装自体は上記で完結しているのですが、上記の結果の妥当性についてちょっと確認したいことが出来たので、ここに追記として書いていきます。

（書きながら実験しています）

まず、差がそれなりに出つつ「2データ間の平均値の差は有意でない」という結果が出ているわけですが、仮にt検定を実施する際のサンプル数を増やしたらどうなるのでしょうか？つまりサンプル数が増えれば、データの分散は対応して小さくなってくれるのでしょうか？

以下に四半期ごとにデータの集計を変えたものを記載します。

サンプル数は十分そうですが、やはり平均値の差に統計的有意さは出ません。

少し原点に立ち返り、実際のデータの分布をみてみます（例として2018年を表示）。

Box Plotと平均値を表示しています。

分布が似通っているので確かに有意さは無さそうですが、それ以上に分布が広がっていることが気になります。これは例えば顧客Segmentなどによる影響でしょうか？少し見てみましょう。

（ちなみに今回はt検定の実装という文脈から入ったので最初にデータを確認しませんでしたが、こういうことをする際は最初にデータの分布を見ないといけませんね）

あまり違いが見えませんが、それ以上に「以上に合計Salesが大きいデータ」に引っ張られているようにも見えます。

このような分布から外れたデータを異常値などと呼ぶのですが、サンプル数が決して多くはないデータにおける平均値（と分散）は異常値の影響を受けやすく、もちろんt検定の結果にも影響します。

異常値処理の方法はいくつかあるのですが、今回はTableauで実装しやすい四分位：Quartileを使用します。

つまりBox Plotの上端と下端を超えたデータを省くような実装を考えます。

Inter Quatile Range

WINDOW_PERCENTILE(SUM([Sales]),0.75)

- WINDOW_PERCENTILE(SUM([Sales]),0.25)

Outlier Label

SUM([Sales]) > WINDOW_PERCENTILE(SUM([Sales]),0.75) + 1.5 * [Inter Quatile Range]

OR

SUM([Sales]) < WINDOW_PERCENTILE(SUM([Sales]),0.25) - 1.5 * [Inter Quatile Range]

図中の赤い点を省いたらOrderごとの分布はどうなりそうでしょうか？見てみましょう。

いくつかの四半期では、もしかしたら有意差が出そうですね。

この異常値処理をt検定の式の中に組み込む方法を以降で考えていきます。

まずはベクトルサイズを確認し、実装が正しそうか確認します。

vector size (east) (outlier) (Westも同様)

SCRIPT_REAL(

'length(subset(.arg1,!is.na(.arg1)))'

, IF NOT [Outlier Label] THEN SUM([East Sales]) END

)

各ラベルの上側が異常値を含めたベクトルサイズ、下が異常値除外したベクトルサイズです。異常値の個数も確かめましたが合っていました。

Tableauの操作の順序によると表計算フィルターは表計算の後に走ってしまうので、表計算の中にIF文でのフィルタリングを実装する形になります。

このIF文込みの計算式が入力ベクトルになり、subset()の中で欠損値は除外されるので、結果的に異常値抜きの入力ベクトルがインプットになる、という算段です。

ということで、やっとt検定の式を改良できます。

F-test Result (outlier)

//H1: true difference in means is not equal to 0

IF [p-value of F_test (outlier)] < 0.05

THEN 'F' //H0 is rejected

ELSE 'T' //H0 is not rejected

END

p-value of F_test (outlier)

SCRIPT_REAL(

'var.test(

subset(.arg1,!is.na(.arg1))

, subset(.arg2,!is.na(.arg2))

, var.equal=T)$p.value'

, IF NOT [Outlier Label] THEN SUM([East Sales]) END

, IF NOT [Outlier Label] THEN SUM([West Sales]) END

)

p-value of t_test (outlier)

SCRIPT_REAL(

't.test(

subset(.arg1,!is.na(.arg1))

, subset(.arg2,!is.na(.arg2))

, var.equal=' + [F-test Result (outlier)] +')$p.value'

, IF NOT [Outlier Label] THEN SUM([East Sales]) END

, IF NOT [Outlier Label] THEN SUM([West Sales]) END

)

p_value < 0.05 (t_test) (outlier)

// H1: true difference in means is not equal to 0

IF [p-value of t_test (outlier)] < 0.05

THEN '★'

ELSE ''

END

ビューの細かい作りはWorkbookを見て頂くとして、結果を以下に記載します。

特に四半期ごとの（Segmentで分けない）平均Sales差は統計的有意性が２つの四半期で出るようになりました。もちろん統計的有意性を出すことが目的では無いのですが。

欠損値を除外する処理はやはり実装しておいた方がいいのかなと思いつつ、このあたりで追記を着地させることにします。

---

最後に

今回はt検定について取り上げました。そして追記では欠損値処理にも少々触れました。

本当は正規分布を仮定しない場合の実装についても考書こうか考えたのですが、F検定あたりでスタミナが切れました。気になる方はお調べ下さい。

実装はそんなに難しくはならないと思います。

ご質問等はTwitter、Linkedinへよろしくお願いします。それでは。