# Estimating the Coefficient of Heritability Using Analysis of Variance 

# linear model
# y = mi + ai + eij

# Precalculated values from the database
Y2 <- 746983561
SY2n <- 18691786.63
SSy2 <- 18773473
p <- 5
n <- 40
n0 <- 8

# Sums of the squares of the deviations from the mean
  # SS between groups of half-sibs (by fathers)
SSa <- SY2n - (Y2/n)
SSa

  # SS within groups of half-sibs
SSe <- SSy2 - SY2n
SSe

# mean od squares ~ variance
  # sire variance
MSa <- SSa/(p-1)
MSa

  # residual variance 
MSe <- SSe/(n-p)
MSe

# genetics and environmental variance
 Vg <- (MSa - MSe)/n0
 Vg
Ve <- MSe
# intraclass correlation coefficient
r <- Vg/(Vg + Ve)
r


# heritability
h2 <- 4*r
h2


# standard error of heritability estimation

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



### ANOVA calculation using functions in R: lm()
### first it is necessary to read data from the file  data-anova-cz.xlsx
 
 ANOVA1 <- lm(y ~ sire, data = data1)
 anova(ANOVA1)



# MSa <- 4299.4   # sire variance  (between groups of half-sibs)
# MSe <- 2333.9   # residual/environmental variance (within groups of half-sibs)



# REML - package lme4
# install the package lme4 

library(lme4)

REML1 <- lmer(y ~ 1 + (1|sire), data=data1)
summary(REML1)