<p>Landslides are among the most devastating natural hazards globally, and those triggered by earthquakes pose serious risks due to their sudden onset, and high concentration within affected zones. The landslide community mostly focuses on landslide susceptibility and size models, where these two approaches are typically treated independently, limiting their usefulness for comprehensive hazard assessments. Moreover, most regional landslide modeling frameworks lack integration between occurrence and size predictions and often disregard spatial dependencies that are characteristic of co-seismic landslides. To address this gap, we developed a unified, spatially dependent framework to jointly estimate the probability of landslide occurrence and the probability that a landslide exceeds a given size threshold within each slope unit. Here we show that integrating binomial and log-Gaussian generalized additive models (GAMs) with spatial smoothing significantly improves the accuracy and interpretability of co-seismic landslide hazard prediction. We applied this approach to a comprehensive landslide inventory from the 2022 Mw 6.8 Luding earthquake in Sichuan, China, and used a Markov Chain Monte Carlo (MCMC) method to estimate size exceedance probabilities. Results indicate that accounting for spatial dependence improves prediction performance for both susceptibility (AUC = 0.975) and size (MAE = 0.844, RMSE = 1.077, <i>r</i> = 0.781) prediction, and enables robust hazard mapping with uncertainty quantification. These findings highlight the value of a joint, spatially aware modeling approach in understanding and forecasting co-seismic landslide hazards, offering a scalable tool for post-earthquake risk assessment and future regional applications.</p>