#######################################################################################
################ Statistická laboratoř - Zimní škola ##################################
################   ANOVA, Kontingenční tabulky  ###################################
############### Adéla Vrtková, Martina Litschmannová ##################################
#######################################################################################
library(moments)

######################################################################################
############ ANOVA ##############
# Jak jsou na tom výrobci A, C a D co se týče kapacity? Jsou na tom výrobci stejně? 

## -> Vícevýběrové testy
data=read.csv2("http://am-nas.vsb.cz/lit40/DATA/aku_new.csv", header=TRUE)
dataS=stack(data)
colnames(dataS)=c("kap","vyrobce")
dataS=na.omit(dataS)

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

## samotnému testování předchází explorační analýza
tapply(data.acd$kap,data.acd$vyrobce,summary,na.rm=TRUE)
tapply(data.acd$kap,data.acd$vyrobce,sd,na.rm=TRUE)

tapply(data.acd$kap,data.acd$vyrobce,skewness,na.rm=TRUE)
tapply(data.acd$kap,data.acd$vyrobce,moments::kurtosis,na.rm=TRUE)-3 
# napíšeme-li název funkce ve formátu "název balíčku::funkce" (např. moments::kurtosis(a)),
# bude volána funkce z dané knihovny (používáme k zajištění volání funkce kurtosis
# z balíčku moments)

## 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$kap~data.acd$vyrobce,
        ylab="kapacita (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 pro posouzení normality do jednoho obrázku -> funkce layout
par(mfcol = c(3,2))  # pořadí vykreslovaných grafů - po sloupcích
hist(data.acd$kap[data.acd$vyrobce=="A"],
     main="A",
     xlab="kapacita (mAh)",
     ylab="četnost",
     xlim=c(1550,1900),
     ylim=c(0,25))
hist(data.acd$kap[data.acd$vyrobce=="C"],
     main="C",
     xlab="kapacita (mAh)",
     ylab="četnost",
     xlim=c(1550,1900),
     ylim=c(0,25))
hist(data.acd$kap[data.acd$vyrobce=="D"],
     main="D",
     xlab="kapacita (mAh)",
     ylab="četnost",
     xlim=c(1550,1900),
     ylim=c(0,25))
qqnorm(data.acd$kap[data.acd$vyrobce=="A"],
       main="")
qqline(data.acd$kap[data.acd$vyrobce=="A"])
qqnorm(data.acd$kap[data.acd$vyrobce=="C"],
       main="")
qqline(data.acd$kap[data.acd$vyrobce=="C"])
qqnorm(data.acd$kap[data.acd$vyrobce=="D"],
       main="")
