Discussion :

[Julien H]

Bonjour,

Je travaille sur les poissons, je cherche à savoir si il existe des différences de biomasse entre mes différentes stations de prélèvements.
J'ai donc plusieurs stations de prélèvements et un certain nombre de valeurs de masse pour chacune de ces stations (nombre variable en fonction des stations).
Cf PJ

J'ai utilisé un test de KW pour savoir si les différences observées étaient significatives (distribution non normale), elle le sont.
Je voudrais maintenant savoir quelles sont ces différences (entre quelles stations).

Quel test post-hoc dois je utiliser ?
J'ai essayé le test de Mann Whitney mais les résultats sont à mon sens assez étranges (la différence la plus significative est observée entre deux stations pour lesquelles la différence graphique est très faible ...)


Merci d'avance

Julien

[Julien H]

Personne n'a de réponse ?

[A. D.]

Bonjour,

Votre question me semble plus de l'ordre de la statistique et non de la programmation avec le logiciel R, qui est le thème de ce forum.
De plus, difficile de vous répondre sans savoir à quoi ressemblent vos données, quelles tests vous avez essayés et quels résultats vous avez obtenus.
Vous pourriez par exemple poster un cours extrait de vos données et le script utilisé ainsi que les résultats obtenus, peut-être aurez-vous plus de chance d'obtenir une réponse. De plus, les membres du forum sont bénévoles donc peut-être que personne n'a encore eu le temps de se pencher sur votre problème.

Bonne continuation


Cordialement,


A.D.

[Julien H]

Quelques précisions donc.

Voilà le script :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
####Loading database####
data2 <- read.table("Biomasse2.csv", header=T, sep=";", dec = ",", strip.white=T)
data2$Marae <- as.numeric(levels(data2$Marae))[data2$Marae]
data2$RDA2 <- as.numeric(levels(data2$RDA2))[data2$RDA2]
data2$RDA1 <- as.numeric(levels(data2$RDA1))[data2$RDA1]
data2$Pont2 <- as.numeric(levels(data2$Pont2))[data2$Pont2]
data2$Pont1 <- as.numeric(levels(data2$Pont1))[data2$Pont1]
 
####Test KW####
a=data2$Marae
b=data2$RDA2
c=data2$RDA1
d=data2$Pont2
e=data2$Pont1
 
kruskal.test(list(a,b,c,d,e))

Les données sont en PJ de mon premier post. Je les remets ici.

Les résultats obtenus sont :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
data:  list(a, b, c, d, e)
Kruskal-Wallis chi-squared = 123.8981, df = 4, p-value <2.2e-16

Et ma question est : quel test post-hoc utiliser pour savoir où se situent mes différences ? Sachant que les résultats donnés par Man-Whitney me semblent surprenants (pas de différences entre des données graphiquement très différentes, cf. graphique premier post)

Merci d'avance !

[bordi]

Mon niveau n'est pas encore très élevé en statistique, mais essayez de regarder ca si cela peut répondre.

Multiple comparison test after Kruskal-Wallis

http://www.inside-r.org/packages/cra...docs/kruskalmc

[juliatheric]

La réponse courte est celle de bordi. Une version longue ici-bas.

En plus de la méthode stats::kruskal.test, il y a la méthode coin::kruskal_test qui offre plus de flexibilité pour accomplir le test. Par défaut, la dernière va donner la même statistique que la première. La flexibilité réside dans l'utilisation de coin::independence_test où l'on peut personnaliser comment calculer la statistique. Les sorties ajoutent également plus d'informations.

Code R :

Sélectionner tout - Visualiser dans une fenêtre à part

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
 
invisible(sapply(c(
  "coin", "ggplot2", "ggthemes", "latex2exp", "magrittr", "SuppDists", "pgirmess"),
  function(z) library(z, character.only = TRUE)
))
 
biomasse <- read.delim("biomasse.tsv")
 
biomaSansValManq <- sapply(biomasse, function(z) {
  tmp <- na.exclude(z)
  attr(tmp, "na.action") <- NULL
  return(tmp)
})
 
biomaSansValManq$relevé <- unlist(biomaSansValManq)
biomaSansValManq$effectif <- sapply(biomasse, function(z) na.exclude(z) %>% length)
biomaSansValManq$station <- rep(colnames(biomasse), biomaSansValManq$effectif) %>% factor
 
###stats
sorties$stats <- kruskal.test(x = biomasse, na.action = "na.exclude")
print(sorties$stats)
print(sorties$stats$statistic)
print(sorties$stats$p.value)
print(sorties$stats$parameter)
 
sorties$alpha <- 1e-3
sorties$décision <- ifelse(
  1L - pchisq(sorties$stats$statistic, sorties$stats$parameter) < sorties$alpha,
  paste0("Rejet de H0 avec alpha = ", sorties$alpha, ".\n"),
  "Non-rejet de H0.\n"
)
cat(sorties$décision)
 
###coin
sorties$coin <- kruskal_test(
  formula = biomaSansValManq$relevé ~ biomaSansValManq$station,
  distribution = approximate(B = 1000L)
)
show(sorties$coin)
statistic(sorties$coin, type = "linear")
statistic(sorties$coin, type = "standardized")
statistic(sorties$coin, type = "test")
statistic(sorties$coin)
pvalue(sorties$coin)
expectation(sorties$coin)

On voit déjà qu'avec coin, on peut notamment obtenir les rangs moyens des échantillons dans le test K.-W. et connaître l'intervalle de confiance associé à la p-valeur. Pour se faire une idée sur les relations linéaires qui existeraient entre les échantillons, on peut faire ainsi :

