<p><div><div><div><div data-show-faces="false" data-height="21" data-width="100" data-layout="button_count" data-send="true" data-href="https://www.r-bloggers.com/truncated-poisson-distributions-in-r-and-stan-by-ellis2013nz/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="http://ellisp.github.io/blog/2018/03/20/truncated-poisson"> Peter's stats stuff - R</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<h2>Data that’s been only partially observed</h2>
<p>I’ve been updating my skills in fitting models to truncated data and was pleased to find that, like so much else in statistics, it’s easier than it used to be.</p>
<p>First, some definitions:</p>
<ul>
<li><em>censored</em> data is where some observations are cut off at some maximum or minimum level; those data points are “off the scale” but at least we know they exist and we know which direction they are off the scale. For example, if we were analysing the life span of people born in 1980, anyone who has survived to the end of 2017 has their age at death recorded as “37 or higher”. We know they’re in the data, and their value is at least some minimum amount, but we don’t know with precision what it will end up being.</li>
<li><em>truncated</em> data is where data beyond some maximum or minimum level is just missing. Typically this is because of some feature of the measurement process eg anything smaller than X just doesn’t show up.</li>
</ul>
<p>I’ve got some future blog posts on a more substantive real life issue where I have count data for which, in some situations, I only see the observations with a count higher than some threshold. Let’s imagine, for example, we are looking at deaths per vehicle crash, and are dependent for measurement on newspapers that only report crashes with at least two deaths, even though many crashes have one or zero deaths.</p>
<p>Here’s a greatly simplified example. I generate 1,000 observations of counts, with an average value of 1.3. Then I compare that original distribution with what I’d get if only those of two or higher were observed. It looks like this:</p>
<p><img width="100%" src="http://ellisp.github.io/img/0120-truncated-poisson.svg"></p>
<p>…generated by this code:</p>
<figure>
<pre><code data-lang="r"><span>library</span><span>(</span><span>tidyverse</span><span>)</span><span>
</span><span>library</span><span>(</span><span>scales</span><span>)</span><span>
</span><span>library</span><span>(</span><span>fitdistrplus</span><span>)</span><span>
</span><span>library</span><span>(</span><span>rstan</span><span>)</span><span>
</span><span>library</span><span>(</span><span>truncdist</span><span>)</span><span>
</span><span># original data:
</span><span>set.seed</span><span>(</span><span>321</span><span>)</span><span>
</span><span>a</span><span> </span><span>&lt;-</span><span> </span><span>rpois</span><span>(</span><span>1000</span><span>,</span><span> </span><span>1.3</span><span>)</span><span>
</span><span># truncated version of data:
</span><span>b</span><span> </span><span>&lt;-</span><span> </span><span>a</span><span>[</span><span> </span><span>a</span><span> </span><span>></span><span> </span><span>1</span><span>]</span><span>
</span><span># graphic:
</span><span>data_frame</span><span>(</span><span>value</span><span> </span><span>=</span><span> </span><span>c</span><span>(</span><span>a</span><span>,</span><span> </span><span>b</span><span>),</span><span>
</span><span>variable</span><span> </span><span>=</span><span> </span><span>rep</span><span>(</span><span>c</span><span>(</span><span>&quot;Original data&quot;</span><span>,</span><span> </span><span>&quot;Truncated so only observations of 2 or more show up&quot;</span><span>),</span><span> </span><span>c</span><span>(</span><span>length</span><span>(</span><span>a</span><span>),</span><span> </span><span>length</span><span>(</span><span>b</span><span>))))</span><span> </span><span>%>%</span><span>
</span><span>ggplot</span><span>(</span><span>aes</span><span>(</span><span>x</span><span> </span><span>=</span><span> </span><span>value</span><span>))</span><span> </span><span>+</span><span>
</span><span>geom_histogram</span><span>(</span><span>binwidth</span><span> </span><span>=</span><span> </span><span>1</span><span>,</span><span> </span><span>colour</span><span> </span><span>=</span><span> </span><span>&quot;white&quot;</span><span>)</span><span> </span><span>+</span><span>
</span><span>facet_wrap</span><span>(</span><span>~</span><span>variable</span><span>,</span><span> </span><span>ncol</span><span> </span><span>=</span><span> </span><span>1</span><span>)</span><span> </span><span>+</span><span>
</span><span>ggtitle</span><span>(</span><span>&quot;Comparing full and truncated datasets from a Poisson distribution&quot;</span><span>)</span><span> </span><span>+</span><span>
</span><span>labs</span><span>(</span><span>y</span><span> </span><span>=</span><span> </span><span>&quot;Number of observations&quot;</span><span>)</span><span>
</span><span># fitting a model to original works well:
</span><span>mean</span><span>(</span><span>a</span><span>)</span><span>
</span><span>fitdistr</span><span>(</span><span>a</span><span>,</span><span> </span><span>&quot;Poisson&quot;</span><span>)</span><span>
</span><span># but obviously not if naively done to the truncated version:
</span><span>mean</span><span>(</span><span>b</span><span>)</span><span>
</span><span>fitdistr</span><span>(</span><span>b</span><span>,</span><span> </span><span>&quot;Poisson&quot;</span><span>)</span></code></pre>
</figure>
<p>Estimating the key parameter <code>lambda</code> for the full data (<code>a</code>) works well, giving an estimate of 1.347 that is just over one standard error from the true value of 1.3. The <code>fitdistr</code> function from the <code>MASS</code> package distributed with base R does a nice job in such circumstances.</p>
<p>But the mean value of <code>b</code> is badly biased upwards if used to estimate <code>lambda</code>; at 2.6 the mean of <code>b</code> is roughly twice the correct value of the mean of the underlying distribution. Obviously, removing a whole bunch of data at one end of the distribution is going to make naive estimation methods biased. So we need specialist methods that try to estimate lambda on the assumption that the data come from a Poisson distribution, but only the right-most part of it.</p>
<h2>Maximum likelihood</h2>
<p>The <code>fitdistrplus</code> package by Aurélie Siberchicot, Marie Laure Delignette-Muller and Christophe Dutang in combination with <code>truncdist</code> by Frederick Novomestky and Saralees Nadarajah gives a straightforward way to implement maximum likelihood estimation of a truncated distribution. Methods other than maximum likelihood are also available if required.</p>
<p>You need to make truncated versions of the <code>dpois</code> and <code>ppois</code> functions (or their equivalents for whatever distribution you are modelling) and use these within <code>fitdistrplus::fitdist</code>, which has some added functionality over <code>MASS::fitdistr</code> used in the previous chunk of code.</p>
<figure>
<pre><code data-lang="r"><span>#-------------MLE fitting in R-------------------
# adapted from https://stackoverflow.com/questions/16947799/fitting-a-lognormal-distribution-to-truncated-data-in-r
</span><span>dtruncated_poisson</span><span> </span><span>&lt;-</span><span> </span><span>function</span><span>(</span><span>x</span><span>,</span><span> </span><span>lambda</span><span>)</span><span> </span><span>{</span><span>
</span><span>dtrunc</span><span>(</span><span>x</span><span>,</span><span> </span><span>&quot;pois&quot;</span><span>,</span><span> </span><span>a</span><span> </span><span>=</span><span> </span><span>1.5</span><span>,</span><span> </span><span>b</span><span> </span><span>=</span><span> </span><span>Inf</span><span>,</span><span> </span><span>lambda</span><span>)</span><span>
</span><span>}</span><span>
</span><span>ptruncated_poisson</span><span> </span><span>&lt;-</span><span> </span><span>function</span><span>(</span><span>x</span><span>,</span><span> </span><span>lambda</span><span>)</span><span> </span><span>{</span><span>
</span><span>ptrunc</span><span>(</span><span>x</span><span>,</span><span> </span><span>&quot;pois&quot;</span><span>,</span><span> </span><span>a</span><span> </span><span>=</span><span> </span><span>1.5</span><span>,</span><span> </span><span>b</span><span> </span><span>=</span><span> </span><span>Inf</span><span>,</span><span> </span><span>lambda</span><span>)</span><span>
</span><span>}</span><span>
</span><span>fitdist</span><span>(</span><span>b</span><span>,</span><span> </span><span>&quot;truncated_poisson&quot;</span><span>,</span><span> </span><span>start</span><span> </span><span>=</span><span> </span><span>c</span><span>(</span><span>lambda</span><span> </span><span>=</span><span> </span><span>0.5</span><span>))</span><span> </span></code></pre>
</figure>
<p>Note that to do this I specified the lower threshold as 1.5; as all the data are integers this effectively means we only observe the observations of 2 or more, as is the case. We also needed to specify a reasonably starting value for the estimate of <code>lambda</code> – getting this too far out will lead to errors.</p>
<p>This method gives us an estimate of 1.34 with a standard error of 0.08, which is pretty good given we’ve only got 398 observations now. Of course, we’ve got the luxury of knowing for sure the true data generating process was Poisson.</p>
<h2>Bayesian</h2>
<p>For an alternative Bayesian method, Stan makes it easy to describe data and probability distributions as truncated. The Stan manual has an entire chapter on truncated or censored data. Here’s an example Stan program to estimate the mean of the original Poisson distribution from our truncated data. As well as the original data, which I call <code>x</code> in this program, we need to tell it how many observations (<code>n</code>), the <code>lower_limit</code> that it was truncated by, and whatever is needed to characterise a prior distribution for the parameter we’re estimating.</p>
<p>The key bits of the program below are:</p>
<ul>
<li>In the <code>data</code> chunk, specify that the data for <code>x</code> has a lower limit of <code>lower_limit</code></li>
<li>In the <code>model</code> chunk, specify that distribution of <code>x</code> is truncated via <code>T[lower_limit, ]</code></li>
</ul>
<figure>
<pre><code data-lang="stan">data {
  int n;
  int lower_limit;
  int <null lower="lower_limit"> x[n];
real lambda_start_mu;
real lambda_start_sigma;
}
parameters {
real<null lower="0">lambda;
}
model {
lambda ~ normal(lambda_start_mu, lambda_start_sigma);
for(i in 1:n){
x[i] ~ poisson(lambda) T[lower_limit, ];
}
}

