如何将一个复杂的方程式放到R公式中? [英] How to put a complicated equation into a R formula?

问题描述

我们将树木的直径作为预测变量，并将树木的高度作为因变量.对于此类数据，存在许多不同的方程式，我们尝试对其中的一些进行建模并比较结果.

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.

但是，我们无法弄清楚如何正确地将一个方程式转换为相应的R formula格式.

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

在R中设置的trees数据可以用作示例.

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  

估计

系数a，b和c，这是我们感兴趣的.

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

现在有问题的方程式:

尝试像这样拟合它:

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

出现错误:

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?

这不会给出错误:

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  

外部I()中的整个项仅给出一个系数，这似乎是逻辑.

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?

推荐答案

假设您使用的是nls，则R公式可以使用普通的R函数H(a, b, c, D)，因此该公式可以只是h ~ H(a, b, c, dbh)，并且可以正常工作:

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

以图形方式显示输出:

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

这篇关于如何将一个复杂的方程式放到R公式中?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！