#######################################################################################
###################### Vybrané neparametrické testy ###################################
######################### Martina Litschmannová #######################################
#######################################################################################

################## Přiklad 1. (Kostka) ################################################
#######################################################################################
#H0: Kostka je férová.
#Ha: Kostka není férová.

x = c(1,2,3,4,5,6)
n.obs = c(979,1002,1015,980,1040,984)
p.exp = c(1/6,1/6,1/6,1/6,1/6,1/6)
n.exp = 6000*p.exp
n.exp # nutno zkontrolovat předpoklady testu
# Všechny očekávané četnosti jsou větší než 5.

x.obs = sum(((n.obs-n.exp)^2)/n.exp)
x.obs
p.hodnota = 1-pchisq(x.obs,5)
p.hodnota

# Na hladině významnosti 0,05 nezamítáme HO (p-hodnota = 0,711, Chí-kvadrát test nezávislosti,
# df = 5).

################## Přiklad 2. (Úplně specifikovaný test) ##############################
#######################################################################################
#H0: Počet poruch během 100 hodin provozu lze modelovat Poissonovým rozdělením s parametrem 1,2.
#Ha: Počet poruch během 100 hodin provozu nelze modelovat Poissonovým rozdělením s parametrem 1,2.

x = c(0,1,2,3,4)
n.obs = c(52,48,36,10,4)
p.exp = dpois(x,1.2)
n.exp = 150*p.exp
n.exp # nutno zkontrolovat předpoklady testu
# 4 z 5 očekávaných četnosti, tj. 80%, jsou větší než 5.

x.obs = sum(((n.obs-n.exp)^2)/n.exp)
x.obs
p.hodnota = 1-pchisq(x.obs,4)
p.hodnota

# Na hladině významnosti 0,05 nezamítáme HO (p-hodnota = 0,590, Chí-kvadrát test nezávislosti,
# df = 4).

################## Přiklad 3. (Neúplně specifikovaný test) ##############################
#######################################################################################
#H0: Počet poruch během 100 hodin provozu lze modelovat Poissonovým rozdělením.
#Ha: Počet poruch během 100 hodin provozu nelze modelovat Poissonovým rozdělením.

x = c(0,1,2,3,4)
n.obs = c(52,48,36,10,4)

lambda.t = weighted.mean(x,n.obs)   # odhad parametru Poissonova rozdělení
p.exp = dpois(x,lambda.t)
n.exp = 150*p.exp
n.exp # nutno zkontrolovat předpoklady testu
# 4 z 5 očekávaných četnosti, tj. 80%, jsou větší než 5.

x.obs = sum(((n.obs-n.exp)^2)/n.exp)
x.obs
p.hodnota = 1-pchisq(x.obs,3)
p.hodnota

# Na hladině významnosti 0,05 nezamítáme HO (p-hodnota = 0,491, Chí-kvadrát test nezávislosti,
# df = 3).

################## Přiklad 4. (Dálnice) ###############################################
#######################################################################################

## Načtení dat z xlsx souboru (pomoci balíčku openxlsx)
library(openxlsx)
dalnice = readWorkbook("C:/Users/lit40/Desktop/neparametricke_hypotezy.xlsx",
                    sheet=2,
                    colNames=TRUE,
                    startRow = 1)   
colnames(dalnice)="hodnoty"

#H0: Rozestupy mezi vozidly lze modelovat normálním rozdělením.
#Ha: Rozestupy mezi vozidly nelze modelovat normálním rozdělením.

install.packages("nortest")
library(nortest)
pearson.test(dalnice$hodnoty)

# Určení počtu stupňů volbnosti
pom = pearson.test(dalnice$hodnoty)
attributes(pom)
pom$df

#Na hladině významnosti 0,05 lze zamítnout HO (p-hodnota << 0,001, Chí-kvadrát test dobré shody,
# df = 12).

################## Přiklad 5. (Tetování) ##############################################
#######################################################################################

## Načtení dat z xlsx souboru (pomoci balíčku openxlsx)
library(openxlsx)
tet = readWorkbook("C:/Users/lit40/Desktop/neparametricke_hypotezy.xlsx",
                       sheet=3,
                       colNames=TRUE,
                       startRow = 1)   
tet = tet[,c(6,10)]
colnames(tet)=c("pohlavi","tetovani")

#H0: Mezi pohlavím a pořízením si tetování neexistuje souvislost.
#Ha: Mezi pohlavím a pořízením si tetování existuje souvislost.

data = table(tet$pohlavi,tet$tetovani)
data

mosaicplot(data,
           color = topo.colors(2),
           main = "")

library(lsr)
cramersV(data)

pom = chisq.test(data)
attributes(pom)
pom$expected
# Všechny očekávané četnosti jsou větší než 5.
pom$p.value

#Na hladině významnosti 0,05 lze zamítnout HO (p-hodnota = 0,003, Chí-kvadrát test dobré shody,
# df = 1). Pozorovanou závislost lze hodnotit jako slabou (Cramerovo V = 0,121).

