#######################################################################################
########################### Dvouvýběrové testy ########################################
######################### Martina Litschmannová #######################################
#######################################################################################

################## Přiklad 1. (Cholesterol) ###########################################
#######################################################################################

# Test o shodě středních hodnot / mediánů

## Načtení dat z xlsx souboru (pomoci balíčku readxl)
library(readxl)
chol = read_excel("C:/Martina/Výuka/Pravděpodobnost a statistika/DATA/testy_dvouvyberove.xlsx",
                  sheet = "cholesterol2",
                  skip = 1)   
colnames(chol)=c("mladsi","starsi")

## Převod do standardního datového formátu
chol.s = stack(chol)
chol.s = na.omit(chol.s)
colnames(chol.s) = c ("hodnoty","skupina")

## Explorační analýza
par(mfrow=c(1,1))
boxplot(chol.s$hodnoty~chol.s$skupina)

# Data obsahují odlehlá pozorování.
pom = boxplot(chol.s$hodnoty~chol.s$skupina,plot =F)
pom
# Odstranění odlehlých pozorování:
chol.s$hodnoty.bez = chol.s$hodnoty
chol.s$hodnoty.bez[chol.s$skupina == "mladsi" & chol.s$hodnoty >=4.927]=NA
chol.s$hodnoty.bez[chol.s$skupina == "starsi" & chol.s$hodnoty <=4.116]=NA

boxplot(chol.s$hodnoty.bez~chol.s$skupina)

# Samostatné proměnné
mladsi.bez = chol.s$hodnoty.bez[chol.s$skupina=="mladsi"]
starsi.bez = chol.s$hodnoty.bez[chol.s$skupina=="starsi"]


## Explorační analýza pro data bez odlehlých pozorování
boxplot(chol.s$hodnoty.bez~chol.s$skupina)

tapply(chol.s$hodnoty.bez,chol.s$skupina,summary,na.rm=TRUE)
tapply(chol.s$hodnoty.bez,chol.s$skupina,sd,na.rm=TRUE)

# Ověření normality
par(mfrow=c(1,2))
skupina = c("mladsi","starsi")
for (i in 1:2){
  qqnorm(chol.s$hodnoty[chol.s$skupina == skupina[i]],
         main = "",xlab = "teoretické \nnormované kvantily",ylab = "výběrové kvantily")
  qqline(chol.s$hodnoty[chol.s$skupina == skupina[i]])}

library(moments)
tapply(chol.s$hodnoty.bez,chol.s$skupina,skewness,na.rm=TRUE)
tapply(chol.s$hodnoty.bez,chol.s$skupina,function (x) moments::kurtosis(x,na.rm=TRUE)-3)
# Šikmost i špičatost odpovídá norm. rozdělení! Pro konečné rozhodnutí o normalitě dat použijeme
# test normality.

# Předpoklad normality ověříme Shapirovovým . Wilkovovým testem.
tapply(chol.s$hodnoty.bez,chol.s$skupina,shapiro.test)
# Na hl. významnosti 0.05 nelze předpoklad normality zamítnout 
# (p-hodnota1=0.464, p-hodnota2=0.940, Shapirův-Wilkův test).

# Ověření shody rozptylů

# H0: sigma.mladsi = sigma.starsi
# Ha: sigma.mladsi <> sigma.starsi

var.test(chol.s$hodnoty.bez~chol.s$skupina,ratio=1,conf.level=0.95)
# Na hl. významnosti 0.05 zamítáme předpoklad o shodě rozptylů (p-hodnota<<0.001,test o shodě rozptylů).

# Chceme-li ověření shody rozptylů řešit raději pomocí intervalového odhadu, 
# je pro snazší interpretaci lepší zajistit, aby byl nalezen odhad pro poměr
# max(rozptyl 1, rozptyl 2)/min(rozptyl 1, rozptyl 2). V našem případě tedy poměr 
# rozptyl.starsi/rozptyl.mladsi.

var.test(starsi.bez,mladsi.bez,conf.level=0.95)
# Rozptyl hodnot cholesterolu ve skupině starších mužů je cca 12x větší než 
# rozptyl hodnot cholesterolu ve skupině mladších mužů. S 95% spolehlivostí lze 
# odhadovat, že rozptyl hodnot cholesterolu ve skupině starších mužů je 8 - 19 krát 
# větší než rozptyl hodnot cholesterolu ve skupině mladších mužů. Pozorovanou neshodu 
# mezi rozptyly lze na hladině významnosti 0,05 označit za statisticky významnou (F test, F = 12,2,
# df_num = 83, df_denom = 96, p-hodnota < 0,001).

