Stepwise binary logistic regression, Forward or Backward? I want to run a binary logistic regression to understanding (modeling) factors affecting nest-site selection in a bird species. I have Presence/Absence data and 13 predictors. The sample size is small. 32 Presence points and 64 Absence points.
I want to run a stepwise binary logistic regression in R manually. I have 2 questions about it.
Which method can be better? Forward or Backward?
Which statistic is more important for removing or adding variables at each step? AIC or p-value of ANOVA (Wald test) ?

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 parameters. 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 true 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 or excluding a variable which does impact the outcome.  This is especially damning because all assumptions of the model are met as best as they can be met.  All variables affecting the outcome are eligible to be selected and they are all linear in the log odds.  Nothing is wrong with the model, it's the selection that is causing this apparent mal performance.
There are a host of problems with stepwise regression (I'll link them here for you to read).  It's 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 stepwise.
Code
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)


