13
$\begingroup$

I am working with a process in which I expect my variable to be Poisson distributed. For reasons that have to do with the scale, however, the values I obtain have a minimum of 11.

I have noted than when I estimate the fit of my data to a Poisson with an R function such as fitdistrplus::fitdist my data fits the best to a Poisson distribution (following the AIC values), but only when I subtract the minimum value of 11 to the vector of observations (the rationale for doing so is that values in a Poisson start from zero).

For example, the following gives a terrible fit to a Poisson:

obs <- c(15, 15, 13, 14, 13, 12, 12, 14, 14, 12, 11, 13, 11, 12, 12, 
         13, 14, 14, 13, 11, 12, 14, 17, 12, 14, 12, 12, 12, 12, 12, 
         13, 16, 15, 15, 17, 13, 13, 13, 13, 14, 14, 13, 15)
fitdist(obs, "pois")

But after subtracting the minimum value, the fit is better.

fitdist(obs-11, "pois")

However, this changes the lambda parameter of course, and in the Poisson I believe this changes as well the shape of the distribution. I am not only unsure about whether my transformation here is correct, but if it isn't why the fit improved.

Can anybody point me to a reference on this topic to know if my transformation is correct? i.e. whether a Poisson minus a constant still distributes as a Poisson?

EDIT: For clarification, the variable I am dealing with are rates of de novo mutation in the genome. These are rare events and usually modelled as Poisson processes. In the scale of number of mutations per haploid genome and generation, a few tens of mutations are observed per haploid genome and generation (the obs vector). This is why the minimum value is 11. However, note that there is an effect of the scale here. In the scale of diploid genomes values in obs would double, and in the scale or per nucleotide site (rather than per genome) the values in obs would all be so small (e.g. ~6*10e-8) that would round to zero.

$\endgroup$
7
  • 1
    $\begingroup$ "For reasons that have to do with the scale, however, the values I obtain have a minimum of 11" why? Is this due to some indirect measurement of counts (e.g. a voltage from a counter being hit by electrically charged particles)? $\endgroup$ Commented Jun 12 at 7:05
  • 3
    $\begingroup$ You need to explain a bit more what is going on. Currently this question creates three different answers: 1 interpreting your question as being about a truncated/sensored observation 2 interpreting your question as being about a shifted observation 3 interpreting yiur question being about the final sentence "whether a Poisson minus a constant still distributes as a Poisson?" $\endgroup$ Commented Jun 12 at 7:17
  • 1
    $\begingroup$ You say " I expect my variable to be Poisson distributed" but why is that? It is already not the case because of the cutoff at eleven. Another issue might be that whatever causes that shift or cut-off, is also completely distorting the entire measurement. $\endgroup$ Commented Jun 12 at 7:19
  • 3
    $\begingroup$ Your added context works a bit (it seems like you are not working with raw counts in some cases) but it is still confusing. "a few tens of mutations are observed per haploid genome and generation (the obs vector). This is why the minimum value is 11." I don't see the connection when you say 'this is why'. How come a few tens are observed? Is the number of counts multiplied by ten somehow (e.g. you made ten replicated observations of the 'same' genome giving that rough multiplicity by a factor ten)? Is your example representative, it is not a few tens, but instead ten plus a few singles? $\endgroup$ Commented Jun 12 at 8:07
  • 1
    $\begingroup$ If this is sampled from a bacterial population, you should not expect a Poisson distribution. $\endgroup$ Commented Jun 12 at 18:32

3 Answers 3

19
$\begingroup$

You should incorporate your estimate of the shift into your analysis

As others have pointed out, no, this is not a Poisson distribution (it is actually a shifted Poisson distribution). The bigger issue here (than what essentially amounts to a naming issue) is that you are "estimating" a lower bound of the support of your data by eyeballing the sample and then setting this lower bound as if it were an a priori fixed constant. That is a very bad way to estimate a distribution and it will tend to overstate your confidence in the outcome.

If you would like to fit your data to a shifted Poisson distribution, it would be better to leave the shift as an unknown parameter and estimate this from your data, incorporating the uncertainty in that estimator into your analysis. To do this using maximum likelihood estimation, you would use the log-likelihood function for shifted Poisson data. Taking $x_* \in \mathbb{N}_{0+}$ to be the shift parameter (giving the lower bound of the support) and $\lambda$ to be the rate parameter, your log-likelihood function would be:

$$\begin{align} \ell_\mathbf{x}(x_*, \lambda) &= \sum_{i=1}^n \log \text{ShiftPois}(x_i|x_*, \lambda) \\[6pt] &= \sum_{i=1}^n \log \text{Pois}(x_i-x_*| \lambda) \\[6pt] &= -n \lambda + \log(\lambda) \sum_{i=1}^n (x_i-x_*) - \sum_{i=1}^n \log \Gamma(x_i-x_*+1) \end{align}$$

