R中散点图的回归直线和拟合曲线

Question

我有一组树木的高度和直径的数据。我想找出它们之间的回归关系，并绘制出来。例如，我想尝试a * DIAMETER + b * DIAMETER^2 + C，并以散点图的形式显示其曲线。通过下面的说明，我达到了几条线，但我只想要一条与开发模型相关的趋势线。我该怎么办？

setwd('D:\\PhD\\Data\\Field Measurments\\Data Analysis\\')

dat1 = read.table('Fagus.csv', header = TRUE, sep =',')


# fit a non-linear regression

Height = dat1$Height

Diameter = dat1$Diameter 

plot(Diameter, Height, main="Height Curve", xlab="Diameter", ylab="Height", pch=19)

nls1 <- nls(Height ~ a*(Diameter)^2+b*Diameter+c, data = dat1, start = list(a =a, b=b,c=c), algorithm="port")

lines(fitted(nls1) ~ Diameter, lty = 1, col = "red") # solid red line

上面的说明对我来说是错误的吗？

Answer 1

如上所述，您不应该将系数放入您的公式中。尝试：

nls1 <- nls(Height ~ I(Diameter^2) + Diameter, data = dat1,  algorithm="port")

关于I(Diameter ^2)

“为了避免这种混淆，可以使用函数I()将模型公式中运算符在算术意义上使用的部分括起来。例如，在公式y ~ a + I(b+c)中，术语b+c应解释为b和c的和。”~ formula{stats}文档

我没有运行剩下的部分(在移动设备上)，但是你的代码乍一看看起来还不错。

Answer 2

这里似乎对线性模型和非线性模型存在误解。线性模型的系数是线性的。非线性模型则不是。模型在预测变量(在你的例子中是直径)中是否是线性的是无关的。因此，在您的示例中，是表单的模型：

高度=a*直径+b*直径^2+c

是一个线性模型。您不需要使用nls(...)。您可以通过两种方法之一指定模型公式，这两种方法都会产生相同的结果：

Height~Diameter + I(Diameter^2)

或

Height~poly(Diameter,2,raw=TRUE)

第二种形式使用poly(...)函数创建2阶多项式。raw=T告诉poly(...)生成原始多项式，而不是正交多项式(默认)。第一种形式比较简单，除非你想要阶数大于2的多项式。下面是一个使用这两种形式的例子。

set.seed(1)    # for reproducible example
df <- data.frame(Diameter=sample(1:50,50))
df$Height <- with(df,2*Diameter + .5*Diameter^2 + 4 + rnorm(50,sd=30))

fit <- lm(Height~Diameter + I(Diameter^2),df)
summary(fit)
# ...
# Coefficients:
#               Estimate Std. Error t value Pr(>|t|)    
# (Intercept)   -6.85088   12.26720  -0.558  0.57917    
# Diameter       3.31030    1.10964   2.983  0.00451 ** 
# I(Diameter^2)  0.47717    0.02109  22.622  < 2e-16 ***

fit.poly<- lm(Height~poly(Diameter,2,raw=TRUE),df)
summary(fit.poly)
# Coefficients:
#                                Estimate Std. Error t value Pr(>|t|)    
# (Intercept)                    -6.85088   12.26720  -0.558  0.57917    
# poly(Diameter, 2, raw = TRUE)1  3.31030    1.10964   2.983  0.00451 ** 
# poly(Diameter, 2, raw = TRUE)2  0.47717    0.02109  22.622  < 2e-16 ***

绘制数据和趋势曲线：

df$pred <- predict(fit)
with(df,plot(Height~Diameter))
with(df[order(df$Diameter),],lines(pred~Diameter,col="red",lty=2))

​

​

Answer 3

您的问题出在start=参数上。您需要提供a、b和c参数的实际值。下面是一个可重复使用的示例

#sample data
dat<-data.frame(Diameter = runif(50, 2, 6))
dat<-transform(dat,Height=2*Diameter + .75 * Diameter^2 +4 + rnorm(50))
dat<-dat[order(dat$Diameter), ]

#now fit the model
mynls<-nls(Height ~ a*I(Diameter^2) + b*Diameter + c, dat,
    start=list(a=1, b=1, c=1), algorithm="port")

注意我们如何为每个系数设置默认值1。您可以设置您认为最合适的任何内容。以及如何使用拟合结果绘制原始值

plot(Height~Diameter,dat, main="Height Curve", 
    xlab="Diameter", ylab="Height", pch=19)
lines(fitted(mynls)~ dat$Diameter, col="red")

这给了我们

​

​