【原】【R分享|实战】 科白君教你如何给散点图添加最优拟合曲线和相关性~

   

“ 每期小积累，分析学得美。”   --科白君

"R实战"专题·第3篇

  编辑 | 科白维尼

  5854字 | 10分钟阅读

本期推送内容

只有空荡荡的散点图，似乎不能满足我们。散点图能表示线性回归关系，但需要丰富信息和内容，例如，添加拟合曲线、相关性等等，让散点图可读性更强。然而，拟合的回归模型类型特别多，并不是所有数据都呈线性关系，也有可能是非线性。线性主要包括一次、二次线性关系，对数关系等；非线性包括指数函数，幂函数等。这意味着每组数据时可以选择模型的，我们可以通过计算评价回归模型的指标来选择最优的模型，然后再添加拟合曲线(公式、解释量及显著性)。本期与大家分享以下内容：1)如何添加散点图的拟合曲线;2)如何计算数据的R、p值;3)如何通过计算AIC等指标来选择最优的拟合模型，并计算模型的R^2和p值以及添加最优模型的公式，希望大家有所收获。

本期的推送主要分为三部分：

1）详细讲解如何在散点图添加拟合曲线的一些函数；

2）如何选择最优的拟合模型；

3）如何定义线性与非线性拟合模型，并添加模型的公式。

本期代码有点长，大家可以找个空闲的时候慢慢学习~

首先，加载需要的R包

# 清空环境变量
rm (list = ls ())

#1 加载所需要的R包
# 如果R包还没安装的可以参考上期推文: R包的安装与使用教程进行安装~
pacman : : p_load (ggplot2, ggpubr, ggforce, ggpointdensity, reshape2, scales, grid,
 dplyr, viridis, export)

然后，加载R语言自带的数据集"mpg"，并查看数据结构

data (mpg)
# 了解数据结构
# 该数据集为 234行×11列 tribble数据
# tribble主要用来改进R内置数据框(data.frame)存在的许多问题
str (mpg)

展示：

我们可以用 as.data.frame（）转换数据结构

mpg <- as.data.frame (mpg)
str (mpg)
head (mpg)

展示：

此外，如果不喜欢系统默认的行名可以用 fix（）进行修改

fix (mpg)
head (mpg)

展示修改过程：

直接点任何值进行修改

进入第一部分的讲解：

1.1 绘制基础散点图

# 主要用 geom_point () 函数
p1 <- ggplot (mpg, aes (displ, hwy)) +
  geom_point (size=2)
p1

1.2 我们添加默认的拟合曲线，用 geom_smooth（）/ stat_smooth和（） 函数

# 关于拟合曲线的添加，我们可以先了解一下这个函数 ?geom_smooth
# geom_smooth加入拟合曲线，默认loess, se=TRUE表示加入置信区间
# span控制loess平滑的平滑量，较小的数字产生波动线，较大的数字产生平滑线
p2 <- p1 +
  geom_smooth ()
p2

# 提示出现：默认使用局部加权回归 loess (Locally Weighted Scatterplot Smoothing)  

# loess主要思想是取一定比例的局部数据，在这部分子集中拟合多项式回归曲线

# 这样就可以观察到数据在局部展现出来的规律和趋势

# 曲线的光滑程度与数据比例有关：比例越少，拟合越不光滑，反之越光滑

1.3 使用span(跨度)调整默认平整度的“波动幅度”

p3 <- p1+
  geom_smooth (span = 0.3)
p3

1.4 如果不想用默认loess，可以使用其他回归模型

p4 <- p1 +
  geom_smooth (method = "lm", # lm 表示拟合线性模型
   se = FALSE) # se = FALSE 表示删除置信区间
p4

p5 <- p1+
  geom_smooth (method = "lm",
  formula = y ~ splines : : bs (x, 3),
  se = FALSE)
p5

# splines 包，我们可以了解一下 ?splines

# splines包中的样条函数，df=3，样条具有三个基函数

# 关于样条，常用的有两类，一类是多项式样条，另一类是光滑样条

