## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(
  echo = TRUE,
  rows.print = 5,
  message = FALSE,
  warning = FALSE,
  fig.width = 7, fig.height = 5)
set.seed(177813)

## ----requirement--------------------------------------------------------------
library(PLNmodels)
library(ggplot2)

## ----data_load----------------------------------------------------------------
data(microcosm)
most_abundant <- order(colSums(microcosm$Abundance), decreasing = TRUE)[1:30]
microcosm$Abundance <- microcosm$Abundance[, most_abundant]
mean(microcosm$Abundance == 0)

## ----zi_param, results = FALSE------------------------------------------------
myPLN     <- PLN(Abundance ~ 0 + site + offset(log(Offset)), data = microcosm)
zi_single <- ZIPLN(Abundance ~ 0 + site + offset(log(Offset)), data = microcosm)
zi_row    <- ZIPLN(Abundance ~ 0 + site + offset(log(Offset)), data = microcosm, zi = "row")
zi_col    <- ZIPLN(Abundance ~ 0 + site + offset(log(Offset)), data = microcosm, zi = "col")
zi_site   <- ZIPLN(Abundance ~ 0 + site + offset(log(Offset)) | 0 + site, data = microcosm)

## ----zi_param_table-----------------------------------------------------------
data.frame(
  model  = c("PLN", "ZIPLN (single)", "ZIPLN (row)", "ZIPLN (col)", "ZIPLN (site-dependent)"),
  loglik = c(myPLN$loglik, zi_single$loglik, zi_row$loglik, zi_col$loglik, zi_site$loglik),
  BIC    = c(myPLN$BIC, zi_single$BIC, zi_row$BIC, zi_col$BIC, zi_site$BIC),
  ICL    = c(myPLN$ICL, zi_single$ICL, zi_row$ICL, zi_col$ICL, zi_site$ICL)
) %>% knitr::kable(digits = 1)

## ----fitted, fig.cap = "Fitted vs. observed values"---------------------------
data.frame(
  fitted   = as.vector(fitted(zi_site)),
  observed = as.vector(microcosm$Abundance)
) %>%
  ggplot(aes(x = observed, y = fitted)) +
    geom_point(size = .5, alpha = .25) +
    scale_x_log10(limits = c(1, NA)) +
    scale_y_log10(limits = c(1, NA)) +
    theme_bw() + ggplot2::annotation_logticks()

## ----zi_by_site, fig.cap = "Estimated zero-inflation probability by body site, across the 30 species"----
one_obs_per_site <- !duplicated(microcosm$site)
pi_hat <- zi_site$model_par$Pi[one_obs_per_site, ]
rownames(pi_hat) <- as.character(microcosm$site[one_obs_per_site])

data.frame(site = rep(rownames(pi_hat), ncol(pi_hat)), zi_prob = as.vector(pi_hat)) %>%
  ggplot(aes(x = site, y = zi_prob)) +
    geom_boxplot() + theme_bw() +
    labs(x = "Body site", y = "Zero-inflation probability (per species)")

## ----zi_coef_B0, fig.cap = "Estimated regression coefficients of the zero-inflation component ($B_0$)"----
pheatmap::pheatmap(coefficients(zi_site, "zero"), cluster_rows = FALSE)

## ----zi_coef_B, fig.cap = "Estimated regression coefficients of the count component ($B$)"----
pheatmap::pheatmap(coefficients(zi_site, "count"), cluster_rows = FALSE)

## ----ziplnnetwork, results = FALSE--------------------------------------------
zi_models <- ZIPLNnetwork(Abundance ~ site + offset(log(Offset)), data = microcosm, control = ZIPLNnetwork_param(min_ratio = 0.01))

## ----ziplnnetwork_diagnostic, fig.cap = "Diagnostic of the ZIPLNnetwork fits"----
plot(zi_models, "diagnostic")

## ----ziplnnetwork_criteria, fig.cap = "Evolution of model selection criteria"----
plot(zi_models)

## ----ziplnnetwork_best, fig.cap = "Sparse residual network between the 30 most abundant taxa"----
zi_net <- getBestModel(zi_models, "EBIC")
plot(zi_net)