# Ověření shody středních hodnot (Aspinové-Welchův test)

# Pro snadnou interpretaci intervalového odhadu je vhodné pořadí výběru stanovit tak,
# aby rozdíl jejich průměrů byl kladný.

# H0: mu.starsi = mu.mladsi
# Ha: mu.starsi > mu.mladsi
t.test(starsi.bez,mladsi.bez,mu=0,alternative="greater",var.equal=FALSE,conf.level=0.95)

# Dle výsledků výběrového šetření očekáváme, že střední obsah cholesterolu v krvi straších mužů
# bude cca o 0,52 mmol/l vyšší než střední obsah cholesterolu u mladších mužů.
# S 95% spolehlivostí lze tento rozdíl očekávat větší než 0,46 mmol/l. 
# Na hladině významnosti 0,05 zamítáme předpoklad o shodě středních hodnot cholesterolu ve skupinách 
# mladších a starších mužů ve prospěch alternativy, že starší muži mají vyšší střední hladinu 
# cholesterolu než muži mladší (t-test, t = 13,1, p-hodnota < 0,001).

################## Přiklad 2. (Deprese) ###############################################
#######################################################################################

# Test o shodě středních hodnot / mediánů

## Načtení dat z xlsx souboru (pomoci balíčku readxl)
library(readxl)
deprese = read_excel("C:/Martina/Výuka/Pravděpodobnost a statistika/DATA/testy_dvouvyberove.xlsx",
                     sheet = "deprese")   
colnames(deprese)=c("endo","neuro")

## Převod do standardního datového formátu
deprese.s = stack(deprese)
deprese.s = na.omit(deprese.s)
colnames(deprese.s) = c ("hodnoty","skupina")

## Explorační analýza
par(mfrow = c(1,1))
boxplot(deprese.s$hodnoty~deprese.s$skupina)

# Data obsahují odlehlá pozorování.
pom = boxplot(deprese.s$hodnoty~deprese.s$skupina,plot = F)
pom
# Odstranění odlehlých pozorování
deprese.s$hodnoty.bez = deprese.s$hodnoty
deprese.s$hodnoty.bez[deprese.s$hodnoty>=1547 & deprese.s$skupina == "endo"]=NA

# Samostatné proměnné
endo.bez = deprese.s$hodnoty.bez[deprese.s$skupina=="endo"]
neuro.bez = deprese.s$hodnoty.bez[deprese.s$skupina=="neuro"]

## Explorační analýza pro data bez odlehlých pozorování
par(mfrow = c(1,1))
boxplot(deprese.s$hodnoty.bez~deprese.s$skupina)
# Odlehlá pozorování viditelná na krabicovém grafu již nejsou odlehlými pozorováními z původního
# datového souboru. Tato pozorování bychom neměli již odstraňovat, spíše bychom měli zvolit 
# robustní metody datové analýzy (metody, které nejsou citlivé na porušení předpokladu normality).

tapply(deprese.s$hodnoty.bez,deprese.s$skupina,summary,na.rm=TRUE)
tapply(deprese.s$hodnoty.bez,deprese.s$skupina,sd,na.rm=TRUE)

# Ověření normality
par(mfrow = c(1,2))
skup = c("endo","neuro")
for (i in 1:2){
  qqnorm(deprese.s$hodnoty[deprese.s$skupina == skup[i]],
         main = "",xlab = "teoretické \nnormované kvantily",ylab = "výběrové kvantily")
  qqline(deprese.s$hodnoty[deprese.s$skupina == skup[i]])
}

library(moments)
tapply(deprese.s$hodnoty.bez,deprese.s$skupina,skewness,na.rm=TRUE)
tapply(deprese.s$hodnoty.bez,deprese.s$skupina,function(x) moments::kurtosis(x,na.rm=TRUE)-3)
# Šikmost i špičatost odpovídá norm. rozdělení. Pro konečné rozhodnutí o normalitě dat použijeme
# test normality.

# Předpoklad normality ověříme Shapirovovým - Wilkovovým testem.
tapply(deprese.s$hodnoty.bez,deprese.s$skupina,shapiro.test)
# Na hl. významnosti 0.05 zamítáme předpoklad normality 
# (Shapirův-Wilkův test, p-hodnota_endo<<0,001, p-hodnota_neuro=0,281).


