用 R 进行线性回归模拟后绘图

Question

我正在用 R 做线性回归的模拟。

我考虑的回归模型是y_i = a + b_1 * x_1i + b_2 * x_2i + e_i.

参数设计如下：

x_1i ~ N(2,1), x_2i ~ 泊松 (4), e_i ~ N(0, 1), theta = (a, b_1, b_2)

我正在做的以下代码是我想使用上面提到的分布生成 100 个 (y, x_1, x_2) 的独立随机样本 1000 次，并且我还想估计 theta_hat （theta 的估计量）。得到 theta_hat 后，我想绘制 (a_hat)、b_1 (b_1_hat)、b_2 (b_2_hat), 分别.

## Construct 1000 x_1
x_1_1000 <- as.data.frame(replicate(n = 1000,expr = rnorm(n = 100, 
  mean = 2, sd = 1)))
colnames(x_1_1000) <- paste("x_1", 1:1000, sep = "_")

x_2_1000 <- as.data.frame(replicate(n = 1000,expr = rpois(n = 100, 
  lambda = 4)))
colnames(x_2_1000) <- paste("x_2", 1:1000, sep = "_")

error_1000 <- as.data.frame(replicate(n = 1000, expr = rnorm(n = 100,
  mean = 0, sd = 1)))
colnames(error_1000) <- paste("e", 1:1000, sep = "_")

y_1000 <- as.data.frame(matrix(data = 0, nrow = 100, ncol = 1000))
y_1000 = 1 + x_1_1000 * 1 + x_2_1000*(-2) + error_1000
colnames(y_1000) <- paste("y", 1:1000, sep = "_")
######################################################################
lms <- lapply(1:1000, function(x) lm(y_1000[,x] ~ x_1_1000[,x] + x_2_1000[,x]))
theta_hat_1000 <- as.data.frame(sapply(lms, coef))

进行线性回归后，我只是将结果存储到 lms 中，这是一个列表。因为我只是想要系数的数据，所以我将这些模拟系数存储到“theta_hat_1000”但是，当我想绘制分布图时，我最终无法得到我想要的。我尝试了两种方法解决问题，但仍然很困惑。

我尝试的第一种方法是将数据框重命名为“theta_hat_1000”。我已成功重命名 column_i，其中 i 从 1 到 1000。但是，我无法成功重命名行。

rownames(theta_hat_1000[1,]) <- "ahat"
rownames(theta_hat_1000[2,]) <- "x1hat"
rownames(theta_hat_1000[3,]) <- "x2hat"

上面列出的代码没有显示任何错误消息，但最终未能更改行名称。因此，我尝试了以下代码

rownames(theta_hat_1000) <- c("ahat", "x1hat", "x2hat")

已成功重命名。但是，当我想检查数据框中是否存储了任何内容时，它会报告“NULL”

theta_hat_1000$ahat

NULL

因此，我注意到有些奇怪的地方。因此，我尝试了如下第二种方式。

我试图取消列出“theta_hat_1000”，这是一个存储在我的全局环境中的列表。但是，在做了这样的事情之后，我没有得到我想要的。预期结果只是得到三行，每行有 1000 个值，但实际情况是我得到了 3000 个 obs 和 1 列。

理想的结果是得到三列，每列有 1000 个值，并将它们放入数据框中以进行进一步的处理，例如使用 ggplot 来演示估计系数的分布。

我坚持了好几个小时。如果有人能帮助我并给我一些建议，我将不胜感激。

Answer 1

代码中的这一行 theta_hat_1000$ahat 不起作用，因为“ahat”是数据框中的行名而不是列名。您将通过调用 theta_hat_1000["ahat",].

获得结果

不过，据我所知，您想要的结果实际上是一个包含 3 列（和 1000 行）的数据框，代表回归模型的 3 个参数（截距、x1、x2）。您代码中的这一行 as.data.frame(sapply(lms, coef)) 生成一个包含 3 行和 1000 列的数据框。例如，您可以在将矩阵更改为数据框之前转置矩阵以获得 1000 行和 3 列。

theta_hat_1000 <- sapply(lms, coef)
theta_hat_1000 <- as.data.frame(t(theta_hat_1000))
colnames(theta_hat_1000) <- c("ahat", "x1hat", "x2hat")

head(theta_hat_1000)
       ahat     x1hat     x2hat
1 2.0259326 0.7417404 -2.111874
2 0.7827929 0.9437324 -1.944320
3 1.1034906 1.0091594 -2.035405
4 0.9677150 0.8168757 -1.905367
5 1.0518646 0.9616123 -1.985357
6 0.8600449 1.0781489 -2.017061

现在您还可以使用 theta_hat_1000$ahat.

调用变量

用 R 进行线性回归模拟后绘图

Plotting after Doing Simulation of Linear Regression with R

simulation

r

linear-regression