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

[Arne Henningsen]

Hi,

I don't know much about non-linear models, but there is another possibility to 
fit these models: 

1) get some starting values for the parameters
2) take the derivatives of the model with respect to the parameters at the 
point of the starting values of the parameters
3) perform a linear estimation of this linearized model (using systemfit) to 
get new parameter estimates
4) got to step 2) and take these new parameter estimates in place of the 
starting values
5) iterate this until the parameters stay stable from one to the next 
iteration

This has three advantages:
1) It is not much work to write these function since systemfit already exists
2) If the model is linear in parameters, it is identical to the linearized 
model and, thus, the first iteration leads directly to the optimum
3) You get get the SEs from the last iteration of systemfit

Does this approach also have disadvantages (e.g. non-convergence of parameters 
in many cases)?

Best wishes,
Arne


On Saturday 18 October 2003 01:54, Jeff D. Hamann wrote:
> 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
>
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://www.stat.math.ethz.ch/mailman/listinfo/r-help

-- 
Arne Henningsen
Department of Agricultural Economics
Christian-Albrechts-University Kiel
24098 Kiel, Germany
Tel: +49-431-880-4445
Fax: +49-431-880-1397 
ahenningsen at email.uni-kiel.de
http://www.uni-kiel.de/agrarpol/ahenningsen/