Biology 708, McMaster University

Multivariate modeling via mixed models

This tutorial follows from our introduction to multivariate linear models, extending it by using multivariate linear mixed models.

Useful packages:

library(MCMCglmm)
library(lme4)
library(brms)
library(tidyr)
library(dplyr)
library(corrplot)
library(broom.mixed)
library(dotwhisker)
library(ggplot2); theme_set(theme_bw())

data

This is an old Drosophila melanogaster data set from Ian Dworkin’s PhD work (Dworkin 2005). This was from a study that was meant to test predictions of a model on how mutational and environmental variation can influence the overall structure of phenotypic variation. For this study Dworkin measured several traits (lengths) on the first leg as well as the number of sex comb teeth (a structure used to clasp females during copulation) for different wild type strains (line) reared at different developmental temperatures (temp), with and without a mutation that effects proximal-distal axis development in limbs (genotype).

Get data from https://doi.org/10.5061/dryad.8376 or this repo:

dll_data <- read.csv("../data/dll.csv")

## make temp a factor (25 vs 30 degrees)
dll_data$temp <- factor(dll_data$temp)
## scale relevant variables (fancier but less repetition than previously)
morph_vars <- c("femur","tibia","tarsus","SCT")
morph_vars_sc <- paste(morph_vars,"s",sep="_")
dll_data2 <- dll_data
## c() drops unwanted structure from the results of scale()
for (i in 1:length(morph_vars)) {
    dll_data2[[morph_vars_sc[i]]] <- c(scale(dll_data[[morph_vars[i]]]))
}

Mixed models

The previous expression for the multivariate model was

\[ \mathbf{Y} = \mathbf{XB} + \mathbf{E} \]

In contrast, the expression for the mixed model is

\[ y = \mathbf{X\beta} + \mathbf{Zb} + \epsilon \]

where \(\mathbf{b}\) is a set of Gaussian variables with a variance-covariance matrix \(\mathbf{\Sigma}\) which we estimate.

---

Suppose we have observations of \(m\) individuals, with \(p\) different observations (traits, or time points, or types of measurement, or …) for each individual. The way we’re going to make this work is to expand (“melt”, or gather() or pivot_longer() if we’re using tidyr) the data set to be \(mp\) observations long, then treat each individual (which was previously a single row of the data set but now comprises \(p\) rows) as a group (we’ll call this units):

A trick to do multivariate mixed models using lme4

lme4 does not (currently) have a natural syntax for multivariate responses, however, as I eluded to in class, there is an important relationship between multivariate response models and so called “repeated” measures (or longitudinal) models. As such we can use a few tricks to estimate the model in lme4. Below this, I will go through the same model using MCMCglmm, which is a library which has a more natural syntax for such multivariate responses, but is explicitly Bayesian, so you need to provide prior distributions.

melting code for lme4

What we need to first do (for lme4) is to generate a single column that represents the numeric values for our response traits, and then a second variable that stores the trait type.

dll_melt <- (dll_data2
    %>% select(-c(femur,tibia,tarsus,SCT))  ## drop unscaled vars
    %>% mutate(units=factor(1:n()))
    %>% gather(trait,value, -c(units,replicate,line,genotype,temp))
    %>% drop_na()
    %>% arrange(units)
)
## saveRDS(dll_melt, file="dll_melt.rds")

(in case you’re familiar with current-day tidyverse: gather() is an earlier version of pivot_longer(), which is now recommended).

---

Look at how this has changed the structure of the data from a wide format:

head(dll_data) # original wide

##   replicate   line genotype temp femur tibia tarsus SCT
## 1         1 line-1      Dll   25 0.590 0.499  0.219   9
## 2         1 line-1      Dll   25 0.550 0.501  0.214  13
## 3         1 line-1      Dll   25 0.588 0.488  0.211  11
## 4         1 line-1      Dll   25 0.588 0.515  0.211  NA
## 5         1 line-1      Dll   25 0.596 0.502  0.207  12
## 6         1 line-1      Dll   25 0.577 0.499  0.207  14

to a long format:

