QGparams: Quantitative Genetics parameters from GLMM estimates.
In QGglmm: Estimate Quantitative Genetics Parameters from Generalised Linear Mixed Models

Description Usage Arguments Details Value Author(s) See Also Examples

Function to estimate the quantitative genetics parameters on the observed scale based on estimates on the latent scale. The function yields a data.frame which includes the phenotypic mean and variance, as well as the additive genetic variance and heritability, on the observed scale.

1
2
3

    QGparams(mu, var.a, var.p, model = "", width = 10, predict = NULL,
             closed.form = TRUE, custom.model = NULL, n.obs = NULL, 
             cut.points = NULL, theta = NULL, verbose = TRUE)

`mu`	Latent intercept estimated from a GLMM (ignored if predict is not `NULL`). (numeric of length 1)
`var.a`	Latent additive genetic variance estimated from a GLMM. (numeric of length 1)
`var.p`	Latent total phenotypic variance estimated from a GLMM. Usually, the sum of the estimated variances of the random effects, plus the "residual" variance. (numeric of length 1)
`model`	Name of the used model, i.e. distribution.link. Ignored if `custom.model` is not NULL. (character) Available models are : "Gaussian" Gaussian distribution with identity link (e.g. LMM) "binom1.probit" Binomial with 1 trial (binary data) with a probit link "binomN.probit" Binomial with N tria with a probit link (require the parameter `n.obs`) "binom1.logit" Binomial with 1 trial (binary) with a logit link "binomN.logit" Binomial with N trial with a logit link (require the parameter `n.obs`) "Poisson.log" Poisson distribution wiht a log link "Poisson.sqrt" Poisson distribution with a square - root link "negbin.log" Negative - Binomial distribution wiht a log link (require the parameter `theta`) "negbin.sqrt" Negative - Binomial distribution with a square - root link (require the parameter `theta`) "ordinal" Multiple threshold model for ordinal categorical traits (require the parameter `cut.points`)
`width`	Parameter for the integral computation. The integral is evaluated from `mu` - `width * sqrt(var.p)` to `mu` + `width * sqrt(var.p)`. The default value is 10, which should be sensible for most models. (numeric)
`predict`	Optional vector of predicted values on the latent scale (i.e. matrix product Xb). The latent predicted values must be computed while only accounting for the fixed effects (marginal to the random effects). (numeric)
`closed.form`	When available, should closed forms be used instead of integral computations? (boolean, ignored if `model = "ordinal"`)
`custom.model`	If the model used is not available using the `model` argument, a list of functions describing the model can be provided. (list of functions, see Details)
`n.obs`	Number of "trials" for the "binomN" distribution. (numeric)
`cut.points`	Values for the "cut points" in the multiple threshold model ("ordinal"). Should be a vector of length equal to the number of categories plus one, starting with the value -Inf and ending with the value Inf. (numeric)
`theta`	Dispersion parameter for the Negative Binomial distribution. The parameter `theta` should be such as the variance of the distribution is `mean + mean^2 / theta`. (numeric)
`verbose`	Should the function be verbose? (boolean)

The function typically uses precise integral numerical approximation to compute quantitative genetics parameters on the observed scale, from latent estimates yielded by a GLMM. If closed form solutions for the integrals are available, it uses them if closed.form = TRUE.

Only the most typical distribution/link function couples are implemented in the function. If you used an "exotic" GLMM, you can use the custom.model argument. It should take the form of a list of functions. The first function should be the inverse of the link function named inv.link, the second function should be the "distribution variance" function named var.func and the third function should be the derivative of the inverse link function named d.inv.link (see Example below).

Some distributions require extra-arguments. 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).

If fixed effects (apart from the intercept) have been included in the GLMM, they can be included as marginal predicted values, i.e. predicted values excluding the random effects, which can be calculated as the matrix product Xb where X is the design matrix and b is the vector of fixed effects estimates. To do so, provide the vector of marginal predicted values using the argument predict. Note this can considerably slow down the algorithm, especially when no closed form is used.

