library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.2     ✔ tibble    3.3.0
## ✔ lubridate 1.9.4     ✔ tidyr     1.3.1
## ✔ purrr     1.1.0     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(dplyr)
library(readxl)
library(knitr)

my_data <- read_excel("Project wormen.xlsx", sheet = "R")

kable(my_data)
Mutatie kwijt Geen paralyse Gedeeltelijke paralyse volledige paralyse Conditie herhaling
2 8 0 0 WT16 1
2 8 0 0 WT16 2
0 7 2 0 WT16 3
0 7 3 0 WT25 1
4 4 2 0 WT25 2
4 1 5 0 WT25 3
2 9 0 0 SNA16 1
0 7 3 0 SNA16 2
1 5 3 0 SNA16 3
2 6 1 1 SNA25 1
0 3 5 0 SNA25 2
4 3 3 0 SNA25 3

Maak de data tidy op een manier die gebruikt kan worden voor de kruskal wallis test

long_data <- my_data %>%
  pivot_longer(
    cols = c(`Mutatie kwijt`,
             `Geen paralyse`,
             `Gedeeltelijke paralyse`,
             `volledige paralyse`),
    names_to = "Gedrag",
    values_to = "Observaties"
  )

kable(long_data)
Conditie herhaling Gedrag Observaties
WT16 1 Mutatie kwijt 2
WT16 1 Geen paralyse 8
WT16 1 Gedeeltelijke paralyse 0
WT16 1 volledige paralyse 0
WT16 2 Mutatie kwijt 2
WT16 2 Geen paralyse 8
WT16 2 Gedeeltelijke paralyse 0
WT16 2 volledige paralyse 0
WT16 3 Mutatie kwijt 0
WT16 3 Geen paralyse 7
WT16 3 Gedeeltelijke paralyse 2
WT16 3 volledige paralyse 0
WT25 1 Mutatie kwijt 0
WT25 1 Geen paralyse 7
WT25 1 Gedeeltelijke paralyse 3
WT25 1 volledige paralyse 0
WT25 2 Mutatie kwijt 4
WT25 2 Geen paralyse 4
WT25 2 Gedeeltelijke paralyse 2
WT25 2 volledige paralyse 0
WT25 3 Mutatie kwijt 4
WT25 3 Geen paralyse 1
WT25 3 Gedeeltelijke paralyse 5
WT25 3 volledige paralyse 0
SNA16 1 Mutatie kwijt 2
SNA16 1 Geen paralyse 9
SNA16 1 Gedeeltelijke paralyse 0
SNA16 1 volledige paralyse 0
SNA16 2 Mutatie kwijt 0
SNA16 2 Geen paralyse 7
SNA16 2 Gedeeltelijke paralyse 3
SNA16 2 volledige paralyse 0
SNA16 3 Mutatie kwijt 1
SNA16 3 Geen paralyse 5
SNA16 3 Gedeeltelijke paralyse 3
SNA16 3 volledige paralyse 0
SNA25 1 Mutatie kwijt 2
SNA25 1 Geen paralyse 6
SNA25 1 Gedeeltelijke paralyse 1
SNA25 1 volledige paralyse 1
SNA25 2 Mutatie kwijt 0
SNA25 2 Geen paralyse 3
SNA25 2 Gedeeltelijke paralyse 5
SNA25 2 volledige paralyse 0
SNA25 3 Mutatie kwijt 4
SNA25 3 Geen paralyse 3
SNA25 3 Gedeeltelijke paralyse 3
SNA25 3 volledige paralyse 0

Maak een safe_kw functie, zodat er geen warning komt over eventuele ‘gepaarde metingen’.

safe_kw <- function(x, g) { #x: hoevaak is iets geteld. G: groeperings variabele, dus bijv. conditie
  if (length(unique(g)) < 2 #Als G minder dan 2 unieke groepen heeft kan er niks vergeleken worden
      || length(unique(x[!is.na(x)])) < 2 #Als x minder dan 2 unieke values heeft zal de test niet uitgevoerd worden. bijv. Alles is 0 -> geef NA
      ) { 
    return(NA_real_)
  }
  kruskal.test(x ~ g)$p.value #Runned de test waarbij de spreiding over X tussen groepen G wordt vergeleken. Geeft alleen de P-value terug.
}

x: numeric variable (e.g., counts for a particular outcome)

g: grouping variable (e.g., Conditie)

Run de kruskal test waarbij vergeleken wordt hoevaak een bepaald gedrag geobserveerd is per conditie.
Hieruit komen dan de P-values voor het gedrag vergeleken tussen de vier condities: WT16 vs WT25 vs SNA16 vs SNA25.

long_data |>
  group_by(Gedrag) |>
  summarise(
    p_value = safe_kw(Observaties, Conditie), 
    .groups = "drop")
