setwd('C:/Disk D/ND-Skola/02-Vyuka/04-Aplikovana statistika 2016/Sbirka')
#Pr. 9.1
source('AS-funkce.R')                                            #nacteni skriptu s funkce Scheffe -- potreba stahnout novou verzi skriptu ze studijnich materialu
data <- read.delim('cas_delniku.txt', sep='', dec=',', header=T) #nacteni dat
head(data)
X  <- data$X
ID <- data$ID
nazvy <- c('delnik1', 'delnik2', 'delnik3')

#Pruzkumova analyza
prum <- sm <- N <- NULL                   # vytvoreni prazdnych vektoru, do kterych budeme nasledne doplnovat hodnoty pomoci cyklu
for(i in 1:3){
  prum[i] <- mean(X[ID==i])               # do vektoru prum na prvni pozici vloz prumer z hodnot X, pro ktere bylo ID=1 (prumer casu prvniho delnika), nasledne totez pro druheho a tretiho delnika
  N[i]    <- length(X[ID==i])             # do vektoru N na prvni pozici vloz pocet hodnot X, pro ktere bylo ID=1 (pocet namerenych casu prvniho delnika), nasledne totez pro druheho a tretiho delnika
  sm[i]   <- sd(X[ID==i])                 # do vektoru sm na prvni pozici vloz sm.odch. z hodnot X, pro ktere bylo ID=1 (sm.odchylka casu prvniho delnika), nasledne totez pro druheho a tretiho delnika
}
prum[4] <- mean(X)                        # do vektoru prum na ctvrou pozici vloz prumer vsech namerenych casu ukonu
N[4] <- length(X)                         # do vektoru N na ctvrou pozivi vloz pocet vsech namerenych casu
sm[4] <- sd(X)                            # do vektoru sm na ctvrtou pozici vloz sm. odchylku vsech namerenych casu ukonu

tab <- data.frame(prumer= prum, N=N, sm.odch=sm,
                  row.names=c(nazvy, 'vsichni')) #prehledna tabulka vysledku
round(tab, 4)

# Krabicovy graf
boxplot(X[ID==1], X[ID==2], X[ID==3],
        names=nazvy, ylab='cas (v min)',
        main='Cas k provedeni ukonu',
        col='khaki1', border='orange4')
mtext('Krabicovy graf', line=0.4, cex=0.8)
points(prum[1:3], pch=15, col='brown4')

# Testy normality casu zvlast pro kazdeho delnika
for(i in 1:3){
  print(shapiro.test(X[ID==i])$p.value)
}
# Pro vsechny tri vybery byla p-hodnota > 0.05 -> H0 o normalite vsech tri vyberu nezamitame na hl. vyzn. 0.05.

#Test homogenity dat
# H0:sigma1^2=sigma2^2=sigma3^2
library(lawstat)                         # knihovna obsahujici funkci levenee.test
levene.test(X, ID, location='mean')
# P-hodnota Levenova testu je vetsi nez 0.05, tedy H0 o shode rozptylu vsech tri vyberu nezamitame na hladine vyznamnosti 0.05.


# ANOVA - analyza rozptylu jednoducheho trideni
# H0:mu_1=mu_2=mu_3
#H1: alespon jedna dvojice mu_i a mu_j se vyznamne lisi.

r <- length(unique(ID)) # pocet skupin
n <- length(X)          # celkovy pocet pozorovani

Xi. <- ni <- NULL       # priprava prazdnych vektoru
for(i in 1:3){
  Xi.[i] <- sum(X[ID==i])   # do vektoru Xi. na prvni pozici vloz soucet namerenych casu prvniho delnika; totez proved pro druheho a tretiho delnika
  ni[i]  <- length(X[ID==i]) # do vektori ni na prvni pozivi vloz pocet namerenych casu prvniho delnika; totez proved pro druheho a tretiho delnika
}
Mi. <- Xi./ni                # aritm. prumery namerenych casu pro kazdeho delnika
X.. <- sum(X)                # soucet vsech nymerenych hodnot
M.. <- mean(X)               # aritmeticky prumer vsech namerenych hodnot

SA  <- sum(ni*(Mi.-M..)^2)   # skupinovy soucet ctvercu
fA  <- r-1                   # stupne volnosti skupinoveho souctu ctvercu
ST  <- sum((X-M..)^2)        # celkovy soucet ctvercu
fT  <- n-1                   # stupne volnosti celkoveho souctu ctvercu
SE  <- ST-SA                 # rezidualni soucet ctvercu
fE  <- n-r                   # stupne volnosti rezidualniho souctu ctvercu
FA  <- (SA/fA)/(SE/fE)       # testovaci statistika pro ANOVU

alpha <- 0.05                # hladina vyznamnosti alpha
qf(1-alpha, r-1, n-r)        # dolni hranice kritickeho oboru (u ANOVY ma kriticky obor vzdy tvar (dh, nekonecno))

p.val <- 1-pf(FA, r-1, n-r)  # p-hodnota pro ANOVU

# Protoze FA nalezi do kritickeho oboru, H0 o shode strednich hdonot mu_1=mu_2=mu_3 zamitame na hladine vyznamnosti 0.05.
# Protoze p-hodnota < 0.05, H0 o shode strednch hodnot zamitame na hl. vyznamnosti 0.05.

# Prehledna tabulka vsech mezivysledku
tab <- data.frame(SA=SA, fA=fA, SE=SE, fE=fE, ST=ST, fT=fT,
                  FA=FA, p.val=p.val)

round(tab,4)

# H0 o shode strednich hodnot jsme zamitli -> Scheffeho metoda mnohonasobneho porovnavani
# Zajima nas, ktere dvojice strednich hodnot se od sebe vyznamne lisi
# H01: mu_1=mu_2
# H02: mu_2=mu_3
# H03: mu_1=mu_3
Scheffe(X, ID, names=nazvy, alpha=0.05)
# Porovnavame pravou a levou stranu nerovnosti, kde funkce Scheffe() nam vyhodi vysledky pravych a levych stran pro vsechny hypotezy najednou
# H0i zamitame, pokud je leva strana vetsi nez prava strana:
# 1. a 2. delnik: L=0.583333 > R=0.4282131 -> H01   zamitame na hl. vyzn. 0.05
# 2. a 3. delnik: L=0.05     < R=0.3830054 -> H02 nezamitame na hl. vyzn. 0.05
# 1. a 3. delnik: L=0.633333 > R=0.4282131 -> H03   zamitame na hl. vyzn. 0.05
# -> Zaver: na hl. vyzn. 0.05 jsme zamitli nulovou hypotezu o shode strednich hodnot casu vykonu prvniho a tretiho delnika a prvniho a druheho delnika
# -> prvni delnik pracuje vyrazne odlisne od zbylych dvou delniku.
# Z krab. grafu vidime, ze casy ukonu prvniho delnika jsou mensi, nez casy zbylych dvou, tedy prvni delnik pracuje vyznamne rychleji nez zbyli dva delnici
# Zaver analyzy pro vedouciho podniku je, ze prvnimu denikovi by mel zvysit plat, nebo udelit premii :).