Linear Mixed Model

Estimate of fixed effects and prediction of random effects is given by:

A brief history of ASReml

• The core algorithm of ASReml by Gilmour et al. (2002) is written in Fortran.
• asreml R package by Butler et al. (2009) wraps the core algorithm.

• ASReml uses REML (Patterson and Thompson, 1971) and AI algorithm and sparse matrix methods (Gilmour et al. 1995) to fit the specified (generalised) linear mixed models.

Illustrative Example^*: Multi-environmental Field Trial

Model Fitting in asreml

asreml(fixed=response ~ fixed formula,
       random= ~ random formula,
       rcov= ~ rcov formula,
       data=data.frame)

Model 1: Simple Model

Model 1 continued

fit1 is equivalent to

fit1a <- asreml(Yield ~ Site, data=data1
                random= ~ idv(Site):id(Geno))
fit1b <- asreml(Yield ~ Site, data=data1,
                random= ~ id(Site):idv(Geno))

where the variance structures of random effects are made explicit. Additionally we can define the residual error variance structure explicit:

fit1c <- asreml(Yield ~ Site, data=data1,
                random= ~ Site:Geno,
                rcov= ~ idv(units))

where units is a special reserve word in asreml equivalent to factor(1:nrow(data1)).

Variance vs. Correlation Structure

• id() is a correlation structure \(\boldsymbol{I}_n\)
• idv() is a variance structure with common (homogeneous) variance \(\sigma^2\boldsymbol{I}_n\)
• idh() (which is the same as diag()) is a variance structure with a separate (heterogeneous) variance for each level of factor \(\text{diag}(\sigma^2_1, ..., \sigma^2_n)\).

Appending the name of the correlation function with:

• v results in adding a homogeneous variance or
• h results in heterogeneous variance.

Separable Structure

• Given a random term A:B where A and B are factors then the generic form is f(A):g(B) where

  • f(A) generates a variance matrix G_A say, and
  • g(B) generates a variance matrix G_B say.

  Then the corresponding random effects, , has properties:

  • ; and
  • .

• Note assumption of separability

  • differs to assuming a completely unstructured variance ; and
  • greatly reduces computational load.

Model 2: MET Unstructured Model

• Same genotype across sites
• General idea: borrow strength across sites to get more accurate Site:Geno prediction.

summary(fit2, all=T) %>% str()

  List of 9
   $ call        : language asreml(fixed = Yield ~ Site, random = ~us(Site):id(Geno), data = data1)
   $ loglik      : num -13.4
   $ nedf        : int 15
   $ sigma       : num 1.18
   $ varcomp     :'data.frame':    7 obs. of  5 variables:
    ..$ gamma     : num [1:7] 0.15 0.221 0.462 0.663 0.156 ...
    ..$ component : num [1:7] 0.209 0.308 0.644 0.925 0.217 ...
    ..$ std.error : num [1:7] 0.118 0.593 0.733 0.846 0.592 ...
    ..$ z.ratio   : num [1:7] 1.774 0.52 0.878 1.093 0.366 ...
    ..$ constraint: Factor w/ 3 levels "?","Positive",..: 3 2 2 2 2 1 2
    .. ..- attr(*, "names")= chr [1:7] "Site:Geno!Site.A:A" "Site:Geno!Site.B:A" "Site:Geno!Site.B:B" "Site:Geno!Site.C:A" ...
   $ distribution: chr "gaussian"
   $ link        : chr "identity"
   $ coef.fixed  : num [1:4, 1:3] 0 -1.04 0.156 5.618 NA ...
    ..- attr(*, "dimnames")=List of 2
    .. ..$ : chr [1:4] "Site_A" "Site_B" "Site_C" "(Intercept)"
    .. ..$ : chr [1:3] "solution" "std error" "z ratio"
   $ coef.random : num [1:9, 1:3] -0.238 -0.301 0.539 -0.121 -0.385 ...
    ..- attr(*, "dimnames")=List of 2
    .. ..$ : chr [1:9] "Site_A:Geno_G1" "Site_A:Geno_G2" "Site_A:Geno_G3" "Site_B:Geno_G1" ...
    .. ..$ : chr [1:3] "solution" "std error" "z ratio"
   - attr(*, "class")= chr "summary.asreml"

Model 3: G-BLUP Model

• Genotypes are often related with its relatedness inferred from pedigree or genetic markers.

• This known structure (denoted \(\textbf{G}_g\) here) is fitted in the variance structure as follows:
  

Say,

Model 4: Multi-Trait MET Model

• It is perhaps more realistic to fit a model where the residual error variance differ at different trials.

• at(f,l) conditions on level l of factor f. If l is not specified, conditions on each level of f.

Analysis on Subset of Data

fit4A1 <- asreml(Yield ~ Site, random=~Site:Geno, 
                 data=data1, subset=Site=="A")

fit4A2 <- asreml(Yield ~ 1, random=~ Geno,
                 data=data1, subset=Site=="A")

