I'm trying to simulate 2 data sets for testing some variable selection methods for logistic regression models.
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.
Specifictly,
-
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 ...