Predicting osteoporosis risk using nighttime eating habits

This cross-sectional study aimed to investigate the association between nighttime eating exposure (NEE, exposure) and osteoporosis (outcome). The sample size was calculated using the formula for cross-sectional surveys (α = 0.05, Z₂/α = 1.96). As shown in Equation (1).

Based on an osteoporosis patient prevalence p = 39.19% (16), and an allowable error δ = 0.2p = 0.07838, the minimum sample size was calculated as 150. The final sample size was determined to be 150 participants. This study employed consecutive sampling from patients with nighttime eating habits at the Neck, Shoulder, Low Back, and Leg Pain Center of Luoyang Orthopedic Hospital (Henan Provincial Orthopedic Hospital) from January 2024 to May 2025. Inclusion criteria: (1) Patients with nighttime eating habits; (2) Able to undergo both lumbar spine and femoral neck bone mineral density examinations within an interval not exceeding 7 days to minimize external influences and ensure consistency; (3) Age 18 years or older; clear consciousness, possessing normal reading and communication abilities. Exclusion criteria: (1) Presence of congenital scoliosis, spinal fracture, spinal tumor, ankylosing spondylitis; (2) Suffering from diseases severely affecting bone metabolism, such as chronic kidney/liver failure, hyperparathyroidism, rheumatoid arthritis, and other autoimmune diseases; (3) Long-term use of medications affecting bone metabolism, such as glucocorticoids, antiepileptic drugs, etc.; (4) Patients with incomplete questionnaires or clinical data. Osteoporosis diagnosis in enrolled patients followed the Guidelines for Diagnosis and Treatment of Primary Osteoporosis (2022) (17). Informed consent was obtained from all participants prior to any procedures. Before any procedure, all participants provided written informed consent. The study was conducted in accordance with the principles and guidelines of the Declaration of Helsinki and was approved by the Ethics Committee of Henan Luoyang Orthopedic Hospital (Henan Orthopedic Hospital) (approval No. 2023KYKT0021-02).

During the dietary assessment, rigorously trained researchers employed the Retrospective Dietary Survey Auxiliary Reference Food Atlas, developed by Professor Wang Zhixu, alongside food models to assist patients, family members, or caregivers in recalling dietary intake over a period of three days. The average daily nutrient intake was calculated using Nutrition Calculator v2.8.0.5, provided by the Institute of Nutrition and Food Safety, Chinese Center for Disease Control and Prevention, Beijing, China. This study utilized a 24-h dietary recall form, where dietary intake was evaluated through face-to-face interviews conducted over three consecutive days. All patients exhibiting nighttime eating habits were interviewed to describe the types and quantities of foods consumed within the past 24 h. In cases where patient recall proved challenging, family members supplemented the dietary information. All dietary surveys were conducted following patient admission. The timing of food intake was incorporated into the standard 24-h dietary recall form. Due to the absence of a standardized definition in prior research, the NEE window was provisionally defined as occurring between 8:00 p.m. and 6:00 a.m., based on relevant literature and dietary habits. It is important to note that this study is a preliminary exploration of the relationship between NEE and obesity prevalence (OP); therefore, special circumstances, such as night-shift workers, were not considered separately, and only food consumption within this defined time window was included. The calculation formula for NEE is as shown in: Equation (2)

Investigators personally distributed questionnaires and provided detailed explanations to participants regarding the survey's purpose, importance, requirements, and completion methods. All questionnaires were collected immediately upon completion by patients and meticulously checked by investigators for completeness and accuracy. A total of 194 questionnaires were distributed, with 186 valid questionnaires ultimately collected, yielding an effective response rate of 95.6%. Patient disease data were collected and entered by investigators via the hospital's nursing information system, with each entry verified individually.

Data Preparation and Splitting: The osteoporosis status was converted into a binary factor variable (positive/negative). A stratified random sampling strategy (strata = "osteoporosis") was employed to split the dataset into training (70%) and testing (30%) sets at a 7:3 ratio. Reproducibility was ensured using set.seed(123). Stratified sampling guaranteed consistent outcome variable distribution between groups.

Feature Engineering: Predictor variables in the training set were converted to a design matrix via the model.matrix function, with the intercept term removed ([, −1]). All continuous variables were centered and scaled, and zero-variance features (zv) were eliminated.

