How to put a complicated equation into a R formula?

Question

We have the diameter of trees as the predictor and tree height as the dependent variable. A number of different equations exist for this kind of data and we try to model some of them and compare the results.

However, we we can't figure out how to correctly put one equation into the corresponding R formula format.

The trees data set in R can be used as an example.

data(trees)
df <- trees
df$h <- df$Height * 0.3048   #transform to metric system
df$dbh <- (trees$Girth * 0.3048) / pi   #transform tree girth to diameter

First, the example of an equation that seems to work well:

form1 <- h ~ I(dbh ^ -1) + I( dbh ^ 2)  
m1 <- lm(form1, data = df)
m1

Call:
lm(formula = form1, data = df)

Coefficients:
(Intercept)    I(dbh^-1)     I(dbh^2)  
27.1147      -5.0553       0.1124  

Coefficients a, b and c are estimated, which is what we are interested in.

Now the problematic equation:

Trying to fit it like this:

form2 <- h ~ I(dbh ^ 2) / dbh + I(dbh ^ 2) + 1.3

gives an error:

m1 <- lm(form2, data = df)
Error in terms.formula(formula, data = data) 
invalid model formula in ExtractVars

I guess this is because / is interpreted as a nested model and not an arithmetic operator?

This doesn't give an error:

form2 <- h ~ I(I(dbh ^ 2) / dbh + I(dbh ^ 2) + 1.3)
m1 <- lm(form2, data = df)

But the result is not the one we want:

m1
Call:
lm(formula = form2, data = df)

Coefficients:
(Intercept)  I(I(dbh^2)/dbh + I(dbh^2) + 1.3)  
19.3883                            0.8727  

Only one coefficient is given for the whole term within the outer I(), which seems to be logic.

How can we fit the second equation to our data?

Answer 1

Assuming you are using nls the R formula can use an ordinary R function, H(a, b, c, D), so the formula can be just h ~ H(a, b, c, dbh) and this works:

# use lm to get startingf values
lm1 <- lm(1/(h - 1.3) ~ I(1/dbh) + I(1/dbh^2), df)
start <- rev(setNames(coef(lm1), c("c", "b", "a")))

# run nls
H <- function(a, b, c, D) 1.3 + D^2 / (a + b * D + c * D^2)
nls1 <- nls(h ~ H(a, b, c, dbh), df, start = start)

nls1 # display result

Graphing the output:

plot(h ~ dbh, df)
lines(fitted(nls1) ~ dbh, df)

Answer 2

You've got a couple problems. (1) You're missing parentheses for the denominator of form2 (and R has no way to know that you want to add a constant a in the denominator, or where to put any of the parameters, really), and much more problematic: (2) your 2nd model isn't linear, so lm won't work.

Fixing (1) is easy:

form2 <- h ~ 1.3 + I(dbh^2) / (a + b * dbh + c * I(dbh^2))

Fixing (2), though there are many ways to estimate parameters for a nonlinear model, the nls (nonlinear least squares) is a good place to start:

m2 <- nls(form2, data = df, start = list(a = 1, b = 1, c = 1))

You need to provide starting guesses for the parameters in nls. I just picked 1's, but you should use better guesses that ballpark what the parameters might be.