Lots of Seth fixes, change paper name, change repo name #55

Author: athowes

README.md

@@ -1,14 +1,14 @@
-# elgm-inf
+# naomi-aghq
 
-Code for the manuscript Howes, Stringer, Flaxman and Eaton "Fast approximate Bayesian inference of HIV indicators using the Naomi small-area estimation model" (in preparation).
+Code for the manuscript Howes, Stringer, Flaxman and Eaton "Fast approximate Bayesian inference of HIV indicators using PCA adaptive Gauss-Hermite quadrature" (in preparation).
 
 [Naomi](https://github.com/mrc-ide/naomi) ([Eaton et al, 2021](https://onlinelibrary.wiley.com/doi/10.1002/jia2.25788)) is a spatial evidence synthesis model used to produce district-level HIV epidemic indicators in sub-Saharan Africa.
 Multiple outcomes of interest, including HIV prevalence, HIV incidence and treatment coverage are jointly modelled using both household survey data and routinely reported health system data.
 The model is provided as a [tool](https://naomi.unaids.org/) for countries to input their data to and generate estimates during a yearly process supported by UNAIDS.
 Currently, inference is conducted using empirical Bayes and a Gaussian approximation via the [`TMB`](https://kaskr.github.io/adcomp/_book/Introduction.html) R package.
-We propose a new inference method which extends adaptive Gauss-Hermite quadrature to deal with >20 hyperparameters, enabling fast and accurate inference for Naomi and other [extended latent Gaussian](https://www.tandfonline.com/doi/full/10.1080/10618600.2022.2099403) models.
+We propose a new inference method extending adaptive Gauss-Hermite quadrature to deal with >20 hyperparameters, enabling fast and accurate inference for Naomi and other [extended latent Gaussian](https://www.tandfonline.com/doi/full/10.1080/10618600.2022.2099403) models.
 Using data from Malawi, our method improves the accuracy of inferences across a range of metrics, while being substantially faster to run than Hamiltonian Monte Carlo with the No-U-Turn sampler.
-By extending the [`aghq`](https://github.com/awstringer1/aghq) package ([Stringer, 2021](https://arxiv.org/abs/2101.04468)) we facilitate easy, flexible use of our method when provided a [`TMB`](https://kaskr.github.io/adcomp/_book/Introduction.html) C++ template for the model's log-posterior.
+Our implementation uses the [`aghq`](https://github.com/awstringer1/aghq) package ([Stringer, 2021](https://arxiv.org/abs/2101.04468)) facilitating easy, flexible use of the method when provided a [`TMB`](https://kaskr.github.io/adcomp/_book/Introduction.html) C++ template for the model's log-posterior.
 
 ![Example district-level Naomi model outputs for adults aged 15-49.](naomi_results.png)
 

elgm-inf.Rproj renamed to naomi-aghq.Rproj

@@ -1,5 +1,5 @@
 ---
-title: Appendix to "Fast approximate Bayesian inference of HIV indicators using the Naomi small-area estimation model"
+title: Appendix to "Fast approximate Bayesian inference of HIV indicators using principal components analysis adaptive Gauss-Hermite quadrature"
 author:
   - Adam Howes^[Department of Mathematics, Imperial College London]
   - Alex Stringer^[Department of Statistics and Actuarial Science, University of Waterloo]

src/docs_paper/appendix.Rmd

@@ -671,3 +671,17 @@ @article{karatzoglou2019package
   year={2019}
 }
 
+  @Article{brooks2017glmmTMB,
+    author = {Mollie E. Brooks and Kasper Kristensen and Koen J. {van
+      Benthem} and Arni Magnusson and Casper W. Berg and Anders Nielsen
+      and Hans J. Skaug and Martin Maechler and Benjamin M. Bolker},
+    title = {{glmmTMB} Balances Speed and Flexibility Among Packages
+      for Zero-inflated Generalized Linear Mixed Modeling},
+    year = {2017},
+    journal = {The R Journal},
+    doi = {10.32614/RJ-2017-066},
+    pages = {378--400},
+    volume = {9},
+    number = {2},
+  }
+

src/docs_paper/citations.bib

@@ -93,7 +93,9 @@ mvQuad::rescale(gg5, m = mu, C = cov, dec.type = 1)
 figA5 <- add_points(figA0, gg5) +
   labs(size = "")
 
-figA <- (figA1 + figA2) / (figA3 + figA4 + figA5)
+figA <- (figA1 + figA2) / (figA3 + figA4 + figA5) +
+  plot_annotation(tag_levels = "A") &
+  theme(plot.tag.position  = c(0.15, 0.95))
 
 #' Fig B
 

src/docs_paper/figures.R

@@ -60,7 +60,7 @@ depends:
       id: latest
       use:
         depends/posterior-contraction.png: posterior-contraction.png
-        depends/mean-sd.png: mean-sd.png
+        depends/mean-sd.png: mean-sd-alt.png
         depends/mean-sd.csv: mean-sd.csv
   - check_hyper-marginals:
       id: latest

src/docs_paper/orderly.yml

@@ -58,10 +58,18 @@ posterior_contraction_plot <- bind_rows(df_hyper, df_latent) %>%
   coord_flip() +
   labs(x = "", y = "Posterior contraction", col = "Type", shape = "Type") +
   scale_color_manual(values = c("#56B4E9", "#009E73")) +
-  scale_y_continuous(limits = c(-0.3, 1), breaks = c(-0.2, 0, 0.2, 0.4, 0.6, 0.8, 1)) +
+  scale_y_continuous(
+    limits = c(-1, 1),
+    breaks = c(-1.0, -0.8, -0.6, -0.4, -0.2, 0, 0.2, 0.4, 0.6, 0.8, 1.0)
+  ) +
+  geom_hline(yintercept = 0, linetype = "dashed", size = 0.25) +
+  annotate("text", x = 36, y = -0.6, size = 3, label = "Prior tighter") +
+  annotate("text", x = 36, y = 0.4, size = 3, label = "Posterior tighter") +
+  annotate("segment", x = 35, xend = 35, y = 0, yend = -1.0, size = 0.25, arrow = arrow(length = unit(0.2, "cm"))) +
+  annotate("segment", x = 35, xend = 35, y = 0, yend = 1.0, size = 0.25, arrow = arrow(length = unit(0.2, "cm"))) +
   theme_minimal()
 
-ggsave("posterior-contraction.png", posterior_contraction_plot, h = 6.5, w = 6.25)
+ggsave("posterior-contraction.png", posterior_contraction_plot, h = 6.5, w = 6.25, bg = "white")
 
 #' Which have lower amounts of posterior contraction?
 names(subset(posterior_contraction, posterior_contraction < 0.5))
@@ -126,47 +134,68 @@ df_metrics <- df_metrics %>%
 
 write_csv(df_metrics, "mean-sd.csv")
 
-mean_sd_plot <- ggplot(df_plot, aes(x = truth, y = approximate)) +
+mean_sd_plot <- ggplot(df_plot, aes(x = truth, y = approximate - truth)) +
   geom_point(shape = 1, alpha = 0.4) +
   facet_grid(indicator ~ method) +
-  coord_fixed(ratio = 1) +
-  geom_abline(slope = 1, intercept = 0, linetype = "dashed") +
+  geom_abline(slope = 0, intercept = 0, linetype = "dashed") +
   geom_text(data = df_metrics, aes(x = -Inf, y = Inf, label = label), size = 3, hjust = 0, vjust = 1.5) +
-  labs(x = "NUTS", y = "") +
+  labs(x = "NUTS", y = "Approximation - NUTS") +
   theme_minimal()
 
 ggsave("mean-sd.png", mean_sd_plot, h = 6, w = 6.25)
 
 #' Split into two plots for presentations etc.
 
+jitter_amount <- 0.02
+
 mean_plot <- df_plot %>%
   filter(indicator == "Posterior mean estimate") %>%
-  ggplot(aes(x = truth, y = approximate)) +
-  geom_point(shape = 1, alpha = 0.4) +
+  ggplot(aes(x = truth, y = approximate - truth)) +
+  geom_jitter(shape = 1, alpha = 0.4, width = jitter_amount, height = jitter_amount) +
+  lims(y = c(-0.4, 0.4)) +
   facet_grid(indicator ~ method) +
-  coord_fixed(ratio = 1) +
-  geom_abline(slope = 1, intercept = 0, linetype = "dashed") +
+  geom_abline(slope = 0, intercept = 0, linetype = "dashed", size = 0.25) +
   geom_text(
     data = filter(df_metrics, indicator == "Posterior mean estimate"),
     aes(x = -Inf, y = Inf, label = label), size = 3, hjust = 0, vjust = 1.5
   ) +
-  labs(x = "NUTS", y = "") +
+  labs(x = "NUTS", y = "Approximation - NUTS") +
   theme_minimal()
 
 ggsave("mean.png", mean_plot, h = 4, w = 6.25)
 
 sd_plot <- df_plot %>%
   filter(indicator == "Posterior SD estimate") %>%
-  ggplot(aes(x = truth, y = approximate)) +
-  geom_point(shape = 1, alpha = 0.4) +
+  ggplot(aes(x = truth, y = approximate - truth)) +
+  geom_jitter(shape = 1, alpha = 0.4, width = jitter_amount, height = jitter_amount) +
+  lims(y = c(-0.6, 0.6)) +
   facet_grid(indicator ~ method) +
-  coord_fixed(ratio = 1) +
-  geom_abline(slope = 1, intercept = 0, linetype = "dashed") +
+  geom_abline(slope = 0, intercept = 0, linetype = "dashed", size = 0.25) +
   geom_text(
     data = filter(df_metrics, indicator == "Posterior SD estimate"),
     aes(x = -Inf, y = Inf, label = label), size = 3, hjust = 0, vjust = 1.5
   ) +
-  labs(x = "NUTS", y = "") +
+  labs(x = "NUTS", y = "Approximation - NUTS") +
   theme_minimal()
 
 ggsave("sd.png", sd_plot, h = 4, w = 6.25)
+
+sd_plot_alt <- sd_plot +
+  theme(
+    strip.background = element_blank(),
+    strip.text.x = element_blank()
+  )
+
+y_axis <- ggplot(data.frame(l = mean_plot$labels$y, x = 1, y = 1)) +
+  geom_text(aes(x, y, label = l), angle = 90) +
+  theme_void() +
+  coord_cartesian(clip = "off")
+
+mean_plot_alt <- mean_plot
+mean_plot_alt$labels$x <- ""
+mean_plot_alt$labels$y <- sd_plot_alt$labels$y <- ""
+
+mean_sd_plot_alt <- y_axis + (mean_plot_alt / sd_plot_alt) +
+  plot_layout(widths = c(1, 30))
+
+ggsave("mean-sd-alt.png", mean_sd_plot_alt, h = 6, w = 6.25)