计算给定分布中的经验 Beta

Question

我必须计算经验β，并通过测量差异将其与理论β进行比较。这是开发新代码的模型：

if (!require("pwr")) install.packages("pwr")

# parameters
N=1000
n = 25
k = 500
a_teo = 0.05

#Here we create the universe 
poblacion <- rnorm(N, 10, 10)
m.pob <- mean(poblacion)
sd.pob <- sd(poblacion)

#Data simulation and store of p.values
p <- vector(length=k)
for (i in 1:k){
 muestra <- poblacion[sample(1:N, n)]
 p[i] <- t.test(muestra, mu=15)$p.value
}

#Empiric beta
beta_emp <- length(p[p>a_teo])/k

#Theoretical beta
d=(10 - 15) / sd.pob
beta_teo <- 1 - pwr.t.test(n,
                           d=d,
                           sig.level=a_teo,
                            type="one.sample")$power
#Print theoretical vs empirical bet
sprintf("beta_teo = %.3f -- beta_emp = %.3f -- Diff = %.3f", beta_teo, beta_emp, beta_teo-beta_emp)

使用前面的内容作为参考，我必须运行一个非“正态”分布的模拟，而只是具有范围（2、3、4、5）内的值，并测量经验和理论 beta 之间的差异。这是我的代码：

if (!require("pwr")) install.packages("pwr")

N=1000
n = 25
k = 500
a_teo = .05

poblacion <- sample(2:5, N, replace = T)
m.pob <- mean(poblacion)
sd.pob <- sd(poblacion)

p <- vector(length=k)
for (i in seq_along(p)){
  muestra <- sample(2:5, N, replace = T)
  p[i] <- t.test(muestra, mu=m.pob)$p.value
}

#Empiric beta
beta_emp <- length(p[p>a_teo])/k

#Theoretical beta
d=(10 - m.pob) / sd.pob
beta_teo <- 1 - pwr.t.test(n,
                           d=d,
                           sig.level=a_teo,
                           type="one.sample")$power
sprintf("beta_teo = %.3f -- beta_emp = %.3f -- Diff = %.3f", beta_teo, beta_emp, beta_teo-beta_emp)

我得到“零”作为理论贝塔值。有人能告诉我问题出在哪里吗？ “10”值是示例中的平均值，因此可能应该对此进行更改。

Answer 1

有几个容易被误解的地方。

• 您必须使用

  rnorm

  模拟高斯（=正常）样本，但您使用

  sample

• 效应大小为

  (mu - mu1) / sd

  ，其中

  mu

  是模拟样本的实际总体平均值，

  mu1

  是备择假设下的总体平均值

library(pwr)

n = 20   # sample size
k = 5000 # number of simulated samples
a_teo = .05

m.pob <- 2 # population mean under H1
mu <- 1    # actual population mean
stdev <- 1.5

p <- numeric(k)
for (i in seq_along(p)){
  muestra <- rnorm(n, mu, stdev)
  p[i] <- t.test(muestra, mu = m.pob)$p.value
}

#Empiric beta
beta_emp <- length(p[p>a_teo])/k

#Theoretical beta
d <- (mu - m.pob) / stdev
beta_teo <- 1 - pwr.t.test(n,
                           d = d,
                           sig.level = a_teo,
                           type = "one.sample")$power
sprintf("beta_teo = %.3f -- beta_emp = %.3f -- Diff = %.3f", beta_teo, beta_emp, beta_teo-beta_emp)