Feature Selection and Model Construction: A multi-method fusion strategy was used for feature selection to enhance robustness. Model construction: LASSO Regression: The glmnet package was used to fit a LASSO-regularized logistic regression model (family = "binomial," alpha = 1). The λ value was optimized via 10-fold cross-validation. This step was limited to hyperparameter tuning for the LASSO model, aiming to prevent overfitting. Logistic Regression: A standard logistic regression model was fitted using identical preprocessing and 10-fold cross-validation. Random Forest: The randomForest package was used to fit an ensemble model with 500 decision trees (ntree = 500). Predictor importance was calculated based on the mean decrease in Gini impurity. The top 5 important variables from each of the three methods were selected, and their union was taken,finally incorporating 7 variables to construct the final nomogram prediction model.

Model validation: A multi-level validation strategy was employed to assess model performance: the model's performance was evaluated on an independent test set, the Receiver Operating Characteristic curve was calculated to assess discriminative ability, and Decision Curve Analysis and Clinical Impact Curve were used to evaluate clinical utility; additionally, 1,000 Bootstrap resamplings were performed to calculate the model's calibration C-index and Optimism, quantifying the assessment of the model's stability and internal validity.

We selected the NHANES database cycles 2013–2014 and 2017–2018; data were merged initially, but the resulting dataset size was small after direct processing, necessitating the use of random forest imputation for handling missing data. Based on literature review and the focus of this study investigating the relationship between NEE and OP (18), with selection results shown in Supplementary Figure 1. Dietary data were obtained using two non-consecutive 24-h dietary recall interviews following USDA Food and Nutrition Database guidelines. The first interview was conducted in-person, and the second was administered via telephone 3 to 10 days later. Dietary intake data estimated the types and quantities of all foods and beverages (including all types of water) consumed during the 24-h period (midnight to midnight) preceding each interview. Total nutrient intake profiles were generated for each participant, representing daily total energy and nutrient intake from foods and beverages. Time quantification variables were processed using the lubridate package to convert eating time strings (DR1_020) into POSIXct time objects, extracting hour and minute numerical variables. Energy intake statistics were calculated by grouping by subject ID (SEQN): total_kcal represented total daily energy intake; night_kcal represented the sum of energy intake during the nighttime period (night = 1). Data quality control measures included: excluding extreme energy intake (retaining only records with 500 ≤ total_kcal ≤ 5,000 kcal/day); handling zero values (when total_kcal = 0, night_ratio was forced to NA); and coding missing time data (original code 99) as NA. Covariates considered included age, race/ethnicity, education level, household income, marital status, smoking status, alcohol consumption, body mass index (BMI), diabetes status, hypertension status, lumbar spine bone mineral density (BMD), femoral neck BMD, and white blood cell count. The outcome variable was defined as self-reported physician-diagnosed osteoporosis/bone fragility.

NEE served as the exposure variable in this study. As prior research lacks directly relevant exposure variables, no direct genetic instrument for NEE was identified; therefore, indirect exposure variables (sleep deprivation, 24-h dietary patterns, and BMI) were selected to investigate the relationship between NEE and OP using multivariable Mendelian randomization (MVMR). Genetic variants for dietary variables (ukb-b-14294) and BMI (ukb-b-19953) were obtained from European population studies within the UK Biobank; BMI data (ebi-a-GCST006686) were sourced from the European Bioinformatics Institute. Instrumental variable (IV) selection adhered to three core MR assumptions: First, single nucleotide polymorphisms (SNPs) must be strongly associated with the exposure (p < 5 × 10⁻⁶) and exhibit an F-statistic >10. Second, independent SNPs were ensured using linkage disequilibrium (LD) clumping (R² < 0.001, distance = 10,000 kb). Finally, SNPs associated with the outcome or confounding variables were excluded (p < 5 × 10⁻⁶). The OP outcome data (ukb-b-12141) were derived from a European population study within the UK Biobank by Ben Elsworth et al., comprising 462,933 individuals of European ancestry, including 455.386 with OP. MR-PRESSO analysis detected no outlier SNPs.

