如何多次重复模拟？

Question

我对 Rstudio 真的很陌生，所以我希望有人可以帮助我。 所以我有这个代码：

x = 1:5
alpha = 1
beta = 1.5
betaD = 0.1
s = 1
sa = 0.2
sb = 0.2
N = 10

grp = factor(rep(c("Control", "Treatment"), c(N,N)))

for(i in 1:(2*N)) {
  ai = rnorm(1, 0, sa)
  bi = rnorm(1, 0, sb)
  intercept = alpha+ai
  slope = beta + bi + ifelse(grp[i]=="Treatment", betaD, 0.0)

  y = intercept+ slope*x + rnorm(length(x), 0, s)

  tmp = data.frame(subject=i, x=x, y=y, a=ai, b=bi, group=grp[i])
  if(i==1) dataset = tmp
  else dataset = rbind(dataset, tmp)
}

require(lme4)

fitAll= lmList(y~x|subject, data=dataset)
slopes = coef(fitAll)$x
boxplot(slopes~grp)
t.test(slopes~grp, var.equal=TRUE)

fit0 = lmer(y~ x +(x|subject), data=dataset, REML=FALSE)
fit1 = lmer(y~ group*x +(x|subject), data=dataset, REML=FALSE)
anova(fit0, fit1)

当我运行它时，它会生成这个：

Two Sample t-test

data:  slopes by grp
t = -2.2495, df = 18, p-value = 0.03723
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.66690111 -0.02277686
sample estimates:
  mean in group Control mean in group Treatment 
               1.362975                1.707814

还有这个：

Data: dataset
Models:
fit0: y ~ x + (x | subject)
fit1: y ~ group * x + (x | subject)
     Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)  
fit0  6 326.65 342.28 -157.32   314.65                           
fit1  8 324.34 345.18 -154.17   308.34 6.3072      2     0.0427 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

基本上我想要做的是在代码中重复，这样当我点击运行时，它会生成这个，无论我指定多少次。然后我希望它将它生成的 p 值分为两组，一组 p 值高于 0.05，另一组低于 0.05

正如我所说，我对此真的很陌生，所以如果有人可以简单地向我解释，那将不胜感激。

Answer 1

要多次运行代码，请使用replicate。类似的东西

replicate(
  100,
  {
     # Your code that creates the random dataset and runs ANOVA
  }
)

Answer 2

为了简单起见，我从t.test 中获取了 p 值，它可能不是您要的 p 值。但是，它适用于演示目的。

只需将您的代码包装在一个函数中，并根据需要多次使用replicate：

do_once <- function()
{
  x = 1:5
  alpha = 1
  beta = 1.5
  betaD = 0.1
  s = 1
  sa = 0.2
  sb = 0.2
  N = 10

  grp = factor(rep(c("Control", "Treatment"), c(N,N)))

  for(i in 1:(2*N)) {
    ai = rnorm(1, 0, sa)
    bi = rnorm(1, 0, sb)
    intercept = alpha+ai
    slope = beta + bi + ifelse(grp[i]=="Treatment", betaD, 0.0)

    y = intercept+ slope*x + rnorm(length(x), 0, s)

    tmp = data.frame(subject=i, x=x, y=y, a=ai, b=bi, group=grp[i])
    if(i==1) dataset = tmp
    else dataset = rbind(dataset, tmp)
  }

  require(lme4)

  fitAll= lmList(y~x|subject, data=dataset)
  slopes = coef(fitAll)$x
  boxplot(slopes~grp)
  t.test(slopes~grp, var.equal=TRUE)$p.value  
}
p_vals <- replicate(10, do_once())

要获得低于 0.05 的 p 值，只需

p_vals[p_vals < 0.05]

是的，这与 Rstudio 无关，R 代码可以在任何 IDE 和纯 R 控制台中运行。