【问题标题】:Bad plot when plotting 95% confidence interval of linear model prediction绘制线性模型预测的 95% 置信区间时的错误图
【发布时间】:2017-02-12 09:27:00
【问题描述】:

我需要一些帮助来绘制带有置信区间的预测。考虑下面的例子

library(Hmisc)
data("mtcars") 

mfit = lm(mpg ~ vs + disp + cyl, data = mtcars)

#disp and cyl at their mean
newcar = data.frame(vs = c(0,1), disp = 230, cyl = 6.188)

pmodel <- predict(mfit, newcar, se.fit=TRUE)

当所有其他变量保持不变(平均值/众数)时,我想绘制vs(0 和 1 时)的效果。

为此,我在下面运行此代码:

plot(1:2, pmodel$fit[1:2], ylim=c(0,1), pch=19, xlim=c(.5,2.5), xlab="X",
     ylab = "Predicted values", xaxt = "n", main = "Figure1")
arrows(1:2, (pmodel$fit[1:2] - 1.96 * pmodel$fit[1:2]),
       1:2, (pmodel$fit[1,1] + 1.96 * pmodel$fit[1:2]),
       length=0.05, angle=90, code=3)
axis(1, at=c(1,2), labels=c("0","1"))

我在这里做错了什么?谢谢!

【问题讨论】:

  • 您也可以使用专门设计的函数,例如 visreg 包:library(visreg); visreg(mfit, xvar="vs")

标签: r plot regression linear-regression lm


【解决方案1】:

注意你有ylim = c(0, 1),这是不正确的。在绘制置信区间时,我们必须确保ylim 覆盖CI 的下限和上限。

## lower and upper bound of CI
lower <- with(pmodel, fit - 1.96 * se.fit)
upper <- with(pmodel, fit + 1.96 * se.fit)

## x-location to plot
xx <- 0:1

## set `xlim` and `ylim`
xlim <- range(xx) + c(-0.5, 0.5)  ## extends an addition 0.5 on both sides
ylim <- range(c(lower, upper))

## produce figure
plot(xx, pmodel$fit, pch = 19, xlim = xlim, ylim = ylim, xaxt = "n",
     xlab = "X", ylab = "Predicted values", main = "Figure1")
arrows(xx, lower, xx, upper, length = 0.05, angle = 90, code = 3)
axis(1, at = xx)

您的代码中还有其他几个 cmets:

  • fit - 1.96 * se.fit 不是fit - 1.96 * fit;
  • 您可以直接在 x-location 0:1 上绘图,而不是 1:2
  • fit[1:2] 而不是fit[1,1]

【讨论】:

    【解决方案2】:

    在ggplot中:

    df <- data.frame(x=1:2, pred=pmodel$fit[1:2])
    df$lower <- df$pred - 1.96 * pmodel$se.fit[1:2] 
    df$upper <- df$pred + 1.96 * pmodel$se.fit[1:2]
    ggplot(df, aes(x, pred, group=x, col=as.factor(x))) + geom_point(size=2) + 
      geom_errorbar(aes(ymin=lower, ymax=upper), width=0.1)
    

    【讨论】:

      猜你喜欢
      • 2022-01-13
      • 2021-07-05
      • 2017-11-20
      • 2020-09-16
      • 1970-01-01
      • 1970-01-01
      • 2018-03-09
      • 2013-04-21
      • 2021-07-06
      相关资源
      最近更新 更多