Ich möchte Dunnetts Test in Python machen

Einführung

Dies ist mein erster Beitrag. Ich hoffe, Sie werden es mit einem warmen Gefühl lesen. Hintergrund ist ein Grundlagenforscher, der in vielen Pharmaunternehmen arbeitet, daher sind statistische Kenntnisse und Programmierung nur Selbststudium ...

Nachtrag (14.02.2020)

Es gab verschiedene Druckfehler ... Es tut mir leid, wenn Sie versucht haben, es zu verschieben, und es hat nicht funktioniert ... Ich hätte es reparieren können, also danke. Dunnett's test In der biologischen Forschung möchten Sie häufig einen ** Student-t-Test durchführen und einen signifikanten Unterschied zeigen ** ([Der Ausdruck "signifikanter Unterschied" sollte beseitigt werden]) (https://www.nature.com/articles/d41586-). 019-00857-9) ist eine weitere Gelegenheit ...).

Bei der Bewertung zwischen mehreren Gruppen besteht ein Problem der Multiplizität, sodass es viele Möglichkeiten gibt, den Tukey-HSD-Test zu verwenden. Der Dunnett-Test ist auch einer der Mehrgruppen-Paarvergleichstests. Diese Methode wird verwendet, wenn Sie nicht alle Paarvergleiche durchführen müssen und ** Kontrollgruppe (oder Kontrollgruppe) mit anderen Gruppen ** vergleichen möchten. Ich benutze es ziemlich oft, weil es ** eine höhere Erkennungsleistung ** hat als der Vergleich aller Gruppen.

Auf Wikipedia

In der Statistik ist der Dannett-Test (Dunnett-Test) eines der zahlreichen Vergleichsverfahren [1]. Entwickelt vom kanadischen Statistiker Charles Danette, um jede der vielen Behandlungsgruppen mit einer einzigen Kontrollgruppe zu vergleichen. Mehrfachvergleiche mit der Kontrollgruppe werden auch als Viele-zu-Eins-Vergleiche bezeichnet. Wikipedia

Es wurde beschrieben als.

Ich verstehe, dass Sie ** verwenden können, wenn die Bevölkerung normal und gleichmäßig verteilt ist **.

Ich möchte alles mit Python machen!

Ich habe Excel verwendet, um die Daten zu analysieren, aber jetzt mache ich eine Datenanalyse mit Python (mit Colaboratory), weil ich Python persönlich üben möchte. Für die statistische Analyse brachte ich jedoch Daten ein und verwendete R mit EZR und SAS mit EXSUS.

