#####################################################################################
################ BIOSTATISTIKA PRO MATEMATICKOU BIOLOGII - CVICENI 3 ################
#####################################################################################


############### Nacteni dat ############### 

kardio <- read.table(file="Biostatistika_data_kardio.txt", sep="\t", dec=".", header=T) # nacteni dat



############### Popisna statistika spojitych promennych ###############

# popisna statistika jedne promenne:
summary(kardio$vek)
summary(kardio$vek_chybny)
summary(kardio$vek[300:350])
summary(kardio$vek_chybny[300:350])

# smerodatna odchylka:
sd(kardio$vek)
sd(kardio$vek_chybny)
sd(kardio$vek[300:350])
sd(kardio$vek_chybny[300:350])

# popisna statistika spojite promenne podle kategorialni promenne:
summary(kardio$vek[kardio$pohlavi=="M"])
summary(kardio$vek[kardio$pohlavi=="Z"])

# vytvareni hezcich tabulek:
library(Epi)
stat.table(pohlavi_rek,contents=list(count(),mean(vek),median(vek),min(vek),max(vek)),data=kardio) # N, prumer, median, minimum a maximum veku podle pohlavi

# uprava popisku sloupcu, aby tabulka vypadala lepe:
stat.table(index=list("Pohlavi"=pohlavi_rek),contents=list(N=count(),"Prumer"=mean(vek),"Median"=median(vek),"Minimum"=min(vek),"Maximum"=max(vek)),data=kardio)

# ulozeni tabulky do promenne, prevod na typ "data frame" a ulozeni do textoveho souboru:
tab1 <- stat.table(index=list("Pohlavi"=pohlavi_rek),contents=list(N=count(),"Prumer"=mean(vek),"Median"=median(vek),"Minimum"=min(vek),"Maximum"=max(vek)),data=kardio) # ulozeni do promenne tab1
tab1 <- as.data.frame(t(tab1))    # prevod na typ "data frame"; je nutne tab1 transponovat pomoci t(tab1)
tab1[,2]<-round(tab1[,2],1)       # zaokrouhleni prumeru na jedno desetinne misto
write.table(tab1,file="Tab1.txt",quote=FALSE, sep="\t",eol="\n", na="NA",dec=".",row.names=TRUE) # ulozeni tabulky do textoveho souboru

# dalsi funkce slouzici k sumarizaci dat: 
# "aggregate" z knihovny "stats" - defaultni knihovna, nemusi se instalovat
# "cast" z knihovny "reshape" - knihovnu "reshape" je treba nainstalovat a nacist:
install.packages("reshape")       # instalace pomoci prikazove radky (je tez mozno knihovnu nainstalovat kliknutim na: Packages -> Install package(s) -> zvoleni statu, napr. Austria -> vyber knihovny -> OK)
library(reshape) # nacteni knihovny "reshape" 

# seznam funkci pro popisnou statistiku:
x <- kardio$vek
length(x)                         # pocet prvku promenne x
mean(x)                           # prumer
median(x)                         # median
min(x)                            # minimum
max(x)                            # maximum
range(x)                          # rozsah hodnot promenne x
c(min(x),max(x))                  # stejny vysledek jako range(x)
var(x)                            # rozptyl
sd(x)                             # smerodatna odchylka
quantile(x,0.25)                  # 1. kvartil (tzn. 25. percentil)
quantile(x,0.5)                   # 2. kvartil (tzn. 50. percentil) - totez co median(x) 
quantile(x,0.75)                  # 3. kvartil (tzn. 75. percentil)
sum(x)                            # suma hodnot x
prod(x)                           # soucin hodnot x

x <- kardio$vaha
mean(x)
x <- kardio$vaha[!is.na(kardio$vaha)]
mean(x)

## Ukol 1. Vytvorte tabulku, ktera bude venovana srovnani celkoveho cholesterolu podle koureni.
#          Tabulka bude obsahovat pocet validnich zaznamu, prumer, median, 5. percentil, 95. percentil a smerodatnou odchylku.



rm(list=ls())


#####################################################################################

############### Vyberovy prumer a stredni chyba prumeru (SE) ###############

kardio <- read.table(file="Biostatistika_data_kardio.txt", sep="\t", dec=".", header=T) # nacteni dat

