rm(list=ls())
library(xlsx)
## Loading required package: rJava
## Loading required package: xlsxjars
library(synbreed)
library(regress)
library(gwaR)

#only run these two lines if you need to compute A
#if you will load A, just comment this lines and ignore them
setwd("C:/Users/steibelj/OneDrive/Documents/job/ans824/2017")
load("MSUPRP_meat.RData")

#load Andres data
setwd("C:\\Users\\steibelj\\OneDrive\\Documents\\GWA\\")
andy_dts<-read.xlsx("Datos.xls",sheetIndex=1)
head(andy_dts)
##     ID sex      TITAH      TITAS       TITAS2 lrf_10wk
## 1 1001   M  0.0833082  0.0247694 0.0006135232     6.35
## 2 1002   F -0.0270106  0.1880564 0.0353652096     9.14
## 3 1003   M  0.1188594 -0.0396142 0.0015692848     6.35
## 4 1004   F  0.0093282 -0.1463548 0.0214197275     5.08
## 5 1006   F -0.0070914  0.0105274 0.0001108262     5.59
## 6 1007   M -0.0932860 -0.0781200 0.0061027344     5.84
#compute A matrix and save it or load it if it has been previously saved
#(select by commenting the appropriate lines)

A_matrix<-kin(MSUPRP_meat,ret="add")
save(A_matrix,file="A_MSU.Rdata")
#load("A_MSU.Rdata")

#take a peak at a random diagonal block
A_matrix[100:105,100:105]
##        1032   1033   1034   1035   1036   1037
## 1032 1.0000 0.0625 0.5000 0.0625 0.5000 0.0625
## 1033 0.0625 1.0000 0.0625 0.5000 0.0625 0.5000
## 1034 0.5000 0.0625 1.0000 0.0625 0.5000 0.0625
## 1035 0.0625 0.5000 0.0625 1.0000 0.0625 0.5000
## 1036 0.5000 0.0625 0.5000 0.0625 1.0000 0.0625
## 1037 0.0625 0.5000 0.0625 0.5000 0.0625 1.0000
## attr(,"class")
## [1] "relationshipMatrix" "matrix"
#set animal ID as rownames of input data
#this is required by gwaR functions to work correctly
rownames(andy_dts)<-andy_dts$ID
head(andy_dts)
##        ID sex      TITAH      TITAS       TITAS2 lrf_10wk
## 1001 1001   M  0.0833082  0.0247694 0.0006135232     6.35
## 1002 1002   F -0.0270106  0.1880564 0.0353652096     9.14
## 1003 1003   M  0.1188594 -0.0396142 0.0015692848     6.35
## 1004 1004   F  0.0093282 -0.1463548 0.0214197275     5.08
## 1006 1006   F -0.0070914  0.0105274 0.0001108262     5.59
## 1007 1007   M -0.0932860 -0.0781200 0.0061027344     5.84

Now we have all the elements needed to use our gwaR package. This package can fit gaussian mixed models. By default it includes an additive genetic effect with arbitrary variance-covariance matrix. In our example we have no extra random effects. But we will have fixed effects of sex, TITAH, TITAS and TITAS2.

#fit animal model
mod_res<-gblup(rsp = "lrf_10wk",data = andy_dts,      #response and dataset
               design = c(y~sex+ TITAH+TITAS+TITAS2), #fixed effectsmodel
               G=A_matrix)                          #var-cov for additive eff

once the mixed model object is obtained we can look into different summaries for variance component estimation and fixed effect estimates and tests

#check estimation of variance components and h2
varcomp(mod_res) #notice h2=0.3017
##     Estimate   StdError  prop.var         se
## G  0.3551757 0.12229010 0.3017387 0.09582012
## In 0.8219214 0.08108272 0.6982613 0.11731924
#fixed effects
print(mod_res)
## gblup analysis of trait: lrf_10wk 
## 
## fixed effects equation:
## y ~ sex + TITAH + TITAS + TITAS2
## 
## random effects equation:
## ~G + In
## 
## log-likelihood: -479.824 converged in: 5 iterations 
## 
## estimated fixed effects:
##              Estimate   StdError   test.st      p.value
## (Intercept) 5.9708793 0.20356791 29.331141 0.0000000000
## sexM        0.1825175 0.06702912  2.722958 0.0064700241
## TITAH       0.5578666 0.33563218  1.662137 0.0964853512
## TITAS       1.2686541 0.37147047  3.415222 0.0006373007
## TITAS2      4.8124005 2.42955482  1.980775 0.0476165577
## 
## estimated variance components:
##     Estimate   StdError  prop.var
## G  0.3551757 0.12229010 0.3017387
## In 0.8219214 0.08108272 0.6982613
#TitaH not sig. TITAS signif TITAS2 marginally significance
#these are type 3 tests.

Notice: our program DELETES any column and row of A that does not have a corresponding row in the phenotypic data. Consequently this is similar to what SAS should be doing. If SAS is obtaining a different result, perhaps there is an error somewhere in the input data.

Alternative fixed effects

I run this one to use the same parameterization used by Sebastian

#fit animal model
mod_res<-gblup(rsp = "lrf_10wk",data = andy_dts,      #response and dataset
               design = c(y~sex+ TITAS-1), #fixed effectsmodel
               G=A_matrix)                          #var-cov for additive eff
print(mod_res)
## gblup analysis of trait: lrf_10wk 
## 
## fixed effects equation:
## y ~ sex + TITAS - 1
## 
## random effects equation:
## ~G + In
## 
## log-likelihood: -482.5874 converged in: 5 iterations 
## 
## estimated fixed effects:
##       Estimate  StdError  test.st     p.value
## sexF  6.007795 0.2022734 29.70137 0.000000000
## sexM  6.195052 0.2016487 30.72201 0.000000000
## TITAS 1.147982 0.3669057  3.12882 0.001755097
## 
## estimated variance components:
##     Estimate   StdError  prop.var
## G  0.3528689 0.12178873 0.2990339
## In 0.8271611 0.08097719 0.7009661
varcomp(mod_res)
##     Estimate   StdError  prop.var         se
## G  0.3528689 0.12178873 0.2990339 0.09522211
## In 0.8271611 0.08097719 0.7009661 0.11684917