Function to estimate the variancecovariance matrix of a variance component on the observed scale based on estimates on the latent scale. Contrary to the univariate function, this one cannot use the analytical closed forms and yields a list of paramaters instead of a data.frame.
1 2 
mu 
Vector of latent intercepts estimated from a GLMM (ignored if predict is not 
vcv.comp 
Component variancecovariance matrix (Gmatrixlike). (numeric) 
vcv.P 
Total phenotypic variancecovariance matrix. Usually, the sum of all the estimated variancecovariance matrices. (numeric) 
models 
A vector containing the names of the model used or a list which elements contain the list of the functions needed (inverselink, distribution variance and derivative of the inverselink, as stated in the output of

rel.acc 
Relative accuracy of the integral approximation. (numeric) 
width 
Parameter for the integral computation. The default value is 10, which should be sensible for most models. (numeric) 
predict 
Optional matrix of predicted values on the latent scale (each trait in each column). The latent predicted values must be computed while only accounting for the fixed effects (marginal to the random effects). (numeric) 
n.obs 
Number of "trials" for the "binomN" distribution. (numeric) 
theta 
Dispersion parameter for the Negative Binomial distribution. The parameter 
verbose 
Should the function be verbose? (boolean) 
The function typically uses integral numerical approximation provided by the R2Cuba package to compute multivariate quantitative genetics parameters on the observed scale, from latent estimates yielded by a GLMM. It cannot use closed form solutions.
Only the most typical distribution/link function couples are implemented through the models
argument. If you used an "exotic" GLMM, you can provide a list containg lists of functions corresponding to the model. The list of functions should be implemented as is the output of QGlink.funcs()
, i.e. three elements: the inverse link functions named inv.link
, the derivative of this function named d.inv.link
and the distribution variance named var.func
(see Example below).
Some distributions require extraarguments. This is the case for "binomN", which require the number of trials N, passed with the argument n.obs
. The distribution "negbin" requires a dispersion parameter theta
, such as the variance of the distribution is mean + mean**2/theta
(mean/dispersion parametrisation). For now, the arguments n.obs
and theta
can be used for ONE distribution only.
If fixed effects (apart from the intercept) have been included in the GLMM, they can be included through the argument predict
as a matrix of the marginal predicted values, i.e. predicted values excluding the random effects, for each trait (one trait per column of the matrix, see Example below).Note that computation can be extremely slow in that case.
The function yields a list containing the following values:
mean.obs 
Vector of phenotypic means on the observed scale. 
vcv.P.obs 
Phenotypic variancecovariance matrix on the observed scale. 
vcv.comp.obs 
Component variancecovariance (Gmatrixlike, but broadsense) on the observed scale. 
Pierre de Villemereuil & Michael B. Morrissey
QGmvparams
, QGlink.funcs
, QGmvmean
, QGvcov
, QGmvpsi
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34  ##Example using a bivariate model (Binary trait/Gaussian trait)
#Parameters
mu=c(0,1)
G=diag(c(0.5,2))
M=diag(c(0.2,1)) #Maternal effect VCV matrix
P=diag(c(1,4))
#Broadsense "Gmatrix" on observed data scale
## Not run: QGmvicc(mu=mu,vcv.comp=G,vcv.P=P,models=c("binom1.probit","Gaussian"))
#Maternal effect VCV matrix on observed data scale
## Not run: QGmvicc(mu=mu,vcv.comp=M,vcv.P=P,models=c("binom1.probit","Gaussian"))
#Reminder: the results are the same here because we have no correlation between the two traits
#Defining the model "by hand" using the list
list.models=list(
model1=list(inv.link=function(x){pnorm(x)},
d.inv.link=function(x){dnorm(x)},
var.func=function(x){pnorm(x)*(1pnorm(x))}),
model2=list(inv.link=function(x){x},
d.inv.link=function(x){1},
var.func=function(x){0})
)
#Running the same analysis than above
## Not run: QGmvicc(mu=mu,vcv.comp=M,vcv.P=P,models=list.models)
#Using predicted values
#Say we have 100 individuals
n=100
#Let's simulate predicted values
p<matrix(c(runif(n),runif(n)),ncol=2)
#Note that p has as many as columns as we have traits (i.e. two)
#Multivariate analysis with predicted values
## Not run: QGmvicc(predict=p,vcv.comp=M,vcv.P=P,models=c("binom1.probit","Gaussian"))
#That can be extremely long to run!!

Questions? Problems? Suggestions? Tweet to @rdrrHQ or email at ian@mutexlabs.com.
All documentation is copyright its authors; we didn't write any of that.