diff --git a/R/psmelt.R b/R/psmelt.R
index 5953714..abfe77d 100644
--- a/R/psmelt.R
+++ b/R/psmelt.R
@@ -63,7 +63,6 @@
 #' p = ggplot(mdf, aes(x=SampleType, y=Abundance, fill=Genus))
 #' p = p + geom_bar(color="black", stat="identity", position="stack")
 #' print(p)
-# TODO: test behavior vs. phyloseq's
 psmelt = function(physeq){
   # Access covariate names from object, if present
   if(!inherits(physeq, "phyloseq")){
@@ -152,9 +151,9 @@ psmelt = function(physeq){
       tb <- tb %>%
           dplyr::left_join(tax, by = "OTU")
   }
-  # Arrange by Abundance; TODO: verify this sorting behavior matches phyloseq's
+  # Arrange by Abundance, then OTU names (to approx. phyloseq behavior)
   tb <- tb %>%
-      dplyr::arrange(desc(Abundance))
+      dplyr::arrange(desc(Abundance), OTU)
   # Return as a data.frame for phyloseq compatibility
   tb %>% as.data.frame
 }
diff --git a/tests/testthat/test-psmelt.R b/tests/testthat/test-psmelt.R
new file mode 100644
index 0000000..b09c63e
--- /dev/null
+++ b/tests/testthat/test-psmelt.R
@@ -0,0 +1,35 @@
+context("Equivalence of psmelt() with phyloseq::psmelt()")
+
+library(dplyr)
+
+# Test on a subset of GlobalPatterns
+data(GlobalPatterns)
+set.seed(20190421)
+ps <- GlobalPatterns %>%
+    {prune_taxa(sample(taxa_names(.), 200), .)} %>%
+    tax_glom("Genus")
+
+test_that("psmelt() functionally matches phyloseq::psmelt()", {
+    # In this example, phyloseq::psmelt() drops the "Species" column from the
+    # tax_table (which is full of NAs). But speedyseq::psmelt() does not
+    # discard columns even if they contain no non-missing data. Other
+    # differences: phyloseq's output has rownames (which seem not meaningful);
+    # speedyseq's output does not; both data frames are sorted by Abundance,
+    # but the row order differs in cases where Abundance and OTU names differ.
+    options(stringsAsFactors = TRUE)
+    tb1 <- phyloseq::psmelt(ps)
+    tb2 <- psmelt(ps)
+    expect_true(is.factor(tb1$Kingdom))
+    expect_true(all_equal(tb1, select(tb2, -Species),
+        ignore_col_order = FALSE, ignore_row_order = TRUE))
+    expect_equal(tb1$Abundance, tb2$Abundance)
+    expect_equal(tb1$OTU, tb2$OTU)
+    options(stringsAsFactors = FALSE)
+    tb1 <- phyloseq::psmelt(ps)
+    tb2 <- psmelt(ps)
+    expect_false(is.factor(tb1$Kingdom))
+    expect_true(all_equal(tb1, select(tb2, -Species),
+        ignore_col_order = FALSE, ignore_row_order = TRUE))
+    expect_equal(tb1$Abundance, tb2$Abundance)
+    expect_equal(tb1$OTU, tb2$OTU)
+})