#######################################################################################
################ Statistická laboratoř - 10. workshop #################################
############### Adéla Vrtková, Martina Litschmannová ##################################
#######################################################################################

### Skript je určen pro účastníky workshopů pořádaných Statistickou laboratoří.
### V případě, že jej, popř. jeho části, chcete v nezměněné či mírně upravené podobě
### používat v rámci vlastních publikací nebo ve výukových kurzech realizovaných 
### na vašem pracovišti, budeme rády, když nás o tom budete informovat a v příslušných
### materiálech uvedete poděkování Statistické laboratoři (http://k470.vsb.cz/statlab/).
### Za akceptování této prosby předem děkujeme!

##################################################################################
#### Akumulátory - dvouvýběrové testy #############################################

## Už umíme - načtení dat, vynechání chybějících hodnot
data=read.csv("http://am-nas.vsb.cz/lit40/DATA/aku2.csv",sep=";", dec=".", header=TRUE)
data=na.omit(data)

# zajímá nás pouze výrobce A a C
data.ac=subset(data, data$Výrobce %in% c("A","C"))
data.ac$Výrobce=factor(data.ac$Výrobce)

## Testování hypotéz předchází explorační analýza
library(moments)

tapply(data.ac$pokles,data.ac$Výrobce,summary,na.rm=TRUE)
tapply(data.ac$pokles,data.ac$Výrobce,sd,na.rm=TRUE)

tapply(data.ac$pokles,data.ac$Výrobce,skewness,na.rm=TRUE)
tapply(data.ac$pokles,data.ac$Výrobce,kurtosis,na.rm=TRUE)-3

# prvotní náhled poskytnou krabicové grafy, histogramy, QQ-grafy
pom=layout(matrix(c(1,2,3,1,4,5),3,2))
layout.show(pom)

boxplot(split(data.ac$pokles,data.ac$Výrobce),
        ylab="pokles kapacity baterií (mAh)")
hist(data.ac$pokles[data.ac$Výrobce=="A"],
     ylab="četnost",
     xlab="pokles kapacity baterií (mAh)",
     xlim=c(100,400),
     ylim=c(0,40),
     main="Výrobce A")
hist(data.ac$pokles[data.ac$Výrobce=="C"],
     ylab="četnost",
     xlab="pokles kapacity baterií (mAh)",
     xlim=c(100,400),
     ylim=c(0,70),
     main="Výrobce C")
qqnorm(data.ac$pokles[data.ac$Výrobce=="A"],
       xlab="teoretické kvantily",
       ylab="výběrové kvantily",
       main="Výrobce A")
qqline(data.ac$pokles[data.ac$Výrobce=="A"])
qqnorm(data.ac$pokles[data.ac$Výrobce=="C"],
       xlab="teoretické kvantily",
       ylab="výběrové kvantily",
       main="Výrobce C")
qqline(data.ac$pokles[data.ac$Výrobce=="C"])

## u poklesů výrobce C se nachází odlehlé pozorování -> odstraníme a znovu se na data podíváme
pom = boxplot(data.ac$pokles[data.ac$Výrobce=="C"], plot=FALSE)
data.ac.bez = data.ac
data.ac.bez$pokles[data.ac$pokles %in% pom$out] = NA

pom=layout(matrix(c(1,2,3,1,4,5),3,2))
layout.show(pom)

boxplot(split(data.ac.bez$pokles,data.ac.bez$Výrobce),
        ylab="pokles kapacity baterií (mAh)")
hist(data.ac.bez$pokles[data.ac.bez$Výrobce=="A"],
     ylab="četnost",
     xlab="pokles kapacity baterií (mAh)",
     xlim=c(100,400),
     ylim=c(0,40),
     main="Výrobce A")
hist(data.ac.bez$pokles[data.ac.bez$Výrobce=="C"],
     ylab="četnost",
     xlab="pokles kapacity baterií (mAh)",
     xlim=c(100,400),
     ylim=c(0,70),
     main="Výrobce C")
qqnorm(data.ac.bez$pokles[data.ac.bez$Výrobce=="A"],
       xlab="teoretické kvantily",
       ylab="výběrové kvantily",
       main="Výrobce A")
qqline(data.ac.bez$pokles[data.ac.bez$Výrobce=="A"])
qqnorm(data.ac.bez$pokles[data.ac.bez$Výrobce=="C"],
       xlab="teoretické kvantily",
       ylab="výběrové kvantily",
       main="Výrobce C")
qqline(data.ac.bez$pokles[data.ac.bez$Výrobce=="C"])

## přejdeme k testování normality
tapply(data.ac.bez$pokles,data.ac.bez$Výrobce,shapiro.test)  # usnadnění práce pomocí tapply, aplikujeme Shapirův-Wilkův test -> normalita se nezamítá

