如何把一个复杂的公式转换成一个R公式？(How to put a complicated equat

我们拥有的树木作为预测结果和树高为因变量的直径。许多不同的方程存在此类型的数据，我们试图模拟一些人并比较结果。

但是，我们也不能弄清楚如何正确地把一个方程到相应的R formula格式。

该trees数据中设定R可以用来作为一个例子。

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()这似乎是逻辑。

我们怎样才能符合第二个方程我们的数据？

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)

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分的，但你应该用更好的猜测称，看球的参数可能是什么。

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 ，但需要正确的预测考虑到我们减去一个常数从LHS的事实：

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

因为我们已经使用了恒定的LHS（这几乎，但不相当，适合使用偏移的框架-我们只能用如果我们的配方为偏移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)