We demonstrate the power of Gibbs point process models from the spatial statistics literature when applied to studies of resolved galaxies. We conduct a rigorous analysis of the spatial distributions of objects in the star formation complexes of M33, including giant molecular clouds (GMCs) and young stellar cluster candidates (YSCCs). We choose a hierarchical model structure from GMCs to YSCCs based on the natural formation hierarchy between them. This approach circumvents the limitations of the empirical two-point correlation function analysis by naturally accounting for the inhomogeneity present in the distribution of YSCCs. We also investigate the effects of GMCs’ properties on their spatial distributions. We confirm that the distribution of GMCs and YSCCs are highly correlated. We found that the spatial distributions of YSCCs reaches a peak of clustering pattern at around 250 pc scale compared to a Poisson process. This clustering mainly occurs in regions where the galactocentric distance is greater than 4.5 kpc. Furthermore, the galactocentric distance of GMCs and their mass have strong positive effects on the correlation strength between GMCs and YSCCs. We outline some possible implications of these findings for our understanding of the cluster formation process.