[R] selfStart function problem
Lauri Mehtätalo
Lauri.Mehtatalo at metla.fi
Fri Jul 19 11:23:34 CEST 2002
Hello, list,
I am making a self-starting nonlinear function to model the relation of tree height (H)
and diameter (D) in a forest stand. The function I am trying is
H=a*exp[b*(D+5.8)^(-c)].
To calculate the initial estimates of the parameters, I linearized the formula by taking
logarithms and fixing the parameter c=1. Then I calculated the initial estimates of a
and b using lm() on the linearized form. The function works well if I am giving a single
variable as an dependent variable, but if I use some transformation of height (e.g.
(H-1.3)) as dependent variable, the expression in lm()-model is not right. That is
perhaps because of the mode of the object LHS is now "call", but I have not found
how to get the program handle it correctly. Has someone any idea? I am using R1.5.1
on Windows NT.
Lauri
My code is
> dhexp <- deriv(~ a*exp(b*(D+5.8)^(-c)), c("a","b","c"), function(D,a,b,c) {} )
>
> dhexpInit <- function(mCall, LHS, data) {
+ equation <- paste('log(',LHS,') ~ I((',mCall[["D"]],'+5.8)^-1.0)')
+ print(equation)
+ linear <- lm(as.formula(equation),data=data)
+ value <- c(exp(coef(linear)[1]), coef(linear)[2], 1.0)
+ names(value) <- mCall[c("a", "b", "c")]
+ value
+ }
>
> dhexp <- selfStart(dhexp, initial=dhexpInit)
>
> getInitial(height ~ dhexp(diameter,a,b,c), spruce)
[1] "log( height ) ~ I(( diameter +5.8)^-1.0)"
a b c
39.78098 -25.34938 1.00000
>
> getInitial(I(height-1.3) ~ dhexp(diameter,a,b,c), spruce)
[1] "log( I ) ~ I(( diameter +5.8)^-1.0)"
[2] "log( height - 1.3 ) ~ I(( diameter +5.8)^-1.0)"
Error in log(x) : Non-numeric argument to mathematical function
>
tree.nr time diameter height
3486 2 0 14.3 9.9
3487 2 1 14.6 10.3
3488 2 2 14.8 10.5
3489 2 3 15.2 11.0
3490 4 0 18.4 13.7
3491 4 1 19.1 14.3
3492 4 2 19.5 14.3
3493 4 3 20.5 15.1
3494 6 0 16.0 13.9
3495 6 1 16.4 14.4
- - - - - - - - - - - - -
Lauri Mehtätalo
Metsäntutkimuslaitos
Joensuun tutkimuskeskus
PL 68 (Yliopistokatu 7)
80101 JOENSUU
lauri.mehtatalo at metla.fi
p. 013-251 4110
- - - - - - - - - - - - -
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !) To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
More information about the R-help
mailing list