# Ověření shody mediánů (Mannův - Whitneyův test)
# POZOR! V nulové hypotéze testujeme, zda se rozdělení populací, z nichž
# byly výběry pořízeny liší pouze svou polohou (konkrétně mediánem).

par(mfrow = c(2,1))
hist(endo.bez,xlim = c(0,1400))
hist(neuro.bez,xlim = c(0,1400))

# H0: med.endo = med.neuro
# Ha: med.endo < med.neuro
wilcox.test(endo.bez,neuro.bez,
            mu=0,
            exact=FALSE,
            alternative="less",
            conf.level=0.95)

# Na hladině významnosti 0,05 nelze zamítnout hypotézu o shodě 
# mediánů dob do remise onemocnění pro obě skupiny pacientů (Wilcoxonův test s korekcí na spojitost,
# W = 39, p-hodnota = 0,088).

################## Přiklad 3. (Osmolalita - párová data) ##############################
#######################################################################################

# Načtení dat 

library(readxl)
osmolalita = read_excel("C:/Martina/Výuka/Pravděpodobnost a statistika/DATA/testy_dvouvyberove.xlsx", 
                        sheet = "osmolalita", skip = 1)  
osmolalita = osmolalita[,c(2,3)]
colnames(osmolalita)=c("o8","o11")

# Výpočet nárůstu osmolality
osmolalita$narust = osmolalita$o11 -  osmolalita$o8

## Explorační analýza 
par(mfrow = c(1,1))
boxplot(osmolalita$narust)
# Data obsahují odlehlá pozorování.

# Odstranění odlehlé hodnoty
boxplot(osmolalita$narust,plot = F)
osmolalita$narust.bez = osmolalita$narust
osmolalita$narust.bez[osmolalita$narust==150]=NA

## Explorační analýza pro data bez odlehlých pozorování
par(mfrow = c(1,1))
boxplot(osmolalita$narust.bez, ylab = "nárůst osmolality mezi 8 h a 11 h (mmol/kg)")
abline(h=0,lty=2) # přímka definující nulový efekt

summary(osmolalita$narust.bez,na.rm = T)
sd(osmolalita$narust.bez,na.rm = T)

# Ověření normality
qqnorm(osmolalita$narust.bez)
qqline(osmolalita$narust.bez)
hist(osmolalita$narust.bez)
library(moments)
skewness(osmolalita$narust.bez,na.rm = T)
moments::kurtosis(osmolalita$narust.bez,na.rm = T)-3
# Šikmost i špičatost odpovídají normálnímu rozdělení.

# Předpoklad normality ověříme Shapirovovým - Wilkovovým testem.
shapiro.test(osmolalita$narust.bez)
# Na hl. významnosti 0.05 nelze předpoklad normality zamítnout 
# (Shapirův-Wilkův test, W = 0,966, p-hodnota=0,800).

# Párový t-test
# H0: mu.narust = 0 mm
# Ha: mu.narust > 0 mm

t.test(osmolalita$narust.bez,mu=0,alternative="greater")

# Dle výběrového šetření lze očekávat, že osmolalita moči se
# mezi 8 a 11 hodinou zvýší o cca 18 mmol/kg. S 95% spolehlivostí
# lze očekávat, že dojde k navýšení nejméně o 1 mmol/kg.
# Na hladině významnosti 0,05 lze tento nárůst označit za statisticky
# významný (párový t-test, t = 1,9, df = 14, p-hodnota = 0,042).

################## Přiklad 4. (Porovnání kvality) #####################################
#######################################################################################

n.MM = 200  		# rozsah výběru pro firmu MM
x.MM = 14      	# počet vadných součástek pro firmu MM
p.MM = x.MM/n.MM 	# relativní četnost vadných součástek firmy MM (bodový odhad pravděpodobnosti,
# že firma MM vyrobí vadnou součástku)
p.MM

n.PP = 100  		# rozsah výběru pro firmu PP
x.PP = 10      	# počet vadných součástek pro firmu PP
p.PP = x.PP/n.PP 	# relativní četnost vadných součástek firmy PP (bodový odhad pravděpodobnosti,
# že firma PP vyrobí vadnou součástku)
p.PP

