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 ...
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 **.
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#)**.
C'est merveilleux! !! Merci au créateur. Avec cela, il semblait que je pourrais vivre de python seul.
Pourtant…
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 ~~)
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
En ce moment, ** la liberté v ** est
Par rapport à cette statistique de test T ...
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 ...
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
(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)
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
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
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 $)
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é
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]
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]
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 ...?
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
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.
[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)
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