#
## KONTINGENCNI TABULKY

# data: pastelky velke (150 pastelek), planovani

# JEDEN VYBER, JEDEN FAKTOR: krevni skupiny
# mam nahodny vyber, pozorovani jsou roztrizena podle 1 faktoru (zde krevni skupina)
vyber<-c(25,7,1,17) # cetnosti pro krevni skupinu A, B, AB, 0.
psti<-c(0.43, 0.19, 0.10, 0.28)  # pravdepodobnosti, ktere chci testovat
# H0: muj nahodny vyber odpovída temto pravdepodobnostem 
# presneji: cetnosti ocekavane na zaklade pravdepodobnosti jsou srovnatelne 
# s cetnostmi pozorovanymi
# chi-kvadrat test dobre shody:
chisq.test(vyber, p=psti)   # vektor pravdepodobnosti musim oznacit jako "p=..."
                            # jinak program predpoklada, ze je to druhy vyber "y"
# test dobre shody ma predpoklad: ocekavane cetnosti >= 5
chisq.test(vyber, p=psti)$expected  # zde jsou nachystane
# vysledek: p=0.17 ... nezamitam hypotezu, ze pozorovane cetnosti odpovidaji
                     # danym pravdepodobnostem


# VYTVORIT KONTINGENCNI TABULKU z datove tabulky
load("data_planovani.RData")
summary(planovani)  # vsimni si, jak jsou data zapsana: co jedinec, to radek
                    # tvar sumarni tabulky musim teprve vytvorit
xtabs(~plan+vzdelani, data=planovani)   # funkce xtabs tvori tabulku
ktplan<-xtabs(~plan+vzdelani, data=planovani)  # uloz tabulku do objektu "ktplan"
addmargins(ktplan)    # k tabulce dopln marginalni soucty (jen pro tisk)


# DVA FAKTORY: NEZAVISLOST?
# tady chi-squared test potrebuje dostat tu tabulku:
chisq.test(ktplan)
# kontrola nejmensi OCEKAVANE cetnosti:
chisq.test(ktplan)$expected    # moc desetinnych mist, zaokrouhlime to pro lepsi cteni
round(chisq.test(ktplan)$expected,1)  # zaokrouhli vysledky na 1 desetinne misto (digits)
  # jake dalsi polozky obsahuje objekt vysledku testu?
  names(chisq.test(ktplan))
  # [1] "statistic" "parameter" "p.value"   "method"    "data.name"
  # [6] "observed"  "expected"  "residuals" "stdres"  


# DVA VYBERY: STEJNA STRUKTURA?
xtabs(~porodnice+vzdelani, data=planovani)  # jak vypadá tabulka?
ktporod<-xtabs(~porodnice+vzdelani, data=planovani) # ulož tabulku do objektu
chisq.test(ktporod) 
round(chisq.test(ktporod)$expected, 1) # zaokrouhli vysledky na 1 desetinné místo (digits)

# CTYRPOLNI TABULKA
# data hrabos
kthrabos <- xtabs(pocty~ Frenkelia+Sarcocystis, data=Hrabos)
chisq.test(kthrabos, correct = F) # Yatesova oprava vypnuta
chisq.test(kthrabos)$expected     # jenže podmínka o_i>5 nesplněna
chisq.test(kthrabos, correct = T) # Yatesova oprava zapnuta

# mira tesnosti: phi koeficient a Cramerovo V
# vlastni funkce na vypocet Cramerova V pro 2x2 tabulky
cramer<-function(tab){a<-tab[1,1];b<-tab[1,2];c<-tab[2,1];d<-tab[2,2];
(a*d-b*c)/sqrt((a+b)*(c+d)*(a+c)*(b+d))}
cramer(kthrabos)

# vypocet phi koeficientu podle vzorecku
# potrebuju chi-kvadrat statistiku (bez Yatesovy opravy)
chisq.test(kthrabos, correct=F)$stat
# toto cislo podelit n=sum(kthrabos) a cele odmocnit sqrt()
phi<-function(tab){sqrt(chisq.test(tab, correct=F)$stat/sum(tab))} # phi koeficient
phi(kthrabos) # warovani se opet tyka nizkych ocekavanych cetnosti (<5)

# vypocty schovane v baliku vcd
install.packages("vcd")  # stahni balik vcd
library(vcd)  # otevri knihovnu novych prikazu
assocstats(kthrabos)  # pouzij funkci assocstats
  # bohuzel Cramerovo V je pocitano bez znamenka, vzdy kladne, to je chyba/skoda.

# modifikovany hrabos: jen parazitovani jedinci
hrabos2<-matrix(c(4,27,11,0),nrow = 2,ncol=2, byrow = T)
hrabos2
cramer(hrabos2) # tady uz je jasna NEGATIVNI asociace
phi(hrabos2)  # phi koeficient nenese informaci o pozitivite nebo negativite
assocstats(hrabos2)

# Fisheruv presny test
  # nulove hypotezy zustavaji stejne: Jsou vlastnosti nezavisle? 
  #                              nebo Maji vybery stejnou pravdepodobnostni strukturu?
fisher.test(kthrabos)
fisher.test(hrabos2)
  # hypoteza se presneji tyka pomeru sanci..., ale interpretace se nemeni

# McNEMARUV TEST SYMETRIE (BOWKERŮV TEST)
ktles<-xtabs(pocty~r1994+r1995, data=Les)
addmargins(ktles)
# nulova hypoteza: tabulka (pravdepodobnosti) je symetricka;
  # pro nas pochopitelneji: (relativni) cetnosti stupnu nakazy jsou srovnatelne
  # v prvnim i ve druhem mereni. Spravne se vyrok tyka marginalnich pravdepodobnosti.
mcnemar.test(ktles)