使用 ggplot 从其系数绘制回归
Plotting regression from its coefficients with ggplot
我正在尝试用 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
或类似的,包括在每个级别中相应的估计值,但无法弄清楚如何。
我通常通过生成预测数据帧来做到这一点。
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.
我认为查看 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()
我正在尝试用 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
或类似的,包括在每个级别中相应的估计值,但无法弄清楚如何。
我通常通过生成预测数据帧来做到这一点。
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.
我认为查看 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()