Simultaneous model for the joint and marginal distributions of multivariate product-multinomial categorical responses