poisson.Rmd

---
title: "Poisson Log-Normal Mixed Model (Simplified for Binder)"
author:
  name: Christof Seiler
  affiliation: Department of Statistics, Stanford University
output:
  BiocStyle::html_document:
    toc_float: true
bibliography: bibliography.bib
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE)
```

# Goal

Reanalysis of mass cytometry data from @aghaeepour2017immune using the Poisson Log-Normal Mixed Model.

# Prerequisites

Parse input parameters.

```{r parse_input_parameters}
ncells = Inf
zenodo_url = "https://zenodo.org/record/2652578/files/"
prefit = paste0("cytoeffect_plmm_ncells_",ncells,".Rdata")
prefit
```

Load packages.

```{r load_packages}
library("cytoeffect")
library("tidyverse")
library("magrittr")
library("ggthemes")
library("cowplot")
library("ggcorrplot")
theme_set(theme_few())
```

# Load Data

Download preprocessed data from Zenodo.

```{r download_zenodo}
rdata_filenames = prefit
for(filename in rdata_filenames) {
  if(!file.exists(filename)) {
    download.file(url = paste0(zenodo_url, filename), 
                  destfile = filename,
                  mode = "wb")
  }
}
```

## HMC Diagnostics

Postprocessing of posterior samples. Traceplot of posterior samples.

```{r post_sampling, fig.wide=TRUE}
load(file = prefit)
pars_str = "beta"
rstan::traceplot(obj$fit_mcmc, inc_warmup = FALSE, pars = pars_str)
```

Some more MCMC diagnostics. According to empirically findings, Rhat > 1.1 is usually indicative of problems in the fit.

```{r mcmc_diagnostics}
pars = c("beta",
         "sigma","sigma_term","sigma_donor",
         "Cor","Cor_term","Cor_donor")
tb = rstan::summary(obj$fit_mcmc, 
                    pars = pars)$summary %>% 
  as.tibble(rownames = "pars", .before = 1) %>% 
  dplyr::select(pars, n_eff, Rhat)
tb %<>% na.omit() # Stan fills upper triangle with zeros
tb %>% arrange(n_eff)
tb %>% arrange(desc(Rhat))
tb %>% summarize(min = min(n_eff), max = max(n_eff))
tb %>% summarize(min = min(Rhat), max = max(Rhat))
```

# Results

Plot posterior regression coefficients.

```{r plot_beta}
p1 = plot(obj, type = "beta") + 
  ggtitle(expression("Fixed Effects"~beta)) +
  theme(legend.position = "bottom") +
  guides(col = guide_legend(ncol = 1)) + 
  scale_color_manual(values=c("#5DA5DA", "#FAA43A", "#60BD68"))
p1
plot(obj, type = "beta") + 
  facet_wrap(~condition, scales = "free_x") +
  theme(legend.position = "bottom") +
  guides(col = guide_legend(ncol = 1)) + 
  scale_color_manual(values=c("#5DA5DA", "#FAA43A", "#60BD68"))
```

Extract expected count difference for pSTAT1.

```{r fixed_effects_pSTAT1}
post_beta = rstan::extract(obj$fit_mcmc, pars = "beta")[[1]]
first_index = which(levels(pull(obj$df_samples_subset, obj$condition)) 
                    == "1st trimester")
third_index = which(levels(pull(obj$df_samples_subset, obj$condition)) 
                    == "3rd trimester")
pstat1_index = which(obj$protein_names == "pSTAT1")
first_log_count = quantile(post_beta[,pstat1_index,first_index], 
                           probs = c(0.025, 0.5, 0.975))
first_log_count
exp(first_log_count)
third_log_count = quantile(post_beta[,pstat1_index,third_index], 
                           probs = c(0.025, 0.5, 0.975))
third_log_count
exp(third_log_count)
diff_log_count = quantile(
  post_beta[,pstat1_index,third_index] - post_beta[,pstat1_index,first_index], 
  probs = c(0.025, 0.5, 0.975))
diff_log_count
exp(diff_log_count)
```

Posterior multivariate pairs plot.

```{r posterior_pair_plot}
pSTAT1_index = which(obj$protein_names == "pSTAT1")
pSTAT3_index = which(obj$protein_names == "pSTAT3")
pSTAT5_index = which(obj$protein_names == "pSTAT5")
post_beta = rstan::extract(obj$fit_mcmc, pars = "beta")[[1]]
tb_log_count = bind_rows(
  tibble(
    term = levels(pull(obj$df_samples_subset, obj$condition))[1],
    pSTAT1 = post_beta[,pSTAT1_index,1],
    pSTAT3 = post_beta[,pSTAT3_index,1],
    pSTAT5 = post_beta[,pSTAT5_index,1]
  ),
  tibble(
    term = levels(pull(obj$df_samples_subset, obj$condition))[2],
    pSTAT1 = post_beta[,pSTAT1_index,2],
    pSTAT3 = post_beta[,pSTAT3_index,2],
    pSTAT5 = post_beta[,pSTAT5_index,2]
  )
)
plot_diag = function(marker) {
  ggplot(tb_log_count, aes_string(marker, fill = "term")) + 
    geom_histogram(bins = 40, position = "identity", alpha = 0.5) +
    scale_fill_manual(values=c("#5DA5DA", "#FAA43A"))
}
plot_off_diag = function(marker1, marker2) {
  ggplot(tb_log_count, aes_string(marker1, marker2, color = "term")) +
    geom_density2d() + 
    scale_color_manual(values=c("#5DA5DA", "#FAA43A"))
}
ppair = plot_grid(
  plot_diag("pSTAT1") + theme(legend.position = "none"),
  NULL, 
  NULL,
  plot_off_diag("pSTAT1","pSTAT3") + theme(legend.position = "none"), 
  plot_diag("pSTAT3") + theme(legend.position = "none"), 
  NULL,
  plot_off_diag("pSTAT1","pSTAT5") + theme(legend.position = "none"), 
  plot_off_diag("pSTAT3","pSTAT5") + theme(legend.position = "none"), 
  plot_diag("pSTAT5") + theme(legend.position = "none"),
  ncol = 3
)
plot_grid(ppair,
          get_legend(plot_diag("pSTAT1") + theme(legend.position = "bottom")),
          ncol = 1,
          rel_heights = c(1, .1))