<p>With a little more effort it’s possible to extend this by making Stan estimate <code>lower_limit</code> from the data; not necessary in my hypothetical example because I know where the minimum cut-off point of observed data lies.</p>
<p>Here’s how the data are fed to Stan from R:</p>
<figure>
<pre><code data-lang="r"><span>#-------------Calling Stan from R--------------
</span><span>data</span><span> </span><span>&lt;-</span><span> </span><span>list</span><span>(</span><span>
</span><span>x</span><span> </span><span>=</span><span> </span><span>b</span><span>,</span><span>
</span><span>lower_limit</span><span> </span><span>=</span><span> </span><span>2</span><span>,</span><span>
</span><span>n</span><span> </span><span>=</span><span> </span><span>length</span><span>(</span><span>b</span><span>),</span><span>
</span><span>lambda_start_mu</span><span> </span><span>=</span><span> </span><span>2</span><span>,</span><span>
</span><span>lambda_start_sigma</span><span> </span><span>=</span><span> </span><span>1</span><span>
</span><span>)</span><span>
</span><span>fit</span><span> </span><span>&lt;-</span><span> </span><span>stan</span><span>(</span><span>&quot;0120-trunc.stan&quot;</span><span>,</span><span> </span><span>data</span><span> </span><span>=</span><span> </span><span>data</span><span>,</span><span> </span><span>cores</span><span> </span><span>=</span><span> </span><span>4</span><span>)</span><span>
</span><span>fit</span><span>
</span><span>plot</span><span>(</span><span>fit</span><span>)</span><span> </span><span>+</span><span> </span><span>ggtitle</span><span>(</span><span>&quot;Credibility interval for lambda, estimated by Stan from truncated data&quot;</span><span>,</span><span>
</span><span>&quot;(Correct value is 1.3)&quot;</span><span>)</span><span> </span><span>+</span><span>
</span><span>labs</span><span>(</span><span>y</span><span> </span><span>=</span><span> </span><span>&quot;Estimated parameters&quot;</span><span>)</span><span> </span><span>+</span><span>
</span><span>theme_minimal</span><span>(</span><span>base_family</span><span> </span><span>=</span><span> </span><span>&quot;myfont&quot;</span><span>)</span></code></pre>
</figure>
<p>This gives us a posterior distribution for <code>lambda</code> that matches that from the <code>fitdistrplus</code> method: 1.35 with a standard deviation of 0.08. The <code>rstan</code> package automatically turns this into a <code>ggplot2</code> image of a credibility interval:</p>
<p><img width="100%" src="http://ellisp.github.io/img/0120-stan-fit.svg"></p>
<p>So, nice. Two simple ways to estimate the original distribution from truncated data.</p>
	<p>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="http://ellisp.github.io/blog/2018/03/20/truncated-poisson"> Peter's stats stuff - R</a></strong>.</div>
