# deltaSE: Delta Method to Calculate Standard Errors for Functions of... In gremlin: Mixed-Effects REML Incorporating Generalized Inverses

## Description

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

## Usage

 ``` 1 2 3 4 5 6 7 8 9 10``` ```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

## 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

`deriv`
 ``` 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 40 41 42 43 44``` ``` # 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) #<-- variance (do sqrt below) stopifnot(abs(dOut[, "Std. Error"] - sqrt(aiFnOut)) < 1e-10) ```