在图上添加回归线方程和 R^2
Posted
技术标签:
【中文标题】在图上添加回归线方程和 R^2【英文标题】:Add regression line equation and R^2 on graph 【发布时间】:2011-11-24 20:50:47 【问题描述】:我想知道如何在ggplot
上添加回归线方程和 R^2。我的代码是:
library(ggplot2)
df <- data.frame(x = c(1:100))
df$y <- 2 + 3 * df$x + rnorm(100, sd = 40)
p <- ggplot(data = df, aes(x = x, y = y)) +
geom_smooth(method = "lm", se=FALSE, color="black", formula = y ~ x) +
geom_point()
p
我们将不胜感激。
【问题讨论】:
对于 lattice 图形,请参阅latticeExtra::lmlineq()
。
@JoshO'Brien Error: 'lmlineq' is not an exported object from 'namespace:latticeExtra'
【参考方案1】:
另一种选择是创建一个自定义函数,使用 dplyr
和 broom
库生成方程:
get_formula <- function(model)
broom::tidy(model)[, 1:2] %>%
mutate(sign = ifelse(sign(estimate) == 1, ' + ', ' - ')) %>% #coeff signs
mutate_if(is.numeric, ~ abs(round(., 2))) %>% #for improving formatting
mutate(a = ifelse(term == '(Intercept)', paste0('y ~ ', estimate), paste0(sign, estimate, ' * ', term))) %>%
summarise(formula = paste(a, collapse = '')) %>%
as.character
lm(y ~ x, data = df) -> model
get_formula(model)
#"y ~ 6.22 + 3.16 * x"
scales::percent(summary(model)$r.squared, accuracy = 0.01) -> r_squared
现在我们需要将文本添加到绘图中:
p +
geom_text(x = 20, y = 300,
label = get_formula(model),
color = 'red') +
geom_text(x = 20, y = 285,
label = r_squared,
color = 'blue')
【讨论】:
【参考方案2】:这里是大家最简单的代码
注意:显示 Pearson 的 Rho 和 not R^2。
library(ggplot2)
library(ggpubr)
df <- data.frame(x = c(1:100)
df$y <- 2 + 3 * df$x + rnorm(100, sd = 40)
p <- ggplot(data = df, aes(x = x, y = y)) +
geom_smooth(method = "lm", se=FALSE, color="black", formula = y ~ x) +
geom_point()+
stat_cor(label.y = 35)+ #this means at 35th unit in the y axis, the r squared and p value will be shown
stat_regline_equation(label.y = 30) #this means at 30th unit regresion line equation will be shown
p
【讨论】:
与上述相同的问题,在您的情节中显示为 rho 而不是 R²! 实际上你可以只添加 R2:stat_cor(aes(label = ..rr.label..))
我发现这是最简单的解决方案,可以最好地控制标签的位置(我找不到使用 stat_poly_eq 将 R^2 置于方程下方的简单方法)并且可以结合stat_regline_equation()
绘制回归方程
'ggpubr' 似乎没有被积极维护;因为它在 GitHub 中有许多未解决的问题。无论如何,stat_regline_equation()
和stat_cor()
中的大部分代码都是在我的包“ggpmisc”中未经确认而复制的。它取自stat_poly_eq()
,它得到了积极维护,并在复制后获得了一些新功能。示例代码需要最少的编辑才能使用“ggpmisc”。【参考方案3】:
使用ggpubr:
library(ggpubr)
# reproducible data
set.seed(1)
df <- data.frame(x = c(1:100))
df$y <- 2 + 3 * df$x + rnorm(100, sd = 40)
# By default showing Pearson R
ggscatter(df, x = "x", y = "y", add = "reg.line") +
stat_cor(label.y = 300) +
stat_regline_equation(label.y = 280)
# Use R2 instead of R
ggscatter(df, x = "x", y = "y", add = "reg.line") +
stat_cor(label.y = 300,
aes(label = paste(..rr.label.., ..p.label.., sep = "~`,`~"))) +
stat_regline_equation(label.y = 280)
## compare R2 with accepted answer
# m <- lm(y ~ x, df)
# round(summary(m)$r.squared, 2)
# [1] 0.85
【讨论】:
您是否见过一种简洁的编程方式来指定label.y
的数字?
@MarkNeal 可能会得到 y 的最大值,然后乘以 0.8。 label.y = max(df$y) * 0.8
@MarkNeal 好点子,也许可以在 GitHub ggpubr 上将问题作为功能请求提交。
提交的自动定位问题here
@zx8754 ,在你的情节中,它显示为 rho 而不是 R² ,有什么简单的方法可以显示 R² 吗?【参考方案4】:
我在我的包ggpmisc
中包含了一个统计信息stat_poly_eq()
,它允许这个答案:
library(ggplot2)
library(ggpmisc)
df <- data.frame(x = c(1:100))
df$y <- 2 + 3 * df$x + rnorm(100, sd = 40)
my.formula <- y ~ x
p <- ggplot(data = df, aes(x = x, y = y)) +
geom_smooth(method = "lm", se=FALSE, color="black", formula = my.formula) +
stat_poly_eq(formula = my.formula,
aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")),
parse = TRUE) +
geom_point()
p
此统计数据适用于没有缺失项的任何多项式,并且希望具有足够的灵活性以普遍有用。 R^2 或调整后的 R^2 标签可以与任何带有 lm() 的模型公式一起使用。作为一个 ggplot 统计数据,它在组和方面的表现都符合预期。
“ggpmisc”包可通过 CRAN 获得。
版本 0.2.6 刚刚被 CRAN 接受。
它通过 @shabbychef 和 @MYaseen208 解决 cmets。
@MYaseen208 这显示了如何添加一个帽子。
library(ggplot2)
library(ggpmisc)
df <- data.frame(x = c(1:100))
df$y <- 2 + 3 * df$x + rnorm(100, sd = 40)
my.formula <- y ~ x
p <- ggplot(data = df, aes(x = x, y = y)) +
geom_smooth(method = "lm", se=FALSE, color="black", formula = my.formula) +
stat_poly_eq(formula = my.formula,
eq.with.lhs = "italic(hat(y))~`=`~",
aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")),
parse = TRUE) +
geom_point()
p
@shabbychef 现在可以将方程中的变量与用于轴标签的变量进行匹配。要将 x 替换为 z 并将 y 替换为 h 可以使用:
p <- ggplot(data = df, aes(x = x, y = y)) +
geom_smooth(method = "lm", se=FALSE, color="black", formula = my.formula) +
stat_poly_eq(formula = my.formula,
eq.with.lhs = "italic(h)~`=`~",
eq.x.rhs = "~italic(z)",
aes(label = ..eq.label..),
parse = TRUE) +
labs(x = expression(italic(z)), y = expression(italic(h))) +
geom_point()
p
作为这些正常的 R 解析表达式,希腊字母现在也可以在等式的 lhs 和 rhs 中使用。
[2017-03-08] @elarry 编辑以更准确地解决原始问题,展示如何在方程式标签和 R2 标签之间添加逗号。
p <- ggplot(data = df, aes(x = x, y = y)) +
geom_smooth(method = "lm", se=FALSE, color="black", formula = my.formula) +
stat_poly_eq(formula = my.formula,
eq.with.lhs = "italic(hat(y))~`=`~",
aes(label = paste(..eq.label.., ..rr.label.., sep = "*plain(\",\")~")),
parse = TRUE) +
geom_point()
p
[2019-10-20] @helen.h 我在下面给出了 stat_poly_eq()
与分组一起使用的示例。
library(ggpmisc)
df <- data.frame(x = c(1:100))
df$y <- 20 * c(0, 1) + 3 * df$x + rnorm(100, sd = 40)
df$group <- factor(rep(c("A", "B"), 50))
my.formula <- y ~ x
p <- ggplot(data = df, aes(x = x, y = y, colour = group)) +
geom_smooth(method = "lm", se=FALSE, formula = my.formula) +
stat_poly_eq(formula = my.formula,
aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")),
parse = TRUE) +
geom_point()
p
p <- ggplot(data = df, aes(x = x, y = y, linetype = group)) +
geom_smooth(method = "lm", se=FALSE, formula = my.formula) +
stat_poly_eq(formula = my.formula,
aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")),
parse = TRUE) +
geom_point()
p
[2020-01-21] @Herman 乍一看可能有点反直觉,但是在使用分组时要获得单个方程需要遵循图形的语法。要么将创建分组的映射限制为单个图层(如下所示),要么保留默认映射并在您不希望分组的图层中使用常量值覆盖它(例如colour = "black"
)。
继续前面的例子。
p <- ggplot(data = df, aes(x = x, y = y)) +
geom_smooth(method = "lm", se=FALSE, formula = my.formula) +
stat_poly_eq(formula = my.formula,
aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")),
parse = TRUE) +
geom_point(aes(colour = group))
p
[2020-01-22] 为了完整起见,一个带有分面的示例,证明在这种情况下也满足了图形语法的期望。
library(ggpmisc)
df <- data.frame(x = c(1:100))
df$y <- 20 * c(0, 1) + 3 * df$x + rnorm(100, sd = 40)
df$group <- factor(rep(c("A", "B"), 50))
my.formula <- y ~ x
p <- ggplot(data = df, aes(x = x, y = y)) +
geom_smooth(method = "lm", se=FALSE, formula = my.formula) +
stat_poly_eq(formula = my.formula,
aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")),
parse = TRUE) +
geom_point() +
facet_wrap(~group)
p
【讨论】:
需要注意的是,公式中的x
和y
是指绘图层中的x
和y
数据,不一定是在范围内my.formula
的构建时间。因此,公式应该总是使用 x 和 y 变量吗?
好点@elarry!这与 R 的 parse() 函数的工作方式有关。通过反复试验,我发现aes(label = paste(..eq.label.., ..rr.label.., sep = "*plain(\",\")~"))
可以完成这项工作。
@HermanToothrot 通常 R2 是回归的首选,因此stat_poly_eq()
返回的数据中没有预定义的 r.label。你可以使用stat_fit_glance()
,同样来自包'ggpmisc',它以数值形式返回R2。请参阅帮助页面中的示例,并将stat(r.squared)
替换为sqrt(stat(r.squared))
。
@PedroAphalo 如果我使用像公式 = y~x+z 这样的多元模型,是否可以重命名第三个变量?
我才知道,显然,我们不能在 plotly 中使用 ggpmisc::stat_poly_eq,它没有在 plotly 中实现。【参考方案5】:
这是一种解决方案
# GET EQUATION AND R-SQUARED AS STRING
# SOURCE: https://groups.google.com/forum/#!topic/ggplot2/1TgH-kG5XMA
lm_eqn <- function(df)
m <- lm(y ~ x, df);
eq <- substitute(italic(y) == a + b %.% italic(x)*","~~italic(r)^2~"="~r2,
list(a = format(unname(coef(m)[1]), digits = 2),
b = format(unname(coef(m)[2]), digits = 2),
r2 = format(summary(m)$r.squared, digits = 3)))
as.character(as.expression(eq));
p1 <- p + geom_text(x = 25, y = 300, label = lm_eqn(df), parse = TRUE)
编辑。我找出了我选择此代码的来源。这是 ggplot2 google 群组中原始帖子的link
【讨论】:
@JonasRaedle 关于使用annotate
获得更好看的文本的评论在我的机器上是正确的。
这看起来不像我机器上发布的输出,其中标签被覆盖的次数与调用数据的次数一样多,导致标签文本厚而模糊。首先将标签传递给 data.frame (请参阅下面的评论中的我的建议。
@PatrickT:删除aes(
和对应的)
。 aes
用于将数据帧变量映射到可视变量——这里不需要,因为只有一个实例,所以你可以把它全部放在主要的 geom_text
调用中。我会将其编辑为答案。
这个解决方案的问题似乎是,如果数据集更大(我的数据集是 370000 个观察值),该函数似乎会失败。我会推荐来自@kdauria 的解决方案,它的作用相同,但要快得多。
对于那些想要 r 和 p 值而不是 R2 和方程的人:eq
【参考方案6】:
受this answer 中提供的方程式风格的启发,一种更通用的方法(多个预测器 + 乳胶输出作为选项)可以是:
print_equation= function(model, latex= FALSE, ...)
dots <- list(...)
cc= model$coefficients
var_sign= as.character(sign(cc[-1]))%>%gsub("1","",.)%>%gsub("-"," - ",.)
var_sign[var_sign==""]= ' + '
f_args_abs= f_args= dots
f_args$x= cc
f_args_abs$x= abs(cc)
cc_= do.call(format, args= f_args)
cc_abs= do.call(format, args= f_args_abs)
pred_vars=
cc_abs%>%
paste(., x_vars, sep= star)%>%
paste(var_sign,.)%>%paste(., collapse= "")
if(latex)
star= " \\cdot "
y_var= strsplit(as.character(model$call$formula), "~")[[2]]%>%
paste0("\\hat",.,"_i")
x_vars= names(cc_)[-1]%>%paste0(.,"_i")
else
star= " * "
y_var= strsplit(as.character(model$call$formula), "~")[[2]]
x_vars= names(cc_)[-1]
equ= paste(y_var,"=",cc_[1],pred_vars)
if(latex)
equ= paste0(equ," + \\hat\\varepsilon_i \\quad where \\quad \\varepsilon \\sim \\mathcalN(0,",
summary(MetamodelKdifEryth)$sigma,")")%>%paste0("$",.,"$")
cat(equ)
model
参数需要一个 lm
对象,latex
参数是一个布尔值,用于请求简单字符或乳胶格式的方程,...
参数将其值传递给 format
功能。
我还添加了一个将其输出为乳胶的选项,因此您可以在 rmarkdown 中使用此功能,如下所示:
```r echo=FALSE, results='asis'
print_equation(model = lm_mod, latex = TRUE)
```
现在使用它:
df <- data.frame(x = c(1:100))
df$y <- 2 + 3 * df$x + rnorm(100, sd = 40)
df$z <- 8 + 3 * df$x + rnorm(100, sd = 40)
lm_mod= lm(y~x+z, data = df)
print_equation(model = lm_mod, latex = FALSE)
此代码产生:
y = 11.3382963933174 + 2.5893419 * x + 0.1002227 * z
如果我们要求乳胶方程,将参数四舍五入为 3 位:
print_equation(model = lm_mod, latex = TRUE, digits= 3)
这会产生:
【讨论】:
【参考方案7】:真的很喜欢@Ramnath 解决方案。为了允许使用自定义回归公式(而不是固定为 y 和 x 作为文字变量名称),并将 p 值添加到打印输出中(正如@Jerry T 评论的那样),这里是 mod:
lm_eqn <- function(df, y, x)
formula = as.formula(sprintf('%s ~ %s', y, x))
m <- lm(formula, data=df);
# formating the values into a summary string to print out
# ~ give some space, but equal size and comma need to be quoted
eq <- substitute(italic(target) == a + b %.% italic(input)*","~~italic(r)^2~"="~r2*","~~p~"="~italic(pvalue),
list(target = y,
input = x,
a = format(as.vector(coef(m)[1]), digits = 2),
b = format(as.vector(coef(m)[2]), digits = 2),
r2 = format(summary(m)$r.squared, digits = 3),
# getting the pvalue is painful
pvalue = format(summary(m)$coefficients[2,'Pr(>|t|)'], digits=1)
)
)
as.character(as.expression(eq));
geom_point() +
ggrepel::geom_text_repel(label=rownames(mtcars)) +
geom_text(x=3,y=300,label=lm_eqn(mtcars, 'hp','wt'),color='red',parse=T) +
geom_smooth(method='lm')
不幸的是,这不适用于 facet_wrap 或 facet_grid。
【讨论】:
非常简洁,我引用了here。澄清一下 - 您的代码在 geom_point() 之前是否缺少ggplot(mtcars, aes(x = wt, y = mpg, group=cyl))+
?一个半相关的问题 - 如果我们在 aes()
中为 ggplot 引用 hp 和 wt,那么我们可以 grab 将它们用于调用lm_eqn
,那么我们只需要在一个地方编码?我知道我们可以在调用 ggplot() 之前设置 xvar = "hp"
,并在两个位置都使用 xvar 来替换 hp,但是这感觉应该是不必要的。
非常好的解决方案!感谢分享!【参考方案8】:
我更改了stat_smooth
和相关函数的源代码的几行,以创建一个添加拟合方程和R 平方值的新函数。这也适用于构面图!
library(devtools)
source_gist("524eade46135f6348140")
df = data.frame(x = c(1:100))
df$y = 2 + 5 * df$x + rnorm(100, sd = 40)
df$class = rep(1:2,50)
ggplot(data = df, aes(x = x, y = y, label=y)) +
stat_smooth_func(geom="text",method="lm",hjust=0,parse=TRUE) +
geom_smooth(method="lm",se=FALSE) +
geom_point() + facet_wrap(~class)
我使用@Ramnath 答案中的代码来格式化等式。 stat_smooth_func
函数不是很健壮,但使用它应该不难。
https://gist.github.com/kdauria/524eade46135f6348140。如果出现错误,请尝试更新 ggplot2
。
【讨论】:
非常感谢。这不仅适用于方面,甚至适用于组。我发现它对于分段回归非常有用,例如stat_smooth_func(mapping=aes(group=cut(x.val,c(-70,-20,0,20,50,130))),geom="text",method="lm",hjust=0,parse=TRUE)
,结合 ***.com/questions/19735149/… 的 EvaluateSmooths
@aelwan,随意更改这些行:gist.github.com/kdauria/…。然后source
脚本中的整个文件。
@kdauria 如果我在每个 facet_wraps 中有几个方程并且我在每个 facet_wrap 中有不同的 y_values 怎么办。有什么建议如何固定方程的位置?我使用此示例dropbox.com/s/9lk9lug2nwgno2l/R2_facet_wrap.docx?dl=0 尝试了 hjust、vjust 和 angle 的几个选项,但我无法在每个 facet_wrap 中将所有方程置于同一级别
@aelwan,等式的位置由以下几行确定:gist.github.com/kdauria/…。我在 Gist 中制作了函数的 xpos
和 ypos
参数。因此,如果您希望所有方程重叠,只需设置xpos
和ypos
。否则,xpos
和 ypos
根据数据计算。如果你想要更高级的东西,在函数中添加一些逻辑应该不会太难。例如,也许您可以编写一个函数来确定图形的哪一部分空间最多,然后将函数放在那里。
我在 source_gist 中遇到了一个错误:r_files[[which]] 中的错误:无效的下标类型“闭包”。请参阅此帖子以获取解决方案:***.com/questions/38345894/r-source-gist-not-working【参考方案9】:
我已将 Ramnath 的帖子修改为 a) 使其更通用,因此它接受线性模型作为参数而不是数据框,并且 b) 更恰当地显示负数。
lm_eqn = function(m)
l <- list(a = format(coef(m)[1], digits = 2),
b = format(abs(coef(m)[2]), digits = 2),
r2 = format(summary(m)$r.squared, digits = 3));
if (coef(m)[2] >= 0)
eq <- substitute(italic(y) == a + b %.% italic(x)*","~~italic(r)^2~"="~r2,l)
else
eq <- substitute(italic(y) == a - b %.% italic(x)*","~~italic(r)^2~"="~r2,l)
as.character(as.expression(eq));
用法将更改为:
p1 = p + geom_text(aes(x = 25, y = 300, label = lm_eqn(lm(y ~ x, df))), parse = TRUE)
【讨论】:
这看起来很棒!但我在多个方面绘制 geom_points,其中 df 根据 facet 变量而有所不同。我该怎么做? Jayden 的解决方案效果很好,但是字体看起来很丑。我建议将用法更改为:p1 = p + annotate("text", x = 25, y = 300, label = lm_eqn(lm(y ~ x, df)), colour="black", size = 5, parse=TRUE)
编辑:这也解决了您在图例中显示字母时可能遇到的任何问题。
@Jonas,出于某种原因,我收到了"cannot coerce class "lm" to a data.frame"
。这种替代方法有效:df.labs <- data.frame(x = 25, y = 300, label = lm_eqn(df))
和 p <- p + geom_text(data = df.labs, aes(x = x, y = y, label = label), parse = TRUE)
@PatrickT - 如果您使用 Ramnath 的解决方案调用 lm_eqn(lm(...))
,就会收到错误消息。您可能在尝试了一个之后又尝试了这个,但忘记确保您已重新定义 lm_eqn
@PatrickT:您能否将您的答案单独作为一个答案?我很乐意投票!以上是关于在图上添加回归线方程和 R^2的主要内容,如果未能解决你的问题,请参考以下文章