<hr>
<a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="https://feedburner.google.com/fb/a/mailverify?uri=RBloggers">daily e-mail updates</a></strong> about <a rel="nofollow" href="https://www.r-project.org/" title="The R Project for Statistical Computing">R</a> news and <a rel="nofollow" href="https://www.r-bloggers.com/search/tutorial" title="R tutorials">tutorials</a> on topics such as: <a rel="nofollow" href="https://www.r-bloggers.com/search/data%20science" title="Data science">Data science</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Big%20Data" title="Big Data">Big Data, <a rel="nofollow" href="https://www.r-users.com/" title="R jobs">R jobs</a>, visualization (<a rel="nofollow" href="https://www.r-bloggers.com/search/ggplot2" title="ggplot and ggplot2 tutorials">ggplot2</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/boxplot" title="Boxplots using lattice and ggplot2 tutorials">Boxplots</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/map" title="Maps and gis">maps</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/animation" title="Animation in R">animation</a>), programming (<a rel="nofollow" href="https://www.r-bloggers.com/search/RStudio" title="RStudio IDE for R">RStudio</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/sweave" title="Sweave and literate programming">Sweave</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX" title="LaTeX in R">LaTeX</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/SQL" title="SQL and databases">SQL</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/eclipse" title="Eclipse IDE for R">Eclipse</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/git" title="git and github, Version Control System">git</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/hadoop" title="Large data in R using Hadoop">hadoop</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping" title="Web Scraping of google, facebook, yahoo, twitter and more using R">Web Scraping</a>) statistics (<a rel="nofollow" href="https://www.r-bloggers.com/search/regression" title="Regressions and ANOVA analysis tutorials">regression</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/PCA" title="principal component analysis tutorial">PCA</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/time+series" title="Time series">time series</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/trading" title="finance trading">trading</a>) and more...
<img width="1" alt="" src="http://feeds.feedburner.com/~r/RBloggers/~4/b06-kiNu0co" height="1"></a></div></null></null></code></pre></figure>