## ověření shody rozptylů -> funkce v R - var.test
var.test(data.ac.bez$pokles~data.ac.bez$Výrobce)  # hypotézu o shodě rozptylů nelze zamítnout

## nezávislost + normalita + shoda rozptylů -> dvouvýběrový Studentův t-test
t.test(data.ac.bez$pokles~data.ac.bez$Výrobce,
       var.equal=TRUE, # tímto parametrem rozlišujeme mezi dvouvýběrovým t testem a Aspinové - Welchovým testem
       mu=0,
       alternative="less",
       conf.level=0.95)

# trochu jinak - máme dva sloupce zvlášť - využijeme např. nemáme-li vytvořený standartní datový formát
poklesA=data.ac.bez$pokles[data.ac.bez$Výrobce=="A"]
poklesC=data.ac.bez$pokles[data.ac.bez$Výrobce=="C"]

t.test(poklesA,poklesC,
       var.equal=TRUE, # tímto parametrem rozlišujeme mezi dvouvýběrovým t testem a Aspinové - Welchovým testem
       mu=0,
       alternative="less",
       conf.level=0.95)

## Závěr: Na hladině významnosti 5% zamítáme nulovou hypotézu ve prospěch alternativy. 
## Lze tedy předpokládat, že tvrzení výrobce A je pravdivé - tj. že akumulátory výrobce A mají menší pokles kapacity.

############################################################################################################
############ Postoupíme ve zobecnění o krok dále, aneb přidáme si do srovnání dalšího výrobce ##############
# Jak jsou na tom výrobci A, C a D co se týče poklesu kapacity? Jsou na tom výrobci stejně? ################

## -> Vícevýběrové testy
data=read.csv("http://am-nas.vsb.cz/lit40/DATA/aku2.csv",sep=";", dec=".", header=TRUE)
data=na.omit(data)

# zajímají nás pouze výrobci A, C a D
data.acd=subset(data, data$Výrobce %in% c("A","C","D"))
data.acd$Výrobce=factor(data.acd$Výrobce)  # nutno u daného faktoru odstranit úrovně, které se v dané proměnné již nevyskytují

## opět předchází explorační analýza
tapply(data.acd$pokles,data.acd$Výrobce,summary,na.rm=TRUE)
tapply(data.acd$pokles,data.acd$Výrobce,sd,na.rm=TRUE)

tapply(data.acd$pokles,data.acd$Výrobce,skewness,na.rm=TRUE)
tapply(data.acd$pokles,data.acd$Výrobce,kurtosis,na.rm=TRUE)-3

## pro grafické zobrazení více výběrů je nejvhodnější vícenásobný boxplot - poskytuje přehledné srovnání
par(mfrow=c(1,1))
boxplot(data.acd$pokles~data.acd$Výrobce,
        ylab="poklesy kapacit (mAh)")

# tentokráte se rozhodneme, že odlehlá pozorování v datech ponecháme -> pokračujeme dále v analýze
# chceme-li opět různé grafické nástroje do jednoho obrázku -> funkce layout
pom = layout(matrix(c(1,2,3,4,5,6),3,2))
layout.show(pom)
hist(data.acd$pokles[data.acd$Výrobce=="A"],
     main="A",
     xlab="pokles kapacit (mAh)",
     ylab="četnost",
     xlim=c(100,380))
hist(data.acd$pokles[data.acd$Výrobce=="C"],
     main="C",
     xlab="pokles kapacit (mAh)",
     ylab="četnost",
     xlim=c(100,380))
hist(data.acd$pokles[data.acd$Výrobce=="D"],
     main="D",
     xlab="pokles kapacit (mAh)",
     ylab="četnost",
     xlim=c(100,380))
qqnorm(data.acd$pokles[data.acd$Výrobce=="A"],
       main="",
       xlab="teoretické kvantily",
       ylab="výběrové kvantily")
qqline(data.acd$pokles[data.acd$Výrobce=="A"])
qqnorm(data.acd$pokles[data.acd$Výrobce=="C"],
       main="",
       xlab="teoretické kvantily",
       ylab="výběrové kvantily")
qqline(data.acd$pokles[data.acd$Výrobce=="C"])
qqnorm(data.acd$pokles[data.acd$Výrobce=="D"],
       main="",
       xlab="teoretické kvantily",
       ylab="výběrové kvantily")
qqline(data$pokles[data$Výrobce=="D"])

## přichází na řadu ověření normality - viděli jsme šikmost, špičatost i QQ-grafy,...
# ...nakonec normalitu ověříme exaktně pomocí Shapirova-Wilkova testu
tapply(data.acd$pokles,data.acd$Výrobce,shapiro.test)

## jelikož je předpoklad normality splněn, bude použit pro test o shodě rozptylů Bartlettův test
install.packages("stats")
library(stats)
bartlett.test(data.acd$pokles~data.acd$Výrobce)

