模拟的争论结果

编程问答 2022-05-16

问题描述

我对 R 很陌生，我正在尝试学习如何进行一些模拟。目前我有一个程序可以执行以下操作：

  1. 在一个函数中，使用 DGP 创建假数据，该数据以 tibble 形式返回
  2. 在另一个函数中，将假观察随机分配给治疗
  3. 在最终函数中，将随机分配结果与假数据合并并运行回归。我使用以下代码返回一个包含估计值和 p 值的列表
    tauhat <- trobust %>% filter(term=="TTRUE") %>% pull(estimate)
    pvalue <- trobust %>% filter(term=="TTRUE") %>% pull(p.value)
    return(list(tauhat,pvalue))

如果我运行这些函数一次，我会得到类似下面的内容

> finitepop(finiteN=20)
[[1]]
[1] 0.3730686

[[2]]
[1] 0.03445962

然后我使用复制来重复这个过程，比如 100 次。我最终得到了一个 2X100 的东西 - 也许它是一个数组？ - 我想把它变成一个 100X2 的小标题，也就是说，一个小标题，其中包含估计值和 p 值以及作为观察值存储的模拟结果的列。模拟结果看起来像

> finitesim <- (replicate(n=reps,finitepop(finiteN=20)))
> finitesim
     [,1]        [,2]      [,3]      [,4]       [,5]      [,6]      [,7]     
[1,] -0.03096849 0.206797  0.2386698 0.09374408 0.1462773 0.2479394 0.2177207
[2,] 0.8850678   0.2622687 0.2105784 0.5990369  0.3279901 0.1063231 0.2489028
     [,8]      [,9]       [,10]       [,11]     [,12]      [,13]     
[1,] 0.1661424 0.00977172 -0.08761129 0.1170922 -0.1559203 0.278062  
[2,] 0.2086819 0.9390261  0.6071284   0.472165  0.4214389  0.05973561

我应该如何将结果转换为漂亮的小标题？

编辑：下面是一个 MWE，为了方便起见，我将右侧变量更改为 x，并且我没有为 lm_robust 创建聚类结构

library(tidyverse)
library(lmertest)                       #for lmer
library(merTools)                       #for lmer
library(estimatr)                       #for cluster robust se


finitepop <- function(finiteN){
    
    fakedata <- tibble(
        id=1:finiteN,x=rnorm(n=finiteN),y=rnorm(n=finiteN)
        )
    robust <- lm_robust(data=fakedata,y~x,cluster=id)
    trobust <- tidy(robust)
    tauhat <- trobust %>% filter(term=="x") %>% pull(estimate)
    pvalue <- trobust %>% filter(term=="x") %>% pull(p.value)
    return(list(tauhat,pvalue))
    }


finitesim <- (replicate(n=10,finitepop(finiteN=20),simplify=FALSE))

finitesim

解决方法

使用可以使用 map_df 包（purrr 的一部分）中的 tidyverse：

finitepop <- function(finiteN){
  
  fakedata <- tibble(
    id=1:finiteN,x=rnorm(n=finiteN),y=rnorm(n=finiteN)
  )
  robust <- lm_robust(data=fakedata,y~x,cluster=id)
  trobust <- tidy(robust)
  tauhat <- trobust %>% filter(term=="x") %>% pull(estimate)
  pvalue <- trobust %>% filter(term=="x") %>% pull(p.value)
  print(pvalue)
  return(tibble(tauhat,pvalue))
}

finitesim <- replicate(n=10,finitepop(finiteN=20),simplify=FALSE) %>%
  purrr::map_df(as.data.frame)

> finitesim
         tauhat     pvalue
1   0.035057186 0.89818890
2  -0.248569087 0.24159959
3   0.111054217 0.75700470
4   0.596779950 0.00223398
5  -0.004052686 0.98418837
6  -0.105390590 0.67410417
7  -0.107913504 0.54778478
8  -0.021681712 0.89834059
9  -0.161811559 0.49091499
10  0.241477999 0.21281508
rrsimulationtibble