Supervized classification of multivariate count table with the Poisson discriminant Analysis

Covariates and offsets

Just like PLN, PLN-LDA generalizes to a formulation close to a multivariate generalized linear model where the main effect is due to a linear combination of the discrete group structure, \(d\) covariates \(\mathbf{x}_i\) and a vector \(\mathbf{o}_i\) of \(p\) offsets in sample \(i\). The latent layer then reads \[\begin{equation} \mathbf{Z}_i \sim \mathcal{N}({\mathbf{o}_i + \mathbf{g}_i^\top \mathbf{M} + \mathbf{x}_i^\top\mathbf{B}},\boldsymbol\Sigma) \end{equation}\] where \(\mathbf{B}\) is a \(d\times p\) matrix of regression parameters.

Prediction

Given:

• a new observation \(\mathbf{Y}\) with associated offset \(\mathbf{o}\) and covariates \(\mathbf{x}\)
• a model with estimated parameters \(\hat{\boldsymbol{\Sigma}}\), \(\hat{\mathbf{B}}\), \(\hat{\mathbf{M}}\) and group counts \((n_1, \dots, n_K)\)

We can predict the observation’s group using Bayes rule as follows: for \(k \in {1, \dots, K}\), compute \[\begin{equation} \begin{aligned} f_k(\mathbf{Y}) & = p(\mathbf{Y} | \mathbf{g} = k, \mathbf{o}, \mathbf{x}, \hat{\mathbf{B}}, \hat{\boldsymbol{\Sigma}}) \\ & = \boldsymbol{\Phi}_{PLN}(\mathbf{Y}; \mathbf{o} + \boldsymbol{\mu}_k + \mathbf{x}^\top \hat{\mathbf{B}}, \hat{\boldsymbol{\Sigma}}) \\ p_k & = \frac{n_k}{\sum_{k' = 1}^K n_{k'}} \end{aligned} \end{equation}\] where \(\boldsymbol{\Phi}_{PLN}(\bullet; \boldsymbol{\mu}, \boldsymbol{\Sigma})\) is the density function of a PLN distribution with parameters \((\boldsymbol{\mu}, \boldsymbol{\Sigma})\). \(f_k(\mathbf{Y})\) and \(p_k\) are respectively plug-in estimates of (i) the probability of observing counts \(\mathbf{Y}\) in a sample from group \(k\) and (ii) the probability that a sample originates from group \(k\).

