Testovanie
Simona Mikovčíková
5.12.2020
library(readxl)
library(randtests)
library(ggplot2)
library(cowplot)
library(nortest)
library(moments)
library(RColorBrewer)
library(dunn.test)
library(outliers)
library(EnvStats)##
## Attaching package: 'EnvStats'## The following objects are masked from 'package:moments':
##
## kurtosis, skewness## The following objects are masked from 'package:stats':
##
## predict, predict.lm## The following object is masked from 'package:base':
##
## print.defaultlibrary(car)## Loading required package: carData##
## Attaching package: 'car'## The following object is masked from 'package:EnvStats':
##
## qqPlotlibrary(BSDA)## Loading required package: lattice##
## Attaching package: 'BSDA'## The following objects are masked from 'package:carData':
##
## Vocab, Wool## The following object is masked from 'package:datasets':
##
## Orange1 COVID(I)
Máme reálne dáta pre Slovensko, kde by ma zaujímalo, či je stredná hodnota počtu vykonaných testov za deň za sledované obdobie rovná 2000 testom, testujte na hladine významnosti 0.01.
data <- read_excel("C:/Users/Simi/Desktop/R_ko(3.semester)/3_zadanie/data_testovanie_zadanie.xlsx", sheet = 1)
head(data)## # A tibble: 6 x 8
## Datum `Pocet potvrden~ `Pocet vyliecen~ `Pocet aktivnyc~
## <dttm> <dbl> <dbl> <dbl>
## 1 2020-03-06 00:00:00 1 0 1
## 2 2020-03-07 00:00:00 3 0 3
## 3 2020-03-08 00:00:00 5 0 5
## 4 2020-03-09 00:00:00 7 0 7
## 5 2020-03-10 00:00:00 7 0 7
## 6 2020-03-11 00:00:00 10 0 10
## # ... with 4 more variables: `Dennych testov` <dbl>, `Dennych
## # prirastkov` <dbl>, `Pocet umrti` <dbl>, `Denny prirastok umrti` <dbl>denne<-data$`Dennych testov`
head(denne)## [1] 378 66 76 69 111 97denne## [1] 378 66 76 69 111 97 95 58 197 289 159 318 225 301 368
## [16] 440 242 451 476 335 913 747 720 401 688 877 1191 1454 1889 1524
## [31] 1036 1448 2042 1690 2301 2174 1580 1324 1302 1439 2967 3351 3144 3323 2458
## [46] 2694 3468 4525 3840 4828 4839 3171 1767 5472 4584 5150 3698 1450 1584 2060
## [61] 4742 5161 4694 3910 1488 786 2063 4326 4876 3992 4084 2476 971 2041 3371
## [76] 2933 2751 2236 1649 645 1464 2839 2352 1848 3433 1606 274 6418 2336 2135
## [91] 1832 2639 1180 160 851 1545 1500 1262 1511 479 47 847 1163 787 806
## [106] 1278 301 41 661 1257 936 1515 1611 931 62 784 2063 1708 1801 2216
## [121] 808 50 873 2225 2284 2172 2879 960 279 1163 2205 2336 1862 2161 410
## [136] 24 3333 2571 2251 2049 2275 767 216 1548 2296 1851 2176 2884 585 766
## [151] 1320 2538 2667 2473 3099 1068 564 1454 3131 2741 2738 3235 2013 481 1583
## [166] 3684 3435 3245 3833 1723 929 2103 4090 3636 4360 4453 1951 588 2763 2428
## [181] 3519 4772 5947 2462 922 2891 5309 5021 4266 6191 3080 1425 4323 3235 4027
## [196] 5542 5750 3443 1952 2664 6231 5213 5540 6483 5655 1954 2561 7350 9170 7824
## [211] 9492 8565 3382 3750 8679 9518alfa<-0.01
df1<-data.frame(denne)Testujeme náhodnosť hodôt. Použijem Wald-Wolfowitz runs test. \[H_0:výber je náhodný\] \[H_1: ¬H_0\]
p<-runs.test(denne,plot=T)
p$p.value<alfa## [1] TRUEturning.point.test(denne)##
## Turning Point Test
##
## data: denne
## statistic = -4.1594, n = 216, p-value = 3.19e-05
## alternative hypothesis: non randomnessNulovú hypotézu o náhodnosti zamietame na hladine významnosti 0.01, čiže dáta nie sú náhodné.Môže to byť spôbené tým,že priebeh pandémie ma exponenciálny charakter a pretože počet nakazených stúpa covid sa šíri čo raz rýchlejšie a preto sa denne vykona viac testov. Testujeme normalitu premennej Denné testy \[H_0: N(\mu,\sigma)\] \[H_1:¬ H_0\]
hist(denne, # histogram
col="light blue",
border="black",
prob = TRUE,
main = "Histogram Dennych testov")
lines(density(denne),
lwd = 2,
col = "dark blue")
Z grafu vidíme,že stredná hodnota je približne 2000 čo zodpovedá zadaní z úlohy. Testujeme normalitu Shapiro-Wilkovým aj Lillieforsovým testom
test1<-lillie.test(denne)
test1$p.value<alfa## [1] TRUEtest2<-shapiro.test(denne)
test2$p.value<alfa## [1] TRUEKeďže p-hodnota je menšia ako hladina významnosti 0.01, \(H_0\) zamietame, Prírastok denných testov nie je daná normálnym rozdelením.
qqnorm(denne,pch=1)
qqline(denne,col="blue",lwd=2)
Z grafu vidíme, že body sa vychylujú od priamky z čoho môžme usúdiť v tomto pripade sa nejedná o normálne rozdelenie.
qqPlot(denne)
## [1] 216 211ggplot(df1, aes(sample=denne)) +
stat_qq() + stat_qq_line(colour= "#FF6666")+
labs(x= "Teoretický", y = "Nameraný", title = "QQ-plot")+
theme(plot.title = element_text(hjust=0.5))
Aj z kvantil-kvantilového grafu môžeme vidieť, že dáta zrejme nebudú normálne rozdelené. Aj keď väčšina bodov sa nachádza na priamke.
Testujeme prítomnosť extremálnych hodnôt pomocou metódy založenej na medzikvantilovom rozpätí.
k1<-1.5 # pre vybočujúce hodnoty
k2<-3 # pre extremálne hodnoty
Q1<-quantile(denne, probs=0.25) # dolný kvartil
Q3<-quantile(denne, probs=0.75) # horný kvartil
IQR<-IQR(denne)
denne[denne<(Q1-k1*IQR)]## numeric(0)denne[denne>(Q3+k1*IQR)]## [1] 7350 9170 7824 9492 8565 8679 9518data$Datum[data$`Dennych testov`==7350]## [1] "2020-09-29 UTC"data$Datum[data$`Dennych testov`==9170]## [1] "2020-09-30 UTC"data$Datum[data$`Dennych testov`==7824]## [1] "2020-10-01 UTC"data$Datum[data$`Dennych testov`==9492]## [1] "2020-10-02 UTC"data$Datum[data$`Dennych testov`==8565]## [1] "2020-10-03 UTC"data$Datum[data$`Dennych testov`==8679]## [1] "2020-10-06 UTC"data$Datum[data$`Dennych testov`==9518]## [1] "2020-10-07 UTC"denne[denne>(Q3+k2*IQR)]## numeric(0)denne[denne<(Q1-k2*IQR)]## numeric(0)Nízke vybočujúce hodnoty v dátach nie sú. Vysoké vybočujúce hodnoty, ktoré sú extremálymi je 7. Všetky hodnoty sa nachádzajú v intervale od 29.9.2000 po 7.10.2020 s vynímkou 5.10.2020 tam nastal pokles.
boxplot(denne, col ="red")
hodnoty sú vyšikmené, čo vidíme aj z boxplotu, a preto použijeme znamienkový test
skewness(denne)## [1] 1.308421SIGN.test(denne,md=2000,conf.level = 1-alfa)##
## One-sample Sign-Test
##
## data: denne
## s = 111, p-value = 0.7338
## alternative hypothesis: true median is not equal to 2000
## 99 percent confidence interval:
## 1596.551 2316.033
## sample estimates:
## median of x
## 2045.5
##
## Achieved and Interpolated Confidence Intervals:
##
## Conf.Level L.E.pt U.E.pt
## Lower Achieved CI 0.9884 1606.000 2301.000
## Interpolated CI 0.9900 1596.551 2316.033
## Upper Achieved CI 0.9922 1584.000 2336.000Medián je štatisticky významne odlišný od strednej hodnoty 2000. H_o zamietame na hladine významnosti 0.01.Stredná hodnota počtu vykonaných testov za deň v priebehu sledovaného obdobia nie je rovná 2000. Pri parametrickej verzii by som použila test pri neznámom rozptyle (smean-mu0)/(s/sqrt(n))
2 COVID(II)
V tejo časti by ma zaujímalo, či sa počet denných prírastkov za obdobie apríl až máj štatisticky významne líši od počtu denných prírastkov za obdobie jún až júl. Testujte na hladine významnosti 0.05.
jar<-data$`Dennych prirastkov`[27:87]
leto<-data$`Dennych prirastkov`[88:148]
alfa <-0.05
df_jar<-data.frame(jar)
df_leto<-data.frame(leto)Najprvn testujeme náhodnosť dát. \[H_0:výber je náhodný\] \[H_1: ¬H_0\]
runs.test(jar)##
## Runs Test
##
## data: jar
## statistic = -6.2782, runs = 6, n1 = 27, n2 = 30, n = 57, p-value =
## 3.424e-10
## alternative hypothesis: nonrandomnessruns.test(leto)##
## Runs Test
##
## data: leto
## statistic = -4.4271, runs = 14, n1 = 30, n2 = 30, n = 60, p-value =
## 9.553e-06
## alternative hypothesis: nonrandomnessturning.point.test(jar)##
## Turning Point Test
##
## data: jar
## statistic = -0.32219, n = 56, p-value = 0.7473
## alternative hypothesis: non randomnessturning.point.test(leto)##
## Turning Point Test
##
## data: leto
## statistic = -1.0547, n = 58, p-value = 0.2916
## alternative hypothesis: non randomnessNemôžme zamietnuť nulovú hypotézu o náhodnoti.Dáta sú náhodné. Testujeme normalitu
g1<-ggplot(df_jar, aes(jar)) +
geom_histogram(aes(y= ..density..), bins=7, fill="#4472c4")+
labs(x="Namerané hodnoty", y="Pravdepodobnosť", title = "Histogram")+
geom_density(alpha=.4, fill="#FF6666")+
theme(plot.title = element_text(hjust=0.5))
g2<-ggplot(df_leto, aes(leto)) +
geom_histogram(aes(y= ..density..), bins=7, fill="#4472c4")+
labs(x="Namerané hodnoty", y="Pravdepodobnosť", title = "Histogram")+
geom_density(alpha=.4, fill="#FF6666")+
theme(plot.title = element_text(hjust=0.5))
plot_grid(g1,g2,labels = "AUTO")
Z histogramov nie je symetrický a zdá že dáta nebuduu normálne rozdelené.
Otestujeme normalitu dát.Použijeme Lillieforsov aj Shapiro-Wilkov test. P hodnoty porovnáme \[H_0: N(\mu,\sigma)\] \[H_1: ¬ H_0\]
lillie.test(jar)##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: jar
## D = 0.25629, p-value = 8.26e-11shapiro.test(jar)##
## Shapiro-Wilk normality test
##
## data: jar
## W = 0.7089, p-value = 1.072e-09lillie.test(leto)##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: leto
## D = 0.18628, p-value = 1.617e-05shapiro.test(leto)##
## Shapiro-Wilk normality test
##
## data: leto
## W = 0.86283, p-value = 6.409e-06Na základe p-hodnôt (menšie ako hladina významnosti) sme oboma testami zamietli H0 o normalite výberu. Dáta nepochádzajú z normálneho rozdelenia preto použijeme neparametrický test.
df<-data.frame(prírastok=c(jar,leto),
Obdobie= rep(c("apríl-máj","jún_júl"),c(length(jar),length(leto))))
ggplot(df,aes(x=Obdobie,y=prírastok))+
geom_boxplot(fill=c("#4472c4", "#FF6666"))+
labs(x= "Obdobie", y= "Denné prírastky" , title = "Porovnanie apríla-maja a juna-jula")+
scale_fill_brewer(palette = "Set1")+
theme(plot.title = element_text(hjust=0.5), legend.position = "none")
Z boxplotu vidime že tieto dve obdobia sa význačne nelíšia. Významnosť budeme testovať vhodným testom.Test vyberieme podľa toho či vyšikmenie bude štatisticky významné alebo nevýznamné.Môžeme testovať Wilcoxonovým 2-výberovým testom alebo Kolmogorov- Smirnov testom.
skewness(jar)## [1] 2.078117skewness(leto)## [1] 0.9950504agostino.test((jar))##
## D'Agostino skewness test
##
## data: (jar)
## skew = 2.0267, z = 5.0227, p-value = 5.094e-07
## alternative hypothesis: data have a skewnessagostino.test(leto)##
## D'Agostino skewness test
##
## data: leto
## skew = 0.97041, z = 2.99001, p-value = 0.00279
## alternative hypothesis: data have a skewness\[H_0: F_x = F_y\] \[H_1:¬ H_0\]
fligner.test(prírastok~ Obdobie,df)##
## Fligner-Killeen test of homogeneity of variances
##
## data: prírastok by Obdobie
## Fligner-Killeen:med chi-squared = 0.12036, df = 1, p-value = 0.7286nulovú hypotézu nemôžem zamietnuť, rozptyly sú rovnaké a preto môžem použiť 2-výberový Wilcoxonov test.
wilcox.test(jar,leto, paired = F)##
## Wilcoxon rank sum test with continuity correction
##
## data: jar and leto
## W = 1910.5, p-value = 0.7994
## alternative hypothesis: true location shift is not equal to 0wilcox.test(prírastok~ Obdobie,df)##
## Wilcoxon rank sum test with continuity correction
##
## data: prírastok by Obdobie
## W = 1910.5, p-value = 0.7994
## alternative hypothesis: true location shift is not equal to 0H0 nemôžeme zamietnuť na hladine významnosti 0.05. Neexistuje štatisticky významný rozdiel v dennom prirastku medzi obdobím apríl-máj a jún-júl.
Pri parametrickej verzii by som ako prvé testovala rovnosť rozptylov F-testom a podlľa toho či by sa rozpyl rovnal alebo nie, by som použila dvojvýberový párový alebo nepárový test.
3 COVID(III)
Rozdeľme pozorované obdobie na tri obdobia. Prvé obdobie nech je od začiatku pozorovaní do 15.5.2020, druhé od 16.5.2020 do 25.8.2020 a posledné obdobie zvyšné pozorovania. Má obdobie v ktorom sledujeme počet úmrtí štatisticky významný vplyv na túto premennú? Aký? Testujte na hladine významnosti 0.1.
prve_obdobie<-data$`Denny prirastok umrti`[1:71]
druhe_obdobie<-data$`Denny prirastok umrti`[72:84]
tretie_obdobie<-data$`Denny prirastok umrti`[85:216]
alfa<-0.1
df1<- data.frame(data)
df2<-data.frame(prve_obdobie)
df3<-data.frame(druhe_obdobie)
df4<-data.frame(tretie_obdobie)
umrtia<-as.factor(df1$Denny.prirastok.umrti) #faktor
df1$Pocet.umrti<-as.factor(df1$Denny.prirastok.umrti)Testujeme či je výber náhodný \[H_0:výber je náhodný\] \[H_1: ¬H_0\]
turning.point.test(prve_obdobie)##
## Turning Point Test
##
## data: prve_obdobie
## statistic = 2.1628, n = 28, p-value = 0.03055
## alternative hypothesis: non randomnessturning.point.test(tretie_obdobie)##
## Turning Point Test
##
## data: tretie_obdobie
## statistic = 3.2835, n = 25, p-value = 0.001025
## alternative hypothesis: non randomnessNulovú hypotézu o náhodnosti nemôžem zamietnuť na hladine významnosti 0.05 v prípade prvom a tretiom období.No v druhom ondobí nulovú hypotézu o náhodnosti zamietam.Čiže V prvom aj tretiom období boli dáta nahodné. V druhom období dáta nie sú náhodné pretože vrámci 16.5-25.8.2020 počet denných prípadov sa rapídne nezvyšoval a možno preto sa aj úmrtosť prudko nezvyšovala ako to bolo vrámci prvého a tretieho obdobia.
g1<-ggplot(df2, aes(prve_obdobie)) +
geom_histogram(aes(y= ..density..), bins=7, fill="#4472c4")+
labs(x="Namerané hodnoty", y="Pravdepodobnosť", title = "Prvé obdobie")+
geom_density(alpha=.4, fill="#FF6666")+
theme(plot.title = element_text(hjust=0.5))
g3<-ggplot(df4, aes(tretie_obdobie)) +
geom_histogram(aes(y= ..density..), bins=7, fill="#4472c4")+
labs(x="Namerané hodnoty", y="Pravdepodobnosť", title = "Druhé obdobie")+
geom_density(alpha=.4, fill="#FF6666")+
theme(plot.title = element_text(hjust=0.5))
plot_grid(g1,g3,labels = "AUTO")
Z histogramu pre prvé obdobie a tretie obdobie vidíme, že hodnoty sú “vynechné” a nie sú ani symetrické
Testujeme či výber pochádza z normálneho rozdelenia
lillie.test(prve_obdobie)##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: prve_obdobie
## D = 0.41851, p-value < 2.2e-16shapiro.test(prve_obdobie)##
## Shapiro-Wilk normality test
##
## data: prve_obdobie
## W = 0.57672, p-value = 5.38e-13lillie.test(tretie_obdobie)##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: tretie_obdobie
## D = 0.50291, p-value < 2.2e-16shapiro.test(tretie_obdobie)##
## Shapiro-Wilk normality test
##
## data: tretie_obdobie
## W = 0.31168, p-value < 2.2e-16Hypotézu o normalite na základe p-hodnôt zamietame pre oba výbery, preto použijeme neparametrický test. Keďže dáta nie sú normálne rozdelené, odľahlé hodnoty hľadáme pomocou medzikvantilovej metódy. Na základe neho vidíme, že náhodný výber neobsahuje žiadne vybočujúce hodnoty a v tretom období obsahuje extremálne hodnoty.
Q1<-quantile(prve_obdobie, probs=0.25) # dolný kvartil
Q3<-quantile(prve_obdobie, probs=0.75) # horný kvartil
IQR<-IQR(prve_obdobie)
k1<-1.5 # pre vybočujúce hodnoty
k2<-3 # pre extremálne hodnoty
prve_obdobie[prve_obdobie<(Q1-k1*IQR)]## numeric(0)prve_obdobie[prve_obdobie>(Q3+k1*IQR)]## [1] 4V prvom období sa nachádzajú 4 vybočujúce hodnoty.
Q1<-quantile(druhe_obdobie, probs=0.25) # dolný kvartil
Q3<-quantile(druhe_obdobie, probs=0.75) # horný kvartil
IQR<-IQR(druhe_obdobie)
k1<-1.5 # pre vybočujúce hodnoty
k2<-3 # pre extremálne hodnoty
druhe_obdobie[druhe_obdobie<(Q1-k1*IQR)]## numeric(0)druhe_obdobie[druhe_obdobie>(Q3+k1*IQR)]## numeric(0)V druhom období sa nachádzajú žiadne vybočujúce hodnoty, pretože toto obdobie bolo celkom kludné vhľadom na počet úmrtí.
Q1<-quantile(tretie_obdobie, probs=0.25) # dolný kvartil
Q3<-quantile(tretie_obdobie, probs=0.75) # horný kvartil
IQR<-IQR(tretie_obdobie)
k1<-1.5 # pre vybočujúce hodnoty
k2<-3 # pre extremálne hodnoty
tretie_obdobie[tretie_obdobie<(Q1-k1*IQR)]## numeric(0)tretie_obdobie[tretie_obdobie>(Q3+k1*IQR)]## [1] 1 2 2 4 1 1 1 1 3 1 3 6 1 2V tretom období sa nachádza najviac extremálych hodnôt. V tomto štádiu úmrtnosť prudko vzrástla. Testujem vplyv jedného faktora(počet úmrtí) na premennú obdobie. Požijeme neparametrickú analýzu rozptylu Kruskal_Wallisov test \[H_0:F_1=F_2=F_3 \] \[H_1: F_1\not= F_2\not= F_3\]
df <- data.frame(umrtia=c(prve_obdobie,druhe_obdobie,tretie_obdobie),
obdobie=rep(c("1","2","3"),c(length(prve_obdobie),length(druhe_obdobie),length(tretie_obdobie))))
ggplot(df, aes(x= obdobie, y=umrtia))+
geom_boxplot(aes(fill=obdobie))+
labs(x= "Obdobia", y= "Počet úmrtí" , title = "Porovnanie 3 období")+
scale_fill_brewer(palette = "Blues", name="Obdobia", labels= c("6.3.2020-15.5.2020","16.5.2020-25.8.2020","26.8.2020-7.10.2020"))+
theme(plot.title = element_text(hjust=0.5))
Sledujeme jeden znak na troch výberoch. Ide o časový rad.Z grafu vidíme, že prvé obdobie má väščí počet denných prírastkov úmrtia ako druhé obdobie. A v tretom období sa vyskytujú najčastejšie extrémne hodnoty ale podobá sa druhému obdobiu. Použitím analýzy rozptylu zistíme, či ide o štatisticky významný rodiel.
Rozhodli sme sa pre neparametrický Kruskal-Wallisov test.
kruskal.test( obdobie~ umrtia, df)##
## Kruskal-Wallis rank sum test
##
## data: obdobie by umrtia
## Kruskal-Wallis chi-squared = 14.714, df = 5, p-value = 0.01166Na základe p-hodnoty Kruskal-Wallisovho testu, ktorá je menšia ako hladina významnosti alfa,nulovú hypotézu zamietame. Existuje teda na hladine významnosti 0.05 štatisticky významný rozdiel v úmrtiami medzi obdobiami. Aby sme zistili medzi ktorými, použujeme Dunnov test.
dunn.test(df$umrtia, df$obdobie) ## Kruskal-Wallis rank sum test
##
## data: x and group
## Kruskal-Wallis chi-squared = 12.1962, df = 2, p-value = 0
##
##
## Comparison of x by group
## (No adjustment)
## Col Mean-|
## Row Mean | 1 2
## ---------+----------------------
## 2 | 2.500656
## | 0.0062*
## |
## 3 | 3.105166 -1.023039
## | 0.0010* 0.1531
##
## alpha = 0.05
## Reject Ho if p <= alpha/2Porovnaním jednotlivých úmrtí v tabuľke vidíme, že štatisticky významný rozdiel je medzi druhým a tretím obdobím, teda medzi 16.5.2020-25.8.2020 a 26.8.2020-7.10.2020.Ako sme už videli z boxplotu Tento rozdiel je zapríčinený vybočujúcimi hodnotami.
Ak by bola normalita splnená použili by sme parametricky test jednofaktorovej analýzy rozptylu.
4 Železo
U 123 pacientov s nízkou hladinou železa v krvi sa sledovala efektivita doplnku stravy so zvýšeným obsahom železa. Na hladine významnosti 0,05 testujte, či je doplnok účinný, teda zvyšuje hladinu železa. Testujte, či pr prvom pozorovaní ide naozaj o nízke hodnoty železa, teda či je stredná hodnota aspoň 4. Pri druhom meraní zistite, či je smerodajná odchýlka nameraných hodnôt väčšia ako 0.5.
data <- read_excel("C:/Users/Simi/Desktop/R_ko(3.semester)/3_zadanie/data_testovanie_zadanie.xlsx", sheet = 2)
fe_1<-data$Fe_1
fe_2<-data$Fe_2
n<-123
alfa<-0.05
mu<-4
sigma<-0.05
df1<-data.frame(data$Fe_1)
df2<-data.frame(data$Fe_2)Testujeme náhodnosť hodôt. Použijem Wald-Wolfowitz runs test. \[H_0:výber je náhodný\] \[H_1: ¬H_0\]
runs.test(fe_1,plot=T)
##
## Runs Test
##
## data: fe_1
## statistic = 1.0909, runs = 68, n1 = 61, n2 = 61, n = 122, p-value =
## 0.2753
## alternative hypothesis: nonrandomnessruns.test(fe_2,plot=T)
##
## Runs Test
##
## data: fe_2
## statistic = -0.18182, runs = 61, n1 = 61, n2 = 61, n = 122, p-value =
## 0.8557
## alternative hypothesis: nonrandomnessNulovú hypotézu o náhodnosti nemôžem zamietnuť na hladine významnosti 0.05, čiže oba súbory dát sú náhodné.
turning.point.test(fe_1)##
## Turning Point Test
##
## data: fe_1
## statistic = 0.71814, n = 123, p-value = 0.4727
## alternative hypothesis: non randomnessturning.point.test(fe_2)##
## Turning Point Test
##
## data: fe_2
## statistic = -1.0054, n = 123, p-value = 0.3147
## alternative hypothesis: non randomnessAni na základe testu kritických bodov nemôžem zamietnuť nulovú hypotézu o náhodnosti.
Ttestujeme normalitu vo všeobecnosti bez použitia parametrov.Použijeme Lillieforsov aj Shapiro-Wilkov test. P hodnoty porovnáme. \[H_0:X ~ N(\mu,\sigma)\] \[H_1: ¬ H_0\]
graf1<-ggplot(df1, aes(data$Fe_1)) +
geom_histogram(aes(y= ..density..), bins=7, fill="#4472c4")+
labs(x="Namerané hodnoty", y="Pravdepodobnosť", title = "Histogram")+
geom_density(alpha=.4, fill="#FF6666")+
theme(plot.title = element_text(hjust=0.5))
graf2<-ggplot(df2, aes(data$Fe_2)) +
geom_histogram(aes(y= ..density..), bins=7, fill="#4472c4")+
labs(x="Namerané hodnoty", y="Pravdepodobnosť", title = "Histogram")+
geom_density(alpha=.4, fill="#FF6666")+
theme(plot.title = element_text(hjust=0.5))
plot_grid(graf1,graf2, labels = "AUTO")
Z histogramov, ktoré u celkom symetrické usudzujeme, že dáta by mohli byť normálne rozdelené.
graf1<-ggplot(df1, aes(sample=data$Fe_1)) +
stat_qq() + stat_qq_line(colour= "#FF6666")+
labs(x= "Teoretický", y = "Nameraný", title = "QQ-plot")+
theme(plot.title = element_text(hjust=0.5))
graf2<-ggplot(df2, aes(sample=data$Fe_2)) +
stat_qq() + stat_qq_line(colour= "#FF6666")+
labs(x= "Teoretický", y = "Nameraný", title = "QQ-plot")+
theme(plot.title = element_text(hjust=0.5))
plot_grid(graf1,graf2,labels = "AUTO")
Aj z kvantil-kvantilového grafu vidíme že, vačšina bodov sa nachádza na priamke, preto usudzujeme že ide naozaj o normálne rozdelenie. Hypotézu otestujeme pomocou Lillieforsovho a Shapiro-Wilkovho testu.
lillie.test(fe_1)##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: fe_1
## D = 0.082024, p-value = 0.04115shapiro.test(fe_1)##
## Shapiro-Wilk normality test
##
## data: fe_1
## W = 0.98579, p-value = 0.227lillie.test(fe_2)##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: fe_2
## D = 0.064082, p-value = 0.2474shapiro.test(fe_2)##
## Shapiro-Wilk normality test
##
## data: fe_2
## W = 0.9853, p-value = 0.2046Na základe p-hodôt ktoré sú väčšie ako alfa, nezamietame nulovú hypotézu o normalite oboch výverov aj napriek tomu, že pre prvý výber vyšiel podľa Lillieforsov testu nie normálne rozdelenný ale kedže Shapiro_Wilkov test má väčšiu váhu, použujem ten.
Kedže máme splnený predpoklad normality, testujeme, či či najmenšia/najväčšia hodnota v dátach je vybočujúca. Extrémalne hodnoty testujem Grubsovým testom.
grubbs.test(fe_1)##
## Grubbs test for one outlier
##
## data: fe_1
## G = 3.15120, U = 0.91794, p-value = 0.08101
## alternative hypothesis: lowest value -0.355089575052261 is an outliergrubbs.test(fe_1,opposite = T)##
## Grubbs test for one outlier
##
## data: fe_1
## G = 3.01338, U = 0.92496, p-value = 0.1335
## alternative hypothesis: highest value 5.89059244096279 is an outliergrubbs.test(fe_2)##
## Grubbs test for one outlier
##
## data: fe_2
## G = 2.90526, U = 0.93025, p-value = 0.1943
## alternative hypothesis: lowest value 1.03126377426088 is an outliergrubbs.test(fe_2,opposite = T)##
## Grubbs test for one outlier
##
## data: fe_2
## G = 2.86879, U = 0.93199, p-value = 0.2199
## alternative hypothesis: highest value 6.97499354928732 is an outlierNa základe p-hodnoty vidíme, že maximalna ani minimálna hodnota zo súborov sa nepotvrdila ako extremálna.
ggplot(data = df1, aes(x = "", y = fe_1)) +
geom_boxplot(fill="#4472c4", outlier.colour = "red")+labs(x="", y="Obsah železa v krvi", title = "Boxplot pre hladina železa v krvi")+
theme(plot.title = element_text(hjust=0.5))+
geom_hline(aes(yintercept=mean(fe_1)), colour="#ff9999")+
geom_hline(aes(yintercept=4), colour="#e60000")
Z boxplotu vidíme, že odhad strednej hodnoty je menší ako testovaná hodnota 4. Otestujeme, či ide o významný rozdiel.
Testujeme, či stredná hodnota je aspoň 4. \[H_0:\mu=4\] \[H_1: \mu > 4\]
t<-(mean(fe_1)-mu)/sigma/sqrt(n)
q<-qt(alfa,df=n-1)
t<q## [1] TRUENulovú hypotézu zamietame na hladine významnosti 0.05.
p<- pt(t, df=n-1)
p<alfa## [1] TRUEAj pomocou p-hodnoty sme dosiahli rovnakého záveru. To znamená, že stredná hodnota hladiny železa v krvi je dosahuje aspoň 4.
Testujeme, či smerodajná odchylka nameraných hodnôt je väčšia ako 0.5. \[H_0:\sigma=0.5\] \[H_1: \sigma>0.5\]
chi2<-varTest(fe_2,alternative = "greater",sigma.squared = 0.5)
chi2$p.value<alfa## [1] TRUENulovú hypotézu zamietame na hladine významnsoti 0.05. smerodajná odchylka pri nameraní hodnôt obsahu železa je väčšia ako 0.5.
Ak by nebola splnená normalita použila by som alternatívu k jednovýberovému parametrickému t-testu (=jednovýberový test o strednej hodnote pri neznámom rozptyle), najprvn by som porovnala zošikmenia a na základe toho, či by boli štatisticky významné alebo nevýznamné by som použila jednovýberový Wilcoxonov test alebo znamienkový test.