predictnl: Estimation of standard errors using the numerical delta...

In mclements/numdelta: Numerical delta method

Description

A simple, yet exceedingly useful, approach to estimate the variance of a function using the numerical delta method. A number of packages provide functions that analytically calculate the gradients; we use numerical derivatives, which generalises to models that do not offer analytical derivatives (e.g. ordinary differential equations, integration), or to examples that are tedious or error-prone to calculate (e.g. sums of predictions from GLMs).

Usage

## Default S3 method:
predictnl(object, fun, newdata=NULL, ...)
## S3 method for class 'lm'
predictnl(object, fun, newdata=NULL, ...)
## S4 method for signature 'mle'
predictnl(object, fun, ...)
## S4 method for signature 'mle2'
predictnl(object, fun, ...)

Arguments

object	An object with coef, vcov and `coef<-` methods (required).
fun	A function that takes object as the first argument, possibly with newdata and other arguments (required). See notes for why it is often useful to include newdata as an argument to the function.
newdata	An optional argument that defines newdata to be passed to fun.
...	Other arguments that are passed to fun.

Details

The signature for fun is either fun(object, ...) or fun(object, newdata=NULL, ...).

The different predictnl methods call the utility function numDeltaMethod, which in turn calls the grad function for numerical differentiation. The numDeltaMethod function calls the standard coef and vcov methods, and the non-standard `coef<-` method for changing the coefficients in a regression object. This non-standard method has been provided for several regression objects and essentially mirrors the coef method.

One potential issue is that some predict methods do not re-calculate their predictions for the fitted dataset (i.e. when newdata=NULL). As the predictnl function changes the fitted coefficients, it is required that the predictions are re-calculated. One solution is to pass newdata as an argument to both predictnl and fun; alternatively, newdata can be specified in fun or a specialised numdelta:::predict method can be called. These approaches are described in the examples below. The numDeltaMethod method called by predictnl provides a warning when the variance estimates are zero, which may be due to this cause.

For completeness, it is worth discussing why the example predictnl(fit,predict) does not work for when fit is a glm object. First, predict.glm does not update the predictions for the fitted data. Second, the default predict method has a signature predict(object, ...), which does not include a newdata argument. We could then either (i) require that a newdata argument be passed to the fun function for all examples, which would make this corner case work, or (ii) only pass the newdata argument if it is non-null or in the formals for the fun function, which would fail for this corner case. The current API defaults to the latter case (ii). To support this approach, the predictnl.lm method replaces a null newdata with object$data. We also provide a revised numdelta:::predict.lm method that performs the same operation, although its use is not encouraged due to its clumsiness.

Value

Returns an object of class predictnl with methods print and confint.

Examples

df <- data.frame(x=0:1, y=c(10, 20))
fit <- glm(y ~ x, df, family=poisson)

## Some of the prediction methods do not recalculate the predictions
## if the predictions are within the fitted data (that is, newdata=NULL). 
predictnl(fit,
          function(obj)
          diff(predict(obj,type="response"))) ## WRONG - with warning
## We can fix this in three ways.
## First, we can pass a newdata argument
predictnl(fit, function(obj,newdata)
          diff(predict(obj,newdata,type="response")))
## Second, we can specify the newdata argument in the prediction call
predictnl(fit,
          function(obj)
          diff(predict(obj,newdata=obj$data,type="response")))
## Third, we can use a specialised (unexported) predict method from numdelta
predictnl(fit,
          function(obj)
          diff(numdelta:::predict.lm(obj,type="response")))

## Simpler example
predict(fit,se.fit=TRUE)                # what we should get
predictnl(fit,predict)                  # WRONG - with warning 
predictnl(fit,predict.glm)              # okay 
predictnl(fit,predict,newdata=fit$data) # okay
predictnl(fit,numdelta:::predict.lm)    # okay

## A more complex example
## Not run: 
## log-hazard survival with natural splines
require(bbmle)
require(splines)
require(survival)
set.seed(12345)
y <- rweibull(100,2,10)
loghX <- ns(y,df=5,intercept=TRUE)
h <- function(coef,x) exp(drop(predict(loghX,x) %*% coef))
S <- function(coef,time)
     sapply(time, function(yi)
            exp(-integrate(function(x) h(coef,x),0,yi)$value))
negloglike <- function(coef) -sum(log(S(coef,y)) + log(h(coef,y)))
init <- structure(c(-1.4696, -2.1932, 1.2183, -5.1613, 0.607806),
            .Names = c("beta1", "beta2", "beta3", "beta4", "beta5"))
parnames(negloglike) <- names(init)
fit <- mle2(negloglike,init,vecpar=TRUE)

x <- seq(0,max(y),length=301)
pred1 <- predictnl(fit,function(obj) log(-log(S(coef(obj),x))))
plot(x,S(coef(fit),x),type="l")
lines(survfit(Surv(y,rep(TRUE,length(y)))~1),lty=2)
matlines(x,exp(-exp(confint(pred1))),col="blue",lty=1)
lines(x,pweibull(x,2,10,lower.tail=FALSE),col="red")

pred2 <- predictnl(fit,function(obj) log(h(coef(obj),x)))
plot(x,h(coef(fit),x),type="l")
matlines(x,exp(confint(pred2)),col="blue",lty=1)
lines(x,dweibull(x,2,10)/pweibull(x,2,10,lower.tail=FALSE),col="red")

## End(Not run)
 

mclements/numdelta documentation built on May 22, 2019, 3:10 p.m.