Poisson PCA: Poisson Measurement Error corrected PCA, with Application to Microbiome Data

TL;DR

This paper introduces Poisson PCA, a semiparametric method for principal component analysis of Poisson-distributed data, effectively correcting bias and accounting for sequencing depth, with applications to microbiome data.

Contribution

We develop a novel semiparametric Poisson PCA approach that corrects bias in variance estimation and computes principal scores without relying on parametric models.

Findings

Our method outperforms PLN in identifying main principal components.

It is faster and more robust to outliers than existing parametric approaches.

Applications to microbiome data demonstrate practical effectiveness.

Abstract

In this paper, we study the problem of computing a Principal Component Analysis of data affected by Poisson noise. We assume samples are drawn from independent Poisson distributions. We want to estimate principle components of a fixed transformation of the latent Poisson means. Our motivating example is microbiome data, though the methods apply to many other situations. We develop a semiparametric approach to correct the bias of variance estimators, both for untransformed and transformed (with particular attention to log-transformation) Poisson means. Furthermore, we incorporate methods for correcting different exposure or sequencing depth in the data. In addition to identifying the principal components, we also address the non-trivial problem of computing the principal scores in this semiparametric framework. Most previous approaches tend to take a more parametric line. For example the Poisson-log-normal (PLN) model, approach. We compare our method with the PLN approach and find that our method is better at identifying the main principal components of the latent log-transformed Poisson means, and as a further major advantage, takes far less time to compute. Comparing methods on real data, we see that our method also appears to be more robust to outliers than the parametric method.

Figures

Equations

$$\mathrm{Var}(\mathbf{X})$$

$$=\mathbb{E}\left(\mathrm{diag}(\mathbf{\Lambda })\right)+\mathrm{Var}\left(\mathbf{\Lambda }\right)$$

$$\mathrm{Var}\left(\mathbf{\Lambda }\right)$$

$$=\mathrm{Var}(\mathbf{X})-\mathrm{diag}\left(\mathbb{E}(\mathbf{X})\right)$$

$$\widehat{\mathbb{E}(\mathbf{X})}$$

$$\widehat{\mathrm{Var}(\mathbf{X})}$$

$$\widehat{\mathrm{Var}\left(\mathbf{\Lambda }\right)}=\frac{1}{n-1}{X}^{T}X-\frac{1}{n(n-1)}{X}^{T}11^{T}X-\mathrm{diag}\left(\frac{1}{n}{X}^{T}1\right)$$

$$\mathrm{Var}(\mathbf{X}|S)$$

$$=\mathbb{E}\left(\mathrm{diag}(\mathbf{\Lambda }S)|S\right)+\mathrm{Var}\left(\mathbf{\Lambda }S|S\right)$$

$$=S\mathbb{E}\left(\mathrm{diag}(\mathbf{\Lambda })|S\right)+S^2\mathrm{Var}\left(\mathbf{\Lambda }|S\right)$$

$$\mathrm{Var}\left(\mathbf{\Lambda }\right)=\mathbb{E}\mathrm{Var}\left(\mathbf{\Lambda }|S\right)$$

$$=\mathbb{E}\mathrm{Var}\left(\frac{\mathbf{X}}{S}|S\right)-\mathbb{E}\mathrm{diag}\left(\mathbb{E}\left(\frac{\mathbf{X}}{S^2}\right)|S\right)$$

$$=\mathbb{E}\mathrm{Var}\left(\frac{\mathbf{X}}{S}|S\right)-\mathrm{diag}\left(\mathbb{E}\left(\frac{\mathbf{X}}{S^2}\right)\right)$$

$$\mathrm{Var}\left(\frac{\mathbf{X}}{S}\right)$$

$$=\mathbb{E}\left(\mathrm{Var}\left(\frac{\mathbf{X}}{S}|S\right)\right)+\mathrm{Var}\left(\frac{\mathbb{E}(\mathbf{\Lambda })S}{S}|S\right)$$

$$=\mathbb{E}\left(\mathrm{Var}\left(\frac{\mathbf{X}}{S}|S\right)\right)$$

$$\mathrm{Var}\left(\mathbf{\Lambda }\right)$$

$$\widehat{\mathrm{Var}\left(\frac{\mathbf{X}}{S}\right)}$$

$$\widehat{\mathbb{E}\left(\frac{\mathbf{X}}{S^2}\right)}$$

$$\widehat{\mathrm{Var}\left(\mathbf{\Lambda }\right)}$$

$$\mathrm{Var}(g(\mathbf{X}))=\mathbb{E}\left(\mathrm{Var}(g(\mathbf{X})|\Lambda )\right)+\mathrm{Var}\left(\mathbb{E}(g(\mathbf{X})|\Lambda )\right)$$

$$\mathrm{Var}(f(\mathbf{\Lambda }))=\mathrm{Var}(g(\mathbf{X}))-\mathbb{E}\left(\mathrm{Var}(g(\mathbf{X})|\mathbf{\Lambda })\right)$$

$$\mathbb{E}(g(\mathbf{X})|\mathbf{\Lambda }=\lambda )$$

$$e^{\lambda }f(\lambda )$$

$$e^{\lambda }{\lambda }^k$$

$$\mathrm{Var}(g(\mathbf{X})|\mathbf{\Lambda }=\lambda )$$

$$=e^{-\lambda }\sum _{n=0}^{\infty }\frac{{\lambda }^n}{n!}g(n{)}^2-e^{-2\lambda }{\left(\sum _{n=0}^{\infty }\frac{{\lambda }^n}{n!}g(n)\right)}^2$$

$$=\sum _{n=0}^{\infty }g(n{)}^2{e}^{-\lambda }\frac{{\lambda }^n}{n!}-e^{-2\lambda }\left(\sum _{n=0}^{\infty }\sum _{m=0}^{\infty }\frac{{\lambda }^n}{n!}\frac{{\lambda }^m}{m!}g(m)g(n)\right)$$

$$=\sum _{n=0}^{\infty }g(n{)}^2{e}^{-\lambda }\frac{{\lambda }^n}{n!}-\sum _{n=0}^{\infty }h(n)e^{-2\lambda }\frac{{\lambda }^n}{n!}$$