The posterior probability \(\hat{\pi}_k(\mathbf{Y})\) that observation \(\mathbf{Y}\) belongs to group \(k\) and most likely group \(\hat{k}(\mathbf{Y})\) can thus be defined as \[\begin{equation} \begin{aligned} \hat{\pi}_k(\mathbf{Y}) & = \frac{p_k f_k(\mathbf{Y})}{\sum_{k' = 1}^K p_{k'} f_{k'}(\mathbf{Y})} \\ \hat{k}(\mathbf{Y}) & = \underset{k \in \{1, \dots, K\}}{\arg\max} \hat{\pi}_k(\mathbf{Y}) \end{aligned} \end{equation}\]

Optimization by Variational inference

Classification and prediction are the main objectives in (PLN-)LDA. To reach this goal, we first need to estimate the model parameters. Inference in PLN-LDA focuses on the group-specific main effects \(\mathbf{M}\), the regression parameters \(\mathbf{B}\) and the covariance matrix \(\boldsymbol\Sigma\). Technically speaking, we can treat \(\mathbf{g}_i\) as a discrete covariate and estimate \([\mathbf{M}, \mathbf{B}]\) using the same strategy as for the standard PLN model. Briefly, we adopt a variational strategy to approximate the log-likelihood function and optimize the consecutive variational surrogate of the log-likelihood with a gradient-ascent-based approach. To this end, we rely on the CCSA algorithm of Svanberg (2002) implemented in the C++ library (Johnson 2011), which we link to the package.

Structure of PLNLDAfit

The myLDA_nocov variable is an R6 object with class PLNLDAfit, which comes with a couple of methods. The most basic is the show/print method, which sends a brief summary of the estimation process and available methods:

myLDA_nocov

## Linear Discriminant Analysis for Poisson Lognormal distribution
## ==================================================================
##  nb_param   loglik       BIC       ICL
##       357 -799.824 -1494.514 -1082.492
## ==================================================================
## * Useful fields
##     $model_par, $latent, $latent_pos, $var_par, $optim_par
##     $loglik, $BIC, $ICL, $loglik_vec, $nb_param, $criteria
## * Useful S3 methods
##     print(), coef(), sigma(), vcov(), fitted()
##     predict(), predict_cond(), standard_error()
## * Additional fields for LDA
##     $percent_var, $corr_map, $scores, $group_means
## * Additional S3 methods for LDA
##     plot.PLNLDAfit(), predict.PLNLDAfit()

Comprehensive information about PLNLDAfit is available via ?PLNLDAfit.

Specific fields

The user can easily access several fields of the PLNLDAfit object using S3 methods, the most interesting ones are

• the \(p \times p\) covariance matrix \(\hat{\boldsymbol{\Sigma}}\):

sigma(myLDA_nocov) %>% corrplot::corrplot(is.corr = FALSE)

• the regression coefficient matrix \(\hat{\mathbf{B}}\) (in this case NULL as there are no covariates)

coef(myLDA_nocov)

## NULL

• the \(p \times K\) matrix of group means \(\mathbf{M}\)

myLDA_nocov$group_means %>% head() %>% knitr::kable(digits = 2)

The PLNLDAfit class also benefits from two important methods: plot and predict.

plot method

The plot methods provides easy to interpret graphics which reveals here that the groups are well separated:

plot(myLDA_nocov)

By default, plot shows the first 3 axis of the LDA when there are 4 or more groups and uses special representations for the edge cases of 3 or less groups.

ggplot2-savvy users who want to make their own representations can extracts the \(n \times (K-1)\) matrix of sample scores from the PLNLDAfit object …

myLDA_nocov$scores %>% head %>% knitr::kable(digits = 2)

…or the \(p \times (K-1)\) matrix of correlations between scores and (latent) variables

myLDA_nocov$corr_map %>% head %>% knitr::kable(digits = 2)

predict method

The predict method has a slightly different behavior than its siblings in other models of the PLNmodels. The goal of predict is to predict the discrete class based on observed species counts (rather than predicting counts from known covariates).

By default, the predict use the argument type = "posterior" to output the matrix of log-posterior probabilities \(\log(\hat{\pi})_k\)

predicted.class <- predict(myLDA_nocov, newdata = trichoptera)
## equivalent to 
## predicted.class <- predict(myLDA_nocov, newdata = trichoptera,  type = "posterior")
predicted.class %>% head() %>% knitr::kable(digits = 2)

You can also show them in the standard (and human-friendly) \([0, 1]\) scale with scale = "prob" to get the matrix \(\hat{\pi}_k\)

predicted.class <- predict(myLDA_nocov, newdata = trichoptera, scale = "prob")
predicted.class %>% head() %>% knitr::kable(digits = 3)

Setting type = "response", we can predict the most likely group \(\hat{k}\) instead:

predicted.class <- predict(myLDA_nocov, newdata = trichoptera,  type = "response")
predicted.class

##  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 
##  1  1  1  1  1  1  9  2  1  1  1  1  2  3  2  2  2  3  3  3  9  3  4  1  4  4 
## 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 
## 12  5  4  5  6  7  7  7  8  8  8  8  8  1  9  9  9 10 10 10 10 11 12 
## Levels: 1 2 3 4 5 6 7 8 9 10 11 12

We can assess that the predictions are quite similar to the real group (this is not a proper validation of the method as we used data set for both model fitting and prediction and are thus at risk of overfitting).

table(predicted.class, trichoptera$Group, dnn = c("predicted", "true"))

##          true
## predicted  1  2  3  4  5  6  7  8  9 10 11 12
##        1  10  0  0  1  0  0  0  0  1  0  0  0
##        2   1  4  0  0  0  0  0  0  0  0  0  0
##        3   0  1  4  0  0  0  0  0  0  0  0  0
##        4   0  0  0  3  1  0  0  0  0  0  0  0
##        5   0  0  0  0  2  0  0  0  0  0  0  0
##        6   0  0  0  0  0  1  0  0  0  0  0  0
##        7   0  0  0  0  0  0  3  0  0  0  0  0
##        8   0  0  0  0  0  0  0  4  1  0  0  0
##        9   1  0  1  0  0  0  0  0  3  0  0  0
##        10  0  0  0  0  0  0  0  0  0  4  0  0
##        11  0  0  0  0  0  0  0  0  0  0  1  0
##        12  0  0  0  0  1  0  0  0  0  0  0  1

Finally, we can get the coordinates of the new data on the same graph at the original ones with type = "scores". This is done by averaging the latent positions \(\hat{\mathbf{Z}}_i + \boldsymbol{\mu}_k\) (found when the sample is assumed to come from group \(k\)) and weighting them with the \(\hat{\pi}_k\). Some samples, have compositions that put them very far from their group mean.

library(ggplot2)
predicted.scores <- predict(myLDA_nocov, newdata = trichoptera,  type = "scores")
colnames(predicted.scores) <- paste0("Axis.", 1:ncol(predicted.scores))
predicted.scores <- as.data.frame(predicted.scores)
predicted.scores$group <- trichoptera$Group
plot(myLDA_nocov, map = "individual", nb_axes = 2, plot = FALSE) + 
  geom_point(data = predicted.scores, 
             aes(x = Axis.1, y = Axis.2, color = group, label = NULL))

Specific fields

All fields of our new PLNLDA fit can be accessed as before with similar results. The only important difference is the result of coef: since we included a covariate in the model, coef now returns a 1-column matrix for \(\hat{\mathbf{B}}\) instead of NULL

coef(myLDA_cov) %>% head %>% knitr::kable()

The group-specific main effects can still be accessed with $group_means

myLDA_cov$group_means %>% head %>% knitr::kable(digits = 2)

plot method

Once again, the plot method is very useful to get a quick look at the results.

plot(myLDA_cov)

predict method

We can again predict the most likely group for each sample :

predicted.class_cov <- predict(myLDA_cov, newdata = trichoptera, type = "response")

and check that we recover the correct class in most cases (again, we used the same data set for model fitting and group prediction only for ease of exposition):

table(predicted.class_cov, trichoptera$Group, dnn = c("predicted", "true"))

##          true
## predicted  1  2  3  4  5  6  7  8  9 10 11 12
##        1  11  0  0  1  0  0  0  0  1  0  0  0
##        2   1  4  0  0  0  0  0  0  0  0  0  0
##        3   0  1  4  0  0  0  0  0  0  0  0  0
##        4   0  0  0  3  1  0  0  0  0  0  0  0
##        5   0  0  0  0  2  0  0  0  0  0  0  0
##        6   0  0  0  0  0  1  0  0  0  0  0  0
##        7   0  0  0  0  0  0  3  0  0  0  0  0
##        8   0  0  0  0  0  0  0  4  1  0  0  0
##        9   0  0  1  0  0  0  0  0  3  0  0  0
##        10  0  0  0  0  0  0  0  0  0  4  0  0
##        11  0  0  0  0  0  0  0  0  0  0  1  0
##        12  0  0  0  0  1  0  0  0  0  0  0  1