ASRemlR Snippets: Constraining Variance Components
Often, after fitting a complex linear mixed model (LMM), we are interested in constraining some of the variance components (VCs). There are many reasons why we might want to do this, such as when we combine data from several sources and we want to have a common parameter, or when we want to simplify a model to evaluate a specific hypothesis relating to these components.
In this snippet we will briefly explore an example where we constrain two VCs to be identical. The example comes from Becker (1984) and corresponds to data from 28 white pine (Pinus monticola) families formed by crossing 4 male to 7 female parents. The progeny was evaluated in a randomized complete block design with 4 replicates and planted in large plots. Our response of interest was epicotyl length, the part of the embryonic seedling stem from which the plant’s shoot system develops.
Using ASRemlR, the reference model \((\texttt{mod0})\) fitted to this data is:

In the above model we have a fixed effect of replicate (\(\texttt{rep}\)), and random effects for the parents and their interaction (\(\texttt{female}\), \(\texttt{male}\) and \(\texttt{female:male}\)). This LMM produces the following output in ASRemlR:
component  std.error  z.ratio  bound  %ch  

male  0.1203614  0.1203723  0.9999094  P  0 
female  0.0209448  0.0420690  0.4978693  P  0 
female:male  0.1368513  0.0628149  2.1786443  P  0 
units!R  0.2004263  0.0314941  6.3639416  P  0 
In this example, the variance of males (\(\hat{{\sigma_m}}^2\) = 0.120) is approximately six times larger than the variance for females (\(\hat{{\sigma_f}}^2\) = 0.021). Therefore, the additive variance from males and females are contrastingly different: \(\texttt{V}_{Am} = 4 \space \sigma_m^2 \neq 4 \space \sigma_f^2 = \texttt{V}_{Af}\).
The above results may seem biologically unreasonable for pines, where the genetic contribution of male and female parents (i.e., \(\texttt{V}_{A}\)) is often very similar. Hence, we expect their genetic variances to be almost identical. That is, we think \(\sigma_f^2 = \sigma_m^2\). The inconsistent VC estimates from our reference model could be the result of random sampling of families, or to the limited number of parents considered in this dataset (recall that we have only 4 male and 7 female parents).
In ASRemlR, the constraint \(\sigma_f^2 = \sigma_m^2\) is easy to implement. This is done by creating a matrix, say \(\texttt{M}\), that defines the grouping of the VCs. All VCs in the same group (i.e., with the same number) will be constrained to be identical. For our problem, the matrix looks like:

The first column, \(\texttt{V1}\), indicates the grouping; hence, there are three VCs two of which (\(\texttt{female}\) and \(\texttt{male}\)) will be identical. The second column, \(\texttt{V2}\), is used to specify the multiplicative relationships between components, but this is a topic of another ASRemlR snippet.
The construction of the above matrix \(\texttt{M}\) can be facilitated with the following R code. This starts by requesting the table of VC parameters (\(\texttt{gam}\)) from ASRemlR:

It then proceeds to generate our matrix \(\texttt{M}\) of constraints using only the column names from the table \(\texttt{gam}\):

Now, we are ready to fit the model in ASRemlR with the constraint \(\sigma_f^2 = \sigma_m^2\). The matrix \(\texttt{M}\) is incorporated using the option \(\texttt{vcc=M}\):

The output from fitting this constrained model in ASRemlR is:
component  std.error  z.ratio  bound  %ch  

male  0.0639041  NA  NA  C  0.5 
female  0.0639041  0.0369027  1.731692  P  0.5 
female:male  0.1327523  0.0602139  2.204676  P  0.0 
units!R  0.2005090  0.0312342  6.419535  P  0.0 
In the above table see that \(\hat{{\sigma_m}}^2 = \hat{{\sigma_m}}^2\) = 0.064. Also note that the zratio for this variance component (1.731) is much higher than the zratios for both \(\texttt{male}\) (1.000) or \(\texttt{female}\) (0.498) obtained from our reference model (from the first table), indicating better use of the available information.
We compare our constrained model (\(\texttt{modC}\)) against the reference model (\(\texttt{mod0}\)) statistically, using a likelihood ratio test (LRT). This is implemented using the function (\(\texttt{lrt.asreml}\)):

df  LRstatistic  Pr(Chisq)  

mod0/modC  1  1.03068  0.309999 
The LRT indicates that our null hypothesis of \(\texttt{H}_0: \sigma_f^2 = \sigma_m^2\) cannot be rejected in favor of our alternative hypothesis \(\texttt{H}_1: \sigma_f^2 \neq \sigma_m^2\) (pvalue = 0.31). That is, despite the estimated variance for \(\texttt{males}\) being approximatelly six times larger than the variance for \(\texttt{females}\) under \(\texttt{mod0}\), there is no sufficient statistical evidence that they are different. Probably, this result is due to the limited number of parents and crosses considered in this study (i.e., low statistical power).
And as a final extension, for the constrained analysis we also can calculate a narrowsense heritability using the function \(\texttt{vpredict}\), as:

Estimate  SE  

h2  0.5543991  0.2486712 
The above output provides us with an estimate of \(\texttt{h}^2\), 0.554 (SE = 0.249), a reasonably high heritability estimate but with a large amount of uncertainty.
The appeal of this ASRemlR approach for constraining variance components is that, for any complex linear mixed models, constraints can be easily specified in column \(\texttt{V1}\) of matrix \(\texttt{M}\), for most hypothesis of interest, or to evaluate models with specific biological or mathematical constraints.
Author
Salvador A. Gezan
References
Becker (1984). Manual of Quantitative Genetics, 4th Edition. Academic Enterprises, Pullman, WA.
Files to download
Notes: SAG May2020