Fix variance components in sommer's mmer2?

Viewed 115

When fitting a model with mmer2() of the sommer package - is it possible to fix variance components to a certain value similar to what we can do e.g. with the PARMS statement in SAS or the G.param and R.param arguments in ASReml-R v3?

I am specifically interested in fixing individual components irrespective of whether they are on the G-side or R-side of my mixed model.

1 Answers

If you have sommer >= 3.7 forcing specific variance or covariance components can be done using the Gt (initial values) and Gtc (contraints) arguments from the vs() function that is used to specify the variance model for a random effect.

For example assume you fit the following multivariate mixed model for two traits:

library(sommer)
data(DT_cpdata)
#### create the variance-covariance matrix for id levels
A <- A.mat(GT) # additive relationship matrix
D <- D.mat(GT) # additive relationship matrix
#### look at the data and fit the model
DT$idd <- DT$id
head(DT)
ans.m <- mmer(cbind(Yield,color)~1,
               random=~ vs(id, Gu=A)
               + vs(Rowf,Gtc=diag(2)),
               rcov=~ vs(units),
               data=DT)

In this model there's 2 random effects (id and Rowf) and a residual term for units. By default an unstructured variance model is fitted (i.e. units and id) but you can see how for Rowf the Gtc argument is a constraint matrix that makes a diagonal model for Rowf:

> diag(2)
     [,1] [,2]
[1,]    1    0
[2,]    0    1 

which follows the rules: 0: not to be estimated 1: estimated and constrained to be positive 2: estimated and unconstrained 3: fixed variance-covariance component provided in Gt

After that is easy to see that if you want to force certain variance components you have to provide initial values (Gt argument) and specify a constraint matrix with values 3's (the function fixm() can creates such matrix).

ans.mf <- mmer(cbind(Yield,color)~1,
              random=~ vs(id, Gu=A, Gt=ans.m$sigma_scaled$id,Gtc=fixm(2))
              + vs(Rowf,Gt=ans.m$sigma_scaled$Rowf, Gtc=fixm(2))
              + vs(idd, Gu=D),
              rcov=~ vs(units, Gt=ans.m$sigma_scaled$units, Gtc=fixm(2)),
              data=DT)

ans.m$sigma
ans.mf$sigma

If you check the results you can see that you have forced the variance components for id, Rowf and units and estimated properly the components for idd.

Related