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.

T{1i}=\frac{|μ_1-μ_i|}{\sqrt{(\frac{1}{n_1}+\frac{1}{n_i})\sigma^2}}

pourtant

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

argument Contenu
distribution Choisissez si Dunnett est un côté ou les deux côtés. Il semble qu'il puisse être utilisé pour des distributions autres que Dunnett.
q, prob q est le point de division (ici, la statistique de test T est entrée), et prob est la probabilité. prob est la probabilité que la variable de probabilité soit inférieure à q. Par conséquent, la valeur p peut être calculée comme 1 - prob. Par exemple, niveau de signification 5%Pour calculer la valeur critique de prob= 0.Réglez sur 95. L'un ou l'autre est manquant et utilisé comme argumentL'argument manquant est la valeur de retourCe sera. Cette fois, je veux connaître la valeur p, donc je vais supprimer le prob.
df Degré de liberté v
nparms Dans le cas de Danette, le nombre total de groupes k-1
parameters λ12,λ13,...λ1i

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

ダウンロード (1).png

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

ダウンロード (2).png

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]
g_v(x)=\frac{ν^{\frac{v}{2}}}{\Gamma(\frac{v}{2})2^{\frac{v}{2}-1}}x^{v-1}exp(-\frac{vx^2}{2})

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

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

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

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

Recommended Posts

Je veux faire le test de Dunnett en Python
Je veux écrire en Python! (2) Écrivons un test
Je veux faire quelque chose avec Python à la fin
Je veux faire quelque chose comme sort uniq en Python
[Python] Ce que j'ai fait pour faire un test unitaire
Je veux créer une fenêtre avec Python
Je veux fusionner des dictionnaires imbriqués en Python
Je veux afficher la progression en Python!
Je veux faire un patch monkey seulement en partie en toute sécurité avec Python
Je veux écrire en Python! (1) Vérification du format de code
Je souhaite intégrer une variable dans une chaîne Python
Je veux facilement implémenter le délai d'expiration en python
Même avec JavaScript, je veux voir Python `range ()`!
Je veux échantillonner au hasard un fichier avec Python
Je veux travailler avec un robot en python.
Je veux écrire en Python! (3) Utiliser des simulacres
Je veux utiliser le jeu de données R avec python
Je veux manipuler des chaînes dans Kotlin comme Python!
[Python] Comment faire PCA avec Python
Je veux faire ○○ avec les Pandas
Je veux déboguer avec Python
Je veux pouvoir exécuter Python avec VS Code
Je veux ajouter un joli complément à input () en python
J'ai essayé d'implémenter la permutation en Python
Je veux imprimer dans la notation d'inclusion
Comment faire R chartr () en Python
J'ai essayé d'implémenter PLSA dans Python 2
Je veux utiliser jar de python
Je veux créer un environnement Python
Je veux analyser les journaux avec Python
Comment faire un test de sac avec python
Je veux jouer avec aws avec python
J'ai essayé d'implémenter ADALINE en Python
Je voulais résoudre ABC159 avec Python
J'ai essayé d'implémenter PPO en Python
Je veux intégrer Matplotlib dans PySimpleGUI
Je veux résoudre APG4b avec Python (seulement 4.01 et 4.04 au chapitre 4)
Je souhaite rechercher le texte intégral avec elasticsearch + python
[Couches Python / AWS Lambda] Je souhaite réutiliser uniquement le module dans AWS Lambda Layers
Je voulais faire quelque chose comme la pipe d'Elixir en Python
Je veux réussir le test G dans un mois Jour 1
Je veux le faire avec Python lambda Django, mais je vais m'arrêter
Je veux utiliser MATLAB feval avec python
Je veux corriger Datetime.now dans le test de Django
Je veux mémoriser, y compris les arguments de mots clés de Python
Python: j'ai pu récurer en lambda
Je souhaite envoyer un e-mail depuis Gmail en utilisant Python.
[Python] Je veux gérer 7DaysToDie depuis Discord! 1/3
Implémentation minimale d'Union Find en Python
Je veux faire un jeu avec Python
J'ai écrit "Introduction à la vérification des effets" en Python
Je ne veux pas passer un test de codage
Je souhaite stocker les informations de la base de données dans la liste
Je souhaite convertir une table convertie en PDF en Python en CSV
Je veux convertir par lots le résultat de "chaîne de caractères" .split () en Python
Je veux expliquer en détail la classe abstraite (ABCmeta) de Python
J'ai essayé d'implémenter TOPIC MODEL en Python
Je souhaite utiliser le répertoire temporaire avec Python2
Maintenance de l'environnement réalisée avec Docker (je souhaite post-traiter GrADS en Python
Je veux colorer une partie de la chaîne Excel avec Python