Je veux faire le test de Dunnett en Python

introduction

Ceci est mon premier message posté. J'espère que vous le lirez avec un sentiment chaleureux. Le fond est un chercheur de base travaillant dans de nombreuses sociétés pharmaceutiques, donc les connaissances statistiques et la programmation ne sont que des auto-études ...

Postscript (14/02/2020)

Il y a eu plusieurs erreurs d'impression ... Je suis désolé si vous avez essayé de le déplacer et que cela n'a pas fonctionné ... J'aurais dû être en mesure de le réparer, alors merci. Dunnett's test En recherche biologique, vous voulez souvent faire ** le test t de Student et montrer une différence significative ** (L'expression «différence significative» doit être éliminée. 019-00857-9) est une autre opportunité ...).

Lors de l'évaluation entre plusieurs groupes, il y a un problème de multiplicité, il existe donc de nombreuses possibilités d'utiliser le test Tukey-HSD. Le test de Dunnett est également l'un des tests de comparaison de paires multi-groupes, et c'est une méthode utilisée lorsque vous n'avez pas besoin d'effectuer toutes les comparaisons de paires et que vous souhaitez comparer ** le groupe témoin (ou le groupe témoin) avec d'autres groupes **. Je l'utilise beaucoup car il ** a une puissance de détection plus élevée ** que de comparer tous les groupes.

Sur Wikipedia

En statistique, le test de Dannett (test de Dunnett) est l'une des multiples procédures de comparaison [1]. Développé par le statisticien canadien Charles Danette pour comparer chacun des nombreux groupes de traitement avec un seul groupe témoin. Les comparaisons multiples avec le groupe témoin sont également appelées comparaisons plusieurs à un. Wikipédia

Il a été décrit comme.

Je comprends que vous pouvez utiliser ** lorsque la population est normalement répartie et uniformément répartie **.

Je veux tout faire avec python!

J'ai utilisé Excel pour analyser les données, mais maintenant je fais une analyse de données avec python (en utilisant Colaboratory) parce que je veux pratiquer python personnellement. Cependant, pour l'analyse statistique, j'ai apporté des données et utilisé R avec EZR et SAS avec EXSUS.

Veuillez lire l'article de @ issakuss (Python Statistical Techniques-Statistical Analysis Against Python-), où ** [pingwouin](https: / J'ai trouvé un paquet appelé /pingouin-stats.org/api.html#)**.

Vous pouvez utiliser la plupart des méthodes de test! !! !!

C'est merveilleux! !! Merci au créateur. Avec cela, il semblait que je pourrais vivre de python seul.

Pourtant…

Il n'y a pas de test de Dunnett dans le pingouin ...

Alors, implémentons-le nous-mêmes tout en étudiant! J'ai décidé de faire ça. (~~ Eh bien, si vous le cherchez, il peut y avoir un paquet quelque part ~~)

politique

Je devais connaître la théorie pour la mettre en œuvre moi-même, si légèrement! !! étudié.

La méthode de calcul de la statistique de test T est similaire à celle de Tukey. En supposant qu'il existe k groupes de 1 à i, la ** statistique de test T ** est 1 avec le groupe de contrôle 1.

pourtant

• $ μ $: Moyenne
• $ n $: Nombre d'échantillons
• $ σ ^ 2 $: Covariance (Je pense que c'est de notoriété publique car je n'ai souvent pas compris la définition des personnages, mais laissez-moi la décrire).

En ce moment, ** la liberté v ** est $ v = \sum_{i=1}^{k}n_i - k$ (Autrement dit, (nombre total d'échantillons) - (nombre total de groupes)).

Par rapport à cette statistique de test T ...

La table de Danette! !! !! ………?

table? C'était comme (rires) En parlant de cela, il y avait beaucoup de tableaux au dos des manuels universitaires ... Le tableau actuel est [this](https://www.weblio.jp/content/%E3%83%80%E3%83%8D%] E3% 83% 83% E3% 83% 88% E3% 81% AE% E8% A1% A8). En pingouin, il semblait mettre dans une table Tukey et y faire référence, mais je n'aimais pas mettre dans une table ...

SAS fait quelque chose d'extraordinaire

Il semble que SAS utilise une ** fonction probmc ** et n'utilise pas de table. J'ai donc décidé d'implémenter cette fonction (la plupart du temps j'ai eu du mal). Heureusement, la formule a été répertoriée dans ici, donc je vais l'implémenter tel quel. Je vais.

La fonction probmc fonctionne avec les arguments suivants (simple à écrire, je ne connais pas SAS).

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

est. Qu'est-ce que λ, mais c'est une valeur calculée à partir du nombre d'échantillons à comparer $ \lambda_{1i} = \sqrt{\frac{n_i}{n_1+n_i}} $ Ce sera. Si le nombre d'échantillons dans chaque groupe est le même, les paramètres ne seront pas passés, mais dans le calcul que si vous implémentez celui qui semble être le plus compliqué, vous pouvez faire d'autres choses, ** Si le nombre d'échantillons dans chaque groupe est différent ** Mettre en place.

(Je confirmerai plus tard si la formule sera la même si le nombre d'échantillons est le même pour la formule lorsque le nombre d'échantillons dans chaque groupe est différent)

fonction probmc (tour de cou des deux côtés avec différents numéros d'échantillons dans chaque groupe)

Le prob obtenu par la fonction probmc ressemble à ceci:

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

 Et la formule qui fonctionne à l'intérieur est la suivante.

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

Si vous pouvez mettre en œuvre cette formule, vous avez terminé! (Je n'avais pas l'énergie pour comprendre le contenu ...)

ici $ \ phi (x) $: fonction de densité de probabilité (pdf) $ \ Phi (x) $: Fonction de densité cumulée (cdf) Ce sera. Ceux-ci peuvent être implémentés en utilisant scipy et 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])

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

Ce sera.

De plus, $ d \ mu_v (x) $ est défini dans SAS comme suit:

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

ici $ \ Gamma (x) $: fonction Gamma Et il y a aussi une fonction dans scipy (également présente en maths).

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

** Un problème ici ... ** L'intervalle d'intégration de l'équation d'origine est $\int_{0}^{∞}d\mu_v(x)$ C'était, mais ça $\int_{0}^{∞}d\mu_v(x)=\int_{?}^{?}g_v(x)dx $ En regardant le graphique en essayant de convertir en ...

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): #changer 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…? ** Est-il inversé uniformément du côté négatif? ** C'était un problème pour remplacer l'intervalle d'intégration ... Cela ressemble à une fonction sigmoïde, donc je l'utilise probablement de cette façon. Malheureusement, je ne pouvais pas comprendre la signification de la formule, alors j'ai décidé de continuer avec la plage de ** x de 0 à ∞ ** pour le moment.

