Fixed_Effect_model_for_SNP_effects.R

# Linear Models for the Prediction of Animal Breeding Values, Raphael Mrode, 3rd Edition
# Fixed Effect Model for SNP effects
# Date : 2020-11-11
# Julius Mugambe 


# Clear workspace
rm(list = ls())

#set working direction
setwd('K:/R-MME')
# Paste working directory
getwd()

## [1] "K:/R-MME"

# Load data
# The observations are Daughter Yield Deviations (DYDs) for fat yield
# Effective Daughter Contributions (EDC) for each bull are also given

# read data 
data = read.table("data2.txt", header = TRUE, sep = "", stringsAsFactors = FALSE)

# variance component and ratio
sigma_a = 35.241
sigma_e = 245
alpha = sigma_e/sigma_a

#pedigree 
Ped <- read.table("Ped_2.txt", header = TRUE, sep = "", stringsAsFactors = FALSE)
# view Pedigree

# Observation records for only the reference population
y <- data[1:8,6]
# View y
y

## [1]  9.0 13.4 12.7 15.4  5.9  7.7 10.2  4.8

# function to make inverse of numerator relationship matrix
ainv = function(ped) {
  n = nrow(ped)
  Ainv = matrix(c(0), nrow = n, ncol = n)
    for (i in 1:n) {
    animal = ped[i, 1]
    sire = ped[i, 2]
    dam = ped[i, 3]
    
    if (sire == 0 & dam == 0) {
      # both parents unknown
      alpha = 1
      Ainv[animal, animal] = alpha + Ainv[animal, animal]
    } else if (sire != 0 & dam == 0) {
      # sire known
      alpha = 4/3
      Ainv[animal, animal] = alpha + Ainv[animal, animal]
      Ainv[animal, sire] = -alpha/2 + Ainv[animal, sire]
      Ainv[sire, animal] = -alpha/2 + Ainv[sire, animal]
      Ainv[sire, sire] = alpha/4 + Ainv[sire, sire]
    } else if (sire == 0 & dam != 0) {
      # dam known
      alpha = 4/3
      Ainv[animal, animal] = alpha + Ainv[animal, animal]
      Ainv[animal, dam] = -alpha/2 + Ainv[animal, dam]
      Ainv[dam, animal] = -alpha/2 + Ainv[dam, animal]
      Ainv[dam, dam] = alpha/4 + Ainv[dam, dam]
    } else {
      # both parents known
      alpha = 2
      Ainv[animal, animal] = alpha + Ainv[animal, animal]
      Ainv[animal, sire] = -alpha/2 + Ainv[animal, sire]
      Ainv[sire, animal] = -alpha/2 + Ainv[sire, animal]
      Ainv[animal, dam] = -alpha/2 + Ainv[animal, dam]
      Ainv[dam, animal] = -alpha/2 + Ainv[dam, animal]
      Ainv[sire, sire] = alpha/4 + Ainv[sire, sire]
      Ainv[sire, dam] = alpha/4 + Ainv[sire, dam]
      Ainv[dam, sire] = alpha/4 + Ainv[dam, sire]
      Ainv[dam, dam] = alpha/4 + Ainv[dam, dam]
    }
  }
  return(Ainv)
}
# A-inverse from the full pedigree from all the animals
A <- ainv(Ped)

#Constructing the Z- matrix for only the 8 animals to be considered for the first three SNPS
# Z = M-P. 

M <-as.matrix(data[,7:9]) # considering only the 3 first SNPs

P <- (colMeans(M)/2)

# Use sweep function to make matrix Z from M and the 2nd allele frequencies for all the SNPs of the 8 animals
Z_all <- sweep(M,2,P*2)

#For the first 8 animals (reference population)
Z <- Z_all[1:8,]

D <- (round(diag(A),3))