$$s_{n}(\mathbf{X})=\left\{\begin{array}[]{ll}1&\textrm{if }\mathbf{X}=n\\ 0&\textrm{otherwise}\end{array}\right.$$

$$t_n(\mathbf{X})=(-1{)}^{\mathbf{X}-n}(\genfrac{}{}{0ex}{}{\mathbf{X}}{n})$$

$$\widehat{\mathrm{Var}\left(f({\mathbf{\Lambda }}_{\mathbf{l}})\right)}$$

$$=\frac{1}{n-1}\sum _{i=1}^{n}g(X_{il}{)}^2-\frac{1}{n(n-1)}{\left(\sum _{i=1}^{n}g(X_{il})\right)}^2-\frac{1}{n}\sum _{i=1}^n\left(\sum _{k=0}^{\infty }g(k{)}^2{s}_k(X_{il})-\sum _{k=0}^{\infty }h(k)t_k(X_{il})\right)$$

$$=\frac{1}{n}\sum _{i=1}^n\sum _{k=0}^{X_{il}}(-1{)}^{X_{il}-k}(\genfrac{}{}{0ex}{}{X_{il}}{k})h(k)-\frac{1}{n(n-1)}\left(\sum _{\stackrel{i\ne j}{i,j=1}}^{n}g(X_{il})g(X_{jl})\right)$$

$$\log (\lambda )=\log (a)+\sum _{n=1}^{\infty }\frac{(-1{)}^n}{n}{\left(\frac{\lambda }{a}-1\right)}^n$$

$$g(n)=\left\{\begin{array}[]{ll}\log(a)+\sum_{i=1}^{N}\frac{1}{i}&\textrm{if }n=0\\ \log(a)+\sum_{i=1}^{N}\frac{1}{i}+\sum_{i=1}^{n\land N}\frac{(-1)^{i}n!}{a^{i}(n-i)!}\sum_{m=i}^{N}\frac{\binom{m}{i}}{m}&\textrm{if }1\leqslant n<l_{0}\\ \log(n)&\textrm{if }n\geqslant l_{0}\end{array}\right.$$

$$\tilde{\Sigma }=\Sigma -a{1}^T-1a^T$$

$$\left(\Sigma -a{1}^T-1a^T\right)1$$

$$\Sigma 1-a{1}^{T}1-1a^{T}1$$

Full text

Keywords: Dimension reduction; Compositional data; Sequencing depth correction; Count data; Log-normal Poisson

1 Introduction

Principal component analysis (PCA) is a widely used dimension reduction and data exploration technique. It is particularly suited for real valued multivariate data analysis with the underlying assumption that the measurement error in the data is additive. This can be seen from the fact that one way to define the PCA is by minimizing the reconstruction error of the data by the first $k$ principal components with the objective function defined as squared error loss. For other types of data, this loss function can be natually extended to be the negative loglikelihood function. This can be most natually developed to the whole exponential family distribution, for example, Collins et al. (2001) and Chiquet et al. (2018) maximized the exponential family loglikelihood of the data by assuming a lower dimensional representation of the natural parameters through the canonical link functions. Such an approach can be deemed as totally parametric and thus is sensitive to the model assumptions, outliers and the choice of the dimension for the latent parameter space. To some extent, this is more a parametric modelling method instead of data exploration.

Our focus in this paper is to develop a new PCA method for count data by correcting Poisson measurement errors, especially for those data observed from the high throughput sequencing technologies. Measurement error is an important issue in statistical analysis, that occurs when we want to analyse some latent variable, but our attempts to measure the variable are subject to random error. Our framework is conceptually different from that of Collins et al. (2001) and Chiquet et al. (2018). A closely related approach is the work of Liu et al. (2018) who develop the same estimator in the untransformed case without sequencing depth noise, as we do, but then proceed to develop shrinkage estimators to improve performance, rather than deal with sequencing depth noise or transformations. We assume that the latent Poisson means do not lie exactly on a low-dimensional space, but instead follow some unknown distribution, and we wish to estimate the principal components of this distribution or of a non-linear transformation of this distribution. When no transformation is needed, a PCA on an unbiased estimate of the covariance matrix of the latent Poisson means is straightforward. When a non-linear transformation is applied on the latent Poisson means, our idea is to derive a non-linear transformation of the data from which an unbiased estimator of the covariance matrix of the transformed latent Poisson means can be derived and thus the PCA can be applied on this unbiased estimate of the covariance. Our motivation for extending our work to transformed latent Poisson means is that previous research suggests that a logarithmic scale is an appropriate way to study microbial abundance. This also makes sense from the point of view of differential equation models of species abundance, which tend to imply exponential-type growth. We have avoided any parametric constraints on this latent distribution. Although our research is driven by the application to microbiome data which we will describe in more detail in the following subsection, the method is geneally applicable to high dimensional count data which are popularly observed in many different areas including social science, ecology, microbiology and bioscience in general.

The structure of the paper is as follows. In Section 2, we introduce the features of microbiome data which motivated this research. In Section 3, we develop the bias corrected estimators for the covariance matrices of functions of the latent Poisson means. Then in Section 4 we develop methods for handling the sequencing depth issue, particularly in the $\log$-transformed case. In Section 5, we provide a method for projecting the data onto the principal component space. We then demonstrate the performence of our methods on simulated data in Section 6 and on a real microbiome data set (Caporaso et al, 2010) in Section 7. Finally in Section 8, we summarise our work and suggest promising directions for future research.

The methods in this paper are all implemented in the R package PoissonPCA, which is available from the first author’s website www.mathstat.dal.ca/~tkenney/PoissonPCA/, and will be submitted to CRAN.

2 Features of microbiome data

The microbiome is the collection of all microorganisms in a particular environment. This collection mostly consists of bacteria, archaea, eukaryotes and viruses. Recent improvements in technology have allowed the bacteria and archaea components to be measured en-masse. That is, identifying DNA from all bacteria and archaea present in a given sample can be collected and counted, with the objective of giving a complete picture of the whole microbial community. This picture can then be compared to metadata, such as whether the microbial community was from a healthy individual or from an individual suffering from a particular illness. The idea is that the interaction of the thousands of microbes present in the whole community may be more informative than simple analyses that try to attribute complex conditions to the actions of a single species. There have been a large number of papers that have identified associations between the microbiome and various environmental conditions such as nutrient abundance (Arrigo 2005), and human health (Fujimura et al. 2010, Sekirov et al. 2010) including obesity (Turnbaugh et al. 2009), Crohn’s disease (Quince et al. 2013), diabetes (Kostic et al. 2015).

The microbiome data present a number of interesting data analysis challenges that make the application of PCA difficult. The feature of microbiome data that is the focus of this paper is measurement error. While each environment has a community with a certain total abundance of each Operational Taxonomic Unit (OTU), we only observe a small sample from that community, which for this paper, we will model as following a Poisson distribution with mean given by the underlying total abundance. In practice the sequencing procedure is more complicated (there are several different sequencing procedures that can be used, each of which involves their own noise and bias — Gorzelak et al. (2015) give more details about how the sequencing procedure influences the data) but the Poisson noise error seems a reasonable starting point. The aim of this paper is to adapt PCA to account for this Poisson measurement error.

A further feature of microbiome data is sequencing depth. Essentially, the issue is that the Poisson means of our observed sample are subject to large multiplicative noise. For this reason, microbiome data is usually treated as compositional data, with only the proportions of each OTU considered, for example Li (2015). However, this loses a lot of information about the Poisson distribution. For the Poisson distribution, the relative noise in the observed proportions will decrease as sequencing depth increases, so samples at higher sequencing depth should be given more weight. Another commonly used approach to the sequencing depth problem is rarefaction, where large parts of the data are thrown away, so that all samples have the same total OTU count. This practice is mainly due to the use of arbitrary indices which are sensitive to sample size. McMurdie and Holmes (2014) argue strongly against the use of rarefaction.

The other features of microbiome data include the high dimensionality and sparsity. There is substantial literature on ways to deal with these issues by imposing sparsity constraints on the PCA. For example Zou et al. (2006), and Jolliffe et al. (2003). Salmon et al. (2014) combine the Poisson PCA of Collins et al. (2001) with an $l^{1}$ penalty to create a sparse Poisson PCA, and apply this method to Photon-limited image denoising. However, the sparseness issue is not the focus for the current paper, and we will avoid it by grouping our OTUs at a higher taxonomic level, i.e. instead of considering variables that roughly correspond to the abundance of a particular species of microbe, we will look at higher-level variables that correspond to the total abundance of a whole genus or higher taxonomic level. This will reduce the dimension enough to avoid the curse of dimensionality for the datasets we consider. In principle the $l^{1}$ penalty could be applied to our method to provide a sparse Poisson-error corrected PCA, but that is a topic for future work.

There has been other work incorporating the Poisson error structure of microbiome data into the analysis, but it has been focused on parametric modelling, rather than Exploratory data analysis. For example Holmes et al. (2012) model microbiome data as a Dirichlet multinomial distribution; Knights et al. (2011) use a Bayesian mixture model; Shafiei et al. (2014) and Shafiei et al. (2015) develop hierarchical Bayesian multinomial mixture models; and Cai et al. (2017) apply non-negative matrix factorisation with Poisson error.

3 Method

We express the situation more formally as follows: We have $n$ samples of a $p$-dimensional random vector $\mathbf{X}$ which follows a Poisson distribution with conditional mean given by a latent random vector $\mathbf{\Lambda}$ and a “sequencing depth” $S$ (this comes from our specific application to microbiome data). That is, we have $\mathbf{X}_{i}|\mathbf{\Lambda}_{i},S\sim Po(\mathbf{\Lambda}_{i}S)$, $\mathbf{X}_{i}$’s $(i=1,\cdots,p)$ are conditionally independent given the random vector $\mathbf{\Lambda}$ and the sequencing depth $S$. $\mathbf{X}$ can be thought of as $\mathbf{\Lambda}$ with some Poisson measurement error. Our samples form an $n\times p$ matrix $X$ of observed data, where the $i$th row is an independent sample of the random vector $\mathbf{X}$. We want to estimate the distribution of $\mathbf{\Lambda}$ and its relationship to other variables. For this purpose, we develop a method of estimating a PCA on this latent $\mathbf{\Lambda}$, or on some transformation of $\mathbf{\Lambda}$, $f(\mathbf{\Lambda})$. Our approach to this problem is moment-based — we aim to correct the bias in the variance estimates, and use the unbiased variance estimates to perform PCA.

We start with the linear case, where we want to estimate the PCA of $\mathbf{\Lambda}$. There are two cases to consider: the case without sequencing depth noise (where $S_{i}=1$ for all $i=1,\cdots,n$), and the case where sequencing depth correction is needed. We then develop the methods for estimating the PCA of a nonlinear transformation $f(\mathbf{\Lambda})$.

3.1 Poisson Noise Correction in PCA without Sequencing Depth Corrections

We first consider the simplest case, where $S_{i}=1$ for all $i=1,\cdots,n$. That is, any differences in total abundance between samples is a result of differences in underlying latent abundances. The model here is fairly straightforward. Latent Poisson means are a random vector $\mathbf{\Lambda}$, and observed counts are a random vector $\mathbf{X}$ with $\mathbf{X}_{i}|\mathbf{\Lambda}_{i}\sim Po(\mathbf{\Lambda}_{i})$ conditionally independent given $\mathbf{\Lambda}$. We want to estimate the principal components of the random vector $\mathbf{\Lambda}$ from a matrix $X$ whose rows are independent samples of the random vector $\mathbf{X}$. The principal components of $\mathbf{\Lambda}$ are the eigenvectors of the variance of the vector $\mathbf{\Lambda}$, so we only need to provide an estimate for the variance of the random vector $\mathbf{\Lambda}$. We have

[TABLE]

where $\operatorname{diag}(\mathbf{\Lambda})$ is a diagonal matrix with entries $\mathbf{\Lambda}$. Now we have the following straightforward estimators for $\operatorname{Var}(\mathbf{X})$ and $\operatorname{diag}\left({\mathbb{E}}(\mathbf{X})\right)$:

[TABLE]

plugging these in gives the estimator

[TABLE]

3.2 Poisson Noise Correction in PCA with Sequencing Depth Corrections

Suppose we now add sequencing depth noise to the problem, so that the Poisson mean is $S\mathbf{\Lambda}$, where $S$ is a random scalar and $\mathbf{\Lambda}$ is a random vector. To avoid the identifiability issue we assume $\mathbf{\Lambda}^{T}1=1$, and that $\mathbf{\Lambda}$ and $S$ are independent. We can apply the same procedure

[TABLE]

Thus

[TABLE]

Now since $\mathbf{X}_{i}|\mathbf{\Lambda}_{i},S\sim Po(\mathbf{\Lambda}_{i}S)$, $(i=1,\cdots,p)$, we have

[TABLE]

Thus

[TABLE]

We can now empirically estimate this quantity from a sample. Suppose we have a matrix $X$, where each row is a sample from the underlying distribution of $\mathbf{X}$. If the sequencing depth $S$ is known, then we can compute $\frac{X}{S}$ for each sample, and we have empirical estimators

[TABLE]

where $D$ is a diagonal matrix with entries the sequencing depths of each observation. Plugging this into our formula for $\operatorname{Var}(\mathbf{\Lambda})$ gives us the estimator

[TABLE]

The correction here is similar to the correction in (Li, Palta and Shao, 2004), where they added a correction term $\frac{{\beta_{1}}^{2}}{2{\sigma_{\epsilon}}^{2}}\sum_{i=1}^{N}\frac{W_{i}}{{L_{i}}^{2}}$ to the likelihood term in a regression model, where $W_{i}$’s are the Poisson samples, and $L_{i}$’s are the sequencing depths. Their case is a regression model, and in their case the sequencing depth is more directly observed (being number of hours of sleep).

3.3 Poisson PCA on Transformed Poisson Means without Sequencing Depth Corrections

Suppose now that we are interested in the covariance matrix of some transformation $f(\mathbf{\Lambda})$, where $f$ is some univariate function applied elementwise to $\mathbf{\Lambda}$. For example, we might set $f(\mathbf{\Lambda})=\log(\mathbf{\Lambda})$. We can use the same technique as in Section 3.1, but with some univariate transformation $g(\mathbf{X})$ of the original data. Using the law of total variance, we have:

[TABLE]

Now if we can find $g(\mathbf{X})$ so that ${\mathbb{E}}(g(\mathbf{X})|\mathbf{\Lambda})=f(\mathbf{\Lambda})$ then we get

[TABLE]

The first step, therefore is finding $g(\mathbf{X})$ such that ${\mathbb{E}}(g(\mathbf{X})|\mathbf{\Lambda})=f(\mathbf{\Lambda})$. Without loss of generality, we take $\mathbf{\Lambda}$ as univariate. We can directly calculate ${\mathbb{E}}(g(\mathbf{X})|\mathbf{\Lambda})$ as follows:

[TABLE]

That is, we can take $\frac{g(n)}{n!}$ to be the coefficients of the Taylor series of $e^{\lambda}f(\lambda)$. As an example, suppose we want $g(\mathbf{X})$ to give an unbiased estimator for $f(\mathbf{\Lambda})=\mathbf{\Lambda}^{k}$ for some fixed $k$. We have

[TABLE]

We therefore get our estimator $g(\mathbf{X})=\mathbf{X}(\mathbf{X}-1)\cdots(\mathbf{X}-k+1)$, which is a well-known unbiassed estimator for $\lambda^{k}$.

Next we need to calculate the conditional variance $\operatorname{Var}(g(\mathbf{X})|\mathbf{\Lambda})$. Since the elements of $\mathbf{X}$ are conditionally independent given $\mathbf{\Lambda}$, this matrix is a diagonal matrix and it will suffice to take $\mathbf{\Lambda}$ as univariate for the estimation of each element of this matrix. We can estimate this from first principles:

[TABLE]

where $h(n)=\sum_{m=0}^{n}\binom{n}{m}g(m)g(n-m)$ is the coefficient of the Taylor series $e^{2\lambda}f(\lambda)^{2}=\sum_{n=0}^{\infty}h(n)\frac{\lambda^{n}}{n!}$. Since $\mathbf{\Lambda}$ is latent, we need to replace the expressions $e^{-\lambda}\frac{\lambda^{n}}{n!}$ and $e^{-2\lambda}\frac{\lambda^{n}}{n!}$ by estimators which are functions of $\mathbf{X}$. Since we will be taking the expectation of our estimators of the conditional variance, it makes sense to focus on the bias of these estimators. For $e^{-\lambda}\frac{\lambda^{n}}{n!}$, there is an unbiassed estimator

[TABLE]

which can reasonably be used. This also has a computational advantage because it results in only needing to compute a single term in the sum. For the term $e^{-2\lambda}\frac{\lambda^{n}}{n!}$, the unbiassed estimator is

[TABLE]

This gives us the overall estimator for the variance of $f(\mathbf{\Lambda_{l}})$, where $\mathbf{\Lambda_{l}}$ is an element of the vector $\mathbf{\Lambda}$.

[TABLE]

The estimates of the off-diagonal elements of the covariance matrix can be computed as the corresponding elements of the sample covariance matrix of $g(X)$, where $g$ is applied elementwise to the observed data matrix $X$.

In Appendix A, we show the consistency of this estimator:

Theorem 3.1**.**

If $e^{2\lambda}f(\lambda)^{2}=\sum_{n=0}^{\infty}\frac{h(n)\lambda^{n}}{n!}$ is a globally convergent Taylor series, $\mathbf{\Lambda}$ is a non-negative random variable with finite raw moments $\mu_{n}$ and $\sum_{n=0}^{\infty}\frac{\mu_{n}|h_{n}|}{n!}$ converges, then the estimator (1) is consistent.

In practice, this estimator can have high variance, leading to poor finite-sample performance for certain transformations. Sometimes this issue can be alleviated by using suitable approximations, and simplifying the algebra in equation (1). We show how this can be done for the commonly used log transformation in Appendix B. We also suggest a possible approach to the more general case, that could be given further study in Appendix C.

3.4 The Special Case of Log Transformation

In the special case when $f(\lambda)=\log(\lambda)$, the approach from Section 3.3 does not work, because $f$ does not have a globally convergent Taylor series. We can compute a Taylor series for $\log(\lambda)=\log\left(a\left(1+\frac{\lambda-a}{a}\right)\right)$ for some value of $a$. That is, we get

[TABLE]

This Taylor series is convergent for $0<\lambda<2a$, however, when we expand the terms and try to reorganise it as a power series in $\lambda$, the coefficients do not converge. Therefore, we can only use a truncated power series. That is, we truncate the above series at some chosen value $N$, to get a polynomial. This produces a reasonable approximation to $\log(\lambda)$ on the interval $(0,2a)$. For larger values of $\lambda$, we have that ${\mathbb{E}}(\log(\mathbf{X})|\mathbf{\Lambda}=\lambda)\approx\log(\lambda)$. Therefore, we set $g(X)=\log(X)$ instead for larger values of $X$. The result of this is that ${\mathbb{E}}(g(\mathbf{X})|\mathbf{\Lambda})\approx\log(\mathbf{\Lambda})$ for $\mathbf{\Lambda}$ not too close to zero. In summary our function $g(n)$ is given by

[TABLE]

There is some difficulty choosing the values of $a$ and $N$. From experiments, we see that $a=3,N=4$ and $g(X)=\log(X)$ whenever $X\geqslant 7$, gives a fairly good estimator. We discuss the corresponding estimator for conditional variance in Appendix B.

4 Sequencing Depth Correction

The methods described above for variance estimators of transformed latent Poisson means cover the situation without sequencing depth correction. In this section, we give two approaches to correcting the resulting covariance matrix for sequencing depth.

4.1 Method 1: Compositional Covariance Matrix for Sequencing Depth Correction

The first method is to get a compositional covariance matrix. That is, suppose the estimated covariance matrix is $\Sigma$. We want to find the matrix $\tilde{\Sigma}$ with the following properties:

For any $u,v$ with $u^{T}1=v^{T}1=0$, we have $u^{T}\tilde{\Sigma}v=u^{T}\Sigma v$ 2. 2.

$\tilde{\Sigma}$ is symmetric. 3. 3.

$\tilde{\Sigma}1=0$

Conditions 2 and 3 are necessary for $\tilde{\Sigma}$ to be the covariance matrix of compositional data. Condition 1 says that the covariance matrices $\tilde{\Sigma}$ and $\Sigma$ give the same covariance for the projection of $f(\mathbf{\Lambda})$ onto the space orthogonal to the vector $1$. Now to solve this, we have the following.

Proposition 4.1**.**

The matrix $\tilde{\Sigma}_{p\times p}$ that satisfies Conditions 1–3 is given by

[TABLE]

where $a=(pI+11^{T})^{-1}\Sigma 1$.

Proof.

We first note that for any $v$ with $v^{T}1=0$, we have $\tilde{\Sigma}v-\Sigma v=\alpha(v)1$ for some scalar $\alpha(v)$. This is because $u^{T}\left(\tilde{\Sigma}v-\Sigma v\right)=0$ for any $u$ with $u^{T}1=0$. By linearity, $\alpha(v)$ is a linear function of $v$, so we have $\alpha(v)=a^{T}v$ for some fixed vector $a$. Now we have $\left(\tilde{\Sigma}-\Sigma-1a^{T}\right)v=0$ for all $v$ with $v^{T}1=0$. If we let $\Sigma_{r}=\tilde{\Sigma}-\Sigma-1a^{T}$, then this gives $\Sigma_{r}x=\Sigma_{r}\frac{1}{p}11^{T}x$ because $\frac{1}{p}11^{T}x$ is the orthogonal projection of $x$ onto the 1-dimensional space spanned by 1. Therefore we have $\Sigma_{r}=\Sigma_{r}\frac{11^{T}}{p}=b1^{T}$ where $b=\frac{\Sigma_{r}1}{p}$. Thus we have $\tilde{\Sigma}-\Sigma-1a^{T}=b1^{T}$ so $\tilde{\Sigma}=\Sigma+1a^{T}+b1^{T}$ Since $\tilde{\Sigma}$ and $\Sigma$ are both symmetric, we get $a=b$, which gives the form in the proposition. All that now remains is to find the value of $a$. This is forced by $\tilde{\Sigma}1=0$ which gives

[TABLE]

∎

We have that $(pI+11^{T})(cI+d11^{T})=cpI+(c+pd+pd)11^{T}$, so solving for $cp=1$ and $c+2pd=0$ we get that $(pI+11^{T})^{-1}=\frac{I}{p}-\frac{11^{T}}{2p^{2}}$, so

[TABLE]

4.2 Method 2: Minimum Variance by Correcting for the Maximum Sequencing Depth Noise

The second method is specific to the $\log(\mathbf{\Lambda})$ case. In this case, the sequencing depth becomes additive noise on $\log(\mathbf{\Lambda})$, proportional to the vector 1. This means that $\Sigma=\Sigma^{*}+c11^{T}$. The value of $c$ is not identifiable, but it is limited by the constraint that $\Sigma^{*}$ should be non-negative definite. One natural choice is therefore to select the largest value of $c$ such that $\Sigma^{*}$ remains non-negative definite. In the case where $1$ is an eigenvector of $\Sigma$, this will give exactly the same result as the compositional method. In other cases, the results may differ slightly.

To solve for the largest value of $c$ such that $\Sigma-c11^{T}$ is non-negative definite, we consider the determinant $\left|\Sigma-c11^{T}\right|$ as a continuous function of $c$. While $\Sigma-c11^{T}$ is positive definite, the determinant is positive. We therefore set $c$ to be the smallest solution to $\left|\Sigma-c11^{T}\right|=0$. We solve this with a change of basis $\Sigma=ASA^{T}$ where $A$ is an orthogonal matrix with $Ae_{1}=\frac{1}{\sqrt{p}}1$, where $e_{1}=(1,0,\cdots,0)^{T}$. Now we have that

[TABLE]

where $S_{\hat{1}\hat{1}}$ is the matrix $S$ with the first row and column removed and the first entry replaced by $1$. For example, if

[TABLE]

then

[TABLE]

More formally $S_{\hat{1}\hat{1}}=S-e_{1}s^{T}-s{e_{1}}^{T}+(S_{11}+1)e_{1}{e_{1}}^{T}$ where $s=Se_{1}$ is the first column of $S$. This gives

[TABLE]

This gives us that

[TABLE]

where $\widetilde{\Sigma}$ is the compositional estimate from Section 4.1. So we are trying to solve for $c$ such that

[TABLE]

Thus

[TABLE]

In practice, this method does not perform well, because the small eigenvalues and eigenvectors are not well estimated. Indeed some of the estimates are often negative. We therefore improve practical performance by introducing some thresholding. Let $\Sigma=VDV^{T}$ where

[TABLE]

with $\lambda_{1}>\lambda_{2}>\cdots>\lambda_{p}$. Since $\Sigma$ should be non-negative definite, we first replace any negative eigenvalues by 0. Next we choose the smallest $j$ such that $\frac{\lambda_{2}+\cdots+\lambda_{j}}{\lambda_{2}+\cdots+\lambda_{p}}>0.9$. The idea here is that we assume that $\lambda_{1}$ is influenced mostly by sequencing depth and take eigenvalues accounting for 90% of the variance we are trying to estimate. We now replace all smaller eigenvalues $\lambda_{k}$ for $k>j$ by $\lambda_{j}$. We now have a new diagonal matrix

[TABLE]

Then we apply the above minimum variance sequencing depth correction procedure to $V\widetilde{D}V^{T}$.

5 Projection onto Principal Component Space

Having estimated the principal components of the latent Poisson means, the remaining problem is to “project” the observations onto this principal component space. The methods in this paper are based on a semiparametric framework. We are assuming that the latent Poisson means follow some unknown distribution, and are seeking to use non-parametric methods to estimate the principal components of this latent distribution. In this spirit, our projection method should be based on the same assumptions. We will assume for each observed point, that there is some latent random vector $\mathbf{\Lambda}$ of Poisson means, and that we want to project these Poisson means onto the principal component space. In the non-latent variance of principal components, we have an observation $X$, and we seek the projection $X^{*}$ in principal component space that minimises $(X-X^{*})^{T}(X-X^{*})$. In our case, we have a latent random vector $\mathbf{\Lambda}$, and have calculated the principal component space of $f(\mathbf{\Lambda})$. If we knew $\mathbf{\Lambda}$, we would seek the projection $R$ in the principal component space that minimises $\left(f(\Lambda)-R\right)^{T}\Sigma^{-1}\left(f(\Lambda)-R\right)$, where $\Sigma$ is the covariance matrix of $f(\mathbf{\Lambda})$. (Because the principal component space is spanned by eigenvectors of $\Sigma$, it turns out that it is equivalent to take the projection $R$ that minimises $\left(f(\Lambda)-R\right)^{T}\left(f(\Lambda)-R\right)$.)

However, we do not observe $\mathbf{\Lambda}$ directly, instead, we observe $\mathbf{X}\sim Po(\mathbf{\Lambda})$. We therefore simultaneously expect $\mathbf{\Lambda}$ to be such that $\mathbf{X}$ has high likelihood, and such that $f(\mathbf{\Lambda})$ is close to the principal component space. We therefore create the objective function

[TABLE]

where $\log(\Lambda)$ is the vector whose $i$th element is $\log(\Lambda_{i})$ — that is it is the logarithm function applied elementwise to the vector $\Lambda$. Note that because we have computed the covariance matrix of $f(\Lambda)$, there is a canonical choice of how to balance the two constraints, and we do not need a tuning parameter. If the latent variable $f(\Lambda)$ is normal, then $L(\Lambda)$ is the negative log-likelihood of the joint distribution of $\mathbf{X}$ and $\mathbf{\Lambda}$ (treating $R$ as parameters).

To solve this, we substitute the linear model $f(\Lambda)=\mu+Va$ where $V$ is the matrix of principal component vectors. That is $\Sigma=VDV^{T}$. Suppose that our projection is onto the first $r$ principal components and let

[TABLE]

That is, $a_{[r]}$ is the vector of the first $r$ elements of $a$. Then we have $R=\mu+Va_{[r]}$. In this basis, the objective function becomes

[TABLE]

where $D$ is the diagonal matrix of eigenvalues of $\Sigma$, in the order corresponding to $V$. The estimate for $a$ is therefore given by setting the derivative of this to zero. That is

[TABLE]

In the particular case where $f(\Lambda)=\log(\Lambda)$ (evaluated elementwise) this equation becomes

[TABLE]

Differentiating again, the Hessian matrix for the case where $f(\Lambda)=\log(\Lambda)$ is

[TABLE]

where $S=VD^{-1}I_{i\geqslant k}V^{T}$ Newton’s method now gives

[TABLE]

6 Simulation

6.1 Simulation Design

We assess the performance of the method under three different transformations — the default untransformed Poisson means, a polynomial transformation $f(x)=0.004x^{3}-0.4x^{2}+16x$, and a logarithmic transformation. In each case, we set the eigenvalues of the covariance matrix to a fixed sequence $d$, and simulate the eigenvectors $V$ uniformly. We now simulate the true transformed covariance matrix as $M=VDV^{T}$, where $D$ is a diagonal matrix with diagonal entries $d$. We simulate an $n\times p$ data matrix $T$ of transformed Poisson means, where each row of $T$ is i.i.d. with a chosen distribution (usually multivariate normal) with covariance matrix $M$. We compute the Poisson means from $T$ by $\Lambda_{ij}=f^{-1}(T_{ij})$. Finally, we simulate $X_{ij}\sim Po(\Lambda_{ij})$. For each scenario, we simulate 100 datasets. We use the same eigenvalues, but different eigenvectors in the 100 simulated covariance matrices. We also simulate different means for $T$ for each simulated dataset, using a normal distribution with mean 100 and standard deviation 15.

For simulations with sequencing depth correction, in the compositional case, we modify the above procedure by setting $d_{p}=0$, and $v_{p}=1$, with other eigenvectors in $V$ simulated orthogonally to 1. In the non-compositional case, simulate $\Lambda$ as above. In both cases, we then simulate sequencing depth noise following a gamma distribution with shape parameter 4 and scale parameter 100.

6.2 Simulation Results

6.2.1 Assessing Simulation Results

We want to devise a method to assess the performance of our corrected PCA on the simulated data. In this case, we have a true covariance matrix $\Sigma=VDV^{T}$, where $V$ is orthogonal and $D$ is diagonal. The principal component matrix $V$ is defined by the property of simultaneously maximising $\sum_{i=1}^{k}(V^{T}\Sigma V)_{ii}$ for all $k$, subject to being orthogonal. We have obtained an estimated covariance matrix $\widehat{\Sigma}=\widehat{V}\widehat{D}\widehat{V}^{T}$ with $\widehat{V}$ orthogonal and $\widehat{D}$ diagonal. We want to assess the performance of $\widehat{V}$ as an estimate for $V$. Since our objective is to maximise the variance explained by the first $k$ principal components, we will measure this by comparing cumulative true variance explained by the first $k$ estimated principal components, i.e. the cumulative sums $\sum_{i=1}^{k}(\widehat{V}^{T}\Sigma\widehat{V})_{ii}$, with the true cumulative sum $\sum_{i=1}^{k}(V^{T}\Sigma V)_{ii}$.

6.2.2 Simulation results for untransformed PCA without sequencing depth correction

The simulation results for normally distributed Poisson means without sequencing depth noise with data dimensions equal to 80 are shown in Figure 1. We compare the percentages of the total variance explained by the true eigenvectors, the eigenvectors found by our Poisson PCA method, and the eigenvectors found by “naive PCA” i.e. applying PCA directly to the observed data, without any correction for Poisson error. We see that when $n=100$, our method is comparable with naive PCA, but as sample size increases, our method performs much better, as expected. For large $n$, our method appears to be converging to the correct eigenvalues as $n\rightarrow\infty$, while uncorrected PCA is not improving with larger sample sizes, because of the bias. Similar results when $p=30$ are shown in Supplementary Figure LABEL:linearsimresultsp30.

6.2.3 Simulation results for Polynomial Transformed PCA without sequencing depth correction

Figure 2 shows the simulation results for polynomial-transformed Poisson means. Again, our method is comparable to other methods for small sample sizes, and converges towards the truth as sample size increases. In this case, the larger variance of the estimators means that covergence of our method is slower than in the untransformed case. However, Poisson PCA begins to outperform other methods for sample size around 5,000, which is a plausible sample size for some real data sets. Supplemental Figure LABEL:polynomialsimresultsp30 shows similar results for $p=30$.

6.2.4 Simulation results for Log Transformed PCA without sequencing depth correction

Figure 3 compares Poisson PCA to PLNPCA (Chiquet et. al, 2018), which is a parametric method based on a Poisson log-normal assumption using a variational approximation to estimate the parameters by MLE, and to PCA directly performed on the log-transformed data (with 0 counts arbitrarily replaced by 0.1). We compare the methods on both log-normal and non log-normal Poisson means. For the non log-normal case, we simulated some data in polar coordinates with direction uniformly distributed over the sphere and radius distributed following a centred gamma distribution with shape $\alpha=5$ and scale $\theta=\frac{1}{\sqrt{5}}$. The data were then stretched to have the proscribed variance. Similar results with $p=30$ are shown in Supplemental Figure LABEL:logsimresultsp30.

When sample size is small, Poisson PCA clearly outperforms both PLNPCA methods and the naive PCA on log-tranformed data in the first few dimensions which capture about 70% to 80% of the total variance. After this point, the total variance explained by naive PCA on log-tranformed data starts to get higher. Between the two PLNPCA methods, the full rank log-Normal Poisson PCA seems to outperform the log-normal Poisson model with rank chosen by BIC. When sample size is large ($n=1,000$) the Poisson PCA and two PLNPCA methods show comparable results and perform better than the naive PCA on log-tranformed data. In practice, when we use PCA, we aim in particular to identify the first few principal components, so better performance at identifying the earlier prinicipal components is a very strong asset of our method. All the methods are robust to changing the underlying distribution of the Poisson means. Results for $n=5,000$ are shown in Supplemental Figure LABEL:logsimresultsn5000; we do not give results for larger $n$ because of the computation time needed by PLNPCA.

In addition to the better identification of the early principal components, Poisson PCA also has a major advantage in computation time. Figure 4 compares the logarithm of the average processor time (amount of computation performed) over 100 simulated datasets, for the three methods. For the limited rank PLNPCA, with rank chosen by BIC, we compared all possible ranks. Limiting the rank would improve the runtime of this method, to make it more comparable to the direct PLNPCA estimate. Poisson PCA is faster than the likelihood-based method by an order of magnitude. The runtimes are approximately linear in $log(p)$, so we have shown regression lines. Our Poisson PCA method is not quite linear, probably because the singular value decomposition, performed after variance estimation, takes nonlinear time. Because Poisson PCA is much faster than the other methods, this nonlinearity is more apparent for our method.

6.2.5 Simulation results for untransformed PCA and log-transformed PCA with sequencing depth correction

Figure 5 compares Poisson PCA with compositional PCA (i.e. data are rescaled so each row has total 1, and PCA is applied to the rescaled data). The scaled Poisson means are compositional, and follow a multivariate normal distribution; but we have multiplied each Poisson mean by a random sequencing depth. Our semiparametric method has done better than compositional PCA at identifying the principal components.

Figure 6 shows the results for log-normal simulated data with sequencing depth noise. In the compositional log Poisson means simulation, we use the log of observed total counts as an offset in the parametric PLN model. We see that this is a challenging problem, but that our compositional sequencing depth correction method (from Section 4.1) offers a slight improvement over the PLNPCA method. We do not include larger sample sizes because of the computation time needed to fit the PLNPCA model. We also compare our two methods for correcting sequencing depth noise when the latent covariance matrix is not compositional. We see that the minimum variance approach (from Section 4.2) has not done well at identifying the first eigenvector, because it has not removed enough sequencing depth noise. However, retaining the signal in the direction 1 allows it to get improved results for later principal components. Overall it seems that in most circumstances it would be better to use the compositional method, but in rare cases where we are aiming to find a large number of principal components, it may be preferable to use the minimum variance approach. Additional simulation results are shown in Supplemental Figure LABEL:logsimresultsseqdepthsupplemental.

6.3 Projection onto Principal Component Space

We now test the performance of the projection method from Section 5 on simulated data. We use our standard log-normal Poisson simulation without sequencing depth correction. We compare two projection approaches: (a) pure normal error — we take the log of the observed data and project this onto the principal component space (substituting $-3$ for $\log(0)$); (b) the combined likelihood method from Section 5. We apply each of these projections both using the true covariance matrix, and using the estimated covariance matrix from Poisson PCA.

For assessing the quality of our projection, we have the true transformed Poisson mean $\log(\Lambda)$, so we can project this Poisson mean onto the principal component space, and measure the distance between the correct projection and the estimated projection. Since the projection is orthogonal, Euclidean distance suffices for this measure.

Figure 7 compares the results on both the true principal component space, and the estimated principal component space. We see that when we project onto the true principal component space, our projection method outperforms the naïve projection method for all numbers of principal components. When we project onto the estimated principal component space, we see that our method performs excellently for any projection that might ever be used in practice. However, when we try to project onto high dimensional component spaces, our method starts to perform less well. In this case, because there is very little signal in the remaining dimensions, our method may soak up some of the noise. Optimisation is also more challenging in higher dimensions, so inaccuracies in the optimisation may cause deterioration in performance. Since the typical use of PCA is to obtain a low-dimensional representation of the data, performance on small numbers of principal components is most important. Therefore, the results of our method are good. Similar results for $p=30$ are in Supplementary Figure LABEL:ProjectionSimulationResultsp30.

7 Real Data Analysis

7.1 The Moving Picture Dataset

The moving picture dataset (Caporaso et al. 2011) is a large temporal microbiome dataset comprising samples from 4 body sites of two individuals over a period of approximately two years. There are over 300 samples from each site for individual 2 and over 100 samples from each site for individual 1. Previous analysis of this dataset has shown that it is relatively easy to distinguish the two individuals from the species-level OTUs, but using higher-level OTUs to distinguish them is more difficult. We apply our Poisson PCA method across a range of levels from phylum to genus from this dataset. Besides making the problem challenging compared to the results of previous studies, the focus on higher level OTUs also ensures that the number of variables is less than the number of observations, which is a necessary condition for the consistency of classical PCA. We present the results for the first two principal components at the order level (225 variables) in this section. Results at other levels are in the supplemenary materials.

There is substantial variation in sequencing depth between samples, with number of total reads varying between less than 10,000 and over 60,000. The distributions of sequencing depths were different for the two individuals. Since sequencing depth is largely considered noise due to experimental procedure, a method which separates samples based on sequencing depth is not a good method. To prevent methods that fail to correctly compensate for sequencing depth from wrongly seeming good, we used subsampling to ensure that the sequencing depths had the same distribution for each individual. (Details of the subsampling procedure are in Supplemental Appendix LABEL:RealDataSubsamplingProcedure.) We retain a large variability of sequencing depth within each individual’s samples, but remove the ability to use sequencing depth as a distinguishing feature, so that methods that better correct for sequencing depth noise are expected to perform better on the subsampled data.

7.2 Data Analysis

We compare three methods for dimension reduction of the log-transformed data: our compositional sequencing depth corrected Poisson PCA method using the projection method from Section 5; PLNPCA with rank chosen by BIC in the range 1–30, using total count as an offset; and log compositional PCA, i.e. naive PCA applied to the log of compositional data with zero counts replaced by 0.01 before converting to compositional data (similar results are obtained using other values to replace 0). Since the differences between the two individuals’ microbial communities are likely to be a large source of variability in the data, we expect to see a separation between samples from the two individuals after performing PCA on the data from each body site.

Results for the gut data are in Figure 8. We see that all three methods show a reasonable separation between the two individuals. PLNPCA shows several outliers, which can be explained by the difficulty of attempting to fit data to a prescribed parametric distribution if the distribution does not fit.

Figure 9 shows the projection of the tongue data. For this dataset, Poisson PCA shows a better separation between two individuals. PLNPCA has again been heavily influenced by a few outliers. Log compositional PCA shows a clear clustering, but the clustering is not related to the individual. Further examination reveals that the first principal component for both individuals is heavily influenced by a single order from phylum Cyanobacteria (see Supplementary Figure LABEL:Cyanobacteria). This applies to both Poisson PCA and log compositional PCA. However, it is more obvious in the case of log compositional PCA, and there is a clear gap between cases where the observed count of cyanobacteria is zero and where it is non-zero, which leads to the two clusters seen in Figure 9(a). The abundance of Cyanobacteria does not vary much between the two individuals, so because this abundance dominates the first principal component for log compositional PCA, we do not distinguish between the individuals using this method. Because Poisson PCA does not attach so much weight to the zero values, it is able to identify some signal, even with the Cyanobacteria abundance dominating the first principal component.

Figure 10 shows the results for the left palm. Previous research (e.g. Shafiei et al. 2015, Cai et al. 2016) has found these samples harder to classify than the gut or tongue data, so it is not surprising that all PCA methods show more mixing of the two groups. However, we see the same patterns as for the tongue data, with Poisson PCA producing the best separation, PLNPCA being adversely influenced by outliers and log compositional PCA producing a false clustering. This time, the false clustering can be seen to be closely related to sequencing depth (see Supplementary Figure LABEL:MovPicLPalmOrderSeqDepth).

Figure 11 shows the results for the right palm. Previous studies have found the right palms the most difficult locations to distinguish between the individuals, and the PCA results confirm this, with substantial mixing of the groups using all PCA methods, with log compositional PCA perhaps showing the best separation (certainly the best linear separation).

8 Conclusions and Future Work

8.1 Conclusions

We have developed methods for correcting for Poisson measurement error in estimating the variance of a latent variable. Our methods can estimate the variance of the latent Poisson means, or the variance of some transformation of the latent means. We have also incorporated ways to adjust our methods to deal with additional multiplicative noise in the Poisson means (sequencing depth). We have tested our methods across a range of scenarios and compared them with existing methods. We found that our method performs well on both simulated and real data, and is not so heavily influenced by outliers as parametric methods, and successfully adjusts for sequencing depth unlike naïve PCA on the logarithm of the compositional data.

While our method has focused on Poisson noise, so is less broad than some parametric methods that extend to exponential family distributions, our work is broader in that it can be applied to arbitrary transformations of the data, not just the logarithm. The ideas behind our work should easily be applicable to other families of noise distribution, particularly discrete. We have also developed methods for dealing with sequencing depth, that are not available in other methods. Furthermore, our method is an order of magnitude faster than parametric methods, because it does not require optimisation of all parameters.

8.2 Future Work

There are a number of avenues for extending this research. Firstly, for the log transformation, we were not able to find a perfect unbiased estimator, and instead used an approximation. Further work to find a better approximation and corresponding estimators should lead to better performance of the method. Secondly, while we have focussed on Poisson measurement error correction, the same ideas could be applicable to other measurement error correction for PCA. Our method was based on the fact that the Poisson distribution has one parameter, so the measurement error variance can be directly estimated from the data. This technique could be applied to other single parameter measurement error distributions. Thirdly, Liu et al. (2018) used shrinkage methods to improve the finite-sample performance of the untransformed Poisson PCA estimators; similar ideas could be applied to improve the finite sample performance of the sequencing depth corrected estimators and the transformed estimators. Finally, a large area of research in PCA is the implementation of sparseness constraints. Since our method provides an unbiassed variance estimator for the Poisson means, it should be straightforward to implement a standard regularisation method to ensure consistency in the $p\gg n$ case.

Appendix A Proof of Consistency of Transformed Variance Estimator

Theorem** (3.1).**

If $e^{2\lambda}f(\lambda)^{2}=\sum_{n=0}^{\infty}\frac{h(n)\lambda^{n}}{n!}$ is a globally convergent Taylor series, $\Lambda$ is a non-negative random variable with finite raw moments $\mu_{n}$ and $\sum_{n=0}^{\infty}\frac{\mu_{n}|h_{n}|}{n!}$ converges, then the estimator (1) is consistent.

Proof.

Since the estimator (1) is unbiassed, the strong law of large numbers tells us that this is a consistent estimator provided ${\mathbb{E}}\left(g(\mathbf{X})\right)$ and ${\mathbb{E}}\left(\sum_{k=0}^{\mathbf{X}}(-1)^{\mathbf{X}-k}{\mathbf{X}\choose k}h(k)\right)$ are both finite. We have

[TABLE]

which must be finite since otherwise the variance of $f(\Lambda)$ which we are trying to estimate will not be finite.

We have that

[TABLE]

We are given that $\sum_{k=0}^{\infty}\frac{{\mathbb{E}}\left(\Lambda^{k}\right)}{k!}|h(k)|$ converges, so by Fubini’s theorem,

[TABLE]

is also finite. ∎

Appendix B Estimating Conditional Variance

We are interested in the logarithm case, so for large $\mathbf{X}$, we have $g(\mathbf{X})=\log(\mathbf{X})$. We now want to estimate $\operatorname{Var}(\log(\mathbf{X})|\Lambda)$. We will estimate it by taking the average of some function $\tilde{h}(X)$ for each $X$. That is, we want to design a function $\tilde{h}$ such that

[TABLE]

We will assume that $\mathbf{\Lambda}$ is relatively large. We note that

[TABLE]

For large $\mathbf{\Lambda}$, we have that $D=\frac{\mathbf{X}-\Lambda}{\Lambda}$ is almost certain to be much less than 1, so we can use the Taylor expansion

[TABLE]

We truncate this at the $T$th term, for some $T$, then we have

[TABLE]

We can truncate this again at $T$ (i.e. remove all powers of $D$ larger than $T$) to get

[TABLE]

Now we compute

[TABLE]

where $\genfrac{\{}{\}}{0.0pt}{}{k}{j}$ is a Stirling number of the second kind. The second last equation comes from setting $l=k-j$, and the last equation comes from setting $m=n-k$. We have that $\binom{n}{m}\genfrac{\{}{\}}{0.0pt}{}{n-m}{n-l-m}$ is the number of partitions of a set of size $n$ into $n-l$ blocks with $m$ chosen singleton blocks. This gives us that

[TABLE]

is the number of partitions of $n$ with $n-l$ blocks and no singletons, since for any partition into $n-l$ blocks, of which $M$ are singletons, this partition is included $\binom{M}{m}$ times in the $m$th term in the series (once for each subset of $m$ singletons). Therefore the total term assigned to this partition is

[TABLE]

We will denote the number of partitions with $n-l$ blocks and no singletons by $\genfrac{\{}{\}}{0.0pt}{}{n}{n-l}_{\textrm{sf}}$. We therefore have

[TABLE]

We now substitute this into Equation (2) (noting that for $n>0$, we have $\genfrac{\{}{\}}{0.0pt}{}{n}{0}_{\textrm{sf}}=0$, and for $k>\frac{n}{2}$, we have $\genfrac{\{}{\}}{0.0pt}{}{n}{k}_{\textrm{sf}}=0$) to get:

[TABLE]

Recall that $\frac{1}{(\mathbf{X}+1)\cdots(\mathbf{X}+l)}$ is an estimator for $\Lambda^{-l}$, so substituting in this estimator gives us the estimator $\tilde{h}(\mathbf{X})$ as

[TABLE]

where

[TABLE]

Our Poisson PCA package truncates at $T=9$. The terms $a_{l}$ can be precalculated once for all data points, so computation of $\tilde{h}$ is quick.

Appendix C Alternative Estimator for Conditional Variance

Recall that the general estimator for conditional variance of a transformation $g(X)$ depends on computing an estimator for $e^{-2\Lambda}\frac{\Lambda^{n}}{n!}$ for an integer $n$. The unbiased estimator for this quantity is $t_{n}=(-1)^{X-n}\binom{X}{n}$. However, this estimator can have very high variance, resulting in poor finite-sample performance. For the special case $g(X)=\log(X)$, we computed a useable approximation in Appendix B. For the general transformation case, we do not have a good general method. In this section, we suggest one possible approach that could be further studied.

The idea is to try to find two observations, $X_{1}$ and $X_{2}$ with the same value of $\mathbf{\Lambda}$. Then we have

[TABLE]

as an estimator for $e^{-2\lambda}\frac{(2\lambda)^{n}}{n!}$. We can then average this estimator over all pairs of values of $X_{1}$ and $X_{2}$ that come from approximately the same value of $\mathbf{\Lambda}$.

In practice we don’t know whether two values $X_{1}$ and $X_{2}$ come from the same value of $\Lambda$. We therefore weight all pairs by how plausible it is that they have the same value of $\mathbf{\Lambda}$. We use the weight derived from the likelihood ratio statistic:

[TABLE]

and take the weighted average of $t^{\prime}_{n}(X_{1},X_{2})$ divided by $2^{n}$ to get the estimator of $e^{-2\lambda}\frac{\lambda^{n}}{n!}$.

That is, our estimator is

[TABLE]