# Background

The following random regression (RR) model was used for all subsequent analyses

\begin{align} \label{RR2} \text{PSA}_{tjk} =\mu + \sum_{k=0}^{2}\phi(t)_{jk}\beta_k + \sum_{k=0}^{2}\phi(t)_{jk} u_{jk} + \sum_{k=0}^{1}\phi(t)_{jk} s_{jk} + e_{tjk} \end{align}

As mentioned previously, $$\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. The random additive genetic effects are described using a second-order Legendre polynomial, while a first-order Legendre polynomial is used to describe the experiment effects across time points.

The random regrression model was compared to a conventional, single time point model (TP). The model is as follows \begin{align} \label{TP} \mathbf{y} &= \mathbf{Zu} + \mathbf{Qs} + \mathbf{e} , \end{align}

Here, $$\mathbf{y}$$ is the PSA at time $$t$$; $$\mathbf{Z}$$ and $$\mathbf{Q}$$ are incidence matrices corresponding to the random additive genetic effect ($$\mathbf{u}$$), and random experimental effect ($$\mathbf{s}$$), respectively; and $$\mathbf{e}$$ is the random residual error. The variances were based on the following assumptions $$\mathbf{u} \sim N(0, \mathbf{G}\sigma_g^2)$$, $$\mathbf{s} \sim N(0, \mathbf{I}\sigma_s^2)$$, and $$\mathbf{e} \sim N(0, \mathbf{I}\sigma_e^2)$$. Here, $$\sigma_g^2$$ is the additive genetic variance; $$\sigma_s^2$$ is an environmental variance associated with experiment; and $$\sigma_e^2$$ is the residual variance.

