如何在 R 中为多元回归运行蒙特卡罗模拟？答案

【问题标题】：How to run a montecarlo simulation for multple regression in R?如何在 R 中为多元回归运行蒙特卡罗模拟？
【发布时间】：2021-02-04 19:25:20
【问题描述】：

我想为预测 mpg 的多元回归模型运行蒙特卡罗模拟，然后评估每辆车的性能比另一辆车更好（更低 mpg）的次数。这就是我目前所得到的

library(pacman)
pacman::p_load(data.table, fixest, stargazer, dplyr, magrittr)

df <- mtcars
fit <- lm(mpg~cyl + hp, data = df)
fit$coefficients[1]

beta_0 = fit$coefficients[1] # Intercept 
beta_1 = fit$coefficients[2] # Slope (cyl)
beta_2 = fit$coefficients[3] # slope (hp)
set.seed(1)  # Seed
n = 1000     # Sample size
M = 500      # Number of experiments/iterations

## Storage 
slope_DT <- rep(0,M)
slope_DT_2 <- rep(0,M)
intercept_DT <- rep(0,M)

## Begin Monte Carlo

for (i in 1:M){ #  M is the number of iterations
  
  # Generate data
  U_i = rnorm(n, mean = 0, sd = 2) # Error
  X_i = rnorm(n, mean = 5, sd = 5) # Independent variable
  Y_i = beta_0 + beta_1*X_i + beta_2*X_i +U_i  # Dependent variable
  
  # Formulate data.table
  data_i = data.table(Y = Y_i, X = X_i)
  
  # Run regressions
  ols_i <- fixest::feols(data = data_i, Y ~ X)
  
  # Extract slope coefficient and save
  slope_DT_2[i] <- ols_i$coefficients[3]
  slope_DT[i] <- ols_i$coefficients[2]
  intercept_DT[i] <- ols_i$coefficients[1]
  
}


# Summary statistics
estimates_DT <- data.table(beta_2 = slope_DT_2,beta_1 = slope_DT, beta_0 = intercept_DT)

此代码不会为hp 创建任何系数我想知道如何将系数添加到模型中，然后预测结果并测试一辆车的 mpg 比另一辆车低多少次。例如，马自达 RX4 的预测 mpg 比 Datsun 710 低多少倍。关于如何使这项工作的一些想法？谢谢

【问题讨论】：

我认为它只是错字。在每次迭代结束时，您覆盖slope_DT。其中之一必须是slope_DT_2。
嗨。感谢您的观察，我修复了类型，但仍然没有给出我需要的结果。关于如何解决这个问题的任何想法？
您只创建了一个独立变量X_i。我认为您应该创建其中两个，一个用于cyl (X_i_1) 和一个 fpr hp (X_i_2)。然后你可以做data_i = data.table(Y = Y_i, X1 = X_i_1,X2=X_i_2)，公式变成ols_i <- fixest::feols(data = data_i, Y ~ X1 + X2)。然后你也会得到三个 coefficeints 而不是当前代码中的两个。
那行得通，谢谢！。现在我如何使用它来进行预测并评估一辆车的 mpg 比另一辆车低多少次？
那么你想比较每对可用的汽车吗？

标签： r regression simulation montecarlo

【解决方案1】：

就像我在评论中指出的那样，您应该使用两个自变量。此外，我想向您推荐lapply-function，它使代码更短，因为您不需要初始化/存储部分。

estimates_DT <- do.call("rbind",lapply(1:M, function(i) {
  # Generate data
  U_i = rnorm(n, mean = 0, sd = 2) # Error
  X_i_1 = rnorm(n, mean = 5, sd = 5) # First independent variable
  X_i_2 = rnorm(n, mean = 5, sd = 5) #Second ndependent variable
  Y_i = beta_0 + beta_1*X_i_1 + beta_2*X_i_2 + U_i  # Dependent variable

  # Formulate data.table
  data_i = data.table(Y = Y_i, X1 = X_i_1, X2 = X_i_2)
  
  # Run regressions
  ols_i <- fixest::feols(data = data_i, Y ~ X1 + X2)  
  ols_i$coefficients
}))

estimates_DT <- setNames(data.table(estimates_DT),c("beta_0","beta_1","beta_2"))

要比较两辆车的预测，请定义以下函数，将您想要比较的两个车名作为参数：

compareCarEstimations <- function(carname1="Mazda RX4",carname2="Datsun 710") {
  car1data <- mtcars[rownames(mtcars) == carname1,c("cyl","hp")]
  car2data <- mtcars[rownames(mtcars) == carname2,c("cyl","hp")]
  
  predsCar1 <- estimates_DT[["beta_0"]] + car1data$cyl*estimates_DT[["beta_1"]]+car1data$hp*estimates_DT[["beta_2"]]
  predsCar2 <- estimates_DT[["beta_0"]] + car2data$cyl*estimates_DT[["beta_1"]]+car2data$hp*estimates_DT[["beta_2"]]
  
  list(
    car1LowerCar2 = sum(predsCar1 < predsCar2),
    car2LowerCar1 = sum(predsCar1 >= predsCar2)
  )
}

确保作为参数提供的名称是有效名称，例如在rownames(mtcars)。

【讨论】：