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