# Ověření předpokladů
n.MM>30 & n.MM > 9/(p.MM*(1-p.MM))
n.PP>30 & n.PP > 9/(p.PP*(1-p.PP))
# Dále pro obě firmy předpokládáme, že  n/N < 0.05, tj. že daná populace (součástek) má rozsah 
# alespoň 20*n, tj. 20*200 (4 000), resp. 20*100 (2 000) součástek, což je asi vcelku reálný předpoklad.

## Waldův test #############################################################
# "Pořadí výběrů" volíme kvůli snazší interpretaci tak, aby rozdíl "p.1 - p.2" vycházel kladný.
## H0: pi.PP = pi.MM
## Ha: pi.PP > pi.MM

x.obs = (p.PP - p.MM)/sqrt((p.PP*(1-p.PP)/n.PP)+(p.MM*(1-p.MM)/n.MM))
x.obs
p.value = 1-pnorm(x.obs)
p.value

# Pravostranný intervalový odhad
alpha = 0.05
p = (x.PP + x.MM)/(n.PP + n.MM)
s.error = qnorm(1-alpha/2)*sqrt(p*(1-p)*(1/n.MM + 1/n.PP))
dolni.mez = (p.PP - p.MM) - s.error
dolni.mez

# Dle výběrového šetření očekáváme, že pravděpodobnost výroby vadné součástky
# je u firmy PP cca o 3 % větší než u firmy MM. Na základě 95% Waldova pravostranného intervalového odhadu
# (-0,035; inf) lze pozorovaný rozdíl v kvalitě výroby označit za statisticky nevýznamný. Ke stejným závěrům
# můžeme dojít i na základě Waldova pravostranného testu (x_obs = 0,857,
# p-hodnota = 0,196).


## Pearsonův X2 test #######################################################
## H0: pi.PP = pi.MM
## Ha: pi.PP > pi.MM

prop.test(c(x.PP,x.MM),c(n.PP,n.MM),alternative="greater",conf.level=0.95)

# Dle výběrového šetření očekáváme, že pravděpodobnost výroby vadné součástky
# je u firmy PP cca o 3 % větší než u firmy MM. Na základě 95% Clopperova - Pearsonova pravostranného intervalového odhadu
# (-0,035; inf) lze pozorovaný rozdíl v kvalitě výroby označit za statisticky nevýznamný. Ke stejným závěrům
# můžeme dojít i na základě Clopperova - Pearsonova pravostranného testu (x_obs = 0,459,
# df = 1, p-hodnota = 0,249).

################################# Pro zájemce #########################
################# Jiná grafická interpretace k příkladu 3 #############
#######################################################################

# Bodový graf doplněný o přímku y = x
plot(osmolalita$o8,osmolalita$o11,
     xlab = "osmolalita v 8 h (mmol/kg)",
     ylab = "osmolalita v 11 h (mmol/kg)")
abline(0,1,col = "red")

# Krabicový graf pro párová data
s = seq(length(osmolalita$o8)) #posloupnost 1,2,...,rozsah dat 
par(mfrow = c(1,1),bty="l") #bty - nastavení ohraničení grafu (jsou viditelné pouze osy x,y, nikoliv rámeček kolem grafu)
boxplot(osmolalita$o8,osmolalita$o11,
        ylab="osmolalita (mmol/kg)",
        names=c("8 hodin","11 hodin"),
        cex.main = 0.9)
stripchart(list(osmolalita$o8,osmolalita$o11),
           vertical=T,
           pch=16,
           cex=0.7,
           add=T) # přidání bodového grafu do boxplotu
#vykreslení spojnic mezi datovými páry
x0 = rep(1,length(osmolalita$o8))
y0 = osmolalita$o8
x1 = rep(2,length(osmolalita$o11))
y1 = osmolalita$o11
segments(x0,y0,x1,y1,
         col=1,
         lwd=0.5)

#Graf rozdílů párových dat doplněný o bodový a intervalový odhad střední hodnoty,
# resp. mediánu, nárustu osmolality mezi 8 h a 11 h

#Intervalový odhad střední hodnoty (t.test) nebo mediánu (wilcox.test)
res = t.test(osmolalita$narust,conf.level=0.95)
#res = wilcox.test(osmolalita$narust,conf.level=0.95,conf.int = T)

