tentativa de calcular poder num modelo de regressao logistica simples, com reamostragem

poder_logistica_reamostragem

# Calcula a previsão de um modelo. O mesmo que "fitted", 
# mas aceita um modelo com coeficientes modificados
pred <- function(model) {
  pred <- model.matrix(model) %*% model$coefficients
  pred <- exp(pred) / (1+exp(pred))
  return(pred)
}

# Calcula o valor Z do teste de hipóteses para uma variável no modelo glm
zval <- function(model, v) {
  p <- fitted(model)
  V <- diag(p * (1-p))
  X <- model.matrix(model)
  std_beta <- sqrt(diag(solve(t(X) %*% V %*% X)))
  zval <- abs(model$coefficients[v] / std_beta[v])
  return(zval)
}

# Gerando dados e ajustando modelo
N <- 1000
x <- rbinom(N, 1, .5)
y <- ifelse(x==1, rbinom(N, 1, .9), rbinom(N, 1, .4))
d <- data.frame(y, x)
m <- glm(y~x, family=binomial, data=d)

# Aqui eu fixo um valor do efeito de x em y
m$coefficients[['x']] <- log(1)

# Número de simulações
sim <- 1000

# Aqui eu replico a análise SEM reamostragem
a1 <- replicate(sim, {
  m$model[['y']] <- rbinom(N, 1, pred(m))
  m_new <- glm(formula(m), data=m$model, family=binomial)
  zv <- zval(m_new, 'x')
  valor_p <- 2 * pnorm(zv, lower.tail=F)
  return(ifelse(valor_p < .05, 1, 0))   
})
# Aqui eu replico a análise COM reamostragem (com reposição)
a2 <- replicate(sim, {
  am <- sample(1:N, replace=T)
  m$model[['y']][am] <- rbinom(length(am), 1, pred(m)[am])
  m_new <- glm(formula(m), data=m$model[am,], family=binomial)
  zv <- zval(m_new, 'x')
  valor_p <- 2 * pnorm(zv, lower.tail=F)
  return(ifelse(valor_p < .05, 1, 0))   
})

# Deveria dar 5%, pois eu fixei o efeito como nulo
mean(a1)
mean(a2)

# Por que quando eu faço reamostragem a proporção de rejeição aumenta?