# 具体详情参考http://www./article-5298-1.html

1.5 修改添加的拟合曲线和置信区间的颜色等

p6 <- p1 +
  geom_smooth (method = "lm",
     se = TRUE, # T or F 表示是否添加置信区间
  colour = "#00A5FF", # colour 表示拟合曲线的颜色
  fill = "#00A5FF", # fill 表示置信区间的颜色
     alpha = 0.2) # alpha 表示透明度
p6

1.6 平滑曲线自动拟合于每个组(由aes定义)

p7 <- ggplot(mpg, aes(displ, hwy, colour = fl)) +
  geom_point (size=2)+
  geom_smooth (method = "lm",
  se = FALSE)
p7

或者 用facet_wrap（）函数将不同组进行分面

p7_2 <- p7_1 +
  stat_cor (aes (colour = fl),
   method = "pearson", #用pearson相关性分析，可以自行更改方法
   label.x = 3) #在 x = 4 的地方表示相关性信息
p7_2

这里，我们添加相关性分析的信息

# 这里用到 stat_cor () 函数，必须加载 ggpubr 包
# 我们也可以先了解一下该函数?stat_cor
p7_2 <- p7_1 +
  stat_cor (aes (colour = fl),
   method = "pearson", #用pearson相关性分析，可以自行更改方法
   label.x = 3) #在 x = 4 的地方表示相关性信息
p7_2

第二部分，通过计算选择出最优的拟合模型

2.1  利用 "mpg" 数据集中的两列数据构建一个新的数据集

mpg_new <- data.frame (mpg$hwy, mpg$displ) # data.frame 进行构建数据框
head (mpg_new)
colnames (mpg_new)[1]<- "x" # 更改新构建的数据第一列的名称为 "x"
colnames (mpg_new)[2] <- "y" # 更改新构建的数据第二列的名称为 "y"
head (mpg_new) # 当然也可以用上面讲到的fix() 函数直接更改

2.2 我们先把数据描绘出来

p8 <- ggplot (mpg_new, aes (x, y))+
  geom_point (size = 2)+
  stat_smooth (method ="loess",
  formula = y ~ x,
  se = T,
  colour = "#00A5FF",
  fill = "#00A5FF",
  alpha = 0.2,
  span = 0.5) +
  stat_cor (method ="pearson",
   label.x = 30) +
  scale_x_continuous (limits = c(10, 50))+
  scale_y_continuous (limits = c(0, 8))
p8

2.3 设定各类拟合模型，这里主要讲两类：线性模型与非线性模型

models <- list (lm (y ~ x, data = mpg_new), # lm 表示线性模型
   lm (y ~ I(1/x), data = mpg_new),
   lm (y ~ log(x), data = mpg_new),
   nls (y ~ I(x^2 *a)+ b* x + c, data = mpg_new, start = list (a= 1, b = 1, c= 0)), #nls 表示非线性模型
   nls (y ~ I(1/x * a)+ b * x, data = mpg_new, start = list (a = 1, b = 1)),
   nls (y ~(a + b*log(x)), data = mpg_new, start = setNames (coef (lm (y ~ log(x), data = mpg_new)), c("a", "b"))),
   nls (y ~ I(exp(1)^(a + b * x)), data = mpg_new, start = list (a= 0, b= 1)),
   nls (y ~ I(1/x * a)+b, data = mpg_new, start = list (a= 1, b= 1)),
   nls (y ~ a*I(x^b), data = mpg_new, start = list (a= 2, b= 1.5))
) # I() 函数可以避免这种混淆将模型公式中使用运算符的部分括起来
  # 例如，在公式y ~ a + I(b+c)中，b+c项被解释为b和c的和
models[[1]] # 查看第一个公式的信息， [[X]] x为1-9共9个公式，可以直接查看对应的公式信息
# models 查看所有模型的参数

2.4 简单绘制图形，可以先查看不同模型拟合曲线的大致情况