# nahodny vyber z realnych dat:
# Predpokladejme, ze nas datovy soubor je cilova populace, z niz budeme provadet nahodny vyber.

# vybereme 10 lidi a spocitame jejich prumerny vek, smerodatnou odchylku a stredni chybu prumeru
vyb <- sample(kardio$vek)[1:10]   # vyber 10 lidi
mean(vyb)                         # prumer
sd(vyb)                           # smerodatna odchylka
sd(vyb)/length(vyb)               # stedni chyba prumeru

# vybereme jinych 10 lidi a spocitame jejich prumerny vek, smerodatnou odchylku a stredni chybu prumeru
vyb <- sample(kardio$vek)[1:10]   # vyber 10 lidi
mean(vyb)                         # prumer
sd(vyb)                           # smerodatna odchylka
sd(vyb)/length(vyb)               # stedni chyba prumeru

# 100x vybereme 10 lidi a spocitame jejich prumerny vek:
windows()                         # otevreni noveho grafickeho okna
par(mfrow=c(1,2))                 # nastaveni vykresleni 2 grafu do 1 okna
x<-NULL
for (i in 1:100) {                
vyb <- sample(kardio$vek)[1:10]
x[i] <- mean(vyb)}                # cyklus pro výpočet výběrových průměrů pro n=10
hist(x, breaks=10, xlim=c(25,60)) # vykreslení histogramu pro výběrové průměry pro n=10

# vybereme 100 lidi a spocitame jejich prumerny vek, smerodatnou odchylku a stredni chybu prumeru
vyb <- sample(kardio$vek)[1:100]  # vyber 100 lidi
mean(vyb)                         # prumer
sd(vyb)                           # smerodatna odchylka
sd(vyb)/length(vyb)               # stedni chyba prumeru

# vybereme jinych 100 lidi a spocitame jejich prumerny vek, smerodatnou odchylku a stredni chybu prumeru
vyb <- sample(kardio$vek)[1:100]  # vyber 100 lidi
mean(vyb)                         # prumer
sd(vyb)                           # smerodatna odchylka
sd(vyb)/length(vyb)               # stedni chyba prumeru

# 100x vybereme 100 lidi a spocitame jejich prumerny vek:
for (i in 1:100) {                
vyb <- sample(kardio$vek)[1:100]
x[i] <- mean(vyb)}                # cyklus pro výpočet výběrových průměrů pro n=10
hist(x, breaks=10, xlim=c(25,60)) # vykreslení histogramu pro výběrové průměry pro n=10



### ukazka pomoci simulovanych dat:
windows()                         # otevreni noveho grafickeho okna
par(mfrow=c(1,3))                 # nastaveni vykresleni 3 grafu do 1 okna

x <- rep(0, 100)                  # vytvořím si vektor pro ukládání průměrů

# Vyberovy prumer pro n=10:
pom <- rnorm(10, 175, 15)
mean(pom)
sd(pom)
se <- sd(pom)/length(pom)
se
for (i in 1:100) {
pom <- rnorm(10, 175, 15)         # nahodne generovani 10 hodnot se stredni hodnotou 175 a smerodatnou odchylkou 15
x[i] <- mean(pom)}                # cyklus pro výpočet výběrových průměrů pro n=10
hist(x, breaks=10, xlim=c(160,190)) # vykreslení histogramu pro výběrové průměry pro n=10

# Vyberovy prumer pro n=100:
pom <- rnorm(100, 175, 15)
mean(pom)
sd(pom)
se <- sd(pom)/length(pom)
se
for (i in 1:100) {
pom <- rnorm(100, 175, 15)        # nahodne generovani 100 hodnot se stredni hodnotou 175 a smerodatnou odchylkou 15
x[i] <- mean(pom)}                # cyklus pro výpočet výběrových průměrů pro n=100
hist(x, breaks=10, xlim=c(160,190)) # vykreslení histogramu pro výběrové průměry pro n=100

# Vyberovy prumer pro n=1000:
pom <- rnorm(1000, 175, 15)
mean(pom)
sd(pom)
se <- sd(pom)/length(pom)
se
for (i in 1:100) {
pom <- rnorm(1000, 175, 15)       # nahodne generovani 1000 hodnot se stredni hodnotou 175 a smerodatnou odchylkou 15
x[i] <- mean(pom)}                # cyklus pro výpočet výběrových průměrů pro n=1000
hist(x, breaks=10, xlim=c(160,190)) # vykreslení histogramu pro výběrové průměry pro n=1000



