Paarweise Chiquadrat-Test in R

Ich erstelle eine naive Funktion, welche paarweise Chiquadrat-Tests in R durchführt. Die Funktion nimmt zwei Faktoren entgegen, einer für die Werte, und der andere für die Gruppierung.
R
Autor:in

Joe Slam

Veröffentlichungsdatum

24 Mai 2024 - 17:48

Geändert

02 Juni 2024 - 14:16

In R ist keine Funktion implementiert, die einen paarweisen Chiquadrat-Test durchführt. Das mag daran liegen, dass der Chiquadrat-Test als (weniger genaue) Alternative zu dem rechenaufwändigen exakten Fisher-Test entwickelt wurde. Mit der heutigen Computertechnologie ist es jedoch in den meisten Fällen kein Problem, den exakten Test nach Fisher auch auf älterer Hardware anzuwenden, siehe hierzu https://www.produnis.de/R/schliessende-statistik.html#sec-posthocsig.

Paarweise Chiquadrat-Tests

Ich habe aus edukativen Gründen eine Funktion programmiert, die dennoch paarweise Chiquadrat-Tests durchführt.

Die Funktion ist auch in meinem Zusatzpaket jgsbook enthalten und kann mittels jgsbook::pairwise.chisq.test(x, y) aufgerufen werden.

Durch die Mehrfachvergleiche (multiples Testen) entsteht das erhöhte Risiko von falsch-positiven Ergebnissen (Typ-I-Fehler). Bei unserer Funktion müssen wir das Problem der Mehrfachtestung berücksichtigen. In Zeile 50 übergeben wir die berechneten p-Werte der p.adjust()-Funktion, und korrigieren die Werte mit dem gewünschten Verfahren (standardmäßig bonferroni).

pairwise.chisq.test <- function(A, B, p.adjust.method="bonferroni"){
  
  # wir überprüfen, ob die Objekte vom Typ factor sind.
  if (!is.factor(A) | !is.factor(B)) {
    stop("Both objects need to be factors.")
  }

  # erstelle ein Hilfsdatenframe
  df <- data.frame(A=A, B=B)

  # Erzeuge alle möglichen 2er-Kombinationen von B
  kombi <- combn(levels(B), 2)

  # Erzeuge ein Datenframe für die Werterückgabe
  # die erste Zeile wird später wieder gelöscht.
  result <- data.frame(group1 = "A",
                       group2 = "B",
                       xsquare = "9",
                       df = "9999",
                       pvalue = "9",
                       padjust = "9",
                       adjust = p.adjust.method,
                       sig = "")

  # Berechne für jedes Paar
  for (i in 1:(length(kombi)/2)) {
    # füge die aktuelle Paarung in einem Datenframe zusammen
    help <- rbind(df[df$B==kombi[,i][1],],
                  df[df$B==kombi[,i][2],])
    # führe den Chi-Quadrat-Test durch
    pups <- chisq.test(help$A, help$B)
    
    # setze das Signifikanz-Signal zurück
    sig = ""
    
    # setze einen Punkt wenn p < 0.055
    if(p.adjust(as.numeric(pups$p.value), method = p.adjust.method, n=(length(kombi)/2)) < 0.055) {sig="."}
    # setze einen Stern wenn p < 0.05
    if(p.adjust(as.numeric(pups$p.value), method = p.adjust.method, n=(length(kombi)/2)) < 0.05)  {sig="*"}
    # setze zwei Sterne wenn p < 0.01
    if(p.adjust(as.numeric(pups$p.value), method = p.adjust.method, n=(length(kombi)/2)) < 0.01)  {sig="**"}
    # setze drei Sterne wenn p < 0.001
    if(p.adjust(as.numeric(pups$p.value), method = p.adjust.method, n=(length(kombi)/2)) < 0.001) {sig="***"}
    # speichere die Ergebnisse in ein Datenframe
    loop <- data.frame(group1 = kombi[,i][1],
                       group2 = kombi[,i][2],
                       xsquare = as.numeric(pups$statistic),
                       df = as.numeric(pups$parameter),
                       pvalue = as.numeric(pups$p.value),
                       padjust = p.adjust(as.numeric(pups$p.value), method = p.adjust.method, n=(length(kombi)/2)),
                       adjust = p.adjust.method,
                       sig = sig
    )
    # füge die aktuellen Ergebnisse dem "result"-Objekt hinzu
    result <- rbind(result, loop)
  }
  # Fertig mit allen Paarungen
  
  # Entferne den Dummy-Eintrag 
  result <- result[-1,]
  # Korrigiere die Zeilennummerierung
  row.names(result) <- NULL
  # Gebe die Ergebnisse zurück
  return(result)
}
#----------------------------------------------------------------------#

