Dealing with separate sexes in the sommer package

I'm trying to figure out how to specify a particular mixed model in the sommer package for R.

I have measured two traits on a bunch of male individuals, and two traits on a bunch of their female relatives. The aim to estimate genetic (co)variance within and between these 4 traits. However, because individuals cannot be both sexes, I want to fit the model so that it estimates the residual covariance between female trait 1 and female trait 2, and also the residual covariance between male trait 1 and male trait 2, but NOT the residual covariances between any male and female traits (since there should be no information on these covariances in the data). In MCMCglmm, one can accomplish this with code like the following (this assumes the two female traits are in columns 1 and 2, and the two male traits are in columns 3 and 4, of the matrix of response variables):

rcov = ~us(at.level(trait, 1:2)):units + us(at.level(trait, 3:4)):units

But in sommer, it seems there is no equivalent function: I get the error message

Error: On the meantime the only rcov structures available are:
     'rcov=~units' or 'rcov=~at(.):units'.

I next tried assigning each individual a number, and fitting an individual-level random effect like this:

library(sommer)

# Generate some fake data: 
# 100 males and 100 females
# Two traits are measured on each male, and two traits on each female
# 20 individuals per sex are measured for each of 5 different genotypes 
df <- data.frame(
  sex = rep(c("female", "male"), each = 100),
  female_trait_1 = c(rnorm(100), rep(NA, 100)),
  female_trait_2 = c(rnorm(100), rep(NA, 100)),
  male_trait_1 = c(rep(NA, 100), rnorm(100)),
  male_trait_2 = c(rep(NA, 100), rnorm(100)),
  genotype = rep(rep(1:5, each = 20), 2),
  individual = 1:200
)
df$genotype <- as.factor(df$genotype)
df$individual <- as.factor(df$individual)

sommer_test <- mmer2(

  # four traits as multivariate response
  cbind(female_trait_1, 
        female_trait_2,
        male_trait_1, 
        male_trait_2) ~ 1,

  # Fit the random effect of genotype (to estimate genetic covariance within and between sexes)
  # Try to fit US covariance matrices to specific levels of 'trait' (does not work)
  random =~ 
    us(trait):genotype +
    us(at.levels(trait, c("female_trait_1", "female_trait_2"))):individual +
    us(at.levels(trait, c("male_trait_1", "male_trait_2"))):individual,

  data = df
)
summary(sommer_test)

However, the latter also does not work - it runs, but the US matrix for individual is just fitted twice, to all combinations of the traits (including the male-female combinations I was trying to avoid using at.levels). So, it seems like at.levels does not work the same way as at.level in MCMCglmm, since it seems to do nothing as I have used it here.

Am I missing something key, or is this functionality missing in sommer at present?

Many thanks!

Solution

In order to deal with separate sexes in sommer you need to install sommer >=3.7. This model is straightforward manipulating the Gtc argument (in the vs function) which handles the constraints applied for the variance-covariance components. The values 1,2,3 correspond to positive, unconstrained,fixed.

Given said that, the structure of the data makes clear that the covariance among all traits can only be estimated for the "genotype" random effect:

> unsm(4)
     [,1] [,2] [,3] [,4]
[1,]    1    2    2    2
[2,]    0    1    2    2
[3,]    0    0    1    2
[4,]    0    0    0    1

whereas for "individual" is not a full unstructured among traits but instead looks like this:

> mm <- adiag1(unsm(2),unsm(2));mm
     [,1] [,2] [,3] [,4]
[1,]    1    2    0    0
[2,]    0    1    0    0
[3,]    0    0    1    2
[4,]    0    0    0    1

Long story short the model is the following:

# Generate some fake data: 
# 100 males and 100 females
# Two traits are measured on each male, and two traits on each female
# 20 individuals per sex are measured for each of 5 different genotypes 
set.seed(3434)
df <- data.frame(
  sex = rep(c("female", "male"), each = 100),
  female_trait_1 = c(rnorm(100), rep(NA, 100)),
  female_trait_2 = c(rnorm(100), rep(NA, 100)),
  male_trait_1 = c(rep(NA, 100), rnorm(100)),
  male_trait_2 = c(rep(NA, 100), rnorm(100)),
  genotype = rep(rep(1:5, each = 20), 2),
  individual = 1:200
)
df$genotype <- as.factor(df$genotype)
df$individual <- as.factor(df$individual)

library(sommer)
mm <- adiag1(unsm(2),unsm(2));mm
mix <- mmer(cbind(female_trait_1, 
                  female_trait_2,
                  male_trait_1,
                  male_trait_2) ~ 1,
     random=~vs(genotype,Gtc=unsm(4)) + vs(individual,Gtc=mm),
     rcov=~vs(units), na.method.Y = "include",
     data=df)

mix$sigma
cov2cor(mix$sigma$genotype)
cov2cor(mix$sigma$individual)

You may want to use real data to get some meaningful results, but this shows the way. Cheers.