library(dplyr)
library(statforbiology)
#
# Loading data
dataset <- getAgroData("WheatQuality") |>
mutate(across(1:2, .fns = factor))
head(dataset)
## Genotype Block Height TKW Whectol Yield
## 1 arcobaleno 1 90 44.5 83.2 64.40
## 2 arcobaleno 2 90 42.8 82.2 60.58
## 3 arcobaleno 3 88 42.7 83.1 59.42
## 4 baio 1 80 40.6 81.8 51.93
## 5 baio 2 75 42.7 81.3 51.34
## 6 baio 3 76 41.1 81.1 47.78Imagine a field experiment comparing several treatments in a randomised complete block design. In each plot, we usually record several traits of interest, and we may want to investigate how they are associated with one another, without necessarily assuming a causal relationship. So, how can we carry out a sound correlation analysis?
In a previous post, I showed that simply calculating Pearson’s correlation coefficient from the raw plot data is not the best approach (follow this link, if you would like to read that post). In particular, I emphasised two issues.
- First, the observations are not independent: they are grouped by treatment and block. Ignoring this structure can invalidate standard tests and confidence intervals that rely on an independence assumption.
- Second, in a designed experiment, there is no single correlation coefficient that answers every possible question. We might be interested in (i) the correlation between plot measurements, (ii) the correlation between genotype effects, (iii) the correlation between block effects, or (iv) the correlation between residuals. These correlations have different interpretations, and focusing only on the first may lead to an oversimplified interpretation of the results.
In previous post I used the raw data, genotype means, block means and the residuals under simple univariable ANOVA models to calculate the Pearson’s correlation coefficients between yield and TKW.
However, this blog is about modelling, and I started wondering whether I could estimate all four types of correlation within a single modelling framework. An interesting approach was put forward by Piepho (2018) and it is based on multi-trait mixed models, which jointly consider two or more traits as responses and allow us to estimate the variances and covariances associated with different sources of variation.
For a randomised complete block design, we can write the model as:
\[ Y_{ijk} = \mu_i + \beta_{ik} + \tau_{ij} + \epsilon_{ijk}\]
where \(Y_{ijk}\) is the response for the trait \(i\), the treatment \(j\) and the block \(k\), \(\mu_i\) is the mean for the trait \(i\), \(\beta_{ik}\) is the effect of the block \(k\) and trait \(i\), \(\tau_{ij}\) is the effect of treatment \(j\) for the trait \(i\) and \(\epsilon_{ijk}\) is the residual for the trait \(i\), the treatment \(j\) and the block \(k\).
In this model, we make the following assumptions (Piepho, 2018):
- Trait-specific means. The trait effect is treated as fixed, with no intercept, so that each trait has its own mean, (\(\mu_i\)).
- Random genotype and block effects. The genotype and block effects (({ij}) and ({ik})) are treated as random. This is useful here because our main interest is in their variances and covariances, rather than in estimating each individual effect. For each source of variation, we allow the variance to differ between traits and estimate the covariance between traits. In mixed-model terminology, we use an unstructured variance-covariance matrix for the effects associated with each genotype and each block.
- Correlated residuals. The residuals are assumed to follow a multivariate normal distribution, with trait-specific variances and a covariance between traits measured on the same plot. Thus, the residual variance-covariance matrix is unstructured within each plot.
These assumptions allow us to estimate the correlations between genotype effects, between block effects, and between residuals, directly from model fit. Considering that the genotype, block, and residual components are mutually independent, we can also derive the marginal correlation between plot measurements by combining the corresponding variance and covariance components.
Of course, in other applications, genotype effects (and sometimes block effects) might be better treated as fixed. The appropriate choice depends on the research question and the population to which we want to generalise. We will not pursue that alternative here.
Fitting a multi-trait model in R is not entirely straightforward. For this example, I find the ‘sommer’ package (Covarrubias-Pazaran, 2016) particularly useful, because it provides a flexible way to specify mixed models with a range of variance-covariance structures.
Fitting a bivariate mixed model with ‘sommer’
Let’s consider the WheatQuality dataset, available in the statforbiology package. It contains data from a genotype experiment with 27 genotypes arranged in three complete blocks. Suppose that we measured two traits in every plot: grain yield and thousand-kernel weight (TKW). How are these two traits correlated, and does the answer depend on which source of variation we consider?
We can fit a bivariate mixed model using ‘sommer’. The code below is intended for version 4.4.98. Because the package has undergone some recent changes, it may not work with earlier versions.
There are four main steps:
- Reshape the data. We need to convert the dataset to long format, stacking the two traits in a single response column. We also need a factor identifying the original plots, because each plot contributes two observations, one for each trait. The
stackTraits()function in ‘sommer’ does this and creates three variables: ‘record’, identifying the plots; ‘trait’, identifying the traits; and ‘value’, containing the observed responses. - Specify the fixed effects. We use the usual R formula
value ~ trait - 1. Removing the intercept allows us to estimate a separate mean for each trait. - Specify the random effects. In sommer, the general structure
vsm(usm(StructuringFactor), ism(GroupingFactor))allows us to specify the variance-covariance structure of the random effects associated with the levels of ‘GroupingFactor’, using ‘StructuringFactor’ to define the covariance structure across traits. Here,usm()specifies an unstructured covariance matrix (other functions are available in ‘sommer’ to model different structures), whileism()identifies the grouping structure. We use this specification twice: once for genotypes and once for blocks. In both cases, we allow the two traits to have different variances and a non-zero covariance. - Specify the residual structure. Finally, we allow the residuals for the two traits measured on the same plot to be correlated, using
rcov = ~ vsm(usm(trait), ism(record)).
The final model call is reported in the box below:
library(sommer)
# Put the dataset in long format
mdataset <- dataset |>
sommer::stackTraits(c("Yield", "TKW"), keep = 1:2) |>
dplyr::mutate(across(c(1,2), .fns = factor))
head(mdataset)
## Genotype Block record trait value
## 1 arcobaleno 1 r1 Yield 64.40
## 2 arcobaleno 2 r2 Yield 60.58
## 3 arcobaleno 3 r3 Yield 59.42
## 4 baio 1 r4 Yield 51.93
## 5 baio 2 r5 Yield 51.34
## 6 baio 3 r6 Yield 47.78
# Fit the model
mod.bimix <- mmes(value ~ trait - 1,
random = ~ vsm(usm(trait), ism(Genotype)) +
vsm(usm(trait), ism(Block)),
rcov = ~vsm(usm(trait), ism(record)),
data = mdataset)
## BLAS backend: /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRblas.0.dylib.
## BLAS threads configured: unknown (unavailable); allocation-aware CPUs: unknown.
## MME path: AI-first with on-demand EM/PX-GEM recovery [solver=cholmod; latent RR=inactive; Schur=active; chain=inactive; banded=inactive; factor-Schur=inactive; residual=kronecker; grid=inactive; repeated-Kronecker=active]
## Dense block-Schur engine active (27 random-effect blocks, 8 border effects).
## iter LogLik deltaLL update policy trust line notPD time(s)
## 1 -17.135435 - AI auto:em 0.168 1.000 - 0.01
## 2 31.456688 48.5921 AI auto:em 0.648 1.000 - 0.01
## 3 52.818502 21.3618 AI auto:em 1.000 1.000 - 0.01
## 4 53.781853 0.9634 AI auto:em 1.000 1.000 - 0.01
## 5 53.800806 0.0190 AI auto:em 1.000 1.000 - 0.01
## 6 53.800818 0.0000 AI auto:em - - - 0.01
## Fit CPU usage: 0.072 CPU-s / 0.080 elapsed-s = 0.90 average CPUs (not an active-thread count).
## blockSchur usage (6 calls, block extraction, factorization, Schur assembly and log determinant): 0.000 CPU-s / 0.000 elapsed-s = 0.98 average CPUs.
bimix.obj <- summary(mod.bimix)
coefs <- bimix.obj$varcomp
coefs[,-5]
## term parameter estimate StdError
## 1 vsm(usm(trait), ism(Genotype)) variance[Yield] 77.6345915 22.0989241
## 2 vsm(usm(trait), ism(Genotype)) covariance[TKW,Yield] 38.8322623 15.0937768
## 3 vsm(usm(trait), ism(Genotype)) variance[TKW] 53.8615917 15.3547154
## 4 vsm(usm(trait), ism(Block)) variance[Yield] 3.7104778 3.9364268
## 5 vsm(usm(trait), ism(Block)) covariance[TKW,Yield] -0.2428390 1.9074544
## 6 vsm(usm(trait), ism(Block)) variance[TKW] 1.6684568 1.8343662
## 7 vsm(usm(trait), ism(record)) variance[Yield] 6.0939037 1.1951128
## 8 vsm(usm(trait), ism(record)) covariance[TKW,Yield] 0.1635551 0.7242690
## 9 vsm(usm(trait), ism(record)) variance[TKW] 4.4717901 0.8769902The model output contains the estimated variance-covariance parameters. We can use these estimates to calculate the corresponding correlations. Recall that a correlation is the covariance divided by the product of the two standard deviations.
The vpredict() function in sommer allows us to calculate these quantities and obtain their standard errors. In the following code, V1, V2, and so on refer to the variance-covariance parameters in the order in which they appear in the model output. It is important to check this ordering against the output from your fitted model before interpreting the results.
# Correlation between genotype effects
vpredict(mod.bimix, rg ~ V2 / (sqrt(V1)*sqrt(V3) ))
## Estimate SE
## rg 0.6005172 0.1306758
# Correlation between block effects
vpredict(mod.bimix, rt ~ V5 / (sqrt(V4)*sqrt(V6) ))
## Estimate SE
## rt -0.09759916 0.7571335
# Correlation of residuals under the bivariate model
vpredict(mod.bimix, re ~ V8 / (sqrt(V7)*sqrt(V9) ))
## Estimate SE
## re 0.03133109 0.1385389We can also calculate the marginal correlation between plot measurements. Under this model, the marginal variance of each trait is the sum of its genotype, block, and residual variances. Similarly, the marginal covariance between the two traits is the sum of the corresponding covariances. We can therefore calculate the plot-level correlation as follows:
# Correlation of plot values
vpredict(mod.bimix, rp ~ (V2 + V5 + V8) / (sqrt(V1 + V4 + V7)*sqrt(V3 + V6 + V9) ))
## Estimate SE
## rp 0.5350203 0.1281157We can also use covmatrix_mmes() to obtain the estimated correlation matrices and their standard errors.
# Correlation between genotype effects (with standard errors)
covmatrix_mmes(mod.bimix, term = 1)$correlation
## Yield TKW
## Yield 1.0000000 0.6005172
## TKW 0.6005172 1.0000000
covmatrix_mmes(mod.bimix, term = 1)$correlation.se
## Yield TKW
## Yield 0.0000000 0.1306758
## TKW 0.1306758 0.0000000
# Correlation between block effects (with standard errors)
covmatrix_mmes(mod.bimix, term = 2)$correlation
## Yield TKW
## Yield 1.00000000 -0.09759916
## TKW -0.09759916 1.00000000
covmatrix_mmes(mod.bimix, term = 2)$correlation.se
## Yield TKW
## Yield 0.0000000 0.7571335
## TKW 0.7571335 0.0000000
# Correlation between the residuals (with standard errors)
covmatrix_mmes(mod.bimix, term = 3)$correlation
## Yield TKW
## Yield 1.00000000 0.03133109
## TKW 0.03133109 1.00000000
covmatrix_mmes(mod.bimix, term = 3)$correlation.se
## Yield TKW
## Yield 0.0000000 0.1385389
## TKW 0.1385389 0.0000000The estimated correlations are quite close to the Pearson correlation coefficients obtained from the raw data (see my previous post). However, the main advantage of the mixed-model approach is that it separates the different sources of variation and allows us to investigate whether the corresponding covariance components differ from zero.
For example, we can investigate whether the residuals for the two traits are correlated by fitting a reduced model in which the residual variances are trait-specific but the residual covariance is constrained to zero. We can then compare this model with the full model using a likelihood-ratio statistic. In ‘sommer’, we can impose such a restriction by replacing usm() with dsm(), which specifies a diagonal variance-covariance structure. The code is shown below.
One important caveat: likelihood-ratio tests for variance-covariance parameters require care. Testing whether a variance component is zero involves a boundary condition, whereas testing whether a covariance is zero, with both variances positive, is generally an interior test. Nevertheless, the validity of the usual chi-squared reference distribution depends on the model and the assumptions underlying the test. It is therefore important to check how the test is implemented and interpreted for the fitted models.
mod.bimix2 <- mmes(value ~ trait - 1,
random = ~ vsm(usm(trait), ism(Genotype)) +
vsm(usm(trait), ism(Block)),
rcov = ~vsm(dsm(trait), ism(record)),
data = mdataset)
## BLAS backend: /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRblas.0.dylib.
## BLAS threads configured: unknown (unavailable); allocation-aware CPUs: unknown.
## MME path: AI-first with on-demand EM/PX-GEM recovery [solver=cholmod; latent RR=inactive; Schur=active; chain=inactive; banded=inactive; factor-Schur=inactive; residual=diagonal; grid=inactive; repeated-Kronecker=inactive]
## Dense block-Schur engine active (27 random-effect blocks, 8 border effects).
## iter LogLik deltaLL update policy trust line notPD time(s)
## 1 -18.457629 - AI auto:em 0.169 1.000 - 0.00
## 2 29.664739 48.1224 AI auto:em 0.648 1.000 - 0.00
## 3 52.370291 22.7056 AI auto:em 1.000 1.000 - 0.00
## 4 53.735832 1.3655 AI auto:em 1.000 1.000 - 0.00
## 5 53.775232 0.0394 AI auto:em 1.000 1.000 - 0.00
## 6 53.775283 0.0001 AI auto:em - - - 0.00
## Fit CPU usage: 0.015 CPU-s / 0.015 elapsed-s = 1.00 average CPUs (not an active-thread count).
## blockSchur usage (6 calls, block extraction, factorization, Schur assembly and log determinant): 0.000 CPU-s / 0.000 elapsed-s = 0.95 average CPUs.
# Likelihood ratio test
LRT1 <- anova(mod.bimix, mod.bimix2)
## Likelihood ratio test for mixed models
## ==============================================================
## Df AIC BIC loLik Chisq ChiDf PrChisq
## mod1 11 -103.6016 -97.42644 53.80082
## mod2 10 -103.5506 -97.37537 53.77528 0.05107 1 0.82121
## ==============================================================
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Similarly, we can impose restrictions on the variance-covariance matrix for random effects, to test for the hypothesis that the correlation between genotype/block effects is zero.
# Removing correlation between block effects
mod.bimix3 <- mmes(value ~ trait - 1,
random = ~ vsm(usm(trait), ism(Genotype)) +
vsm(dsm(trait), ism(Block)),
rcov = ~vsm(usm(trait), ism(record)),
data = mdataset)
## BLAS backend: /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRblas.0.dylib.
## BLAS threads configured: unknown (unavailable); allocation-aware CPUs: unknown.
## MME path: AI-first with on-demand EM/PX-GEM recovery [solver=cholmod; latent RR=inactive; Schur=active; chain=inactive; banded=inactive; factor-Schur=inactive; residual=kronecker; grid=inactive; repeated-Kronecker=active]
## Dense block-Schur engine active (27 random-effect blocks, 8 border effects).
## iter LogLik deltaLL update policy trust line notPD time(s)
## 1 -17.140136 - AI auto:em 0.168 1.000 - 0.00
## 2 31.427841 48.5680 AI auto:em 0.627 1.000 - 0.00
## 3 52.839209 21.4114 AI auto:em 1.000 1.000 - 0.00
## 4 53.771138 0.9319 AI auto:em 1.000 1.000 - 0.00
## 5 53.792606 0.0215 AI auto:em 1.000 1.000 - 0.00
## 6 53.792621 0.0000 AI auto:em - - - 0.00
## Fit CPU usage: 0.017 CPU-s / 0.017 elapsed-s = 1.00 average CPUs (not an active-thread count).
## blockSchur usage (6 calls, block extraction, factorization, Schur assembly and log determinant): 0.000 CPU-s / 0.000 elapsed-s = 0.97 average CPUs.
LRT2 <- anova(mod.bimix, mod.bimix3)
## Likelihood ratio test for mixed models
## ==============================================================
## Df AIC BIC loLik Chisq ChiDf PrChisq
## mod1 11 -103.6016 -97.42644 53.80082
## mod2 10 -103.5852 -97.41005 53.79262 0.01639 1 0.89812
## ==============================================================
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Removing correlation between genotype effects
mod.bimix4 <- mmes(value ~ trait - 1,
random = ~ vsm(dsm(trait), ism(Genotype)) +
vsm(usm(trait), ism(Block)),
rcov = ~vsm(usm(trait), ism(record)),
data = mdataset)
## BLAS backend: /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRblas.0.dylib.
## BLAS threads configured: unknown (unavailable); allocation-aware CPUs: unknown.
## MME path: AI-first with on-demand EM/PX-GEM recovery [solver=cholmod; latent RR=inactive; Schur=inactive (auto: not selected); chain=inactive; banded=inactive; factor-Schur=inactive; residual=kronecker; grid=inactive; repeated-Kronecker=active]
## iter LogLik deltaLL update policy trust line notPD time(s)
## 1 -18.855128 - AI auto:em 0.176 1.000 - 0.00
## 2 28.611386 47.4665 AI auto:em 0.643 1.000 - 0.00
## 3 48.021734 19.4103 AI auto:em 1.000 1.000 - 0.01
## 4 48.354609 0.3329 AI auto:em 1.000 1.000 - 0.01
## 5 48.355448 0.0008 AI auto:em - - - 0.01
## Fit CPU usage: 0.016 CPU-s / 0.018 elapsed-s = 0.89 average CPUs (not an active-thread count).
## cholmod usage (5 calls, numeric factorization): 0.000 CPU-s / 0.000 elapsed-s = 0.37 average CPUs.
# Likelihood ratio test
LRT3 <- anova(mod.bimix, mod.bimix4)
## Likelihood ratio test for mixed models
## ==============================================================
## Df AIC BIC loLik Chisq ChiDf PrChisq
## mod1 11 -103.6016 -97.42644 53.80082
## mod2 10 -92.7109 -86.53570 48.35545 10.89074 1 0.00097 ***
## ==============================================================
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1In this example, the genotype effects show evidence of a correlation between the two traits, whereas the correlations between block effects and between residuals do not appear to be statistically significant under the tests considered. This does not necessarily mean that the latter correlations are zero; rather, the data did not provide sufficient evidence that they differ from zero.
This illustrates why it can be useful to distinguish among the different sources of correlation in a designed experiment. A single correlation coefficient calculated from the raw data cannot tell us whether the association arises mainly from differences between genotypes, differences between blocks, or the remaining plot-level variation.
I hope this example has been useful. If you know of a better solution, or if you have suggestions for improving this approach, I’d be happy to hear from you. Please drop me a line at the address below.
(P.S. If you enjoyed this post, check out my book “Field Research Methods in Agriculture: An Introduction with R” via the link below!)
Prof. Andrea Onofri
Department of Agricultural, Food and Environmental Sciences
University of Perugia (Italy)
Send comments to: andrea.onofri@unipg.it
References
- Butler, D., Cullis, B.R., Gilmour, A., Gogel, B., Thomson, R., 2018. ASReml-r reference manual, version 4. UK.
- Covarrubias-Pazaran, G., 2016. Genome-Assisted Prediction of Quantitative Traits Using the R Package sommer. PLOS ONE 11, e0156744. https://doi.org/10.1371/journal.pone.0156744
- Piepho, H.-P., 2018. Allowing for the structure of a designed experiment when estimating and testing trait correlations. The Journal of Agricultural Science 156, 59–70.
- Pinheiro, J.C., Bates, D.M., 2000. Mixed-effects models in s and s-plus. Springer-Verlag Inc., New York.
Appendix
Is it possible to fit the above models with other packages in R? In a previous version of this post I have used the ‘asreml’ package (Butler et al., 2018), although this is a paid option. However, here is the relevant code:
mod.asreml <- asreml(cbind(Yield, TKW) ~ trait - 1,
random = ~ Genotype:us(trait, init = c(3.7, -0.25, 1.7)) +
Block:us(trait, init = c(77, 38, 53)),
residual = ~ units:us(trait, init = c(6, 0.16, 4.5)),
data=dataset)
#summary(mod.asreml)$varcomp[,1:3]
## component std.error z.ratio
## Block:trait!trait_Yield:Yield 3.7104778 3.9364268 0.9426005
## Block:trait!trait_TKW:Yield -0.2428390 1.9074544 -0.1273105
## Block:trait!trait_TKW:TKW 1.6684568 1.8343662 0.9095549
## Genotype:trait!trait_Yield:Yield 77.6346623 22.0956257 3.5135761
## Genotype:trait!trait_TKW:Yield 38.8322972 15.0909109 2.5732242
## Genotype:trait!trait_TKW:TKW 53.8616088 15.3520661 3.5084274
## units:trait!R 1.0000000 NA NA
## units:trait!trait_Yield:Yield 6.0939037 1.1951128 5.0990195
## units:trait!trait_TKW:Yield 0.1635551 0.7242690 0.2258209
## units:trait!trait_TKW:TKW 4.4717901 0.8769902 5.0990195
# Reduced model
mod.asreml2 <- asreml(cbind(Yield, TKW) ~ trait - 1,
random = ~ Genotype:us(trait) + Block:us(trait),
residual = ~ units:idh(trait),
data=dataset)
# LRT
lrt.asreml(mod.asreml, mod.asreml2) To the best of my knowledge, there is no way to fit such a complex model with the ‘nlme’ package and related lme() function. However, if we are willing to treat the ‘block’ effect as fixed, we can fit a bivariate mixed model with the lme() function (Pinheiro and Bates, 2000). This alternative retains correlated genotype effects and allows the residual variances to differ between traits, with a correlation between the two observations from the same plot. The code is shown below.
library(nlme)
mod <- lme(value ~ trait*Block,
random = list(Genotype = pdSymm(~trait - 1)),
weight = varIdent(form=~1|trait),
correlation = corCompSymm(form=~1|Genotype/record),
data = mdataset)