使用 R 的 Montecarlo 模拟（如何为 R“翻译”Stata 代码）

Question

我需要从 Stata 到 R 运行相同的练习（蒙特卡罗模拟）。

我在 Stata 中使用的代码是下面的代码。我如何使用 R 做到这一点？（我已经搜索了很多教程，但我仍然没有设法在 R 中做到这一点）。

* Simulations (10, 100 and 1000 sample replications/iterations)

clear
drop _all
set obs 100
set seed 54231
gen x = ((rnormal()))*10 + 40

* Generating true_y, considering Beta = 0,035
 
gen true_y = 5+0.03500*x
save truth, replace
twoway scatter true_y x

program hetero1
version 13
args c
use truth, clear
gen y = true_y + rnormal()
regress y x
end


foreach i in 10 100 1000 {
simulate _b, reps (`i'): hetero1
sum _b_x
twoway histogram _b_x, fraction xline(+0.03500, lcolor(red)) xline(`r(mean)', lcolor(green)) fcolor(none) lcolor(gs0) legend(off) title(`i' Repetições)
graph save graf`i'.gph, replace
}
 gr combine graf10.gph graf100.gph graf1000.gph
 graph export "graf.png", as(png) replace 

最后，考虑到 10、100 和 1000 个样本重复，我需要获得这 3 个直方图（估计的 beta/系数）。红线指的是“真实”系数，绿线是估计系数的平均值 - [见链接中的图片]

Answer 1

这应该这样做：

# Simulations (10, 100 and 1000 sample replications/iterations)
library(ggplot2)
library(dplyr)
library(gridExtra)
n <- 100
set.seed(54231)
x <- rnorm(n)*10 + 40

# Generating true_y, considering Beta = 0,035

true_y <- 5+0.03500*x

plot(x, true_y)

b <- t(replicate(1110, coef(lm(true_y + rnorm(n) ~ x))))

b <- as.data.frame(b) %>% 
  rename("a" = "(Intercept)", 
         "b" = "x") %>% 
  mutate(
    obs = 1:n(), 
    n = case_when(
      obs %in% 1:10 ~ "N = 10", 
      obs %in% 11:110 ~ "N = 100", 
      TRUE ~ "N = 1000"), 
    n = factor(n, levels=c("N = 10", "N = 100", "N = 1000")))
  

b10 <- b %>% filter(n == "N = 10") 
g1 <-   ggplot() + 
  geom_histogram(data = b10, aes(x=b), bins=3, col="white") + 
  geom_vline(xintercept = 0.03500, col="red") + 
  geom_vline(data = b10 %>% summarise(b=mean(b)), aes(xintercept = b), col="green") + 
  facet_wrap(~n) + 
  theme_bw()

b100 <- b %>% filter(n == "N = 100") 
g2 <-   ggplot() + 
  geom_histogram(data = b100, aes(x=b), bins=10, col="white") + 
  geom_vline(xintercept = 0.03500, col="red") + 
  geom_vline(data = b100 %>% summarise(b=mean(b)), aes(xintercept = b), col="green") + 
  facet_wrap(~n) + 
  theme_bw()

b1000 <- b %>% filter(n == "N = 1000") 
g3 <-   ggplot() + 
  geom_histogram(data = b1000, aes(x=b), bins=25, col="white") + 
  geom_vline(xintercept = 0.03500, col="red") + 
  geom_vline(data = b1000 %>% summarise(b=mean(b)), aes(xintercept = b), col="green") + 
  facet_wrap(~n) + 
  theme_bw()

library(gridExtra)
grid.arrange(g1, g2, g3, nrow=2)