Bitte lesen Sie den Artikel von @ issakuss (Python Statistical Techniques-Statistical Analysis Against Python-), in dem ** [pingwouin](https: / Ich habe ein Paket namens /pingouin-stats.org/api.html#)** gefunden.

Sie können die meisten Testmethoden verwenden! !! !!

Es ist wundervoll! !! Danke an den Schöpfer. Damit schien es, als könnte ich allein von Python leben.

Jedoch…

Es gibt keinen Dunnett-Test in Pingouin ...

Lassen Sie es uns also während des Studiums selbst implementieren! Ich habe mich dazu entschlossen. (~~ Nun, wenn Sie danach suchen, gibt es möglicherweise irgendwo ein Paket ~~)

Politik

Ich musste die Theorie kennen, um sie selbst umzusetzen, so leicht! !! studiert.

Die Methode zur Berechnung der Teststatistik T ähnelt Tukey. Unter der Annahme, dass es k Gruppen von 1 bis i gibt, ist die ** Teststatistik T ** 1 mit der Kontrollgruppe als 1.

jedoch

• $ μ $: Durchschnitt
• $ n $: Anzahl der Proben
• $ σ ^ 2 $: Kovarianz (Ich denke, es ist allgemein bekannt, weil ich die Definition von Zeichen oft nicht verstanden habe, aber lassen Sie mich sie beschreiben).

Zu diesem Zeitpunkt ist ** Freiheit v ** $ v = \sum_{i=1}^{k}n_i - k$ (Das heißt (Gesamtzahl der Proben) - (Gesamtzahl der Gruppen)).

Im Vergleich zu dieser Teststatistik T ...

Das ist Danettes Tisch! !! !! ………?

Tabelle? Es war wie (lacht) Apropos, es gab viele Tabellen auf der Rückseite der Lehrbücher der Universität ... Die eigentliche Tabelle ist [this](https://www.weblio.jp/content/%E3%83%80%E3%83%8D% E3% 83% 83% E3% 83% 88% E3% 81% AE% E8% A1% A8). In Pingouin schien es einen Tukey-Tisch einzulegen und sich darauf zu beziehen, aber ich wollte keinen Tisch einlegen ...

SAS macht etwas Erstaunliches

Es scheint, dass SAS eine ** probmc-Funktion ** verwendet und keine Tabelle verwendet. Also habe ich beschlossen, diese Funktion zu implementieren (die meiste Zeit hatte ich Schwierigkeiten). Zum Glück wurde die Formel in hier aufgeführt, sodass ich sie so implementieren werde, wie sie ist. Ich werde.

Die probmc-Funktion arbeitet mit den folgenden Argumenten (einfach zu schreiben, ich kenne SAS nicht).

probmc(distribution, q, prob, df, nparms <, parameters>)

ist. Was ist λ, aber es ist ein Wert, der aus der Anzahl der zu vergleichenden Proben berechnet wird $ \lambda_{1i} = \sqrt{\frac{n_i}{n_1+n_i}} $ Es wird sein. Wenn die Anzahl der Stichproben in jeder Gruppe gleich ist, werden keine Parameter übergeben. Bei der Berechnung können Sie jedoch andere Dinge tun, wenn Sie diejenige implementieren, die am kompliziertesten erscheint. ** Wenn die Anzahl der Stichproben in jeder Gruppe unterschiedlich ist ** Implementieren.

(Ich werde später bestätigen, ob die Formel gleich ist, wenn die Anzahl der Proben mit der Formel übereinstimmt, wenn die Anzahl der Proben in jeder Gruppe unterschiedlich ist.)

Probmc-Funktion (beide Seiten dunnett mit unterschiedlicher Anzahl von Proben in jeder Gruppe)

Das von der probmc-Funktion erhaltene prob sieht folgendermaßen aus:

prob = probmc("dunnet2", T, ., v, k-1, <lambda1,lambda2...>)

 Und die Formel, die darin funktioniert, lautet wie folgt.

```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)

Wenn Sie diese Formel implementieren können, sind Sie fertig! (Ich hatte nicht die Energie, den Inhalt zu verstehen ...)

Hier $ \ phi (x) $: Wahrscheinlichkeitsdichtefunktion (pdf) $ \ Phi (x) $: Funktion der kumulativen Dichte (cdf) Es wird sein. Diese können mit scipy und numpy implementiert werden.

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])

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

Es wird sein.

Außerdem wird $ d \ mu_v (x) $ in SAS wie folgt definiert:

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

Hier $ \ Gamma (x) $: Gammafunktion Und es gibt auch eine Funktion in scipy (auch in Mathe vorhanden).

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))

** Ein Problem hier ... ** Das Integrationsintervall der ursprünglichen Gleichung beträgt $\int_{0}^{∞}d\mu_v(x)$ Es war, aber das $\int_{0}^{∞}d\mu_v(x)=\int_{?}^{?}g_v(x)dx $ Betrachten Sie die Grafik, wenn Sie versuchen, in ...

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): #ändern 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()

Hmm ...? ** Ist es auf der Minus-Seite gleichmäßig umgekehrt? ** Dies war ein Problem beim Ersetzen des Integrationsintervalls ... Es sieht aus wie eine Sigmoid-Funktion, daher verwende ich sie wahrscheinlich auf diese Weise. Leider konnte ich die Bedeutung der Formel nicht verstehen, deshalb habe ich mich entschlossen, vorerst mit dem Bereich von ** x von 0 bis ∞ ** fortzufahren.

Wenn man die Genauigkeit betrachtet, wenn der Wert am größten ist, um die Anzahl der Berechnungen zu verringern, scheint $ \ mu_9 (3) = 0,99999999999998976 $ ausreichend zu sein, sodass der Bereich von x in der Berechnung 0 → 5 beträgt. (Übrigens $ \ mu_9 (2) = 0,99999603534120194 $)

Bereit!

Also die ursprüngliche Formel,

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)

Wird transformiert $ prob=\int_{0}^{∞}\int_{-∞}^{∞}f_{T\lambda}(x,y)g_v(x)dydx $ ist geworden. Implementieren Sie dies

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]

Schauen wir uns nun den Graphen von $ f_ {T \ lambda} (x, y) g_v (x) $ an, um den Rechenaufwand zu reduzieren.

T, v = 10,9 #Stellen Sie den Wert entsprechend ein
lambda_param =[0.5,0.4] #Stellen Sie den Wert entsprechend ein

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()

Es sieht aus wie eine gleichmäßige Funktion! Mit anderen Worten

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

jedoch

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]

Es wird sein.

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]

Komplett! !! ...?

Ich denke, es wird wie folgt zusammengefasst. Ich wollte es ordentlich mit Klassen und so weiter organisieren, aber es tut mir leid, dass es schmutzig ist, weil ich es nicht gut machen kann.

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]

Ich habe es tatsächlich benutzt

Erstellen Sie entsprechende Daten und probieren Sie sie aus. Betrachten Sie das Muster, ob Arzneimittel A und Arzneimittel B wirksame Werte für die Kontrollgruppe ergeben. Zuerst von der Datenerstellung

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)

Wir werden zu diesen Daten gehen. Zunächst aus dem Data Moulding. Es ist notwendig, die Teststatistik T und λ zu berechnen. Ich würde gerne eine ähnliche Methode anwenden, um sie später auf Pigouin anzuwenden.

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'] #Freiheitsgrad
  ng = aov.at[0, 'DF'] + 1 #Gesamtzahl
  grp = data.groupby(between)[dv]
  n_sub = grp.count()
  control_index = n_sub.index.get_loc(control) #Index der Kontrollgruppe abrufen
  n = n_sub.values #Anzahl von Beispielen
  gmeans = grp.mean().values#Jeder Durchschnitt
  gvar = aov.at[1, 'MS'] / n #Jede Distribution

  vs_g = np.delete(np.arange(ng),control_index) #Holen Sie sich einen anderen Index als die Kontrollgruppe

  #Paarweise Kombinationen (Teststatistik T finden)
  mn = np.abs(gmeans[control_index]- gmeans[vs_g])
  se = np.sqrt(gvar[control_index] + gvar[vs_g]) #Die Formel ist etwas anders, aber weil sie eine gleichmäßige Verteilung voraussetzt
  tval = mn / se

  #Finde 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")

ist geworden. Ist es wirklich richtig ...?

Probieren Sie es mit EZR

Ich bin völlig müde ... Schließlich möchte ich die gleiche Analyse mit R unter Verwendung von EZR durchführen. R kann nicht geschrieben werden, bitte verzeihen Sie mir nur die Ausgabe.

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

Das? ?? ?? Ist es ein bisschen anders? ??

Ich bin müde, also werde ich mein Bestes geben, um es noch einmal zu überprüfen. Es scheint, dass R die Monte-Carlo-Methode verwendet, also denke ich, dass die Werte unterschiedlich sind, aber 0,005 ist zu groß, nicht wahr?

Vielen Dank für das bisherige Anschauen. Wie oben erwähnt, haben Amateure ihr Bestes gegeben.

Verweise

[Mehrfachvergleich (Teil 2-1: Spezifische Methode (i)), 3. Armitage-Studiensitzung: Material 2, Masaaki Doi](http://www012.upp.so-net.ne.jp/doi/ Biostat / CT39 / multiple_comparison_2_1.pdf)

Zum Kopieren

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'] #Freiheitsgrad
  ng = aov.at[0, 'DF'] + 1 #Gesamtzahl
  grp = data.groupby(between)[dv]
  n_sub = grp.count()
  control_index = n_sub.index.get_loc(control) #Index der Kontrollgruppe abrufen
  n = n_sub.values #Anzahl von Beispielen
  gmeans = grp.mean().values#Jeder Durchschnitt
  gvar = aov.at[1, 'MS'] / n #Jede Distribution

  vs_g = np.delete(np.arange(ng),control_index) #Holen Sie sich einen anderen Index als die Kontrollgruppe

  #Paarweise Kombinationen (Teststatistik T finden)
  mn = np.abs(gmeans[control_index]- gmeans[vs_g])
  se = np.sqrt(gvar[control_index] + gvar[vs_g]) #Die Formel ist etwas anders, aber weil sie eine gleichmäßige Verteilung voraussetzt
  tval = mn / se

  #Finde 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