#####################################################################################

############### Intervaly spolehlivosti (IS) ###############

##### IS pro vyberovy prumer pri neznamem sigma:

x <- kardio$cel_cholesterol
windows()
hist(x)
prumer <- mean(x)                 # vyberovy prumer
s <- sd(x)                        # vyberova smerodatna odchylka
n <- length(x)                    # velikost vzorku
alpha <- 0.1                      # hladina spolehlivosti - tzn. 90% IS
kvantil <-  qt(1-alpha/2,n-1)     # kvantil Studentova rozdeleni
IS.dolni <- prumer - (s/sqrt(n))*kvantil
IS.horni <- prumer + (s/sqrt(n))*kvantil
IS.dolni
IS.horni
windows()                         # otevreni noveho grafickeho okna
plot(c(IS.dolni,IS.horni),c(0,0),"l",xlim=c(4.8,5.2),ylim=c(-2,0)) # vykresleni IS
points(prumer,0)                  # vykresleni polohy prumeru

# ovlivneni sirky IS spolehlivosti, kterou pozadujeme:
alpha <- 0.05                     # hladina spolehlivosti - tzn. 95% IS
kvantil <-  qt(1-alpha/2,n-1)     # kvantil Studentova rozdeleni
IS.dolni <- prumer - (s/sqrt(n))*kvantil
IS.horni <- prumer + (s/sqrt(n))*kvantil
lines(c(IS.dolni,IS.horni),c(-1,-1))
points(prumer,-1)

# ovlivneni sirky IS velikosti vzorku:
x <- kardio$cel_cholesterol[1:300] # vyberu prvnich 300 lidi
prumer <- mean(x)                 # vyberovy prumer
s <- sd(x)                        # vyberova smerodatna odchylka
n <- length(x)                    # velikost vzorku
alpha <- 0.05                     # hladina spolehlivosti
kvantil <-  qt(1-alpha/2,n-1)     # kvantil Studentova rozdeleni
IS.dolni <- prumer - (s/sqrt(n))*kvantil
IS.horni <- prumer + (s/sqrt(n))*kvantil
lines(c(IS.dolni,IS.horni),c(-2,-2))
points(prumer,-2)

# ovlivneni sirky IS variabilitou nahodne veliciny:
var1 <- rnorm(500, 175, 15)
prumer <- mean(var1)                 # vyberovy prumer
s <- sd(var1)                        # vyberova smerodatna odchylka
n <- length(var1)                    # velikost vzorku
alpha <- 0.1                         # hladina spolehlivosti - tzn. 90% IS
kvantil <-  qt(1-alpha/2,n-1)        # kvantil Studentova rozdeleni
IS.dolni <- prumer - (s/sqrt(n))*kvantil
IS.horni <- prumer + (s/sqrt(n))*kvantil
windows()                            # otevreni noveho grafickeho okna
plot(c(IS.dolni,IS.horni),c(0,0),"l",xlim=c(170,180),ylim=c(-1,0)) # vykresleni IS
points(prumer,0)                     # vykresleni polohy prumeru

var2 <- rnorm(500, 175, 45)
prumer <- mean(var2)                 # vyberovy prumer
s <- sd(var2)                        # vyberova smerodatna odchylka
n <- length(var2)                    # velikost vzorku
alpha <- 0.1                         # hladina spolehlivosti - tzn. 90% IS
kvantil <-  qt(1-alpha/2,n-1)        # kvantil Studentova rozdeleni
IS.dolni <- prumer - (s/sqrt(n))*kvantil
IS.horni <- prumer + (s/sqrt(n))*kvantil
lines(c(IS.dolni,IS.horni),c(-1,-1))
points(prumer,-1)



##### IS pro vyberovy prumer pri znamem sigma:
var <- rnorm(500,175,15)
# predpokladame, ze sigma zname - sigma je 15:
sigma <- 15
prumer <- mean(var)
n <- length(var)
alpha <- 0.05
kvantil <- qnorm(1-alpha/2)
IS.dolni <- prumer - (sigma/sqrt(n))*kvantil
IS.horni <- prumer + (sigma/sqrt(n))*kvantil
IS.dolni
prumer
IS.horni



