# Significance levels of the estimated logistic regression coefficients for artificially generated data sets

I'm trying to simulate 2 data sets for testing some variable selection methods for logistic regression models. The initial step is to fit the logistic regression for all candidate predictors. But something not expected shows.
The data is generated as follows. The response, Y, is from a Bernoulli distribution with the probability $$p_i = { exp(β^TX_i) \over (1 + exp(β^TX_i)},$$ where $$X_i$$ is a 100 ×1 predictor vector and $$\beta$$ is the corresponding parameter vector.

For the ith predictor vector $$X_i = (x_{i,1}, . . . , x_{i,P} )^T$$, $$x_{i,j}$$ are e independently generated from a normal distribution $$N (\mu_j ,1)$$, where $$\mu_j$$ is assumed from the uniform distribution within (−1, 1).

For the two cases, all candidate predictors have the same values for each and every simulation. The coefficients in model 2 are 10 times of those in model 1. Specifically,

• Case 1: The first five parameters are chosen as $$(\beta_1, \beta_2, . . . , \beta_5)^T = (0.5, −2.0, −0.6, 0.5, 1.2)$$, and the other $$\beta_i, i > 5$$, are set as zeros.

• Case 2: The parameter vector is $$(\beta_1, \beta_2, . . . , \beta_5)^T = (5, −20, −6, 5, 12)$$, which is equal to 10 times of the vector in Case 1. The other $$\beta_i, i > 5$$, are set as zeros.

Intuitively, if I run logistic regression as Y~ all of the 100 Xs for the two data sets, the significant levels of the coefficient estimates for case 2 should be higher than those for case 1. However, the estimated standard errors of the coefficients for case 1 are usually much smaller than those in case 2. As a result, the significance levels of the true parameters in case 2 are not higher as I expected. Is it because of the simulation settings or some IRLS estimation issue?
If any expert sees this question, kindly help me with this problem!
My code is as follows.

MyData = function(n, p, param1, param2){
X = matrix(0, nrow = n, ncol = p)
for (i in 1:p){
mu = runif(1, -1, 1)
X[,i] = rnorm(n, mu, 1)
}

para1 = c(param1, rep(0,p-length(param1)))
para2 = c(param2, rep(0,p-length(param2)))
eta1 = X%*%para1
pi1 = exp(eta1)/(1+exp(eta1))
y1 = rbinom(n, 1, prob=pi1)
eta2 = X%*%para2
pi2 = exp(eta2)/(1+exp(eta2))
y2 = rbinom(n, 1, prob=pi2)

data1 = data.frame(y1,X)
names(data1) = c("Y", paste('X',c(1:p),sep=''))
data2 = data.frame(y2,X)
names(data2) = c("Y", paste('X',c(1:p),sep=''))

return(list(case1 = data1, case2 = data2))
}

syn.data = MyData(20000, 100, c(0.5, -2, -0.6, 0.5, 1.2), c(5, -20, -6, 5,12))
case1 = syn.data$$case1 case2 = syn.data$$case2
model.fit.1 = glm(Y ~ ., data = case1, family = binomial, control = list(maxit = 50))
model.fit.2 = glm(Y ~ ., data = case2, family = binomial, control = list(maxit = 50))


And one simulation result is listed as follows.
summary(model.fit.1)$coeff Estimate Std. Error z value Pr(>|z|) (Intercept) -0.033630247 0.11251447 -0.2988971 7.650186e-01 X1 0.528894792 0.02064601 25.6172945 9.790102e-145 X2 -2.031462066 0.03121973 -65.0698095 0.000000e+00 X3 -0.546228462 0.02057548 -26.5475444 2.741460e-155 X4 0.486106719 0.02063916 23.5526400 1.179548e-122 X5 1.233698734 0.02444974 50.4585719 0.000000e+00 ... summary(model.fit.2)$coeff Estimate Std. Error z value Pr(>|z|) (Intercept) 2.773897e-01 0.32665570 0.849180700 3.957808e-01 X1 5.013686e+00 0.16820019 29.807852085 3.090660e-195 X2 -2.008715e+01 0.63427844 -31.669295410 4.114227e-220 X3 -5.925434e+00 0.19477124 -30.422528791 2.766526e-203 X4 5.027093e+00 0.16797930 29.926859266 8.803500e-197 X5 1.207407e+01 0.38530449 31.336426397 1.489435e-215 ...