-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathscript_norm_graph.R
216 lines (185 loc) · 10.6 KB
/
script_norm_graph.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
### Script pour tester la normalité d'un jeu de données et faire l'analyse statisque
# en fonction de la distribution des données ####
#### Initialisation des variables ############################################################################################
# Définition du répertoire de travail
setwd("~/PATH/TO/DIRECTORY") #mettre le chemin complet du repertoire de travail
# Données
# pas de colonne de n° de ligne.
# les noms de colonnes doivent commencer par une lettre. Si vous avez des lignées avec des noms numériques type 1.1, nomer la L1.1
DATA <- "your_data_file.txt"
# Pour le box plot
# Chacune de ces variables peut prendre la valeur "" si on ne veut rien afficher
BOXPLOT_Xaxis <- "Titre de l'axe x"
BOXPLOT_Yaxis <- "Titre de l'axe y"
BOXPLOT_TITLE <- "Titre du graphique"
# Pour le Bar plot
# Chacune de ces variables peut prendre la valeur "" si on ne veut rien afficher
BAR_TITLE <- "Titre du graphique"
BAR_Yaxis <- "Titre de l'axe y"
###########################################################################################################################################
# Fonctions utilisées dans le script
##########################################################################################################################################
#=======================================================================================================================================
# Réalise le test de Student sur les données contenue dans le dataframe df
#=======================================================================================================================================
student.test <- function (df) {
# Crée un data frame pour stocker les pvalue du test de student
pval <- data.frame()
# Réalise le test de student sur chaque colonne de df en utilisant la première colonne comme référence
for (i in 2 : ncol(df)) {
pval[i,1] <- as.data.frame(t.test(df[,1], df[,i], alternative = "two.sided")$p.value)
}
# Nomme colonnes et lignes du dataframe pval
colnames(pval) <- c("pvalue")
rownames(pval) <- colnames(df)
pval$padjust <- p.adjust(pval$pvalue, method = "BH", n = length(pval$pvalue))
# Renvoie le resultat du test de Student sous la forme d'un dataframe
return(pval)
}
#=======================================================================================================================================
#=======================================================================================================================================
# Identifie la plus grande valeur du jeu de donnée pour établir la hauteur de l'axe y
#=======================================================================================================================================
maximum_y <- function () {
ymax <- max (as.numeric(resume[4,]) + sd_df$SD)
}
#=======================================================================================================================================
#=======================================================================================================================================
# Crée un Barplot représentant les moyennes des données contenues dans df
#=======================================================================================================================================
# La fonction prend en paramètre le titre principal du graphique et le titre de l'axe y
graph.norm <- function(titre_graph, titre_axe_y, couleur = rainbow(ncol(df))) {
ymax <- max (as.numeric(resume[4,]) + sd_df$SD)
# Dessine le barplot
plot <- barplot(c(as.numeric(resume[4, ])),
col = couleur,
ylim = c(0, ymax + ymax*0.2),
main = titre_graph,
ylab = titre_axe_y,
names.arg = colnames(df)
)
# Ajout des barres d'erreurs : en vérité des fléches à double pointe (code=3) et pointes horizontales (angle = 90)
arrows(plot[,1], as.numeric(resume[4,]) - sd_df$SD, plot[,1], as.numeric(resume[4,]) + sd_df$SD, angle = 90, code = 3)
return(plot)
}
#=======================================================================================================================================
#=======================================================================================================================================
# Positonne les étoiles de significativité dans le cas des données normales
#=======================================================================================================================================
star_norm <- function() {
for (i in 1 : ncol(df)) {
x.star <- mean(plot[i, 1]) #permet de sélectionner la colonne voulu dans le plot
y.star <- ymax*0.05 + (as.numeric(resume[4, i]) + sd_df[i, 1]) #position en y des étoiles 5% du ymax au dessus de la barre d'erreur de chaque colonne
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
star = ""
if (is.na(pval[i,2])) { # test le cas ou la pvalue n'est pas calculée (Na)
star = " "
} else if (pval[i, 2] < 0.001) {
star = "***"
} else if (pval[i, 2] < 0.01) {
star = "**"
} else if (pval[i, 2] < 0.05) {
star = "*"
} else if (pval[i, 2] < 0.1) {
star = "."
} else {
star = " "
}
text(star, x = x.star, y = y.star)
}
}
#=======================================================================================================================================
#=======================================================================================================================================
# Fait le test de Kruskal-Wallis
#=======================================================================================================================================
# Retourne la pvalue du test
test_kruskal <- function() {
pval_test <- kruskal.test((as.list(df)))$p.value
if (pval_test < 0.05) {
print ("Le test de Kruskall Wallis compare les médianes, la pvalue est < 0.05 ce qui indique qu’au moins 1 des médianes est différentes des autres, on réalise un test post hoc de Dunn")
test_dunn()
} else {
print ("Le test de Kruskall Wallis compare les médianes, la pvalue est > 0.05 ce qui indique qu’il n'y a pas de différence entre les médianes")
}
return(pval_test)
}
#=======================================================================================================================================
#=======================================================================================================================================
# Fait le test de Dunn
#=======================================================================================================================================
# retourne les pvalue dans un dataframe
test_dunn <- function() {
# installer le package DescTools si nécessaire
if (!requireNamespace("DescTools", quietly = TRUE))
install.packages("DescTools")
library(DescTools)
pval <- as.data.frame(DunnTest(as.list(df), method = "BH", out.list = FALSE)[1])
print(DunnTest(as.list(df), method = "BH"))
return(pval)
}
#=======================================================================================================================================
#=======================================================================================================================================
# Positonne les étoiles de significativité dans le cas des données non paramétriques
#=======================================================================================================================================
star_non_param <- function() {
for (i in 1 : (ncol(df)-1)) {
x.star <- mean(plot[i+1, 1]) #permet de sélectionner la colonne voulu dans le plot. Le +1 permet de commencer à ajouter les étoiles sur la deuxieme colonne
y.star <- ymax*0.05 + (as.numeric(resume[4, i+1]) + sd_df[i+1, 1]) #position en y des étoiles 5% du ymax au dessus de la barre d'erreur de chaque colonne
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
star = ""
if (is.na(pval[i,1])) {
star = " "
} else if (pval[i, 1] < 0.001) {
star = "***"
} else if (pval[i, 1] < 0.01) {
star = "**"
} else if (pval[i, 1] < 0.05) {
star = "*"
} else if (pval[i, 1] < 0.1) {
star = "."
} else {
star = " "
}
text(star, x = x.star, y = y.star)
}
}
######################################################################################################################################
######################################################################################################################################
# Partie principale du script
# Chargement des données
# pas de colonne de n° de ligne.
# les noms de colonnes doivent commencer par une lettre. Si vous avez des lignées avec des noms numériques type 1.1, nomer la L1.1
df <- read.table(DATA, header = TRUE, sep = "\t")
####################################################################################################################################
# I. Statistiques descriptives
# On place dans un tableau les quantiles, la moyenne et la médiane pour toutes les colonnes du fichier de départ
resume <- as.data.frame(do.call(cbind, lapply(df, summary)))
# II. Boxplot des données
boxplot(df,
xlab = BOXPLOT_Xaxis, ylab = BOXPLOT_Yaxis,
main = BOXPLOT_TITLE,
col = rainbow(ncol(df))) # rainbow(ncol(df)) peut être remplacer par un vecteur de couleur (autant que de colonne) du type c("red", "green").
# III. Calcul de la Déviation standard écart type/racine(nb de mesure)
sd_df <- as.data.frame(apply(df, 2, sd, na.rm = TRUE)/sqrt(nrow(df)))
colnames(sd_df) <- c("SD")
# IV. test de normalité de Shapiro sur chaque colonne
shapiro_df <- apply(df, 2, shapiro.test)
if (sum(do.call(rbind,lapply(shapiro_df,function(v){v$p.value})) > 0.05) == ncol(df)) {
print ("Les données de df suivent une loi normale, on réalise un test de Student")
pval <- student.test(df)
ymax <- maximum_y()
plot <- graph.norm("Titre du graphique", "Titre de l'axe y")
star_norm()
} else {
print ("Les données de df ne suivent pas une loi normale, on réalise un test de Kruskall Wallis et un test de Dunn")
pavl_test <- test_kruskal()
pval <- test_dunn()
ymax <- maximum_y()
plot <- graph.norm("Titre du graphique", "Titre de l'axe y")
points(x = plot, y = resume[3, ])
star_non_param()
}
######################################################################################################################################
# R session
if (!require(devtools)) { install.packages("devtools") }
devtools::session_info()