R help - Jul 2002 - selfStart function problem

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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._

"Lauri Meht?talo" <Lauri.Mehtatalo at metla.fi> writes:
> 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
> +     }
Although I am responding quite late to this message I hope this will
help.  I did not see other responses.

The recommended way of handling LHS inside the start function is as

              xy <- sortedXyData(mCall[["D"]], LHS, data)

That will provide you with a data frame with variables 'x' and
'y'
that result from evaluating the expressions for the D argument and the
left hand side.

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._