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