################## Přiklad 6. (Spokojenost) ###########################################
#######################################################################################

data = matrix(c(10,25,50,15,20,10,130,40),nrow=2,byrow=T)
data
rownames(data) = c("Praha","Venkov")
colnames(data) = c("v. nespokojen","s. nespokojen","s. spokojen","v. spokojen")
data

### Explorační analýza
prop.table(data)  # sdružené relativní četnosti
prop.table(data,1)  # řádkové relativní četnosti
prop.table(data,2)  # sloupcové relativní četnosti

# shlukový sloupcový graf
barplot(data,beside = T, legend = rownames(data)) 
barplot(t(data),beside = T, legend = colnames(data),
        col = terrain.colors(4)) 

# skládaný sloupcový graf
barplot(data,beside = F, legend = rownames(data)) 
barplot(t(data),beside = F, legend = colnames(data),
        col = terrain.colors(4))

# mozaikový graf
mosaicplot(data)
mosaicplot(t(data))

par(mfrow = c(1,1),mar = c(2,5,2,2))
mosaicplot(data,
           col = terrain.colors(4),
           main = "",
           las = 1  # změna orientace popisků osy y
           ) 

mosaicplot(data,
           col = terrain.colors(4),
           main = "",
           las = 1,  # změna orientace popisků osy y
           dir = c("h","v") # otočení mozaikového grafu o 90°
           )  

# Cramerovo V
# install.packages("lsr")
library(lsr)
cramersV(data)

#H0: Mezi spokojenosti v práci a umístěním podniku neexistuje souvislost.
#Ha: Mezi spokojenosti v práci a umístěním podniku existuje souvislost.

pom = chisq.test(data)
attributes(pom)
pom$expected
# Všechny očekávané četnosti jsou větší než 5.
pom$p.value

#Na hladině významnosti 0,05 lze zamítnout HO (p-hodnota << 0,001, Chí-kvadrát test dobré shody,
# df = 3). Pozorovanou závislost lze hodnotit jako středně silnou (Cramerovo V = 0,296).

################## Přiklad 7. (Apoplexie) #############################################
#######################################################################################

data = matrix(c(171,3264,117,4320),nrow=2,byrow=T)
data
rownames(data) = c("kuřák","nekuřák")
colnames(data) = c("apoplexie","bez apoplexie")
data

#H0: Mezi kouřením a výskytem apoplexie neexistuje souvislost.
#Ha: Mezi kouřením a výskytem apoplexie existuje souvislost.

mosaicplot(data,
           color = topo.colors(2),
           main = "")

library(lsr)
cramersV(data)
# Dle mozaikového grafu a Cramerova V (0,061) lze souvislost mezi kuřáctvím a výskytem apoplexie
# hodnotit jako velmi slabou.

# adc) riziko = pravděpodobnost
prop.test(171,3435)
# U kuřáku je riziko vzniku apoplexie cca 5,0 %. Se spolehlivostí 95 % lze očekávat, že
# riziko apoplexie bude u kuřáku mezi 4,3 % a 5,8 % (Clopperův - Pearsonův odhad).
prop.test(117,4437)
# U nekuřáku je riziko vzniku apoplexie cca 2,6 %. Se spolehlivostí 95 % lze očekávat, že
# riziko apoplexie bude u nekuřáku mezi 2,2 % a 3,2 % (Clopperův - Pearsonův odhad).

# add) 
install.packages("epiR")
library(epiR)
epi.2by2(data)
# U kuřáků je cca 1,89 krát vyšší riziko apoplexie než u nekuřáků. S 95% spolehlivostí
# lze očekávat, že toto riziko bude 1,50 až 2,38 krát vyšší.

# ade)
# U kuřáku je šance vzniku apoplexie cca 52:1000. Tj. u 1052 kuřáků lze očekávat cca 52 výskytů apoplexie.
# U nekuřáku je šance vzniku apoplexie cca 27:1000. Tj. u 1027 nekuřáků lze očekávat cca 27 výskytů apoplexie.

# adf)
# U kuřáků je cca 1,93 krát vyšší šance apoplexie než u nekuřáků. S 95% spolehlivostí
# lze očekávat, že tato šance bude 1,52 až 2,46 krát vyšší.

# adg)
# Pozor! Příkaz epi.2by2 nemá jako výstup očekávané četnosti pro Chí-kvadrát test nezávislosti.
# Není tak možno ověřit předpoklady testu!

pom = chisq.test(data)
attributes(pom)
pom$expected
# Všechny očekávané četnosti jsou větší než 5.
pom$p.value

#Na hladině významnosti 0,05 lze zamítnout HO (p-hodnota << 0,001, Chí-kvadrát test dobré shody,
# df = 1). Pozorovanou závislost lze hodnotit jako velmi slabou (Cramerovo V = 0,061).