## vzhledem k p-hodnotě Bartlettova testu, na hladině významnosti 0.05 nezamítáme nulovou hypotézu o shodě rozptylů
## Nezávislost + normalita + homoskedasticita -> ANOVA

anova = aov(data.acd$pokles~data.acd$Výrobce)
summary(anova)

# na hladině významnosti 0.05 zamítáme nulovou hypotézu o shodě středních hodnot ve prospěch alternativy,
# lze tedy předpokládat, že se kapacita po 100 cyklech bude mezi výrobci lišit
# víme, že se liší - ale které? -> Post hoc analýza -> mnohonásobné porovnávání

TukeyHSD(anova)

# zajímá nás zejména poslední sloupec s p-hodnotami - je-li p-hodnota menší než hladina významnosti,
# pak považujeme na dané hladině významnosti rozdíl mezi danou dvojicí za statisticky významný

# Závěr: Statisticky významně se liší všichni výrobci, tzn. pokles kapacity závisí na tom, od jakého výrobce
# akumulátor pochází.

############################################################################################################
############ Na závěr se podíváme na všechny výrobce - A, B, C, D ##########################################
# Jak jsou na tom výrobci co se týče poklesu kapacity? Jsou na tom výrobci stejně? #########################
## -> Vícevýběrové testy podruhé

data=read.csv("http://am-nas.vsb.cz/lit40/DATA/aku2.csv",sep=";", dec=".", header=TRUE)
data=na.omit(data)

## opět předchází explorační analýza, náhled na data a rozhodnutí o odlehlých pozorováních
tapply(data$pokles,data$Výrobce,summary,na.rm=TRUE)
tapply(data$pokles,data$Výrobce,sd,na.rm=TRUE)

tapply(data$pokles,data$Výrobce,skewness,na.rm=TRUE)
tapply(data$pokles,data$Výrobce,kurtosis,na.rm=TRUE)-3

# vícenásobný krabicový graf
par(mfrow=c(1,1))
boxplot(data$pokles~data$Výrobce,
        ylab="poklesy kapacit (mAh)")

# odlehlá pozorování tentokrát ponecháme v datovém souboru a podíváme se na histogramy a QQ-grafy
pom = layout(matrix(c(1,2,3,4,5,6,7,8),4,2))
layout.show(pom)
hist(data$pokles[data$Výrobce=="A"],main="A",xlab="pokles kapacit (mAh)",ylab="četnost",xlim=c(100,380))
hist(data$pokles[data$Výrobce=="B"],main="B",xlab="pokles kapacit (mAh)",ylab="četnost",xlim=c(100,380))
hist(data$pokles[data$Výrobce=="C"],main="C",xlab="pokles kapacit (mAh)",ylab="četnost",xlim=c(100,380))
hist(data$pokles[data$Výrobce=="D"],main="D",xlab="pokles kapacit (mAh)",ylab="četnost",xlim=c(100,380))

qqnorm(data$pokles[data$Výrobce=="A"],main="",xlab="teoretické kvantily",ylab="výběrové kvantily")
qqline(data$pokles[data$Výrobce=="A"])
qqnorm(data$pokles[data$Výrobce=="B"],main="",xlab="teoretické kvantily",ylab="výběrové kvantily")
qqline(data$pokles[data$Výrobce=="B"])
qqnorm(data$pokles[data$Výrobce=="C"],main="",xlab="teoretické kvantily",ylab="výběrové kvantily")
qqline(data$pokles[data$Výrobce=="C"])
qqnorm(data$pokles[data$Výrobce=="D"],main="",xlab="teoretické kvantily",ylab="výběrové kvantily")
qqline(data$pokles[data$Výrobce=="D"])

## přichází na řadu ověření normality - viděli jsme šikmost, špičatost i QQ-grafy,...
# ...nakonec normalitu ověříme exaktně pomocí Shapirova-Wilkova testu
tapply(data$pokles,data$Výrobce,shapiro.test)

# tentokráte zamítáme předpoklad normality a proto volíme Kruskalův-Wallisův test alias neparametrickou ANOVU
kruskal.test(data$pokles~data$Výrobce)

# na hladině významnosti 0.05 zamítáme nulovou hypotézu o shodě mediánů ve prospěch alternativy
# lze tedy předpokládat, že se poklesy kapacity budou mezi výrobci lišit

# víme, že se liší - ale které? -> Post hoc analýza -> mnohonásobné porovnávání Dunnové metodou
install.packages("dunn.test")
library(dunn.test)
dunn.test(data$pokles,data$Výrobce)

# Závěr: Statisticky významně se liší všichni výrobci, tzn. pokles kapacity závisí na tom, od jakého výrobce
# akumulátor pochází.