dataA <- subset(data1, Site=="A")
fit4A3 <- asreml(Yield ~ Site, random=~Site:Geno, 
                 data=dataA)

Residuals

resid(fit4A1)

fit4A <- asreml(Yield ~ 1, random=~ Geno,
                data=data1, subset=Site=="A")
fit4B <- asreml(Yield ~ 1, random=~ Geno,
                data=data1, subset=Site=="B")
fit4C <- asreml(Yield ~ 1, random=~ Geno,
                data=data1, subset=Site=="C")
fit4 <- asreml(Yield ~ Site, data=data1, 
      random=~ diag(Site):Geno, # idh(Site):Geno
      rcov=~ at(Site):units)

lucid::vc(fit4)

                 effect component std.error z.ratio constr
   Site:Geno!Site.A.var    0.9034    1.618     0.56      P
   Site:Geno!Site.B.var    1.349     1.669     0.81      P
   Site:Geno!Site.C.var    0.9461    1.54      0.61      P
        Site_A!variance    1.261     1.029     1.2       P
        Site_B!variance    0.6029    0.4922    1.2       P
        Site_C!variance    1.063     0.868     1.2       P

Model 5: Spatial Model

• It is perhaps more realistic to fit a model such that plots that are closer together geographically are more correlated than plots further away from each other.
• We capture this by considering an autoregressive of order 1 structure for row and column direction in the field.

Model 6: MET Factor Analytic Model

• In practice, we replace the unstructured variance models with factor analytic models.

More explanations on these models later today.

Model 7: MPI Model

• We can add multiple (crossed) random effects and multiple givs. Here the model is fitting additive marker effects, residual additive effects (via pedigree), and non-additive effects (capturing interaction effects).

Model 8: Multivariate/Multi-trait Model

• Suppose you record multiple traits for each genotype and you wish to jointly model these traits.
• asreml has a reserve word trait which is a factor with levels corresponding to each trait.
• Notice this model is similar to the Multi-trait MET Model.

Example: Split-Plot Field Trial for Oats

Inference for fixed effects: Wald Test

fit <- asreml(yield ~ Variety*Nitrogen, random= ~ Blocks/Wplots, data=oats)
wald(fit) # same as anova(fit)

  Wald tests for fixed effects
  Response: yield
  Terms added sequentially; adjusted for those above
                   Df Sum of Sq Wald statistic Pr(Chisq)    
  (Intercept)       1     43410        245.141    <2e-16 ***
  Variety           2       526          2.971    0.2264    
  Nitrogen          3     20020        113.057    <2e-16 ***
  Variety:Nitrogen  6       322          1.817    0.9357    
  residual (MS)             177                             
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Kenward-Roger approximation to Wald Test

• Variance parameters are estimated thus the asymptotic results in the previous slide do not hold well.
• Kenward and Roger (1997) presented a scaled Wald statistic with an F approximation to its sampling distribution and an improved estimator for the precision of fixed effects (i.e. an adjusted variance matrix).
• The numerator degrees of freedom for each term is determined as the number of non-singular equations involved in the term.
• The calculation of the denominator degrees of freedom is in general not trivial and is computationally expensive.
• asreml uses a partial implementation of Kenward-Roger approach. More specifically, it uses an unadjusted variance matrix with variance parameters fixed at the REML estimates of the maximal model.

Conditional Tests

• By default, asreml show the incremental tests where the sum of the squares are adjusted for the terms listed above the rows.
• You can switch to conditional tests by using the argument ssType="conditional".

then Nelder (1977, 1994) argues that meaningful and interesting test for terms can only be conducted for those that respect marginality relations.

• For example, consider a factorial model

y ~ 1 + A + B + A:B

then A cannot be adjusted for A:B.

Conditional Tests in asreml

out_wald <- wald(fit, ssType="conditional", denDF="numeric")

out_wald$Wald %>% lucid::lucid()

                   Df denDF   F.inc   F.con Margin        Pr
  (Intercept)       1     5 245     245            0.0000193
  Variety           2    10   1.48    1.48       A 0.272    
  Nitrogen          3    45  37.7    37.7        A 0        
  Variety:Nitrogen  6    45   0.303   0.303      B 0.932

• Variety:Nitrogen is adjusted for other terms;
• Variety is adjusted for Nitrogen and (Intercept); and
• Nitrogen is adjusted for Variety and (Intercept).

Reisdual likelihood ratio test

data("burgueno.rowcol", package="agridat")
fit1 <- asreml(yield ~ gen, random=~rep, trace=F,
               data=burgueno.rowcol)
fit2 <- asreml(yield ~ gen, trace=F,
               data=burgueno.rowcol)
lrt(fit2, fit1)

  Likelihood ratio test(s) assuming nested random models.
  Chisq probability adjusted using Stram & Lee, 1994.
            Df LR statistic Pr(Chisq)    
  fit1/fit2  1       21.141 2.133e-06 ***
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

