Counterintuitive coefficients in elastic net logistic regression

Question

In a model run of elastic net logistic regression, I encountered a very counterintuitive coefficient. I have looked into the data, model and script, but, I still cannot seem to wrap my head around the counter-intuitiveness I see regarding the dependent and independent variable. Initially, the V7 should be negative, as it is significantly lower in the dependent variable where the outcome is 1 compared to the outcome of 0, see graph.

Further, the descriptive statistics are:

               Value 0          Value 1
count    749304.000000       402.000000
mean          2.762876         1.618396
std           3.672386         2.488794
min           0.000000         0.000000
25%           0.306000         0.001500
50%           1.662000         0.638250
75%           3.901500         2.338500
max         223.084500        17.217000

But, I end up with coefficients that show the following, here, one should keep an eye on variable number 7 (V7), which I am talking about.

(Intercept) -3.096141e+01
V1           1.436113e-03
V2          -1.774919e-01
V3          -5.586214e-04
V4          -1.763915e-03
V5           6.817795e-03
V6           3.986299e-02
**V7         3.085392e-02**
V8          -1.117509e-02
V9           6.917977e-02

1. Why do I see that coefficient V7 is positive when it clearly is smaller in cases of 1 than cases of 0 in the dependent variable?
2. Do I misinterpret the results of my elastic net regression? I doubt it, as the other variables are intuitively correct?

The script is below:

library(readr)
library(caret)
library(tidyverse)
library(glmnet)
library(ROCR)
library(pROC)
library(doParallel)
registerDoParallel(4, cores = 4)
set.seed(123)
df <- read_csv("df.csv")
View(df)
training.samples <- df$V10 %>% createDataPartition(p = 0.8, list = FALSE)
train <- df[training.samples, ]
test <- df[-training.samples, ]
x.train <- data.frame(train[, names(train) != "V10"])
x.train <- data.matrix(x.train)
y.train <- train$fire
x.test <- data.frame(test[, names(test) != "V10"])
x.test <- data.matrix(x.test)
y.test <- test$fire
nFolds <- 10
foldid <- sample(rep(seq(nFolds), length.out = nrow(train)))
list.of.fits <- list()
for (i in 0:10){
    fit.name <- paste0("alpha", i/10) 
    list.of.fits[[fit.name]] <- cv.glmnet(x.train, y.train, type.measure = "deviance", alpha = i/10, family = "binomial", nfolds = nFolds, foldid = foldid, parallel = TRUE)
}
coef <- coef(list.of.fits[[fit.name]], s = list.of.fits[[fit.name]]$lambda.min)
coef

Answer 1

It is not uncommon for coefficients to change sign when other terms are added to the model, e.g. when x7 is the only term in the model then the coefficient will be negative as you expect, but as you add more terms (either through a stepwise approach, or relaxing the penalties/constraints in the elasticnet) the coefficient on x7 can easily change sign due to relationships between x7, the other predictors, and the response.

See this answer: Binary Logistic Regression: Direction of B's different in multiple than in bivariate cases for one example/explanation that may make this more clear.