IdentifiantMot de passe
Loading...
Mot de passe oublié ?Je m'inscris ! (gratuit)
Navigation

Inscrivez-vous gratuitement
pour pouvoir participer, suivre les réponses en temps réel, voter pour les messages, poser vos propres questions et recevoir la newsletter

Tests d'hypothèse Discussion :

Test post-hoc après Kruskal Wallis


Sujet :

Tests d'hypothèse

  1. #1
    Candidat au Club
    Homme Profil pro
    Étudiant
    Inscrit en
    mai 2015
    Messages
    4
    Détails du profil
    Informations personnelles :
    Sexe : Homme
    Localisation : France, Gironde (Aquitaine)

    Informations professionnelles :
    Activité : Étudiant

    Informations forums :
    Inscription : mai 2015
    Messages : 4
    Points : 2
    Points
    2
    Par défaut Test post-hoc après Kruskal Wallis
    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 ...)
    Nom : Sans titre.png
Affichages : 1081
Taille : 9,3 Ko

    Merci d'avance

    Julien
    Fichiers attachés Fichiers attachés

  2. #2
    Candidat au Club
    Homme Profil pro
    Étudiant
    Inscrit en
    mai 2015
    Messages
    4
    Détails du profil
    Informations personnelles :
    Sexe : Homme
    Localisation : France, Gironde (Aquitaine)

    Informations professionnelles :
    Activité : Étudiant

    Informations forums :
    Inscription : mai 2015
    Messages : 4
    Points : 2
    Points
    2
    Par défaut
    Personne n'a de réponse ?

  3. #3
    Modératrice

    Femme Profil pro
    Statisticienne, Fondatrice de la société DACTA
    Inscrit en
    juin 2010
    Messages
    893
    Détails du profil
    Informations personnelles :
    Sexe : Femme
    Âge : 34
    Localisation : France, Loire Atlantique (Pays de la Loire)

    Informations professionnelles :
    Activité : Statisticienne, Fondatrice de la société DACTA

    Informations forums :
    Inscription : juin 2010
    Messages : 893
    Points : 2 673
    Points
    2 673
    Par défaut
    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.

    Forum R
    Fournir le code utilisé (pensez aux balises code !), les packages nécessaires, ainsi qu'un court mais représentatif extrait du jeu de données et les éventuels messages d'erreur.
    Recherche d'informations concernant R : RSiteSearch / tutoriels : http://r.developpez.com/cours/ .

    Pensez également au bouton "Résolu" et à voter (en bas à droite des messages) lorsque vous avez obtenu une réponse satisfaisante.

  4. #4
    Candidat au Club
    Homme Profil pro
    Étudiant
    Inscrit en
    mai 2015
    Messages
    4
    Détails du profil
    Informations personnelles :
    Sexe : Homme
    Localisation : France, Gironde (Aquitaine)

    Informations professionnelles :
    Activité : Étudiant

    Informations forums :
    Inscription : mai 2015
    Messages : 4
    Points : 2
    Points
    2
    Par défaut
    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 !
    Fichiers attachés Fichiers attachés

  5. #5
    Membre habitué
    Homme Profil pro
    Inscrit en
    octobre 2007
    Messages
    190
    Détails du profil
    Informations personnelles :
    Sexe : Homme
    Localisation : France

    Informations forums :
    Inscription : octobre 2007
    Messages : 190
    Points : 182
    Points
    182
    Par défaut
    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

  6. #6
    Membre actif
    Homme Profil pro
    Bioinformaticien
    Inscrit en
    octobre 2008
    Messages
    126
    Détails du profil
    Informations personnelles :
    Sexe : Homme
    Localisation : Autre

    Informations professionnelles :
    Activité : Bioinformaticien
    Secteur : Enseignement

    Informations forums :
    Inscription : octobre 2008
    Messages : 126
    Points : 296
    Points
    296
    Par défaut
    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)

    Nom : biomasse.png
Affichages : 1098
Taille : 30,1 Ko

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

Discussions similaires

  1. Réponses: 4
    Dernier message: 10/06/2015, 10h12
  2. Réponses: 3
    Dernier message: 24/05/2012, 15h58
  3. Friedman test post-hoc
    Par julieh47 dans le forum SAS STAT
    Réponses: 1
    Dernier message: 08/04/2010, 16h53
  4. Réponses: 0
    Dernier message: 11/12/2009, 20h58
  5. test post-hoc non-parametriques?
    Par ganod dans le forum SAS STAT
    Réponses: 0
    Dernier message: 31/07/2009, 12h07

Partager

Partager
  • Envoyer la discussion sur Viadeo
  • Envoyer la discussion sur Twitter
  • Envoyer la discussion sur Google
  • Envoyer la discussion sur Facebook
  • Envoyer la discussion sur Digg
  • Envoyer la discussion sur Delicious
  • Envoyer la discussion sur MySpace
  • Envoyer la discussion sur Yahoo