head(dll_melt)

##   replicate   line genotype temp units    trait  value
## 1         1 line-1      Dll   25     1  femur_s  1.514
## 2         1 line-1      Dll   25     1  tibia_s  0.597
## 3         1 line-1      Dll   25     1 tarsus_s  1.751
## 4         1 line-1      Dll   25     1    SCT_s -1.222
## 5         1 line-1      Dll   25     2  femur_s  0.138
## 6         1 line-1      Dll   25     2  tibia_s  0.654

lmer fit

We can now go ahead and fit the model where we need to include trait as a predictor variable (where each trait is now a “repeated measure” from a particular subject/individual)

t1 <- system.time(
    lmer1 <- lmer(value ~ trait:(genotype*temp) - 1 +
                      (trait-1|line) + (trait-1|units),
                  data=dll_melt,
                  control=lmerControl(
                      optimizer="bobyqa",
                      check.nobs.vs.nlev="ignore",
                      check.nobs.vs.nRE="ignore"))
)

lme4: fixed-effects formula

• the system.time() business is just to keep track of how long the command takes to run (we could have kept the HTML output clean by hiding this in the .Rmd file, at the cost of making the .Rmd more complicated/confusing for those who were following along there)
• it doesn’t make sense to consider effects that apply equally to all traits
• so, let trait interact with all the other variables, but nothing else - use -1 to drop intercept
• specification is equivalent to
  • trait:(1+genotype+temp+genotype:temp)
  • or trait + trait:genotype + trait:temp + trait:genotype:temp

lme4: random-effects formula

• (trait-1|line) means “variance/covariances of traits among lines”
• -1 so we consider traits, not differences among traits
• (trait-1|units) specifies “var/cov of traits among units (individuals)”

lme4: control specificatons

• lmer always includes a residual variance term. This is redundant because we have only one data point per individual per trait: lmerControl(...) tells lmer not complain
• As a result, our within-individual correlation estimates will be slightly overestimated (not so for GLMMs)
• optimizer="bobyqa" is the only choice of optimizers that does not give a convergence warning (ask about this if you feel like going down a rabbit hole)

Notes on lme4 results

