hen_data

load libraries

library(pedigree)

## Loading required package: Matrix

library(pedigreemm)

## Loading required package: lme4

library(lme4)
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(gwaR)

## 
## Attaching package: 'gwaR'

## The following object is masked from 'package:pedigree':
## 
##     gblup

library(regress)
library(Matrix)

load health phenotypes data in July:

load("hen_social_data.RData")
health_data_July_all<-hen_social_data[["health_data_July_all"]]
dim(health_data_July_all)

## [1] 1067   19

colnames(health_data_July_all)

##  [1] "date"           "ID"             "assessor"       "pen"           
##  [5] "neck"           "back"           "wing"           "tail"          
##  [9] "cloaca"         "breast"         "plumage"        "footpad"       
## [13] "toes"           "mass"           "age"            "sire"          
## [17] "barn1"          "transform_toes" "transform_tail"

load A matrix: Am_reduced (after filtering sire)

Am_reduced <- hen_social_data[["Am_reduced"]]
dim(Am_reduced)

## [1] 1070 1070

load G matrix: match all phenotypes

GRMs_match <- hen_social_data[["GRMs_match"]]

load association matrix in week 1 (association_matrix_week1), file name week 1, by time order it is week 10 (with all 5 pens)

association_matrix_week1 <- hen_social_data[["association_matrix_week1"]]
dim(association_matrix_week1)

## [1] 1063 1063

run gblup and EM-REML

# match ID in phenotype and association matrix
common_rows <- intersect(rownames(association_matrix_week1), rownames(health_data_July_all))
association_matrix_week1 <- association_matrix_week1[common_rows, common_rows]
health_data_July <- health_data_July_all[common_rows,]

# GBLUP
July_model_tail <-gblup(rsp="tail",data=health_data_July,design = c(y~1,~pen),GRMs_match,pos = c(TRUE,TRUE,TRUE))
varcomp(July_model_tail)

##      Estimate  StdError   prop.var         se
## G    54.13590 12.488245 0.23353643 0.05147107
## pen  13.17127  9.976046 0.05681942 0.04080168
## In  164.50207 10.851562 0.70964415 0.06918962

# EM-REML function, need to run the isge function in the isge file first.
# use initial value for social genetic effects at 60
source("./isge_pen_random_initial.R")
July_model_tail.cov.isge <-igest(July_model_tail,association_matrix_week1,tol = 10^-5,k_iter = 500)

## [1] "The Tolerance criterion has been met"

July_model_tail.cov.isge$sigma

##            Estimate
## sigma u   52.662078
## sigma c   59.473578
## sigma uc   6.469835
## sigma p    3.371895
## sigma e  162.377534