Debonding of orthodontic brackets is an interface-dominated failure in the tooth–adhesive–bracket system, driven by mixed-mode peeling and shear. We develop a stochastic analytical–numerical framework (high-order adhesive model, mixed-mode cohesive interfaces, and in-house finite elements) with stochastic perturbation and pseudo arc-length continuation to quantify how uncertainty affects post-peak response and instability. Simulations show a clear sensitivity hierarchy: 10% uncertainty in adhesive elastic modulus has only a minor effect on peak debonding load and does not change the failure mode, whereas 10% uncertainty in adhesive thickness produces ~10% scatter in peak load and can switch the response from stable separation to unstable debonding. The framework enables probabilistic identification of dominant uncertainty sources and supports strategies to reduce unintended bracket failure.