p8_1 <- ggplot (mpg_new, aes (x, y)) + geom_point (size = 3.5) +
  stat_smooth (method = "lm", formula = as.formula(models[[1]]), size = 1, se = FALSE, colour ="black") +
  stat_smooth (method = "lm", formula = as.formula(models[[2]]), size = 1, se = FALSE, colour ="blue") +
  stat_smooth (method = "lm", formula = as.formula(models[[3]]), size = 1, se = FALSE, colour ="yellow") +
  stat_smooth (method = "nls", formula = as.formula(models[[4]]), size = 1, se = FALSE, colour ="purple", linetype = 2) +# linetype=2 为 dashed 虚线
  stat_smooth (method = "nls", formula = as.formula(models[[5]]), size = 1, se = FALSE, colour = "red", linetype = 2) +
  stat_smooth (method = "nls", formula = as.formula(models[[6]]), size = 1, se = FALSE, colour ="green", linetype = 2) +
  stat_smooth (method = "nls", formula = as.formula(models[[7]]), size = 1, se = FALSE, colour ="violet")+
  stat_smooth (method = "nls", formula = as.formula(models[[8]]), size = 1, se = FALSE, colour ="orange", linetype = 2) +
  stat_smooth (method ="nls", formula = as.formula(models[[9]]), size = 1, se = FALSE, colour ="grey", linetype = 2)
p8_1 

2.5 通过计算AIC值，选出最优模型

# AIC值越小，拟合模型更好
library (AICcmodavg); library (plyr); library (stringr) #读取所需要用的R包 ，library (pkg1) + ; library (pkg2)  但是只能放在同一行才不会出错
ldply (models, function(mod){ data.frame (AICc = AICc(mod), AIC = AIC(mod), model = deparse(formula(mod))) })

第三部分，选出最优方程后，我们定义拟合模型

3.1 根据AIC值，我们选出了最优模型为 models[[4]]：y ~ a*I(x^2) + b*x+c

# 接下来我们定义函数，计算得出后面需要要添加的公式
# 先定义 lm 线性函数
lm_eqn = function (mpg_new){
  m <- lm (y ~ I(x^2)+ x,mpg_new)
  eq <- substitute(italic(y) == ~c + ~b %.% italic(x)+ ~a %.% italic(x)^2*","~~italic(r)^2~"="~r2,
    list (a = format (coef(m)[2], digits = 3), # ~ 表示一个空格
   b = format (coef(m)[3], digits = 3), # ","表示，；"=" 表示 =
   c = format (coef(m)[1], digits = 3), # digits 表示保留三位有效数字
   r2 = format (summary(m)$r.squared, digits = 3))) # substitute() 函数能够解释模型的表达式
  as.character (as.expression (eq))
}  # 这里需要注意的是，我发现如果按a b c (即123)的顺序出来的公式与上面模型计算出来的参数有差异
# 因此我建议以模型计算出来的参数为准，防止定义时，系统计算出错


# 非线性模型函数的定义方法

# 非线性模型R^2计算 参考http://blog.sina.com.cn/s/blog_78c5f0530101c7w1.html
# 我已经把这部分内容用# 符号隐藏了，删除 # 尝试
# lm_eqn = function(mpg_new){
# mod <- nls (y~ a*I(x^b),start = list(a=2, b= 1.5), data = mpg_new)
# a=0
# b=0
# for(i in 1:100){ # 按理这个循环要循环234次，但是100次的时候就已经是负值了
# ymean=mean(y) # 结合上面计算的AIC值，意味着这个模型不可用
# a=a+(y[i]-fitted(mod)[i])^2
# b=b+(y[i]-ymean)^2
# }
# r2=1-a/b # 由于非线性的R^2 不能够直接用summary访问， 只能用循环公式来计算
# eq <- substitute(italic(y) == ~a %.% italic(x) ^ b*","~~italic(r)^2~"="~r2,
# list(a = format(coef(mod)[1], digits = 3),
# b = format(coef(mod)[2], digits = 3),
# r2 = format(r2,digits = 3)))
# as.character(as.expression(eq))
# }

3.2 然后把上面定义的函数添加到图里