Ordinal model is different from the other models, because it yields multivariate inference on the observed data scale, even though the latent scale is not multivariate. As a consequence, this model can only be accessed using the function QGparams and has an output similar to the one of QGmvparams.

The function yields a data.frame containing the following values:

`mean.obs`	Phenotypic mean on the observed scale.
`var.obs`	Phenotypic variance on the observed scale.
`var.a.obs`	Additive genetic variance on the observed scale.
`h2.obs`	Heritability on the observed scale.

Pierre de Villemereuil & Michael B. Morrissey

QGmvparams, QGpred, QGlink.funcs, QGmean, QGvar.dist, QGvar.exp, QGpsi

## Example using binary data
# Parameters
mu <- 0
va <- 1
vp <- 2

# Simulating data l = mu + a + e
lat<- mu + rnorm(1000, 0, sqrt(va)) + rnorm(1000, 0, sqrt(vp - va))
y<- rbinom(1000, 1, pnorm(lat))

# Expected results
QGparams(mu = 0, var.p = 2, var.a = 1, model = "binom1.probit")
# Simulated results for mean and variance
mean(y)
var(y)

# Using integral approximations
QGparams(mu = 0, var.p = 2, var.a = 1, model = "binom1.probit", closed.form = FALSE)
# Note that the approximation is exactly equal to the results obtained with the closed form

# Let's create a custom model
custom <- list(inv.link = function(x){pnorm(x)},
    var.func = function(x){pnorm(x) * (1 - pnorm(x))},
    d.inv.link = function(x){dnorm(x)})
    
QGparams(mu = 0, var.p = 2, var.a = 1, custom.model = custom)

# Using an ordinal model (with 4 categories)
QGparams(mu = 0.1, var.a = 1, var.p = 2, cut.points = c( - Inf, 0, 0.5, 1, Inf), model = "ordinal")
# Note the slightly different output (see QGmvparams)

# Integrating over a posterior distribution
# e.g. output from MCMCglmm named "model"
# df <- data.frame(mu = model$Sol[, 'intercept'], 
#                  va = model$VCV[, 'animal'], 
#                  vp = rowSums(model$VCV))
# params <- apply(df, 1, function(row){
#        QGparams(mu = row$mu, var.a = row$va, var.p = row$vp, model = "Poisson.log")
# })

Loading required package: mvtnorm
Loading required package: R2Cuba
[1] "Using the closed forms for a Binomial1-probit model."
  mean.obs var.obs  var.a.obs    h2.obs
1      0.5    0.25 0.05305165 0.2122066
[1] 0.529
[1] 0.2494084
[1] "Computing observed mean..."
[1] "Computing variances..."
[1] "Computing Psi..."
  mean.obs var.obs  var.a.obs    h2.obs
1      0.5    0.25 0.05305165 0.2122066
[1] "Computing observed mean..."
[1] "Computing variances..."
[1] "Computing Psi..."
  mean.obs var.obs  var.a.obs    h2.obs
1      0.5    0.25 0.05305165 0.2122066
[1] "Using the closed forms for an ordinal model (ignoring the closed.form argument)"
$mean.obs
[1] 0.4769798 0.1143395 0.1070148 0.3016659

$vcv.P.obs
            [,1]        [,2]        [,3]        [,4]
[1,]  0.24947007 -0.05453763 -0.05104389 -0.14388855
[2,] -0.05453763  0.10126597 -0.01223602 -0.03449232
[3,] -0.05104389 -0.01223602  0.09556262 -0.03228271
[4,] -0.14388855 -0.03449232 -0.03228271  0.21066358

$vcv.G.obs
             [,1]          [,2]          [,3]         [,4]
[1,]  0.052875103 -1.305491e-03 -0.0052947327 -0.046274880
[2,] -0.001305491  3.223269e-05  0.0001307274  0.001142531
[3,] -0.005294733  1.307274e-04  0.0005301965  0.004633809
[4,] -0.046274880  1.142531e-03  0.0046338088  0.040498540

$h2.obs
[1] 0.2119496849 0.0003182973 0.0055481580 0.1922427224