Estimate and Standard Error of a Nonlinear Function of Estimated Regression Coefficients

deltaMethod is a generic function that uses the delta method to get a first-order approximate standard error for a nonlinear function of a vector of random variables with known or estimated covariance matrix.

Usage

deltaMethod(object, ...)

# Default S3 method
deltaMethod(object, g., vcov., func=g., constants, level=0.95, 
  rhs, singular.ok=FALSE, ..., envir=parent.frame())
# S3 method for class 'lm'
deltaMethod (object, g., vcov.=vcov(object, complete=FALSE), 
           parameterNames=names(coef(object)), ..., envir=parent.frame())
# S3 method for class 'nls'
deltaMethod(object, g., vcov.=vcov(object, complete=FALSE), ..., envir=parent.frame())
# S3 method for class 'multinom'
deltaMethod (object, g., vcov. = vcov(object, complete=FALSE), 
           parameterNames = if (is.matrix(coef(object))) 
           colnames(coef(object)) else names(coef(object)), ..., envir=parent.frame()) 
# S3 method for class 'polr'
deltaMethod (object, g., vcov.=vcov(object, complete=FALSE), ..., envir=parent.frame())
# S3 method for class 'survreg'
deltaMethod (object, g., vcov. = vcov(object, complete=FALSE), 
           parameterNames = names(coef(object)), ..., envir=parent.frame())
# S3 method for class 'coxph'
deltaMethod (object, g., vcov. = vcov(object, complete=FALSE), 
           parameterNames = names(coef(object)), ..., envir=parent.frame())
# S3 method for class 'mer'
deltaMethod (object, g., vcov. = vcov(object, complete=FALSE),
           parameterNames = names(fixef(object)), ..., envir=parent.frame())
# S3 method for class 'merMod'
deltaMethod (object, g., vcov. = vcov(object, complete=FALSE),
           parameterNames = names(fixef(object)), ..., envir=parent.frame())
# S3 method for class 'lme'
deltaMethod (object, g., vcov. = vcov(object, complete=FALSE),
           parameterNames = names(fixef(object)), ..., envir=parent.frame())
# S3 method for class 'lmList'
deltaMethod (object, g.,  ..., envir=parent.frame())

Arguments

object

For the default method, object is either (1) a vector of p named elements, so names(object) returns a list of p character strings that are the names of the elements of object; or (2) a model object for which there are coef and vcov methods, and for which the named coefficient vector returned by coef is asymptotically normally distributed with asymptotic covariance matrix returned by vcov. For the other methods, object is a regression object for which coef(object) or fixef(object) returns a vector of parameter estimates.

g.

A quoted string that is the function of the parameter estimates to be evaluated; see the details below.

vcov.

The (estimated) covariance matrix of the coefficient estimates. For the default method, this argument is required. For all other methods, this argument must either provide the estimated covariance matrix or a function that when applied to object returns a covariance matrix. The default is to use the function vcov.

Note that arguments supplied to ... are not passed to vcov. when it's a function; in this case either use an anonymous function in which the additional arguments are set, or supply the coefficient covariance matrix directly.

func

A quoted string used to annotate output. The default of func = g. is usually appropriate.

parameterNames

A character vector of length p that gives the names of the parameters in the same order as they appear in the vector of estimates. This argument will be useful if some of the names in the vector of estimates include special characters, like I(x2^2), or x1:x2 that will confuse the numerical differentiation function. See details below.

constants

This argument is a named vector whose elements are constants that are used in the f argument. It isn't generally necessary to specify this argument but it may be convenient to do so.

level

level for confidence interval, default 0.95.

rhs

hypothesized value for the specified function of parameters; if absent no hypothesis test is performed.

singular.ok

if TRUE (the default is FALSE), aliased coefficients are eliminated from the coefficient vector and covariance matrix; otherwise, aliased coefficients produce an error.

...

Used to pass arguments to the generic method.