ggsave(filename = "posterior_multivariate_plmm.pdf", width = 8, height = 6)
```

Plot posterior standard deviation.

```{r posterior_sigma}
p2 = plot(obj, type = "sigma") + 
  ggtitle("Marker Standard Deviation"~sigma) +
  theme(legend.position = "bottom") +
  guides(col = guide_legend(ncol = 1)) +
  scale_color_manual(values=c("#5DA5DA", "#FAA43A", "#F17CB0"))
p2
```

Plot posterior correlations.

```{r posterior_cor}
plist = plot(obj, type = "Cor")
plist
```

Pairwise correlation change between conditions.

```{r correlation_uncertainty}
marker_pair = c("pSTAT3","pSTAT5")
Cor = rstan::extract(obj$fit_mcmc, pars = "Cor")[[1]]
Cor_term = rstan::extract(obj$fit_mcmc, pars = "Cor_term")[[1]]
Cor_diff = Cor_term - Cor
tb_cor = Cor_diff[,
                  which(obj$protein_names == marker_pair[1]),
                  which(obj$protein_names == marker_pair[2])] %>% as.tibble
tb_cor %<>% mutate(
  side = if_else(tb_cor$value > 0, 
                 true = paste0("positive (", 100*mean(tb_cor$value > 0), "%)"),
                 false = paste0("negative (", 100*mean(tb_cor$value <= 0), "%)"))
)
# keep colors consistent
if(mean(tb_cor$value > 0) == 1) {
  fill_colors = "#E46726"
} else {
  fill_colors = c("#6D9EC1","#E46726")
}
ggplot(tb_cor, aes(value, fill = side)) + 
  geom_histogram(bins = 50, alpha = 0.7) +
  xlab(paste0("Cor_term(", paste(marker_pair, collapse = ", "),")" )) +
  ggtitle("Posterior Distribution") + 
  scale_fill_manual(values = fill_colors)
```

Check if overall correlation structure changes between conditions.

```{r compare_covariance}
value = sapply(1:nrow(Cor_diff), function(i) {
  mask = which(upper.tri(Cor_diff[i,,]), arr.ind = T)
  cord = Cor_diff[i,,]
  mean(cord[lower.tri(cord)] > 0)
})
tb_cor = tibble(value = value)
tb_cor %<>% mutate(
  side = if_else(tb_cor$value > 0.5, 
                 true = paste0("> 1/2 (", 100*mean(tb_cor$value > 0.5), "%)"),
                 false = paste0("<= 1/2 (", 100*mean(tb_cor$value <= 0.5), "%)"))
)
p_global = ggplot(tb_cor, aes(value, fill = side)) + 
  geom_histogram(bins = 25, alpha = 0.7) +
  ggtitle(expression("Overall P(Corr"~Omega~"(3rd) > Corr"~Omega~"(1st))")) +
  scale_fill_manual(values = fill_colors) +
  theme(legend.position = "bottom") +
  xlab("probability")
p_global
```

Plot differential correlations.

```{r plot_differential_cor}
cor_increase = apply(X = Cor_diff, MARGIN = c(2,3), FUN = function(x) mean(x > 0))
colnames(cor_increase) = rownames(cor_increase) = obj$protein_names
p_local = ggcorrplot(cor_increase, hc.order = TRUE, type = "lower",
           outline.col = "lightgray",
           colors = c("#6D9EC1", "white", "#E46726")) +
  ggtitle(expression("P(Corr"~Omega~"(3rd) > Corr"~Omega~"(1st))")) +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank()) +
  scale_fill_gradient2(limit = c(0, 1), midpoint = 0.5,
                       low = "#6D9EC1", mid =  "white", high = "#E46726",
                       name = "probability")
p_local
```

Combine plot for paper.

```{r compbine_plot}
pall = plot_grid(
  p1, p2, 
  plist[[1]] + ggtitle(expression("Marker Corr"~Omega~"(1st trimester)")),
  plist[[2]] + ggtitle(expression("Marker Corr"~Omega~"(3rd trimester)")),
  p_global, p_local, 
  rel_heights = c(0.38,0.31,0.31),
  nrow = 3, labels = "AUTO"
)
ggsave(plot = pall, 
       filename = "posterior_summary_plmm.pdf", 
       width = 8, height = 11)
```

# Session Info {.unnumbered}

```{r session_info}
sessionInfo()
```

# References {.unnumbered}