qqline(data.acd$kap[data.acd$vyrobce=="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$kap,data.acd$vyrobce,shapiro.test)

## jelikož předpoklad normality nelze zamítnout, bude použit pro test o shodě rozptylů Bartlettův test
par(mfrow = c(1,1))
boxplot(data.acd$kap~data.acd$vyrobce,
        ylab="kapacita (mAh)")
s2 = tapply(data.acd$kap,data.acd$vyrobce,var)
max(s2)/min(s2)

install.packages("stats")
library(stats)
bartlett.test(data.acd$kap~data.acd$vyrobce)

## 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$kap~data.acd$vyrobce)
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ýrobci C a D od výrobce A, 
# výrobci C a D se statisticky významně neliší.
# Lze usoudit, že očekávaná kapacita akumulátoru záleží na tom, od jakého výrobce akumulátor pochází.

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

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

tapply(dataS$kap,dataS$vyrobce,skewness,na.rm=TRUE)
tapply(dataS$kap,dataS$vyrobce,moments::kurtosis,na.rm=TRUE)-3

# vícenásobný krabicový graf
par(mfrow=c(1,1))
boxplot(dataS$kap~dataS$vyrobce,
        ylab="kapacita (mAh)")

# odlehlá pozorování tentokrát ponecháme v datovém souboru a podíváme se na histogramy a QQ-grafy
par(mfcol = c(4,2))  # pořadí vykreslovaných grafů - po sloupcích
hist(dataS$kap[dataS$vyrobce=="A"],main="A",xlab="kapacita (mAh)",ylab="četnost",xlim=c(1550,1900),ylim=c(0,25))
hist(dataS$kap[dataS$vyrobce=="B"],main="B",xlab="kapacita (mAh)",ylab="četnost",xlim=c(1550,1900),ylim=c(0,25))
hist(dataS$kap[dataS$vyrobce=="C"],main="C",xlab="kapacita (mAh)",ylab="četnost",xlim=c(1550,1900),ylim=c(0,25))
hist(dataS$kap[dataS$vyrobce=="D"],main="D",xlab="kapacita (mAh)",ylab="četnost",xlim=c(1550,1900),ylim=c(0,25))

qqnorm(dataS$kap[dataS$vyrobce=="A"],main="")
qqline(dataS$kap[dataS$vyrobce=="A"])
qqnorm(dataS$kap[dataS$vyrobce=="B"],main="")
qqline(dataS$kap[dataS$vyrobce=="B"])
qqnorm(dataS$kap[dataS$vyrobce=="C"],main="")
qqline(dataS$kap[dataS$vyrobce=="C"])
qqnorm(dataS$kap[dataS$vyrobce=="D"],main="")
qqline(dataS$kap[dataS$vyrobce=="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(dataS$kap,dataS$vyrobce,shapiro.test)

# Předpoklad normality nelze zamítnout - opět otestujeme shodu rozptylů.
bartlett.test(dataS$kap,dataS$vyrobce)

# Tentokráte zamítáme shodu rozptylů, a proto volíme Kruskalův-Wallisův test alias neparametrickou ANOVU.
kruskal.test(dataS$kap~dataS$vyrobce)

# Na hladině významnosti 0.05 zamítáme nulovou hypotézu o shodě mediánů, 
# lze tedy předpokládat, že se mediány kapacit 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(dataS$kap,dataS$vyrobce)

# Závěr: Mediány kapacit akumulátorů výrobce A se statisticky významně 
# liší od mediánů kapacit akumulátorů výrobců B, C i D. 
# Mediány kapacit akumulátorů výrobců B, C a D se statisticky významně neliší.


######################################################################################
#### Akumulátory - kontingenční tabulky ##############################################

# Podívejme se pouze na výrobce B a D.
data.bd=subset(dataS, dataS$vyrobce %in% c("B","D"))
data.bd$vyrobce=factor(data.bd$vyrobce)

# Zajímá nás, zda kapacita akumulátoru přesáhla 1650 mAh.
# definování nové dichotmické proměnné pomocí podmíněného příkazu - velmi podobné excelovskému "když"
data.bd$ind=ifelse(data.bd$kap>1650, # podmínka
                   "ANO", # co se má stát, je-li podmínka splněna
                   "NE") # není-li podmínka splněna

### Prezentace kategoriální proměnné pomocí známé funkce table ################################
## funkci table už známe z explorační analýzy kvalitativní proměnné - používali jsme ji pro jednu proměnnou
table(data.bd$ind)
table(data.bd$vyrobce)

## pomocí funkce table získáme kontingenční tabulku
tab=table(data.bd$vyrobce,data.bd$ind)  # Pozor na pořadí!!! Dodržujme "příčina"-"důsledek"
tab

prop.table(tab)          # relativní četnosti
prop.table(tab)*100      # relativní četnosti v %

prop.table(tab,1)        # řádkové relativní četnosti
prop.table(tab,2)        # sloupcové relativní četnosti

## grafická prezentace pomocí mozaikového grafu
par(mfrow=c(1,1))
mosaicplot(tab)

mosaicplot(tab,                     
           color = terrain.colors(2),
           main = "",
           cex=1.5)                      # nastavení velikosti písma                    

install.packages("lsr")
library(lsr)
cramersV(tab) # Míra kontingence Cramerovo V naznačuje slabou závislost

### Test nezávislosti v kont. tabulce a míra kontingence
# H0: Výrobce a výskyt kapacity větší než 1650 mAh jsou statisticky nezávislé veličiny.
# HA: Výrobce a výskyt kapacity větší než 1650 mAh nejsou statisticky nezávislé veličiny.

chisq.test(tab)

# Pozor - funkce sama nezkontroluje předpoklady testu!
pom=chisq.test(tab)
pom$expected 
# očekávané četnosti>=2, 
# alespoň 80% pak >5 
# - splněno -> lze testovat korektně -> dívám se na p-hodnotu
pom$p.value

## Na hladině významnosti 0.05 nezamítám H0, tzn. výrobce a výskyt kapacity přesahující 1650 mAh  
# jsou statisticky nezávislé veličiny. Tj. z daného hlediska nezáleží na tom, od jakého výrobce akumulátory pocházejí.


##### Relativní riziko a poměr šancí
# Z grafu to vypadá, že výrobce B je na tom lépe - lze kvantifikovat kolikrát a jestli vůbec stat. významně?
tab
install.packages("epiR")
library(epiR)
epi.2by2(tab)

## Interpretace RR: Výrobce B má cca 1.3x větší riziko (pravděpodobnost), že jeho akumulátory 
##                  budou mít kapacity přes 1650 mAh, než výrobce D.

## Intervalový odhad RR napoví statistickou významnost.
## S 95% spolehlivostí lze očekávat, že výrobce B má 0.93 až 1.70 krát větší riziko, že jeho akumulátory 
## budou mít kapacity přes 1650 mAh, než výrobce D.

## Interpretace OR: Výrobce B má cca 1.6x větší šanci, že jeho akumulátory 
##                  budou mít kapacity přes 1650 mAh, než výrobce D.

## Intervalový odhad OR napoví statistickou významnost.
## S 95% spolehlivostí lze očekávat, že výrobce B má 0.87 až 2.95 krát větší šanci, že jeho akumulátory 
## budou mít kapacity přes 1650 mAh, než výrobce D.

# Jelikož intervalové odhady RR i OR obsahují jedničku - jsou stat. nevýznamné
# Nelze říci, že výrobce B (popř. D) má statisticky významně zvýšenou/sníženou šanci/riziko (pravděpodobnost) na přesáhnutí kapacity 1650 mAh.