p8_2 <- ggplot (mpg_new, aes (x, y)) +
  geom_point (shape = 19,alpha = 0.3,stroke = 0,size = mpg_new$y)+ # 这里选择了16号形状 可以用show_point_shapes() 访问形状类别
  stat_smooth (method = "lm",
  formula = y ~ I(x^2) + x,
  se = T,
  colour = "#00A5FF",
  fill = "#00A5FF",
  alpha = 0.2) +
  scale_x_continuous (limits = c(10, 50)) +
  scale_y_continuous (limits = c(0, 8)) +
  theme_bw () +
  theme (panel.border = element_blank (),
  panel.grid = element_blank (),
  axis.line = element_line (colour = "black")) +
  theme (axis.text = element_text (colour = 'black',size = 12),
   axis.title = element_text(size = 14,colour = 'black')) + # ' ' 与 "" 用法相同
  annotate ("text",  # annotate () 函数为添加文本信息
   x = 10, # 表示公式添加的位置 x轴
   y = 7.5, # 表示公式添加的位置 y轴
   label = lm_eqn (mpg_new), # 添加的公式为上面定义的内容 (包括计算)
   hjust = 0,
   size = 6,
   parse = TRUE)+ # parse = T 表示能够解析出模型公式，如果为F 则识别的公式与正常的表达不同
  labs (x ="Soil depth (cm)", y ="Soil organic carbon (g/kg)")
p8_2
graph2ppt (file="p9_2.ppt", append = FALSE)

最后经过调整的结果：

额外福利，我找到了另外两种能够计算回归模型参数的R包，1) basicTrendline包

library("basicTrendline") # 数据还是用的新构建的dataset
trendline (mpg_new$x, mpg_new$y, model="line2P", ePos.x = "topleft", summary=TRUE, eDigit=5)
trendline (mpg_new$x, mpg_new$y, model="line3P", CI.fill = FALSE, CI.color = "black", CI.lty = 2)
trendline (mpg_new$x, mpg_new$y, model="log2P", ePos.x= "top", linecolor = "red", CI.color = NA)
trendline (mpg_new$x, mpg_new$y, model="exp2P", CI.fill = FALSE, CI.color = "black", CI.lty = 2)
trendline (mpg_new$x, mpg_new$y, model="exp3P", xname="T", yname = paste(delta^15,"N") ,
 yhat=FALSE, Rname=1, Pname=0, ePos.x = "bottom")
# 如果出现报错，很可能是因为这个数据用这个模型时不符合条件

简单介绍其中一种情况：

2)  ggpmisc包，利用 stat_poly_eq（）函数来完成

# 先定义一个拟合模型
my_formular <- y ~ I(x^2) + x

# 添加R^2
ggplot(mpg_new, aes(x, y)) +
  geom_point() +
  geom_smooth(method = "lm", formula = my_formular) +
  stat_poly_eq(formula = my_formular, rr.digits = 4, parse = TRUE)

# 添加各种标签 包括拟合模型的公式，R^2，p值，AIC值等
ggplot(mpg_new, aes(x, y)) +
  geom_point() +
  geom_smooth(method = "lm", formula = my_formular) +
  stat_poly_eq(aes(label = paste(stat(eq.label), # 可以? stat_poly_eq 里面有对应的解释
 stat(adj.rr.label), # 可以根据自己需要添加对应的信息
 stat(p.value.label),
 stat(AIC.label),
 sep = "*\", \"*")), # 这行代码是识别和分割x，如果删除则无法正常识别
 coef.digits = 4, #保留4位有效数字
   formula = my_formular, parse = TRUE)
# 但是能够发现，系统计算出来的拟合公式是错误的
# 正确的应该是 y ~ 0.003*x^2 -0.335*x +9.310

# 更改以下代码，添加文本信息，添加上面我们定义的函数
ggplot(mpg_new, aes(x, y)) +
  geom_point() +
  geom_smooth(method = "lm", formula = my_formular) +
  stat_poly_eq(geom = "text", label.x = 20, label.y = 0, hjust = 0,
   formula = my_formular, parse = TRUE, label = lm_eqn(mpg_new))