Aussi, en regardant la précision lorsque la valeur est la plus grande afin de réduire le nombre de calculs, $ \ mu_9 (3) = 0.99999999999998976 $ semble être suffisant, donc la plage de x est de 0 → 5 dans le calcul. (Au fait, $ \ mu_9 (2) = 0.99999603534120194 $)

Prêt!

Donc la formule originale,

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)

Est transformé $ prob=\int_{0}^{∞}\int_{-∞}^{∞}f_{T\lambda}(x,y)g_v(x)dydx $ est devenu. Mettre en œuvre ceci

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]

Regardons maintenant le graphique de $ f_ {T \ lambda} (x, y) g_v (x) $ pour réduire la quantité de calcul.

T, v = 10,9 #Définissez la valeur de manière appropriée
lambda_param =[0.5,0.4] #Définissez la valeur de manière appropriée

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

Cela ressemble à une fonction uniforme! En d'autres termes

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

pourtant

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]

Ce sera.

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]

Achevée! !! ...?

Je pense donc que ce sera résumé comme suit. Je voulais l'organiser proprement en utilisant des classes et ainsi de suite, mais je suis désolé que ce soit sale parce que je ne peux pas le faire bien.

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]

Je l'ai utilisé

Créez des données appropriées et essayez-les. Considérez le schéma selon lequel le médicament A et le médicament B donnent des valeurs efficaces pour le groupe témoin. D'abord de la création des données

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)

Nous irons à ces données. Tout d'abord, du moulage des données. Il est nécessaire de calculer les statistiques de test T et λ. Je voudrais prendre une méthode similaire, en pensant à la mettre sur pigouin plus tard.

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'] #Degré de liberté
  ng = aov.at[0, 'DF'] + 1 #Nombre total
  grp = data.groupby(between)[dv]
  n_sub = grp.count()
  control_index = n_sub.index.get_loc(control) #Obtenir l'index du groupe de contrôle
  n = n_sub.values #Le nombre d'échantillons
  gmeans = grp.mean().values#Chaque moyenne
  gvar = aov.at[1, 'MS'] / n #Chaque distribution

  vs_g = np.delete(np.arange(ng),control_index) #Obtenir un index autre que le groupe de contrôle

  #Combinaisons par paires (trouver la statistique de test T)
  mn = np.abs(gmeans[control_index]- gmeans[vs_g])
  se = np.sqrt(gvar[control_index] + gvar[vs_g]) #La formule est un peu différente, mais parce qu'elle suppose une distribution égale
  tval = mn / se

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

est devenu. Est-ce vraiment vrai ...?

Essayez-le avec EZR

Je suis complètement fatigué ... Enfin, j'aimerais faire la même analyse avec R en utilisant EZR. R ne peut pas être écrit, alors pardonnez-moi uniquement pour la sortie.

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

cette? ?? ?? Est-ce un peu différent? ??

Je suis fatigué, alors je ferai de mon mieux pour le revoir. Il semble que R utilise la méthode de Monte Carlo, donc je pense que les valeurs sont différentes, mais 0,005 est trop grand, n'est-ce pas?

Merci d'avoir regardé jusqu'à présent. Comme mentionné ci-dessus, les amateurs ont fait de leur mieux.

Les références

[Comparaison multiple (Partie 2-1: Méthode spécifique (i)), 3e session d'étude Armitage: Matériel 2, Masaaki Doi](http://www012.upp.so-net.ne.jp/doi/ biostat / CT39 / multiple_comparison_2_1.pdf)

Pour copie

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'] #Degré de liberté
  ng = aov.at[0, 'DF'] + 1 #Nombre total
  grp = data.groupby(between)[dv]
  n_sub = grp.count()
  control_index = n_sub.index.get_loc(control) #Obtenir l'index du groupe de contrôle
  n = n_sub.values #Le nombre d'échantillons
  gmeans = grp.mean().values#Chaque moyenne
  gvar = aov.at[1, 'MS'] / n #Chaque distribution

  vs_g = np.delete(np.arange(ng),control_index) #Obtenir un index autre que le groupe de contrôle

  #Combinaisons par paires (trouver la statistique de test T)
  mn = np.abs(gmeans[control_index]- gmeans[vs_g])
  se = np.sqrt(gvar[control_index] + gvar[vs_g]) #La formule est un peu différente, mais parce qu'elle suppose une distribution égale
  tval = mn / se

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