Population genetic data from multiple taxa can address comparative phylogeographic questions about community-scale response to environmental shifts, and a useful strategy to this end is to employ hierarchical co-demographic models that directly test multi-taxa hypotheses within a single, unified analysis while benefiting in statistical power from aggregating datasets. This approach has been applied to classical phylogeographic datasets such as mitochondrial barcodes as well as reduced-genome polymorphism datasets that can yield 10,000s of SNPs, produced by emergent technologies such as RAD-seq and GBS. A strategy for the latter had been accomplished by adapting the site frequency spectrum to a novel summarization of population genomic data across multiple taxa called the aggregate site frequency spectrum (aSFS), which potentially can be deployed under various inferential frameworks including approximate Bayesian computation, random forest, and composite likelihood optimization. Here, we introduce the R package Multi-DICE, a wrapper program that exploits existing simulation software for straight-forward and flexible execution of hierarchical model-based inference using the aSFS, which is derived from genomic-scale data, as well as mitochondrial data. We validate several novel software features such as applying alternative inferential frameworks, enforcing a minimal threshold of time surrounding event pulses, and specifying flexible hyperprior distributions. In sum, Multi-DICE provides comparative analysis within the familiar R environment while allowing a high degree of user customization, and will thus serve as a valuable tool for comparative phylogeography and population genomics.