模拟逻辑回归二元结果数据的正确方法？

Question

我正在尝试进行一项涉及 3 个协变量和基于这三个协变量的二元结果的模拟研究。我所做的是基于正态分布模拟三个连续变量 x1,x2,x3，并为每个协变量分配系数 $gamma$。然后我使用函数 $exp(x^T\gamma)/(1+exp(x^T\gamma))$ 生成 $\pi$ 的值，然后使用基于 $\pi 的二项式分布生成二进制结果$。我的任务是测试在 R 中使用

glm()

和族“二项式”是否可以很好地估计系数。然而，估计系数与真实系数并不是很接近。但是，如果我的数据模拟过程正确，我认为这种情况不应该发生。

我已经在R中附加了我的代码。我使用没有拦截的逻辑模型的原因是因为我在生成

t

时没有假设拦截。

r <- c(0.33,3.75,2.82)
# specifying x when there's no correlation

x1 <- rnorm(500,1,1)
x2 <- rnorm(500,2,1)
x3 <- rnorm(500,5,1)

df2 <- data.frame(x1=x1,x2=x2,x3=x3)
x <- as.matrix(df2)

df2 <- df2 %>%
  mutate(l = x %*% r,
         p = 1/(1+l),
         t = rbinom(n(),1,t))

coef_corr_2 <- glm(t ~ x1+x2+x3-1,data=df2,family="binomial") %>% coef()

估计的系数值为

> glm(t ~ x1+x2+x3-1,data=df2,family="binomial") %>% coef()
          x1           x2           x3 
-0.008725804 -0.022818767  0.083695370 

这与真实值相去甚远。

我想知道我在数据生成阶段是否做错了什么？我希望有人能够发现我的代码或我的推理过程中的问题，这实际上只是通过如此简单的生成过程，

glm()

应该给出非常好的系数估计。

Answer 1

当我必须这样做时，我更喜欢这样做：

set.seed(42)

n <- 10000

x1 = rnorm(n)
x2 = rnorm(n)
x3 = rnorm(n)
z = 0 + 1 * x1 + 2 * x2 + 3 * x3
p = 1 / (1 + exp(-z))
y = rbinom(n, 1, p)

data = data.frame(y, x1, x2, x3)

glm(y ~ -1 + x1 + x2 + x3, data = data, family = "binomial") |>
  coef()
#>       x1       x2       x3 
#> 1.000824 2.002293 3.005233

^{创建于 2024-04-16，使用 reprex v2.1.0}