envir

Environment in which g. is evaluated; not normally specified by the user.

Details

Suppose \(x\) is a random vector of length \(p\) that is at least approximately normally distributed with mean \(\beta\) and estimated covariance matrix \(C\). Then any function \(g(\beta)\) of \(\beta\), is estimated by \(g(x)\), which is in large samples normally distributed with mean \(g(\beta)\) and estimated variance \(h'Ch\), where \(h\) is the first derivative of \(g(\beta)\) with respect to \(\beta\) evaluated at \(x\). This function returns both \(g(x)\) and its standard error, the square root of the estimated variance.

The default method requires that you provide \(x\) in the argument object, \(C\) in the argument vcov., and a text expression in argument g. that when evaluated gives the function \(g\). The call names(object) must return the names of the elements of x that are used in the expression g..

Since the delta method is often applied to functions of regression parameter estimates, the argument object may be the name of a regression object from which the estimates and their estimated variance matrix can be extracted. In most regression models, estimates are returned by the coef(object) and the variance matrix from vcov(object). You can provide an alternative function for computing the sample variance matrix, for example to use a sandwich estimator.

For mixed models using lme4 or nlme, the coefficient estimates are returned by the fixef function, while for multinom, lmList and nlsList coefficient estimates are returned by coef as a matrix. Methods for these models are provided to get the correct estimates and variance matrix.

The argument g. must be a quoted character string that gives the function of interest. For example, if you set m2 <- lm(Y ~ X1 + X2 + X1:X2), then deltaMethod(m2,"X1/X2") applies the delta method to the ratio of the coefficient estimates for X1 and X2. The argument g. can consist of constants and names associated with the elements of the vector of coefficient estimates.

In some cases the names may include characters such as the colon : used in interactions, or mathematical symbols like + or - signs that would confuse the function that computes numerical derivatives, and for this case you can replace the names of the estimates with the parameterNames argument. For example, the ratio of the X2 main effect to the interaction term could be computed using deltaMethod(m2, "b1/b3", parameterNames=c("b0", "b1", "b2", "b3")). The name “(Intercept)” used for the intercept in linear and generalized linear models is an exception, and it will be correctly interpreted by deltaMethod. Another option is to use back-ticks to quote nonstandard R names, as in deltaMethod(m2,"X1/`X1:X2`").

For multinom objects, the coef function returns a matrix of coefficients, with each row giving the estimates for comparisons of one category to the baseline. The deltaMethod function applies the delta method to each row of this matrix. Similarly, for lmList and nlsList objects, the delta method is computed for each element of the list of models fit.

For nonlinear regression objects produced by the nls function, the call coef(object) returns the estimated coefficient vectors with names corresponding to parameter names. For example, m2 <- nls(y ~ theta/(1 + gamma * x), start = list(theta=2, gamma=3)) will have parameters named c("theta", "gamma"). In many other familiar regression models, such as those produced by lm and glm, the names of the coefficient estimates are the corresponding regressor names, not parameter names.

For mixed-effects models fit with lmer and glmer from the lme4 package or lme and nlme from the nlme package, only fixed-effect coefficients are considered.

For regression models for which methods are not provided, you can extract the named vector of coefficient estimates and and estimate of its covariance matrix and then apply the default deltaMethod function.

Note: Earlier versions of deltaMethod included an argument parameterPrefix that implemented the same functionality as the parameterNames argument, but which caused several problems that were not easily fixed without the change in syntax.

Value

An object of class "deltaMethod", inheriting from "data.frame", for which a print method is provided. The object contains columns named Estimate for the estimate, SE for its standard error, and columns for confidence limits and possibly a hypothesis test. The value of g. is given as a row label.

References

Fox, J. (2016) Applied Regression Analysis and Generalized Linear Models, Third Edition. Sage.

Fox, J. and Weisberg, S. (2019) An R Companion to Applied Regression, Third Edition, Sage.

Weisberg, S. (2014) Applied Linear Regression, Fourth Edition, Wiley, Section 6.1.2.

Author

Sanford Weisberg, sandy@umn.edu, John Fox jfox@mcmaster.ca, and Pavel Krivitsky.

Examples

m1 <- lm(time ~ t1 + t2, data = Transact) 
deltaMethod(m1, "b1/b2", parameterNames= paste("b", 0:2, sep="")) 
#>       Estimate      SE   2.5 % 97.5 %
#> b1/b2  2.68465 0.31899 2.05945 3.3099
deltaMethod(m1, "t1/t2", rhs=1) # use names of preds. rather than coefs.
#>       Estimate      SE   2.5 %  97.5 % Hypothesis z value  Pr(>|z|)    
#> t1/t2  2.68465 0.31899 2.05945 3.30985    1.00000  5.2813 1.283e-07 ***
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
deltaMethod(m1, "t1/t2", vcov=hccm) # use hccm function to est. vars.
#>       Estimate      SE   2.5 % 97.5 %
#> t1/t2  2.68465 0.55836 1.59030  3.779
deltaMethod(m1, "1/(Intercept)")
#>               Estimate         SE      2.5 % 97.5 %
#> 1/Intercept  0.0069267  0.0081825 -0.0091107  0.023
# The next example invokes the default method by extracting the
# vector of estimates and covariance matrix explicitly
deltaMethod(coef(m1), "t1/t2", vcov.=vcov(m1))
#>       Estimate      SE   2.5 % 97.5 %
#> t1/t2  2.68465 0.31899 2.05945 3.3099
# Mroz example, in carData, logistic regression using 2 predictors
mroz.mod1 <- glm(lfp ~ k5 + age, family=binomial, data=Mroz)
# The S() method in car returns exponentiated logistic regression coefficients by default
S(mroz.mod1)
#> Call: glm(formula = lfp ~ k5 + age, family = binomial, data = Mroz)
#> 
#> Coefficients:
#>             Estimate Std. Error z value Pr(>|z|)    
#> (Intercept)  3.08578    0.49712   6.207 5.39e-10 ***
#> k5          -1.32042    0.18638  -7.085 1.39e-12 ***
#> age         -0.05847    0.01090  -5.364 8.13e-08 ***
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> 
#> (Dispersion parameter for binomial family taken to be 1)
#> 
#>     Null deviance: 1029.75  on 752  degrees of freedom
#> Residual deviance:  964.48  on 750  degrees of freedom
#> 
#>  logLik      df     AIC     BIC 
#> -482.24       3  970.48  984.36 
#> 
#> Number of Fisher Scoring iterations: 4
#> 
#> Exponentiated Coefficients and Confidence Bounds
#>               Estimate     2.5 %     97.5 %
#> (Intercept) 21.8845911 8.3698492 58.8728998
#> k5           0.2670230 0.1832502  0.3808541
#> age          0.9432093 0.9230190  0.9633534
#> 
# Exponentiated coefficient estimates can also be computed one at a time by deltaMethod
deltaMethod(mroz.mod1, "exp(k5)")
#>         Estimate       SE    2.5 % 97.5 %
#> exp(k5) 0.267023 0.049767 0.169482 0.3646
# The point estimate is the same as in the S() output but the confidence interval is 
# different because S() simply exponentiated the ends of the confidence interval before 
# exponentiation.
# A more complex case is to get the estimate and standard error for the inverse logit
# of the coefficient estimate.
deltaMethod(mroz.mod1, "exp(b1)/(1 + exp(b1))", parameterNames= paste("b", 0:2, sep=""))
#>                       Estimate       SE    2.5 % 97.5 %
#> exp(b1)/(1 + exp(b1)) 0.210748 0.031001 0.149988 0.2715
# the parameterNames argument was used to change the name of the coefficient from
# k5 to b1.