3

数据：

df <- structure(list(x = c(9.5638945103927, 13.7767187698566, 6.0019477258207, 
10.1897072092089, 15.4019854273531, 10.9746646056535, 12.9429073949468, 
20.7513493525379, 18.5764146937149, 2.91302077116471, 13.6523222711501, 
10.0920467755108), y = c(83.949498880077, 18.066881289085, 71.3052196358606, 
39.8975644317452, 57.2933166677927, 87.8484256883889, 92.6818329896141, 
49.8297961197214, 56.3650103496898, 14.7950650020996, 37.9271392096266, 
50.4357237591891), z = c("a", "c", "e", "f", "b", "a", "b", "a", 
"b", "a", "c", "d")), .Names = c("x", "y", "z"), row.names = c(NA, 
-12L), class = "data.frame")

我的模型：

mod <- glm(y ~ x + I(x^2) + z, family=quasipoisson, data = df)
summary(mod)

我想绘制这样的东西：

ggplot(df, aes(x=x,y=y)) + 
  geom_point() +
  stat_smooth(method="lm",se=FALSE,
              formula=y~x+I(x^2),fill="transparent",
              colour="black") +
  stat_smooth(method="lm",geom="ribbon",
              formula=y~x+I(x^2),fill="transparent",
              colour="red",linetype="dashed",fullrange=TRUE) +
          scale_x_continuous(limits=c(-2,35)) +
          coord_cartesian(xlim=c(2,25),
                          ylim=range(pretty(df$y)))

在此处输入图像描述

但是，我绘制的模型显然与 不同的模型mod，没有z也没有quasiposson。

如何绘制类似 ggplot 但使用我的实际模型的东西？然而，我看过，predict当有多个解释时，我完全不知道该怎么做，就像我的情况一样。我不在乎这样做ggplot2

4

2 回答 2

7

似乎您可以通过使用 轻松地将您的示例调整到您的新模型stat_smooth(method='glm', family=quasipoisson, ...)，但是将 添加z到公式中会导致错误。从ggplot2 docs看，您可以看到predictdf用于生成间隔限制的内容。查看该函数的代码，看起来它仅适用于跨 x 维度的预测。但是我们可以很容易地编写我们自己的在多个维度上工作的版本，然后将限制绘制为单独的层。

mypredictdf <- function (model, newdata, level=0.95){
  pred <- stats::predict(model, newdata = newdata, se =TRUE, type = "link")
  std <- qnorm(level/2 + 0.5)
  data.frame(newdata,
             y = model$family$linkinv(as.vector(pred$fit)),
             ymin = model$family$linkinv(as.vector(pred$fit - std * pred$se)),
             ymax = model$family$linkinv(as.vector(pred$fit + std * pred$se)), 
             se = as.vector(pred$se))
}
px <- with(df, seq(from=min(x), to=max(x), length=100))
pdf <- expand.grid(x=px, z=unique(df$z))
pdf <- mypredictdf(mod, newdata=pdf)
g <- ggplot(data=pdf, aes(group=z))
g <- g + geom_point(data=df, aes(x=x, y=y, color=z))
g <- g + geom_ribbon(aes(x=x, ymin=ymin, ymax=ymax),
                     alpha=0.2)
g <- g + geom_line(aes(x=x, y=y, color=z))

一个面板

看起来刻面是个好主意：

g <- g + facet_wrap(~z)

刻面版

于 2014-03-15T00:11:58.570 回答
3

Here is one approach that deals with multiple variables ( y = f(x,z) in your case).

mod <- glm(y ~ x + I(x^2) + z, family=quasipoisson, data = df)
pred <- predict(mod, type="response",se.fit=T)
df$pred <- pred$fit
df$se   <- pred$se.fit

ggplot(df, aes(x=y))+
  geom_point(aes(y=pred, color=z),size=3)+
  geom_errorbar(aes(ymin=pred-se, ymax=pred+se, color=z),width=1.5)+
  geom_abline(intercept=0, slope=1, color="blue", linetype=2)+
  labs(x="Actual", y="Predicted")

This plots predicted y vs. actual y, grouped by z, with error bars = ±1 × se. To get 95% CL on the predicted, you would need to use ±1.96 × se. The dotted line is a reference (actual = predicted) which would represent a perfect fit. You can see from this that z=b and z=c are problematic, but z in (a,d,e,f) all fit the data rather well.

If you have more than 2 variables, the grouping becomes problematic, but the asic approach of plotting predicted y vs. actual y still applies.

于 2014-03-06T23:05:24.997 回答