diff --git a/analysis/assess_fits_68k_pbmc.Rmd b/analysis/assess_fits_68k_pbmc.Rmd
index 0d6ad10..9363cbe 100644
--- a/analysis/assess_fits_68k_pbmc.Rmd
+++ b/analysis/assess_fits_68k_pbmc.Rmd
@@ -23,6 +23,7 @@ library(ggplot2)
 library(cowplot)
 library(ggrepel)
 set.seed(1)
+source("../code/plot_functions.R")
 ```
 
 Load the 68k PBMC data, and the results of running `fit_poisson_nmf`
@@ -232,3 +233,25 @@ ggsave("../plots/loadings_scatterplot_pbmc68k_k12.png",
 ggsave("../plots/structure_plot_pbmc68k.eps",p3,height = 2,width = 8)
 ggsave("../plots/genes_scatterplot_pbmc68k.eps",p4,height = 4,width = 4.25)
 ```
+
+```{r elbo-plots, fig.height=6, fig.width=9}
+load("../output/pbmc68k/lda-pbmc68k.RData")
+p <- vector("list",12)
+for (i in 2:12) {
+  runs   <- which(dat$k == i)
+  p[[i]] <- create_elbo_plot(fits[runs],dat[runs,"runtime"],i)
+}
+plot_grid(p[[2]],p[[3]],p[[4]],p[[5]],
+          p[[6]],p[[7]],p[[8]],p[[9]],
+          p[[10]],p[[11]],p[[12]],
+          nrow = 3,ncol = 4)
+```
+
+```{r plots-for-paper-4, message=FALSE, echo=FALSE}
+ggsave("../plots/lda_elbo_pbmc68k.eps",
+       plot_grid(p[[2]],p[[3]],p[[4]],p[[5]],
+                 p[[6]],p[[7]],p[[8]],p[[9]],
+                 p[[10]],p[[11]],p[[12]],
+                 nrow = 3,ncol = 4),
+   	   height = 6,width = 9)
+```