First, LASSO regression (glmnet package), logistic regression (method = "glm"), and random forest (randomForest package) were applied to the training set data for feature selection to identify core predictors of osteoporosis. Subsequently, selected variables were integrated to construct a visual nomogram (rms:nomogram) for individualized risk quantification. Further, model discriminative ability was evaluated using receiver operating characteristic (ROC) curves (pROC package calculating AUC and 95% CI), while calibration curve analysis (rms:calibrate) supplemented by the Hosmer-Lemeshow goodness-of-fit test was performed on both training and validation sets to assess agreement between predicted risk and observed frequencies. Finally, decision curve analysis (DCA) (rmda:decision_curve) quantified the clinical net benefit across risk thresholds of 10–80%, comprehensively evaluating clinical utility. Analyses used R software (version 4.2.1). All NHANES analyses incorporated survey weights to account for the complex multi-stage probability sampling design. Categorical variables are presented as weighted frequencies (%), with group comparisons performed using weighted chi-square tests. Continuous variables are expressed as mean ± standard error (SE), and group comparisons used weighted linear regression models. Three multivariable logistic regression models assessed the association between NEE and OP: Crude model (unadjusted); Model 1 (adjusted for age, sex, race/ethnicity, education level, and marital status); Model 2 (further adjusted for smoking status, alcohol consumption, BMI, diabetes, hypertension, lumbar spine BMD, femoral neck BMD, and white blood cell count). Inverse-variance weighted (IVW) served as the primary method for Mendelian randomization (MR) analysis, with MR-Egger and weighted median as supplementary methods. MR findings are presented as odds ratios (ORs) per unit increase in the genetically predicted exposure. The MR-Egger intercept tested for horizontal pleiotropy (intercept close to 0 and p > 0.05 indicating absence). Heterogeneity among SNPs was assessed using Cochran's Q test (p > 0.05 and I² < 25% considered acceptable). MR reliability was evaluated via funnel plot symmetry and leave-one-out analysis. MR analyses used the "MRPRESSO" and "TwoSampleMR" R packages. p < 0.05 was considered statistically significant.

This study included 186 participants. The prevalence of osteoporosis in the overall cohort was 14.6%, with balanced distribution between the training and validation sets (14.5% vs. 14.8%). No statistically significant differences (p > 0.05) were observed between the training and validation sets regarding clinical characteristics including osteoporosis status, age, hypertension, diabetes, education level, height, BMI, etc. (Table 1).

Analysis of factors influencing osteoporosis based on the training set, using three methods (LASSO regression, logistic regression, random forest), all identified NEE as one of the influencing factors for osteoporosis. NEE showed a moderate effect in LASSO regression and random forest, while its importance score was lower in logistic regression. From the figure, it can be seen that the importance of NEE in model validation was higher than traditional social-behavioral factors like smoking and alcohol consumption in two models, but weaker than the gold standard BMD and age, indicating its value as an independent risk factor for supplementary screening, potentially providing early warning, especially in individuals with borderline BMD values. (Figure 2).

The influencing factors screened by the three models from the training set were used for nomogram and ROC curve analysis, constructing an osteoporosis prediction model. The nomogram model integrated key predictive variables, allowing for visual quantification of individual OP risk. Among them, NEE contributed significantly as an independent behavioral risk factor, forming a core predictive framework together with BMD indicators (lumbar/femoral). In ROC curve comparisons, the combined model including NEE (NEE + Femoral shaft density + Lumbar bone density + Age + Drinking) demonstrated high predictive performance (AUC = 0.808), significantly superior to single-variable models. This confirmed a synergistic effect between NEE and BMD indicators, and that NEE is an indispensable independent variable in the osteoporosis prediction model. (Figures 3A–C).

In the validation set, this prediction model maintained good discriminative ability, with an AUC of 0.898, confirming that variable synergy enhances model discrimination. DCA further confirmed that this NEE-containing model (red curve) yielded consistently higher net benefit than the "all-variables model" (green curve) and the "null model" (blue dashed line) within the 0.1–0.6 risk threshold range, with the most substantial net benefit advantage observed specifically within the clinically critical 0.2–0.4 threshold interval (Figures 4A,B). To evaluate the stability of the nomogram's predictive performance, Bootstrap resampling (1,000 times) showed the model's corrected C-index was 0.795 (95% CI: 0.752–0.838), indicating low model optimism and good internal stability (Supplementary Table 1).

Osteopenia was defined as a T-score between −1.0 and −2.5 standard deviations at either the lumbar spine or femoral neck. The worst-site principle was applied, whereby meeting the criterion at any site classified an individual as osteopenic. Three feature learning methods (LASSO regression, logistic regression, and random forest) were used to screen for influencing factors of osteopenia. The union of the top five influencing factors from these three methods was taken, finally resulting in seven important influencing factors. (Figure 5).

Subgroup analysis revealed that the effect of NEE is prominent in the osteopenic population. Full cohort analysis showed that for each 1-unit increase in NEE, osteoporosis risk increased by 20% (OR = 1.20, 95%CI:1.16–1.44, p < 0.05). In the osteopenic population, the effect size of NEE was significantly higher than in the normal bone mass population (OR = 1.12 vs. 1.06), and statistical significance jumped to the p < 0.001 level. This between-group difference was directly supported by interaction analysis (interaction p = 0.0297), indicating that osteopenic status may amplify the pathological effect of nighttime eating. At a deeper mechanistic level, in the osteopenic subgroup, NEE may enhance bone metabolic sensitivity leading to abnormally elevated nocturnal bone resorption peaks, disrupting the bone resorption balance. In the subgroup, the protective effect of BMD was significantly weakened (femoral density OR = 0.10 vs. full cohort 0.33). (Table 2).

