PythonでDunnett's testがしたい

#はじめに
初めての投稿となります。
生暖かい気持ちで読んでいただけると幸いです。
Backgroundはそこらへんの製薬企業で働いている基礎研究者なので統計学的知識もプログラミングも独学のみです…
###追記(2020/02/14)
色々と誤植がありました…万が一動かそうとして動かなかった方は申し訳ありません…
修正できたはずですのでよろしくお願いします。
##Dunnett's test
生物系研究ではよく**Student's t testを行い、有意差を示したい**ということがありますね（[有意差という表現をなくすべき](https://www.nature.com/articles/d41586-019-00857-9)というのはまた別の機会に…）。

多群間で評価を行う場合は多重性の問題が起こるため、Tukey-HSD testなどを使う機会が多いです。
Dunnett's testも多群の対比較の検定のひとつで、すべての対比較を行う必要がなく、**コントロール群（や対照群）と他の群らで比較したい場合**に用いる手法です。すべての群の対比較をするよりも**検出力が高くなる**ので結構使います。

Wikipediaでは
>統計学において、ダネットの検定（ダネットのけんてい、英: Dunnett's test）は、多重比較手順の一つである[1]。カナダの統計学者チャールズ・ダネットによって、多くの処理群のそれぞれと単一の対照群を比較するために開発された。対照群に対する多重比較は多対一（many-to-one）比較とも呼ばれる。Wikipedia

と記載されています。

使用していいのは、**母集団が正規分布かつ等分散の時**、と理解しています。

##pythonで全部したい！
データの解析を行う上でexcelを使っていたのですが、個人的にpythonを練習したいというのもあり今はデータ解析をpythonで行っています（Colaboratoryを用いて）。ただし、統計解析についてはデータを持ち込んでEZRでRや、EXSUSでSASを使用したりしていました


