[R] names of parameters from nonlinear model?

Jeff D. Hamann jeff.hamann at forestinformatics.com
Tue Dec 2 09:40:25 CET 2003


I've been trying to figure out how to build a list of terms from a nonlinear
model (terms() returns a error). I need to compute and evaluate the partial
derivatives (Jacobian) for each equaiton in  a set of equations.

For example:

> eqn <- q ~ s0 + s1 * p + s2 * f + s3 * a
> sv2 <- c(d0=3,d1=4.234,d2=4,s0=-2.123,s1=0.234,s2=2.123,s3=4.234)
> names( sv2 )
[1] "d0" "d2" "d1" "s0" "s2" "s3" "s1"
...after some processing....
> sv2
        d0         d2         d1         s0         s2         s3         s1
99.8954229  0.3346356 -0.3162988 58.2754311  0.2481333  0.2483023  0.1603666

# error becuase of the gradient
> nlols.derivs <- eval( deriv( eqn, names( sv2 ), hessian=T ) )
> J <- attr( nlols.derivs, "gradient" )

in this case J is

> J
      d0 d2 d1 s0    s2 s3      s1
 [1,]  0  0  0  1  98.0  1 100.323
 [2,]  0  0  0  1  99.1  2 104.264
...blah, blah, blah...
[6,]  0  0  0  1 108.2  6  99.456
 19,]  0  0  0  1  89.3 19 105.769
[20,]  0  0  0  1  93.0 20 113.490
> se <- sqrt( nlols$eq[[2]]$mse * diag( solve( crossprod( J ) ) ) )
Error in solve.default(crossprod(J)) : Lapack routine dgesv: system is
exactly singular
> print( se )

## gives the correct results becuase the J matrix is full rank
nlols.derivs <- eval( deriv( nlols$eq[[2]]$formula, c( "s0", "s1", "s2",
"s3" ), hessian=T ) )
J <- attr( nlols.derivs, "gradient" )
se <- sqrt( nlols$eq[[2]]$mse * diag( solve( crossprod( J ) ) ) )
print( se )


So I can do the following:

1) parse the matrix, column by column, and perform some operation to test to
see if the column is all zeros and tally up a character vector with the
names of the columns to obtain the terms in equation i

        eqn.terms <- vector()
        for( v in 1:length( est$estimate ) ) {
          j <- attr( eval( deriv( as.formula( eqns[[i]] ), names(
startvals ) ) ), "gradient" )
          if( qr( j[,v] )$rank > 0 ) {
            eqn.terms <- rbind( eqn.terms,
                               name <- names( est$estimate )[v] )
          }
        }

        derivs[[i]] <- deriv( as.formula( eqns[[i]] ), eqn.terms,
hessian=T )
        jacobian <- attr( eval( derivs[[i]] ), "gradient" )
        se <- sqrt( mse[i] * diag( solve( crossprod( jacobian ) ) ) )


2) decomopse the matrix (svd(crossprod(J)) or qr(crossprod(J)) ?
3) ask for help and see if there's a magic R function I can't find in the
code/faq/help docs/etc...

Sorry for the lame question but I'm not seeing an "elegant" solution...

Jeff.


---
Jeff D. Hamann
Forest Informatics, Inc.
PO Box 1421
Corvallis, Oregon USA 97339-1421
(office) 541-754-1428
(cell) 541-740-5988
jeff.hamann at forestinformatics.com
www.forestinformatics.com




More information about the R-help mailing list