#######################################################################################
################ Statistická laboratoř - Zimní škola ##################################
################   Máme data - a co dál? (2. část)  ###################################
############### Adéla Vrtková, Martina Litschmannová ##################################
#######################################################################################

##  Máme data, a co dál? - 1. část
#   1. Spustíme potřebné balíčky, které obsahují další statistické funkce
#   2. Nastavíme pracovní adresář, odkud importujeme data, popř. kam chceme ukládat výstupy
#   3. Importujeme data (z pracovního adresáře, z internetu)
#   4. Pre-processing -> a) Podíváme se na data
#                        b) uložíme si data ve více formátech (každá funkce má "radši" jiný formát)
#   5. Analýza kvalitativních proměnných

##  Máme data, a co dál? - 2. část
#   6. Analýza kvantitativních proměnných
#   7. Identifikace a rozhodnutí o vyloučení/ponechání odlehlých pozorování

#######################################################################################
## Body 1.-4. jsou hotovy  ####################################

data=read.csv2(file="http://am-nas.vsb.cz/lit40/DATA/aku_new.csv")
dataS=stack(data)
colnames(dataS)=c("kap","vyrobce")
dataS=na.omit(dataS)
a=dataS$kap[dataS$vyrobce=="A"]
b=dataS$kap[dataS$vyrobce=="B"]
c=dataS$kap[dataS$vyrobce=="C"]
d=dataS$kap[dataS$vyrobce=="D"]

#######################################################################################
## 6. Explorační analýza a vizualizace kvantitativní proměnné #########################

## Popisná statistika #################################################################
summary(dataS$kap)

# Výpočet průměru jedné proměnné
mean(dataS$kap)
mean(a)

# Pozor na chybějící hodnoty
mean(data$C)
mean(data$C,na.rm=TRUE)
mean(na.omit(data$C))

# Výpočet mediánu jedné proměnné
quantile(dataS$kap,probs=0.5)
quantile(a,probs=0.5)

# Určení rozsahu
length(na.omit(dataS$kap))

#######################################################################################
## funkce, které umožní aplikaci vybrané funkce na sloupce datového rámce, na hodnoty vektoru apod.
x=1:4
#     - lapply - na vstup aplikuje zvolenou funkci, vrátí list (seznam) hodnot
lapply(x,sqrt)
lapply(x,sqrt)[[2]]

#     - sapply - na vstup aplikuje volenou funkci, vrátí vektor hodnot
sapply(x,sqrt)
sapply(x,sqrt)[2]

#     - vapply - stejné jako sapply, akorát s parametrem navíc, který specifikuje výstup
vapply(x,sqrt,numeric(1))
vapply(x,sqrt,numeric(1))[2]

#     - tapply - lze aplikovat na vektor, jemuž je přiřazen faktor rozlišující hodnoty do skupin
#              - použijeme dále 
#######################################################################################

# Výpočet průměrů kapacit dle výrobců
head(dataS)
tapply(dataS$kap, dataS$vyrobce, mean)

# Výpočet mediánu kapacit dle výrobců
tapply(dataS$kap, dataS$vyrobce, quantile, probs=0.5)

# Obdobně lze určit další charakteristiky.

# Pozor! Funkce pro výpočet šikmosti (skewness) a špičatosti (kurtosis) nejsou součástí základního R, najdete je v balíčku moments
install.packages("moments")
library(moments)
tapply(dataS$kap, dataS$vyrobce, skewness)
tapply(dataS$kap, dataS$vyrobce, kurtosis)

# další charakteristiky -> var(), sd(), min(), max(), skewness(), kurtosis() 


## Krabicový graf - angl. boxplot #####################################################
boxplot(a)

# Úprava grafu
boxplot(a,
        main="Nepoužívejme nadpis, když používáme titulek", 
        xlab="Výrobce A",
        ylab="kapacita (mAh)")

# Zabarvení boxplotu, zobrazení průměru
boxplot(a, 
        xlab="Výrobce A",
        ylab="kapacita (mAh)",
        col="grey")
points(1, mean(a,na.rm=TRUE), pch=3)

# Horizontální boxplot
boxplot(a,
        main="Výrobce A",
        horizontal=TRUE,
        xlab="kapacita (mAh)")
points(mean(a,na.rm=TRUE),1, pch=3)

# Další parametry
boxplot(a, 
        xlab="Výrobce A",
        ylab="kapacita (mAh)",
        col="grey",
        outline=FALSE,			    # nezobrazí odlehlá pozorování
        boxwex=0.5)                         # změní šířku krabice na 1/2

# Přidání přímky y=mean(a)
abline(h=mean(a,na.rm=TRUE), col="red",lty=2)

# Funkci boxplot lze využít nejen k vykreslení grafu
boxplot(a,plot=FALSE)

boxplot(a)$out

## Vícenásobný krabicový graf - angl. multiple boxplot ################################
boxplot(data)

boxplot(data, 
        boxwex=0.5,
        col="grey")

boxplot(data, 
        boxwex=c(0.5,0.5,1.5,1.5),
        col=c("red","yellow","grey","blue"))

boxplot(data, 
        boxwex=c(0.5,0.5,1.5,1.5),
        col=terrain.colors(4))
abline(h=1700,col="red",lty=2)

# Mezera v boxplotu - myšlenka - třetí pozice bude vynechána
boxplot(data,
        at=c(1,2,4,5))

