7 Neparametrija – Podatkovna znanost: R u znanosti i obrazovanju

7.1 Obrada i prikaz neparametrijskih podataka

Neparametrijska statistika koristi se kod podataka koji imaju neke od sljedećih karakteristika:

mali broj ispitanika ili nekih drugih entiteta u istraživanju (n<20 ili n<30, treba biti oprezan pri korištenju točnih brojeva jer navedene vrijednosti ovise o predmetu istraživanja i primijenjenoj metodologiji)
podaci su izraženi na nominalnim ili ordinalnim mjernim ljestvicama
raspodjela se značajno razlikuje od normalne (U-raspodjela, multimodalna raspodjela, Poissonova raspodjela)

Neke karakteristike korištenja neparametrijske statistike mogu se sažeti na slijedeći način:

Može se koristiti za sve vrste ljestvica
Jednostavan izračun
Ne traži populacijske podatke (nije i razina zaključivanja) #Nejasno je što se misli…
Rezultati mogu izgledati gotovo jednako kao i kod parametrijske statistike
U društvenim znanostima, pri ispitivanju stavova, ekstremne vrijednosti – često je obrnuta normalna raspodjela (U)

No, postoje i određena ograničenja u odnosu na druge procedure:

Sadrži nepotpune modele / podatke
Višestruke usporedbe – post hoc testovi
Snaga zaključivanja je manja od neparametrijskih testova

Snaga nekog testa odnosi se na vjerojatnost odbacivanja nul-hipoteze ako je ta hipoteza zaista pogrešna. Sposobnost otkrivanja razlike ako ta razlika zaista postoji.

U slijedećem tabličnom prikazu možemo vidjeti usporedbe između nekih svojstava neparametrijske i parametrijske statistike.

parametrijski		neparametrijski
Intervalna, omjerna	Ljestvica mjerenja	Nominalna, ordinalna
Normalna ili približno normalna	Oblik raspodjele	Nije relevantno
Veća	Statistička snaga	Manja
Podjednake	Varijance	Nije relevantno
Veća	Osjetljivost	Manja
Veća	Razumljivost za korisnike	Manja

Kako provjeriti normalnost raspodjele dobivenih vrijednosti nekog istraživanja? Jedan od načina je prikaz raspodjele i odstupanje od normalne, zamišljene raspodjele. Drugi način je pomoću Q-Q prikaza.

Tako npr. u R sučelju možemo simulirati normalnu raspodjelu s N=100 te prikazati pomoću histograma i Q-Q slikovnog prikaza raspodjelu podataka.

normal.results <- rnorm (n=100)
hist(x=normal.results)

qqnorm(y=normal.results)

Ukoliko Q-Q slikovni prkaz prati zamišljnu linearnu (dijagnonalnu) liniju tada je dobivena raspodjela u okvirima normalne raspodjele.

Pomoću Shapiro-Wilk testa možemo numerički provjeriti da li navedena raspodjela podataka odstupa od normalne.

shapiro.test( x = normal.results )


    Shapiro-Wilk normality test

data:  normal.results
W = 0.98419, p-value = 0.277

Kako vidimo p vrijednost je veća od 0.05 te je statistički ne značajna razini od 95%. U tom slučaju smatramo da se oblik distribucije ne razlikuje značajno od normalne.

U sljedećim dijelovima ćemo prikazati neparametrijske postupke ovisno o vrsti podataka koje imamo te specifičnim pitanjima na koja želimo odgovoriti.

7.2 Jedan uzorak - test hipoteze

Primjena neparametrijskog testa za jedan uzorak može se pronaći od znanosti do obrazovanja. Tako, primjerice želimo vidjeti da li je profesor ili ocjenjivač sklon davanju ocjena koje su različite od sredine tj. ocjene 3. Kako bi provjerili hipotezu, možemo se poslužiti primjerom sljedećih ocjena profesora Petra.