附上完整的代码：

#######################################
# 添加最优拟合模型曲线和相关性的散点图
#######################################

# 清空环境变量
rm(list=ls())

#1 加载所需要的R包
# p_load 需要下载 pacman 包才能运行
# 如果R包还没安装的可以参考上期推文: R包的安装与使用教程进行安装~
pacman::p_load (ggplot2, ggpubr, ggforce, ggpointdensity,reshape2,scales,grid,
 dplyr, viridis, export)

#2 加载R语言自带的数据集“mpg”
data (mpg)

# 了解数据结构
# 该数据集为 234行×11列 tribble数据
# tribble主要用来改进R内置数据框(data.frame)存在的许多问题
str (mpg)

# 不习惯的话，可转换数据结构，as.data.frame()进行转化
mpg <- as.data.frame(mpg)
str (mpg)
head (mpg)

# 如果不喜欢系统默认的行名可以用fix()进行修改
fix (mpg)
head (mpg)

#3 添加散点图的拟合曲线，包括R、p及拟合曲线的类型和公式
#1) 绘制基础散点图
# 主要用 geom_point() 函数
p1 <- ggplot (mpg, aes (displ, hwy)) +
  geom_point (size=2) 
p1 

#2) 添加默认的拟合曲线，用 geom_smooth() 函数
# 关于拟合曲线的添加，我们可以先了解一下这个函数 ?geom_smooth
# geom_smooth加入拟合曲线，默认loess,se=TRUE表示加入置信区间
# span控制loess平滑的平滑量，较小的数字产生波动线，较大的数字产生平滑线
p2 <- p1 +
  geom_smooth()
p2
# 提示出现表示：默认使用局部加权回归 loess (Locally Weighted Scatterplot Smoothing)
# loess主要思想是取一定比例的局部数据，在这部分子集中拟合多项式回归曲线
# 这样就可以观察到数据在局部展现出来的规律和趋势
# 曲线的光滑程度与选取数据比例有关：比例越少，拟合越不光滑，反之越光滑

#3) 使用span(跨度)调整默认平整度的“波动幅度”
p3 <- p1+
  geom_smooth (span = 0.3)
p3

#4) 如果不想用默认loess，可以使用其他模型公式
p4 <- p1 +
  geom_smooth (method = "lm", # lm 表示拟合线性模型
   se = FALSE) # se = FALSE 表示删除置信区间
p4

p5 <- p1+
  geom_smooth(method = "lm",
  formula = y ~ splines::bs(x, 3),
  se = FALSE)
p5
# splines 包，我们可以了解一下 ?splines
# splines包中的样条函数，df=3，样条具有三个基函数
# 关于样条，常用的有两类，一类是多项式样条，另一类是光滑样条
# 具体详情参考http://www./article-5298-1.html

#5) 修改添加的拟合曲线和置信区间的颜色等
p6 <- p1 +
  geom_smooth(method = "lm",
  se = TRUE, # T or F 表示是否添加置信区间
  colour = "#00A5FF", # colour 表示拟合曲线的颜色
  fill = "#00A5FF", # fill 表示置信区间的颜色
  alpha = 0.2) # alpha 表示透明度
p6

#6) 平滑曲线自动拟合于每个组(由aes定义)
p7 <- ggplot(mpg, aes(displ, hwy, colour = fl)) +
  geom_point (size=2) +
  geom_smooth (method = "lm",
  se = FALSE)
p7

# 可以利用facet_wrap()函数将不同组进行分面
p7_1 <- ggplot(mpg, aes(displ, hwy, colour = fl)) +
  geom_point (size=2) +
  geom_smooth (method = "lm",
   se = FALSE) +
  facet_wrap(~fl,scales = "free")
p7_1

