deltaSE {gremlin}

Delta Method to Calculate Standard Errors for Functions of (Co)variances.

Description

Calculates the standard error for results of simple mathematical functions of (co)variance parameters using the delta method.

Usage

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"))

Arguments

calc

A character expression, formula, or list (of formula or expression) expressing the mathematical function of (co)variance component for which to calculate standard errors.

object

A fitted model object of class ‘gremlin’.

scale

A character indicating whether to calculate the function and standard error on the original data scale (“theta”) or on the underlying scale to which (co)variance components are transformed for the model fitting calculations (“nu”). Defaults to “theta” if not specified.

Details

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.

Value

A data.frame containing the “Estimate” and “Std. Error” for the mathematical function(s) of (co)variance components.

Methods (by class)

• default: Default method

• formula: Formula method

• list: List method

Author(s)

[e-mail]

References

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

See Also

deriv

Examples

  # 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 
    dOut <- deltaSE(thetaV1 ~ V1*V2, grS, "nu")  #<-- V2 is sigma2e
    aiFnOut <- nuVar2thetaVar_lambda(grS)[1]  #<-- variance (do sqrt below)
    stopifnot(abs(dOut[, "Std. Error"] - sqrt(aiFnOut)) < 1e-10)