Understanding log odds ratio from quartiles in logistic regression Data is found by running data(PimaIndiansDiabetes) from the mlbench package.
PimaIndiansDiabetes <- PimaIndiansDiabetes %>% na.omit()

logmod4 <- glm(diabetes ~ pregnant + glucose + mass + pedigree, 
     family="binomial", data=PimaIndiansDiabetes)

# Extracting 1st & 3rd quantiles from PIMA dataset.
BMI_1Q <- quantile(PimaIndiansDiabetes$mass, 0.25)
BMI_3Q <- quantile(PimaIndiansDiabetes$mass, 0.75)

# Getting the 95% confidence interval for the bmi parameter.
BMI_Conf <- confint(logmod4, 'mass', level = 0.95)

BMI_Conf

# Getting the 95% confidence interval for our previous log-odds ratio.
odds_ratio <- (exp(BMI_Conf * (BMI_1Q - BMI_3Q)))

# Assuming that all other variables remain constant, we get a probability of:
odds_ratio/(1+odds_ratio)

##### Output:
####     2.5 %  ,    97.5 % 
#### 0.4155208  , 0.2599249 

Can someone tell me the reason why the following statement is true?

Keeping other factors constant, women with a BMI at the 1st quartile are between 26% to 42% less prone to obtain diabetes than women with a BMI at the 3rd quartile.

mass = BMI readings, by the way.

The interquartile range for BMI is $36.6-27.3 = 9.3$, so a correct interpretation of the confidence interval would be:

Keeping other variables constant, the data are compatible at the 95% level with an
increase of the odds for diabetes by a factor of $1.62$ to $2.67$
comparing a women with a BMI at the first quartile with one at the
third quartile (i.e. for a difference in BMI of $9.3$ units).

The most plausible value (the point estimate) for the odds ratio is $\exp(0.07810\cdot 9.3) = 2.07$, so roughly a doubling of the odds for diabetes.
If you are interested in the difference for the predicted probabilities and a corresponding confidence interval, you could apply a technique described in this thread. What makes this not straightforward is the fact the you need to assume values for the other variables in the model. A common choice is to fix them at their mean values. Here is an example using the emmeans package:
em <- emmeans(logmod4, "mass", regrid = "response", at = list(mass = c(BMI_3Q, BMI_1Q)))

em

mass  prob     SE  df asymp.LCL asymp.UCL
 36.6 0.378 0.0245 Inf     0.330     0.426
 27.3 0.227 0.0220 Inf     0.184     0.270

contrast(em, "pairwise", infer = c(TRUE, TRUE))

contrast            estimate     SE  df asymp.LCL asymp.UCL z.ratio p.value
mass36.6 - mass27.3    0.151 0.0258 Inf       0.1     0.201   5.847  <.0001

The predicted probabilities at the first and third quartile of BMI are $0.227$ and $0.378$, respectively. The estimated difference is thus $0.151$ with an asymptotic 95% confidence interval of $(0.100;\,0.201)$.
