[R] nlm, hessian, and derivatives in obj function?

[Jeff D. Hamann]

I've been working on a new package and I have a few questions regarding the
behaviour of the nlm function. I've been (for better or worse) using the nlm
function to fit a linear model without suppling the hessian or gradient
attributes in the objective function. I'm curious as to why the nlm requires
31 iterations (for the linear model), and then it doesn't work when I try to
add the derivative information. I know using nlm for a linear model isn't
the "optimal" method, but I would like to make sure the parameter estimates
and the se's are matching before I attempt more difficult problems.

rm(list=ls(all=TRUE))
print( "running nlsystemfit models test at end...")
data( kmenta )
attach( kmenta )
##demand2 <- q ~ d0 + d1 * p + d2 * d
supply2 <- q ~ s0 + s1 * p + s2 * f + s3 * a
##system2 <- list( demand2, supply2 )
##labels <- list( "Demand", "Supply" )
##inst <- ~ d + f + a
##sv2 <- c(d0=3,s2=2.123,d2=4,s0=-2.123,s3=4.234,d1=4.234,s1=0.234)
sv2 <- c(s0=-2.123,s1=0.234,s2=2.123,s3=4.234)

obj <- function( s, eqn, data, parmnames )
{
  ## get the values of the parameters
  for( i in 1:length( parmnames ) )
    {
      name <- names( parmnames )[i]
      val <- s[i]
      storage.mode( val ) <-  "double"
      assign( name, val )
    }

  lhs <- as.matrix( eval( as.formula( eqn )[[2]] ) )
  rhs <- as.matrix( eval( as.formula( eqn )[[3]] ) )
  resid <- crossprod( lhs - rhs )

  ## just how does this work...
  attr( obj, "value" ) <- resid
  attr( obj, "gradient" ) <- attr( eval( deriv3( eqn, names(
parmnames ) ) ), "gradient" )

}

res <- nlm( obj, sv2, hessian=T, eqn=supply2, data=kmenta, parmnames=sv2,
check.analyticals=T)

I haven't been able to get nlm to function as I keep getting the following
error message:

Error in nlm(obj, sv2, hessian = T, eqn = supply2, data = kmenta, parmnames
= sv2,  :
 invalid function value in 'nlm' optimizer


If I perform the fit without the derivative information, I get the correct
estimates,

$minimum
[1] 92.55106

$estimate
[1] 58.2754312  0.1603666  0.2481333  0.2483023

$gradient
[1] 8.552542e-08 9.087699e-06 5.716032e-06 2.163105e-06

$hessian
         [,1]       [,2]     [,3]     [,4]
[1,]   40.000   4000.762   3865.0   420.00
[2,] 4000.762 401486.918 386045.8 42007.76
[3,] 3865.000 386045.812 379593.1 39762.40
[4,]  420.000  42007.764  39762.4  5740.00

$code
[1] 1

$iterations
[1] 31

I was under the impression that you could also obtain the se of the
parameter estimates using the sqrt( diag( res$hessian ) ), but I haven't
been able to reproduce the se computed by the Jacobian

se <- sqrt( mse * diag( solve( crossprod( J ) ) ) )    # gives the correct
results...
hse <- sqrt( ( res$minimum / 8 ) * diag( solve( res$hessian ) ) )   # gives
similar results, but why 8?

I've tried to put the functionality to include the jacobian and hessian in
the objective function for nlm without success as I don't know what the form
of the functions will be ahead of time.

and get the se from the sqrt( diag( hessian ) ), but it's nowhere close?

Jeff.

---
Jeff D. Hamann
Hamann, Donald and Associates, Inc.
PO Box 1421
Corvallis, Oregon USA 97339-1421
(office) 541-754-1428
(cell) 541-740-5988
jeff_hamann at hamanndonald.com
www.hamanndonald.com