par(mfrow = c(1,1),mar = c(2,7,2,2))
stripchart(osmolalita$narust,
           vertical=T,
           pch=16,
           method="jitter",
           main="Změna osmolality mezi 8 h a 11 h (mmol/kg)",
           ylab="Změna osmolality (mmol/kg):\n11 h - 8 h",
           cex.main = 0.8)
points(1,res$estimate,col="red",pch=16,cex=1.2) # bod definující medián
arrows(1,res$conf.int[1],
       1,res$conf.int[2],
       col="red",
       code=3, # definuje styl "šipek"
       lwd=2,
       angle=90,
       length = 0.1) # definuje délku krajních úseček
abline(h=0,lty=2) # přímka definující nulový efekt

########################################################
######## Blandův - Altmanův graf #######################
########################################################

BA_plot = function(x,y,x.lab = "(x + y)/2",y.lab = "y - x",pomer = F){
  # x ... měření 1
  # y ... měření 2
  # x.lab, y.lab ... popisky os
  # pomer ... TRUE/FALSE určuje, zda má být na vertikální osu vynášen poměr y/x
  
  if (pomer & y.lab == "y - x") {
    y.lab = "y/x"
  }
  
  diff = y - x
  diff_ln = log(y)-log(x)
  
  prumer_diff = mean(diff)
  prumer_ln = mean(diff_ln)
  sd_diff = sd(diff)
  sd_ln = sd(diff_ln)
  
  if (pomer) {
    prumer = exp(prumer_ln) # geometrický průměr poměrů
  } else {
    prumer = prumer_diff  # průměr rozdílů
  }
  
  # Interval, v němž rozdíly leží s pravděpodobností 95 %, tj. limity shody
  if (pomer) {
    lb = exp(prumer_ln - qnorm(0.975)*sd_ln)
  } else {
    lb = prumer_diff - qnorm(0.975)*sd_diff
  }
  if (pomer) {
    ub = exp(prumer_ln + qnorm(0.975)*sd_ln)
  } else {
    ub = prumer_diff + qnorm(0.975)*sd_diff
  }
  
  # Definice vykreslovaných hodnot dle pomer
  xx = (x + y)/2
  if (pomer) {
    yy = y/x
  } else {
    yy = diff
  }
  
  # Definice rozsahu osy y
  if (lb > 0){
    lo = min(0.9*lb,min(yy))
  } else {
    lo = min(1.1*lb,min(yy))
  }
  
  if (ub > 0){
    up = max(0.9*lb,max(yy))
  } else {
    up = max(1.1*lb,max(yy))
  }
  
  plot(xx,yy,
       xlab = x.lab,
       ylab = y.lab,
       pch = 16,
       ylim = c(lo,up))
  
  # Definice nulového efektu
  if (pomer) {
    e0 = 1
  } else {
    e0 = 0
  }
  
  abline(h=e0,lty = 2, col = "red") # přímka definující nulový efekt
  abline(h=prumer,lty = 2) # přímka definující průměr
  
  # Limity shody
  # přímka definující horní mez intervalu, v němž rozdíly leží s 95% p-stí
  abline(h = ub,lty = 3) 
  # přímka definující dolní mez intervalu, v němž rozdíly leží s 95% p-stí
  abline(h=lb,lty = 3) 
  
  print(paste0("Limity shody: (",round(lb,2),"; ",round(ub,2),")"))
  print(paste0("Průměr: ",round(prumer,2)))
}

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

# Aplikace funkce BA_plot
par(mfrow = c(1,1),mar = c(5,5,2,2))
BA_plot(osmolalita$o8,osmolalita$o11,pomer = F,
        x.lab = "průměr osmolality v 8 a v 11 hodin (mmol/kg)",
        y.lab = "nárůst osmolality \nmezi 8 a 11 hodinou (mmol/kg)")

# Aplikace funkce BA_plot na logaritmovaná data
BA_plot(log(osmolalita$o8),log(osmolalita$o11),
        x.lab = "průměr přirozených logaritmů osmolality v 8 a v 11 hodin (mmol/kg)",
        y.lab = "rozdíl logaritmů osmolality \nv 8 a v 11 hodin (mmol/kg)")

# Aplikace funkce BA_plot pro závislost poměrů měření ku geometrickému průměru měření
BA_plot(osmolalita$o8,osmolalita$o11,pomer = T,
        x.lab = "geometrický průměr osmolality v 8 a v 11 hodin (mmol/kg)",
        y.lab = "poměr osmolality v 11 hodin \nku osmolalitě v 8 hodin")