| Title: | Mixed-Effects REML Incorporating Generalized Inverses |
|---|---|
| Description: | Fit linear mixed-effects models using restricted (or residual) maximum likelihood (REML) and with generalized inverse matrices to specify covariance structures for random effects. In particular, the package is suited to fit quantitative genetic mixed models, often referred to as 'animal models'. Implements the average information algorithm as the main tool to maximize the restricted log-likelihood, but with other algorithms available. |
| Authors: | Matthew Wolak [cre, aut] |
| Maintainer: | Matthew Wolak <[email protected]> |
| License: | GPL-3 | file LICENSE |
| Version: | 1.1.0 |
| Built: | 2025-02-05 05:41:25 UTC |
| Source: | https://github.com/matthewwolak/gremlin |
Fit linear mixed-effects models using restricted (or residual) maximum likelihood (REML) and with generalized inverse matrices to specify covariance structures for random effects. In particular, the package is suited to fit quantitative genetic mixed models, often referred to as 'animal models' (Henderson 1973). Implements the average information algorithm (Johnson & Thompson 1995; Gilmour et al. 1995; Meyer & Smith 1996) as the main tool to maximize the restricted log-likelihood, but with other algorithms available.
The average information algorithm combined with sparse matrix techniques can potentially make model fitting very efficient.
Maintainer: Matthew Wolak [email protected] (ORCID)
Henderson, C.R. 1973. Sire evaluation and genetic trends. Journal of Animal Science 1973:10-41.
Johnson, D.L. and R. Thompson. 1995. Restricted maximum likelihood estimation of variance components for univariate animal models using sparse matrix techniques and average information. Journal of Dairy Science 78:449-456.
Gilmour, A.R., R. Thompson, and B.R. Cullis. 1995. Average information REML: An efficient algorithm for variance parameter estimation in linear mixed models. Biometrics 51:1440-1450.
Meyer, K. and S.P. Smith. 1996. Restricted maximum likelihood estimation for animal models using derivatives of the likelihood. Genetics, Selection, and Evolution 28:23-49.
Mrode, R.A. 2005. Linear Models for the Prediction of Animal Breeding Values, 2nd ed. CABI Publishing, Cambridge, MA, USA.
Useful links:
Report bugs at https://github.com/matthewwolak/gremlin/issues
## Not run: # Following the example from Mrode 2005, chapter 11. library(nadiv) #<-- to construct inverse of the numerator relatedness matrix pedMrode11 <- prepPed(Mrode11[, 1:3]) Ainv <- makeAinv(pedMrode11)$Ainv gr11lmm <- gremlin(WWG11 ~ sex - 1, random = ~ calf, data = Mrode11, ginverse = list(calf = Ainv), Gstart = matrix(0.2), Rstart = matrix(0.4), #<-- specify starting values maxit = 15, #<-- maximum iterations v = 2, vit = 1, #<-- moderate screen output (`v`) every iteration (`vit`) algit = "AI") #<-- only use Average Information algorithm iterations summary(gr11lmm) # Compare the model to a Linear Model with no random effects ## Use `update()` to update the model gr11lm <- update(gr11lmm, random = ~ 1) #<-- `~1`=drop all random effects summary(gr11lm) # Do analysis of variance between the two models ## See AIC or evaluate likelihood ratio against a Chi-squared distribution anova(gr11lm, gr11lmm) ## End(Not run)## Not run: # Following the example from Mrode 2005, chapter 11. library(nadiv) #<-- to construct inverse of the numerator relatedness matrix pedMrode11 <- prepPed(Mrode11[, 1:3]) Ainv <- makeAinv(pedMrode11)$Ainv gr11lmm <- gremlin(WWG11 ~ sex - 1, random = ~ calf, data = Mrode11, ginverse = list(calf = Ainv), Gstart = matrix(0.2), Rstart = matrix(0.4), #<-- specify starting values maxit = 15, #<-- maximum iterations v = 2, vit = 1, #<-- moderate screen output (`v`) every iteration (`vit`) algit = "AI") #<-- only use Average Information algorithm iterations summary(gr11lmm) # Compare the model to a Linear Model with no random effects ## Use `update()` to update the model gr11lm <- update(gr11lmm, random = ~ 1) #<-- `~1`=drop all random effects summary(gr11lm) # Do analysis of variance between the two models ## See AIC or evaluate likelihood ratio against a Chi-squared distribution anova(gr11lm, gr11lmm) ## End(Not run)
Check and/or set optimization algorithms to use. Intended for internal use within gremlin
algChk(algit, maxit, ctrl, mc)algChk(algit, maxit, ctrl, mc)
algit |
A character |
maxit |
An |
ctrl |
A |
mc |
A model |
A list containing a vector algit specifying the currently
implemented optimization algorithms for each iteration along with a vector
each containing the type of first (fdit) and second (sdit)
derivatives for each iteration (else NA if either is not applicable).
REML Likelihood Ratio Tests for gremlin models using anova()
## S3 method for class 'gremlin' anova(object, ..., model.names = NULL)## S3 method for class 'gremlin' anova(object, ..., model.names = NULL)
object |
An object of |
... |
Additional objects of |
model.names |
Optional character vector with model names to be used in the anova table |
A data.frame containing the nested comparison of model
objects via a REML likelihood ratio test.
mod11 <- gremlin(WWG11 ~ sex - 1, random = ~ sire, data = Mrode11) mod11red <- gremlinR(WWG11 ~ sex - 1, data = Mrode11) anova(mod11, mod11red)mod11 <- gremlin(WWG11 ~ sex - 1, random = ~ sire, data = Mrode11) mod11red <- gremlinR(WWG11 ~ sex - 1, data = Mrode11) anova(mod11, mod11red)
Determine whether the optimization has converged on a maximum of the log-likelihood function
ccFun(obj = NULL) ccFun1(obj = NULL) ccFun2(obj = NULL) ccFun3(obj = NULL) ccFun4(obj = NULL)ccFun(obj = NULL) ccFun1(obj = NULL) ccFun2(obj = NULL) ccFun3(obj = NULL) ccFun4(obj = NULL)
obj |
Optional gremlin model object. If |
A logical value whether the current REML iteration has passed
the convergence criteria
Meyer, K. 2019. WOMBAT A program for mixed model analyses by restricted maximum likelihood. User Notes. 27 September 2019.
grS2 <- gremlinR(WWG11 ~ sex, random = ~ sire, data = Mrode11, maxit = 2) ccFun1(grS2) ccFun2(grS2) grS <- gremlinR(WWG11 ~ sex, random = ~ sire, data = Mrode11) ccFun(grS)grS2 <- gremlinR(WWG11 ~ sex, random = ~ sire, data = Mrode11, maxit = 2) ccFun1(grS2) ccFun2(grS2) grS <- gremlinR(WWG11 ~ sex, random = ~ sire, data = Mrode11) ccFun(grS)
Only calculate values of a sparse matrix inverse corresponding to non-zero locations for the Cholesky factorization.
chol2inv_ii(L, Z = NULL)chol2inv_ii(L, Z = NULL)
L |
A lower-triangle Cholesky factorization ($L L' = C$). |
Z |
A sparse matrix containing the partial inverse of $L L'$ from a previous call to the function. Must contain the “Zdiagp” attribute. |
If $L L' = C$, function efficiently gives diag(Cinv) by only calculating elements of Cinv based on non-zero elements of $L$ and $L$. Follows the method and equations by Takahashi et al. (1973).
A sparse matrix containing the partial inverse of C ($L L'$)
along with attribute “Zdiagp” indicating the location for diagonals
of Z in slot x of the object Z.
Takahashi, Fagan, & Chin. 1973. Formation of a sparse bus impedance matrix and its application to short circuit study. 8th PICA Conference Proceedings, Minneapolis, MN.
Converts lists of (co)variance parameters either between list and
vector format or between the theta and nu scales.
stTrans(x) conTrans(Gcon, Rcon) start2theta(Gstart, Rstart, name = NULL) matlist2vech(theta) vech2matlist(vech, skeleton) theta2nu_trans(theta) nu2theta_trans(nu) theta2nu_lambda(theta, thetaG, thetaR) nu2theta_lambda(nu, sigma2e, thetaG, thetaR) nuVar2thetaVar_lambda(object) nuAI2thetaAIinv_lambda(object) nu2theta_noTrans(nu, thetaG, thetaR)stTrans(x) conTrans(Gcon, Rcon) start2theta(Gstart, Rstart, name = NULL) matlist2vech(theta) vech2matlist(vech, skeleton) theta2nu_trans(theta) nu2theta_trans(nu) theta2nu_lambda(theta, thetaG, thetaR) nu2theta_lambda(nu, sigma2e, thetaG, thetaR) nuVar2thetaVar_lambda(object) nuAI2thetaAIinv_lambda(object) nu2theta_noTrans(nu, thetaG, thetaR)
x, theta, nu
|
A |
Gcon, Rcon
|
A |
Gstart, Rstart
|
A |
name |
An (optional) character |
vech |
A |
skeleton |
An example structure to map |
thetaG, thetaR
|
A |
sigma2e |
A |
object |
An object of |
stTrans Transform start parameters into lower triangle
matrices of class dsCMatrix.
conTrans Transformation of starting constraints to correct format.
start2theta Converts lists of starting values for (co)variance parameters to a theta object used to structure the (co)variance components within gremlin.
matlist2vech Converts a list of (co)variance parameter
matrices to a vector with a “skel” attribute.
vech2matlist Converts a vector of (co)variance parameters to a list of covariance matrices.
theta2nu_trans Transforms theta to nu scale by taking the Cholesky factor of each covariance matrix and then replacing the diagonals with their (natural) logarithms. Done to ensure matrices are positive definite.
nu2theta_trans Back transformation from
theta2nu_trans: exponentiates the diagonal elements of each matrix
then calculates the cross-product.
theta2nu_lambda Transformation that factors out a residual
variance so that nu contains the ‘lambda’ parameterization:
ratios of variance parameters with the residual variance.
nu2theta_lambda Back transformation from
theta2nu_lambda.
nuVar2thetaVar_lambda Transformation of Sampling Variances
from lambda Scale for theta.
nuAI2thetaAIinv_lambda Transform AI matrix from lambda
Scale to AI-inverse of theta.
nu2theta_noTrans Structures theta when not
transformed.
Functions are specified to mostly return either a list of
matrices (structure as defined by the “skel” attribute or in
the skeleton object) or a vector containing the (co)variance
parameters of the model. Additional list elements returned can be:
A vector indexing the G-structure components.
A vector indexing the R-structure components.
Alternatively, nuVar2thetaVar_lambda and nuAI2thetaAIinv_lambda
return a vector and matrix, respectively, holding the sampling
(co)variances of the model (co)variance parameters both on the theta
scale. These are elements of the inverse Average Information matrix.
# User-specified starting parameters thetaOut <- start2theta(Gstart = list(matrix(1), matrix(2)), Rstart = matrix(3)) ## convert to a vector and then back into a matrix list thetav <- matlist2vech(thetaOut$theta) theta <- vech2matlist(thetav, attr(thetav, "skel")) identical(thetaOut$theta, theta) #<-- should be TRUE # lambda parameterization transformation nu <- theta2nu_lambda(theta, thetaOut$thetaG, thetaOut$thetaR) # back-transform from (lambda scale) nu to theta ## For example, when the sigma2e estimate=0.5 theta2 <- nu2theta_lambda(nu, sigma2e = 0.5, thetaOut$thetaG, thetaOut$thetaR)# User-specified starting parameters thetaOut <- start2theta(Gstart = list(matrix(1), matrix(2)), Rstart = matrix(3)) ## convert to a vector and then back into a matrix list thetav <- matlist2vech(thetaOut$theta) theta <- vech2matlist(thetav, attr(thetav, "skel")) identical(thetaOut$theta, theta) #<-- should be TRUE # lambda parameterization transformation nu <- theta2nu_lambda(theta, thetaOut$thetaG, thetaOut$thetaR) # back-transform from (lambda scale) nu to theta ## For example, when the sigma2e estimate=0.5 theta2 <- nu2theta_lambda(nu, sigma2e = 0.5, thetaOut$thetaG, thetaOut$thetaR)
Calculates the standard error for results of simple mathematical functions of (co)variance parameters using the delta method.
deltaSE(calc, object, scale = c("theta", "nu")) ## Default S3 method: deltaSE(calc, object, scale = c("theta", "nu")) ## S3 method for class 'formula' deltaSE(calc, object, scale = c("theta", "nu")) ## S3 method for class 'list' deltaSE(calc, object, scale = c("theta", "nu"))deltaSE(calc, object, scale = c("theta", "nu")) ## Default S3 method: deltaSE(calc, object, scale = c("theta", "nu")) ## S3 method for class 'formula' deltaSE(calc, object, scale = c("theta", "nu")) ## S3 method for class 'list' deltaSE(calc, object, scale = c("theta", "nu"))
calc |
A character |
object |
A fitted model object of |
scale |
A |
The delta method (e.g., Lynch and Walsh 1998, Appendix 1; Ver Hoef 2012) uses
a Taylor series expansion to approximate the moments of a function of
parameters. Here, a second-order Taylor series expansion is implemented to
approximate the standard error for a function of (co)variance parameters.
Partial first derivatives of the function are calculated by algorithmic
differentiation with deriv.
Though deltaSE can calculate standard errors for non-linear functions
of (co)variance parameters from a fitted gremlin model, it is limited
to non-linear functions constructed by mathematical operations such as the
arithmetic operators +, -, *, / and ^,
and single-variable functions such as exp and log. See
deriv for more information.
A data.frame containing the “Estimate” and
“Std. Error” for the mathematical function(s) of (co)variance
components.
deltaSE(default): Default method
deltaSE(formula): Formula method
deltaSE(list): List method
Lynch, M. and B. Walsh 1998. Genetics and Analysis of Quantitative Traits. Sinauer Associates, Inc., Sunderland, MA, USA.
Ver Hoef, J.M. 2012. Who invented the delta method? The American Statistician 66:124-127. DOI: 10.1080/00031305.2012.687494
# Calculate the sum of the variance components grS <- gremlin(WWG11 ~ sex - 1, random = ~ sire, data = Mrode11) deltaSE(Vsum ~ V1 + V2, grS) deltaSE("V1 + V2", grS) #<-- alternative # Calculate standard deviations (with standard errors) from variances ## Uses a `list` as the first (`calc`) argument ### All 3 below: different formats to calculate the same values deltaSE(list(SD1 ~ sqrt(V1), SDresid ~ sqrt(V2)), grS) #<-- formulas deltaSE(list(SD1 ~ sqrt(G.sire), SDresid ~ sqrt(ResVar1)), grS) deltaSE(list("sqrt(V1)", "sqrt(V2)"), grS) #<-- list of character expressions # Additive Genetic Variance calculated from observed Sire Variance ## First simulate Full-sib data set.seed(359) noff <- 5 #<-- number of offspring in each full-sib family ns <- 100 #<-- number of sires/full-sib families VA <- 1 #<-- additive genetic variance VR <- 1 #<-- residual variance datFS <- data.frame(id = paste0("o", seq(ns*noff)), sire = rep(paste0("s", seq(ns)), each = noff)) ## simulate mid-parent breeding value (i.e., average of sire and dam BV) ### mid-parent breeding value = 0.5 BV_sire + 0.5 BV_dam #### var(mid-parent BV) = 0.25 var(BV_sire) + 0.25 var(BV_dam) = 0.5 var(BV) datFS$midParBV <- rep(rnorm(ns, 0, sqrt(0.5*VA)), each = noff) ## add to this a Mendelian sampling deviation to get each offspring BV datFS$BV <- rnorm(nrow(datFS), datFS$midParBV, sqrt(0.5*VA)) datFS$r <- rnorm(nrow(datFS), 0, sqrt(VR)) #<-- residual deviation datFS$pheno <- rowSums(datFS[, c("BV", "r")]) # Analyze with a sire model grFS <- gremlin(pheno ~ 1, random = ~ sire, data = datFS) # calculate VA as 2 times the full-sib/sire variance deltaSE(VAest ~ 2*V1, grFS) # compare to expected value and simulated value VA #<-- expected var(datFS$BV) #<-- simulated (includes Monte Carlo error) # Example with `deltaSE(..., scale = "nu") ## use to demonstrate alternative way to do same calculation of inverse ## Average Information matrix of theta scale parameters when lambda = TRUE ### what is done inside gremlin::nuVar2thetaVar_lambda grSlambda <- gremlin(WWG11 ~ sex - 1, random = ~ sire, data = Mrode11, control = gremlinControl(lambda = TRUE)) dOut <- deltaSE(thetaV1 ~ V1*V2, grSlambda, "nu") #<-- V2 is sigma2e aiFnOut <- nuVar2thetaVar_lambda(grSlambda)[1] #<--variance (sqrt below) stopifnot(abs(dOut[, "Std. Error"] - sqrt(aiFnOut)) < 1e-10)# Calculate the sum of the variance components grS <- gremlin(WWG11 ~ sex - 1, random = ~ sire, data = Mrode11) deltaSE(Vsum ~ V1 + V2, grS) deltaSE("V1 + V2", grS) #<-- alternative # Calculate standard deviations (with standard errors) from variances ## Uses a `list` as the first (`calc`) argument ### All 3 below: different formats to calculate the same values deltaSE(list(SD1 ~ sqrt(V1), SDresid ~ sqrt(V2)), grS) #<-- formulas deltaSE(list(SD1 ~ sqrt(G.sire), SDresid ~ sqrt(ResVar1)), grS) deltaSE(list("sqrt(V1)", "sqrt(V2)"), grS) #<-- list of character expressions # Additive Genetic Variance calculated from observed Sire Variance ## First simulate Full-sib data set.seed(359) noff <- 5 #<-- number of offspring in each full-sib family ns <- 100 #<-- number of sires/full-sib families VA <- 1 #<-- additive genetic variance VR <- 1 #<-- residual variance datFS <- data.frame(id = paste0("o", seq(ns*noff)), sire = rep(paste0("s", seq(ns)), each = noff)) ## simulate mid-parent breeding value (i.e., average of sire and dam BV) ### mid-parent breeding value = 0.5 BV_sire + 0.5 BV_dam #### var(mid-parent BV) = 0.25 var(BV_sire) + 0.25 var(BV_dam) = 0.5 var(BV) datFS$midParBV <- rep(rnorm(ns, 0, sqrt(0.5*VA)), each = noff) ## add to this a Mendelian sampling deviation to get each offspring BV datFS$BV <- rnorm(nrow(datFS), datFS$midParBV, sqrt(0.5*VA)) datFS$r <- rnorm(nrow(datFS), 0, sqrt(VR)) #<-- residual deviation datFS$pheno <- rowSums(datFS[, c("BV", "r")]) # Analyze with a sire model grFS <- gremlin(pheno ~ 1, random = ~ sire, data = datFS) # calculate VA as 2 times the full-sib/sire variance deltaSE(VAest ~ 2*V1, grFS) # compare to expected value and simulated value VA #<-- expected var(datFS$BV) #<-- simulated (includes Monte Carlo error) # Example with `deltaSE(..., scale = "nu") ## use to demonstrate alternative way to do same calculation of inverse ## Average Information matrix of theta scale parameters when lambda = TRUE ### what is done inside gremlin::nuVar2thetaVar_lambda grSlambda <- gremlin(WWG11 ~ sex - 1, random = ~ sire, data = Mrode11, control = gremlinControl(lambda = TRUE)) dOut <- deltaSE(thetaV1 ~ V1*V2, grSlambda, "nu") #<-- V2 is sigma2e aiFnOut <- nuVar2thetaVar_lambda(grSlambda)[1] #<--variance (sqrt below) stopifnot(abs(dOut[, "Std. Error"] - sqrt(aiFnOut)) < 1e-10)
class ‘gremlin’Extracts the fixed effect estimates from a model of class ‘gremlin’.
## S3 method for class 'gremlin' fixef(object, add.dropped = FALSE, ...)## S3 method for class 'gremlin' fixef(object, add.dropped = FALSE, ...)
object |
An object of |
add.dropped |
A |
... |
Additional arguments. |
A numeric vector of fixed effect estimates.
fixef(gremlin(WWG11 ~ sex - 1, random = ~ sire, data = Mrode11))fixef(gremlin(WWG11 ~ sex - 1, random = ~ sire, data = Mrode11))
Fit and setup functions for linear mixed-effect model (Gaussian responses).
gremlin( formula, random = NULL, rcov = ~units, data = NULL, ginverse = NULL, Gstart = NULL, Rstart = NULL, Bp = NULL, Gcon = NULL, Rcon = NULL, maxit = 20, algit = NULL, vit = 1, v = 1, control = gremlinControl(), ... ) gremlinR( formula, random = NULL, rcov = ~units, data = NULL, ginverse = NULL, Gstart = NULL, Rstart = NULL, Bp = NULL, Gcon = NULL, Rcon = NULL, maxit = 20, algit = NULL, vit = 1, v = 1, control = gremlinControl(), ... ) ## S3 method for class 'gremlin' getCall(x, ...) ## S3 method for class 'gremlin' update(object, ...) gremlinSetup( formula, random = NULL, rcov = ~units, data = NULL, ginverse = NULL, Gstart = NULL, Rstart = NULL, Bp = NULL, Gcon = NULL, Rcon = NULL, maxit = 20, algit = NULL, vit = 1, v = 1, control = gremlinControl(), ... ) mkModMats( formula, random = NULL, rcov = ~units, data = NULL, subset = NULL, ginverse = NULL, na.action = na.pass, offset = NULL, contrasts = NULL, Xsparse = TRUE, ... )gremlin( formula, random = NULL, rcov = ~units, data = NULL, ginverse = NULL, Gstart = NULL, Rstart = NULL, Bp = NULL, Gcon = NULL, Rcon = NULL, maxit = 20, algit = NULL, vit = 1, v = 1, control = gremlinControl(), ... ) gremlinR( formula, random = NULL, rcov = ~units, data = NULL, ginverse = NULL, Gstart = NULL, Rstart = NULL, Bp = NULL, Gcon = NULL, Rcon = NULL, maxit = 20, algit = NULL, vit = 1, v = 1, control = gremlinControl(), ... ) ## S3 method for class 'gremlin' getCall(x, ...) ## S3 method for class 'gremlin' update(object, ...) gremlinSetup( formula, random = NULL, rcov = ~units, data = NULL, ginverse = NULL, Gstart = NULL, Rstart = NULL, Bp = NULL, Gcon = NULL, Rcon = NULL, maxit = 20, algit = NULL, vit = 1, v = 1, control = gremlinControl(), ... ) mkModMats( formula, random = NULL, rcov = ~units, data = NULL, subset = NULL, ginverse = NULL, na.action = na.pass, offset = NULL, contrasts = NULL, Xsparse = TRUE, ... )
formula |
A |
random |
A |
rcov |
A |
data |
A |
ginverse |
A |
Gstart |
A |
Rstart |
A |
Bp |
A prior specification for fixed effects. |
Gcon, Rcon
|
A |
maxit |
An |
algit |
A |
vit |
An |
v |
An |
control |
A |
... |
Additional arguments to be passed to control the model fitting. |
x, object
|
An object of class |
subset |
An expression for the subset of |
na.action |
What to do with NAs. |
offset |
Should an offset be specified. |
contrasts |
Specify the type of contrasts for the fixed effects. |
Xsparse |
Should sparse matrices be used for the fixed effects design matrix. |
A list containing an object of class grMod and, if a
model was fit (gremlin or gremlinR) then an object containing
details of the REML iterations (object itMat). An object of class
grMod contains:
The model call.
A list of the model matrices used to construct the
mixed model equations.
The response vector.
The number of responses.
The number of columns of the original response.
The fixed effects design matrix.
The number of columns in X.
The residual design matrix.
A list of the design matrices for each random term.
The number of parameters in the G structure.
A list of generalized inverse matrices.
The log-determinants of the generalized inverse matrices - necessary to calculate the log-likelihood.
A numeric value containing the sum
of the log determinants of the random effects that do not change between
log-likelihood iterations (i.e., the part of the log determinants of
(co)variance matrices to be estimated that have been factored out).
A logical indicating which terms in the random
formula have generalized inverses associated with them (non-diagonal
matrices in the Kronecker product).
Numeric vectors or scalars
describing the numbers of random effects or some function of random and
fixed effects.
(Co)variance component constraints and boundaries of the allowable parameter space for each component.
A vector of the (co)variance parameters to
be estimated by REML withe the attribute “skel” giving the
skeleton to recreate a list of matrices from this vector.
Vectors indexing the random and residual
(co)variances, respectively, in a list of (co)variance matrices (i.e.,
theta).
A list of transformed (co)variance matrices
to be fit by REML. If a residual variance has been factored out of the
mixed model equations, nu contains the ‘lambda’
parameterization with expresses the (co)variance components as ratios
of variance parameters with the residual variance. The ‘nu’ scale
(co)variances are the ones actually fit by REML.
The estimate of the factored out residual variance from
the mixed model equations (i.e., the ‘lambda’ scale)
.
An integer for the total number of (co)variances to be
estimated.
A logical indicating whether the ‘lambda’
scale parameterization has been used.
A logical to indicate if the model is univariate or not.
Sparse matrices of class Matrix that
form the mixed model equations and do not change between iterations of
REML. These are the column bound ‘X’ and ‘Z’ design
matrices for fixed and random effects, the cross-product of W,
the Right-Hand Side of the mixed model equations, and the inverse of
the fixed effect prior matrix (zeroes on the diagonal if no priors have
been specified). Note, these may be NULL if lambda=FALSE,
because the NULL objects are not used or do change between
REML iterations.
A Matrix containing the symbolic Cholesky factorization
of the coefficient matrix of the Mixed Model Equations.
A one column matrix of solutions in the mixed model
equations.
A one column matrix of variances for the solutions
to the mixed model equations. These are obtained from the diagonal of
the inverse Coefficient matrix in the mixed model equations. If lambda
is TRUE then these are on the lambda scale and must be
multiplied by sigma2e to be converted to the original data scale.
A one column matrix of residual deviations, response minus
the values expected based on the solutions, corresponding to the order
in modMats$y.
A matrix of values containing the Average Information
matrix, or second partial derivatives of the likelihood with respect to
the transformed (co)variance components (‘nu’). The inverse of
this matrix gives the sampling (co)variances of these transformed
(co)variance components.
A single column matrix of first derivatives of
the transformed (co)variance parameters (‘nu’) with respect to
the log-Likelihood.
See the parameter described above.
A character vector of REML algorithms to use in each
iteration.
A character vector of which first derivative/gradient
algorithm to use each iteration (if first derivatives of the likelihood
function are to be calculated).
A character vector of which second derivative/hessian
algorithm to use each iteration (if second derivatives of the likelihood
function are to be calculated).
See the parameter described above.
See the parameter described above.
A numeric vector of convergence criteria thresholds.
See gremlinControl and ccFun for details.
numeric values for the effective numbers to
use as “zero” and maximum negative or positive numbers. Values
less than ezero are treated as zero and fixed to this value.
Values less than -1*einf or greater than einf are
restricted to these values. See gremlinControl for more
details.
A numeric value indicating the step-reduction to use.
See gremlinControl for more details.
A numeric value indicating the increment used in the
first derivative (gradient) finite difference method.
A matrix of details about each iteration. Rows
indicate each REML iteration (rownames reflect the REML algorithm used)
and columns contain:
(Co)variance parameters.
See ‘sigma2e’ described above.
Estimates for two these two components of the log of the REML likelihoods. These are obtained from Cholesky factorization of the coefficient matrix of the mixed model equations.
The REML log-likelihood.
Time elapsed for each REML iteration.
mkModMats(): Generates model matrices.
grSire <- gremlin(WWG11 ~ sex, random = ~ sire, data = Mrode11) # Now drop sire random effects and use the `anova` method to compare models grLM <- update(grSire, random = ~ 1) #<-- use `~1` to drop all random effects anova(grSire, grLM) # Modular functions ## get model matrices for a mixed model mM11 <- mkModMats(WWG11 ~ sex - 1, random = ~ sire, data = Mrode11) ## setup model, but do not evaluate the log-likelihood grSetup <- gremlinSetup(WWG11 ~ sex - 1, random = ~ sire, data = Mrode11) ## maximize the restricted maximum likelihood grOut <- remlIt(grSetup) summary(grOut) ## Not run: # Following the example from Mrode 2005, chapter 11. library(nadiv) #<-- to construct inverse of the numerator relatedness matrix pedMrode11 <- prepPed(Mrode11[, 1:3]) Ainv <- makeAinv(pedMrode11)$Ainv gr11lmm <- gremlin(WWG11 ~ sex - 1, random = ~ calf, data = Mrode11, ginverse = list(calf = Ainv), Gstart = matrix(0.2), Rstart = matrix(0.4), #<-- specify starting values maxit = 15, #<-- maximum iterations v = 2, vit = 1, #<-- moderate screen output (`v`) every iteration (`vit`) algit = "AI") #<-- only use Average Information algorithm iterations summary(gr11lmm) # Compare the model to a Linear Model with no random effects ## Use `update()` to update the model gr11lm <- update(gr11lmm, random = ~ 1) #<-- `~1`=drop all random effects summary(gr11lm) # Do analysis of variance between the two models ## See AIC or evaluate likelihood ratio against a Chi-squared distribution anova(gr11lm, gr11lmm) ## End(Not run)grSire <- gremlin(WWG11 ~ sex, random = ~ sire, data = Mrode11) # Now drop sire random effects and use the `anova` method to compare models grLM <- update(grSire, random = ~ 1) #<-- use `~1` to drop all random effects anova(grSire, grLM) # Modular functions ## get model matrices for a mixed model mM11 <- mkModMats(WWG11 ~ sex - 1, random = ~ sire, data = Mrode11) ## setup model, but do not evaluate the log-likelihood grSetup <- gremlinSetup(WWG11 ~ sex - 1, random = ~ sire, data = Mrode11) ## maximize the restricted maximum likelihood grOut <- remlIt(grSetup) summary(grOut) ## Not run: # Following the example from Mrode 2005, chapter 11. library(nadiv) #<-- to construct inverse of the numerator relatedness matrix pedMrode11 <- prepPed(Mrode11[, 1:3]) Ainv <- makeAinv(pedMrode11)$Ainv gr11lmm <- gremlin(WWG11 ~ sex - 1, random = ~ calf, data = Mrode11, ginverse = list(calf = Ainv), Gstart = matrix(0.2), Rstart = matrix(0.4), #<-- specify starting values maxit = 15, #<-- maximum iterations v = 2, vit = 1, #<-- moderate screen output (`v`) every iteration (`vit`) algit = "AI") #<-- only use Average Information algorithm iterations summary(gr11lmm) # Compare the model to a Linear Model with no random effects ## Use `update()` to update the model gr11lm <- update(gr11lmm, random = ~ 1) #<-- `~1`=drop all random effects summary(gr11lm) # Do analysis of variance between the two models ## See AIC or evaluate likelihood ratio against a Chi-squared distribution anova(gr11lm, gr11lmm) ## End(Not run)
Change default settings for gremlin models.
gremlinControl( cctol = c(5e-04, 1e-05, 0.01, NULL), ezero = 1e-08, einf = 1e+30, step = 0.3, h = .Machine$double.eps^(1/3), lambda = FALSE, algorithm = NULL, algArgs = list() )gremlinControl( cctol = c(5e-04, 1e-05, 0.01, NULL), ezero = 1e-08, einf = 1e+30, step = 0.3, h = .Machine$double.eps^(1/3), lambda = FALSE, algorithm = NULL, algArgs = list() )
cctol |
Convergence criteria tolerances (Meyer 2007, 2019). See details
about the convergence checks in |
ezero |
Effective zero to be used, values less than this number are treated as zero and fixed to this value. |
einf |
Effective infinite value to be used, values are limited to a to this variable as a maximum. |
step |
A |
h |
A |
lambda |
A |
algorithm |
A |
algArgs |
A |
A list of class gremlinControl to be used by
gremlinSetup and later functions when fitting the model.
Meyer, K. 2007. WOMBAT - a tool for mixed model analyses in quantitative genetics by restricted maximum likelihood (REML). Journal of Zhejiang University SCIENCE B 8(11):815-821.
Meyer, K. 2019. WOMBAT A program for mixed model analyses by restricted maximum likelihood. User Notes. 27 September 2019.
str(gremlinControl())str(gremlinControl())
Extracts the log-likelihood or AIC from a gremlin model fit.
## S3 method for class 'gremlin' logLik(object, ...) npar.gremlin(object) ## S3 method for class 'gremlin' AIC(object, ..., k = 2, fxdDf = FALSE)## S3 method for class 'gremlin' logLik(object, ...) npar.gremlin(object) ## S3 method for class 'gremlin' AIC(object, ..., k = 2, fxdDf = FALSE)
object |
An object of |
... |
Additional arguments. |
k |
A numeric value for the penalty per parameter. Default is 2, as in classic AIC. |
fxdDf |
A logical indicating whether to penalize according to the number
of fixed effect parameters. Since only models fit by REML can be compared,
these must always be the same and so become a constant. Hence, the default
is |
Function npar.gremlin returns an object with attributes n.fxd
and n.bndry which give additional information about the parameters
estimated and contributing to the overall df of the model. n.fxd
returns the total number of parameters (No. fixed effects + No. (co)variance
components) minus the number of parameters constrained to a certain value. Thus,
n.fxd represents the number of parameters that can vary and, as a
consequence, affect the log-likelihood.
The attribute n.bndry reports the number of parameters that were
restrained to stay inside the boundaries of allowable parameter space (e.g.,
a variance that was not allowed to be negative).
numeric values for the log-likelihood, the number of
parameters estimated by the model (sum of fixed effects and random effect
(co)variance components), and Akaike's Information Criterion.
grS <- gremlin(WWG11 ~ sex - 1, random = ~ sire, data = Mrode11) logLik(grS) AIC(grS)grS <- gremlin(WWG11 ~ sex - 1, random = ~ sire, data = Mrode11) logLik(grS) AIC(grS)
Data from chapter 11 in Mrode 2005. The variables are as follows:
Mrode11Mrode11
An object of class data.frame with 5 rows and 5 columns.
calf. a factor with levels 4 5 6 7 8
dam. a factor with levels 2 5 6
sire. a factor with levels 1 3 4
sex. a factor with levels male female
WWG11. a numeric vector
Mrode, R.A. 2005. Linear Models for the Prediction of Animal Breeding Values, 2nd ed. CABI Publishing, Cambridge, MA, USA.
data(Mrode11)data(Mrode11)
Extract the number of 'observations' in a gremlin model fit.
## S3 method for class 'gremlin' nobs(object, use.fallback = FALSE, ...)## S3 method for class 'gremlin' nobs(object, use.fallback = FALSE, ...)
object |
An object of |
use.fallback |
logical: should fallback methods be used to try to guess the value? Included for compatibility. |
... |
Further arguments to be passed to the methods. |
A single number, usually an integer, but can be NA.
grS <- gremlin(WWG11 ~ sex - 1, random = ~ sire, data = Mrode11) nobs(grS)grS <- gremlin(WWG11 ~ sex - 1, random = ~ sire, data = Mrode11) nobs(grS)
Evaluate the REML likelihood and algorithms for iterating to find maximum REML estimates.
reml( nu, skel, thetaG, sLc, modMats, W, Bpinv, nminffx, nminfrfx, rfxlvls, rfxIncContrib2loglik, thetaR = NULL, tWW = NULL, RHS = NULL ) em( nuvin, thetaG, thetaR, conv, rfxlvls, tugug, trace, y = NULL, r = NULL, nminffx = NULL ) ai(nuvin, skel, thetaG, modMats, W, sLc, sln, r, thetaR = NULL, sigma2e = NULL) tugug_trace(thetaG, nb, rfxlvls, listGeninv, Cinv, sln, pinv = NULL) gradFun( nuvin, thetaG, rfxlvls, sln, tugug, trace, sigma2e = NULL, r = NULL, nminfrfx = NULL ) gradFun_fd(nuvin, grObj, lL, fd = c("fdiff", "cdiff", "bdiff"))reml( nu, skel, thetaG, sLc, modMats, W, Bpinv, nminffx, nminfrfx, rfxlvls, rfxIncContrib2loglik, thetaR = NULL, tWW = NULL, RHS = NULL ) em( nuvin, thetaG, thetaR, conv, rfxlvls, tugug, trace, y = NULL, r = NULL, nminffx = NULL ) ai(nuvin, skel, thetaG, modMats, W, sLc, sln, r, thetaR = NULL, sigma2e = NULL) tugug_trace(thetaG, nb, rfxlvls, listGeninv, Cinv, sln, pinv = NULL) gradFun( nuvin, thetaG, rfxlvls, sln, tugug, trace, sigma2e = NULL, r = NULL, nminfrfx = NULL ) gradFun_fd(nuvin, grObj, lL, fd = c("fdiff", "cdiff", "bdiff"))
nu, nuvin
|
A |
skel |
A skeleton for reconstructing the list of (co)variance parameters. |
thetaG, thetaR
|
|
sLc |
A sparse |
modMats |
A |
W, tWW
|
A sparse |
Bpinv |
A matrix inverse of the matrix containing the prior specification for fixed effects. |
nminffx, nminfrfx, rfxlvls
|
|
rfxIncContrib2loglik |
A |
RHS |
A sparse |
conv |
A |
tugug |
A list of numeric values for the $u_g' u_g$ products of the
solution vector for each of the |
trace |
A list of traces of the inverse coefficient matrix for each of
of the |
y |
The response vector. |
sln, r
|
Sparse |
sigma2e |
A |
nb |
The number of columns in X. |
listGeninv |
A list of generalized inverse matrices. |
Cinv |
A sparse |
pinv |
An integer vector of the matrix permutation. |
grObj |
An list of class |
lL |
A numeric value for REML log-likelihood value. |
fd |
A character indicating whether forward, combined, or backward finite differences (“fdiff”, “cdiff”, or “bdiff”, respectively) are to be calculated. |
A list or matrix containing any of the previous
parameters described above, or the following that are in addition to or
instead of parameters above:
The REML log-likelihood.
Components of the REML log-likelihood derived from the Cholesky factor of the Coefficient matrix to the Mixed Model Equations.
A vector containing the diagonal elements of the inverse
of the Coefficient matrix to the Mixed Model Equations (i.e., the
diagonal entries of Cinv).
A matrix of values containing the Average Information
matrix, or second partial derivatives of the likelihood with respect to
the transformed (co)variance components (nu). The inverse of this matrix
gives the sampling variances of these transformed (co)variance components.
A single column matrix of first derivatives of
the transformed (co)variance parameters (nu) with respect to the
log-Likelihood.
Conduct REML iterations to estimate (co)variance parameters of a linear mixed-effect model (Gaussian responses).
remlIt(grMod, ...) ## Default S3 method: remlIt(grMod, ...) ## S3 method for class 'gremlinR' remlIt(grMod, ...)remlIt(grMod, ...) ## Default S3 method: remlIt(grMod, ...) ## S3 method for class 'gremlinR' remlIt(grMod, ...)
grMod |
A gremlin model of class |
... |
Additional arguments to be passed to control the model fitting. |
A list containing an object of class grMod and
matrix containing details of the REML iterations (object
itMat). See gremlin for descriptions of grMod
and itMat objects.
remlIt(default): Default method
remlIt(gremlinR): gremlinR method
grSsetp <- gremlinSetup(WWG11 ~ sex - 1, random = ~ sire, data = Mrode11) grS <- remlIt(grSsetp)grSsetp <- gremlinSetup(WWG11 ~ sex - 1, random = ~ sire, data = Mrode11) grS <- remlIt(grSsetp)
class ‘gremlin’Residuals of class ‘gremlin’.
## S3 method for class 'gremlin' residuals(object, type = "response", scaled = FALSE, ...)## S3 method for class 'gremlin' residuals(object, type = "response", scaled = FALSE, ...)
object |
An object of |
type |
The type of residuals which should be returned. Only implement “response” currently. Can be abbreviated. |
scaled |
Logical value indicating whether to scale residuals by the residual standard deviation. |
... |
Additional arguments. |
A numeric vector of residuals.
grS <- gremlin(WWG11 ~ sex - 1, random = ~ sire, data = Mrode11) residuals(grS)grS <- gremlin(WWG11 ~ sex - 1, random = ~ sire, data = Mrode11) residuals(grS)
Extract the length of time to fit the model.
runtime(object, ...)runtime(object, ...)
object |
An object of |
... |
Further arguments to be passed to the methods. |
A numeric of class ‘difftime’ with an attribute of units (e.g., seconds or minutes).
Summarize and print results of linear mixed model fitted with gremlin.
## S3 method for class 'gremlin' summary(object, ...) ## S3 method for class 'summary.gremlin' print(x, digits = max(3, getOption("digits") - 3), ...)## S3 method for class 'gremlin' summary(object, ...) ## S3 method for class 'summary.gremlin' print(x, digits = max(3, getOption("digits") - 3), ...)
object, x
|
An object of |
... |
Additional arguments to be passed to control the output. |
digits |
An |
A list of class summary.gremlin or a printed value
to the screen with no return values.
Model log-likelihood.
Function call and model fixed, random, and residual formulae.
A numeric of class ‘difftime’ containing
the length of time to run the model. See how this is handled in
update.gremlin.
A logical indicating if the model was transformed
to the variance ratio, or lambda scale.
A named vector listing summary output for the
model residuals.
Table of variance components and approximate
standard errors (calculated from the inverse of the average information
matrix). If a (co)variance component is fixed or at the boundary of
its parameter space then an NA is returned for the standard error
and a column with constraint types is added to the table. Alternative
methods (e.g., profile likelihood CIs) should be pursued for obtaining
uncertainties associated with fixed or boundary parameters.
A matrix containing the sampling correlations
of the (co)variance components. Note this is on the underlying nu
scale that the model is fitting.
Table of fixed effects and standard errors (calculated from the corresponding diagonal elements of the inverse of the coefficient matrix, transformed where necessary).
grS <- gremlin(WWG11 ~ sex - 1, random = ~ sire, data = Mrode11) summary(grS)grS <- gremlin(WWG11 ~ sex - 1, random = ~ sire, data = Mrode11) summary(grS)
Methods to efficiently calculate a matrix trace depending on the class of matrix.
tr(X, ...) ## Default S3 method: tr(X, ...) ## S3 method for class 'dgCMatrix' tr(X, ...) ## S3 method for class 'dsCMatrix' tr(X, ...)tr(X, ...) ## Default S3 method: tr(X, ...) ## S3 method for class 'dgCMatrix' tr(X, ...) ## S3 method for class 'dsCMatrix' tr(X, ...)
X |
A matrix. |
... |
Additional arguments. |
A numeric value for the sum of the diagonal elements.
tr(default): Default method
tr(dgCMatrix): Method for matrix X of class Matrix:::dgCMatrix
tr(dsCMatrix): Method for matrix X of class Matrix:::dsCMatrix