for $x_* \leqslant x_{(1)}$ (with log-likelihood of negative infinity otherwise). For a given shift $x_* \leqslant x_{(1)}$ we get the rate estimator $\hat{\lambda}(x_*) = \sum_{i=1}^n (x_i-x_*) / n$, which yields the profile log-likelihood:

$$\begin{align} \ell_\mathbf{x}(x_*) &\equiv \ell_\mathbf{x}(x_*, \hat{\lambda}(x_*)) \\[12pt] &= - \sum_{i=1}^n (x_i-x_*) + \log \bigg( \frac{1}{n} \sum_{i=1}^n (x_i-x_*) \bigg) \sum_{i=1}^n (x_i-x_*) - \sum_{i=1}^n \log \Gamma(x_i-x_*+1) \\[6pt] &= \bigg[ \log \bigg( \frac{1}{n} \sum_{i=1}^n (x_i-x_*) \bigg) -1 \bigg] \sum_{i=1}^n (x_i-x_*) - \sum_{i=1}^n \log \Gamma(x_i-x_*+1) \\[6pt] &= \bigg[ \log \bigg( \sum_{i=1}^n (x_i-x_*) \bigg) - \log(n) -1 \bigg] \sum_{i=1}^n (x_i-x_*) - \sum_{i=1}^n \log \Gamma(x_i-x_*+1). \\[6pt] \end{align}$$

We can estimate the shift $x_*$ using the MLE $\hat{x}_* = \text{arg max}_{0 \leqslant r \leqslant x_{(1)}} \ell_\mathbf{x}(r)$, which can be computed as a simple maxima over a finite set. With some further analytical work we can derive the properties of these estimators, including their standard errors, etc.


Approximation of behaviour of estimators: The joint behaviour of the estimators in this problem is complicated by the fact that one parameter is continuous and one is discrete. As a simplifying approximation we can treat both parameters as if they were continuous and derive the resulting Fisher information to determine the (approximating) asymptotic standard error. If we were to treat the shift $x_*$ as a real parameter then we get the partial derivatives:

$$\begin{align} \frac{\partial \ell_\mathbf{x}}{\partial x_*}(x_*, \lambda) &= -n\log(\lambda) + \sum_{i=1}^n \psi (x_i-x_*+1), \\[12pt] \frac{\partial \ell_\mathbf{x}}{\partial \lambda}(x_*, \lambda) &= - n + \frac{1}{\lambda} \sum_{i=1}^n (x_i-x_*), \\[12pt] \frac{\partial^2 \ell_\mathbf{x}}{\partial x_*^2}(x_*, \lambda) &= - \sum_{i=1}^n \psi' (x_i-x_*+1), \\[12pt] \frac{\partial^2 \ell_\mathbf{x}}{\partial x_* \partial \lambda}(x_*, \lambda) &= - \frac{n}{\lambda}, \\[12pt] \frac{\partial^2 \ell_\mathbf{x}}{\partial \lambda^2}(x_*, \lambda) &= - \frac{1}{\lambda^2} \sum_{i=1}^n (x_i-x_*), \\[12pt] \end{align}$$

where $\psi$ is the digamma function and $\psi'$ is the trigamma function. The second-order partial derivatives can be encapsulated in the Hessian matrix:

$$\begin{align} H(x_*, \lambda) &= - \frac{1}{\lambda} \begin{bmatrix} \lambda \sum_{i=1}^n \psi' (x_i-x_*+1) & n \\[6pt] n & \lambda^{-1} \sum_{i=1}^n (x_i-x_*) \\[6pt] \end{bmatrix}. \end{align}$$