##### IS pro podil rozptylu sigma1^2 a sigma2^2:

# realny priklad:
x <- kardio$cel_cholesterol
install.packages("plotrix")
library(plotrix)
multhist(list(x[kardio$pohlavi_rek=="Z"],x[kardio$pohlavi_rek=="M"]))
x1 <- x[kardio$pohlavi_rek=="Z"]  # celkovy cholesterol u zen
n1 <- length(x1)
s1 <- sd(x1)
x2 <- x[kardio$pohlavi_rek=="M"]  # celkovy cholesterol u muzu
n2 <- length(x2)
s2 <- sd(x2)
alpha <- 0.05
IS.dolni <- (s1^2/s2^2)/qf(1-alpha/2,n1-1,n2-1) 
IS.horni <- (s1^2/s2^2)/qf(alpha/2,n1-1,n2-1) 
IS.dolni
IS.horni
windows()                         # otevreni noveho grafickeho okna
plot(c(IS.dolni,IS.horni),c(0,0),"l") # vykresleni IS
points(s1^2/s2^2,0)               # vykresleni polohy prumeru

# simulace:
x1 <- rnorm(100,185,10)
x2 <- rnorm(120,145,40)
multhist(list(x1,x2))
n1 <- length(x1)
n2 <- length(x2)
s1 <- sd(x1)
s2 <- sd(x2)
alpha <- 0.05
IS.dolni <- (s1^2/s2^2)/qf(1-alpha/2,n1-1,n2-1) 
IS.horni <- (s1^2/s2^2)/qf(alpha/2,n1-1,n2-1) 
IS.dolni
s1^2/s2^2
IS.horni



##### IS pro odhad stredni hodnoty rozdilu prumeru:

# realny priklad:
x <- kardio$vek
multhist(list(x[kardio$pohlavi_rek=="Z"],x[kardio$pohlavi_rek=="M"]))
x1 <- x[kardio$pohlavi_rek=="Z"]  # vek u zen
prumer1 <- mean(x1)
n1 <- length(x1)
s1 <- sd(x1)
x2 <- x[kardio$pohlavi_rek=="M"]  # vek u muzu
prumer2 <- mean(x2)
n2 <- length(x2)
s2 <- sd(x2)
alpha <- 0.05
kvantil <- qt(1-alpha/2,n1+n2-2)  # kvantil Studentova rozdeleni
IS.dolni <- prumer1-prumer2-sqrt(s1^2/n1+s2^2/n2)*kvantil
IS.horni <- prumer1-prumer2+sqrt(s1^2/n1+s2^2/n2)*kvantil
IS.dolni
IS.horni
windows()                         # otevreni noveho grafickeho okna
plot(c(IS.dolni,IS.horni),c(0,0),"l") # vykresleni IS
points(prumer1-prumer2,0)         # vykresleni polohy prumeru

# simulace:
x1 <- rnorm(100,185,15)
x2 <- rnorm(120,145,15)
multhist(list(x1,x2))
prumer1 <- mean(x1)
n1 <- length(x1)
s1 <- sd(x1)
prumer2 <- mean(x2)
n2 <- length(x2)
s2 <- sd(x2)
alpha <- 0.05
kvantil <- qt(1-alpha/2,n1+n2-2)  # kvantil Studentova rozdeleni
IS.dolni <- prumer1-prumer2-sqrt(s1^2/n1+s2^2/n2)*kvantil
IS.horni <- prumer1-prumer2+sqrt(s1^2/n1+s2^2/n2)*kvantil
IS.dolni
prumer1-prumer2
IS.horni



##### IS pro velicinu s binomickym rozdelenim:

# Vypoctete IS pro relativni cetnost pacientek s rakovinou prsu v danem souboru 1000 zen, kdyz vime, ze pacientek s rakovinou prsu bylo 120.
r <- 120
n <- 1000
pi <- r/n
alpha <- 0.05
# aby bylo mozne pouzit kvantil, standardizovaneho normalniho rozdeleni, musi byt splneno: n>30, r>5, (n-r)>5
kvantil <- qnorm(1-alpha/2)
IS.dolni <- pi - kvantil*sqrt(pi*(1-pi)/n) 
IS.horni <- pi + kvantil*sqrt(pi*(1-pi)/n) 
IS.dolni
pi
IS.horni
