模拟一个简单的线性模型

Question

我正在尝试模拟一个简单的线性模型 100 次，并从线性模型中找到 B1 的 LS 估计。

set.seed(123498)
x<-rnorm(z, 0, 1)
e<-rnorm(z, 0 ,2)
y<-0.5 + 2*x + e
model<- lm(y~x)
simulaten=100
z=10

for (i in 1:simulaten){
  e<-rnorm(n, 0 ,2)
  x<-rnorm(n, 0, 1)
  y<-0.5 + 2*x + e
  model<- lm(y~x)}
  summary(model)

这是我的 for 循环实现的目标还是我错过了目标？

Answer 1

model 在每次迭代中更新。所以 summary 返回最后一个 'model' 的汇总输出。我们可以将其存储在 list.

中

# // initialize empty list of length equals length of simulaten
modellst <- vector('list', simulaten)
 for(i in seq_len(simulaten)) {
      e <- rnorm(n, 0 ,2)
      x <- rnorm(n, 0, 1)
      y <- 0.5 + 2*x + e
      # // assign the model output to the corresponding list element
      modellst[[i]] <- lm(y~x)  
 }

Answer 2

这是一个replicate解决方案。我已将 n（在问题中忘记）和 simulaten 设置为较小的值。

n <- 100
simulaten <- 4

set.seed(123498)

model_list <- replicate(simulaten, {
      e <- rnorm(n, 0, 2)
      x <- rnorm(n, 0, 1)
      y <- 0.5 + 2*x + e
      lm(y ~ x)
}, simplify = FALSE)

model_list

编辑

可以从模型列表中获得一些统计数据。使用应用于每个模型的函数 coef 提取系数。
使用 sapply 完成，返回的对象是一个 2 行矩阵。

betas <- sapply(model_list, coef)

str(betas)
# num [1:2, 1:1000] 0.671 1.875 0.374 2.019 0.758 ...
# - attr(*, "dimnames")=List of 2
#  ..$ : chr [1:2] "(Intercept)" "x"
#  ..$ : NULL

至于图表，这里有一个例子。请注意，为了使 x 轴达到所有 x 值，在第一次调用 hist 时，参数 xlim 设置为 range(betas).

lgd <- c(expression(beta[0]), expression(beta[1]))
hist(betas[1, ], freq = FALSE, col = "lightblue", xlim = range(betas), ylim = c(0, 2.5), xlab = "betas", main = "")
hist(betas[2, ], freq = FALSE, col = "blue", add = TRUE)
legend("top", legend = lgd, fill = c("lightblue", "blue"), horiz = TRUE)

模拟一个简单的线性模型

simulating a simple linear model

simulation

r

linear-regression

编辑