Taking the expected value of the random variables under the shifted Poisson model and defining the function $h(\lambda) \equiv \mathbb{E}(\psi' (X+1) | X \sim \text{Pois}(\lambda))$ (and noting that $\mathbb{E}(X_i-x_*) = \lambda$) then gives the Fisher information:

$$\begin{align} \mathcal{I}(x_*, \lambda) &= \frac{n}{\lambda} \begin{bmatrix} \lambda h(\lambda) & 1 \\[6pt] 1 & 1 \\[6pt] \end{bmatrix}. \end{align}$$

Inverting this matrix gives the (approximate) asymptotic variance of the estimators:

$$\hat{\mathbb{V}}(\hat{x}_*, \hat{\lambda}) = \mathcal{I}(\hat{x}_*, \hat{\lambda})^{-1} = \frac{1}{n} \cdot \frac{\hat{\lambda}}{\hat{\lambda} h(\hat{\lambda}) - 1} \begin{bmatrix} 1 & -1 \\[6pt] -1 & \hat{\lambda} h(\hat{\lambda}) \\[6pt] \end{bmatrix}.$$

This should give you a reasonable approximation to the asymptotic variance for large $n$. The approximation involves treating the discrete estimator for the shift (lower bound of the support) as if it were an estimator for a continuous parameter.


Implementation in R: We can implement this estimator in R using the function below. This function has been programmed to compute the MLE of the shifted Poisson model and return a fitdist object of the form used in the package you are using.

fitdist.shiftpois <- function(data) {
  
  #Check input data
  if(!is.vector(data))              stop('Error: Input data should be a vector')
  if(!is.numeric(data))             stop('Error: Input data should be a numeric vector')
  if(any(data != as.integer(data))) stop('Error: Input data should contain integers')
  if(min(data) < 0)                 stop('Error: Input data should be non-negative')

  #Set preliminary values and profile log-likelihood
  PAR <- c('shift', 'lambda')
  n   <- length(data)
  MIN <- min(data)
  LL  <- rep(-Inf, MIN+1)
  for (r in 0:MIN) { 
     RATE <- mean(data-r)
     LL[r+1] <- sum(dpois(data-r, lambda = RATE, log = TRUE)) }
  
  #Estimate the parameters
  MAXLOGLIKE <- max(LL)
  EST.XSTAR  <- which.max(LL)-1
  EST.LAMBDA <- mean(data-EST.XSTAR)
  ESTIMATE   <- c(EST.XSTAR, EST.LAMBDA)
  names(ESTIMATE) <- PAR
  
  #Estimate the asymptotic variance matrix
  HH <- function(lambda) { 
    UPPER <- ceiling(lambda + 5*sqrt(lambda))   #Five SD above mean
    DD    <- dpois(0:UPPER, lambda = lambda)
    TT    <- trigamma(0:UPPER+1)
    sum(TT*DD) }
  EST.HH <- HH(EST.LAMBDA)
  VAR <- matrix(0, nrow = 2, ncol = 2)
  rownames(VAR) <- PAR
  colnames(VAR) <- PAR
  CONST <- (EST.LAMBDA/(EST.LAMBDA*EST.HH-1))/n
  VAR[1, 1] <- CONST
  VAR[1, 2] <- -CONST
  VAR[2, 1] <- -CONST
  VAR[2, 2] <- CONST*EST.LAMBDA*EST.HH
  SD <- c(sqrt(VAR[1, 1]), sqrt(VAR[2, 2]))
  names(SD) <- PAR
  COR <- VAR[1, 2, drop = FALSE]/sqrt(VAR[1, 1]*VAR[2, 2])
  
  #Generate the output
  OUT <- list(estimate = ESTIMATE, method = 'mle', 
              sd = SD, cor = COR, vcov = VAR, loglik = MAXLOGLIKE,
              aic = 4-2*MAXLOGLIKE, bic = 2*log(n)-2*MAXLOGLIKE,
              n = n, data = data, distname = 'shiftpois',
              fix.arg = NULL, fix.arg.fun = NULL, dots = NULL,
              convergence = 0, discrete = TRUE, weights = NULL)
  class(OUT) <- 'fitdist'
  
  #Return the output
  OUT }

Using your data we obtain the fitted distribution object:

#Fit the shifted Poisson distribution
(FIT <- fitdist.shiftpois(obs))
    
Fitting of the distribution ' shiftpois ' by maximum likelihood 
Parameters:
        estimate Std. Error
shift   11.00000  0.5337348
lambda   2.27907  0.5812697
$\endgroup$
2
  • 5
    $\begingroup$ Sometimes this shift parameter is a property of a measurement apparatus and can be reasonably estimated from multiple measurements combined, or by using some initial calibration. In such a case, adding it as a parameter to the analysis will make the outcomes of estimates more variable and potentially biased (e.g. smaller samples may get estimated with a larger shift). $\endgroup$ Commented Jun 12 at 7:11
  • 3
    $\begingroup$ This is an excellent answer, thank you very much. Actually, the assumption that the shifting variable is continuous is correct. In the question, I was giving the values of "obs" and the constant 11 as rounded to integers under the assumption that this is the correct way to work with a Poisson, since it is a discrete distribution. $\endgroup$ Commented Jun 12 at 8:12
17
$\begingroup$

Is a Poisson minus a constant still a Poisson?

No. To begin with, the support of a Poisson is the non-negative integers.

But you didn't subtract a constant from a Poisson, you subtracted a constant from the data and compared that to a Poisson; that corresponds to adding a constant to the Poisson-distributed variable.

If you move along by a positive constant, it's no longer got support on (and only on) all of the non-negative integers.

There's more to it than that, but that's sufficient to make it not-Poisson.

However, if you add or subtract a constant from a Poisson-distributed variable, the distribution is that of a shifted Poisson.

[It looks like you may have chosen your shift with reference to the data; if that was the case, its then an estimated parameter (albeit perhaps inefficiently), not a known constant. Beware doing one but treating it like the other. If that's not the case, and the 11 came from considerations external to the data, then ignore this.]

You ask about whether subtracting 11 and fitting a Poisson is "correct" but I am not sure what you're asking for there. It may in some cases be a pretty reasonable activity but I can't tell what led you to this choice.

Your actual variable almost certainly not actually a shifted-Poisson, though the approximation looks to be fairly reasonable. Is it even a count-variable? Can you describe more about what the variable is measuring/counting?

Why do you need to choose a distributional model?


Further response after seeing some edits to the question and reading some of your comments:

You now add crucial information that you should have begun with; omitting that information left us answering questions that are of little use.

You claim that your original (unrounded/unbinned) variable is actually continuous. It isn't. It's scaled count, and hence still strictly discrete. If I take a set of integers and divide each one by a different set of integers, the set of possible outcomes is still countable (and duplication of observations - ties - has positive probability).

You should certainly not round values. However, neither should you scale for the exposures (those scaling constants you divided by). You should keep the counts and the exposures separated; you use them both but not like this. I will take it that your assumption of a constant rate per unit of the scaling factor is reasonable, but the thing you avoided using - per nucleotide - might perhaps be better (however, I'll leave that matter to your judgement on the basis of scientific knowledge you should know better than I do). It doesn't matter that the nucleotide counts are very large, because you aren't scaling the data and you aren't rounding.

The problem with scaling the counts is that - while it gets the means on a comparable basis - you screw up the variances.

A better way to deal with such scaled counts is to not scale them, but to use the scale factor as an exposure type measure. That is, to scale the Poisson rate in the model, not scaling the observed count itself. This way you avoid screwing up the variance.

If the original counts are $Y_1, Y_2, ...$ and the exposures (what you scale by) for each observation are $h_1$, $h_2$,... then rather than transform $Z_i = Y_i / h_i$ to make it so that $E(Z_i)$ are similar, you would instead write the model as $Y_i\sim \text{Pois}(h_i\cdot\mu)$.

An overall estimate of $\mu$ (across observations that have the same rate $\mu$) can then be attained by adding all the counts and all the exposures that should have the same $\mu$; $\hat{\mu} = \sum_i Y_i/\sum_i h_i$.

If you have predictors (IVs/covariates) you would usually use a Poisson GLM; with a log-link the logs of the $h_i$ become an offset term in that model.

If you need to compare rates across two groups of values (which are assumed to have the same rate within group but not necessarily across it) that can be done using standard methods (e.g. poisson.test in R). By conditioning on the total count it essentially converts into a standard test of proportions. Multiple groups just become a test of multinomial proportions known under $H_0$ (they are obtained from the known exposure sets), and so you don't necessarily need a specialized routine. More complicated situations are typically handled with standard tests in the GLM framework.

[It might well be that there is heterogeneity in rates (due to factors other than your scaling factor); if some factors can be accounted for they may be modelled, but if there are substantial sources of unmeasured heterogeneity you may end up with overdispersion. Then a simple Poisson model may not be suitable.]

$\endgroup$
3
  • $\begingroup$ Thank you for the useful reply, I added more context to my question in a recent edit. Answering the question about why I decided to subtract 11 from the data: I noticed that the shape of my distribution was "about right" for a Poisson even when the fit of "obs" was terrible, and the reason for that was that with a mean higher than 11 the shape of the Poisson changes in a way that does not fit the data. $\endgroup$ Commented Jun 12 at 8:02
  • $\begingroup$ Knowing that the Poisson is usually used for modelling processed that are rare in frequency, it was intuitive to me that the minimum value of "obs" should be zero. This way, I am really assuming that all observations have a minimum value of 11 and the variation on top of that is a Poisson starting from zero. Doing so the fit was much better, and agrees with my expectation for this process. But as you and others mention, this seems not to be strictly a Poisson anymore $\endgroup$ Commented Jun 12 at 8:04
  • $\begingroup$ I made some edits to my answer. It's not been made clear what you're using your model to do, so I addressed both parameter estimation and basic testing, but of course there may well be other things you're actually trying to do. $\endgroup$
    – Glen_b
    Commented Jun 12 at 14:50
13
$\begingroup$

You know that with a Poisson, the mean equals the variance. Adding a constant to a random variable does not alter the variance but changes the mean. Therefore, a Poisson minus a constant cannot still be Poisson-distributed.

Your problem sounds like either a left-truncated or a left-censored Poisson. If observations $\le 10$ are missing entirely, that is truncated. If they are merely recoded to 11, that is censored.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.