diff --git a/lbm-per-taxa.R b/lbm-per-taxa.R
new file mode 100644
index 0000000..09af5c7
--- /dev/null
+++ b/lbm-per-taxa.R
@@ -0,0 +1,35 @@
+library(microbiome)
+library(sbm)
+
+source("load-europe.R")
+
+base_otu <- cpd_phyloseq_europe
+
+taxoEstimateBipartiteSBM <- function(phyloseq_obj, taxonomic_level = "kingdom") {
+    taxolevel_otu <- microbiome::aggregate_taxa(phyloseq_obj, taxonomic_level)
+    taxolevel_mat <- phyloseq::otu_table(taxolevel_otu) |> as.matrix()
+
+    taxolevel_lbm <- estimateBipartiteSBM(netMat = taxolevel_mat, model = "poisson", dimLabels = c(row = taxonomic_level, col = "sample"))
+
+    rownames(taxolevel_lbm$probMemberships[[1]]) <- rownames(taxolevel_mat)
+    rownames(taxolevel_lbm$probMemberships[[2]]) <- colnames(taxolevel_mat)
+
+    return(taxolevel_lbm)
+}
+
+kingdom_lbm <- taxoEstimateBipartiteSBM(
+    phyloseq_obj = base_otu,
+    taxonomic_level = "kingdom"
+)
+
+kingdom_lbm$memberships
+
+phylum_lbm <- taxoEstimateBipartiteSBM(
+    phyloseq_obj = base_otu,
+    taxonomic_level = "phylum"
+)
+
+phylum_lbm$memberships
+
+mydf <- otu_table(aggregate_taxa(cpd_phyloseq_europe, "phylum"))
+library(ggplot2)