# Boxplot konstruován z datové matice
boxplot(dataS$kap~dataS$vyrobce)

# Pořadí krabic dle přání s využitím funkce list
pom=list(a, b, c, d)
boxplot(pom)

pom=list(d, c, b, a)
boxplot(pom, names=c("D","C","B","A"))

## Histogram ##########################################################################
hist(a)
hist(a, breaks=15) # Co dělají různé hodnoty parametru breaks s grafem?
hist(a, main="")

hist(a, 
     main="Histogram pro kapacitu akumulátorů, výrobce A", 
     xlab="kapacita (mAh)",
     ylab="počet akumulátorů",
     col="blue", 
     border="grey",
     labels=TRUE)         # přidá absolutní četnosti daných kategorií ve formě popisků  

hist(a,
     main="Histogram pro kapacitu akumulátorů, výrobce A", 
     xlab="kapacita (mAh)",
     ylab="počet akumulátorů",
     col="blue", 
     border="grey",
     labels=TRUE,
     ylim=c(0,30)) # nastavení rozsahu osy y

hist(a, 
     main="Histogram pro kapacitu akumulátorů, výrobce A", 
     xlab="kapacita (mAh)",
     ylab="f(x)",
     col="blue", 
     border="grey",
     freq=FALSE,        # změna měřítka na ose y --> f(x)
     ylim=c(0,0.018))	       
lines(density(a))        # připojí graf odhadu hustoty pravděpodobnosti

#######################################################################################
## Více grafů do jednoho obrázku -> funkce layout nebo par ############################
## Podívejte se na možnosti kombinování grafů - http://www.statmethods.net/advgraphs/layout.html

## Kombinace histogramu a boxplotu ####################################################
pom=layout(mat = matrix(1:8,2,4, byrow=FALSE), height = c(2.5,1))
layout.show(pom)
par(oma=c(2,2,3,2),mar=c(2,2,3,2))

hist(a, 
     main="Výrobce A",
     xlab="kapacita (mAh)", 
     ylab="četnost", 
     ylim=c(0,32), 
     xlim=c(min(dataS$kap),max(dataS$kap)))
boxplot(a, 
        horizontal=TRUE, 
        ylim=c(min(dataS$kap),max(dataS$kap)), 
        boxwex=1.5)

hist(b, 
     main="Výrobce B", 
     xlab="kapacita (mAh)", 
     ylab="četnost", 
     ylim=c(0,32), 
     xlim=c(min(dataS$kap),max(dataS$kap)))
boxplot(b, 
        horizontal=TRUE, 
        ylim=c(min(dataS$kap),max(dataS$kap)), 
        boxwex=1.5)

hist(c, 
     main="Výrobce C", 
     xlab="kapacita (mAh)", 
     ylab="četnost", 
     ylim=c(0,32), 
     xlim=c(min(dataS$kap),max(dataS$kap)))
boxplot(c, 
        horizontal=TRUE, 
        ylim=c(min(dataS$kap),max(dataS$kap)), 
        boxwex=1.5)

hist(d, 
     main="Výrobce D", 
     xlab="kapacita (mAh)", 
     ylab="četnost", 
     ylim=c(0,32), 
     xlim=c(min(dataS$kap),max(dataS$kap)))
boxplot(d, 
        horizontal=TRUE, 
        ylim=c(min(dataS$kap),max(dataS$kap)), 
        boxwex=1.5)
mtext("Histogramy a boxploty pro jednotlivé výrobce", cex = 1.5, outer=TRUE, side=3)

# Pro pokročilé - pomocí for-cyklu
pom=layout(mat = matrix(1:8,2,4, byrow=FALSE), height = c(2.5,1))
layout.show(pom)
par(oma=c(2,2,3,2), mar=c(2,2,3,2)) 
# parametry oma, mar - viz https://www.r-graph-gallery.com/74-margin-and-oma-cheatsheet/

for (i in 1:4){
  hist(data[,i], 
       main=paste("Výrobce",colnames(data)[i]), 
       xlab="kapacita (mAh)", 
       ylab="četnost", 
       xlim=c(min(data,na.rm=TRUE), max(data,na.rm=TRUE)), 
       ylim=c(0,32))
  boxplot(data[,i], 
          horizontal=TRUE, 
          ylim=c(min(data,na.rm=TRUE), max(data,na.rm=TRUE)), 
          boxwex=1.5)
}
mtext("Histogramy a boxploty pro jednotlivé výrobce", cex = 1.5, outer=TRUE, side=3)

#######################################################################################
## 7. Identifikace odlehlých pozorování (a jejich odstranění z dat. rámce)#############

## S individuálním posouzením, jak naložit s odlehlými hodnotami ######################
# Seřázení datového rámce dataS podle proměnnýchdataS$vyrobce a dataS$kap
par(mfrow = c(1,1)) # jednoduché nastavení grafického okna - 1 řádek, 1 sloupec
boxplot(data,plot=FALSE)
boxplot(data)

boxplot(dataS$kap~dataS$vyrobce,plot=FALSE)
boxplot(dataS$kap~dataS$vyrobce)

dataS$kap.bez=dataS$kap # není dobré si přepisovat původní data
dataS$kap.bez[dataS$kap>=1866.9 & dataS$vyrobce=="A"]=NA

dataS$kap.bez