@issakussさんの記事（[Python統計テクニック集 ～Pythonで抗う統計解析～](https://qiita.com/issakuss/items/3a3fe667a8863add99a3)）を読ましていただき、そこで**[pingwouin](https://pingouin-stats.org/api.html#)**というパッケージを知りました。

###なんとほとんどの検定法が使える！！！

すばらしいですね！！
作成者に感謝です。これでpythonだけで生きていけそうでした。

ところが…
###pingouinにDunnett's testがない…

ということになり、ならもう勉強ついでに自分で実装してみよう！ということにしました。
（~~まあ探せばどっかにパッケージあるのかもしれないけど~~）

#方針
自分で実装する上で理論を知る必要がありましたので、軽く！！勉強しました。

検定統計量Tを算出する方法はTukeyと同じようです。
群が1~iまでのk個あるとすると、対照群を1として**検定統計量T**は


$$ T{1i}=\frac{|μ_1-μ_i|}{\sqrt{(\frac{1}{n_1}+\frac{1}{n_i})\sigma^2}}$$
ただし
- $μ$: 平均
- $n$: サンプル数
- $σ^2$: 共分散
とします（私自身文字の定義が分からないのが多かったので常識かと思いますが記述させてください）。

このとき**自由度v**は
$$ v = \sum_{i=1}^{k}n_i - k$$ 
となります（つまり（サンプル数の総和）ー（全群数））。

この検定統計量Tと比較するのが…
###ダネットの表です！！！………？
表？って感じでした（笑）
そういえば大学の教科書の裏にいっぱい表が乗ってたなあ…実際の表は[こんなん](https://www.weblio.jp/content/%E3%83%80%E3%83%8D%E3%83%83%E3%83%88%E3%81%AE%E8%A1%A8)です。pingouinではTukeyの表を入れて参照する形で行っているようでしたが、表を入れるのは嫌でした…

##SASはなんかすごいことしてる
SASでは**probmc関数**というものを用いていて、表を使用していないそうです。なのでこの関数を実装することにしました（苦労したのはほとんどこの部分）。
[こちら](https://support.sas.com/documentation/cdl_alternate/ja/lefunctionsref/67960/HTML/default/p06idazwtugz0an1rlks4swmutms.htm)にありがたいことに式が載っておりましたので、そのまま実装していきます。

probmc関数は以下のような引数で動きます（簡単に書きます、SASはわかりません）。

```PROBMC(distribution, q, prob, df, nparms <, parameters>)```

|引数|内容| 
|:---|---|
|distribution|Dunnettが片側か両側かを選びます。Dunnett以外の分布にも使えるようです。|
|q, prob|qは分位点（ここでは検定統計量Tを入れる）、probは確率です。probは確率変数がq未満となる確率です。したがって、p値はは1– probのように計算できます。たとえば、有意水準5%の臨界値を計算するには、prob= 0.95と設定します。どちらかを欠損させて引数とし、**欠損した引数が返り値**となります。今回はp値が知りたいのでprobを欠損させることになります。|
|df|自由度v|
|nparms|ダネットの時は全群数k-1|
|parameters|λ12,λ13,...λ1i|

です。
λってなんだってことになりますが、比較するサンプル数から計算される値で
$$ \lambda_{1i} = \sqrt{\frac{n_i}{n_1+n_i}} $$
となります。各群のサンプル数が同じ場合はparametersは渡されないのですが、一番ややこしそうなやつを実装すれば他のもできるだろうという算段で、**各群のサンプル数が異なる場合**について実装します。

（各群のサンプル数が異なる場合の式に対してサンプル数が同じになれば式が同じになるのかは後日確かめます）

##probmc関数(各群のサンプル数が異なる両側dunnett）
probmc関数によって得られるprobは次のようになります
```prob = probmc("dunnet2", T, ., v, k-1, <lambda1,lambda2...>)```
そしてこの内部で動いている式は以下のようになります。

```math
prob = \int_{0}^{∞}\int_{-∞}^{∞} ϕ(y)∏_{i=1}^{k-1}\biggl[Φ\biggl(\frac{\lambda_iy+Tx}{\sqrt{1-\lambda_i^2}}\biggr)- Φ\biggl(\frac{\lambda_iy-Tx}{\sqrt{1-\lambda_i^2}}\biggr)\biggr]dy d\mu_v(x)
```

この式が実装できれば完成です！
（中身の理解をする元気はありませんでした…）

ここで
$\phi(x)$: 正規分布の確率密度関数（pdf: Probability density function）
$\Phi(x)$: 累積分布関数（cdf: Cumulative density function）
となります。
これらはscipy, numpyを使用することで実装可能です。

 ```python
import numpy as np
from scipy.stats import norm

lambda_params = [lambda_1,lambda_2]
T=T


def pdf(x):
  return norm.pdf(x)

def cdf(x):
  return norm.cdf(x)

def f(x,y):
  return pdf(y)*np.prod([cdf((lambda_i*y+T*x)/(np.sqrt(1-lambda_i**2)))-cdf((lambda_i*y-T*x)/(np.sqrt(1-lambda_i**2))) for lambda_i in lambda_params])
```

これで
$$ prob=\int_{0}^{∞}\int_{-∞}^{∞}f_{T\lambda}(x,y)dyd\mu_v(x) $$

となります。

また、$d\mu_v(x)$はSASで次のように定義されています。

```math
d\mu_v(x)=\frac{ν^{\frac{v}{2}}}{\Gamma(\frac{v}{2})2^{\frac{v}{2}-1}}x^{v-1}exp(-\frac{vx^2}{2})dx
```

ここで
$\Gamma(x)$: ガンマ関数
であり、これもscipyで関数があります(mathでも存在）。

```python
from scipy.special import gamma

v=v

def du(v,x):
  return (v**(v/2)*x**(v-1)*np.exp(-v*x**2/2))/(gamma(v/2)*2**(v/2-1))
```

**ここで一つ問題が…**
元の式の積分区間が
$$\int_{0}^{∞}d\mu_v(x)$$ 
だったのですが、これを
$$\int_{0}^{∞}d\mu_v(x)=\int_{?}^{?}g_v(x)dx $$
に変換しようとしたときにグラフを見てみると…

```python
from scipy import integrate
import matplotlib.pyplot as plt

x=np.arange(-3,3,0.01)

plt.figure()
for v in range(9,30,3): #vを変化させる

  def duv(x):
    return du(v,x)

  y = [integrate.quad(duv,0,xi)[0] for xi in x]

  plt.plot(x,y,label=i)
plt.xlabel("x")
plt.ylabel("$\mu_v(x)$")
plt.legend()
plt.show()
```
![ダウンロード (1).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/330449/a3cc2eb3-fd72-c51b-4208-77155906ac6d.png)

ん…？**偶奇でマイナス側で反転する？**これは積分区間を置換する上で困った…形を見る感じシグモイド関数っぽいのでおそらくそういう形で使用しているんでしょう。
残念ながら式の意味まで理解ができなかったので、**xの範囲は0→∞**でとりあえず進めることにします。

また計算を少なくするため最も値が大きくなる時の精度を見てみると、$\mu_9(3)=0.9999999999998976$で十分そうなので計算上はxの範囲は0→5ということにします。
(ちなみに$\mu_9(2)=0.9999603534120194$）

###準備完了！

ということで元の式、

```math
prob = \int_{0}^{∞}\int_{-∞}^{∞} ϕ(y)∏_{i=1}^{k-1}\biggl[Φ\biggl(\frac{\lambda_iy+Tx}{\sqrt{1-\lambda_i^2}}\biggr)- Φ\biggl(\frac{\lambda_iy-Tx}{\sqrt{1-\lambda_i^2}}\biggr)\biggr]dy d\mu_v(x)
```

は変形して
$$ prob=\int_{0}^{∞}\int_{-∞}^{∞}f_{T\lambda}(x,y)g_v(x)dydx $$
となりました。
これを実装して

```python
def f_g(x,y):
  return f(x,y) * duv(x)

def dunnett_prob(T,v,lambda_params):
  T=T
  v=v
  lambda_params=lambda_params
  return integrate.dblquad(f_g,0,np.inf,lambda x:-np.inf,lambda x:np.inf)[0]
```

ここで計算量を減らすために$f_{T\lambda}(x,y)g_v(x)$のグラフをみてみます。

```python
T, v = 10,9 #適当に値を設定
lambda_param =[0.5,0.4] #適当に値を設定

y = np.arange(-10,10,0.1)

plt.figure()
for x in np.arange(0,2,0.5):
  z = [f_g(x,yi) for yi in y]
  plt.plot(y,z,label=x)
plt.legend()
plt.xlabel("y")
plt.ylabel("$f(x,y)g(x)$")
#plt.ylim(0,0.0001)
plt.show()
```
![ダウンロード (2).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/330449/78318ad8-a2c0-a615-5332-e456d1e3a370.png)

どうやら偶関数のようですね！つまり

```math
prob ≒ 2 \int_{0}^{5}\int_{0}^{∞}f_{T\lambda}(x,y)g_v(x)dx
```

ただし

```math
f_{T\lambda}(x,y) = ϕ(y)∏_{i=1}^{k-1}\biggl[Φ\biggl(\frac{\lambda_iy+Tx}{\sqrt{1-\lambda_i^2}}\biggr)- Φ\biggl(\frac{\lambda_iy-Tx}{\sqrt{1-\lambda_i^2}}\biggr)\biggr]
```
$$g_v(x)=\frac{ν^{\frac{v}{2}}}{\Gamma(\frac{v}{2})2^{\frac{v}{2}-1}}x^{v-1}exp(-\frac{vx^2}{2})$$
となります。

```python
def dunnett_prob(T,v,lambda_params):
  T=T
  v=v
  lambda_params=lambda_params
  return 2*integrate.dblquad(f_g,0,3,lambda x:0,lambda x:np.inf)[0]

```

##完成！！…？
というわけでまとめると以下のようになるかと思います。
クラスやらなんやらを使って綺麗にまとめたかったのですが、うまくできないので汚くて申し訳ありません。

```python:dunnett.py
import numpy as np
from scipy.stats import norm
from scipy import integrate
from scipy.special import gamma

def dunnett_prob(T,v,lambda_params):

  def pdf(x):
    return norm.pdf(x)

  def cdf(x):
    return norm.cdf(x)

  def f(x,y):
    return pdf(y)*np.prod([cdf((lambda_i*y+T*x)/(np.sqrt(1-lambda_i**2)))-cdf((lambda_i*y-T*x)/(np.sqrt(1-lambda_i**2))) for lambda_i in lambda_params])

  def duv(x):
    return (v**(v/2)*x**(v-1)*np.exp(-v*x**2/2))/(gamma(v/2)*2**(v/2-1))

  def f_g(x,y):
    return f(x,y) * duv(x)

  return 2*integrate.dblquad(f_g,0,5,lambda x:0,lambda x:np.inf)[0]
```

##実際に使ってみた
適当なデータを作成してみて、実際に試してみます。
薬剤Ａと薬剤Ｂがcontrol群にたいして有効な値をだしているかというパターンを考えます。
まずデータ作成から

```python
import pandas as pd

control = np.random.normal(10,scale=2,size=4)
drug_A = np.random.normal(10.5,scale=2,size=6)
drug_B = np.random.normal(14,scale=2,size=5)


df_c = pd.DataFrame(data=control,columns=["value"])
df_c["group"] = "control"
df_A = pd.DataFrame(data=drug_A,columns=["value"])
df_A["group"] = "drug_A"
df_B = pd.DataFrame(data=drug_B,columns=["value"])
df_B["group"] = "drug_B"
df = pd.concat([df_c,df_A,df_B],ignore_index=True)
```

|value | group |
|-----:|-------|
|11.389|control|
| 7.154|control|
| 7.932|control|
|14.729|control|
|10.507|drug_A |
| 6.578|drug_A |
| 6.580|drug_A |
|13.098|drug_A |
|10.131|drug_A |
|13.748|drug_A |
|15.245|drug_B |
|14.078|drug_B |
|17.533|drug_B |
|11.082|drug_B |
|15.413|drug_B |


このデータに対して行っていきます。
まずはデータの成型から。検定統計量Ｔとλを算出する必要があります。のちのちpigouinに載せることを考えて同じような手法を取りたいと思っています。

```python
import pingouin as pg

def dunnett_two_side_unbalanced(data, dv, between,control):

  #First compute the ANOVA
  aov = pg.anova(dv=dv, data=data, between=between, detailed=True)
  v = aov.at[1, 'DF'] #自由度
  ng = aov.at[0, 'DF'] + 1 #全群数
  grp = data.groupby(between)[dv]
  n_sub = grp.count()
  control_index = n_sub.index.get_loc(control) #対照群のindexを取得
  n = n_sub.values #サンプル数
  gmeans = grp.mean().values#各平均
  gvar = aov.at[1, 'MS'] / n #各分散

  vs_g = np.delete(np.arange(ng),control_index) #対照群以外のindexを取得

  # Pairwise combinations（検定統計量Ｔを求める）
  mn = np.abs(gmeans[control_index]- gmeans[vs_g])
  se = np.sqrt(gvar[control_index] + gvar[vs_g]) #式は少し違うが等分散を仮定しているため
  tval = mn / se

  #lambdaを求める
  lambda_params = n[vs_g]/(n[control_index]+n[vs_g])

  pval = [1-dunnett_prob(t,v,lambda_params) for t in tval]

  stats = pd.DataFrame({
                      'A': [control]*len(vs_g),
                      'B': np.unique(data[between])[vs_g],
                      'mean(A)': np.round(gmeans[control_index], 3)*len(vs_g),
                      'mean(B)': np.round(gmeans[vs_g], 3),
                      'diff': np.round(mn, 3),
                      'se': np.round(se, 3),
                      'T': np.round(tval, 3),
                      'p-dunnett': pval
                      })
  return stats

dunnett_two_side_unbalanced(df,"value","group","control")

```

|   A   |  B   |mean(A)|mean(B)|diff | se  |  T  |p-dunnett|
|-------|------|------:|------:|----:|----:|----:|--------:|
|control|drug_A|  10.30|  10.11|0.194|1.917|0.101|  0.99312|
|control|drug_B|  10.30|  14.67|4.369|1.993|2.193|  0.08992|

となりました。果たしてあっているのか…？

##EZRでやってみる
もう完全に疲れてきました…最後にEZRを用いてＲで同じ解析を行ってみたいと思います。Ｒは書けないので出力だけでご勘弁ください。

```R

Multiple Comparisons of Means: Dunnett Contrasts


Fit: aov(formula = value ~ group.factor, data = Dataset)

Linear Hypotheses:
                      Estimate Std. Error t value Pr(>|t|)  
drug_A - control == 0  -0.1939     1.9174  -0.101    0.992  
drug_B - control == 0   4.3693     1.9926   2.193    0.084
```

###あれ？？？ちょっと違うぞ？？
もう疲れたのでまた今度頑張って見直してみます。Ｒではモンテカルロ法を用いているらしいので値は異なるかと思いますが、0.005は大きすぎますよね？

ここまで見ていただいた方がもしいらっしゃればありがとうございます。
以上素人が頑張ってみたでした。



##参考文献
[多重比較について（その２－１：具体的手法（i））, 第３回 Armitage 勉強会：資料２, 土居正明](http://www012.upp.so-net.ne.jp/doi/biostat/CT39/multiple_comparison_2_1.pdf)


##コピペ用
```python
import numpy as np
from scipy.stats import norm
from scipy import integrate
from scipy.special import gamma
import pingouin as pg


def dunnett_two_side_unbalanced(data, dv, between,control):

  def dunnett_prob(T,v,lambda_params):

    def pdf(x):
      return norm.pdf(x)

    def cdf(x):
      return norm.cdf(x)

    def f(x,y):
      return pdf(y)*np.prod([cdf((lambda_i*y+T*x)/(np.sqrt(1-lambda_i**2)))-cdf((lambda_i*y-T*x)/(np.sqrt(1-lambda_i**2))) for lambda_i in lambda_params])

    def duv(x):
      return (v**(v/2)*x**(v-1)*np.exp(-v*x**2/2))/(gamma(v/2)*2**(v/2-1))

    def f_g(x,y):
      return f(x,y) * duv(x)

    return 2*integrate.dblquad(f_g,0,5,lambda x:0,lambda x:np.inf)[0]

  #First compute the ANOVA
  aov = pg.anova(dv=dv, data=data, between=between, detailed=True)
  v = aov.at[1, 'DF'] #自由度
  ng = aov.at[0, 'DF'] + 1 #全群数
  grp = data.groupby(between)[dv]
  n_sub = grp.count()
  control_index = n_sub.index.get_loc(control) #対照群のindexを取得
  n = n_sub.values #サンプル数
  gmeans = grp.mean().values#各平均
  gvar = aov.at[1, 'MS'] / n #各分散

  vs_g = np.delete(np.arange(ng),control_index) #対照群以外のindexを取得

  # Pairwise combinations（検定統計量Ｔを求める）
  mn = np.abs(gmeans[control_index]- gmeans[vs_g])
  se = np.sqrt(gvar[control_index] + gvar[vs_g]) #式は少し違うが等分散を仮定しているため
  tval = mn / se

  #lambdaを求める
  lambda_params = n[vs_g]/(n[control_index]+n[vs_g])

  pval = [1-dunnett_prob(t,v,lambda_params) for t in tval]

  stats = pd.DataFrame({
                      'A': [control]*len(vs_g),
                      'B': np.unique(data[between])[vs_g],
                      'mean(A)': np.round(gmeans[control_index], 3)*len(vs_g),
                      'mean(B)': np.round(gmeans[vs_g], 3),
                      'diff': np.round(mn, 3),
                      'se': np.round(se, 3),
                      'T': np.round(tval, 3),
                      'p-dunnett': pval
                      })
  return stats
```