LoginSignup
0
0

More than 1 year has passed since last update.

@mosomoso_1910

pythonでDunnett’s test

Last updated at Posted at 2022-08-10

はじめに

本記事は以前の記事(https://qiita.com/mosomoso_1910/items/101783b4e59446d00d42)の改訂版となります。
（凡ミスで計算違いしてました。。。あと無限大の積分でもそこまで計算速度は変わらなさそうなので式のままにしています）

Backgroundはそこらへんの製薬企業で働いている基礎研究者なので統計学的知識もプログラミングも独学のみです…
生暖かい気持ちで読んでいただけると幸いです。

Dunnett's test

生物系研究ではよくStudent's t testを行い、有意差を示したいということがありますね（有意差という表現をなくすべきというのはまた別の機会に…）。

多群間で評価を行う場合は多重性の問題が起こるため、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で抗う統計解析～）を読ましていただき、そこで**pingwouin**というパッケージを知りました。

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

すばらしいですね！！
作成者に感謝です。これで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と比較するのが…

ダネットの表です！！！………？

表？って感じでした（笑）
そういえば大学の教科書の裏にいっぱい表が乗ってたなあ…実際の表はこんなんです。pingouinではTukeyの表を入れて参照する形で行っているようでしたが、表を入れるのは嫌でした…

SASはなんかすごいことしてる

SASではprobmc関数というものを用いていて、表を使用していないそうです。なのでこの関数を実装することにしました（苦労したのはほとんどこの部分）。
こちらにありがたいことに式が載っておりましたので、そのまま実装していきます。

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...>)
そしてこの内部で動いている式は以下のようになります。

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を使用することで実装可能です。

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で次のように定義されています。

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でも存在）。

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 $$
に変換しようとしたときにグラフを見てみると…

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

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

準備完了！

ということで元の式、

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 $$
となりました。
これを実装して

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

完成！！…？

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

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 integrate.dblquad(f_g,-np.inf,np.inf,lambda x:0,lambda x:np.inf)[0]

実際に使ってみた

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

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に載せることを考えて同じような手法を取りたいと思っています。

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 = np.sqrt(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.992065
control drug_B 10.30 14.67 4.369 1.993 2.193 0.083971

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

EZRでやってみる

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


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

ということで無事一致しました！！！！やったーーー！！！

参考文献

多重比較について（その２－１：具体的手法（i））, 第３回 Armitage 勉強会：資料２, 土居正明

コピペ用

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 = np.sqrt(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

最後に

問題として、実行時間が非常に長いです。
例題では25秒かかっています。
ここを改善していただける方、よろしくお願いいたします。
scikit-posthocsのissuesに投げておいたので、もしかしたら入れてくれるかもしれません。
(https://github.com/maximtrp/scikit-posthocs/issues/24)

0
0
1

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
0
0