使用 ggplot 从其系数绘制回归答案

【问题标题】：Plotting regression from its coefficients with ggplot使用 ggplot 从其系数绘制回归
【发布时间】：2018-02-24 13:22:09
【问题描述】：

我正在尝试使用ggplot 绘制一些线性和多项式回归。这在估计geom_smoot 函数内部的回归系数时非常简单：

ggplot (mtcars, aes(x=wt, y=mpg, fill=factor(cyl), colour=factor(cyl))) + geom_smooth(method='lm', formula = y ~ poly(x,2)) + geom_point()

但是，在这里我只想根据之前关于回归参数的知识绘制一个预测（或更多，如上例所示）。

所以这里我的回归估计可以通过以下方式打印：

dlply(mtcars,.(cyl), lm, formula=mpg ~ poly(wt,2)) %>%
llply(summary) %>%
ldply(coefficients)

现在我想以相反的方式构建绘图，从估计到绘图。或者更好的是，根据这些估计值的其他值构建预测（例如，Intercept=20、poly(wt,2)1=-15 和 poly(wt,2)2=4 用于cyl=4），然后获得如上图所示的图。

但这里是我不知道如何进行的地方。我想我需要为cyl 的每个级别使用不同的geom_smooth、geom_line 或类似的，包括在每个级别中的相应估计值，但无法弄清楚如何。

【问题讨论】：

标签： r ggplot2 regression

【解决方案1】：

我通常通过生成预测数据框来做到这一点。

mod <- lm(mpg ~ poly(wt,2), mtcars)
pred <- data.frame(wt = seq(0,6,0.01))
pred$mpg <- predict(mod, pred)
ggplot() +
  geom_line(data = pred, aes(x=wt, y=mpg)) +
  geom_point(data = mtcars, aes(x=wt, y=mpg, colour=factor(cyl)))

当然，您可以将参数更改为您喜欢的任何值。

pred$mpg <- 57 - pred$wt * 21 + pred$wt^2 * 3.3

或者，您可以使用stat_function：

ggplot(pred, aes(x=wt)) +
  stat_function(fun = function(x) 57 - 21*x + 3.3*x^2) +
  geom_point(data = mtcars, aes(y=mpg, colour=factor(cyl)))

最后一点：you can't interpret the coefficients of a poly() fit the way you think you would。

【讨论】：

此预测基于lm 函数估计的回归参数。我如何使用它来根据问题中指出的自定义参数来男性化一条线？
根据您的建议，我实际上可以手动修改估算值（ej：mod$coefficients[3] = 4），然后使用prediction。没有比这更直接的方法了，不是吗？
编辑了我的答案以提供解决方案。

【解决方案2】：

我认为这可能是查看broom 包的一个很好的练习。我不太确定你想走哪条路，所以这里有一些我发现的例子：

绘制回归：

我不知道你是如何绘制多项式函数的，所以这是一个练习，但这里有一些代码可以将多项式回归到数据框中：

library(dplyr)
library(broom)
library(tidyr)

mtcars %>% group_by(cyl) %>% do(tidy(lm(mpg ~ poly(wt, 2), data=.))) %>%
     select(cyl, term, estimate) %>%
     spread(term, estimate)
# Source: local data frame [3 x 4]
# Groups: cyl [3]
# 
# cyl `(Intercept)` `poly(wt, 2)1` `poly(wt, 2)2`
# * <dbl>         <dbl>          <dbl>          <dbl>
# 1     4      26.66364     -10.170962       3.003872
# 2     6      19.74286      -2.426656      -1.589859
# 3     8      15.10000      -6.003055      -1.933630

但这里有一个线性回归：

fit <- mtcars %>% group_by(cyl) %>% do(tidy(lm(mpg ~ wt, data=.))) %>%
  select(cyl, term, estimate) %>%
  spread(term, estimate)
ggplot(mtcars, aes(x=wt, y=mpg, colour=cyl)) + geom_point() +
  geom_abline(data=fit, aes(slope=wt, intercept=`(Intercept)`, colour=cyl))

您不能只绘制拟合图，因为您需要提供 x 和 y 值，因此可能需要一些预测值：

wt <- c(2:5)
mtcars %>% group_by(cyl) %>% do(augment(lm(mpg ~ poly(wt, 2), data=.), newdata=data.frame(wt=wt))) %>%
  ggplot(aes(x=wt, y=.fitted, group=cyl, colour=cyl)) + geom_line()

【讨论】：