---
title: "Random regression model selection"
author: "Malachy Campbell"
date: "8/7/2017"
output:
  rmdformats::html_clean:
    fig_width: 6
    fig_height: 6
    highlight: kate
    thumbnails: true
    lightbox: true
    gallery: true
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

#Background
The purpose of this script is to determine which RR model best fits the data. We will first fit the simplest model in asreml and gradually move to a more complex model. The general structure of the RR model is 
$$ PSA_{tjk} = \mu + \sum_{k=0}^{2}\phi(t)_{jk}\beta_k + \sum_{k=0}^{nr}\phi(t)_{jk} u_{jk} + \sum_{k=0}^{nr}\phi(t)_{jk} s_{jk} + e_{tjk}$$. Where $\beta$ is the fixed second-order Legendre polynomial to model the overall trend in the trait overtime, $u_{jk}$ and $s_{jk}$ are the $k^{th}$ random regression coefficients for additive genetic effect and random experiment of line $j$, $nr$ is the order of polynomial for the random effects, and $e_{tjk}$ is the random residual. $\beta$ was selected based on visual inspection of the mean trend in PSA. Various polynomial functions and residual variance structures are evaluated for line and experiment, and residuals, respectively. For each trait, the models were ranked based on goodness-of-prediction using Akaike's information criterion (AIC) scores. 

Since asreml-R does not allow the use of Legendre polynomials all analysis were done with the standalone version of ASREML. Below the ".as" file is provided as well as the call from the commandline.

```{r load packages, echo=F}
library(knitr)
```
#Running the models
##as file for random regression model selection
```{asreml .as file}
RR model selection
 NID !A
 Exp !A
 Rep 2
 DayOfImaging 20 !I
 PSA !/100000
G2.grm
PSA.cleaned.csv !SKIP 1 !MAXITER 1000 !EXTRA 100 !WORKSPACE 6144 !ASUV !DOPATH 2

!PATH 1 #Model1
!ASSIGN USg !< !INIT
0.254496
!>
PSA ~ mu leg(DayOfImaging,2) !r leg(DayOfImaging,0).Exp leg(DayOfImaging,0).grm(NID) !f mv

!PATH 1
PSA ~ mu pol(DayOfImaging,2) !r Exp leg(DayOfImaging,2).NID !f mv
residual id(2142).us(DayOfImaging)

!PATH 2 #Model2: 19358.83 -38673.65 -38495.60
!ASSIGN USe !< !INIT
1.60E-03 1.52E-03 2.06E-03 3.00E-03 4.59E-03 7.12E-03 1.13E-02 1.82E-02 2.86E-02 4.11E-02 6.09E-02 9.26E-02 1.33E-01 1.87E-01 2.69E-01 0.349805 0.503757 0.740224 1.06494 1.52276
!>
!ASSIGN USg !< !INIT
0.254496
!>
PSA ~ mu leg(DayOfImaging,2) !r leg(DayOfImaging,0).Exp us(leg(DayOfImaging,0) $USg).grm(NID) !f mv
residual id(2142).idh(DayOfImaging $USe)

!PATH 3 #Model3: 
!ASSIGN USg !< !INIT
0.254496
0.200363 0.159373
!>
PSA ~ mu leg(DayOfImaging,2) !r leg(DayOfImaging,0).Exp us(leg(DayOfImaging,1) $USg).grm(NID) !f mv

!PATH 4 #Model4: 23273.62 -46499.24 -46305.01
!ASSIGN USe !< !INIT
1.60E-03 1.52E-03 2.06E-03 3.00E-03 4.59E-03 7.12E-03 1.13E-02 1.82E-02 2.86E-02 4.11E-02 6.09E-02 9.26E-02 1.33E-01 1.87E-01 2.69E-01 0.349805 0.503757 0.740224 1.06494 1.52276
!>
!ASSIGN USg !< !INIT
0.254496
0.200363 0.159373
!>
PSA ~ mu leg(DayOfImaging,2) !r leg(DayOfImaging,0).Exp us(leg(DayOfImaging,1) $USg).grm(NID) !f mv
residual id(2142).idh(DayOfImaging $USe)

!PATH 5 #Model5: 
!ASSIGN USg !< !INIT
0.667162
0.587027 0.519362
0.185506 0.164992 0.526885E-01
!>
PSA ~ mu leg(DayOfImaging,2) !r leg(DayOfImaging,0).Exp us(leg(DayOfImaging,2) $USg).grm(NID) !f mv

!PATH 6 #Model6: 24718.93 -49383.86 -49165.35
!ASSIGN USe !< !INIT
1.60E-03 1.52E-03 2.06E-03 3.00E-03 4.59E-03 7.12E-03 1.13E-02 1.82E-02 2.86E-02 4.11E-02 6.09E-02 9.26E-02 1.33E-01 1.87E-01 2.69E-01 0.349805 0.503757 0.740224 1.06494 1.52276
!>
!ASSIGN USg !< !INIT
0.254496
0.200363 0.159373
0.514966E-01 0.414476E-01 0.109234E-01
!>
PSA ~ mu leg(DayOfImaging,2) !r leg(DayOfImaging,0).Exp us(leg(DayOfImaging,2) $USg).grm(NID) !f mv
residual id(2142).idh(DayOfImaging $USe)

!PATH 7 #Model7: 
!ASSIGN USg !< !INIT
0.685896
0.589393 0.509177
0.190862 0.166014 0.545985E-01
!>
!ASSIGN USp !< !INIT
0.162107
0.151803 0.147320
!>
PSA ~ mu leg(DayOfImaging,2) !r us(leg(DayOfImaging,1) $USp).Exp us(leg(DayOfImaging,2) $USg).grm(NID) !f mv

!PATH 8 #Model8: 27537.59 -55017.19 -54782.49
!ASSIGN USe !< !INIT
1.37E-03 1.22E-03 2.07E-03 3.32E-03 5.32E-03 7.88E-03 1.18E-02 1.53E-02 2.04E-02 2.76E-02 4.03E-02 5.80E-02 8.19E-02 0.114627 0.164452 0.224128 0.336555 0.522506 0.790597 1.17392
!>
!ASSIGN USg !< !INIT
0.629183
0.495574 0.394557
0.129607 0.104540 0.281278E-01
!>
!ASSIGN USp !< !INIT
0.594161E-01
0.392295E-01 0.260640E-01
!>
PSA ~ mu leg(DayOfImaging,2) !r us(leg(DayOfImaging,1) $USp).Exp us(leg(DayOfImaging,2) $USg).grm(NID) !f mv
residual id(2142).idh(DayOfImaging $USe)
```

#Run .as file with asreml
Here, is what a typical run looks like. This was done for eight models. The results are summarised in the table I of the manuscript.

```{bash, eval = F}
cd ~/Desktop/RR/RR/ModelSelection/
asreml ModelSel.as
```

The "best" model is model 8 with an $AIC$ of -55,017.19, $BIC$ of -54,782.49, and log likelihood of 27,537.59.