# 添加相关性信息
# 这里用到 stat_cor() 函数，必须加载 ggpubr 包
# 我们也可以先了解一下该函数?stat_cor
p7_2 <- p7_1 +
  stat_cor (aes(colour = fl), 
   method = "pearson", #用pearson相关性分析，可以自行更改方法
   label.x = 3) #在 x = 4 的地方表示相关性信息
p7_2

#7) 利用数据对各类模型进行拟合，再通过计算AIC值选择最优的拟合模型
# 为什么要计算AIC来选择最优模型？这里参考 https://www.baidu.com/link?url=TYvKoyxcrhDDN3nSR6hP6NdgvveMvGi5mH9cP5m2M7bmfCEbLB01rnzbnRIWR1EKJ6pjMUVdPlezPNJjQHVcK_&wd=&eqid=90468e34000e3b880000000360729ac2
# 利用 "mpg" 数据集中的两列数据构建一个新的数据集
mpg_new <- data.frame(mpg$hwy, mpg$displ, mpg$fl) # data.frame 进行构建数据框
head (mpg_new)
colnames (mpg_new)[1] <- "x" # 更改新构建的数据第一列的名称为 "x"
colnames (mpg_new)[2] <- "y" # 更改新构建的数据第二列的名称为 "y"
colnames (mpg_new)[3] <- "group" # 更改新构建的数据第二列的名称为 "group"

head (mpg_new) # 当然也可以用上面讲到的fix() 函数直接更改

# 我们先把数据描绘出来
p8 <- ggplot (mpg_new, aes(x, y))+
  geom_point (size=2)+
  stat_smooth (method = "loess",
  formula = y ~ x,
  se = T,
  colour = "#00A5FF",
  fill = "#00A5FF",
  alpha = 0.2,
  span = 0.5) +
  stat_cor (method = "pearson",
   label.x = 30) +
  scale_x_continuous (limits = c(10,50))+
  scale_y_continuous (limits = c(0,8))
p8

# 设定各类拟合模型，这里主要讲两类：线性模型与非线性模型
models <- list(lm(y ~ x, data = mpg_new), # lm 表示线性模型
   lm(y ~ I(1/x), data = mpg_new),
   lm(y ~ log(x), data = mpg_new),
   nls(y ~ I(x^2 *a) + b*x+c, data = mpg_new, start = list (a=1, b=1, c=0)), #nls 表示非线性模型
   nls(y ~ I(1/x * a) + b * x, data = mpg_new, start = list(a = 1, b = 1)), 
   nls(y ~ (a + b*log(x)), data = mpg_new, start = setNames(coef(lm(y ~ log(x), data=mpg_new)), c("a", "b"))),
   nls(y ~ I(exp(1)^(a + b * x)), data = mpg_new, start = list(a=0,b=1)),
   nls(y ~ I(1/x * a)+b, data = mpg_new, start = list(a=1,b=1)),
   nls(y ~ a*I(x^b), data = mpg_new, start = list(a=2,b=1.5))
) # I() 函数可以避免这种混淆将模型公式中使用运算符的部分括起来
  # 例如，在公式y ~ a + I(b+c)中，b+c项被解释为b和c的和
models[[1]] # 查看第一个公式的信息， [[X]] x为1-9共9个公式，可以直接查看对应的公式信息
# models 查看所有模型的参数

# 可以先查看不同模型拟合曲线的大致情况
p8_1 <- ggplot(mpg_new, aes(x, y)) + geom_point(size = 3.5) +
  stat_smooth(method = "lm", formula = as.formula(models[[1]]), size = 1, se = FALSE, colour = "black") +
  stat_smooth(method = "lm", formula = as.formula(models[[2]]), size = 1, se = FALSE, colour = "blue") +
  stat_smooth(method = "lm", formula = as.formula(models[[3]]), size = 1, se = FALSE, colour = "yellow") +
  stat_smooth(method = "nls", formula = as.formula(models[[5]]), size = 1, se = FALSE, colour = "purple",linetype = 2) +
  stat_smooth(method = "nls", formula = as.formula(models[[5]]), size = 1, se = FALSE, colour = "red",linetype = 2) +
  stat_smooth(method = "nls", formula = as.formula(models[[6]]), size = 1, se = FALSE, colour = "green",linetype = 2) +
  stat_smooth(method = "nls", formula = as.formula(models[[7]]), size = 1, se = FALSE, colour = "violet") +
  stat_smooth(method = "nls", formula = as.formula(models[[8]]), size = 1, se = FALSE, colour = "orange",linetype = 2,level = 0.95) +
  stat_smooth(method = "nls", formula = as.formula(models[[9]]), size = 1, se = FALSE, colour = "grey",linetype = 2,level = 0.95)