PodaciOcjene =("
  Profesor  Ocjena
 'Petar'  3
 'Petar'  4
 'Petar'  5
 'Petar'  4
 'Petar'  4
 'Petar'  5
 'Petar'  5
 'Petar'  3
 'Petar'  2
 'Petar'  5     
")

PodaciOcjeneData = read.table(textConnection(PodaciOcjene),header=TRUE)

#PodaciOcjeneData$Ocjena = factor(PodaciOcjeneData$Ocjena,
#                       ordered = TRUE)

xtabs( ~ Profesor + Ocjena,
       data = PodaciOcjeneData)

        Ocjena
Profesor 2 3 4 5
   Petar 1 2 3 4

Navedeni primjer napravljen je prema uzoru poznatog paketa rcompanion (Mangiafico, 2023) te online knjige (Mangiafico2016?).

U slijedećem koraku napravit ćemo križanje tablice (crosstabulation) koje pokazuje učestalosti i postotke.

xtabs( ~ Profesor + Ocjena,
       data = PodaciOcjeneData)

        Ocjena
Profesor 2 3 4 5
   Petar 1 2 3 4

#proporcija u slijedećem koraku
#prvo dodjelimo podatke varijabli i zatim prikažemo proporcije

proporcije = xtabs( ~ Profesor + Ocjena,
       data = PodaciOcjeneData)

prop.table(proporcije,margin = 1)

        Ocjena
Profesor   2   3   4   5
   Petar 0.1 0.2 0.3 0.4

Slikovni prikaz možemo napraviti također koristeći bazične barplot() funkcije.

bar_ocjena = xtabs(~ Ocjena,
          data=PodaciOcjeneData)

barplot(bar_ocjena,  
        col="dark gray",
        xlab="Petar - ocjene",
        ylab="Učestalost")

Iz slikovnog prikaza je vidljivo kako je raspodjela negativno asimetrična tj. veći broj vrijednosti ocjena je pomaknut prema vrijednosti 5. Deskriptivna statistika može se prikazati i pomoću FSA paketa na slijedeći način.

library(FSA)

## FSA v0.9.5. See citation('FSA') if used in publication.
## Run fishR() for related website and fishR('IFAR') for related book.

Summarize(Ocjena ~ Profesor,                                                     
          data=PodaciOcjeneData,
          digits=3)

  Profesor  n mean    sd min   Q1 median Q3 max
1    Petar 10    4 1.054   2 3.25      4  5   5

Vratimo se na početak primjera jer trebamo odgovoriti na pitanje da li je profesor ili ocjenjivač sklon davanju većih ocjena tj. trebamo ispitati razlikuju li se ocjenjivanje statistički značajno od 3? Za dobiti odgovor na to pitanje, osim slikovnog prikaza, potrebno je provesti test - Wilcoxon tekst predznaka (Wilcoxon signed-rank test).

Wilcoxon test pozivamo funkcijom wilcox.test() gdje je u argumentima funkcije potrebno odrediti vrijednost “mu” tj. referentnu vrijednost od koje se određuje odstupanje distribucije ocjena. U našem primjeru to je vrijednost 3.

#prije primjene Wilcox testa - pretvaranje varijable u čistu numeričku varijablu
PodaciOcjeneData$Ocjena <- as.numeric(PodaciOcjeneData$Ocjena)
wilcox.test(PodaciOcjeneData$Ocjena,
            mu=3,
            conf.int=TRUE,
            conf.level=0.95)


    Wilcoxon signed rank test with continuity correction

data:  PodaciOcjeneData$Ocjena
V = 33.5, p-value = 0.03125
alternative hypothesis: true location is not equal to 3
90 percent confidence interval:
 3.49997 5.00000
sample estimates:
(pseudo)median 
      4.499989

Iz rezultata primjene Wilcoxon testa vidimo i raspon intervala pouzdanosti. Wilcoxon test je značajan na razini od 95%. Pomoću rcompanion paketa može se izračunati i efekt učinka za Wilcoxon test (Mangiafico, 2023).

library(rcompanion)
wilcoxonOneSampleR(PodaciOcjeneData$Ocjena,mu=3)

    r 
0.788

7.3 Dva nezavisna uzorka

Dva vrlo često korištena testa za usporedbu dva nezavisna uzorka su medijan test i test sume rangova (Wilcoxon t test, Mann-Whitney U test).

Medijan test:

Svodi se na hi-kvadrat test
Izračunavanje centralne vrijednosti svih rezultata
Učestalost vrijednosti iznad i ispod medijana
Izračun hi-kvadrat testa

Test sume rangova (Wilcoxon t-test, Mann-Whitney U test) koristi drugačiju proceduru od medijan testa.

Koriste rangove - više informacija od medijan testa
Veća je snaga, snažniji test
Izračun Z vrijednosti iz sume rangova jednog i drugog uzorka (jer su dva nezavisna uzorka)

Wilcoxon test dolazi u dvije inačice, test za nezavisne uzorke (two sample Wilcoxon test) i za zavisne uzorke, kada je riječ o samo jednom uzorku ali u dva ponovljena mjerenja (one sample Wilcoxon test).

Dio neparametrijskih procedura ugrađen je u temeljni R paket stats. Tako, Wilcoxonov test sume rangova poziva se pomoću funkcije wilcox.test.

Zamislimo slijedeći istraživački primjer. U jednom natjecanju kušanja vina sudjelovali su Petar i Ivan koji su dobili sljedeće ocjene (od 1 do 5) od strane kušača.

UlazniPodaci =("
 Vinar Likert
 Petar   4
 Petar   5
 Petar   4
 Petar   3
 Petar   4
 Petar   4
 Petar   4
 Petar   3
 Petar   5
 Petar   5
 Ivan    2
 Ivan    3
 Ivan    2
 Ivan    2
 Ivan    1
 Ivan    2
 Ivan    3
 Ivan    4
 Ivan    2
 Ivan    3
")

PodaciIvanPetar = read.table(textConnection(UlazniPodaci),header=TRUE) #Prebacivanje rezultata u tablicu
# Prebacivanje varijable Likert u factor
#PodaciIvanPetar$Likert = factor(PodaciIvanPetar$Likert, ordered = TRUE)
#PodaciIvanPetar$Vinar = factor(PodaciIvanPetar$Vinar)

Pomoću paketa lattice (Sarkar, 2008) možemo prikazati raspodjelu vrijednosti za Petra i Ivana.

library(lattice)
histogram(~ Likert | Vinar,
          data=PodaciIvanPetar,
          layout=c(1,2)      #  stupci i redovi na odvojenim prikazima
          )

Primjer uporabe neparametrijske statistike može vidjeti koristeći neke od paketa u R-u koji imaju ugrađene funkcije za korištenje neparametrijskih procedura i to: BSDA (R-BSDA?), DescTools (Andri et mult. al., 2022).

Opisna statistika se jednostavno može prikazati pomoću paketa psych (William Revelle, 2023).

library(psych)
describe.by(PodaciIvanPetar$Likert,PodaciIvanPetar$Vinar)


 Descriptive statistics by group 
group: Ivan
   vars  n mean   sd median trimmed  mad min max range skew kurtosis   se
X1    1 10  2.4 0.84      2    2.38 0.74   1   4     3 0.28    -0.84 0.27
------------------------------------------------------------ 
group: Petar
   vars  n mean   sd median trimmed  mad min max range  skew kurtosis   se
X1    1 10  4.1 0.74      4    4.12 0.74   3   5     2 -0.12    -1.35 0.23

Primjena Wilcoxonovog testa za dva nezavisna uzorka u našem primjeru je vrlo jednostavna.

wilcox.test(formula = Likert ~ Vinar, data = PodaciIvanPetar)

Warning in wilcox.test.default(x = DATA[[1L]], y = DATA[[2L]], ...): cannot
compute exact p-value with ties


    Wilcoxon rank sum test with continuity correction

data:  Likert by Vinar
W = 7.5, p-value = 0.001062
alternative hypothesis: true location shift is not equal to 0

Rezultati primjene Wilcoxonovog testa pokazuju kako se procjene značajno razlikuju te kako Petar ima značajno veće vrijednosti tj. dobio je bolje ocjene za svoj proizvod od Ivana.

7.4 Jedan uzorak - dva mjerenja (zavisni test)

Kada imamo model u istraživanju jednog uzorka ali u dva ponovljena mjerenja, tada je test izbora Wilcoxonov test rangova (Wilcoxon Signed-Rank Test). Pozivanje testa je s istovjetnom funkcijom wilcox.test ali je potrebno uključiti opciju paired=TRUE.

Za demonstraciju primjene testa za zavisne uzorke, koristit ćemo primjer iz paketa MASS (Venables & Ripley, 2002) tj. podatke u datoteci immer. Ovaj primjer se odnosi na istraživanje prinosa žitarica u dvije godine (1931 i 1932) tako da prvo trebamo učitati paket MASS a zatim vidjeti zaglavlje datoteke immer.

library(MASS)
head(immer)

  Loc Var    Y1    Y2
1  UF   M  81.0  80.7
2  UF   S 105.4  82.3
3  UF   V 119.7  80.4
4  UF   T 109.7  87.2
5  UF   P  98.3  84.2
6   W   M 146.6 100.4

Prije primjene Wilcoxonovog testa za zavisne uzorke, možemo prikazati slikovno raspodjelu vrijednosti i provjeru značajnosti odstupanja od normaliteta.

#slikovni prikaz raspodjele varijable Y1
hist(x=immer$Y1)

#opisna statistika za varijable Y1 i Y2
describe(immer$Y1)

   vars  n   mean    sd median trimmed   mad  min   max range skew kurtosis
X1    1 30 109.05 28.67 102.95  106.63 30.91 69.1 191.5 122.4 0.79     0.15
     se
X1 5.24

describe(immer$Y2)

   vars  n  mean    sd median trimmed   mad  min   max range skew kurtosis   se
X1    1 30 93.13 24.28  92.95    91.8 24.02 49.9 147.7  97.8 0.37    -0.63 4.43

#Provjera odstupanja od normaliteta
shapiro.test(x=immer$Y1)


    Shapiro-Wilk normality test

data:  immer$Y1
W = 0.92973, p-value = 0.04829

Usporedba da li se značajno razlikuju dvije godine u prinosu žita pomoću wilcox.text() funkcije. Za zavisne uzorke obavezno mora biti uključena funkcija ‘paired=TRUE’.

wilcox.test(immer$Y1, immer$Y2, paired=TRUE)

Warning in wilcox.test.default(immer$Y1, immer$Y2, paired = TRUE): cannot
compute exact p-value with ties


    Wilcoxon signed rank test with continuity correction

data:  immer$Y1 and immer$Y2
V = 368.5, p-value = 0.005318
alternative hypothesis: true location shift is not equal to 0

Rezultati primjene neparametrijskog testa za zavisne uzorke pokazuju kako se prinos žita između dvije godine značajno razlikuje (p<0.01).

7.5 Tri i više nezavisnih uzoraka

Neparametrijska inačica analize varijance (ANOVA) koja se najčešće koristi u praksi je Kruskal-Wallis test. Ovaj test je ekstenzija Wilcoxonovog testa za nezavisne uzorke u situaciji kada imamo više od dva nezavisna uzorka a jednu zavisnu varijablu tj. predmet mjerenja.

Za simulaciju primjene K-W testa možemo upotrijebiti već dostupne podatke u sustavu tj. primjer PlantGrowth gdje su tri eksperimentalne manipulacije i vrijednosti rasta.

Pomoću funkcija provjeravamo zaglavlje datoteke i koliko razina ili skupina ima varijabla koja je u osnovi nezavisna varijabla u ovom modelu.

head(PlantGrowth)

  weight group
1   4.17  ctrl
2   5.58  ctrl
3   5.18  ctrl
4   6.11  ctrl
5   4.50  ctrl
6   4.61  ctrl

levels(PlantGrowth$group)

[1] "ctrl" "trt1" "trt2"

Primjena Kruskal Wallis testa je jednostavna, pomoću funkcije kruskal.test().

kruskal.test(weight ~ group, data = PlantGrowth)


    Kruskal-Wallis rank sum test

data:  weight by group
Kruskal-Wallis chi-squared = 7.9882, df = 2, p-value = 0.01842

Rezultati primjene K-W testa pokazuju kako se grupe značajno razlikuju na razini od 95% (p<0.05) ali ne znamo u kojem smjeru tj. nemamo informaciju između kojih skupina je statistički značajna razlika.

Za post-hoc test koristit ćemo funkciju pairwise.wilcox.test().

pairwise.wilcox.test(PlantGrowth$weight, PlantGrowth$group,
                 p.adjust.method = "BH")

Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
exact p-value with ties


    Pairwise comparisons using Wilcoxon rank sum test with continuity correction 

data:  PlantGrowth$weight and PlantGrowth$group 

     ctrl  trt1 
trt1 0.199 -    
trt2 0.095 0.027

P value adjustment method: BH

Nakon primjene post-hoc testa vidimo da između grupa trt1 i trt2 postoji značajna razlika na razini od 95% (p<0.05).

7.6 Tri i više zavisnih uzoraka

Neparametrijska inačica analize varijance za zavisne uzorke je - Friedmanov test. Koristi se za slučajeve kada zbog malih uzoraka ili nejednakih varijanci nismo u mogućnosti napraviti analizu varijance za ponovljena mjerenja. Za Friedmanov test će nam trebati više paketa, ali nam je većina poznata.

library(tidyverse)
library(ggpubr)
library(rstatix)
library(datarium)

U prvom koraku ćemo pogledati na koji kako izgledaju podatci koje imamo. Za ilustraciju ćemo koristiti podatke o mjerenju samopoštovanja iz paketa datarium.

Važno je primjetiti da kao i kod analize varijance za ponovljena mjerenja postoji “id” sudionika kao poseban stupac. Ta varijabla nam je važna zbog kasnijeg rada s podatcima. U ovom slučaju se samopoštovanje mjerilo tri puta.

data("selfesteem", package = "datarium")
head(selfesteem, 5)

# A tibble: 5 × 4
     id    t1    t2    t3
  <int> <dbl> <dbl> <dbl>
1     1  4.01  5.18  7.11
2     2  2.56  6.91  6.31
3     3  3.24  4.44  9.78
4     4  3.42  4.71  8.35
5     5  2.87  3.91  6.46

Nažalost, iako se unos zavisnih podataka temelji na unosu u različite stupce jedan pored drugog, koji se naziva i “široki” (engl. wide) format, taj način nije kompatibilan s funkcijama za analizu u R-u. Stoga, prije analize moramo presložiti skup podataka iz “širokog” u “dugi” (engl. long) format, gdje sva mjerenja idu jedan ispod drugog. Princip je isti kao i kod analize varijance za ponovljena mjerenja. Koristimo naredbu gather() u kojoj moramo stvoriti i dvije nove varijable: time i score, a koje se odnose na redni broj mjerenja i vrijednost za svako mjerenje. Iz ovog razloga je važno da imamo id sudionika jer id služi kao poveznica između različitih mjerenja za svakog ispitanika, kada su u dugom formatu.

selfesteem <- selfesteem %>%
  pivot_longer(cols=c("t1","t2","t3"), names_to="time", values_to = "score") %>%
convert_as_factor(vars = c("id", "time"))
head(selfesteem)

# A tibble: 6 × 3
  id    time  score
  <fct> <fct> <dbl>
1 1     t1     4.01
2 1     t2     5.18
3 1     t3     7.11
4 2     t1     2.56
5 2     t2     6.91
6 2     t3     6.31

U opisnoj statistici ovih podataka možemo primjetiti kako se vrijednosti povaćavaju nakon svakog mjerenja, ali nismo sigurni je li promjena značajna. Podatke ćemo prikazati i vizualno.

selfesteem %>%
  group_by(time) %>%
  get_summary_stats(score, type = "common")

# A tibble: 3 × 11
  time  variable     n   min   max median   iqr  mean    sd    se    ci
  <fct> <fct>    <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 t1    score       10  2.05  4.00   3.21 0.571  3.14 0.552 0.174 0.395
2 t2    score       10  3.91  6.91   4.60 0.89   4.93 0.863 0.273 0.617
3 t3    score       10  6.31  9.78   7.46 1.74   7.64 1.14  0.361 0.817

ggboxplot(selfesteem, x = "time", y = "score", add = "jitter")

U sljedećem koraku možemo napraviti Friedmanov test, koji nije kompliciran za provesti. Možete primjetiti kako uspoređujemo rezultate ovisno o vremenskoj dimenziji, ali kontrolirajući identifikacijsku oznaku sudionika.

selfesteem_comparison <- selfesteem %>% friedman_test(score ~ time |id)
selfesteem_comparison

# A tibble: 1 × 6
  .y.       n statistic    df        p method       
* <chr> <int>     <dbl> <dbl>    <dbl> <chr>        
1 score    10      18.2     2 0.000112 Friedman test

Rezultat Friedmanovog testa pokazuje kako postoji statistički značajna razlika između skupina, ali, kao i kod analize varijance i K-W testa, ne znamo između kojih skupina. Stoga moramo napraviti post hoc test, a u ovom slučaju možemo koristiti Wilcoxon test za zavisne uzorke. Za razliku od korištenja kod K-W testa, sada ćemo dodati naredbu paired=TRUE.

selfesteem_subgroup_c<-pairwise.wilcox.test(selfesteem$score, selfesteem$time,p.adjust.method = "BH", paired = TRUE)
selfesteem_subgroup_c


    Pairwise comparisons using Wilcoxon signed rank exact test 

data:  selfesteem$score and selfesteem$time 

   t1     t2    
t2 0.0029 -     
t3 0.0029 0.0039

P value adjustment method: BH

U ovom slučaju, postoje statistički značajne razlike između svih mjerenja. Dakle, možemo zaključiti kako je postojala statistički značajna promjena u razini samopoštovanja nakon svakog mjerenja, ili preciznije, razina samopoštovanja se povećavala sa svakim novim mjerenjem.