如何将复杂的方程式放入 R 公式中？

Question

我们以树的直径为自变量，以树高为因变量。这类数据有多种不同的方程式可用于建模和比较结果。但是，我们无法确定如何将其中一种方程正确地放入相应的R公式格式中。可以使用R中的“trees”数据集作为示例。

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

首先，这是一个看起来运作良好的方程示例：

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是我们感兴趣的估计值。

现在来看一下问题所在的方程式：

尝试这样适配：

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

我猜这是因为/被解释为嵌套模型而不是算术运算符？

这不会产生错误：

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

但结果并非我们想要的：

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

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

在外部的I()中，只给出了整个项的一个系数，这似乎是合理的。

我们如何将第二个方程拟合到我们的数据中？

- donodarazao

Answer 1

假设您正在使用 nls，R公式可以使用普通的R函数H(a, b, c, D)，因此公式可以只写为h ~ H(a, b, c, dbh)，这样就可以实现：

# 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)

- G. Grothendieck

---

我会将这个答案标记为正确的，因为a）它包括如何估计起始值，b）使用普通的R函数使我们能够非常容易地拟合其他非线性函数，c）它绘制了结果。谢谢！ - donodarazao

Answer 2

你有几个问题。(1) 在form2的分母中缺少括号（R无法知道您想要在分母中添加常数a，或者在哪里放置任何参数），而且更为严重的是：(2) 您的第二个模型不是线性的，因此lm无法工作。

修复（1）很容易：

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

修复（2），虽然有很多方法可以估计非线性模型的参数，但nls（非线性最小二乘）是一个不错的起点：

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

你需要为nls中的参数提供起始猜测值。我只是选了1，但你应该使用更好的猜测值，以大致确定参数可能的取值范围。

- Gregor Thomas

---

感谢您的回答！要发现那些问题，甚至要找到解决方案，我们可能需要花费很长时间。 - donodarazao

Answer 3

编辑: 已修复，不再错误地使用偏移量...

一个补充@shujaa的答案:

您可以将问题转换为

H = 1.3 + D^2/(a+b*D+c*D^2)

到

1/(H-1.3) = a/D^2+b/D+c

这通常会破坏模型的假设（即，如果H服从具有恒定方差的正态分布，则1/（H-1.3）不会）。但是，我们还是尝试一下：

data(trees)
df <- transform(trees,
            h=Height * 0.3048,   #transform to metric system
            dbh=Girth * 0.3048 / pi   #transform tree girth to diameter
            )
lm(1/(h-1.3) ~ poly(I(1/dbh),2,raw=TRUE),data=df)

## Coefficients:
##                    (Intercept)  poly(I(1/dbh), 2, raw = TRUE)1  
##                       0.043502                       -0.006136  
## poly(I(1/dbh), 2, raw = TRUE)2  
##                       0.010792  

这些结果通常足以获得nls拟合的良好起始值。但是，您可以通过使用链接函数来实现某些非线性形式，从而比这更好，glm可以做到这一点。具体地说，

(fit2 <- glm(h-1.3 ~ poly(I(1/dbh),2,raw=TRUE),
             family=gaussian(link="inverse"),data=df))

## Coefficients:
##                    (Intercept)  poly(I(1/dbh), 2, raw = TRUE)1  
##                       0.041795                       -0.002119  
## poly(I(1/dbh), 2, raw = TRUE)2  
##                       0.008175  
## 
## Degrees of Freedom: 30 Total (i.e. Null);  28 Residual
## Null Deviance:       113.2 
## Residual Deviance: 80.05     AIC: 125.4 
## 

你可以看到结果与线性拟合大致相同，但并非完全相同。

pframe <- data.frame(dbh=seq(0.8,2,length=51))

我们使用 predict，但需要修正预测结果，以考虑从左侧减去常数的事实。

pframe$h <- predict(fit2,newdata=pframe,type="response")+1.3
p2 <- predict(fit2,newdata=pframe,se.fit=TRUE) ## predict on link scale
pframe$h_lwr <- with(p2,1/(fit+1.96*se.fit))+1.3
pframe$h_upr <- with(p2,1/(fit-1.96*se.fit))+1.3
png("dbh_tmp1.png",height=4,width=6,units="in",res=150)
par(las=1,bty="l")
plot(h~dbh,data=df)
with(pframe,lines(dbh,h,col=2))
with(pframe,polygon(c(dbh,rev(dbh)),c(h_lwr,rev(h_upr)),
      border=NA,col=adjustcolor("black",alpha=0.3)))
dev.off()

由于我们在等式左侧使用了常量（这几乎符合使用“偏移”的框架，但并不完全符合—如果我们的公式是 1/H-1.3=a/D^2 +...，也就是说常量调整是在链接（反转）比例尺上而不是原始比例尺上），这与ggplot 的geom_smooth 框架并不完全匹配。

library("ggplot2")
ggplot(df,aes(dbh,h))+geom_point()+theme_bw()+
   geom_line(data=pframe,colour="red")+
   geom_ribbon(data=pframe,colour=NA,alpha=0.3,
             aes(ymin=h_lwr,ymax=h_upr))

ggsave("dbh_tmp2.png",height=4,width=6)

- Ben Bolker