summarise-freundlicher wilcoxon test

[phish108]

Der wilcoxon test ist ein wichtiger Test, um nicht normalverteilte Skalenniveaus zu vergleichen.

Leider ist in infer keine entsprechende Funktion vorhanden und die Funktion wilcox.test aus stats ist mit dplyr super umständlich zu bedienen.

Hier wäre eine summarise-freundliche Funktion wünschenswert.

Die Funktion ist recht einfach und intuitiver als die standard Funktion. Sie hält sich auch an die Namenskonventionen von infer:

wilcox_test = function(response, explanatory, subset = c(), ...) {
    if (length(subset) != 2)
        subset = unique(explanatory)

    tibble(response, explanatory) %>% 
        wilcox.test(response ~ explanatory, ., explanatory %in% subset, ...) %>% 
        tidy
}

Wir können damit in summarise() mehrere Wilcoxon Tests in einem Rutsch durchführen.

daten %>% 
    summarise(
        u00_1_kurs= q00_1 %>% wilcox_test(explanatory = kurs,  alternative = "two.sided") %>% pull(p.value),
        u00_1_q00_2 =  q00_1 %>% wilcox_test(explanatory = kurs,  subset = c("M", "W"), alternative = "two.sided") %>% pull(p.value)
    )

Achtung: list() oder pull(p.value) wurden absichtlich nicht in die Funktion eingebunden, damit wir die Freiheit haben, die entsprechende Manipulation selbst durchzuführen.

Im Gegensatz zur COIN bibliothek ist dieser Test keine Neuimplementierung, sondern kapselt nur die originale Funktionalität.

[phish108]

vielleicht wäre etwas in der Art noch besser:

wilcox_test = function(data, response = NULL, explanatory, subset = c(), ...) {
    if (length(subset) != 2)
        subset = unique(explanatory)

    if (is.null(response)) 
        response = data

    tibble(response, explanatory) %>% 
        wilcox.test(response ~ explanatory, ., explanatory %in% subset, ...) %>% 
        tidy
}

Diese Funktion dürfte aber nicht funktionieren, wie wir es von dplyr funktionen erwarten würden.

Mit dieser Funktion wären dann Operationen wie mit t_test() möglich:

daten %>%
    wilcox_test(response = foo, explanatory = bar)

[phish108]

fast noch eleganter ist diese funktion:

wilcoxon_test = function(daten, order = c(), ...) {
    if (length(order) !=  2) 
        order = daten %>% pull(e) %>% unique()

    daten %>% 
        filter(e %in% order) %>% 
        mutate( e = e %>% fct_relevel(order) ) %>% 
        wilcox.test(x ~ e, ., subset = e %in% order, ...) %>%
        pluck("p.value")
} 

Allerdings macht sie Annahmen über die Variablennamen.

[phish108]

Die folgende Funktion macht jetzt alles richtig:

wilcoxon_test = function(data, response, explanatory, order = c(), ...) {
    if (!is_tibble(data))
        data = tibble(r = data, e = explanatory)
    else 
        data = data %>% select(r = {{ response }}, e = {{ explanatory }} )

    if (length(order) != 2)
        order = data %>% pull(e) %>% unique()

    data %>% 
        filter(e %in% order) %>%
        mutate(e = e %>% fct_relevel(order)) %>%
        wilcox.test(r ~ e, ., subset = e %in% order, ...) %>%
        pluck("p.value")
}