## # A tibble: 4 × 2
##   Gedrag                 p_value
##   <chr>                    <dbl>
## 1 Gedeeltelijke paralyse   0.236
## 2 Geen paralyse            0.104
## 3 Mutatie kwijt            0.696
## 4 volledige paralyse       0.392

Omdat wij ook graag willen weten of er een significant verschil te vinden is tussen het WT en de SNA variant, is er gekeken naar de individuele p-values van de observaties op gedrag tussen de verschillende condities.

Dit is gedaan a.d.h.v. een pairwise wilcoxon rank-sum test ook wel genoemd Mann-Whitney U test.

Hieronder is er als eerst geprobeerd om alleen een Wilcox test uit te voeren tussen de condities SNA 25 en WT 25. Waarbij er gekeken werd naar of er een significant verschil zat tussen deze twee bij het vertonen van volledige paralyse.

pairwise_data <- long_data %>%
  filter(Gedrag == "volledige paralyse", Conditie %in% c("SNA25", "WT25"))

wilcox.test(Observaties ~ Conditie, data = pairwise_data, exact = FALSE) 
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  Observaties by Conditie
## W = 6, p-value = 0.505
## alternative hypothesis: true location shift is not equal to 0

exact = FALSE is hier gebruikt om een warning uit te zetten.
Deze warning zegt normaal gesproken iets over dat er een approximate p-value is uitgerekend ipv een exacte berekening.

De functie hieronder voert automatisch een wilcoxon rank-sum test (Mann-whitney U) uit voor elk mogelijk paar condities.
De functie zal teruggeven worden in een tabel waarin een overzicht staat van alle paren met p-waarden.

pairwise_wilcox <- function(df) { #Er wordt een nieuwe functie gemaakt waarbij de input een data frame is met de kolommen conditie en observaties en de output alle p-values zijn van de wilcoxon rank-sum test.
  condities <- unique(df$Conditie) #Hier wordt er een vector gemaakt voor de verschillende condities
  combs <- combn(condities, 2, simplify = FALSE) #Maak alle mogelijke paren van twee condities. simplify = FALSE zodat het een lijst blijft
  

  map_dfr(combs, function(pair) {#De map_dfr(combs, function(pair) zorgt ervoor dat alle paren in de lijst worden gevonden, waarbij map_dfr zelf de resultaten in één data frame zet.
    subset_df <- df %>% filter(Conditie %in% pair) #Alleen rijen worden geselecteerd waarbij de conditie gelijk is aan één van twee condities huidige paar.
    test_result <- wilcox.test(Observaties ~ Conditie, data = subset_df, exact = FALSE) #Hier wordt de wilcoxon test uitgevoerd.
    tibble( #Resultaten worden opgeslagen in een tibbe.
      Conditie1 = pair[1], #Eerste conditie van het paar bijv. WT25
      Conditie2 = pair[2], #Tweede conditie van het paar bijv. SNA25
      p_value = test_result$p.value
    )
  })
}

Je gebruikt de gemaakte pairwise_wilcox functie om nu voor alle combinaties de p-waarden te berekenen.

pairwise_results <- long_data %>%
  group_by(Gedrag) %>% #Groepeer per gedrag, omdat je het statistische verschil in gedrag tussen twee condities wil weten.
  group_modify(~ pairwise_wilcox(.x)) %>% #Voer voor elke groep de wilxoc test uit a.d.h.v. de eigen gemaakte functie
  ungroup() #Zorgt voor het verwijderen van groeperingen, zodat er een overzichtelijk dataframe uitkomt

Toon de resultaten

kable(pairwise_results)
Gedrag Conditie1 Conditie2 p_value
Gedeeltelijke paralyse WT16 WT25 0.1156880
Gedeeltelijke paralyse WT16 SNA16 0.3457786
Gedeeltelijke paralyse WT16 SNA25 0.1840386
Gedeeltelijke paralyse WT25 SNA16 0.6428348
Gedeeltelijke paralyse WT25 SNA25 1.0000000
Gedeeltelijke paralyse SNA16 SNA25 0.6428348
Geen paralyse WT16 WT25 0.1156880
Geen paralyse WT16 SNA16 0.8221868
Geen paralyse WT16 SNA25 0.0721982
Geen paralyse WT25 SNA16 0.2682859
Geen paralyse WT25 SNA25 1.0000000
Geen paralyse SNA16 SNA25 0.1840386
Mutatie kwijt WT16 WT25 0.4935628
Mutatie kwijt WT16 SNA16 0.8136637
Mutatie kwijt WT16 SNA25 0.8136637
Mutatie kwijt WT25 SNA16 0.5001843
Mutatie kwijt WT25 SNA25 0.8136637
Mutatie kwijt SNA16 SNA25 0.6530951
volledige paralyse WT16 WT25 NaN
volledige paralyse WT16 SNA16 NaN
volledige paralyse WT16 SNA25 0.5049851
volledige paralyse WT25 SNA16 NaN
volledige paralyse WT25 SNA25 0.5049851
volledige paralyse SNA16 SNA25 0.5049851