The model call is: mod.asreml <- asreml(cbind(Yield, TKW) ~ trait - 1, The residual structure, where the observations in the same plot (the term ‘units’ is used in ‘asreml’ to represent the observational units, i.e. the plots) are heteroscedastic and correlated.Therefore, the random model is ~ Genotype:us(trait) + Block:us(trait). For both, I specified a general unstructured variance covariance matrix, so that the traits are heteroscedastic and correlated. The random model, that is composed by the interactions ‘genotype x trait’ and ‘block x trait’.Following Piepho (2018), I removed the intercept, to separately estimate the means for the two traits. The specification is, therefore, ‘cbind(Yield, TKW) ~ trait - 1’.
The fixed model, that only contains the trait effect.When we set a bivariate model with ‘asreml’, we can ‘cbind()’ Yield and TKW and use the name ‘trait’ to refer to them. With the function ‘asreml()’, we need to specify the following components. Therefore, I decided to use the package ‘asreml’ (Butler et al., 2018), although this is not freeware. To the best of my knowledge, there is no way to fit such a complex model with the ‘nlme’ package and related ‘lme()’ function (I’ll gave a hint later on, for a simpler model). Both random effects need to be heteroscedastic (trait-specific variance components) and there must be a correlation between the two traits.
This idea makes sense, because, for this application, we are mainly interested in variances and covariances. Hans-Peter Piepho, in his paper, put forward the idea that the ‘genotype’ and ‘block’ effects for the two traits can be taken as random, even though they might be of fixed nature, especially the genotype effect. Furthermore, residuals belonging to the same plot (the two observed traits) need to be correlated (correlation of errors). \[ Y_\) need to be normally distributed and heteroscedastic, with trait-specific variances. Piepho (2018) showed that, for an experiment like this one, the above correlations can be estimated by coding a multiresponse mixed model, as follows:
In my previous post, I used the above dataset to calculate the Pearson’s correlation coefficient between Yield and TKW for: Head(dataset) # Genotype Block Height TKW Yield Whectol The code below loads the necessary packages, the data and transforms the numeric variable ‘Block’ into a factor. The dataset ‘WheatQuality.csv’ is available on ‘gitHub’ it consists of 81 records (plots) and just as many couples of measures in all. For each plot, the collegue who made the experiment recorded several traits, including yield (Yield) and weight of thousand kernels (TKW).
It is a genotype experiment, laid out in randomised complete blocks, with 27 wheat genotypes and three replicates. I will use the same example as presented in my previous post, which should allow us to compare the results with those obtained by using more traditional methods of data analyses.
0 kommentar(er)
