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,...) }