• the fit is a little slow (about 10 seconds on a fast machine: GLMMs would be even slower)
• print(lmer1), summary(lmer1) useful but awkward (16 fixed effect coefficients, we’ll look at graphical summaries instead
• fixef() (fixed effects), ranef() (line/indiv deviations from pop mean), coef() (line/indiv-level estimates), VarCorr (random effects var/cov)

lme4: random-effects var/cov

• underlying var-cov parameters are labeled \(\theta\): getME(fit,"theta") extracts them if you need them for something.

isSingular(lmer1)

## [1] FALSE

VarCorr(lmer1)

##  Groups   Name          Std.Dev. Corr           
##  units    traitfemur_s  0.652                   
##           traitSCT_s    0.726    0.06           
##           traittarsus_s 0.563    0.58 0.21      
##           traittibia_s  0.668    0.92 0.12 0.48 
##  line     traitfemur_s  0.489                   
##           traitSCT_s    0.392    0.14           
##           traittarsus_s 0.462    0.81 0.37      
##           traittibia_s  0.473    0.87 0.22 0.83 
##  Residual               0.465

lme4: correlation plot

par(mfrow=c(1,2))
vv1 <- VarCorr(lmer1)
## fix unit variance-covariance by adding residual variance:
diag(vv1$units) <- diag(vv1$units)+sigma(lmer1)^2
corrplot.mixed(cov2cor(vv1$line),upper="ellipse")
corrplot.mixed(cov2cor(vv1$units),upper="ellipse")

Correlations of traits across lines are stronger than correlations within individuals. In both cases correlations are all positive (i.e. first axis of variation is size variation?)

lme4: coefficient plot

cc1 <- tidy(lmer1,effect="fixed") %>%
  tidyr::separate(term,into=c("trait","fixeff"),extra="merge",
                  remove=FALSE)
dwplot(cc1)+facet_wrap(~fixeff,scale="free",ncol=2)+
  geom_vline(xintercept=0,lty=2)

## Warning: Using the `size` aesthetic with geom_segment was deprecated in ggplot2 3.4.0.
## ℹ Please use the `linewidth` aesthetic instead.
## ℹ The deprecated feature was likely used in the dotwhisker package.
##   Please report the issue at <https://github.com/fsolt/dotwhisker/issues>.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.

These results tell approximately the same story (coefficients are consistent across traits within a fixed effect, e.g. effect of higher temperature is to reduce scores on all traits).

lme4: random effects

This is very slow, you probably don’t want to evaluate it on your computer

t1_CI <- system.time(cc2 <- tidy(lmer1,
            effect = "ran_pars",
            conf.int = TRUE,
            conf.method = "profile",
            parallel="multicore",
            ncpus=8))
saveRDS(cc2, file="../data/cc2.rds")

(about 40 minutes to run on 8 cores on a fast machine)

MCMCglmm

Another option is the MCMCglmm package, which has a natural interface for general multivariate mixed models. It takes a while to get used to the interface, but here is an example. Please check out here for more information about the package, and here for an overview of how to use it.

---

First I find it easier (given the interface of MCMCglmm) to create a formula with the response variables and predictors. This is only for the fixed effects part of the model.

fmla.MMLM1  <- cbind(femur_s, tibia_s, tarsus_s, SCT_s) ~
    trait:(genotype*temp) - 1

---

Now we need to let MCMCglmm know which family (i.e. distribution) the response variables are. Since all are normal (Gaussian), we can specify it the following way.

fam.test <- rep("gaussian", 4 )

---

Since MCMCglmm is fundamentally a Bayesian approach, it needs a prior. If you provide no prior by default, it tries a “flat” prior, although this rarely works. In this case we provide not-quite-flat priors only for the random effects of line and for the residual variances (could also provide them for the fixed effects). Choosing priors for variances and covariances is tricky (and too technical for right now), so we will use the default prior distributions (inverse-Wishart). Ian Dworkin has some simple tutorials where he draws out the effects of various common priors for variances and covariances.

prior.model.1 <- list( R = list(V=diag(4)/4, nu=0.004),  
                       G = list(G1=list(V=diag(4)/4, nu=0.004)))

---

Let’s take a quick look at the prior

prior.model.1

## $R
## $R$V
##      [,1] [,2] [,3] [,4]
## [1,] 0.25 0.00 0.00 0.00
## [2,] 0.00 0.25 0.00 0.00
## [3,] 0.00 0.00 0.25 0.00
## [4,] 0.00 0.00 0.00 0.25
## 
## $R$nu
## [1] 0.004
## 
## 
## $G
## $G$G1
## $G$G1$V
##      [,1] [,2] [,3] [,4]
## [1,] 0.25 0.00 0.00 0.00
## [2,] 0.00 0.25 0.00 0.00
## [3,] 0.00 0.00 0.25 0.00
## [4,] 0.00 0.00 0.00 0.25
## 
## $G$G1$nu
## [1] 0.004

---

The R matrix prior is for the residuals. This is \(\mathbf{R}_{4,4}\) as we have 4 traits in the VCV matrix. We have have non-zero variances for each trait as the mode for the prior, and 0 for covariances. However this is a weak prior, so even small amounts of data will largely overcome the pull of the prior. Another sensible approach we have used is to make the prior proportional to the overall observed VCV, as there is a large literature on the proportionality of covariance matrices for morphology. This may not be sensible for other types of multivariate response measures.

##,depends.on=c("prior","mod_fm","get_data")}
t2 <- system.time(
    MMLM1.fit <- MCMCglmm(fmla.MMLM1,
                          random=~ us(trait):line, 
                          rcov=~ us(trait):units,
                          prior=  prior.model.1,
                          data= dll_data2, 
                          family = fam.test, 
                          nitt= 10000, burnin= 2000, thin=10)
)

MCMCglmm specification: fixed effects

• fmla.MMLM1 is the formula object we created above
• it looks like the lmer formula
• the trait term is a reserved word in MCMCglmm, letting it know we want to fit a multivariate mixed model
• MCMCglmm automatically melts the data for us (and assigns the name trait the same way we did manually above)

MCMCglmm: random effects

• random=~us(trait):line asks MCMCglmm to fit an unstructured covariance matrix for the line term (i.e the different wild type genotypes we are examining).- “Unstructured” means we are estimating the complete 4 x 4 matrix of covariances (= 4*5/2 = 10 elements total)
• equivalent to (trait-1|line) in lmer model
• MCMCglmm also automatically adds a units
• rcov=~us(trait):units = (trait-1|units)
• MCMCglmm offers a few other options for variance-covariance structures

MCMCglmm: other options

The nitt is how many iterations for the MCMC we want to perform, and the burnin is how many should be ignored at the beginning of the random walk.

MCMCglmm: diagnostics

The fit takes 17 seconds.

Normally we would spend a fair bit of time on diagnostics of the MCMC, but for now we will just quickly check the trace plots and autocorrelation.

In the fitted object $Sol is for solution, which is the term used for fixed effects in MCMCglmm. $VCV is for the variance-covariance matrix.

library(lattice)
xyplot(MMLM1.fit$Sol[,1:4])

xyplot(MMLM1.fit$Sol[,13:16])

acf(MMLM1.fit$Sol[,1:2])

xyplot(MMLM1.fit$VCV[,1:4])

acf(MMLM1.fit$VCV[,1:3])

Nothing terribly worrying.

summary(MMLM1.fit)

## 
##  Iterations = 2001:9991
##  Thinning interval  = 10
##  Sample size  = 800 
## 
##  DIC: 17451 
## 
##  G-structure:  ~us(trait):line
## 
##                                  post.mean l-95% CI u-95% CI eff.samp
## traitfemur_s:traitfemur_s.line      0.3044   0.1510    0.511      800
## traittibia_s:traitfemur_s.line      0.2578   0.1225    0.447      800
## traittarsus_s:traitfemur_s.line     0.2335   0.1025    0.421      800
## traitSCT_s:traitfemur_s.line        0.0335  -0.0847    0.162      800
## traitfemur_s:traittibia_s.line      0.2578   0.1225    0.447      800
## traittibia_s:traittibia_s.line      0.2830   0.1290    0.463      716
## traittarsus_s:traittibia_s.line     0.2289   0.0959    0.408      936
## traitSCT_s:traittibia_s.line        0.0495  -0.0679    0.172      800
## traitfemur_s:traittarsus_s.line     0.2335   0.1025    0.421      800
## traittibia_s:traittarsus_s.line     0.2289   0.0959    0.408      936
## traittarsus_s:traittarsus_s.line    0.2678   0.1359    0.464      800
## traitSCT_s:traittarsus_s.line       0.0824  -0.0234    0.211      800
## traitfemur_s:traitSCT_s.line        0.0335  -0.0847    0.162      800
## traittibia_s:traitSCT_s.line        0.0495  -0.0679    0.172      800
## traittarsus_s:traitSCT_s.line       0.0824  -0.0234    0.211      800
## traitSCT_s:traitSCT_s.line          0.1920   0.0876    0.329      800
## 
##  R-structure:  ~us(trait):units
## 
##                                   post.mean l-95% CI u-95% CI eff.samp
## traitfemur_s:traitfemur_s.units      0.6435   0.6050   0.6837      800
## traittibia_s:traitfemur_s.units      0.4009   0.3690   0.4394      800
## traittarsus_s:traitfemur_s.units     0.2147   0.1870   0.2481      800
## traitSCT_s:traitfemur_s.units        0.0260  -0.0095   0.0540      800
## traitfemur_s:traittibia_s.units      0.4009   0.3690   0.4394      800
## traittibia_s:traittibia_s.units      0.6646   0.6224   0.7083      800
## traittarsus_s:traittibia_s.units     0.1792   0.1467   0.2063      800
## traitSCT_s:traittibia_s.units        0.0599   0.0279   0.0921      800
## traitfemur_s:traittarsus_s.units     0.2147   0.1870   0.2481      800
## traittibia_s:traittarsus_s.units     0.1792   0.1467   0.2063      800
## traittarsus_s:traittarsus_s.units    0.5345   0.5026   0.5707      800
## traitSCT_s:traittarsus_s.units       0.0843   0.0517   0.1135      800
## traitfemur_s:traitSCT_s.units        0.0260  -0.0095   0.0540      800
## traittibia_s:traitSCT_s.units        0.0599   0.0279   0.0921      800
## traittarsus_s:traitSCT_s.units       0.0843   0.0517   0.1135      800
## traitSCT_s:traitSCT_s.units          0.7464   0.7034   0.7972      891
## 
##  Location effects: cbind(femur_s, tibia_s, tarsus_s, SCT_s) ~ trait:(genotype * temp) - 1 
## 
##                                 post.mean l-95% CI u-95% CI eff.samp  pMCMC   
## traitfemur_s:genotypeDll           0.5060   0.3062   0.7329      798 <0.001 **
## traittibia_s:genotypeDll           0.6365   0.4375   0.8516      900 <0.001 **
## traittarsus_s:genotypeDll          0.6190   0.4194   0.8195      800 <0.001 **
## traitSCT_s:genotypeDll             0.6825   0.5097   0.8745      800 <0.001 **
## traitfemur_s:genotypewt            0.1777  -0.0379   0.3749      800  0.097 . 
## traittibia_s:genotypewt            0.0206  -0.1893   0.2135      889  0.860   
## traittarsus_s:genotypewt           0.4481   0.2475   0.6408      800 <0.001 **
## traitSCT_s:genotypewt             -0.0790  -0.2510   0.0814      800  0.365   
## traitfemur_s:temp30               -0.8330  -0.9355  -0.7168      800 <0.001 **
## traittibia_s:temp30               -0.8067  -0.9141  -0.6958      800 <0.001 **
## traittarsus_s:temp30              -1.2364  -1.3373  -1.1290      800 <0.001 **
## traitSCT_s:temp30                 -0.8374  -0.9457  -0.7043      749 <0.001 **
## traitfemur_s:genotypewt:temp30     0.3579   0.2143   0.5007      800 <0.001 **
## traittibia_s:genotypewt:temp30     0.4708   0.3013   0.6004      800 <0.001 **
## traittarsus_s:genotypewt:temp30    0.5128   0.3821   0.6433      800 <0.001 **
## traitSCT_s:genotypewt:temp30       0.7252   0.5747   0.8785      800 <0.001 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Sometimes it is easier to look at the fixed and random effects separately.

summary(MMLM1.fit$Sol)

## 
## Iterations = 2001:9991
## Thinning interval = 10 
## Number of chains = 1 
## Sample size per chain = 800 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##                                    Mean     SD Naive SE Time-series SE
## traitfemur_s:genotypeDll         0.5060 0.1114  0.00394        0.00394
## traittibia_s:genotypeDll         0.6365 0.1058  0.00374        0.00353
## traittarsus_s:genotypeDll        0.6190 0.1024  0.00362        0.00362
## traitSCT_s:genotypeDll           0.6825 0.0933  0.00330        0.00330
## traitfemur_s:genotypewt          0.1777 0.1075  0.00380        0.00380
## traittibia_s:genotypewt          0.0206 0.1036  0.00366        0.00347
## traittarsus_s:genotypewt         0.4481 0.1007  0.00356        0.00356
## traitSCT_s:genotypewt           -0.0790 0.0873  0.00309        0.00309
## traitfemur_s:temp30             -0.8330 0.0557  0.00197        0.00197
## traittibia_s:temp30             -0.8067 0.0570  0.00202        0.00202
## traittarsus_s:temp30            -1.2364 0.0524  0.00185        0.00185
## traitSCT_s:temp30               -0.8374 0.0623  0.00220        0.00228
## traitfemur_s:genotypewt:temp30   0.3579 0.0735  0.00260        0.00260
## traittibia_s:genotypewt:temp30   0.4708 0.0755  0.00267        0.00267
## traittarsus_s:genotypewt:temp30  0.5128 0.0678  0.00240        0.00240
## traitSCT_s:genotypewt:temp30     0.7252 0.0789  0.00279        0.00279
## 
## 2. Quantiles for each variable:
## 
##                                    2.5%     25%     50%     75%   97.5%
## traitfemur_s:genotypeDll         0.2961  0.4328  0.5086  0.5771  0.7236
## traittibia_s:genotypeDll         0.4208  0.5658  0.6328  0.7075  0.8421
## traittarsus_s:genotypeDll        0.4206  0.5513  0.6180  0.6820  0.8237
## traitSCT_s:genotypeDll           0.5025  0.6191  0.6885  0.7445  0.8681
## traitfemur_s:genotypewt         -0.0305  0.1048  0.1769  0.2503  0.3904
## traittibia_s:genotypewt         -0.1844 -0.0463  0.0181  0.0840  0.2351
## traittarsus_s:genotypewt         0.2533  0.3835  0.4458  0.5130  0.6567
## traitSCT_s:genotypewt           -0.2475 -0.1347 -0.0776 -0.0192  0.0888
## traitfemur_s:temp30             -0.9420 -0.8712 -0.8334 -0.7963 -0.7215
## traittibia_s:temp30             -0.9142 -0.8446 -0.8097 -0.7662 -0.6963
## traittarsus_s:temp30            -1.3386 -1.2700 -1.2341 -1.2014 -1.1328
## traitSCT_s:temp30               -0.9584 -0.8779 -0.8366 -0.7925 -0.7168
## traitfemur_s:genotypewt:temp30   0.2183  0.3115  0.3549  0.4052  0.5049
## traittibia_s:genotypewt:temp30   0.3183  0.4219  0.4703  0.5198  0.6209
## traittarsus_s:genotypewt:temp30  0.3843  0.4661  0.5150  0.5580  0.6478
## traitSCT_s:genotypewt:temp30     0.5753  0.6752  0.7274  0.7733  0.8787

And for the random effects.

summary(MMLM1.fit$VCV)

## 
## Iterations = 2001:9991
## Thinning interval = 10 
## Number of chains = 1 
## Sample size per chain = 800 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##                                     Mean     SD Naive SE Time-series SE
## traitfemur_s:traitfemur_s.line    0.3044 0.1038 0.003670       0.003670
## traittibia_s:traitfemur_s.line    0.2578 0.0920 0.003254       0.003254
## traittarsus_s:traitfemur_s.line   0.2335 0.0882 0.003119       0.003119
## traitSCT_s:traitfemur_s.line      0.0335 0.0611 0.002162       0.002162
## traitfemur_s:traittibia_s.line    0.2578 0.0920 0.003254       0.003254
## traittibia_s:traittibia_s.line    0.2830 0.0946 0.003346       0.003535
## traittarsus_s:traittibia_s.line   0.2289 0.0844 0.002984       0.002758
## traitSCT_s:traittibia_s.line      0.0495 0.0598 0.002114       0.002114
## traitfemur_s:traittarsus_s.line   0.2335 0.0882 0.003119       0.003119
## traittibia_s:traittarsus_s.line   0.2289 0.0844 0.002984       0.002758
## traittarsus_s:traittarsus_s.line  0.2678 0.0913 0.003227       0.003227
## traitSCT_s:traittarsus_s.line     0.0824 0.0605 0.002140       0.002140
## traitfemur_s:traitSCT_s.line      0.0335 0.0611 0.002162       0.002162
## traittibia_s:traitSCT_s.line      0.0495 0.0598 0.002114       0.002114
## traittarsus_s:traitSCT_s.line     0.0824 0.0605 0.002140       0.002140
## traitSCT_s:traitSCT_s.line        0.1920 0.0681 0.002406       0.002406
## traitfemur_s:traitfemur_s.units   0.6435 0.0210 0.000743       0.000743
## traittibia_s:traitfemur_s.units   0.4009 0.0182 0.000643       0.000643
## traittarsus_s:traitfemur_s.units  0.2147 0.0155 0.000547       0.000547
## traitSCT_s:traitfemur_s.units     0.0260 0.0161 0.000568       0.000568
## traitfemur_s:traittibia_s.units   0.4009 0.0182 0.000643       0.000643
## traittibia_s:traittibia_s.units   0.6646 0.0228 0.000805       0.000805
## traittarsus_s:traittibia_s.units  0.1792 0.0155 0.000547       0.000547
## traitSCT_s:traittibia_s.units     0.0599 0.0166 0.000586       0.000586
## traitfemur_s:traittarsus_s.units  0.2147 0.0155 0.000547       0.000547
## traittibia_s:traittarsus_s.units  0.1792 0.0155 0.000547       0.000547
## traittarsus_s:traittarsus_s.units 0.5345 0.0180 0.000638       0.000638
## traitSCT_s:traittarsus_s.units    0.0843 0.0160 0.000565       0.000565
## traitfemur_s:traitSCT_s.units     0.0260 0.0161 0.000568       0.000568
## traittibia_s:traitSCT_s.units     0.0599 0.0166 0.000586       0.000586
## traittarsus_s:traitSCT_s.units    0.0843 0.0160 0.000565       0.000565
## traitSCT_s:traitSCT_s.units       0.7464 0.0248 0.000876       0.000830
## 
## 2. Quantiles for each variable:
## 
##                                       2.5%       25%    50%    75%  97.5%
## traitfemur_s:traitfemur_s.line     0.16059  0.229754 0.2858 0.3553 0.5578
## traittibia_s:traitfemur_s.line     0.12820  0.190895 0.2397 0.3049 0.4886
## traittarsus_s:traitfemur_s.line    0.11076  0.170365 0.2184 0.2731 0.4482
## traitSCT_s:traitfemur_s.line      -0.08481 -0.000129 0.0303 0.0670 0.1609
## traitfemur_s:traittibia_s.line     0.12820  0.190895 0.2397 0.3049 0.4886
## traittibia_s:traittibia_s.line     0.14985  0.216015 0.2648 0.3364 0.5153
## traittarsus_s:traittibia_s.line    0.10669  0.169906 0.2132 0.2709 0.4516
## traitSCT_s:traittibia_s.line      -0.05806  0.011198 0.0444 0.0819 0.1919
## traitfemur_s:traittarsus_s.line    0.11076  0.170365 0.2184 0.2731 0.4482
## traittibia_s:traittarsus_s.line    0.10669  0.169906 0.2132 0.2709 0.4516
## traittarsus_s:traittarsus_s.line   0.14568  0.204069 0.2509 0.3071 0.4935
## traitSCT_s:traittarsus_s.line     -0.02199  0.046051 0.0773 0.1123 0.2144
## traitfemur_s:traitSCT_s.line      -0.08481 -0.000129 0.0303 0.0670 0.1609
## traittibia_s:traitSCT_s.line      -0.05806  0.011198 0.0444 0.0819 0.1919
## traittarsus_s:traitSCT_s.line     -0.02199  0.046051 0.0773 0.1123 0.2144
## traitSCT_s:traitSCT_s.line         0.09875  0.143579 0.1803 0.2240 0.3511
## traitfemur_s:traitfemur_s.units    0.60549  0.629174 0.6421 0.6564 0.6853
## traittibia_s:traitfemur_s.units    0.36774  0.387918 0.4001 0.4127 0.4382
## traittarsus_s:traitfemur_s.units   0.18421  0.204318 0.2137 0.2254 0.2459
## traitSCT_s:traitfemur_s.units     -0.00691  0.014989 0.0263 0.0372 0.0592
## traitfemur_s:traittibia_s.units    0.36774  0.387918 0.4001 0.4127 0.4382
## traittibia_s:traittibia_s.units    0.62326  0.648751 0.6641 0.6799 0.7097
## traittarsus_s:traittibia_s.units   0.14940  0.169185 0.1781 0.1893 0.2108
## traitSCT_s:traittibia_s.units      0.02890  0.048569 0.0597 0.0705 0.0931
## traitfemur_s:traittarsus_s.units   0.18421  0.204318 0.2137 0.2254 0.2459
## traittibia_s:traittarsus_s.units   0.14940  0.169185 0.1781 0.1893 0.2108
## traittarsus_s:traittarsus_s.units  0.50268  0.521756 0.5340 0.5463 0.5712
## traitSCT_s:traittarsus_s.units     0.05408  0.073913 0.0840 0.0948 0.1170
## traitfemur_s:traitSCT_s.units     -0.00691  0.014989 0.0263 0.0372 0.0592
## traittibia_s:traitSCT_s.units      0.02890  0.048569 0.0597 0.0705 0.0931
## traittarsus_s:traitSCT_s.units     0.05408  0.073913 0.0840 0.0948 0.1170
## traitSCT_s:traitSCT_s.units        0.70143  0.728177 0.7464 0.7624 0.7958

This is not the most friendly output, and it takes a while to get used to. However, we can see that we still have evidence for something interesting going, and we could extract the vectors of effects as we did above.

let’s look at the VCV matrix for the random effects of line in a slightly clearer way (as a matrix)

##' extract variance-covariance matrices for MCMCglmm objects
##' does not (yet) set cor, stddev, residual-var attributes
##' may be fragile: depends on group vars not having dots in them?
VarCorr.MCMCglmm <- function(object, ...) {
    s <- summary(object$VCV)$statistics[,"Mean"]
    grps <- gsub("^[^.]+\\.([[:alnum:]]+)$","\\1",names(s))
    ss <- split(s,grps)
    getVC <- function(x) {
        nms <- gsub("^([^.]+)\\.[[:alnum:]]+$","\\1",names(x))
        n <- length(nms)
        L <- round(sqrt(n))
        dimnms <- gsub("^([^:]+):.*$","\\1",nms[1:L])
        return(matrix(x,dimnames=list(dimnms,dimnms),
                      nrow=L))
    }
    r <- setNames(lapply(ss,getVC),unique(grps))
    return(r)
}
vv <- VarCorr(MMLM1.fit)

par(mfrow=c(1,2))
corrplot.mixed(cov2cor(vv$line),upper="ellipse")
corrplot.mixed(cov2cor(vv$units),upper="ellipse")

• among-line effects are very similar to lme4 estimates (variables are in a different order)
• among-individual effects a

This is only a taste of what to do, and after this we could start asking questions about this genetic covariance matrix.

Compare fixed effects:

tt <- tidy(MMLM1.fit)
tt2 <- tidy(lmer1)
tt_comb <- bind_rows(lmer=tt,MCMCglmm=tt2,.id="model") %>%
    filter(effect=="fixed")
dwplot(tt_comb)+geom_vline(xintercept=0,lty=2)

brms

see vignette("brms_multivariate")

SLOW

library(brms)
t3 <- system.time(
    brmsfit <- brm(
        mvbind(femur_s, tibia_s, tarsus_s, SCT_s) ~ genotype*temp + (1|p|line) + (1|q|replicate),
        data = dll_data2, chains = 4, cores = 4)
)
saveRDS(brmsfit,file="../data/multiv_brmsfit.rda")

This takes a long time! (\(\approx 40\) minutes, even computing the chains in parallel)

Pros and cons

• MCMCglmm: Bayesian
  • easier to get various kinds of confidence intervals
  • flexible priors
  • more flexible variance structures
  • multiple response types within a unit (e.g. Gaussian + Poisson)
• lme4: frequentist
  • faster
  • simpler interface (??)
  • more convenient outputs
• other options: glmmTMB, brms …

FIXME:

• sort out random effects tidying
• set wildtype as baseline level
• check out convergence warnings with nloptwrap vs bobyqa (why is allFit not working?)

References

Dworkin, Ian. 2005. “Evidence for Canalization of Distal-Less Function in the Leg of Drosophila Melanogaster.” Evolution & Development 7 (2): 89–100. https://doi.org/10.1111/j.1525-142X.2005.05010.x.

---

Course home page           Email the profs           Logo information