p8_1

# 通过计算AIC值，选出最优模型
# AIC值越小，拟合模型更好
library(AICcmodavg); library(plyr); library(stringr) #读取所需要用的R包
ldply(models, function(mod){ data.frame(AICc = AICc(mod), AIC = AIC(mod), model = deparse(formula(mod))) })

#8) 根据AIC值，我们选出了最优模型为 models[[4]]：y ~ a*I(x^2) + b*x+c
# 接下来我们定义函数，计算得出后面需要要添加的公式
# 先定义 lm 线性函数
lm_eqn = function(mpg_new){ 
  m <- lm(y~ I(x^2) + x, mpg_new)
  eq <- substitute(italic(y) == ~c + ~b %.% italic(x) + ~a %.% italic(x)^2*","~~italic(r)^2~"="~r2,
 list(a = format(coef(m)[2], digits = 3), # ~ 表示一个空格
   b = format(coef(m)[3], digits = 3), # ","表示，；"=" 表示 =
   c = format(coef(m)[1], digits = 3), # digits 表示保留三位有效数字
   r2 = format(summary(m)$r.squared, digits = 3))) # substitute() 函数能够解释模型的表达式
  as.character(as.expression(eq))
} # 这里需要注意的是，我发现如果按a b c (即123)的顺序出来的公式与上面模型计算出来的参数有差异
# 因此我建议以模型计算出来的参数为准，防止定义时，系统计算出错

#########################
# 非线性模型函数的定义方法

# 非线性模型R^2计算 参考http://blog.sina.com.cn/s/blog_78c5f0530101c7w1.html
# 我已经把这部分内容用# 符号隐藏了，删除 # 尝试
# lm_eqn = function(mpg_new){
# mod <- nls (y~ a*I(x^b),start = list(a=2, b= 1.5), data = mpg_new)
# a=0
# b=0
# for(i in 1:100){ # 按理这个循环要循环234次，但是100次的时候就已经是负值了
# ymean=mean(y) # 结合上面计算的AIC值，意味着这个模型不可用
# a=a+(y[i]-fitted(mod)[i])^2
# b=b+(y[i]-ymean)^2
# }
# r2=1-a/b # 由于非线性的R^2 不能够直接用summary访问， 只能用循环公式来计算
# eq <- substitute(italic(y) == ~a %.% italic(x) ^ b*","~~italic(r)^2~"="~r2,
# list(a = format(coef(mod)[1], digits = 3),
# b = format(coef(mod)[2], digits = 3),
# r2 = format(r2,digits = 3)))
# as.character(as.expression(eq))
# }

#########################
# 然后把上面定义的函数添加到图里
p9_2 <- ggplot (mpg_new, aes(x, y)) +
  geom_point (shape=19,alpha = 0.3,stroke=0,size=mpg_new$y)+ # 这里选择了16号形状 可以用show_point_shapes() 访问形状类别
  stat_smooth (method = "lm",
  formula = y~ I(x^2) + x,
  se = T,
  colour = "#00A5FF",
  fill = "#00A5FF",
  alpha = 0.2) +
  scale_x_continuous (limits = c(10,50)) +
  scale_y_continuous (limits = c(0,8)) +
  theme_bw() +
  theme(panel.border = element_blank(), 
  panel.grid=element_blank(),
  axis.line = element_line(colour = "black")) +
  theme (axis.text=element_text(colour = 'black',size = 12), 
   axis.title=element_text(size = 14,colour = 'black')) + # ' ' 与 "" 用法相同
  annotate ("text",
   x= 10, # 表示公式添加的位置 x轴
   y= 7.5, # 表示公式添加的位置 y轴
   label =lm_eqn(mpg_new), # 添加的公式为上面定义的内容 (包括计算)
   hjust=0, 
   size=6,
   parse=TRUE)+ # parse = T 表示能够解析出模型公式，如果为F 则识别的公式与正常的表达不同
  labs (x = "Soil depth (cm)", y = "Soil organic carbon (g/kg)")