For both the RR and TP appraches a genomic relationship matrix ($$\mathbf{G}$$) was calculated using VanRaden 2008. \begin{align} \label{G} \mathbf{G} = \frac{\mathbf{Z_{cs} Z_{cs}}' }{m} \end{align}

Here, $$\mathbf{Z_{cs}}$$ is a centered and scaled $$n \times m$$ matrix, where $$m$$ is 33,674 SNPs and $$n$$ is the 357 genotyped rice lines.

To estimate the narrow sense heritability, variance components were obtained for each random term using ASREML for the TP analyses and the RR approach. For the RR approach, additive genetic variance was obtained at each time points using methods described by Mrode (2014). Briefly, for time $$i$$ the genetic variance can be obtained by $$\mathbf{t}_{i} \mathbf{\Omega} \mathbf{t}_{i}'$$, where $$\mathbf{t}_{i} = \phi_{ik}$$, the $$i$$th row vector of the matrix of Legendre polynomials at different time points ($$\mathbf{\phi}$$) for the $$i$$th day of imaging, $$\mathbf{\Omega}$$ is the covariance matrix of RR coefficients for the genetic effects, and $$k$$ is the order of fit. The variance of the experimental effect across time points was calculated using the same approach. For both the single time point analysis $$h^2$$ was estimated as $$\frac{\sigma_g^2}{\sigma_s^2 + \sigma_s^2 + \sigma_e^2}$$.

# Running the models

## as file for random regression

h2 RR model
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 1

!PATH 1 #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

# Estimating $$h^2$$

## Single time point

For the TP analysis, 20 .asr files were created. This chunk of code parses those .asr files, extracts the variance components, and writes those to a file. The r code loads each file and creates an object with the variance components at each time point. ###Parsing asreml output

cd ~/Desktop/RR/TP/h2
for file in *.asr; do less $file | grep -A6 "Model_Term" | sed -n "1,4p" | tr -s " " >$file.varcomp.txt; done

### Estimating $$h^2$$ for single time point approach

library(reshape2)

files <- paste0("TP/h2/TPY", 1:20, ".asr.varcomp.txt")
for(i in 1:20){
tmp\$DOI <- i
if(i == 1){
final <- tmp
}else{
final <- rbind(final, tmp)
}
}
final <- final[c("V1", "V5", "DOI")]
colnames(final)[1:2] <- c("Source", "Sigma")

Var.TP <- dcast(final, Source ~ DOI, value.var = "Sigma")
Var.TP <- Var.TP[2:21]
row.names(Var.TP) <- c("Exp", "G", "E")

Var.TP.per <- rbind(Var.TP[1,]/colSums(Var.TP)*100, Var.TP[2,]/colSums(Var.TP)*100, Var.TP[3,]/colSums(Var.TP)*100)
TP.h2 <- Var.TP[2,]/colSums(Var.TP)

## Random regression

$$h^2$$ for the RR model is straightforward, and is done as described above. The variance components below are obtained from the .asr file.

### Defining some functions

The functions below were adapted from Mrode (2005) by Gota Morota.

##Return coefficient matrix (lambda) of n-th order Legendre polynomials. Scaling method implemented by Gengler et. al. (1999) converts constant Legendre polynomial coefficients into 1
legendre <-
function(n, gengler){

if (nargs()==1){
gengler <- TRUE
}

if (gengler != TRUE & gengler != FALSE){
gengler=TRUE
}

N <- n+1
L <- matrix(0,nrow=N, ncol=N)

for(i in (1:N)){
if(i==1){
L[i,i] <- 1
}
else if(i==2){
L[i,i] <- 1
}
else  {
tmp <- L[i-1,]
tmp2 <- as.numeric()
tmp2 <- c(0,tmp[1:(N-1)])
L[i,] <- (1/(i-2+1))*( (2*(i-2) + 1)*tmp2 -(i-2)*L[i-2,] )
}
}

# Normalize
for (j in (1:N)){
L[j,] <- (sqrt( (2*(j-1)+1)/2)  )*L[j,]
}

# Gengler (1999)
if (gengler==TRUE){
L <- sqrt(2)*L
}

return(L)

}

##Given time points covariate and order of fit for Legendre polynomials, return matrix 'M' containing the polynomials of standardized time. 'M' is order t (number of time points) by k (order of Legendre polynomials)
stdtime <- function(t, n, tmax, tmin){
if(missing(tmax)) {
tmax <- t[which.max(t)]
}
if(missing(tmin)) {
tmin <- t[which.min(t)]
}

N <- n+1
M <- matrix(0, nrow=length(t), ncol=N)
a <- -1 + 2*(t-tmin)/(tmax - tmin)
M[,1] <- 1

for (i in 2:N){
M[,i] <- a^(i-1)
}

return(M)
}

### Estimating $$h^2$$ for random regression approach

#Variance components from the .asr file
##Residual variance. Recall that idh was used, so each time point will have a unique variance
Resvar <- c(1.39E-03, 1.22E-03, 2.06E-03, 3.33E-03,
5.32E-03, 7.88E-03, 1.18E-02, 1.53E-02,
2.04E-02, 2.75E-02, 4.02E-02, 5.80E-02,
8.19E-02, 0.114587, 0.164339, 0.22392,
0.336068, 0.521511, 0.78897, 1.17139)

#Genetic variance. Here a second order Legendre polynomial was used. So G will be a 3 x 3 matrix of
G.mat <- matrix(c(0.629183, 0.495574, 0.129607,
0.495574, 0.394557, 0.104540,
0.129607, 0.104540, 0.281278E-01), 3, 3)

E.mat <- matrix(c(0.594161E-01, 0.392295E-01,
0.392295E-01, 0.260640E-01), 2, 2)

Phi <- stdtime(1:20, 2) %*% t(legendre(2, gengler = F))
G.COV <- Phi %*% G.mat %*% t(Phi)
Gvar <- diag(G.COV)

Phi <- stdtime(1:20, 1) %*% t(legendre(1, gengler = F))
E.COV <- Phi %*% E.mat %*% t(Phi)
Evar <- diag(E.COV)

Var.RR <- rbind(Evar, Gvar, Resvar)
Var.RR.per <- rbind(Evar/colSums(Var.RR)*100, Gvar/colSums(Var.RR)*100, Resvar/colSums(Var.RR)*100)
RR.h2 <- Var.RR[2,]/colSums(Var.RR)

## Plotting $$h^2$$

Here, I’ll plot $$h^2$$ and the variance components for the RR and TP approaches. This chunk of code will generate Figure 4 in the paper.

#single column figure 3.25in width
#pdf("~/Desktop/RR/RR/Figures/Fig4.pdf", h=6, w=3.25, useDingbats = F, pointsize = 10)
par(mar=c(3,3,2,.2), mgp=c(1.8,0.5,0))
nf=layout(rbind(c(1,1), c(2,2), c(3,3)))

plot(13:32, RR.h2, pch=19, cex=0.3, ylab=expression(italic(h^2)), xlab="Days After Transplant", col="black", ylim=c( 0, 1) )
lines(13:32, RR.h2, col="black")

points(13:32, TP.h2, pch=21, cex=0.3, col="grey")
lines(13:32, TP.h2, lty=2, cex=0.3, col="grey")

legend("bottomright", c("TP", "RR"), col=c("grey", "black"), pch=c(21, 19), lty=c(2, 1), cex=1, pt.cex = 0.3)

###Varmat for TP
plot(13:32, Var.TP[2,], pch=19, cex=0.5, ylab=expression(italic(sigma^2)), xlab="Days After Transplant", col="black", ylim=c( 0, max(colSums(Var.TP))*1.05) )
lines(13:32, Var.TP[2,], col="black")

points(13:32, Var.TP[1,], pch=20, col="grey", cex=0.5)
lines(13:32, Var.TP[1,], col="grey")

points(13:32, colSums(Var.TP), pch=21, cex=0.5)
lines(13:32, colSums(Var.TP), col="black", lty=2)

legend("topleft",
c(expression(italic({sigma^{2}}[g])),
expression(italic({sigma^{2}}[Exp])),
expression(italic(sigma^2))),
lty=c(1,1,2),
col=c("black", "grey", "black"),
pch=c(19, 20, 21),
cex=1,
pt.cex = 0.5)

###Varmat for RR
plot(13:32, Var.RR[2,], pch=19, cex=0.5, ylab=expression(italic(sigma^2)), xlab="Days After Transplant", col="black", ylim=c( 0, max(colSums(Var.RR))*1.05) )
lines(13:32, Var.RR[2,], col="black")

points(13:32, Var.RR[1,], pch=20, col="grey", cex=0.5)
lines(13:32, Var.RR[1,], col="grey")

points(13:32, colSums(Var.RR), pch=21, cex=0.5)
lines(13:32, colSums(Var.RR), col="black", lty=2)

legend("topleft",
c(expression(italic({sigma^{2}}[g])),
expression(italic({sigma^{2}}[Exp])),
expression(italic(sigma^2))),
lty=c(1,1,2),
col=c("black", "grey", "black"),
pch=c(19, 20, 21),
cex=1,
pt.cex = 0.5)