Ultracold atoms provide clues to an important many-body problem regarding the
dependence of Bose-Einstein condensation (BEC) transition temperature $T_c$ on
interactions. However, cold atoms are trapped in harmonic potentials and
theoretical evaluations of the $T_c$ shift of trapped interacting Bose gases
are challenging. While previous predictions of the leading-order shift have
been confirmed, more recent experiments exhibit higher-order corrections beyond
available mean-field theories. By implementing two large-N based theories with
the local density approximation (LDA), we extract next-order corrections of the
$T_c$ shift. The leading-order large-N theory produces results quantitatively
different from the latest experimental data. The leading-order auxiliary field
(LOAF) theory containing both normal and anomalous density fields captures the
$T_c$ shift accurately in the weak interaction regime. However, the LOAF theory
shows incompatible behavior with the LDA and forcing the LDA leads to density
discontinuities in the trap profiles. We present a phenomenological model based
on the LOAF theory, which repairs the incompatibility and provides a prediction
of the $T_c$ shift in stronger interaction regime.