p9_2
graph2ppt(file="p9_2.ppt", append=FALSE)

#############
# 举一反三
#############
#这里例举另外两种计算最优模型的方法
#9) 这里用 basicTrendline包
library("basicTrendline")
trendline(mpg_new$x, mpg_new$y, model="line2P", ePos.x = "topleft", summary=TRUE, eDigit=5)
trendline(mpg_new$x, mpg_new$y, model="line3P", CI.fill = FALSE, CI.color = "black", CI.lty = 2)
trendline(mpg_new$x, mpg_new$y, model="log2P", ePos.x= "top", linecolor = "red", CI.color = NA) 
trendline(mpg_new$x, mpg_new$y, model="exp2P", CI.fill = FALSE, CI.color = "black", CI.lty = 2)
trendline(mpg_new$x, mpg_new$y, model="exp3P", xname="T", yname = paste(delta^15,"N") ,
 yhat=FALSE, Rname=1, Pname=0, ePos.x = "bottom")
# 如果出现报错，很可能是因为这个数据用这个模型时不符合条件

#10) 这里用 ggpmisc包
# 先定义一个拟合模型
my_formular <- y ~ I(x^2) + x 

# 添加R^2
ggplot(mpg_new, aes(x, y)) +
  geom_point() +
  geom_smooth(method = "lm", formula = my_formular) +
  stat_poly_eq(formula = my_formular, rr.digits = 4, parse = TRUE)

# 添加各种标签 包括拟合模型的公式，R^2，p值，AIC值等
ggplot(mpg_new, aes(x, y)) +
  geom_point() +
  geom_smooth(method = "lm", formula = my_formular) +
  stat_poly_eq(aes(label = paste(stat(eq.label), # 可以? stat_poly_eq 里面有对应的解释
 stat(adj.rr.label), # 可以根据自己需要添加对应的信息
 stat(p.value.label),
 stat(AIC.label),
 sep = "*\", \"*")), # 这行代码是识别和分割x，如果删除则无法正常识别
 coef.digits = 4, #保留4位有效数字
   formula = my_formular, parse = TRUE)
# 但是能够发现，系统计算出来的拟合公式是错误的
# 正确的应该是 y ~ 0.003*x^2 -0.335*x +9.310 

# 更改以下代码，添加文本信息，添加上面我们定义的函数
ggplot(mpg_new, aes(x, y)) +
  geom_point() +
  geom_smooth(method = "lm", formula = my_formular) +
  stat_poly_eq(geom = "text", label.x = 20, label.y = 0, hjust = 0,
   formula = my_formular, parse = TRUE, label = lm_eqn(mpg_new))

Tips：

1) ######### 代码中的多#号在rstudio里可调节成隐藏或者展开的形式。

2) 选择最优模型时就以AIC值为主，并且在添加公式时候，参数的选择以代码计算出来的为主，系统默认计算出的拟合回归模型的系数记得比对。

3) 其中三种方法可以计算拟合回归模型的参数，个人人为最后一种可能简单，利用 ggpimsc 包里 stat_poly_eq () 函数直接把模型公式拟合出来，然后对比第二种方法计算的参数 (basicTrendline 包)，然后再PPT或者AI里稍微编辑修改即可。

遇到问题欢迎留言~

参考材料，感谢~ 

1. http://blog.sciencenet.cn/home.php?mod=space&uid=651374&do=blog&id=1126673

2. http://www./article/p-fhpazgda-btm.html

3.http://wap.sciencenet.cn/blog-587128-1212579.html?mobile=1