Description Usage Arguments Details Value Author(s) See Also Examples
Function to estimate the Intra - Class Correlation coefficients (ICC, a.k.a. repeatability - like estimates) on the observed scale based on estimates on the latent scale. For a specific variance component, the function yields a data.frame which includes the phenotypic mean and variance, as well as the variance component and associated ICC, on the observed data scale.
1 2 |
mu |
Latent intercept estimated from a GLMM (ignored if predict is not |
var.comp |
Latent variance component for which ICC needs to be computed, 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
|
width |
Parameter for the integral computation. The integral is evaluated from - |
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) |
custom.model |
If the model used is not available using the |
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 precise integral numerical approximation to compute 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.
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.comp.obs |
Component variance on the observed scale. |
icc.obs |
ICC on the observed scale. |
Pierre de Villemereuil & Michael B. Morrissey
QGparams
, QGpred
, QGlink.funcs
, QGmean
, QGvar.dist
, QGvar.exp
, QGpsi
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 35 36 37 38 39 | ## Example using Poisson count data
# Parameters
mu <- 0
va <- 0.5
vm <- 0.2 # Maternal effect
vp <- 1
# Simulating data l = mu + a + e
lat <- mu +
rnorm(1000, 0, sqrt(va)) +
rnorm(1000, 0, sqrt(vm)) +
rnorm(1000, 0, sqrt(vp - (va + vm)))
y <- rpois(1000, exp(lat))
# Computing the broad - sense heritability
QGicc(mu = mu, var.p = vp, var.comp = va, model = "Poisson.log")
# Computing the maternal effect ICC
QGicc(mu = mu, var.p = vp, var.comp = vm, model = "Poisson.log")
# Using integral computation
QGicc(mu = mu, var.p = vp, var.comp = vm, model = "Poisson.log", 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){exp(x)},
var.func = function(x){exp(x)},
d.inv.link = function(x){exp(x)})
QGicc(mu = mu, var.p = vp, var.comp = vm, custom.model = custom)
# Again, exactly equal
# Integrating over a posterior distribution
# e.g. output from MCMCglmm named "model"
# df <- data.frame(mu = model$Sol[, 'intercept'],
# vm = model$VCV[, 'mother'],
# vp = rowSums(model$VCV))
# params <- apply(df, 1, function(row){
# QGicc(mu = row$mu, var.comp = row$vm, var.p = row$vp, model = "Poisson.log")
# })
|
[1] "Using the closed forms for a Poisson-log model."
mean.obs var.obs var.comp.obs icc.obs
1 1.648721 6.319496 1.763407 0.2790424
[1] "Using the closed forms for a Poisson-log model."
mean.obs var.obs var.comp.obs icc.obs
1 1.648721 6.319496 0.6018351 0.09523467
[1] "Computing observed mean..."
[1] "Computing phenotypic variance..."
[1] "Computing component variance..."
mean.obs var.obs var.comp.obs icc.obs
1 1.648721 6.319496 0.6018351 0.09523467
[1] "Computing observed mean..."
[1] "Computing phenotypic variance..."
[1] "Computing component variance..."
mean.obs var.obs var.comp.obs icc.obs
1 1.648721 6.319496 0.6018351 0.09523467
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.