#########################################################
# 7. týden

# Příklad 1: 
# data data-reg1.xlsx
# lineární regrese mezi hodnotami vlastnosti potomků (z) a rodičů - matek (x) > odhad heritability
reg <- lm(z ~ x, data=data_reg1)
summary(reg)

# výpočet heritability - násobíme regresní koeficient 2, protože se jedná o vztah matka-dcera
b <- 0.21505
h2 <- 2*b
h2
# střední chyba odhadu heritability (ze střední chyby odhadu regresního koeficientu * 2)
se_reg <- 0.09903
se_h2 <- 2*se_reg   
se_h2

# Příklad 2
# data data-reg2.xlsx
# odhad heritability na základě regrese potomků na matku uvnitř otců 

reg2 <- lm(PotomciZ ~ Otcové + MatkyX, data=data_reg2)
summary(reg2)
# regresní koeficient b = 0,1003
b <- 0.1003
h2 <- 2*b
h2

se_reg <- 0.3354
se_h2 <- 2*se_reg   
se_h2

#########################################################
# 8. týden
# odhad heritability pomocí lineárních modelů ANOVA a smíšených lineárních modelů REML

# Příklad 1: 1f anova 
# data data-anova1.xlsx

# modelová rovnice:  y = mi + a + e
# můžeme použít funkce aov() nebo lm()

anova1 <- lm(y ~ otec, data=data_anova1)
anova(anova1)

# vybrat hodnoty MS/variancí mezi skupinami a uvnitř skupin z tabulky ANOVA
MSa <- 4299.4
MSe <- 2333.9
n0 <- 8 

# zde končí statistika a začíná genetika
# genetická část výpočtu:
# výpočet variancí
Vg <- (MSa-MSe)/n0
Vg
Ve <- MSe
# intraklasní korelační koeficient
r <- Vg/(Vg+Ve)
r
# výpočet heritability
h2 <- 4*r
h2

# smíšené lineární modely, maximum likelyhood - REML
# nutno si nainstalovat balíček lme4 (v záložce Packages > Install - najít lme4 a nainstalovat)

library(lme4) # načtení balíčku

reml1 <- lmer(y ~ 1 + (1|otec), data=data_anova1)
summary(reml1)

# genetická část - výpočet heritability, je stejný
# porovnejte výsledky

 
# Příklad 2: 2f hierarchická anova 
# data data-anova2.xlsx

#  y = mi + otec + matka(otců) + e
anova2 <- lm(y ~ OTEC/MATKA, data=data_anova2)
anova(anova2)
# vybrat hodnoty variancí z tabulky ANOVA
MSa <- 1426.16
MSe <- 215.48
MSb <- 449.13
p <- 5 # počet otců
m <- 26 # počet matek
n <- 218  # počet potomků

k1 <- n/m   # průměrný počet potomků pro matku
k2 <- m/p   # průměrný počet matek pro otce
k3 <- n/p   # průměrný počet potomků pro otce

# zde končí statistika a začíná genetika
# genetická část výpočtu:
# výpočet variancí
VgM <- (MSb-MSe)/k1
VgM

VgO <- (MSa - MSe - (k2*VgM))/k3
VgO
Vp <- VgO + VgM + MSe

h2o <- 4*(VgO/Vp)
h2o
se1 <- (h2O + (4/(n/p)))*(sqrt(2/p))
se1

h2m <- 4*(VgM/Vp)
h2m
se2 <- (h2M + (4/(n/m)))*(sqrt(2/m))
se2

h2OM <- 2*((VgO+VgM)/Vp)
h2OM
se3 <- 4 * sqrt(h2OM) * sqrt(1/n)
se3

# reml

reml2 <- lmer(y ~ 1 + (1|OTEC/MATKA), data=data_anova2)
summary(reml2)
VgM <- 28.87
VgO <- 21.42
Ve <- 215.02

#... počítání heritabilit je pak stejné jako po ANOVě
# rozdíl mezi odhady pomocí ANOVA a REML je z důvodu nevybalancování dat (různý počet potomků po matkách)
table(data_anova2$MATKA)

