  We investigate the possibility that GW170817 has not been the merger of two
conventional neutron stars (NS) but involved at least one if not two hybrid
stars with a quark matter core which might even belong to a third family of
compact stars. To this end, we develop a Bayesian analysis method for selecting
the most probable equation of state (EoS) under a set of constraints from
compact star physics, which now also include the tidal deformability from
GW170817 and the first result for the mass and radius determination for PSR
J0030+0451 by NICER. We apply this method for the first time to a two-parameter
family of hybrid EoS based on the DD2 model with nucleonic excluded volume for
hadronic matter and the color superconducting generalized nlNJL model for quark
matter. The model has a variable onset of deconfinement and can mimic the
effects of pasta phases with the possibility of a third family of hybrid stars
in the mass-radius ($M-R$) diagram. The main findings of this study are that:
1) the presence of multiple configurations for a given mass (twins) corresponds
to a set of disconnected lines in the diagram of tidal deformabilities for
binary mergers, so that merger events from the same mass range may result in a
probability landscape with different peak positions; 2) the Bayesian analysis
with the above observational constraints favors an early onset of the
deconfinement transition, at masses of $M_{\rm onset}\le 0.8~M_\odot$ with a
$M-R$ relationship that in the range of observed neutron star masses is almost
indistinguishable from that of a soft hadronic APR EoS; 3) a few yet fictitious
measurements of the NICER experiment with a $1\sigma$ range that is half of the
present value and different mass and radius would change the posterior
likelihood so that hybrid EoS with a phase transition onset in the range
$M_{\rm onset} = 1.1 - 1.6~M_\odot $ would be favored.