Die Funktion nimmt folgende Parameter entgegen:

  • A: die erste Variable, ein Faktor, welche die Zieldaten enthält (z.B. Geschlecht)
  • B: die zweite Variable, ein Faktor, welcher die Gruppeneinteilung enthält (z.B. Ort)
  • p.adjust.method: die Methode zur Alphafehler-Korrekur, standardmäßig “bonferroni”, es gehen aber auch “alle anderen”, wie zB holm oder hochberg.

Als Ergebnis liefert uns die Funktion ein Datenframe zurück, mit den Einträgen:

  • Name von Gruppe 1
  • Name von Gruppe 2
  • Teststatistik (xsquare)
  • Freiheitsgrade (df)
  • p-Wert des Chi-Quadrat-Test
  • adjustierter p-Wert (nach gewählter Methode)
  • die gewählte Methode zur p-Wert-Adjustierung
  • Signifikanz-Signal für die adjustierten p-Werte (., *, **, ***)

Beispiel

Als Beispiel untersuchen wir im Datensatz pf8, ob sich die Geschlechter je nach Standort unterscheiden.

# verwende pf8 Datensatz
load(url("https://www.produnis.de/R/data/pf8.RData"))

# führe exakten Test durch
pairwise.chisq.test(pf8$Geschlecht, pf8$Standort)
      group1    group2           xsquare df               pvalue              padjust     adjust sig
1     Rheine   Münster  1.10471703224467  2     0.57559066836014                    1 bonferroni    
2     Rheine      Bahn 0.484729064039409  2    0.784770051592051                    1 bonferroni    
3     Rheine Ladbergen  3.82435691920203  2    0.147758150572077                    1 bonferroni    
4     Rheine  Internet  14.0768553245862  2 0.000877505212656076  0.00877505212656076 bonferroni  **
5    Münster      Bahn                 0  1                    1                    1 bonferroni    
6    Münster Ladbergen  3.97980963867552  1   0.0460487671171598    0.460487671171598 bonferroni    
7    Münster  Internet   8.4878079618588  1  0.00357534323663839   0.0357534323663839 bonferroni   *
8       Bahn Ladbergen  2.85957778936492  1   0.0908313368796546    0.908313368796546 bonferroni    
9       Bahn  Internet  5.76636469124268  1   0.0163357933303218    0.163357933303218 bonferroni    
10 Ladbergen  Internet  16.0075693112762  1 6.30897415058098e-05 0.000630897415058098 bonferroni ***

Es ist zu erkennen, dass sich die Geschlechterverteilungen zwischen den Standortvergleichen Internet-Rheine, Internet-Münster und Internet-Landbergen unterscheiden.

Vergleich mit Fisher Test

Vergleichen wir die Ergebnisse mit den exakten Fisher-Tests:

reporttools::pairwise.fisher.test(pf8$Geschlecht, 
                                  pf8$Standort, 
                                  p.adjust.method = "bonferroni")

    Pairwise comparisons using 

data:  pf8$Geschlecht and pf8$Standort 

          Rheine  Münster Bahn    Ladbergen
Münster   1.00000 -       -       -        
Bahn      1.00000 1.00000 -       -        
Ladbergen 1.00000 0.40161 0.64522 -        
Internet  0.00292 0.02852 0.13207 0.00065  

P value adjustment method: bonferroni

Es ist erkennbar, dass beide Tests zwar in die selbe Richtung tendieren, und auch die gleichen Signifikanzbewertungen ausgeben. Jedoch sind die p-Werte mit unter deutlich abweichend.

Unsere Funktion ist eine nette Übung. Am Besten ist es aber, immer den exakten Fisher-Test zu rechnen!