  Traditional cosmological inference using Type Ia supernovae (SNeIa) have used
stretch- and color-corrected fits of SN Ia light curves and assumed a resulting
fiducial mean and symmetric intrinsic dispersion for the resulting relative
luminosity. As systematics become the main contributors to the error budget, it
has become imperative to expand supernova cosmology analyses to include a more
general likelihood to model systematics to remove biases with losses in
precision. To illustrate an example likelihood analysis, we use a simple model
of two populations with a relative luminosity shift, independent intrinsic
dispersions, and linear redshift evolution of the relative fraction of each
population. Treating observationally viable two-population mock data using a
one-population model results in an inferred dark energy equation of state
parameter $w$ that is biased by roughly 2 times its statistical error for a
sample of N $ \gtrsim$ 2500 SNeIa. Modeling the two-population data with a
two-population model removes this bias at a cost of an approximately $\sim20\%$
increase in the statistical constraint on $w$. These significant biases can be
realized even if the support for two underlying SNeIa populations, in the form
of model selection criteria, is inconclusive. With the current
observationally-estimated difference in the two proposed populations, a sample
of N $ \gtrsim$ 10,000 SNeIa is necessary to yield conclusive evidence of two
populations.
