In a multispecies scenario without interspecific interactions, this could be done in ubms r-package, right? Simply stacking the data for all species into a single matrix, where each row represents a combination of site and species, and each column a sampling event (visit). Then, in the model formula, species identities would be included as random effects: ~(1|sp). What do you think?
If I'm correct, even N-mixture multispecies models could be fitted this way in ubms.