Die Parameternamen orientieren sich an den Konventionen der infer Bibliothek von `tidymodels.

Ich kann die Funktion einzeln in einer Funktionskette aufrufen oder als Teil von summarise.

Damit beide Anwendungskontexte handhabbar werden muss am Anfang der Funktion eine zusätzliche Prüfung durchgeführt werden. Ist nämlich data kein tibble, dann sind wir in einem summarise-Kontext, sonst in einem normalen Datenstrom. Im zweiten Fall müssen wir aber aufpassen, denn unser Ergebnis ist kein tibble mehr, sondern nur ein atomic Vector.

Der zentrale Punkt des order-Parameters ist übrigens, dass wir nicht nur ein subset erstellen, sondern automatisch die Reihenfolge der Level erzwingen. Das ist nämlich ein Problem, wenn wir mit read_csv() o.ä. arbeiten. In diesen Fällen faktorisiert die wilcox.test()-Funktion den erklärenden Vektor und die Reihenfolge ist nicht deterministisch. Diese Situation berücksichtigt diese Funktion, indem es das subset vorweg nimmt und die verbleibenden erklärenden Daten in die richtige Ordnung bringt.

Jedenfalls kann ich mit dieser Version die folgenden Code-Fragmente ausführen:

daten %>%
    summarise(
        u_k = q16_10 %>% wilcoxon_test(explanatory = kurs),
        u_g = q16_10 %>% wilcoxon_test(explanatory = q00_2, order = c("W", "M"))
    )

und

daten %>% wilcoxon_test(response = q16_10, explanatory = kurs)

[phish108]

Die Funktion funktioniert für unabhängige Stichproben super. Bei abhängigen Stichproben liegen die Werte aber in gepaarten Vektoren vor. In diesem Fall haben wir ein A-B Szenario und müssten entsprechend die Vektoren in die Langform bringen.

Das würde ich so ausdrücken:

daten %>%
    wicoxon_test(response = a, explanatory = b, paired = TRUE)

In diesem Fall würde order automatisch gleich c("r", "e") sein und die Daten müssten in die Langform gebracht werden.

Die folgende Funktion implementiert diese Überlegungen.

wilcoxon_test = function(data, response, explanatory, order = c(), paired = FALSE, ...) {
    if (!is_tibble(data))
        data = tibble(r = data, e = explanatory)
    else 
        data = data %>% select(r = {{ response }}, e = {{ explanatory }} )

    if (paired) {
        data = data %>% pivot_longer(c(r, e), names_to = "e", values_to = "r")
        order = c("r", "e")  # wir erzwingen die Reihenfolge
    }

    if (length(order) != 2)
        order = data %>% pull(e) %>% unique()

    data %>% 
        drop_na() %>%                # evtl. Artefakte bereinigen
        filter(e %in% order) %>%
        mutate(e = e %>% fct_relevel(order)) %>%
        wilcox.test(r ~ e, ., subset = e %in% order, paired = paired, ...) %>%
        pluck("p.value")
}

[phish108]

OK, die letzte Version hat viele klassische if Verzweigungen.

Kann man das als Datenstrom ausdrücken?

Jein. Wir brauchen eine zusätzliche Funktion, die klassische R-Akrobatik vor uns verbirgt. Diese Funktion ist eine spezielle if Funktion, die für Datenströme ausgelegt ist. Diese Funktion ist kein Filter, sonder fügt Operationen unter bestimmten Bedingungen in den Datenstrom ein. Diese Bedingungen müssen global TRUE oder FALSE ergeben und dürfen keine logischen Vektoren enthalten.

pipe_if = function(data, expression, when_true, when_false = NULL, id = FALSE) {
    fun = when_false

    if (expression) 
        fun = when_true

    # purrr lambda funktionen in Funktionen umwandeln
    if (is_formula(fun))
        fun = as_mapper(fun)

    # Wenn kein Funktionsparameter übergeben wurde, dann wollen wir einfach die Identitätsfunktion ausführen.
    # D.h. der Datenstrom bleibt unverändert. 
    # Achtung, das kann zu Irritationen führen, wenn keine Lambda-Funktion übergeben wurde. 
    if (!is_function(fun))
        fun = function(x) x   # Identitätsfunktion

    # Wenn wir eine Funktion als Parameter erhalten haben wollen wir sie in den Datenstrom einbetten.
    if (!id)
        return(data %>% fun)

    # Falls keine Veränderung erwünscht ist leiten wir die Daten nur durch. 
    # Wir müssen aber die Funktion trotzdem ausführen.
    data %>% fun
    data
}

Diese Funktion überläd die klassischen if-Bedingungen, so dass wir nun die Funktion wie folgt schreiben können.

wilcoxon_flow = function(data, response, explanatory, order = c(), paired = FALSE, ...)  
    data %>% 
        pipe_if(
            !is_tibble(.), 
            ~ tibble(r = ., e = explanatory), 
            . %>% select(r = {{ response }}, e = {{ explanatory }})
        ) %>%
        pipe_if(
            paired,
            function(d) {
                order <<- c("r", "e") # double arrow für closure zuweisung
                d %>% pivot_longer(c(r, e), names_to = "e", values_to = "r")
            }
        ) %>% 
        pipe_if(
            length(order) != 2,
            function (d) order <<- d %>% pull(e) %>% unique(),
            id = TRUE
        ) %>%
        drop_na() %>%                # evtl. Artefakte bereinigen
        filter(e %in% order) %>%
        mutate(e = e %>% fct_relevel(order)) %>%
        wilcox.test(r ~ e, ., subset = e %in% order, paired = paired, ...) %>%
        pluck("p.value")

Zugegeben. Diese Verkettung ist nicht viel einfacher als das Original. Der Punkt liegt hier beim Konditionalen Einbinden von Operationen in eine Funktionskette.

Diese Logik sollte vielleicht in einem separaten Blog beschrieben werden.