This work targets the identification of a class of models for hybrid dynamical systems characterized by nonlinear autoregressive exogenous (NARX) components, with finite-dimensional polynomial expansions, and by a Markovian switching mechanism. The estimation of the model parameters is performed under a probabilistic framework via Expectation Maximization, including submodel coefficients, hidden state values and transition probabilities. Discrete mode classification and NARX regression tasks are disentangled within the iterations. Soft-labels are assigned to latent states on the trajectories by averaging over the state posteriors and updated using the parametrization obtained from the previous maximization phase. Then, NARXs parameters are repeatedly fitted by solving weighted regression subproblems through a cyclical coordinate descent approach with coordinate-wise minimization. Moreover, we investigate a two stage selection scheme, based on a l1-norm bridge estimation followed by hard-thresholding, to achieve parsimonious models through selection of the polynomial expansion. The proposed approach is demonstrated on a SMNARX problem composed by three nonlinear sub-models with specific regressors.