Code R :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
 
stationCov <- covariance(sorties$coin)
stationEcarts <- variance(sorties$coin) %>% sqrt
 
k <- 0L
stationCor <- matrix(ncol = length(stationEcarts), nrow = length(stationEcarts))
for (i in nrow(stationCov) %>% seq_len) {
  for (j in seq_len(ncol(stationCov) - k) ) {
    stationCor[i, j + k] <- stationCor[j + k, i] <- 
      stationCov[i, j + k] / (stationEcarts[i] * stationEcarts[j + k])
  }
  k <- k + 1L
}
dimnames(stationCor) <- dimnames(stationCov)

La variable stationCor contient les coefficients de corrélation entre les 5 stations. Comme on a rejeté l'hypothèse nulle, on peut passer à pgirmess::kruskalmc.

Code R :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
 
#À combien de comparaisons 2-à-2 s'attend-on ?
choose(ncol(biomasse), 2L)
 
sorties$multicomp <- kruskalmc(
  resp = biomaSansValManq$relevé,
  categ = biomaSansValManq$station,
  probs = sorties$alpha
)

Les résultats.
* Coefficients de corrélation

           Marae      Pont1      Pont2       RDA1       RDA2
Marae  1.0000000 -0.2928710 -0.1382039 -0.1711245 -0.1279573
Pont1 -0.2928710  1.0000000 -0.3517217 -0.4355030 -0.3256448
Pont2 -0.1382039 -0.3517217  1.0000000 -0.2055109 -0.1536696
RDA1  -0.1711245 -0.4355030 -0.2055109  1.0000000 -0.1902742
RDA2  -0.1279573 -0.3256448 -0.1536696 -0.1902742  1.0000000

* Comparaisons 2-à-2 :

              obs.dif critical.dif difference
Marae-Pont1  96.42328     65.41936       TRUE
Marae-Pont2  51.12328     77.10780      FALSE
Marae-RDA1   20.37629     72.11333      FALSE
Marae-RDA2   31.09458     79.38889      FALSE
Pont1-Pont2  45.30000     57.72210      FALSE
Pont1-RDA1  116.79956     50.85824       TRUE
Pont1-RDA2  127.51786     60.73569       TRUE
Pont2-RDA1   71.49956     65.21105       TRUE
Pont2-RDA2   82.21786     73.17612       TRUE
RDA1-RDA2    10.71830     67.89304      FALSE

La méthode de comparaison multiple est expliquée ici : http://giraudoux.pagesperso-orange.f...tellan1988.pdf
En bonus, l'exploration graphique des données initiales et de la distribution K.-W.

Code R :

Sélectionner tout - Visualiser dans une fenêtre à part

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
 
dessBiomasse <- data.frame(
  Biomasse = biomaSansValManq$relevé, 
  Station = biomaSansValManq$station
)
 
biomBatonsIndiv <- ggplot(data = dessBiomasse, aes(x = Biomasse)) + 
  scale_x_log10() + geom_bar(binwidth = 0.5) + facet_grid(. ~ Station) +
  ylab("Effectif") + theme_light()
 
biomBatons <- ggplot(data = dessBiomasse, aes(x = Biomasse, fill = Station)) +
  scale_x_log10() + geom_bar(binwidth = 0.5) + ylab("Effectif") +
  theme_light() + theme(legend.position = "bottom")
 
biomDens <- ggplot(data = dessBiomasse, aes(x = Biomasse, color = Station)) +
  scale_x_log10() + geom_freqpoly(binwidth = 1) + ylab("Effectif") +
  scale_color_brewer(palette = "Set1") + theme_light() +
  theme(legend.position = "bottom")
 
biomBoiteVert <- ggplot(
  data = dessBiomasse,
  aes(x = factor(0L), y = Biomasse, fill = Station, color = Station)) +
  scale_y_log10() + geom_boxplot() + xlab("") + 
  scale_color_brewer(palette = "Dark2") + scale_fill_brewer(palette = "PuOr") +
  theme_light() + theme(legend.position = "bottom")
 
biomBoiteHoriz <- ggplot(
  data = dessBiomasse,
  aes(x = factor(0L), y = Biomasse, fill = Station, color = Station)) +
  scale_y_log10() + geom_boxplot() + xlab("") + coord_flip() +
  scale_fill_hc() + theme_hc()
 
print(biomBatonsIndiv)
print(biomBatons)
print(biomDens)
print(biomBoiteVert)
print(biomBoiteHoriz)
 
sorties$sommeInv <- sum(1L / biomaSansValManq$effectif)
sorties$taille <- sum(biomaSansValManq$effectif)
sorties$nEch <- ncol(biomasse)
 
biomFunc <- function(x = biomaSansValManq$relevé) {
  dKruskalWallis(x, sorties$nEch, sorties$taille, sorties$sommeInv)
}
biomTitreFunc <- latex2exp(paste0(
  "$\\frac{12}{n * (n + 1)} * ", 
  "\\sum_{i = 1}{k}\\frac{R_i ^ 2}{n_i} - 3 * ",
  "(n + 1)$"
))
 
biomKWDistr <- data.frame(x = biomaSansValManq$relevé, y = biomFunc()) %>% 
  ggplot(aes(x, y)) + geom_line() + ylab("") + ggtitle(biomTitreFunc)
 
print(biomKWDistr)

P.S. L'opérateur %>% de magrittr est là seulement pour éviter l'empilement des parenthèses, il n'est pas indispensable.