7. štivo: Testiranje hipoteza kroz primjere - 3. dio (R)
Katarina Kostelić
Ovo je nastavak 7. štiva: Testiranje hipoteza kroz primjere i odnosi se na repliciranje analize koristeći R.
Provedba postupka koristeći R
Prvo ćemo učitati sve podatkovne skupove, pa pogledati kako podaci
izgledaju koristeći head(). JASP koristi GitHub repozitorij
za kod i podatke, tako da podacima možemo pristupiti i direktno putem https://github.com/jasp-stats/jasp-data-library/tree/main,
odabrati odgovarajući skup podataka i direktno učitati u R. Pritom
vodite računa da odabirete .csv format i koristite
raw prikaz (i link).
> heart.rate <- read.csv("https://raw.githubusercontent.com/jasp-stats/jasp-data-library/main/Heart%20Rate/Heart%20Rate.csv")
> head(heart.rate)## Gender Group Heart.Rate
## 1 Female Runners 119
## 2 Female Runners 84
## 3 Female Runners 89
## 4 Female Runners 119
## 5 Female Runners 127
## 6 Female Runners 111> weight.gain <- read.csv("https://raw.githubusercontent.com/jasp-stats/jasp-data-library/refs/heads/main/Weight%20Gain/Weight%20Gain.csv")
> head(weight.gain)## Weight.Before Weight.After Difference
## 1 122.54 135.74 13.20
## 2 120.78 129.36 8.58
## 3 131.12 145.20 14.08
## 4 137.06 145.64 8.58
## 5 163.24 173.80 10.56
## 6 166.32 181.06 14.74> nekretnine <- read.delim("http://sites.williams.edu/rdeveaux/files/2014/09/Saratoga.txt")
> head(nekretnine, 10)## Price Lot.Size Waterfront Age Land.Value New.Construct Central.Air
## 1 132500 0.09 0 42 50000 0 0
## 2 181115 0.92 0 0 22300 0 0
## 3 109000 0.19 0 133 7300 0 0
## 4 155000 0.41 0 13 18700 0 0
## 5 86060 0.11 0 0 15000 1 1
## 6 120000 0.68 0 31 14000 0 0
## 7 153000 0.40 0 33 23300 0 0
## 8 170000 1.21 0 23 14600 0 0
## 9 90000 0.83 0 36 22200 0 0
## 10 122900 1.94 0 4 21200 0 0
## Fuel.Type Heat.Type Sewer.Type Living.Area Pct.College Bedrooms Fireplaces
## 1 3 4 2 906 35 2 1
## 2 2 3 2 1953 51 3 0
## 3 2 3 3 1944 51 4 1
## 4 2 2 2 1944 51 3 1
## 5 2 2 3 840 51 2 0
## 6 2 2 2 1152 22 4 1
## 7 4 3 2 2752 51 4 1
## 8 4 2 2 1662 35 4 1
## 9 3 4 2 1632 51 3 0
## 10 2 2 1 1416 44 3 0
## Bathrooms Rooms
## 1 1.0 5
## 2 2.5 6
## 3 1.0 8
## 4 1.5 5
## 5 1.0 3
## 6 1.0 8
## 7 1.5 8
## 8 1.5 9
## 9 1.5 8
## 10 1.5 6Većina testova spomenutih u štivu dostupna je u osnovnom paketu R-a
stats, dok dodatni paketi mogu pružiti specifične funkcije
ili olakšati interpretaciju rezultata. Evo nekoliko preporučenih
paketa:
stats(ugrađeni paket):- Parametrijski testovi:
- t-test:
t.test() - F-test (test varijanci):
var.test() - ANOVA:
aov(),lm()
- t-test:
- Neparametrijski testovi:
- Wilcoxonov test (za jedan uzorak, nezavisne i uparene uzorke):
wilcox.test() - Kruskal-Wallis test:
kruskal.test() - Friedman test:
friedman.test()
- Wilcoxonov test (za jedan uzorak, nezavisne i uparene uzorke):
- Hi-kvadrat testovi:
- Test sukladnosti, test homogenosti i test neovisnosti:
chisq.test(),fisher.test(),mcnemar.test()
- Test sukladnosti, test homogenosti i test neovisnosti:
- Parametrijski testovi:
BSDA(Basic Statistics and Data Analysis):- Sadrži funkciju
z.test()za provođenje z-testa, što nije izvorno uključeno ustats.
- Sadrži funkciju
rstatix:- Pruža jednostavno sučelje za provođenje i interpretaciju raznih testova (uključujući post-hoc analize) te nudi pregledne izvještaje.
DescTools:- Nudi širok raspon dodatnih funkcija za deskriptivnu statistiku, testove i post-hoc analize.
> install.packages("BSDA")
> install.packages("rstatix")
> install.packages("DescTools")Započinjemo analizom heart.rate skupa podataka. Prvo
koristimo tablicu kontingence kako bi utvrdili strukturu uzorka s
obzirom na spol i vježbanje.
> table(heart.rate$Gender, heart.rate$Group)##
## Control Runners
## Female 200 200
## Male 200 200Sljedeće, dobivamo uvide iz pokazatelja deskriptivne statistike.
> library(psych)
> describe(heart.rate)## vars n mean sd median trimmed mad min max range skew
## Gender* 1 800 1.50 0.5 1.5 1.50 0.74 1 2 1 0.00
## Group* 2 800 1.50 0.5 1.5 1.50 0.74 1 2 1 0.00
## Heart.Rate 3 800 124.49 22.6 124.0 123.77 25.20 69 196 127 0.25
## kurtosis se
## Gender* -2.00 0.02
## Group* -2.00 0.02
## Heart.Rate -0.55 0.80Sljedeće, utvrđujemo pokazatelje desriptivne statistike za varijablu
Heart.Rate za podskupine prema varijabli
Group, što će nam trebati kasnije u analizi.
> describeBy(heart.rate$Heart.Rate, heart.rate$Group)##
## Descriptive statistics by group
## group: Control
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 400 139 18.95 139 139.27 19.27 77 196 119 -0.11 -0.02 0.95
## ------------------------------------------------------------
## group: Runners
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 400 109.98 15.53 109 109.39 16.31 69 164 95 0.37 -0.03 0.78U sljedećem koraku provjeravamo ravnaju li se podaci
Heart.Rate prema normalnoj distribciji.
> shapiro.test(heart.rate$Heart.Rate)##
## Shapiro-Wilk normality test
##
## data: heart.rate$Heart.Rate
## W = 0.98777, p-value = 3.165e-06Ilustrativni prikaz kako bi izgledalo provođenje z-testa, prvo
koristeći paket BSDA, a potom paket
DescTools.
> library(BSDA)## Loading required package: lattice##
## Attaching package: 'BSDA'## The following object is masked from 'package:datasets':
##
## Orange> z.test(heart.rate$Heart.Rate, alternative = "two.sided", mu=120, sigma.x = 22.6)##
## One-sample z-Test
##
## data: heart.rate$Heart.Rate
## z = 5.6193, p-value = 1.917e-08
## alternative hypothesis: true mean is not equal to 120
## 95 percent confidence interval:
## 122.9239 126.0561
## sample estimates:
## mean of x
## 124.49> library(DescTools)##
## Attaching package: 'DescTools'## The following objects are masked from 'package:psych':
##
## AUC, ICC, SD> ZTest(heart.rate$Heart.Rate, alternative = "two.sided", paired = FALSE, mu=120, sd_pop = 22.6)##
## One Sample z-test
##
## data: heart.rate$Heart.Rate
## z = 5.6193, Std. Dev. Population = 22.6, p-value = 1.917e-08
## alternative hypothesis: true mean is not equal to 120
## 95 percent confidence interval:
## 122.9239 126.0561
## sample estimates:
## mean of x
## 124.49Ilustrativni prikaz kako bi izgledalo provođenje klasičnog t-testa
(uz pretpostavku jednakosti varijanci,
var.equal = TRUE).
> t.test(heart.rate$Heart.Rate, alternative = "two.sided", mu=120, paired = FALSE, var.equal = TRUE)##
## One Sample t-test
##
## data: heart.rate$Heart.Rate
## t = 5.6201, df = 799, p-value = 2.637e-08
## alternative hypothesis: true mean is not equal to 120
## 95 percent confidence interval:
## 122.9218 126.0582
## sample estimates:
## mean of x
## 124.49Ilustrativni prikaz kako bi izgledalo provođenje Welchovog t-testa
(uz pretpostavku nejednakosti varijanci,
var.equal = FALSE).
> t.test(heart.rate$Heart.Rate, alternative = "two.sided", mu=120, paired = FALSE, var.equal = FALSE)##
## One Sample t-test
##
## data: heart.rate$Heart.Rate
## t = 5.6201, df = 799, p-value = 2.637e-08
## alternative hypothesis: true mean is not equal to 120
## 95 percent confidence interval:
## 122.9218 126.0582
## sample estimates:
## mean of x
## 124.49Prikaz probedbe Wilcoxonovog testa za jedan uzorak.
> wilcox.test(heart.rate$Heart.Rate, alternative = "two.sided", mu=120, paired = FALSE)##
## Wilcoxon signed rank test with continuity correction
##
## data: heart.rate$Heart.Rate
## V = 186279, p-value = 2.801e-06
## alternative hypothesis: true location is not equal to 120Izračun Cohenovog d, koji mjeri veličinu efekta.
> DescTools::CohenD(heart.rate$Heart.Rate, y=120, pooled = FALSE)## [1] 0.1988191
## attr(,"magnitude")
## [1] "negligible"Također, možemo izračunati i snagu t-testa. Ovdje pretpostavljamo da
je \(\alpha=\beta=0.05\) i to unosimo u
izračun. Dobit ćemo vrijednost delta, koja označava stvarnu
razliku u prosjecima koju ovakav test može ispravno detektirati uz
pretpostavljenu danu razinu pogrešaka tipa I i II.
> power.t.test(n=800, sd=22.6, sig.level = 0.05, power= 0.95, type = "one.sample", alternative = "two.sided")##
## One-sample t test power calculation
##
## n = 800
## delta = 2.883841
## sd = 22.6
## sig.level = 0.05
## power = 0.95
## alternative = two.sidedAlternativno, ako bismo provodili jednosmjerna testiranja, onda bismo na taj način zadali naredbe.
> z.test(heart.rate$Heart.Rate, alternative = "greater", mu=120, sigma.x = 22.6)##
## One-sample z-Test
##
## data: heart.rate$Heart.Rate
## z = 5.6193, p-value = 9.586e-09
## alternative hypothesis: true mean is greater than 120
## 95 percent confidence interval:
## 123.1757 NA
## sample estimates:
## mean of x
## 124.49> ZTest(heart.rate$Heart.Rate, alternative = "greater", paired = FALSE, mu=120, sd_pop = 22.6)##
## One Sample z-test
##
## data: heart.rate$Heart.Rate
## z = 5.6193, Std. Dev. Population = 22.6, p-value = 9.586e-09
## alternative hypothesis: true mean is greater than 120
## 95 percent confidence interval:
## 123.1757 Inf
## sample estimates:
## mean of x
## 124.49> t.test(heart.rate$Heart.Rate, alternative = "greater", mu=120, paired = FALSE, var.equal = TRUE)##
## One Sample t-test
##
## data: heart.rate$Heart.Rate
## t = 5.6201, df = 799, p-value = 1.319e-08
## alternative hypothesis: true mean is greater than 120
## 95 percent confidence interval:
## 123.1744 Inf
## sample estimates:
## mean of x
## 124.49> t.test(heart.rate$Heart.Rate, alternative = "greater", mu=120, paired = FALSE, var.equal = FALSE)##
## One Sample t-test
##
## data: heart.rate$Heart.Rate
## t = 5.6201, df = 799, p-value = 1.319e-08
## alternative hypothesis: true mean is greater than 120
## 95 percent confidence interval:
## 123.1744 Inf
## sample estimates:
## mean of x
## 124.49> wilcox.test(heart.rate$Heart.Rate, alternative = "greater", mu=120, paired = FALSE)##
## Wilcoxon signed rank test with continuity correction
##
## data: heart.rate$Heart.Rate
## V = 186279, p-value = 1.4e-06
## alternative hypothesis: true location is greater than 120> z.test(heart.rate$Heart.Rate, alternative = "less", mu=120, sigma.x = 22.6)##
## One-sample z-Test
##
## data: heart.rate$Heart.Rate
## z = 5.6193, p-value = 1
## alternative hypothesis: true mean is less than 120
## 95 percent confidence interval:
## NA 125.8043
## sample estimates:
## mean of x
## 124.49> ZTest(heart.rate$Heart.Rate, alternative = "less", paired = FALSE, mu=120, sd_pop = 22.6)##
## One Sample z-test
##
## data: heart.rate$Heart.Rate
## z = 5.6193, Std. Dev. Population = 22.6, p-value = 1
## alternative hypothesis: true mean is less than 120
## 95 percent confidence interval:
## -Inf 125.8043
## sample estimates:
## mean of x
## 124.49> t.test(heart.rate$Heart.Rate, alternative = "less", mu=120, paired = FALSE, var.equal = TRUE)##
## One Sample t-test
##
## data: heart.rate$Heart.Rate
## t = 5.6201, df = 799, p-value = 1
## alternative hypothesis: true mean is less than 120
## 95 percent confidence interval:
## -Inf 125.8056
## sample estimates:
## mean of x
## 124.49> t.test(heart.rate$Heart.Rate, alternative = "less", mu=120, paired = FALSE, var.equal = FALSE)##
## One Sample t-test
##
## data: heart.rate$Heart.Rate
## t = 5.6201, df = 799, p-value = 1
## alternative hypothesis: true mean is less than 120
## 95 percent confidence interval:
## -Inf 125.8056
## sample estimates:
## mean of x
## 124.49> wilcox.test(heart.rate$Heart.Rate, alternative = "less", mu=120, paired = FALSE)##
## Wilcoxon signed rank test with continuity correction
##
## data: heart.rate$Heart.Rate
## V = 186279, p-value = 1
## alternative hypothesis: true location is less than 120Za dva uzorka, sad možemo i moramo provjeriti i homogenost varijanci. No, prvo moramo provjeriti jesu li varijable definirane kao odgovarajući tipovi podataka.
> str(heart.rate)## 'data.frame': 800 obs. of 3 variables:
## $ Gender : chr "Female" "Female" "Female" "Female" ...
## $ Group : chr "Runners" "Runners" "Runners" "Runners" ...
## $ Heart.Rate: int 119 84 89 119 127 111 115 109 111 120 ...Za potrebe daljnje analize, Gender i Group
moramo definirati kao faktore. Kad to učinimo na jednostavan način,
koristeći as.factor(), onda će modalitetima varijable biti
pripisane numeričke oznake (bit će prekodirani) po abecednom redu. To
znači da će biti Female = 1 i Male = 2,
odnosno Control = 1 i Runners = 2. Ovo je
korisno zapamtiti (ili zapisati), jer će biti važno za
interpretaciju.
> heart.rate$Group <- as.factor(heart.rate$Group)
> heart.rate$Gender <- as.factor(heart.rate$Gender)
> str(heart.rate)## 'data.frame': 800 obs. of 3 variables:
## $ Gender : Factor w/ 2 levels "Female","Male": 1 1 1 1 1 1 1 1 1 1 ...
## $ Group : Factor w/ 2 levels "Control","Runners": 2 2 2 2 2 2 2 2 2 2 ...
## $ Heart.Rate: int 119 84 89 119 127 111 115 109 111 120 ...Pristupamo provjeri homogenosti varijanci u dvije skupine, koristeći
Leveneov test i Brown-Forsythe test iz paketa DescTools.
Oba testa upućuju na isti zaključak.
> LeveneTest(heart.rate$Heart.Rate, group = heart.rate$Group, center = "mean")## Levene's Test for Homogeneity of Variance (center = "mean")
## Df F value Pr(>F)
## group 1 12.807 0.0003661 ***
## 798
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1> LeveneTest(heart.rate$Heart.Rate, group = heart.rate$Group, center = "median") # Brown-Forsythe-Test## Levene's Test for Homogeneity of Variance (center = "median")
## Df F value Pr(>F)
## group 1 13.238 0.0002919 ***
## 798
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Welchov t-test:
> t.test(heart.rate$Heart.Rate[heart.rate$Group == "Control"], heart.rate$Heart.Rate[heart.rate$Group == "Runners"], alternative = "two.sided", paired = FALSE, var.equal = FALSE)##
## Welch Two Sample t-test
##
## data: heart.rate$Heart.Rate[heart.rate$Group == "Control"] and heart.rate$Heart.Rate[heart.rate$Group == "Runners"]
## t = 23.687, df = 768.42, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 26.61498 31.42502
## sample estimates:
## mean of x mean of y
## 139.00 109.98Mann-Whitney U test (ekstenzija Wilcoxonovog testa za usporedbu dva nezavisna uzorka):
> wilcox.test(heart.rate$Heart.Rate[heart.rate$Group == "Control"], heart.rate$Heart.Rate[heart.rate$Group == "Runners"], alternative = "two.sided", paired = FALSE)##
## Wilcoxon rank sum test with continuity correction
##
## data: heart.rate$Heart.Rate[heart.rate$Group == "Control"] and heart.rate$Heart.Rate[heart.rate$Group == "Runners"]
## W = 140747, p-value < 2.2e-16
## alternative hypothesis: true location shift is not equal to 0Prelazimo na primjenu testova za dva ili više nezavisna uzorka. Ovdje će se koristiti ANOVA i Kruskal-Wallis. Na početku, ispitujemo normalnost rezuduala koristeći Q-Q plot.
> library(ggpubr)## Loading required package: ggplot2##
## Attaching package: 'ggplot2'## The following objects are masked from 'package:psych':
##
## %+%, alpha> model <- lm(Heart.Rate~Group, data = heart.rate)
> ggqqplot(residuals(model))
Osim vizualnog pristupa, možemo koristiti i Shapiro-Wilk test.
> library(rstatix)##
## Attaching package: 'rstatix'## The following object is masked from 'package:stats':
##
## filter> shapiro_test(residuals(model))## # A tibble: 1 × 3
## variable statistic p.value
## <chr> <dbl> <dbl>
## 1 residuals(model) 0.998 0.426Također, provjeravamo razdiobu unutar grupa.
> heart.rate %>%
+ group_by(Group) %>%
+ shapiro_test(Heart.Rate)## # A tibble: 2 × 4
## Group variable statistic p
## <fct> <chr> <dbl> <dbl>
## 1 Control Heart.Rate 0.997 0.772
## 2 Runners Heart.Rate 0.989 0.00409Za velke uzorke, preferiraju se Q-Q plot prikazi naspram testa, jer pri radu s velikim uzorcima Shapiro-Wilk postaje jako osjetljiv i na mala odstupanja od normalnosti.
> ggqqplot(heart.rate, "Heart.Rate", facet.by = "Group")
Postoje tek manja odstupanja od normalnosti (izraženija za trkače), ali moramo voditi računa o tome da nas kod ANOVA-e prvenstveno zanima normalna distribucija reziduala, dok je za podatke važno da dolaze iz normalne distribucije (tolerira se odstupanje od normalnosti). Sljedeći je korak utvrđivanje homogenosti varijanci.
> heart.rate %>% levene_test(Heart.Rate~Group)## # A tibble: 1 × 4
## df1 df2 statistic p
## <int> <int> <dbl> <dbl>
## 1 1 798 13.2 0.000292Temeljem rezultata vidimo da je p-vrijednost < 0,05, što znači da
postoji značajna razlika između varijanci između skupina. Stoga ne
možemo pretpostaviti homogenost varijanci u različitim grupama. U
situaciji kada pretpostavka o homogenosti varijance nije ispunjena,
možemo izračunati Welchov jednosmjerni ANOVA test pomoću funkcije
welch_anova_test() (iz paketa rstatix). Ovaj
test ne zahtijeva pretpostavku o jednakim varijancama.
Za ilustraciju, ovdje se prikazuje kako bi se proveo ANOVA test
(koristeći anova_test() iz paketa
rstatix).
> anova_test(heart.rate, Heart.Rate~Group)## ANOVA Table (type II tests)
##
## Effect DFn DFd F p p<.05 ges
## 1 Group 1 798 561.08 2.38e-94 * 0.413Također, možemo koristiti i funkciju aov() iz paketa
stats. Za ovako dobivene rezultate možemo utvrditi i \(\eta^2\).
> rez <- aov(Heart.Rate~Group, data = heart.rate)
> summary(rez)## Df Sum Sq Mean Sq F value Pr(>F)
## Group 1 168432 168432 561.1 <2e-16 ***
## Residuals 798 239554 300
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1> eta_squared(rez)## Group
## 0.412838> PostHocTest(rez)##
## Posthoc multiple comparisons of means : Tukey HSD
## 95% family-wise confidence level
##
## $Group
## diff lwr.ci upr.ci pval
## Runners-Control -29.02 -31.42487 -26.61513 <2e-16 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Sad prelazimo na rezultate koje možemo koristiti, a to će biti Welchova ANOVA.
> welch_anova_test(heart.rate, heart.rate$Heart.Rate~heart.rate$Group)## # A tibble: 1 × 7
## .y. n statistic DFn DFd p method
## * <chr> <int> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 heart.rate$Heart.Rate 800 561. 1 768. 1.48e-93 Welch ANOVA> ggboxplot(heart.rate, y = "Heart.Rate", facet.by = "Group")
Post-hoc usporedbe:
> heart.rate %>% games_howell_test(Heart.Rate~Group)## # A tibble: 1 × 8
## .y. group1 group2 estimate conf.low conf.high p.adj p.adj.signif
## * <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 Heart.Rate Control Runners -29.0 -31.4 -26.6 0 ****> heart.rate %>% pairwise_t_test(Heart.Rate~Group, pool.sd = FALSE, p.adjust.method = "bonferroni")## # A tibble: 1 × 10
## .y. group1 group2 n1 n2 statistic df p p.adj p.adj.signif
## * <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 Hear… Contr… Runne… 400 400 23.7 768. 1.48e-93 1.48e-93 ****Preostaje još prikazati Kruskall-Wallis test. Taj je test dostupan u
paketu stats i rstatix. Prikazat će se oba,
redom.
> kruskal.test(heart.rate$Heart.Rate~heart.rate$Group)##
## Kruskal-Wallis rank sum test
##
## data: heart.rate$Heart.Rate by heart.rate$Group
## Kruskal-Wallis chi-squared = 345.59, df = 1, p-value < 2.2e-16> kruskal_test(heart.rate, heart.rate$Heart.Rate~heart.rate$Group)## # A tibble: 1 × 6
## .y. n statistic df p method
## * <chr> <int> <dbl> <int> <dbl> <chr>
## 1 heart.rate$Heart.Rate 800 346. 1 3.86e-77 Kruskal-WallisU paketu rstatix na raspolaganju je i funkcija
kruskal_effsize() kojom se utvrđuje veličina efekta (\(\eta^2\)).
> kruskal_effsize(heart.rate, heart.rate$Heart.Rate~heart.rate$Group)## # A tibble: 1 × 5
## .y. n effsize method magnitude
## * <chr> <int> <dbl> <chr> <ord>
## 1 heart.rate$Heart.Rate 800 0.432 eta2[H] largeNa raspolaganju je i Dunnov post-hoc test (prikladan za korištenje kao post-hoc za Kruskal-Wallis test, ako podaci nisu upareni).
> DescTools::DunnTest(heart.rate$Heart.Rate, heart.rate$Group)##
## Dunn's test of multiple comparisons using rank sums : holm
##
## mean.rank.diff pval
## Runners-Control -303.735 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1> boxplot(Heart.Rate ~ Group, data = heart.rate)
Preostaje prikazati postupak za situacije u kojima se u obzir uzimaju
dva faktora. Uz Group, ispitat će se i
Gender.
> ggboxplot(heart.rate, x = "Group", y = "Heart.Rate", color = "Gender")
> model2 <- lm(Heart.Rate ~ Group*Gender, data = heart.rate)
> ggqqplot(residuals(model2))
> shapiro_test(residuals(model2))## # A tibble: 1 × 3
## variable statistic p.value
## <chr> <dbl> <dbl>
## 1 residuals(model2) 0.998 0.338> heart.rate %>%
+ group_by(Group, Gender) %>%
+ shapiro_test(Heart.Rate)## # A tibble: 4 × 5
## Gender Group variable statistic p
## <fct> <fct> <chr> <dbl> <dbl>
## 1 Female Control Heart.Rate 0.995 0.699
## 2 Male Control Heart.Rate 0.994 0.562
## 3 Female Runners Heart.Rate 0.995 0.755
## 4 Male Runners Heart.Rate 0.990 0.196> ggqqplot(heart.rate, "Heart.Rate", ggtheme = theme_bw()) +
+ facet_grid(Group ~ Gender)
> heart.rate %>% levene_test(Heart.Rate ~ Group*Gender)## # A tibble: 1 × 4
## df1 df2 statistic p
## <int> <int> <dbl> <dbl>
## 1 3 796 5.73 0.000697> anova_test(heart.rate, Heart.Rate~Group*Gender)## ANOVA Table (type II tests)
##
## Effect DFn DFd F p p<.05 ges
## 1 Group 1 796 695.647 1.15e-110 * 0.466
## 2 Gender 1 796 185.980 3.29e-38 * 0.189
## 3 Group:Gender 1 796 7.409 7.00e-03 * 0.009Također, možemo koristiti i funkciju aov() iz paketa
stats. Za ovako dobivene rezultate možemo utvrditi i \(\eta^2\).
> rez2 <- aov(Heart.Rate~Group*Gender, data = heart.rate)
> summary(rez2)## Df Sum Sq Mean Sq F value Pr(>F)
## Group 1 168432 168432 695.647 < 2e-16 ***
## Gender 1 45030 45030 185.980 < 2e-16 ***
## Group:Gender 1 1794 1794 7.409 0.00663 **
## Residuals 796 192730 242
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1> eta_squared(rez2)## Group Gender Group:Gender
## 0.412837972 0.110371468 0.004397223> PostHocTest(rez2)##
## Posthoc multiple comparisons of means : Tukey HSD
## 95% family-wise confidence level
##
## $Group
## diff lwr.ci upr.ci pval
## Runners-Control -29.02 -31.17979 -26.86021 <2e-16 ***
##
## $Gender
## diff lwr.ci upr.ci pval
## Male-Female -15.005 -17.16479 -12.84521 <2e-16 ***
##
## $`Group:Gender`
## diff lwr.ci upr.ci pval
## Runners:Female-Control:Female -32.015 -36.02095 -28.009054 <2e-16 ***
## Control:Male-Control:Female -18.000 -22.00595 -13.994054 <2e-16 ***
## Runners:Male-Control:Female -44.025 -48.03095 -40.019054 <2e-16 ***
## Control:Male-Runners:Female 14.015 10.00905 18.020946 <2e-16 ***
## Runners:Male-Runners:Female -12.010 -16.01595 -8.004054 <2e-16 ***
## Runners:Male-Control:Male -26.025 -30.03095 -22.019054 <2e-16 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Sad prelazimo na rezultate koje možemo koristiti, a to će biti Welchova ANOVA.
> heart.rate %>%
+ group_by(Gender) %>%
+ welch_anova_test(Heart.Rate~Group)## # A tibble: 2 × 8
## Gender .y. n statistic DFn DFd p method
## * <fct> <chr> <int> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 Female Heart.Rate 400 394. 1 398. 1.78e-61 Welch ANOVA
## 2 Male Heart.Rate 400 302. 1 364. 1.08e-49 Welch ANOVA> ggboxplot(heart.rate, x= "Gender", y = "Heart.Rate", facet.by = "Group")
Post-hoc usporedbe:
> heart.rate %>%
+ group_by(Gender) %>%
+ games_howell_test(Heart.Rate~Group)## # A tibble: 2 × 9
## Gender .y. group1 group2 estimate conf.low conf.high p.adj p.adj.signif
## * <fct> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 Female Heart.Rate Contr… Runne… -32.0 -35.2 -28.8 0 ****
## 2 Male Heart.Rate Contr… Runne… -26.0 -29.0 -23.1 0 ****> heart.rate %>%
+ group_by(Gender) %>%
+ pairwise_t_test(Heart.Rate~Group, pool.sd = FALSE, p.adjust.method = "bonferroni")## # A tibble: 2 × 11
## Gender .y. group1 group2 n1 n2 statistic df p p.adj
## * <fct> <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <dbl>
## 1 Female Heart.Rate Control Runne… 200 200 19.9 398. 1.78e-61 1.78e-61
## 2 Male Heart.Rate Control Runne… 200 200 17.4 364. 1.08e-49 1.08e-49
## # ℹ 1 more variable: p.adj.signif <chr>Preostaje još prikazati Kruskall-Wallis test. Inačica koja se ovdje
može koristiti dostupna je u paketu rstatix.
> heart.rate %>%
+ group_by(Gender) %>%
+ kruskal_test(Heart.Rate~Group)## # A tibble: 2 × 7
## Gender .y. n statistic df p method
## * <fct> <chr> <int> <dbl> <int> <dbl> <chr>
## 1 Female Heart.Rate 400 210. 1 1.64e-47 Kruskal-Wallis
## 2 Male Heart.Rate 400 181. 1 3.52e-41 Kruskal-Wallis> heart.rate %>%
+ group_by(Group) %>%
+ kruskal_test(Heart.Rate~Gender)## # A tibble: 2 × 7
## Group .y. n statistic df p method
## * <fct> <chr> <int> <dbl> <int> <dbl> <chr>
## 1 Control Heart.Rate 400 93.1 1 4.9 e-22 Kruskal-Wallis
## 2 Runners Heart.Rate 400 60.4 1 7.85e-15 Kruskal-WallisDa bismo dobili detaljnije uvide, možemo kreirati novu varijablu koja će podijeliti skup na četiri podskupa.
> heart.rate$Group_gender <- paste(heart.rate$Group, heart.rate$Gender, sep = "_")
> heart.rate$Group_gender <- factor(heart.rate$Group_gender, levels = c("Control_Female", "Control_Male", "Runners_Female", "Runners_Male"))
>
> str(heart.rate)## 'data.frame': 800 obs. of 4 variables:
## $ Gender : Factor w/ 2 levels "Female","Male": 1 1 1 1 1 1 1 1 1 1 ...
## $ Group : Factor w/ 2 levels "Control","Runners": 2 2 2 2 2 2 2 2 2 2 ...
## $ Heart.Rate : int 119 84 89 119 127 111 115 109 111 120 ...
## $ Group_gender: Factor w/ 4 levels "Control_Female",..: 3 3 3 3 3 3 3 3 3 3 ...> kruskal_test(heart.rate, Heart.Rate~Group_gender)## # A tibble: 1 × 6
## .y. n statistic df p method
## * <chr> <int> <dbl> <int> <dbl> <chr>
## 1 Heart.Rate 800 432. 3 2.79e-93 Kruskal-WallisU paketu rstatix na raspolaganju je i funkcija
kruskal_effsize() kojom se utvrđuje veličina efekta (\(\eta^2\)). Ne može se primijeniti za
interakciju, nego za svaki faktor zasebno.
> kruskal_effsize(heart.rate, Heart.Rate~Group_gender)## # A tibble: 1 × 5
## .y. n effsize method magnitude
## * <chr> <int> <dbl> <chr> <ord>
## 1 Heart.Rate 800 0.539 eta2[H] largeNa raspolaganju je i Dunnov post-hoc test (prikladan za korištenje kao post-hoc za Kruskal-Wallis test, ako podaci nisu upareni).
> DescTools::DunnTest(heart.rate$Heart.Rate, heart.rate$Group_gender)##
## Dunn's test of multiple comparisons using rank sums : holm
##
## mean.rank.diff pval
## Control_Male-Control_Female -167.680 1.2e-12 ***
## Runners_Female-Control_Female -320.615 < 2e-16 ***
## Runners_Male-Control_Female -454.535 < 2e-16 ***
## Runners_Female-Control_Male -152.935 7.2e-11 ***
## Runners_Male-Control_Male -286.855 < 2e-16 ***
## Runners_Male-Runners_Female -133.920 6.8e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1> boxplot(Heart.Rate ~ Group_gender, data = heart.rate, boxwex = 0.5, lex.order = TRUE)
U sljedećem koraku analizira se podatkovni skup
Weight gain.
> str(weight.gain)## 'data.frame': 16 obs. of 3 variables:
## $ Weight.Before: num 123 121 131 137 163 ...
## $ Weight.After : num 136 129 145 146 174 ...
## $ Difference : num 13.2 8.58 14.08 8.58 10.56 ...> describe(weight.gain)## vars n mean sd median trimmed mad min max range
## Weight.Before 1 16 144.64 22.70 138.27 142.51 24.46 117.26 201.74 84.48
## Weight.After 2 16 155.04 21.44 147.84 153.32 22.83 129.36 204.82 75.46
## Difference 3 16 10.41 3.84 11.11 10.56 4.57 3.08 15.62 12.54
## skew kurtosis se
## Weight.Before 0.80 -0.03 5.68
## Weight.After 0.64 -0.58 5.36
## Difference -0.32 -1.34 0.96> shapiro.test(weight.gain$Difference)##
## Shapiro-Wilk normality test
##
## data: weight.gain$Difference
## W = 0.93801, p-value = 0.3254> var.test(weight.gain$Weight.Before, weight.gain$Weight.After, alternative = "two.sided")##
## F test to compare two variances
##
## data: weight.gain$Weight.Before and weight.gain$Weight.After
## F = 1.1217, num df = 15, denom df = 15, p-value = 0.8269
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.3919077 3.2103408
## sample estimates:
## ratio of variances
## 1.121676> t.test(weight.gain$Weight.Before, weight.gain$Weight.After, alternative = "two.sided", paired = TRUE, var.equal = TRUE)##
## Paired t-test
##
## data: weight.gain$Weight.Before and weight.gain$Weight.After
## t = -10.841, df = 15, p-value = 1.71e-08
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
## -12.455282 -8.362218
## sample estimates:
## mean difference
## -10.40875> library(exactRankTests)## Package 'exactRankTests' is no longer under development.
## Please consider using package 'coin' instead.> wilcox.exact(weight.gain$Weight.Before, weight.gain$Weight.After, alternative = "two.sided", paired = TRUE)##
## Exact Wilcoxon signed rank test
##
## data: weight.gain$Weight.Before and weight.gain$Weight.After
## V = 0, p-value = 3.052e-05
## alternative hypothesis: true mu is not equal to 0> data2 <- data.frame(Weight=rep(0,32), time=rep(0,32))
> data2$weight <- c(weight.gain$Weight.Before, weight.gain$Weight.After)
> data2$time <- c(rep(0, length(weight.gain$Weight.Before)), rep(1, length(weight.gain$Weight.After)))
> str(data2)## 'data.frame': 32 obs. of 3 variables:
## $ Weight: num 0 0 0 0 0 0 0 0 0 0 ...
## $ time : num 0 0 0 0 0 0 0 0 0 0 ...
## $ weight: num 123 121 131 137 163 ...> library(effectsize)##
## Attaching package: 'effectsize'## The following objects are masked from 'package:rstatix':
##
## cohens_d, eta_squared## The following object is masked from 'package:psych':
##
## phi> rank_biserial(weight.gain$Weight.Before, weight.gain$Weight.After, paired = TRUE)## r (rank biserial) | 95% CI
## ----------------------------------
## -1.00 | [-1.00, -1.00]Sljedeći po redu je slučaj Nekretnina. Postoji li statistički značajna razlika u cijeni s obzirom na to nalazi li se nekretnina uz obalu?
> library(ggpubr)
> nekretnine$Waterfront<- as.factor(nekretnine$Waterfront)
> model <- lm(nekretnine$Price~nekretnine$Waterfront, data = nekretnine)
> ggqqplot(residuals(model))
Osim vizualnog pristupa, možemo koristiti i Shapiro-Wilk test.
> library(rstatix)
> shapiro_test(residuals(model))## # A tibble: 1 × 3
## variable statistic p.value
## <chr> <dbl> <dbl>
## 1 residuals(model) 0.899 2.60e-32Također, provjeravamo razdiobu unutar grupa.
> nekretnine %>%
+ group_by(Waterfront) %>%
+ shapiro_test(Price)## # A tibble: 2 × 4
## Waterfront variable statistic p
## <fct> <chr> <dbl> <dbl>
## 1 0 Price 0.897 2.20e-32
## 2 1 Price 0.919 1.89e- 1Za velke uzorke, preferiraju se Q-Q plot prikazi naspram testa, jer pri radu s velikim uzorcima Shapiro-Wilk postaje jako osjetljiv i na mala odstupanja od normalnosti. No, ovdje i vizualni prikaz i test upućuju na isti zaključak.
> ggqqplot(nekretnine, "Price", facet.by = "Waterfront")
> nekretnine %>% levene_test(Price~Waterfront)## # A tibble: 1 × 4
## df1 df2 statistic p
## <int> <int> <dbl> <dbl>
## 1 1 1726 4.92 0.0267Temeljem rezultata vidimo da je p-vrijednost < 0,05, što znači da
postoji značajna razlika između varijanci između skupina. Stoga ne
možemo pretpostaviti homogenost varijanci u različitim grupama. U
situaciji kada pretpostavke nisu ispunjene, koristimo Kruskall-Wallis
test. Taj je test dostupan u paketu stats i
rstatix. Prikazat će se oba, redom.
> kruskal.test(nekretnine$Price, nekretnine$Waterfront)##
## Kruskal-Wallis rank sum test
##
## data: nekretnine$Price and nekretnine$Waterfront
## Kruskal-Wallis chi-squared = 20.005, df = 1, p-value = 7.722e-06> kruskal_test(nekretnine, nekretnine$Price ~ nekretnine$Waterfront)## # A tibble: 1 × 6
## .y. n statistic df p method
## * <chr> <int> <dbl> <int> <dbl> <chr>
## 1 nekretnine$Price 1728 20.0 1 0.00000772 Kruskal-WallisU paketu rstatix na raspolaganju je i funkcija
kruskal_effsize() kojom se utvrđuje veličina efekta (\(\eta^2\)).
> kruskal_effsize(nekretnine, nekretnine$Price ~ nekretnine$Waterfront)## # A tibble: 1 × 5
## .y. n effsize method magnitude
## * <chr> <int> <dbl> <chr> <ord>
## 1 nekretnine$Price 1728 0.0110 eta2[H] smallNa raspolaganju je i Dunnov post-hoc test (prikladan za korištenje kao post-hoc za Kruskal-Wallis test, ako podaci nisu upareni).
> DescTools::DunnTest(nekretnine$Price, nekretnine$Waterfront)##
## Dunn's test of multiple comparisons using rank sums : holm
##
## mean.rank.diff pval
## 1-0 578.7573 7.7e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1> boxplot(nekretnine$Price ~ nekretnine$Waterfront, data = nekretnine)
Postoji li statistički značajna razlika u cijeni s obzirom na to koliko kamina ima nekretnina?
> nekretnine$Fireplaces <- as.factor(nekretnine$Fireplaces)
> model <- lm(nekretnine$Price~nekretnine$Fireplaces, data = nekretnine)
> ggqqplot(residuals(model))
> nekretnine %>% levene_test(Price~Fireplaces)## # A tibble: 1 × 4
## df1 df2 statistic p
## <int> <int> <dbl> <dbl>
## 1 4 1723 12.1 0.00000000104> kruskal_test(nekretnine, nekretnine$Price ~ nekretnine$Fireplaces)## # A tibble: 1 × 6
## .y. n statistic df p method
## * <chr> <int> <dbl> <int> <dbl> <chr>
## 1 nekretnine$Price 1728 268. 4 1.07e-56 Kruskal-Wallis> kruskal_effsize(nekretnine, nekretnine$Price ~ nekretnine$Fireplaces)## # A tibble: 1 × 5
## .y. n effsize method magnitude
## * <chr> <int> <dbl> <chr> <ord>
## 1 nekretnine$Price 1728 0.153 eta2[H] large> DescTools::DunnTest(nekretnine$Price, nekretnine$Fireplaces)##
## Dunn's test of multiple comparisons using rank sums : holm
##
## mean.rank.diff pval
## 1-0 365.7473 < 2e-16 ***
## 2-0 665.3022 3.8e-16 ***
## 3-0 928.8736 0.0514 .
## 4-0 1075.6236 0.0163 *
## 2-1 299.5549 0.0011 **
## 3-1 563.1263 0.4434
## 4-1 709.8763 0.2222
## 3-2 263.5714 0.9310
## 4-2 410.3214 0.7676
## 4-3 146.7500 0.9310
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Postoji li statistički značajna razlika u veličini zemljišta s obzirom na to nalazi li se nekretnina uz obalu?
> model <- lm(nekretnine$Lot.Size~nekretnine$Waterfront, data = nekretnine)
> ggqqplot(residuals(model))
> nekretnine %>% levene_test(Lot.Size~Waterfront)## # A tibble: 1 × 4
## df1 df2 statistic p
## <int> <int> <dbl> <dbl>
## 1 1 1726 0.0296 0.863Iako pretpostavka o homogenosti varijanci nije narušena, pretpostavka o normalnosti reziduala jest, pa ćemo koristiti neparametrijske alternative.
> kruskal_test(nekretnine, nekretnine$Lot.Size ~ nekretnine$Waterfront)## # A tibble: 1 × 6
## .y. n statistic df p method
## * <chr> <int> <dbl> <int> <dbl> <chr>
## 1 nekretnine$Lot.Size 1728 0.00645 1 0.936 Kruskal-Wallis> kruskal_effsize(nekretnine, nekretnine$Lot.Size ~ nekretnine$Waterfront)## # A tibble: 1 × 5
## .y. n effsize method magnitude
## * <chr> <int> <dbl> <chr> <ord>
## 1 nekretnine$Lot.Size 1728 -0.000576 eta2[H] small> DescTools::DunnTest(nekretnine$Lot.Size, nekretnine$Waterfront)##
## Dunn's test of multiple comparisons using rank sums : holm
##
## mean.rank.diff pval
## 1-0 10.39019 0.9360
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Postoji li statistički značajna razlika u veličini zemljišta s obzirom na to koliko kamina ima nekretnina?
> model <- lm(nekretnine$Lot.Size~nekretnine$Fireplaces, data = nekretnine)
> ggqqplot(residuals(model))
> nekretnine %>% levene_test(Lot.Size~Fireplaces)## # A tibble: 1 × 4
## df1 df2 statistic p
## <int> <int> <dbl> <dbl>
## 1 4 1723 3.33 0.00998Obje pretpostavke su narušene, pa ćemo koristiti neparametrijske alternative.
> kruskal_test(nekretnine, nekretnine$Lot.Size ~ nekretnine$Fireplaces)## # A tibble: 1 × 6
## .y. n statistic df p method
## * <chr> <int> <dbl> <int> <dbl> <chr>
## 1 nekretnine$Lot.Size 1728 73.1 4 4.95e-15 Kruskal-Wallis> kruskal_effsize(nekretnine, nekretnine$Lot.Size ~ nekretnine$Fireplaces)## # A tibble: 1 × 5
## .y. n effsize method magnitude
## * <chr> <int> <dbl> <chr> <ord>
## 1 nekretnine$Lot.Size 1728 0.0401 eta2[H] small> DescTools::DunnTest(nekretnine$Lot.Size, nekretnine$Fireplaces)##
## Dunn's test of multiple comparisons using rank sums : holm
##
## mean.rank.diff pval
## 1-0 164.8289 1.7e-10 ***
## 2-0 474.9310 1.8e-08 ***
## 3-0 750.1453 0.23586
## 4-0 327.3953 1.00000
## 2-1 310.1021 0.00065 ***
## 3-1 585.3163 0.58445
## 4-1 162.5663 1.00000
## 3-2 275.2143 1.00000
## 4-2 -147.5357 1.00000
## 4-3 -422.7500 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Je li barem 5% promatranih nekretnina novogradnja?
> table(nekretnine$New.Construct)##
## 0 1
## 1647 81> binom_test(x=81, n = 1728, p=0.05, alternative = "less")## # A tibble: 1 × 6
## n estimate conf.low conf.high p p.signif
## * <dbl> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 1728 0.0469 0 0.0561 0.298 nsJe li distribucija tipa goriva (fuel.type) je ista među nekretninama koje imaju i onima koje nemaju središnji sustav za klimatizaciju (central.air)?
> table(nekretnine$Central.Air, nekretnine$Fuel.Type)##
## 2 3 4
## 0 666 248 179
## 1 531 67 37> chisq.test(nekretnine$Central.Air, nekretnine$Fuel.Type)##
## Pearson's Chi-squared test
##
## data: nekretnine$Central.Air and nekretnine$Fuel.Type
## X-squared = 98.079, df = 2, p-value < 2.2e-16Jesu li varijable Fuel.type i Heat.type međusobno neovisne?
> table(nekretnine$Heat.Type, nekretnine$Fuel.Type)##
## 2 3 4
## 2 961 16 144
## 3 230 1 71
## 4 6 298 1> chisq.test(nekretnine$Heat.Type, nekretnine$Fuel.Type)##
## Pearson's Chi-squared test
##
## data: nekretnine$Heat.Type and nekretnine$Fuel.Type
## X-squared = 1594.4, df = 4, p-value < 2.2e-16> CramerV(nekretnine$Heat.Type, nekretnine$Fuel.Type)## [1] 0.6792176> Lambda(nekretnine$Heat.Type, nekretnine$Fuel.Type, direction = "symmetric")## [1] 0.5043937> Lambda(nekretnine$Heat.Type, nekretnine$Fuel.Type, direction = "row")## [1] 0.4645799> Lambda(nekretnine$Heat.Type, nekretnine$Fuel.Type, direction = "column")## [1] 0.5499058Dodatni dobar izvor za učenje inferencijalne statistike uz podršku R-a je udžbenik Summary and Analysis of Extension Program Evaluation in R (Mangiafico, S.S. 2016.).
Korišteni izvori i literatura
Arnholt, A. T., & Evans, B. (2017). Package ‘BSDA’. https://cran.r-project.org/web/packages/BSDA/index.html
Ben-Shachar, M. S., Lüdecke, D., & Makowski, D. (2020). effectsize: Estimation of effect size indices and standardized parameters. Journal of open source software, 5(56), 2815. https://dominiquemakowski.github.io/publication/benshachar2020effectsize/benshachar2020effectsize.pdf, https://CRAN.R-project.org/package=effectsize
Conover, W. J. (1999). Practical nonparametric statistics (Vol. 350). John Wiley & sons.
De Veaux, D. (2015). How much is a Fireplace Worth? Stats 101 Public Library.https://community.amstat.org/stats101/resources/viewdocument?DocumentKey=e4f8d3f1-41a3-4f01-9f8b-f8fbe1562c15&tab=librarydocuments&CommunityKey=5ad27b39-58d0-49e9-9f6f-0c39c82a0401
Excel, M. S. (2007). Microsoft Excel. Denver Co., USA.
Hohenwarter, M., & Hohenwarter, M. (2002). GeoGebra. Available on-line at http://www.geogebra.org/cms/en.
Hollander, M., Wolfe, D. A. & Chicken, E. (2013). Nonparametric statistical methods. John Wiley & Sons Inc.
Holmes, A., Illowsky, B., & Dean, S. (2017). Introductory Business Statistics 2e. OpenStax. https://openstax.org/books/introductory-business-statistics/pages/preface
Hornik, K. (2012). The comprehensive R archive network. Wiley interdisciplinary reviews: Computational statistics, 4(4), 394-398.
Horton, N. J., Baumer, B. S., & Wickham, H. (2015). Setting the stage for data science: integration of data management skills in introductory and second courses in statistics (nycflights13). https://nhorton.people.amherst.edu/precursors/nycflights13.pdf
Hothorn, T., Hornik, K., & Hothorn, M. T. (2022). Package ‘exactRankTests’. https://CRAN.R-project.org/package=exactRankTests
Horvat, J., & Mijoč, J. (2018). Osnove statistike, treće dopunjeno izdanje. Zagreb: Ljevak.
Illowsky, B., & Dean, S. (2018). Introductory statistics. https://openstax.org/books/introductory-statistics-2e/pages/preface
JASP Team (2024). JASP (Version 0.19.3)[Computer software].
Kassambara, A. (2019). rstatix: Pipe-friendly framework for basic statistical tests. CRAN: Contributed Packages. https://CRAN.R-project.org/package=rstatix
Kassambara, A. (2023) ggpubr: ‘ggplot2’ Based Publication Ready Plots. CRAN https://CRAN.R-project.org/package=ggpubr
Kostelić, K. & Etinger, D. (2024). Uvod u R i RStudio. Sveučilište Jurja Dobrile u Puli. https://bookdown.org/kakoste/Uvod_u_R_i_RStudio/
Moore, D. S., McCabe, G. P., and Craig, B. A. (2012). Introduction to the Practice of Statistics (7th ed.). New York: Freeman
Okoye, K., & Hosseini, S. (2024). Analysis of variance (ANOVA) in R: one-way and two-way ANOVA. In R Programming: Statistical Data Analysis in Research (pp. 187-209). Singapore: Springer Nature Singapore. https://link.springer.com/chapter/10.1007/978-981-97-3385-9_9
Revelle, W., & Revelle, M. W. (2015). Package ‘psych’. The comprehensive R archive network, 337(338), 161-165. https://cran.rstudio.org/web/packages/psych/psych.pdf
Signorell, A. (2025) DescTools: Tools for Descriptive Statistics. CRAN. https://CRAN.R-project.org/package=DescTools
Šošic, I. (2004). Primijenjena statistika. Skolska knjiga, Zagreb.
Wickham, H., Francois, R., Henry, L., & Müller, K. (2014). dplyr. A Grammar of Data Manipulation 2020 [Last accessed on 2020 Aug 12] Available from, Rproject.
Yarberry, W., & Yarberry, W. (2021). Dplyr. CRAN recipes: DPLYR, stringr, lubridate, and regex in R, 1-58.