R help - Feb 2011 - Delta method using numerical derivatives

Dear all,

Is there a fairly general R implementation of the delta method that uses
numerical derivatives? 

I realise that the delta method has been implemented using symbolic
derivatives (e.g. alr3::delta.method, emdbook::deltamethod,
msm::deltamethod and survey:::nlcon), however possibly non-linear
estimators using the delta method with numerical derivatives can be
quite useful (e.g. predictnl in Stata and the estimate and predict
statements for proc nlmixed in SAS).

I would like something akin to the following predictnl() function:

## get some data and fit a model 
require(Epi)
data(lungDK)
fit <- glm(D/Y ~ as.factor(A5)+as.factor(A5):P5-1,
           data=lungDK, subset=(P5>=1970), weight=Y,
           family=poisson("identity")) # (Hakulinen and Dyba)
## calculate some summary statistics (e.g. cumulative incidence to 85
years)
require(plyr)
cuminc <- function(object,year=P5,newdata=NULL)
  ddply(newdata,deparse(substitute(year)),function(data)
        sum(5*predict(object,newdata=data)))[,-1,drop=TRUE]
## Now: get summary statistics and SEs using the delta method
predictnl(fit,cuminc,newdata=subset(lungDK,A5<85 & P5>=1970))


My first attempt is below, which depends on the coef() and vcov()
functions and object$coefficients. The approach seems obvious, so I
expect that this has been done before.

Kindly, Mark.


## numerically calculate the gradient (func may return a vector)
grad <- function(func,x,...) # would shadow numDeriv::grad()
  {
    h <- .Machine$double.eps^(1/3)*ifelse(abs(x)>1,abs(x),1)
    temp <- x+h
    h.hi <- temp-x
    temp <- x-h
    h.lo <- x-temp
    twoeps <- h.hi+h.lo
    nx <- length(x)
    ny <- length(func(x,...))
    if (ny==0L) stop("Length of function equals 0")
    df <- if(ny==1L) rep(NA, nx) else matrix(NA, nrow=nx,ncol=ny)
    for (i in 1L:nx) {
      hi <- lo <- x
      hi[i] <- x[i] + h.hi[i]
      lo[i] <- x[i] - h.lo[i]
      if (ny==1L)
        df[i] <- (func(hi, ...) - func(lo, ...))/twoeps[i]
      else df[i,] <- (func(hi, ...) - func(lo, ...))/twoeps[i]
      }
    return(df)
  }
## fun: takes coef as its first argument
## requires: coef() and vcov() on the object
numDeltaMethod <- function(object,fun,...) {
  coef <- coef(object)
  est <- fun(coef,...)
  Sigma <- vcov(object)
  gd <- grad(fun,coef,...)
  se.est <- as.vector(sqrt(diag(t(gd) %*% Sigma %*% gd)))
  data.frame(Estimate = est, SE = se.est)
}
predictnl <- function (object, ...) 
  UseMethod("predictnl")
predictnl.default <- function(object,fun,newdata=NULL,...)
  {
    if (is.null(newdata) && !is.null(object$data))
      newdata <- object$data
    localf <- function(coef,...)
      {
        object$coefficients = coef
        fun(object,...)
      }
    numDeltaMethod(object,localf,newdata=newdata,...)
  }