-2 * (fit2$loglik - fit1$loglik)

  [1] 21.1414

Prediction

fit <- asreml(Yield ~ Site + Geno, data=data1,
              random=~ Site:Geno + at(Site, "A"):Row)
predict(fit, classify="Geno")$predictions$pval

  Notes:
  - The predictions are obtained by averaging across the hypertable
    calculated from model terms constructed solely from factors in
    the averaging and classify sets.
  - Use "average" to move ignored factors into the averaging set.
  - The SIMPLE averaging set:  Site 
  - The ignored set:  Row 
    Geno predicted.value standard.error est.status
  1   G1        5.049309      0.4423711  Estimable
  2   G2        4.529920      0.4253100  Estimable
  3   G3        6.389828      0.4423711  Estimable

Understanding prediction calculation

• classify set: Geno
• averaging set: Site using uniform weights
• ignored set: Row

    Geno predicted.value standard.error est.status
  1   G1        5.049309      0.4423711  Estimable
  2   G2        4.529920      0.4253100  Estimable
  3   G3        6.389828      0.4423711  Estimable

Understanding prediction calculation

predict(fit, classify="Geno:Site")$predictions$pval

    Site Geno predicted.value standard.error est.status
  1    A   G1        5.346355      0.8971775  Estimable
  2    A   G2        4.883918      0.8803011  Estimable
  3    A   G3        6.622785      0.8971775  Estimable
  4    B   G1        4.442765      0.4949729  Estimable
  5    B   G2        3.649645      0.4855925  Estimable
  6    B   G3        5.641703      0.4949729  Estimable
  7    C   G1        5.358807      0.4949729  Estimable
  8    C   G2        5.056196      0.4855925  Estimable
  9    C   G3        6.904995      0.4949729  Estimable

More information on predict

?predict.asreml

And references:

• Lane and Nelder (1982) Analysis of covariance and standardization as instances of prediction. Biometrics 38(3) 613-621
• Gilmour et al. (2004) An efficient computing strategy for prediction in mixed linear models. Computational Statistics and Data Analysis 44(4) 571-586
• Welham et al. (2004) Prediction in linear mixed models. Australian & New Zealand Journal of Statistics 46(3) 325-347

Missing Values in asreml

• Missing values must be encoded as NA.
• By default, missing values in the explanatory variables cause an error: na.method.X="fail".
• If na.method.X="omit" then corresponding records with missing explanatory variables are dicarded.
• By default, missing values in the response are included and estimated from the model fit: na.method.Y="include". In this case a factor mv is created with one level for each missing response.
• If na.method.Y="omit" then missing values in the response are discarded.
• For na.method.X="include", missing values in explanatory variables are replaced with zeros. Thus it is essential that you can scale and center covariates with missing values.

Sparse Format of A^-1

asreml.Ainverse(ped)$ginv

The A Matrix

Convergence

LogLikelihood not converged

By default maximum iteration maxiter is 13 for the AI-REML algorithm. Convergence is determined by REML. You can run a few more iterations by:

asreml_object <- update(asreml_object)

or supplying the argument maxiter with a higher number, say 20

asreml(Yield ~ Site, random=~Geno, 
       maxiter=20)

Changing initial values for asreml

First get a table of initial values:

sv <- asreml(Yield ~ Site, data=data1, 
              random= ~Site:Geno,
              start.values=TRUE)
sv$gammas.table

                 Gamma Value Constraint
  1 Site:Geno!Site.var   0.1          P
  2         R!variance   1.0          P

Modify as you see fit:

sv$gammas.table$Value <- c(0.5, 1)
fit <- asreml(Yield ~ Site, data=data1, 
              random= ~Site:Geno,
              G.param=sv$gammas.table,
              R.param=sv$gammas.table)

Compound Symmetry Model

fit <- asreml(Yield ~ Site, data=data1, 
              random=~ Geno + Site:Geno)
lucid::vc(fit)

               effect component std.error z.ratio constr
        Geno!Geno.var 1.092     1.252        0.87      P
   Site:Geno!Site.var 0.0000004 0.0000002    2.5       B
           R!variance 0.9553    0.3747       2.5       P

sv <- asreml(Yield ~ Site, data=data1, 
    random= ~Geno + Site:Geno, start.values=TRUE) 
sv$gammas.table

                 Gamma Value Constraint
  1      Geno!Geno.var   0.1          P
  2 Site:Geno!Site.var   0.1          P
  3         R!variance   1.0          P

Allow for unconstrained estimate

fit <- asreml(Yield ~ Site, data=data1, 
              random= ~Geno + Site:Geno,
              G.param=sv$gammas.table, 
              R.param=sv$gammas.table) 
vc(fit)

               effect component std.error z.ratio constr
        Geno!Geno.var   1.099      1.255    0.88       P
   Site:Geno!Site.var  -0.03282    0.3954  -0.083      U
           R!variance   0.9755     0.4599   2.1        P

If you used the code "F" instead, it will fix the parameter to the specified initial value.