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 $1sigma$ 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.