# Odhadu koeficientu heritability pomocí analýzy variance 

# model
# y = mi + ai + eij

# předpočítané hodnoty z databáze
Y2 <- 746983561
SY2n <- 18691786.63
SSy2 <- 18773473
p <- 5
n <- 40
n0 <- 8

# součty čtverců odchylek od průměru
  # SS mezi skupinami polosourozenců (podle otců)
SSa <- SY2n - (Y2/n) 
SSa

  # SS uvnitř skupin polosourozenců
SSe <- SSy2 -SY2n
SSe

# průměrné čtverce ~ variance
  # variance otcovská
MSa <- SSa/(p-1)
MSa

  # variance reziduální
MSe <- SSe/(n-p)
MSe

# variance genetická a prostředí
Vg <- (MSa-MSe)/n0 
 Vg
Ve <- MSe
# intraklasní korelační koeficient
r <- Vg/(Vg + Ve)
r



# heritabilita
h2 <- 4*r
h2

# střední chyba odhadu heritability

seh2 <-  4*sqrt((2*((1-r)^2)*(1+(n0-1)*r)^2/(n0*(n0-1)*(p-1))))
seh2



### výpočet ANOVA pomocí funkcí v R: lm()
### nejdříve je nutno načíst data ze souboru data-anova-cz.xlsx

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



# MSa <- 4299.4   # variance otcovská (mezi skupinami polosourozenců)
# MSe <- 2333.9   # variance reziduální/prostřeďová (uvnitř skupin polosourozenců)



# REML - balíček lme4
# doinstalujte si balíček lme4 
library(lme4)

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