Because I'm frankly tired of answering questions about stepwise without something of my own to point to, consider the following.

I'm going to simulate a logistic regression with 10 variables. The variables $x_1, x_2, x_3$ are all independent and have log odds ratios of 0.1, 0.2, and 0.5.  

The variables $x_4, x_5, x_6$ have no effect on the log odds, but are correlated with the variables $x_1, x_2, x_3$ like

$$ \operatorname{Cor}(x_j, x_{j+3}) = 0.3\cdot j $$

So $x_1$ and $x_4$ are correlated, but $x_2$ and $x_4$ are not.

Lastly,the variables $x_7, x_8, x_9$ are independent of all others and have no effect on the log odds.

Those who use stepwise regression seem to think that it can select relevant variables.  So if that were true, then surely stepwise regression could select the right variables for this problem...right?

In 1000 simulations from this process, using 1000 observations...

* $x_1$ is selected 37% of the time
* $x_2$ is selected 70% of the time, and
* $x_3$ is selected 86% of the time\

It would appear that larger effects are being selected with larger frequency.  But how often is the trie model selected?  a whopping 9% of the time.

Let me repeat that.  Almost 90% of the time, you're selecting the wrong model; you're including a variable which actually has 0 impact on the outcome.

There are a host of problems with stepwise regression (I'll link them [here](https://www.stata.com/support/faqs/statistics/stepwise-regression-problems/) for you to read).  Its very clear that there is more nuance to what stepwise regression is doing. It isn't selecting the right model anywhere near enough to justify its use.

Now, I think people may take objection with what I've argued here. "*Demetri, we might not even select the right model even using our scientific judgement.  It isn't fair to criticize stepwise regression on those grounds*".  Ok, maybe, but that is one nail removed from the coffin and some 10 more from that list I've linked.  

Anyway, don't use stewpwise.  


### Code 

```r
library(tidyverse)


N <- 1000
q <- 9
Sigma <- diag(q)
Sigma[1, 4] <- Sigma[4, 1] <- 0.3
Sigma[2, 5] <- Sigma[5, 2] <- 0.6
Sigma[3, 6] <- Sigma[6, 3] <- 0.9

results <- map_dfr(1:1000, ~{
  # Simulate data
  X <- MASS::mvrnorm(N, mu=rep(0, q), Sigma = Sigma)
  beta <- rep(0, 9)
  beta[1:3] <- c( 0.1, 0.2, 0.5)
  y <- rbinom(N, 1, plogis(X%*% beta - 2))
  d <- as_tibble(X) %>% 
       mutate(y=y)
  
  # Fit a model and do stepwise regression
  full.fit <- glm(y~.,, data=d, family = binomial)
  step.fit <- MASS::stepAIC(full.fit, direction = 'both', trace = F)
  
  # Grab which variables were selected
  selected_var <- names(coef(step.fit))
  # Determine if the right model was selected
  correct_model <- c( "(Intercept)","V1", "V2", "V3")
  correct_model_selected <- identical(selected_var,correct_model)
  
  # Determine which of the variables were selected from the sample
  vrs <- str_c("V", 1:10) 
  outcomes<-vrs %in% selected_var
  
  names(outcomes) <- vrs
  
  as.data.frame(t(outcomes)) %>% 
    mutate(correct_model = correct_model_selected)
})

summarise_all(results, mean)

```