在R中创建多个偏差-方差折衷图

Question

我相对来说是R语言的新手。我想知道如何创建以下这个图形。我已经卡在这里已经超过两个小时了。

假设红线是真实关系 y=x^2。假设我想对100个随机样本（蓝线）进行100个线性模型的拟合。

我该怎么做呢？到目前为止，这就是我所拥有的内容：

# create the true relationship
f <- function(x) x^2  # true model
x <- seq(0, 1, by = 0.01)
y <- f(x)

# plot the true function
plot(x, y, type = "l", col = "red", ylim = c(-0.2, 1.2), lwd = 4)

# fit 100 models
set.seed(1)
for (i in 1:100)
{
    errors <- rnorm(n, 0, sigma)       # random errors, have standard deviation sigma
    obs_y <- f(obs_x) + errors         # observed y = true_model + error
    model <- lm(obs_y ~ obs_x)         # fit a linear model to the observed values
    points(obs_x[i], mean(obs_y[i]), col = "green")       # mean values
    abline(model, col = "purple")    # plot the fitted model
}

这将创建以下内容：

在这里，绿点明显关闭了...而且我没有黑点...

谢谢！

- sweetmusicality

---

最好发布你的代码。幸运的是，在你编辑之前我看了一眼。如果是平均线性模型，可以通过将“points”行移出循环来解决绿点问题。 - Sixiang.Hu

Answer 1

这是你的代码经过若干调整后的版本：

f <- function(x) x^2
x <- seq(0, 1, by = 0.05)
n <- length(x)
sigma <- 0.05
y <- f(x)

plot(x, y, type = "l", col = "red", ylim = c(-0.2, 1.2), lwd = 2)

fitted <- ys <- matrix(0, ncol = n, nrow = 100)

set.seed(1)
for (i in 1:100)
{
  errors <- rnorm(n, 0, sigma)
  ys[i, ] <- obs_y <- f(x) + errors
  model <- lm(obs_y ~ x)
  fitted[i, ] <- fitted(model)
  abline(model, col = "purple", lwd = 0.1)
}

points(x = rep(x, each = 100), y = ys, cex = 0.1)
points(x = x, y = colMeans(fitted), col = 'green', cex = 0.3)

- Julius Vainora

---

哇！谢谢你！你能解释一下你的 ys 和 fitted 矩阵代表什么吗？它们是 y-hat 值吗？ - sweetmusicality

@ sweetmusicality，确实，fitted的每一行都包含每个模型的y-hat / fitted值。然后，取该矩阵的列平均值即可得到绿色点。此外，ys的每一行都包含每种情况下观察到的y，这就是黑点所代表的。 - Julius Vainora

等等，fitted和ys有什么不同？ - sweetmusicality

@sweetmusicality，“ys”是我们观察到的（一个带有一些噪声的二次曲线），而“fitted”是我们用模型拟合出来的最佳结果（一条直线）。 - Julius Vainora

Answer 2

使用ggplot，例如：

library(ggplot2)

x <- seq(0, 1, 0.1)
y <- x^2

dat <- as.data.frame(do.call(rbind, lapply(1:100, function(i){
  y_err <- y + rnorm(1, 0, 0.06)
  l <- lm(y_err ~ x)$coefficients
  cbind(samp = i, intercept = l[1], slope = l[2], t(x * l[2] + l[1]), t(y_err))
})), row.names = 1:100)

ggplot() +  
  geom_abline(aes(intercept = dat$intercept, slope = dat$slope)) + 
  geom_point(aes(x = rep(x, each = 100), y = unlist(dat[, 15:25])), alpha = 0.5) +
  geom_line(aes(x = x, y = y), color = "red", lwd = 2) +
  geom_point(aes(x = x, y = colMeans(dat[, 4:14])), color = "green")

- erocoar