How to Simulate Logit

> set.seed(1977)
> n <- 1000
> treat <- sample(c("a", "b", "c"), n, T)
> num.diseases <- sample(0:4, n, T)
> age <- rnorm(n, 50, 10)
> cholesterol <- rnorm(n, 200, 25)
> weight <- rnorm(n, 150, 20)
> sex <- sample(c("female", "male"), n, T)
> X <- as.data.frame(cbind(num.diseases, age, cholesterol, treat))
> X$num.diseases <- as.numeric(X$num.diseases)
> X$age <- as.numeric(X$age)
> X$cholesterol <- as.numeric(X$cholesterol)
> X <- as.data.frame(X[(X$cholesterol - 10) - 5.2 > 0, ])
> dim(X)
[1] 985   4
> head(X)
  num.diseases age cholesterol treat
1            1 349          22     b
2            1  49         738     b
3            1 144         962     c
4            4 702         621     a
5            1 383         214     b
6            5 741         914     c
> L <- 0.1 * (X$num.diseases - 2) + 0.045 * (X$age - 50) + (log(X$cholesterol -
+     10) - 5.2) * (-2 * (X$treat == "a") + 3.5 * (X$treat == "b") +
+     2 * (X$treat == "c"))
> y <- ifelse(runif(985) < plogis(L), 1, 0)
> glm(y ~ num.diseases + age + cholesterol + treat, data = X)
Call:  glm(formula = y ~ num.diseases + age + cholesterol + treat, data = X)

Coefficients:
 (Intercept)  num.diseases           age   cholesterol        treatb
   7.831e-01    -6.689e-03     2.640e-04     3.402e-05     5.350e-02
      treatc
   6.425e-02

Degrees of Freedom: 984 Total (i.e. Null);  979 Residual
Null Deviance:      46.56
Residual Deviance: 39.93        AIC: -348.1