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 ...


1 Answer 1


Your code looks correct. I am not sure to understand your problem. For a given set of true parameters, the size of the standard errors depends on the sample size and the size of the parameters. In your first case, the SE for X1 is ~3.9% the size of the MLE for X1 (0.02064601 / 0.528894792 = 0.039). In the second case, you obtain a comparable value (i.e., 3.4%). This result makes sense because the simulation routine is exactly the same for both cases. The significance level (P-val) will be largely impacted by the size of the true parameters. For a given number of observations, it will be more difficult to show that the estimate differs from 0, that's why P-value for X1 is smaller in 2nd case than in 1st case.

  • $\begingroup$ Thank you so much. I'm not sure I define my question properly. In both cases, the true parameters are both 5 (out of 100 candidate predictors) and the observations are also the same (sample size = 20000). That's why I don't get why the significance levels are not exactly what I expect. $\endgroup$
    – d.lee
    Mar 2, 2019 at 5:18
  • $\begingroup$ Looking at your example, it seems that the true values of the parameters in the 2nd case are multiplied by 10 (true value for X1 = 0.5 in case 1 vs 5 in case 2) - This alone would explain the difference in P-val across the 2 cases - Simulate data for a 3rd case where true values would be 10 times smaller than the ones used in case 1 (so X1 = 0.5 would become 0.05) and again you should obtain larger P-values (Everything else being equal) - This has something to do with the variance of the binomial distribution - By making your true values smaller decrease the signal-to-noise ratio of your model $\endgroup$
    – Nicolas K
    Mar 2, 2019 at 8:42
  • $\begingroup$ I think I sort of understand your point. But if my goal is to test a variable selection method for the logistic regression models, the setting of the original 2 cases does not seem to be good enough for this. Since when all 100 predictors are included in the model at initial step, more predictors tend to be selected for case 2 setting in the simulation procedure. I'm not sure if this is caused by the phenomenon you describe or something else. $\endgroup$
    – d.lee
    Mar 2, 2019 at 12:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.