Many longitudinal studies often require jointly modeling a biomarker and an event outcome, in order to provide more accurate inference and dynamic prediction of disease progression. Cystic fibrosis (CF) studies have illustrated the benefits of these models, primarily examining the joint evolution of lung‐function decline and survival. We propose a novel joint model within the shared‐parameter framework that accommodates nonlinear lung‐function trajectories, in order to provide more accurate inference on lung‐function decline over time and to examine the association between evolution of lung function and risk of a pulmonary exacerbation (PE) event recurrence. Specifically, a two‐level Gaussian process (GP) is used to estimate the nonlinear longitudinal trajectories and a flexible link function is introduced for a more accurate depiction of the binary process on the event outcome. Bayesian model assessment is used to evaluate each component of the joint model in simulation studies and an application to longitudinal data on patients receiving care from a CF center. A nonlinear structure is suggested by both longitudinal continuous and binary evaluations. Including a flexible link function improves model fit to these data. The proposed hierarchical GP model with a flexible power link function where Laplace distribution is the baseline (spep) has the best fit of all joint models considered, characterizing how accelerated lung‐function decline corresponds to increased odds of experiencing another PE.