Analiza Urazowości Zapaśników Stylu Wolnego (2020–2025)
Autor
2026-02-09
0. Instalacja i wczytanie pakietów
Cel i Metodyka: Przygotowanie środowiska pracy. Wczytujemy pakiety niezbędne do analizy, w tym specjalistyczne biblioteki: *
exact2x2: do dokładnych testów Fishera (gdy liczebności są małe). *brglm2: do regresji Firtha (gdy występuje separacja danych, np. w Kadrze). *coin: do obliczania wielkości efektu (Z-score) w testach nieparametrycznych.
required_packages <- c(
"readxl", # wczytywanie plików .xlsx
"dplyr", # manipulacja danymi
"tidyr", # przekształcanie danych
"stringr", # operacje na tekstach
"ggplot2", # wykresy
"scales", # formatowanie osi wykresów
"forcats", # porządkowanie czynników (factor)
"gmodels", # CrossTable (tabele krzyżowe z chi-kwadrat)
"vcd", # testy asocjacji, mosaic plots
"exact2x2", # dokładny test Fishera (rozszerzony)
"dunn.test", # test post-hoc Dunna (po Kruskal-Wallis)
"effsize", # miary wielkości efektu
"coin", # Wilcoxon z prawdziwą statystyką Z (do wyliczenia r)
"brglm2", # regresja Firtha (separacja w regresji logistycznej)
"pwr", # analiza mocy statystycznej post-hoc
"RColorBrewer" # palety kolorów
)
# Zainstaluj brakujące pakiety
new_packages <- required_packages[!(required_packages %in% installed.packages()[, "Package"])]
if (length(new_packages)) install.packages(new_packages, repos = "https://cran.r-project.org")
invisible(lapply(required_packages, library, character.only = TRUE))
# Ustawienia globalne ggplot2
theme_set(theme_minimal(base_size = 13))
options(scipen = 999) # wyłącz notację naukową1. Wczytanie danych
Cel: Import surowych danych z pliku Excel. Weryfikacja poprawności wczytania (liczba wierszy i kolumn).
# Upewnij się, że plik analiza.xlsx znajduje się w tym samym katalogu co plik .Rmd
# Jeśli plik jest w innym miejscu, podaj pełną ścieżkę
df <- read_excel("analiza.xlsx")
cat("Wymiary zbioru danych:", nrow(df), "obserwacji,", ncol(df), "zmiennych\n")## Wymiary zbioru danych: 73 obserwacji, 26 zmiennych2. Czyszczenie i transformacja danych
Cel i Metodyka: Przygotowanie danych do analizy statystycznej. * Konwersja: Zamiana odpowiedzi ankietowych (tekst, np. “1-3 lata”) na liczby. * Faktoryzacja: Uporządkowanie zmiennych kategorialnych (np. Wykształcenie, Poziom) w odpowiedniej kolejności (levels), co jest kluczowe dla czytelności wykresów i wyników testów.
original_names <- names(df)
names(df) <- c(
"ID", "start_time", "completion_time", "email", "name", "last_modified",
"wiek", # Q1
"staz", # Q2
"kat_wagowa", # Q3
"poziom", # Q4
"treningi_tyg", # Q5
"uraz_tak_nie", # Q6
"ile_urazow", # Q7
"lokalizacja", # Q8
"rodzaj_uszkodzenia", # Q9
"absencja", # Q10
"konsultacja", # Q11
"leczenie", # Q12
"okolicznosci", # Q13
"przyczyna", # Q14
"faza_walki", # Q15
"mechanizm", # Q16
"robienie_wagi", # Q17
"kg_zrzucane", # Q18
"metody_redukcji", # Q19
"oslabienie_kontuzja" # Q20
)
# Zachowaj oryginalne etykiety jako atrybut
attr(df, "original_labels") <- setNames(original_names, names(df))
# --- 3a. Staż treningowy ---
df$staz_num <- as.numeric(str_extract(as.character(df$staz), "\\d+\\.?\\d*"))
# --- 3b. Kategoria wagowa ---
df$kat_wagowa_num <- as.numeric(str_extract(as.character(df$kat_wagowa), "\\d+"))
# --- 3c. Poprawka literówki ---
df$okolicznosci <- str_replace(df$okolicznosci, "sparigni", "sparingi")
# --- 3d. Zmienne porządkowe ---
df$treningi_tyg <- factor(df$treningi_tyg,
levels = c("1-3", "4-6", "7-9", "10-12", "Powyżej 12"),
ordered = TRUE)
df$treningi_tyg_num <- as.numeric(df$treningi_tyg)
df$ile_urazow <- factor(df$ile_urazow,
levels = c("1", "2-3", "4-5", "Powyżej 5"),
ordered = TRUE)
df$ile_urazow_num <- as.numeric(df$ile_urazow)
df$absencja <- factor(df$absencja,
levels = c("1-7 dni", "8-21 dni", "1-3 miesiące",
"Powyżej 3 miesięcy", "Uraz zakończył moją karierę"),
ordered = TRUE)
df$absencja_num <- as.numeric(df$absencja)
df$kg_zrzucane <- factor(df$kg_zrzucane,
levels = c("Mniej niż 2 kg", "2-4 kg", "5-7 kg", "Powyżej 7 kg"),
ordered = TRUE)
df$kg_zrzucane_num <- as.numeric(df$kg_zrzucane)
df$oslabienie_kontuzja <- factor(df$oslabienie_kontuzja,
levels = c("Nie", "Trudno powiedzieć", "Raczej tak", "Zdecydowanie tak"),
ordered = TRUE)
df$oslabienie_kontuzja_num <- as.numeric(df$oslabienie_kontuzja)
df$poziom <- factor(df$poziom,
levels = c("Amatorski/Rekreacyjny",
"Zawodniczy (starty w zawodach krajowych)",
"Kadra Narodowa / Międzynarodowy"),
ordered = TRUE)
df$poziom_num <- as.numeric(df$poziom)
# --- 3e. Robienie wagi ---
df$robienie_wagi <- factor(df$robienie_wagi,
levels = c("Nie", "Tak, sporadycznie", "Tak, regularnie"))
# --- 3f. Konwersja na numeryczne ---
df$wiek <- as.numeric(df$wiek)
# --- 3g. Podzbiór: tylko kontuzjowani ---
df_inj <- df %>% filter(uraz_tak_nie == "Tak")3. Statystyka opisowa
Cel: Charakterystyka badanej grupy. Pozwala ocenić strukturę demograficzną (wiek, staż, poziom) i sprawdzić, czy próba jest reprezentatywna. Te dane stanowią tło dla dalszych, bardziej zaawansowanych analiz.
cat("Liczba respondentów:", nrow(df), "\n")## Liczba respondentów: 73cat("Liczba kontuzjowanych:", nrow(df_inj), "\n")## Liczba kontuzjowanych: 61# Zmienne ciągłe
zmienne_ciagle <- c("wiek", "staz_num", "kat_wagowa_num")
etykiety <- c("Wiek", "Staż", "Waga")
opis_ciagle <- data.frame(
Zmienna = etykiety,
n = sapply(zmienne_ciagle, function(v) sum(!is.na(df[[v]]))),
Srednia = sapply(zmienne_ciagle, function(v) round(mean(df[[v]], na.rm = TRUE), 1)),
SD = sapply(zmienne_ciagle, function(v) round(sd(df[[v]], na.rm = TRUE), 1)),
Mediana = sapply(zmienne_ciagle, function(v) median(df[[v]], na.rm = TRUE)),
Min = sapply(zmienne_ciagle, function(v) min(df[[v]], na.rm = TRUE)),
Max = sapply(zmienne_ciagle, function(v) max(df[[v]], na.rm = TRUE))
)
print(opis_ciagle)## Zmienna n Srednia SD Mediana Min Max
## wiek Wiek 73 23.5 5.1 23 14 37
## staz_num Staż 73 9.4 6.4 8 1 27
## kat_wagowa_num Waga 73 77.5 13.1 75 48 125# Zmienne kategorialne - podsumowanie
cat_vars <- c("poziom", "treningi_tyg", "uraz_tak_nie", "robienie_wagi")
for (v in cat_vars) {
cat("\n---", v, "---\n")
print(table(df[[v]]))
}##
## --- poziom ---
##
## Amatorski/Rekreacyjny
## 16
## Zawodniczy (starty w zawodach krajowych)
## 39
## Kadra Narodowa / Międzynarodowy
## 18
##
## --- treningi_tyg ---
##
## 1-3 4-6 7-9 10-12 Powyżej 12
## 10 33 17 9 4
##
## --- uraz_tak_nie ---
##
## Nie Tak
## 12 61
##
## --- robienie_wagi ---
##
## Nie Tak, sporadycznie Tak, regularnie
## 22 27 244. Wykresy
Cel: Wizualizacja najważniejszy trendów. Wykresy pozwalają szybko dostrzec zależności, które później zweryfikujemy testami statystycznymi.
Definicja funkcji rysującej bar_plot, którą można
edytować w locie.
bar_plot <- function(data, var, title, fill_color = "#2c7fb8", flip = FALSE) {
p <- data %>%
filter(!is.na(!!sym(var))) %>%
count(!!sym(var)) %>%
mutate(pct = n / sum(n) * 100) %>%
ggplot(aes(x = reorder(!!sym(var), n), y = pct)) +
geom_col(fill = fill_color, width = 0.7) +
geom_text(aes(label = sprintf("%.1f%%\n(n=%d)", pct, n)),
hjust = -0.1, size = 3.5) +
labs(title = title, x = NULL, y = "Procent (%)") +
coord_flip() +
scale_y_continuous(expand = expansion(mult = c(0, 0.3)))
return(p)
}Wykres 1: Rozkład wieku
ggplot(df, aes(x = wiek)) +
geom_histogram(binwidth = 2, fill = "#2c7fb8", color = "white") +
labs(title = "Rozkład wieku respondentów", x = "Wiek (lata)", y = "Liczba") +
scale_x_continuous(breaks = seq(14, 40, 2))
Wykres 2: Poziom sportowy
bar_plot(df, "poziom", "Poziom sportowy respondentów", "#41b6c4")
Wykres 3: Urazy (Tak/Nie)
df %>%
count(uraz_tak_nie) %>%
mutate(pct = n / sum(n) * 100) %>%
ggplot(aes(x = uraz_tak_nie, y = pct, fill = uraz_tak_nie)) +
geom_col(width = 0.5) +
geom_text(aes(label = sprintf("%.1f%%\n(n=%d)", pct, n)), vjust = -0.5, size = 4) +
scale_fill_manual(values = c("Tak" = "#e31a1c", "Nie" = "#31a354")) +
labs(title = "Odsetek zapaśników z urazem", x = NULL, y = "Procent (%)") +
theme(legend.position = "none") +
scale_y_continuous(expand = expansion(mult = c(0, 0.3)))
Wykresy szczegółowe dla kontuzjowanych
p4 <- bar_plot(df_inj, "lokalizacja", "Lokalizacja najpoważniejszego urazu", "#e31a1c")
p5 <- bar_plot(df_inj, "absencja", "Czas absencji", "#fc8d59")
p6 <- bar_plot(df_inj, "faza_walki", "Faza walki", "#756bb1")
p7 <- bar_plot(df_inj, "mechanizm", "Mechanizm urazu", "#de2d26")
print(p4)
print(p5)
print(p6)
print(p7)
Wykresy dotyczące wagi
p8 <- bar_plot(df, "robienie_wagi", "Redukcja masy ciała", "#feb24c")
p9 <- bar_plot(df %>% filter(!is.na(oslabienie_kontuzja)),
"oslabienie_kontuzja",
"Czy osłabienie przyczyniło się do kontuzji?",
"#d95f0e")
print(p8)
print(p9)
Zaawansowane Wizualizacje
Cel: Głębsza analiza graficzna. * Wykres Mozaikowy (Mosaic Plot): Pokazuje nadreprezentację (kolor niebieski) lub niedoreprezentację (czerwony) w tabeli krzyżowej. Wielkość kafelka odpowiada liczebności grupy. * Boxplot z punktami (Jitter): Przy małych próbach (n=12 zdrowych) zwykły boxplot może mylić. Dodanie punktów (jitter) pokazuje każdą obserwację, co zapewnia pełną transparentność danych.
# 1. Mosaic Plot: Poziom Sportowy vs Uraz
# Używamy pakietu vcd do wizualizacji reszt (residuals)
# Skracamy etykiety poziomu dla czytelności (bo nachodzą na siebie)
df_mosaic <- df
levels(df_mosaic$poziom) <- c("Amatorski", "Zawodniczy", "Kadra")
cat("--- Mosaic Plot: Poziom Sportowy vs Uraz ---\n")## --- Mosaic Plot: Poziom Sportowy vs Uraz ---mosaic(~ poziom + uraz_tak_nie, data = df_mosaic,
shade = TRUE, # Kolorowanie reszt (nad/niedoreprezentacja)
legend = TRUE,
labeling_args = list(set_varnames = c(poziom="Poziom", uraz_tak_nie="Uraz")),
main = "Zależność: Poziom Sportowy a Wystąpienie Urazu")
# 2. Boxplot + Jitter: Wiek vs Uraz
cat("\n--- Boxplot + Jitter: Wiek a Uraz ---\n")##
## --- Boxplot + Jitter: Wiek a Uraz ---ggplot(df, aes(x = uraz_tak_nie, y = wiek, fill = uraz_tak_nie)) +
geom_boxplot(alpha = 0.6, outlier.shape = NA) + # Ukrywamy outliery boxplota, bo pokażemy je jako punkty
geom_jitter(width = 0.2, alpha = 0.8, size = 2) + # Prawdziwe punkty danych
stat_summary(fun = mean, geom = "point", shape = 18, size = 4, color = "black") + # Średnia jako romb
scale_fill_manual(values = c("Tak" = "#e31a1c", "Nie" = "#31a354")) +
labs(title = "Rozkład wieku w grupach z urazem i bez",
subtitle = "Kropki = poszczególni zawodnicy; Romb = średnia",
x = "Czy wystąpił uraz?", y = "Wiek (lata)") +
theme(legend.position = "none")
# 3. Boxplot + Jitter: Staż vs Uraz
cat("\n--- Boxplot + Jitter: Staż a Uraz ---\n")##
## --- Boxplot + Jitter: Staż a Uraz ---ggplot(df, aes(x = uraz_tak_nie, y = staz_num, fill = uraz_tak_nie)) +
geom_boxplot(alpha = 0.6, outlier.shape = NA) +
geom_jitter(width = 0.2, alpha = 0.8, size = 2) +
stat_summary(fun = mean, geom = "point", shape = 18, size = 4, color = "black") +
scale_fill_manual(values = c("Tak" = "#e31a1c", "Nie" = "#31a354")) +
labs(title = "Staż treningowy w grupach z urazem i bez",
subtitle = "Kropki = poszczególni zawodnicy; Romb = średnia",
x = "Czy wystąpił uraz?", y = "Staż (lata)") +
theme(legend.position = "none")
5. Testy Statystyczne
Cel i Metodyka: Weryfikacja hipotez badawczych. * Test asocjacji: Domyślnie używamy testu Fishera (zamiast Chi-kwadrat), ponieważ w wielu podgrupach liczebności są bardzo małe (<5). Test Fishera jest w takiej sytuacji “złotym standardem”. * Test różnic: Dla zmiennych ciągłych (wiek, staż) stosujemy test U Manna-Whitneya (nieparametryczny), ponieważ grupy są bardzo nierówne (61 vs 12) i rozkłady mogą odbiegać od normalnego (test t-Studenta byłby ryzykowny). * Test Kruskala-Wallisa: Odpowiednik ANOVA dla wielu grup (np. porównanie absencji między 3 poziomami sportowymi).
Testy asocjacji (Chi-kwadrat / Fisher)
cat("=== TEST 1: Poziom sportowy x Uraz ===\n")## === TEST 1: Poziom sportowy x Uraz ===tab1 <- table(df$poziom, df$uraz_tak_nie)
print(tab1)##
## Nie Tak
## Amatorski/Rekreacyjny 4 12
## Zawodniczy (starty w zawodach krajowych) 8 31
## Kadra Narodowa / Międzynarodowy 0 18# Jeśli wartości oczekiwane < 5, użyj fisher.test
fisher.test(tab1)##
## Fisher's Exact Test for Count Data
##
## data: tab1
## p-value = 0.05961
## alternative hypothesis: two.sidedcat("\n=== TEST 2: Treningi x Uraz ===\n")##
## === TEST 2: Treningi x Uraz ===tab2 <- table(df$treningi_tyg, df$uraz_tak_nie)
tryCatch(chisq.test(tab2), warning=function(w) fisher.test(tab2, simulate.p.value=TRUE))##
## Fisher's Exact Test for Count Data with simulated p-value (based on
## 2000 replicates)
##
## data: tab2
## p-value = 0.2004
## alternative hypothesis: two.sidedcat("\n=== TEST 3: Robienie wagi x Uraz ===\n")##
## === TEST 3: Robienie wagi x Uraz ===tab3 <- table(df$robienie_wagi, df$uraz_tak_nie)
print(tab3)##
## Nie Tak
## Nie 4 18
## Tak, sporadycznie 7 20
## Tak, regularnie 1 23fisher.test(tab3)##
## Fisher's Exact Test for Count Data
##
## data: tab3
## p-value = 0.1072
## alternative hypothesis: two.sidedTesty różnic w grupach (Mann-Whitney / Kruskal-Wallis)
cat("=== TEST: Wiek (Kontuzjowani vs Nie) ===\n")## === TEST: Wiek (Kontuzjowani vs Nie) ===wilcox.test(wiek ~ uraz_tak_nie, data = df)##
## Wilcoxon rank sum test with continuity correction
##
## data: wiek by uraz_tak_nie
## W = 346.5, p-value = 0.7762
## alternative hypothesis: true location shift is not equal to 0cat("\n=== TEST: Staż (Kontuzjowani vs Nie) ===\n")##
## === TEST: Staż (Kontuzjowani vs Nie) ===wilcox.test(staz_num ~ uraz_tak_nie, data = df)##
## Wilcoxon rank sum test with continuity correction
##
## data: staz_num by uraz_tak_nie
## W = 223.5, p-value = 0.03423
## alternative hypothesis: true location shift is not equal to 0cat("\n=== TEST KW: Absencja vs Poziom ===\n")##
## === TEST KW: Absencja vs Poziom ===kruskal.test(absencja_num ~ poziom, data = df_inj)##
## Kruskal-Wallis rank sum test
##
## data: absencja_num by poziom
## Kruskal-Wallis chi-squared = 5.2362, df = 2, p-value = 0.07294cat("\n=== TEST KW: Liczba urazów vs Poziom ===\n")##
## === TEST KW: Liczba urazów vs Poziom ===kruskal.test(ile_urazow_num ~ poziom, data = df_inj)##
## Kruskal-Wallis rank sum test
##
## data: ile_urazow_num by poziom
## Kruskal-Wallis chi-squared = 4.9227, df = 2, p-value = 0.08532Korelacje Spearmana
cat("Wiek x Liczba urazów:\n")## Wiek x Liczba urazów:print(cor.test(df_inj$wiek, df_inj$ile_urazow_num, method = "spearman"))##
## Spearman's rank correlation rho
##
## data: df_inj$wiek and df_inj$ile_urazow_num
## S = 34452, p-value = 0.4949
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.08905883cat("Staż x Liczba urazów:\n")## Staż x Liczba urazów:print(cor.test(df_inj$staz_num, df_inj$ile_urazow_num, method = "spearman"))##
## Spearman's rank correlation rho
##
## data: df_inj$staz_num and df_inj$ile_urazow_num
## S = 27707, p-value = 0.03721
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.2674083cat("Treningi x Liczba urazów:\n")## Treningi x Liczba urazów:print(cor.test(df_inj$treningi_tyg_num, df_inj$ile_urazow_num, method = "spearman"))##
## Spearman's rank correlation rho
##
## data: df_inj$treningi_tyg_num and df_inj$ile_urazow_num
## S = 29055, p-value = 0.07231
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.23175296. Regresja Logistyczna (Firth)
Cel i Metodyka: Analiza wieloczynnikowa – sprawdzamy, które czynniki (wiek, staż, poziom) zwiększają ryzyko urazu, biorąc pod uwagę wzajemne powiązania. * Dlaczego metoda Firtha? Używamy pakietu
brglm2, ponieważ w danych występuje zjawisko separacji (quasi-complete separation). Przykładowo: 100% zawodników “Kadry Narodowej” miało uraz. Tradycyjna regresja logistyczna wyrzuciłaby błąd (iloraz szans dążący do nieskończoności). Metoda Firtha koryguje to obciążenie, dając wiarygodne wyniki nawet przy takich danych.
Stosujemy metodę Firtha ze względu na możliwą separację danych (np. wszyscy kadrowicze mają kontuzję).
df_model <- df %>%
mutate(
uraz_bin = ifelse(uraz_tak_nie == "Tak", 1, 0),
poziom_f = relevel(factor(as.character(poziom), ordered = FALSE), ref = "Amatorski/Rekreacyjny")
) %>%
filter(!is.na(staz_num), !is.na(treningi_tyg_num))
model1 <- glm(uraz_bin ~ wiek + staz_num + poziom_f + treningi_tyg_num + robienie_wagi,
data = df_model, family = binomial,
method = "brglmFit", type = "AS_mean")
summary(model1)##
## Call:
## glm(formula = uraz_bin ~ wiek + staz_num + poziom_f + treningi_tyg_num +
## robienie_wagi, family = binomial, data = df_model, method = "brglmFit",
## type = "AS_mean")
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.1550 0.1780 0.4507 0.6909 1.3622
##
## Coefficients:
## Estimate Std. Error z value
## (Intercept) 2.62058 2.19258 1.195
## wiek -0.12159 0.09587 -1.268
## staz_num 0.08344 0.09585 0.871
## poziom_fKadra Narodowa / Międzynarodowy 1.23796 1.95078 0.635
## poziom_fZawodniczy (starty w zawodach krajowych) 0.15656 0.86424 0.181
## treningi_tyg_num 0.51241 0.46554 1.101
## robienie_wagiTak, sporadycznie -1.16780 0.82137 -1.422
## robienie_wagiTak, regularnie 0.22922 1.12679 0.203
## Pr(>|z|)
## (Intercept) 0.232
## wiek 0.205
## staz_num 0.384
## poziom_fKadra Narodowa / Międzynarodowy 0.526
## poziom_fZawodniczy (starty w zawodach krajowych) 0.856
## treningi_tyg_num 0.271
## robienie_wagiTak, sporadycznie 0.155
## robienie_wagiTak, regularnie 0.839
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 65.253 on 72 degrees of freedom
## Residual deviance: 51.807 on 65 degrees of freedom
## AIC: 67.807
##
## Type of estimator: AS_mean (mean bias-reducing adjusted score equations)
## Number of Fisher Scoring iterations: 10cat("\nIlorazy szans (OR):\n")##
## Ilorazy szans (OR):print(exp(coef(model1)))## (Intercept)
## 13.7436269
## wiek
## 0.8855112
## staz_num
## 1.0870188
## poziom_fKadra Narodowa / Międzynarodowy
## 3.4485639
## poziom_fZawodniczy (starty w zawodach krajowych)
## 1.1694799
## treningi_tyg_num
## 1.6693136
## robienie_wagiTak, sporadycznie
## 0.3110502
## robienie_wagiTak, regularnie
## 1.25761687. Podsumowanie wyników z korektą Holma
Cel i Metodyka: Zbiorcze zestawienie wszystkich wykonanych testów. * Korekta Holma: Ponieważ wykonaliśmy aż 19 testów, istnieje ryzyko, że któryś wynik wyjdzie “istotny” przez czysty przypadek. Korekta Holma (dla rodzin hipotez) to rygorystyczne podejście, które chroni przed tzw. p-hackingiem i zwiększa wiarygodność naukową pracy.
# Funkcja pomocnicza: V Cramera
cramers_v <- function(tab) {
chi2 <- suppressWarnings(chisq.test(tab))$statistic
n <- sum(tab)
k <- min(nrow(tab), ncol(tab))
v <- sqrt(as.numeric(chi2) / (n * (k - 1)))
return(ifelse(is.nan(v) | is.na(v), NA, round(v, 3)))
}
# --- Odtworzenie testów i pobranie p-wartości ---
# Rodzina A: Predyktory (Chi2/Fisher/U-test)
test1 <- fisher.test(table(df$poziom, df$uraz_tak_nie))
test2 <- fisher.test(table(df$treningi_tyg, df$uraz_tak_nie), simulate.p.value = TRUE, B = 10000)
test3 <- fisher.test(table(df$robienie_wagi, df$uraz_tak_nie))
test6 <- wilcox.test(wiek ~ uraz_tak_nie, data = df, exact = FALSE)
test7 <- wilcox.test(staz_num ~ uraz_tak_nie, data = df, exact = FALSE)
# Rodzina B: Korelacje (Spearman)
cor_a <- cor.test(df_inj$wiek, df_inj$ile_urazow_num, method = "spearman", exact = FALSE)
cor_b <- cor.test(df_inj$staz_num, df_inj$ile_urazow_num, method = "spearman", exact = FALSE)
cor_c <- cor.test(df_inj$treningi_tyg_num, df_inj$ile_urazow_num, method = "spearman", exact = FALSE)
cor_d <- cor.test(df_inj$staz_num, df_inj$absencja_num, method = "spearman", exact = FALSE)
# cor_e requires extra filtering
df_waga <- df %>% filter(!is.na(kg_zrzucane_num), !is.na(oslabienie_kontuzja_num))
if(nrow(df_waga) > 2) {
cor_e <- cor.test(df_waga$kg_zrzucane_num, df_waga$oslabienie_kontuzja_num, method = "spearman", exact = FALSE)
p_cor_e <- cor_e$p.value
} else {
p_cor_e <- NA
}
# Rodzina C: Okoliczności/Mechanizm
df_inj_sub <- df_inj %>% filter(!is.na(okolicznosci), !is.na(lokalizacja))
test4 <- fisher.test(table(df_inj_sub$okolicznosci, df_inj_sub$lokalizacja), simulate.p.value = TRUE, B = 10000)
test5 <- fisher.test(table(df_inj$faza_walki, df_inj$mechanizm), simulate.p.value = TRUE, B = 10000)
test11 <- fisher.test(table(df_inj$okolicznosci, df_inj$faza_walki), simulate.p.value = TRUE, B = 10000)
df_faza <- df_inj %>% filter(!is.na(faza_walki))
# Check if distinct faza_walki values exist > 1 group
if(length(unique(df_faza$faza_walki)) > 1) {
test13 <- kruskal.test(wiek ~ faza_walki, data = df_faza)
p_test13 <- test13$p.value
} else {
p_test13 <- NA
}
# Rodzina D: Leczenie i konsultacje
df_abs <- df_inj %>% filter(!is.na(absencja_num), !is.na(poziom))
test8 <- kruskal.test(absencja_num ~ poziom, data = df_abs)
df_inj_ku <- df_inj %>% filter(!is.na(ile_urazow_num), !is.na(poziom))
test9 <- kruskal.test(ile_urazow_num ~ poziom, data = df_inj_ku)
df_kons <- df_inj %>% filter(!is.na(konsultacja), !is.na(poziom))
test14 <- fisher.test(table(df_kons$poziom, df_kons$konsultacja))
df_lecz <- df_inj %>% filter(!is.na(leczenie), !is.na(poziom))
test15 <- fisher.test(table(df_lecz$poziom, df_lecz$leczenie), simulate.p.value = TRUE, B = 10000)
# Rodzina E: Redukcja
df_waga2 <- df %>% filter(!is.na(robienie_wagi), !is.na(oslabienie_kontuzja))
test12 <- fisher.test(table(df_waga2$robienie_wagi, df_waga2$oslabienie_kontuzja), simulate.p.value = TRUE, B = 10000)
# Zbiór P-values
p_raw <- c(
test1$p.value, test2$p.value, test3$p.value, test4$p.value,
test5$p.value, test6$p.value, test7$p.value, test8$p.value,
test9$p.value, cor_a$p.value, cor_b$p.value, cor_c$p.value,
cor_d$p.value, p_cor_e, test11$p.value, test12$p.value,
p_test13, test14$p.value, test15$p.value
)
rodzina <- c("A", "A", "A", "C", "C", "A", "A", "D", "D",
"B", "B", "B", "B", "B",
"C", "E", "C", "D", "D")
test_names <- c(
"1. Poziom sportowy x Uraz", "2. Treningi/tyg x Uraz", "3. Robienie wagi x Uraz",
"4. Okoliczności x Lokalizacja", "5. Faza walki x Mechanizm",
"6. Wiek: kontuzj. vs nie", "7. Staż: kontuzj. vs nie",
"8. Absencja x Poziom", "9. Liczba urazów x Poziom",
"10a. rho Wiek x Ile urazów", "10b. rho Staż x Ile urazów",
"10c. rho Treningi x Ile urazów", "10d. rho Staż x Absencja",
"10e. rho Kg x Osłabienie",
"11. Okoliczności x Faza walki", "12. Rob.wagi x Osłab.->kont.",
"13. Wiek x Faza walki", "14. Konsultacja x Poziom", "15. Leczenie x Poziom"
)
# Korekta Holma
p_holm_family <- numeric(length(p_raw))
for (fam in unique(rodzina)) {
idx <- which(rodzina == fam)
p_holm_family[idx] <- p.adjust(p_raw[idx], method = "holm")
}
p_holm_global <- p.adjust(p_raw, method = "holm")
# Funkcja formatująca p-wartości
p_val_format <- function(p) {
sapply(p, function(x) {
if(is.na(x)) return("-")
if(x < 0.001) return("<.001")
return(sprintf("%.4f", x))
})
}
# Tabela wyników
results_summary <- data.frame(
Rodzina = rodzina,
Test = test_names,
p_surowe = p_val_format(p_raw),
p_Holm_rodz = paste0(p_val_format(p_holm_family), ifelse(!is.na(p_holm_family) & p_holm_family < 0.05, " *", "")),
p_Holm_global = paste0(p_val_format(p_holm_global), ifelse(!is.na(p_holm_global) & p_holm_global < 0.05, " *", "")),
stringsAsFactors = FALSE
)
results_summary <- results_summary[order(results_summary$Rodzina), ]
knitr::kable(results_summary, caption = "Podsumowanie wyników testów z korektą Holma", align = "llccc")| Rodzina | Test | p_surowe | p_Holm_rodz | p_Holm_global | |
|---|---|---|---|---|---|
| 1 | A | 1. Poziom sportowy x Uraz | 0.0596 | 0.2385 | 0.6557 |
| 2 | A | 2. Treningi/tyg x Uraz | 0.1841 | 0.3682 | 0.9204 |
| 3 | A | 3. Robienie wagi x Uraz | 0.1072 | 0.3216 | 0.7504 |
| 6 | A | 6. Wiek: kontuzj. vs nie | 0.7762 | 0.7762 | 1.0000 |
| 7 | A | 7. Staż: kontuzj. vs nie | 0.0342 | 0.1711 | 0.5549 |
| 10 | B | 10a. rho Wiek x Ile urazów | 0.4949 | 0.4949 | 1.0000 |
| 11 | B | 10b. rho Staż x Ile urazów | 0.0372 | 0.1306 | 0.5549 |
| 12 | B | 10c. rho Treningi x Ile urazów | 0.0723 | 0.1446 | 0.7231 |
| 13 | B | 10d. rho Staż x Absencja | 0.0326 | 0.1306 | 0.5549 |
| 14 | B | 10e. rho Kg x Osłabienie | 0.0014 | 0.0071 * | 0.0270 * |
| 4 | C | 4. Okoliczności x Lokalizacja | 0.4514 | 0.9027 | 1.0000 |
| 5 | C | 5. Faza walki x Mechanizm | 0.0018 | 0.0072 * | 0.0324 * |
| 15 | C | 11. Okoliczności x Faza walki | 0.0340 | 0.1020 | 0.5549 |
| 17 | C | 13. Wiek x Faza walki | 0.4939 | 0.9027 | 1.0000 |
| 8 | D | 8. Absencja x Poziom | 0.0729 | 0.2188 | 0.7231 |
| 9 | D | 9. Liczba urazów x Poziom | 0.0853 | 0.2188 | 0.7231 |
| 18 | D | 14. Konsultacja x Poziom | 0.0496 | 0.1985 | 0.6452 |
| 19 | D | 15. Leczenie x Poziom | 0.1134 | 0.2188 | 0.7504 |
| 16 | E | 12. Rob.wagi x Osłab.->kont. | 0.0531 | 0.0531 | 0.6452 |
cat("\nRodziny hipotez:\n")##
## Rodziny hipotez:cat(" A – Predyktory wystąpienia urazu\n")## A – Predyktory wystąpienia urazucat(" B – Korelacje rang Spearmana\n")## B – Korelacje rang Spearmanacat(" C – Okoliczności i mechanizm urazu\n")## C – Okoliczności i mechanizm urazucat(" D – Absencja, leczenie, konsultacje wg poziomu\n")## D – Absencja, leczenie, konsultacje wg poziomucat(" E – Redukcja masy ciała a percepcja\n")## E – Redukcja masy ciała a percepcja8. Analiza Wrażliwości (Sensitivity Power Analysis)
Cel: Odpowiedź na pytanie: “Dlaczego tak mało wyników jest istotnych statystycznie?”. * Wnioski: Analiza pokazuje, że przy tak małej grupie kontrolnej (n=12) nasze badanie ma moc wykrycia tylko bardzo silnych zależności. Brak istotności statystycznej (p>0.05) w wielu miejscach nie oznacza, że zależności nie ma – oznacza tylko, że przy tej wielkości próby nie udało się jej “udowodnić” matematycznie. To bardzo ważny argument do obrony pracy.
cat("Pytanie: jaki MINIMALNY efekt jesteśmy w stanie wykryć przy mocy 80%?\n")## Pytanie: jaki MINIMALNY efekt jesteśmy w stanie wykryć przy mocy 80%?# a) Moc dla chi-kwadrat
# Zakładamy N = nrow(df) dla testów na całej próbie
pwr_chi <- pwr.chisq.test(w = NULL, N = nrow(df), df = 2, sig.level = 0.05, power = 0.80)
cat(sprintf("a) Chi^2 (df=2): minimalny wykrywalny efekt w = %.3f\n", pwr_chi$w))## a) Chi^2 (df=2): minimalny wykrywalny efekt w = 0.363cat("(0.1=mały, 0.3=średni, 0.5=duży)\n\n")## (0.1=mały, 0.3=średni, 0.5=duży)# b) Moc dla Mann-Whitney
# Liczebności grup
n1 <- sum(df$uraz_tak_nie == "Tak")
n2 <- sum(df$uraz_tak_nie == "Nie")
# Upewniamy się, że są dwie grupy
if(n1 > 0 && n2 > 0) {
n_harm <- 2 * n1 * n2 / (n1 + n2)
pwr_t <- pwr.t.test(n = n_harm / 2, d = NULL, sig.level = 0.05, power = 0.80,
type = "two.sample", alternative = "two.sided")
cat(sprintf("b) Mann-Whitney (przybliżenie t): minimalny wykrywalny efekt d = %.3f\n", pwr_t$d))
cat("(0.2=mały, 0.5=średni, 0.8=duży)\n\n")
} else {
cat("b) Mann-Whitney: brak wystarczającej liczby grup do analizy mocy.\n\n")
}## b) Mann-Whitney (przybliżenie t): minimalny wykrywalny efekt d = 1.323
## (0.2=mały, 0.5=średni, 0.8=duży)# c) Moc dla korelacji Spearmana
pwr_r <- pwr.r.test(n = nrow(df_inj), r = NULL, sig.level = 0.05, power = 0.80)
cat(sprintf("c) Korelacja (n=%d): minimalny wykrywalny r = %.3f\n", nrow(df_inj), pwr_r$r))## c) Korelacja (n=61): minimalny wykrywalny r = 0.350cat("(0.1=mały, 0.3=średni, 0.5=duży)\n")## (0.1=mały, 0.3=średni, 0.5=duży)9. Dodatek Metodologiczny: Uzasadnienie Wyboru Testów i Metodyka
Ta sekcja wyjaśnia zastosowane metody statystyczne. Stanowi ściągawkę do obrony pracy magisterskiej, uzasadniającą dlaczego użyto zaawansowanych technik (np. testu Fishera czy regresji Firtha) zamiast standardowych rozwiązań.
1. Dlaczego Test Dokładny Fishera zamiast Chi-kwadrat?
W standardowej analizie tabel krzyżowych stosuje się test Chi-kwadrat Pearsona. Jednak wymaga on, aby liczebności oczekiwane w komórkach tabeli były odpowiednio duże (zazwyczaj > 5).
- Problem: Nasza grupa kontrolna (osoby bez urazu) jest bardzo mała (n=12). W rezultacie wiele tabel ma liczebności oczekiwane poniżej 5 (a nawet 1).
- Ryzyko: Użycie zwykłego testu Chi-kwadrat w takiej sytuacji prowadziłoby do błędnych wyników (fałszywie istotnych lub nieistotnych) – test ten “łamie się” przy małych liczbach.
- Rozwiązanie: Zastosowano Test Dokładny Fishera (Fisher’s Exact Test). Jest to “złoty standard” przy małych próbach, ponieważ oblicza dokładne prawdopodobieństwo wystąpienia danej konfiguracji danych, a nie polega na przybliżeniach asymptotycznych. Dla większych tabel (np. 4x5) zastosowano wersję z symulacją Monte Carlo, aby uzyskać wiarygodne p-wartości.
2. Dlaczego Regresja Logistyczna Firtha zamiast zwykłej?
- Problem: Wystąpiło zjawisko tzw. separacji
danych (complete or quasi-complete separation). Przykład: 100%
zawodników “Kadry Narodowej” doznało urazu. Zwykła regresja logistyczna
(
glm) nie radzi sobie z taką sytuacją – algorytm nie może zbiec, a błędy standardowe dążą do nieskończoności. - Rozwiązanie: Zastosowano metodę
Logistycznej Regresji Firtha (pakiet
brglm2). Metoda ta wprowadza “karę” (bias reduction) do funkcji wiarygodności, co pozwala na poprawne oszacowanie modelu nawet w przypadku idealnej separacji. - Wynik: Model jest stabilny, a oszacowane ilorazy szans (OR) są wiarygodne.
3. Interpretacja braku istotności statystycznej (p > 0.05)
Wiele wyników w tej pracy ma p-wartość powyżej 0.05 (brak istotności statystycznej).
- Przyczyna: Niekoniecznie oznacza to brak zależności w rzeczywistości. Głównym powodem jest mała liczebność próby (szczególnie grupy bez kontuzji).
- Analiza Mocy (Sensitivity Power Analysis): Wykazała, że przy tej liczebności próby badanie ma moc wykrycia jedynie dużych i bardzo dużych efektów. Efekty słabe i średnie mogą pozostać niewykryte (błąd II rodzaju).
- Wniosek: Brak statystycznej istotności nie dowodzi braku wpływu, ale wskazuje, że przy obecnej wielkości próby nie można go jednoznacznie potwierdzić. Trendy widoczne na wykresach są nadal cenne interpretacyjnie.
4. Korekta Holma dla wielokrotnych porównań
Wykonano w pracy kilkanaście testów statystycznych.
- Ryzyko: Przy wykonywaniu wielu testów rośnie ryzyko, że któryś wyjdzie istotny “przez przypadek” (problem wielokrotnych porównań).
- Rozwiązanie: Zastosowano poprawkę Holma (Holm-Bonferroni method), która jest standardem w rzetelnych badaniach naukowych. Jest ona bezpieczna (kontroluje błąd I rodzaju) i silniejsza (ma większą moc) niż klasyczna poprawka Bonferroniego.
5. Współczynnik Pseudo-R² (McFadden)
W modelu regresji uzyskano Pseudo-R² na poziomie ok. 21%. W naukach społecznych i sporcie, gdzie na wynik wpływa ogromna liczba czynników (technika, genetyka, dieta, psychika, przypadek), jest to wynik satysfakcjonujący, wskazujący, że model ma umiarkowaną siłę wyjaśniającą, choć oczywiście istnieje wiele innych, niemierzonych w ankiecie zmiennych wpływających na ryzyko urazu.