绘制二次回归并显示方程答案

【问题标题】：Plot quadratic regression with equation displayed绘制二次回归并显示方程
【发布时间】：2014-04-09 18:13:53
【问题描述】：

我正在尝试使用 for 循环创建一个包含 6 个图（3 行和 2 列）的 pdf 页面。我能够创建图，但我似乎无法自动为每个图添加回归线。

我正在尝试以下代码。

#Dummy data
Data1 <- data.frame(flow = c(8,8.5,6,7.1,9), SP_elev = c(20,11,5,25,50))
Data2 <- data.frame(flow = c(7,7.2,6.5,8.2,8.5), SP_elev = c(13,15,18,25,19))
Data3 <- data.frame(flow = c(2,3,5,7,9), SP_elev = c(20,25,28,30,35))
Data4 <- data.frame(flow = c(1,4,6,8,9), SP_elev = c(13,15,18,25,19))
Data5 <- data.frame(flow = c(1,4,6,8,9), SP_elev = c(13,15,18,25,19))
Data6 <- data.frame(flow = c(1,4,6,8,9), SP_elev = c(22,23,25,27,29))

#Create Vector list 
dataframes = list("Data1" = Data1, 
                  "Data2" = Data2, 
                  "Data3" = Data3,
                  "Data4" = Data4,
                  "Data5" = Data5,
                  "Data6" = Data6) # I gave up here

# open the PDF device
pdf(file="Dummy_Example.pdf", paper="letter", height=10, width=8)

#Create array of plots 
par(mfrow=c(3,2))

#plot a with regression model
for (i in dataframes) {

plot (i[,c('flow', 'SP_elev')], xlab=expression(paste("Discharge (", ft^3, "/s)",sep = "")), ylab= "Elevation (m)", tck=0.02, adj = 0.5)

#plot regression curve
fit2<-lm(i$SP_elev ~ i$flow + I(i$flow^2), data=i) 
pol2 <- function(x) fit2$coefficient[3]*x^2 + fit2$coefficient[2]*x + fit2$coefficient[1] 
curve(pol2, lwd=1, add=T)
}

# close the PDF device
dev.off()

当我单独制作一个绘图时，我得到了回归曲线代码，但当我尝试自动化它时它似乎不起作用。

另外我还想在图上画出回归曲线的方程。

谢谢，

配音

【问题讨论】：

我重新编辑了我的原始帖子。这就是你要找的东西吗？
此时我可以解决不使用 dput() 输出并且不勾选早期答案的人提出的这个问题......或者我可以观看有关热力学的 CD 光盘。
非常感谢任何帮助。我对 R 很陌生。dput() 到底是什么？
@kevinWright 你能再看看我原来的帖子吗？我用一个可重现的例子添加了一些虚拟数据。
@KevinWright 你能把它作为答案转发吗

标签： r plot batch-processing nonlinear-functions

【解决方案1】：

我编辑了上面的代码以将方程包含在图中。更新：现在更漂亮的方程。

#Dummy data
Data1 <- data.frame(flow = c(8,8.5,6,7.1,9), SP_elev = c(20,11,5,25,50))
Data2 <- data.frame(flow = c(7,7.2,6.5,8.2,8.5), SP_elev = c(13,15,18,25,19))
Data3 <- data.frame(flow = c(2,3,5,7,9), SP_elev = c(20,25,28,30,35))
Data4 <- data.frame(flow = c(1,4,6,8,9), SP_elev = c(13,15,18,25,19))
Data5 <- data.frame(flow = c(1,4,6,8,9), SP_elev = c(13,15,18,25,19))
Data6 <- data.frame(flow = c(1,4,6,8,9), SP_elev = c(22,23,25,27,29))

#Create Vector list 
dataframes = list("Data1" = Data1, 
                  "Data2" = Data2, 
                  "Data3" = Data3,
                  "Data4" = Data4,
                  "Data5" = Data5,
                  "Data6" = Data6) # I gave up here

# open the PDF device
pdf(file="Dummy_Example.pdf", paper="letter", height=10, width=8)

#Create array of plots 
par(mfrow=c(3,2))

#plot a with regression model
for (i in dataframes) {

plot (SP_elev ~ flow, data=i,
      xlab=expression(paste("Discharge (", ft^3, "/s)",sep = "")),
      ylab= "Elevation (m)", tck=0.02, adj = 0.5)

#plot regression curve
fit2<-lm(SP_elev ~ flow + I(flow^2), data=i) 
pol2 <- function(x) fit2$coefficient[3]*x^2 + fit2$coefficient[2]*x + fit2$coefficient[1] 
curve(pol2, lwd=1, add=T, col="blue")

xm <- min(i$flow)
ym <- max(i$SP_elev)

co <- signif(coef(fit2),3)
text(xm, ym,
 bquote(y==.(co[3])*x^2 + .(co[2])*x + .(co[1])),
 adj=c(0,1))


}
dev.off()

【讨论】：

你知道为什么expression(x^2) 不会在上标中显示 2 吗？当我运行代码时，它仍然显示'^'。
更新为真正的上标。