  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 ~250 pc scale compared to a Poisson process.
This clustering mainly occurs in regions where the galactocentric distance
>~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.
