#######################################################################################
################ Statistická laboratoř - 9. workshop ##################################
############### Adéla Vrtková, Martina Litschmannová ##################################
#######################################################################################
library(moments)

#######################################################################################
### Životnost monitorů aneb v praxi vycházíme z dat ###################################

data=read.csv (file="http://am-nas.vsb.cz/lit40/DATA/monitory.csv")

## Testování hypotézy předchází explorační analýza
summary(data$zivotnost) # charakteristiky polohy
sd(data$zivotnost)  # směrodatná odchylka

pom=layout(matrix(c(1,2,3,3), 2, 2), widths=c(1,1), heights=c(1,1))
layout.show(pom)
boxplot(data$zivotnost,   # z krabicového grafu získáme představu o odlehlých pozorováních a symetrii
        main="Krabicový graf pro životnost monitorů",
        horizontal = TRUE,
        ylim=c(200,2200),
        xlab="životnost (hod)")
hist(data$zivotnost,   # histogram rovněž napoví symetrii, zešikmení, díváme se zda připomíná "Gaussovský klobouk"
     main="Histogram pro životnost monitorů",
     xlim=c(200,2200),
     ylab="četnost",
     xlab="životnost (hod)")
qqnorm(data$zivotnost,        # QQ-grafy dají první náhled na možnou normalitu dat
       xlab="teoretické kvantily",
       ylab="výběrové kvantily",
       main="Q-Q graf") 
qqline(data$zivotnost,
       col="red")

## Data obsahují odlehlá pozorování -> vyloučíme je a znovu se na data podíváme
pom = boxplot(data$zivotnost, plot=FALSE)
data$zivotnost.bez = data$zivotnost
data$zivotnost.bez[data$zivotnost %in% pom$out]=NA

summary(data$zivotnost.bez)
sd(data$zivotnost.bez, na.rm=TRUE)

pom=layout(matrix(c(1,2,3,3), 2, 2), widths=c(1,1), heights=c(1,1))
layout.show(pom)
boxplot(data$zivotnost.bez,   # z krabicového grafu získáme představu o odlehlých pozorováních a symetrii
        main="Krabicový graf pro životnost monitorů",
        horizontal = TRUE,
        ylim=c(200,2200),
        xlab="životnost (hod)")
hist(data$zivotnost.bez,   # histogram rovněž napoví symetrii, zešikmení, díváme se zda připomíná "Gaussovský klobouk"
     main="Histogram pro životnost monitorů",
     xlim=c(200,2200),
     ylab="četnost",
     xlab="životnost (hod)")
qqnorm(data$zivotnost.bez,        # QQ-grafy dají první náhled na možnou normalitu dat
       xlab="teoretické kvantily",
       ylab="výběrové kvantily",
       main="Q-Q graf") 
qqline(data$zivotnost.bez,
       col="red")

## pokračujeme dále - potřebujeme ověřit normalitu dat (zatím pouze na základě explorační analýzy)

skewness(data$zivotnost.bez, na.rm=TRUE)   # šikmost
kurtosis(data$zivotnost.bez, na.rm=TRUE)-3 # špičatost

# Ověření normality pomocí exaktního testu -> Shapirův-Wilkův test
# H0: Životnost monitorů má normální rozdělení.
# HA: Životnost monitorů nemá normální rozdělení.

shapiro.test(data$zivotnost.bez)  

# Na hladině významnosti 0,05 nelze předpoklad normality zamítnout
# (p - hodnota = 0,705, Shapirův-Wilkův test).

# normalita dat ověřena -> Studentův t-test

t.test(data$zivotnost.bez,    # data, která chceme otestovat
       mu=1200,               # definování nulové hypotézy
       alternative="greater", # specifikace alternativy
       conf.level=0.95)       # stanovení spolehlivosti = 1-hladina významnosti

## Závěr: Zamítáme nulovou hypotézu ve prospěch alternativy na hladině testu 0.05
## (p - hodnota = 0,006, t-test).
## Můžeme tedy předpokládát statisticky významné zlepšení životnosti monitorů.


#######################################################################
##### Ach, ty volby... ################################################
n.ucast.CVVM=round(0.562*528) # počet rozhodnutých voličů ve výb. šetření CVVM
x.KDUSTAN.CVVM=round(0.095*n.ucast.CVVM)  # odhad počtu voličů koalice KDU-STAN dle CVVM
x.ANO.CVVM=round(0.335*n.ucast.CVVM)      # odhad počtu voličů ANO dle CVVM

## oboustranný 95% intervalový odhad volebních preferencí KDU-ČSL a STAN
binom.test(x.KDUSTAN.CVVM, # počet hlasů pro koalici
           n.ucast.CVVM) # počet všech hlasů

# S 95% spolehlivostí lze volební preference koalice KDU-ČSL a STAN odhadovat 
# mezi 6,4% a 13,3% (Clopperův – Pearsonův odhad).

## Lze tvrdit, že volební preference koalice KDU-ČSL a STAN nepřekročí 10 %? 
# testujeme, H0: pi=0.1, proti HA: pi<0.1
# tj. zda můžeme předpokládat, že koalice bude mít 10% hlasů, proti alternativě, že na ně nedosáhne

binom.test(x.KDUSTAN.CVVM,        # počet hlasů pro koalici
           n.ucast.CVVM,        # počet všech hlasů
           p=0.1,               # nulový hypotéza - "testovaná pravděpodobnost"
           alternative="less")  # upřesnění jednostranné alternativy

# Na hladině významnosti 5% nelze tvrdit, že volební preference koalice 
# KDU-ČSL a STAN jsou nižší než 10% (p - hodnota =0,418, Clopperův – Pearsonův test).

## Liší se odhad volebních preferencí ANO získaný CVVM a SANEP statisticky významně?

n.ucast.SANEP=round(0.562*2407) # počet rozhodnutých voličů ve výb. šetření SANEP
x.ANO.SANEP=round(0.275*n.ucast.SANEP)


# oboustranný intervalový odhad pro ANO - CVVM
binom.test(x.ANO.CVVM,   # počet hlasů pro ANO
           n.ucast.CVVM) # počet všech hlasů

# oboustranný intervalový odhad pro ANO - SANEP
binom.test(x.ANO.SANEP,   # počet hlasů pro ANO
           n.ucast.SANEP) # počet všech hlasů

# intervalový odhad pro rozdíl volebních preferec
prop.test(c(x.ANO.CVVM,x.ANO.SANEP),c(n.ucast.CVVM,n.ucast.SANEP))

# Se spolehlivostí 95 % lze rozdíl mezi volebními preferencemi ANO 
# zjištěnými CVVM a SANEP očekávat mezi -0,2 % a 11,9 %, 
# tj. zjištěné výsledky se statisticky významně neliší (asymptotický interval spolehlivosti).

# Na hladině významnosti 5% nelze zamítnout hypotézu o shodě volebních preferencí 
# zjištěných CVVM a SANEP vůči alternativě o tom, že se zjištěné volební preference 
# liší (p-hodnota = 0,052, asymptotický test s korekcí na spojitost).

######################################################################
#### Zpět k akumulátorům #############################################

# 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"))
boxplot(split(data.ac$pokles,data.ac$ Výrobce))

data.ac$Výrobce=factor(data.ac$Výrobce)  # nutno u daného faktoru odstranit úrovně, které se v dané proměnné již nevyskytují

# 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,
       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,
       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 (p-hodnota << 0,001, dvouvýběrový t-test. 
## Lze tedy předpokládat, že tvrzení výrobce A je pravdivé - tj. že akumulátory výrobce A mají menší pokles kapacity.