1. domača naloga
Nina Cerar
2022-12-01
data(iris)
podatki <- force(iris)Raziskovalno vprašanje
Za vrste rož setosa, versicolor in virginica so izmerili dolžino in širino cvetnega in čašnega lista. Izbrali so 50 rož vsake vrste in so jim izmerili liste. Ali lahko trdimo, da se povprečna dolžina čašnih listov razlikuje od povrečne dolžine cvetnih listov, za vrsto virginica?
head(podatki)## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1 5.1 3.5 1.4 0.2 setosa
## 2 4.9 3.0 1.4 0.2 setosa
## 3 4.7 3.2 1.3 0.2 setosa
## 4 4.6 3.1 1.5 0.2 setosa
## 5 5.0 3.6 1.4 0.2 setosa
## 6 5.4 3.9 1.7 0.4 setosaOpis spremenljivk:
- Sepal.Length: dolžina časnega lista v cm
- Sepal.Width: širina časnega lista v cm
- Petal.Length: dolžina cvetnega lista v cm
- Petal.Width: širina cvetnega lista v cm
- Species: vrsta rože (setosa, versicolor, virginica)
podatki$SpeciesF <- factor(podatki$Species,
levels = c("setosa", "versicolor", "virginica"),
labels = c(1, 2,3))
head(podatki)## Sepal.Length Sepal.Width Petal.Length Petal.Width Species SpeciesF
## 1 5.1 3.5 1.4 0.2 setosa 1
## 2 4.9 3.0 1.4 0.2 setosa 1
## 3 4.7 3.2 1.3 0.2 setosa 1
## 4 4.6 3.1 1.5 0.2 setosa 1
## 5 5.0 3.6 1.4 0.2 setosa 1
## 6 5.4 3.9 1.7 0.4 setosa 1Opisna statistika uporabljenih spremenljivk
podatki1 <- podatki[grepl("virginica", podatki$Species),]
podatki1$diferenca <- podatki1$Sepal.Length - podatki1$Petal.Length
head(podatki1[,c(1,3,7)])## Sepal.Length Petal.Length diferenca
## 101 6.3 6.0 0.3
## 102 5.8 5.1 0.7
## 103 7.1 5.9 1.2
## 104 6.3 5.6 0.7
## 105 6.5 5.8 0.7
## 106 7.6 6.6 1.0library(pastecs)
round(stat.desc(podatki1[,c(1,3,7)]),2)## Sepal.Length Petal.Length diferenca
## nbr.val 50.00 50.00 50.00
## nbr.null 0.00 0.00 0.00
## nbr.na 0.00 0.00 0.00
## min 4.90 4.50 0.30
## max 7.90 6.90 1.80
## range 3.00 2.40 1.50
## sum 329.40 277.60 51.80
## median 6.50 5.55 1.00
## mean 6.59 5.55 1.04
## SE.mean 0.09 0.08 0.05
## CI.mean.0.95 0.18 0.16 0.09
## var 0.40 0.30 0.10
## std.dev 0.64 0.55 0.32
## coef.var 0.10 0.10 0.31Iz zgornjih podatkov lahko vidimo, da je polovica čašnih listov dolgih do vključno 6.5 centimetrov, ostali so daljši. Pri cvetnih listih pa je polovica dolgih do 5.55 centimetrov, ostali so pa daljši. Variacijski razmik za dolžino cvetnih listov znaša 2.4 centimetrov. Največje odstopanje od povprečja se opazi pri čašnih listih.
Preizkus domneve o razliki dveh aritmetičnih sredin za odvisne vzorce
- HO: povprečna dolžina čašnih in cvetnih listov je enaka.
- H1: povprečna dolžina čašnih in cvetnih listov je različna.
library(ggplot2)## Warning: package 'ggplot2' was built under R version 4.2.2ggplot(podatki1, aes(x = diferenca)) +
geom_histogram(binwidth = 0.2 , colour = "pink") +
ylab("frekvenca") +
xlab("diferenca")
Opažam, da je porazdelitev dokaj normalna. Preverimo še s shapiro testom za normalnost.
shapiro.test(podatki1$diferenca)##
## Shapiro-Wilk normality test
##
## data: podatki1$diferenca
## W = 0.98199, p-value = 0.6384Domneva: - Ho: porazdelitev diferenc je normalna. - H1: porazdelitev diferenc ni normalna.
Potrdimo ničelno domnevo pri p vrednost = 0.639. Porazdelitev diferenc je normalna.
Nadaljujemo s t testom: - HO: povprečna dolžina čašnih in cvetnih listov je enaka. - H1: povprečna dolžina čašnih in cvetnih listov je različna.
t.test(podatki1$Sepal.Length, podatki1$Petal.Length,
paired = TRUE,
alternative = "two.sided")##
## Paired t-test
##
## data: podatki1$Sepal.Length and podatki1$Petal.Length
## t = 22.898, df = 49, p-value < 2.2e-16
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
## 0.9450788 1.1269212
## sample estimates:
## mean difference
## 1.036library(effectsize)
effectsize::cohens_d(podatki1$diferenca)## Cohen's d | 95% CI
## ------------------------
## 3.24 | [2.54, 3.93]interpret_cohens_d(3.24, rules="sawilowsky2009")## [1] "huge"
## (Rules: sawilowsky2009)Sklep: Na podlagi vzorčnih podatkov ugotavljamo, da se povprečna dolžina cvetnih in čašnih listov razlikuje (p vrednost < 0.001, zavrnemo ničelno domnevo). Velikost učinka je ogromna (3.24). Ocena povprečne razlike je 1.036 cm, kar pomeni, da je bila povprečna dolžina čašnih listov daljša.
Enofaktorska analiza variance za odvisne vzorce
Za vrste rož setosa, versicolor in virginica so izmerili dolžino in širino cvetnega in čašnega lista. Izbrali so 50 rož vsake vrste in so jim izmerili liste. Ali obstajajo razlike v meritvah za cvet vrste virginica?
podatki <- tibble::rowid_to_column(podatki, "ID")
head(podatki)## ID Sepal.Length Sepal.Width Petal.Length Petal.Width Species SpeciesF
## 1 1 5.1 3.5 1.4 0.2 setosa 1
## 2 2 4.9 3.0 1.4 0.2 setosa 1
## 3 3 4.7 3.2 1.3 0.2 setosa 1
## 4 4 4.6 3.1 1.5 0.2 setosa 1
## 5 5 5.0 3.6 1.4 0.2 setosa 1
## 6 6 5.4 3.9 1.7 0.4 setosa 1Opis spremenljivk:
- ID: identifikator posamezne rože
- Sepal.Length: dolžina časnega lista v cm
- Sepal.Width: širina časnega lista v cm
- Petal.Length: dolžina cvetnega lista v cm
- Petal.Width: širina cvetnega lista v cm
- Species: vrsta rože (setosa, versicolor, virginica)
- SpeciesF: vrsta rože (1 - setosa, 2 - versicolor, 3 - virginica)
library(rstatix)## Warning: package 'rstatix' was built under R version 4.2.2##
## Attaching package: 'rstatix'## The following objects are masked from 'package:effectsize':
##
## cohens_d, eta_squared## The following object is masked from 'package:stats':
##
## filterlibrary(tidyverse)## Warning: package 'tidyverse' was built under R version 4.2.2## ── Attaching packages
## ───────────────────────────────────────
## tidyverse 1.3.2 ──## ✔ tibble 3.1.8 ✔ dplyr 1.0.10
## ✔ tidyr 1.2.1 ✔ stringr 1.4.1
## ✔ readr 2.1.3 ✔ forcats 0.5.2
## ✔ purrr 0.3.5
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ tidyr::extract() masks pastecs::extract()
## ✖ dplyr::filter() masks rstatix::filter(), stats::filter()
## ✖ dplyr::first() masks pastecs::first()
## ✖ dplyr::lag() masks stats::lag()
## ✖ dplyr::last() masks pastecs::last()podatki <- podatki[grepl("virginica", podatki$Species),]
podatki_dolg_format <- podatki[,c(-6,-7)] %>%
gather(key = "cvet", value ="meritve", Sepal.Length, Sepal.Width, Petal.Length, Petal.Width) %>% convert_as_factor(cvet)
tail(podatki_dolg_format, 10)## ID cvet meritve
## 191 141 Petal.Width 2.4
## 192 142 Petal.Width 2.3
## 193 143 Petal.Width 1.9
## 194 144 Petal.Width 2.3
## 195 145 Petal.Width 2.5
## 196 146 Petal.Width 2.3
## 197 147 Petal.Width 1.9
## 198 148 Petal.Width 2.0
## 199 149 Petal.Width 2.3
## 200 150 Petal.Width 1.8Naredimo podatke z dolgih formatom.
library(ggpubr)## Warning: package 'ggpubr' was built under R version 4.2.2ggboxplot(podatki_dolg_format,
x = "cvet",
y = "meritve",
add = "jitter")
Glede na grafikon v skupinah opažamo prisotnost osamelcev. Porazdelitev znotraj skupin zgleda normalno.
Poiščemo osamelce in jih odstranimo. Postopek ponovimo dokler se ne znebimo vseh osamelcev.
podatki_dolg_format %>%
group_by(cvet) %>%
identify_outliers(meritve)## # A tibble: 4 × 5
## cvet ID meritve is.outlier is.extreme
## <fct> <int> <dbl> <lgl> <lgl>
## 1 Sepal.Length 107 4.9 TRUE FALSE
## 2 Sepal.Width 118 3.8 TRUE FALSE
## 3 Sepal.Width 120 2.2 TRUE FALSE
## 4 Sepal.Width 132 3.8 TRUE FALSEpodatki_dolg_format <- podatki_dolg_format %>% filter(!ID == 107, !ID == 118, !ID == 120, !ID == 132)podatki_dolg_format %>%
group_by(cvet) %>%
identify_outliers(meritve)## # A tibble: 2 × 5
## cvet ID meritve is.outlier is.extreme
## <fct> <int> <dbl> <lgl> <lgl>
## 1 Petal.Length 119 6.9 TRUE FALSE
## 2 Sepal.Width 110 3.6 TRUE FALSEpodatki_dolg_format <- podatki_dolg_format %>% filter(!ID == 119, !ID == 110)podatki_dolg_format %>%
group_by(cvet) %>%
identify_outliers(meritve)## # A tibble: 2 × 5
## cvet ID meritve is.outlier is.extreme
## <fct> <int> <dbl> <lgl> <lgl>
## 1 Sepal.Length 123 7.7 TRUE FALSE
## 2 Sepal.Length 136 7.7 TRUE FALSEpodatki_dolg_format <- podatki_dolg_format %>% filter(!ID == 123, !ID == 136)podatki_dolg_format %>%
group_by(cvet) %>%
identify_outliers(meritve)## # A tibble: 1 × 5
## cvet ID meritve is.outlier is.extreme
## <fct> <int> <dbl> <lgl> <lgl>
## 1 Sepal.Length 106 7.6 TRUE FALSEpodatki_dolg_format <- podatki_dolg_format %>% filter(!ID == 106)podatki_dolg_format %>%
group_by(cvet) %>%
identify_outliers(meritve)## [1] cvet ID meritve is.outlier is.extreme
## <0 rows> (or 0-length row.names)S shapiro testom preverimo ali so porazdelitve znotraj posamezne skupine normalne. - Ho: porazdelitev je normalna. - H1: porazdelitev ni normalna.
library(rstatix)
podatki_dolg_format %>% group_by(cvet) %>% shapiro_test(meritve)## # A tibble: 4 × 4
## cvet variable statistic p
## <fct> <chr> <dbl> <dbl>
## 1 Petal.Length meritve 0.962 0.190
## 2 Petal.Width meritve 0.950 0.0689
## 3 Sepal.Length meritve 0.977 0.567
## 4 Sepal.Width meritve 0.962 0.179Glede na dane ugotovitve (visoke p vrednosti: 0.190, 0.069, 0.568 in 0.179) lahko trdimo, da so porazdelitve znotraj posamezne skupine normalne.
library(rstatix)
podatki_dolg_format %>% group_by(cvet) %>% get_summary_stats(meritve, type="common")## # A tibble: 4 × 11
## cvet variable n min max median iqr mean sd se ci
## <fct> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Petal.Length meritve 41 4.8 6.3 5.5 0.6 5.43 0.388 0.061 0.122
## 2 Petal.Width meritve 41 1.4 2.5 2 0.5 2.02 0.269 0.042 0.085
## 3 Sepal.Length meritve 41 5.6 7.4 6.4 0.5 6.46 0.445 0.069 0.14
## 4 Sepal.Width meritve 41 2.5 3.4 3 0.3 2.96 0.242 0.038 0.076Domneve: - Ho: aritmetične sredine za meritve cvetov čašne in cvetne dolžine ter širine so enaka. - H1: vsaj ena aritmetična sredina za meritve cvetov čašne in cvetne dolžine ter širine se razlikuje.
library(rstatix)
ANOVA <- anova_test(dv = meritve,
wid = ID,
within = cvet,
data = podatki_dolg_format)
get_anova_table(ANOVA, correction ="auto")## ANOVA Table (type III tests)
##
## Effect DFn DFd F p p<.05 ges
## 1 cvet 1.92 76.84 2275.956 1.61e-68 * 0.965library(effectsize)
interpret_eta_squared(0.965, rules = "cohen1992")## [1] "large"
## (Rules: cohen1992)Ker se ob ponavljanju preverjanja domnev na istem vzorcu podatkov napihuje verjetnost za napako I. vrste, je kriterij za zavrnitev ničelne domneve običajno strožji. Uporabimo bonferrorijev popravek.
library(rstatix)
popravek <- podatki_dolg_format %>%
pairwise_t_test(meritve ~ cvet, paired = TRUE, p.adjust.method = "bonferroni")
popravek## # A tibble: 6 × 10
## .y. group1 group2 n1 n2 stati…¹ df p p.adj p.adj…²
## * <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 meritve Petal.Leng… Petal… 41 41 50.4 40 7.85e-38 4.71e-37 ****
## 2 meritve Petal.Leng… Sepal… 41 41 -20.9 40 3.84e-23 2.3 e-22 ****
## 3 meritve Petal.Leng… Sepal… 41 41 40.8 40 3.32e-34 1.99e-33 ****
## 4 meritve Petal.Width Sepal… 41 41 -56.9 40 6.87e-40 4.12e-39 ****
## 5 meritve Petal.Width Sepal… 41 41 -24.9 40 5.97e-26 3.58e-25 ****
## 6 meritve Sepal.Leng… Sepal… 41 41 52.2 40 1.97e-38 1.18e-37 ****
## # … with abbreviated variable names ¹statistic, ²p.adj.signifSklep: Ugotavljamo, da se vsaj ena aritmetična sredina razlikuje od ostalih (p<0.001) in da obstaja vsaj ena meritev, ki je različna v skupini meritev za cvet vrste virginica. Ne vemo, katera to je. Velikost učinka je velika (0.965).