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!
I appreciate it in advance.
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
...
