  We theoretically study topological phase transitions in four generalized
versions of the Kane-Mele-Hubbard model with up to $2\times 18^2$ sites. All
models are free of the fermion-sign problem allowing numerically exact quantum
Monte Carlo (QMC) calculations to be performed to extremely low temperatures.
We numerically compute the $\mathbb{Z}_2$ invariant and spin Chern number
$C_\sigma$ directly from the zero-frequency single-particle Green's functions,
and study the topological phase transitions driven by the tight-binding
parameters at different on-site interaction strengths. The $\mathbb{Z}_2$
invariant and spin Chern number, which are complementary to each another,
characterize the topological phases and identify the critical points of
topological phase transitions. Although the numerically determined phase
boundaries are nearly identical for different system sizes, we find strong
system-size dependence of the spin Chern number, where quantized values are
only expected upon approaching the thermodynamic limit. For the Hubbard models
we considered, the QMC results show that correlation effects lead to shifts in
the phase boundaries relative to those in the non-interacting limit, without
any spontaneously symmetry breaking. The interaction-induced shift is
non-perturbative in the interactions and cannot be captured within a "simple"
self-consistent calculation either, such as Hartree-Fock. Furthermore, our QMC
calculations suggest that quantum fluctuations from interactions stabilize
topological phases in systems where the one-body terms preserve the $D_3$
symmetry of the lattice, and destabilize topological phases when the one-body
terms break the $D_3$ symmetry.
