calculate vector valued Hessian in R

Question

I want to calculate a variance-covariance matrix of parameters. The parameters are obtained by a non-linear least squares fit.

library(minpack.lm)
library(numDeriv)

variables

t <- seq(0.1,20,0.3)
a <- 20
b <- 14
c <- 0.4
jitter <- rnorm(length(t),0,0.5)
Hobs <- a+b*exp(-c*t)+jitter

function def

Hhat <- function(parList, t) {parList$a + parList$b*exp(-parL
Hhatde <- function(par, t) {par[1] + par[2]*exp(-par[3]*t)}st$c*t)}
residFun <- function(par, t, observed) observed - Hhat(par,t)

initial conditions

parStart = list(a = 20, b = 10 ,c = 0.5)

nls.lm

library(minpack.lm)
out1 <- nls.lm(par = parStart, fn = residFun, observed = Hobs,
              t = t, control = nls.lm.control(nprint=0))

I wish to calculate manually what is given back via vcov(out1) I tried it with: but sigma and vcov(out1) which don't seem to be the same

J <- jacobian(Hhatde, c(19.9508523,14.6586555,0.4066367 ), method="Richardson", 
method.args=list(),t=t)
sigma <- solve((t(J)%*%J))
vcov(out1)

now trying to do it with the hessian, I can't get it working for error message cf below

hessian

H <- hessian(Hhatde, x = c(19.9508523,14.6586555,0.4066367 ), method="complex", method.args=list(),t=t)

Error in hessian.default(Hhatde, x = c(19.9508523, 14.6586555, 0.4066367),  : 
  Richardson method for hessian assumes a scalar valued function.

How do I do I get my hessian() to work.

I am not very strong on the math here, hence the trial and error approach.

Answer 1

vcov(out1) returns an estimate of the scaled variance-covariance matrix for the parameters in your model. The inverse of the cross product of the gradient, solve(crossprod(J)) returns an estimate of the unscaled variance-covariance matrix. The scaling factor is the estimated variance of the errors. So to calculate the scaled variance-covariance matrix (with some rounding error) using the gradient and the residuals from your model:

df = length(Hobs) - length(out1$par)                # degrees freedom
se_var = sum(out1$fvec^2) / df                      # estimated error variance
var_cov = se_var * solve(crossprod(J))              # scaled variance-covariance 

print(var_cov)                                      
print(vcov(out1))

To brush up on non-linear regression and non-linear least squares, you might wish to check out Seber & Wild's Nonlinear regression, or Bates & Watts' Nonlinear regression analysis and its applications. John Fox also has a short online appendix that you may find helpful.