Among the 15,464 participants included, there were 7,904 males and 7,560 females, with a mean age of (38.69 ± 22.23) years and a mean NEE of (19.96% ± 9.47%). Statistically significant differences (p < 0.05) between groups were observed for age, race, place of birth, education level, household income, smoking history, hypertension history, diabetes history, and other data. Basic characteristic comparisons are shown in Supplementary Table 2.

Weighted logistic regression analysis treating NEE intake level as a continuous variable showed a correlation between NEE intake and OP (OR = 1.67, p = 0.002). After adjusting for variables, Model 1 did not find a correlation, while Model 2 indicated a correlation between NEE and OP. After categorizing NEE and conducting correlation analysis with OP, the crude model and Model 2 showed its correlation. The higher the NEE value, the greater the risk of developing OP. Sensitivity analysis based on the fully-adjusted Model 2, with additional adjustment for total energy intake (Model 3), showed that the NEE–OP association was only slightly attenuated but remained significant (OR = 1.79, p = 0.002), suggesting that the effect of NEE is independent of total energy intake (Table 3). Sensitivity analyses for other thresholds in Model 2 are presented in Supplementary Table 3.

Subgroup analyses stratified by covariates were conducted using logistic regression to explore their relationships with osteoporosis (OP). Significant heterogeneity existed in the association between the exposure factor (NEE) and the outcome; the highest risk increases were observed in elderly individuals (≥60 years), those with low bone mineral density (BMD), and patients with sleep disorders (OR > 1.15). Age (p = 0.001), BMD (p ≤ 0.008), and sleep disorders (p = 0.001) were identified as significant effect modifiers (Figure 6).

A robust nonlinear association exists between NEE and OP: risk increased steeply within the NEE range of 10–25% (p ≤ 0.05), plateauing beyond 25%. This pattern remained highly significant (p ≤ 0.05) across the unadjusted model, demographically adjusted model, and fully adjusted model, confirming NEE as an independent predictor of bone health risk, distinct from age, sex, and metabolic factors. The 10–25% interval represents a critical threshold range for a sharp increase in risk, providing a quantitative basis for chrono-nutritional interventions (Figures 7A–C).

Multivariable Mendelian randomization analysis showed (Figure 8): 24-h dietary pattern (β = 0.0228, p = 0.0033) and BMI were significantly positively associated with osteoporosis risk; while sleep duration was negatively associated with risk. Sensitivity analysis found no horizontal pleiotropy (MR-Egger intercept test p > 0.05) or significant heterogeneity (Cochran's Q test p > 0.05), supporting the robustness of the results. Detailed results are in Supplementary Table 4.

The main strength of this study lies in its comprehensive research strategy: it integrates nationally representative observational data from NHANES, Mendelian randomization analysis for causal inference, and the construction and validation of diagnostic models in independent clinical cohorts, providing multi-level evidence. In the observational analysis, a wide range of covariates, including key bone mineral density indicators, were adjusted. It reveals nonlinear relationships and effect modifications: it clarifies the nonlinear dose–response relationship between NEE and OP and identifies high-risk populations (such as those with low bone mineral density, sleep disorders, etc.). It has certain clinical translatability, as it constructs and validates a practical clinical diagnostic prediction model for NEE.

Limitations include: The cross-sectional design of the clinical component precludes establishing the absolute temporal sequence and long-term effects of NEE on OP. Although Mendelian randomization (MR) was used to infer causality, the indirect exposure variables cannot fully equate to the direct causal effect of NEE. The definition and measurement of NEE (8:00 p.m. to 6:00 a.m.) were based on literature review and habits, lacking an internationally unified standard. Reliance on 24-h dietary recall introduces potential recall bias and measurement error. Specific food types, quality, and special circumstances like shift workers were not differentiated. The clinical cross-sectional study had a relatively small sample size (n = 186) from a single center (China), necessitating caution in extrapolation. While NHANES provided a large, nationally representative sample, its US population has cultural and dietary differences. Despite adjusting for multiple confounders in statistical models, residual confounding from unknown or unmeasured factors (e.g., detailed physical activity levels, specific nutrient intakes, other sleep parameters) may influence the results.