<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/quantifying-the-magnitude-of-a-population-decline-with-bayesian-time-series-modelling/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="http://www.seascapemodels.org/rstats/2017/06/18/estimating-popn-decline.html"> Bluecology blog</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<h2>Quantifying the magnitude of a population decline with Bayesian time-series modelling</h2>
<p>Population abundances tend to vary year to year. This variation can make<br>
it make it hard detect a change and hard to quantify exactly what that<br>
change is.</p>
<p>Bayesian time-series analysis can help us quantify a decline and put<br>
uncertainty bounds on it too. Here I will use the <a rel="nofollow" href="http://www.r-inla.org/" target="_blank">R-INLA<br>
package</a> to fit a time-series model to a<br>
population decline.</p>
<p><img alt="" src="https://i0.wp.com/www.seascapemodels.org/Images/estimating-popn-decline/noisy-abundance-trend-1.png?w=456" data-recalc-dims="1"></p>
<p>For instance, take the pictured time-series. Quantifying change as the<br>
difference between the first and last time-points is obviously<br>
misleading. Doing so would imply that abundance has declined by 77% from<br>
the historical value.</p>
<p>Another approach would be to compare the average of the first and last<br>
decades. Doing so would yield a 72% decline.</p>
<p>A better way might be to model the population trend over time and then<br>
estimate our change from the model. An advantage of a model is we can be<br>
more specific (and transparent) about assumptions. Using Bayesian<br>
time-series analysis we can also pick a model that is appropriate for<br>
the distribution of the data. e.g. here I will use <a rel="nofollow" href="https://en.wikipedia.org/wiki/Negative_binomial_distribution" target="_blank">negative<br>
binomial</a><br>
errors, because the abundances are counts.</p>
<h3>Simulating a time-series</h3>
<p>Let’s simulate a simple time-series of counts that declines at a<br>
constant (exponential) rate:</p>
<div>
<pre><code>set.seed(1969)
n  &lt;-  50
rho  &lt;-  0.8
prec  &lt;-  10
init &lt;- 7
acoef &lt;- c(rep(-0.01, n/2), rep(-0.03, n/2))
eta &lt;- init + (acoef*(1:n))
y &lt;- rnbinom(n, mu  =  exp(eta), size  =  5)
data  &lt;-  list(y =y, z = 1:n)
</code></pre>
</div>
<p>The results are pictured above. In short we calculate <code>eta</code> which is the<br>
linear predictor. Its exponent determines the mean of the negative<br>
binomial process, hence the exponential (not linear) decline. The<br>
<code>size = 5</code> means that the count data will be moderately over-dispersed –<br>
that is quite noisy (the variance will be > than the mean).</p>
<p>Notice that I put a bit of a ‘kink’ in eta by having the rate of decline<br>
(<code>acoef</code>), increase halfway through the time-series like this:</p>
<p><img alt="" src="https://i1.wp.com/www.seascapemodels.org/Images/estimating-popn-decline/expected-abundance-1.png?w=456" data-recalc-dims="1"></p>
<p>We can calculate the magnitude of change like this: <code>1 - (y[n]/y[1])</code> =<br>
0.77 which can be interpreted as what fraction of the intial years<br>
biomass has been lost?</p>
<p>What we want to do is put a smoother over time, so our estimate accounts<br>
for short-term variation. We also want to get credibility intervals (95%<br>
probability intervals) on the estimate of decline.</p>
<h3>Fitting a time-series model with INLA</h3>
<p>R-INLA is a pretty handy package. It is a very<br>
<a rel="nofollow" href="http://www.seascapemodels.org/rstats/2017/04/14/glmm-comparison.html" target="_blank">fast</a><br>
way to fit generalized linear models and can also handle a huge range of<br>
different types of <a rel="nofollow" href="http://www.r-inla.org/models/latent-models" target="_blank">random<br>
effects</a>. Here we will use<br>
an auto-regressive lag one process. That is saying that abundance at<br>
time t depends on abundance at time t-1 with some constant correlation<br>
coefficient.</p>
<p>We could also include other covariates in the model, for instance,<br>
abundance of a predator that might eat our organism of interest.<br>
However, we won’t go into that here.</p>
<p>First up, we specify the formula in INLA’s notation:</p>
<div>
<pre><code>library(INLA)
f1 &lt;- y ~ 1 + f(z, model = &quot;ar1&quot;)
</code></pre>
</div>
<p>Which just says model <code>y</code> as a function of an intercept term (constant<br>
mean) + an autoregressive process that depends on <code>z</code> (just an index for<br>
the year). The <code>f()</code> is INLA’s standard notation for modelling its<br>
different random effects models (AKA latent models), here we chose the<br>
<code>ar1</code> option.</p>
<p>Now the trick is to get INLA to compute the % difference between the<br>
<strong>expectation</strong> for the first and last time points as it goes. That way<br>
we not only get an estimate of the change but also we will get the full<br>
posterior distribution, so we can get our CIs. Note I put emphasis on<br>
<strong>expectation</strong>, because we won’t simply be calculating the difference<br>
between the first and last points (we did that above), but will in fact<br>
be calculating the difference between the model’s estimated mean values<br>
for the first and last points.</p>
<p>To do that we use INLA’s <a rel="nofollow" href="http://www.r-inla.org/faq#TOC-I-have-some-linear-combinations-of-the-nodes-in-the-latent-field-that-I-want-to-compute-the-posterior-marginal-of-is-that-possible-" target="_blank">linear<br>
combination</a><br>
feature.</p>
<p>This feature let’s us compute a linear combination of the estimated<br>
covariates. Here we will ask for the difference between the first and<br>
last values of <code>z</code> (the time-series). <code>inla</code> will multiply the <code>z</code><br>
values by the numbers we give below (here <code>-1</code> and <code>1</code>), then sum them<br>
together.</p>
<div>
<pre><code>lc1 &lt;- inla.make.lincomb(z = c(-1, rep(NA, n-2), 1))
names(lc1) = &quot;lc1&quot;
</code></pre>
</div>
<p>The row of <code>NA</code> in the middle just says to ignore the middle points.<br>
Note that we have also ignored the intercept. More on this in a moment.</p>
<p>Now we can fit the model, using negative binomial errors and specifying<br>
our linear combination:</p>
<div>
<pre><code>m1 &lt;- inla(f1,family = &quot;nbinomial&quot;, data  =  data,
    lincomb = lc1,
    control.inla = list(lincomb.derived.only = TRUE),
    control.predictor = list(compute = TRUE, link = 1)
)
</code></pre>
</div>
<p>And the results summary can be obtained (not shown here):</p>
<div>
<pre><code>summary(m1)
</code></pre>
</div>
<p>Note the large estimate for <code>rho</code>, which is the auto-correlation<br>
parameter. That happens because we have such a strong decline.</p>
<p>We could have also fit just a linear trend to our data, however, because<br>
the rate varies over time (and we wouldn’t know that beforehand if we<br>
had real data, not data we made up), it is nice to use <code>ar1</code> which has a<br>
kind of smoothing effect.</p>
<p>Also worth noting that modelling an trending (non-stationary)<br>
time-series with an <code>ar1</code> is not technically correct if we want to<br>
estimate rho. However, here we use <code>ar1</code> because it acts like a smoother<br>
in a GAM [e.g. <a rel="nofollow" href="https://www.r-bloggers.com/functional-anova-using-inla/" target="_blank">see<br>
here</a>.</p>
<h3>Examining the model’s predictions</h3>
<p>It is pretty straightforward to plot INLA’s predictions for the<br>
time-series. They are all stored in the <code>m1</code> object:</p>
<div>
<pre><code>plot(data$z, data$y, col = 'grey', type = 'l', lwd = 2, xlab = &quot;years&quot;, ylab = &quot;Abundance&quot;)
lines(data$z, m1$summary.fitted.values$mean, col = &quot;tomato&quot;, lwd = 2)
lines(data$z, m1$summary.fitted.values[,3], col = &quot;tomato&quot;, lwd = 1, lty = 2)
lines(data$z, m1$summary.fitted.values[,5], col = &quot;tomato&quot;, lwd = 1, lty = 2)
legend('topright', legend = c(&quot;observed&quot;, &quot;expected&quot;, &quot;95% CIs&quot;), lty = c(1,1,2),
       col = c(&quot;grey&quot;, &quot;tomato&quot;, &quot;tomato&quot;))
</code></pre>
</div>
<p><img alt="" src="https://i0.wp.com/www.seascapemodels.org/Images/estimating-popn-decline/fitted-abundance-1.png?w=456" data-recalc-dims="1"></p>
<p>Now we can extract the change in population size from our linear<br>
combination:</p>
<div>
<pre><code>m1$summary.lincomb.derived$mean

## [1] -1.081196
</code></pre>
</div>
<p>However, this is not quite the number we want. It is the change in <code>eta</code><br>
or the linear predictor. Negative binomial models have a log-link, so<br>
that is the change in log-space. We want the real change. If we<br>
calculate:</p>
<div>
<pre><code>100 * (1 - exp(m1$summary.lincomb.derived$mean))

## [1] 66.08103
</code></pre>
</div>
<p>we get the change as a % from the historical mean The reason being is<br>
that the linear combination is <code>log(eta[n]) - log(eta[1])</code>. Taking the<br>
exponent gets the change and one minus that is the % loss.</p>
<p>It is easy now to get the CIs on the linear combination:</p>
<div>
<pre><code>m1$summary.lincomb.derived

##     ID      mean       sd 0.025quant 0.5quant 0.975quant      mode kld
## lc1  1 -1.081196 0.309715  -1.695217 -1.08053 -0.4716486 -1.079157   0
</code></pre>
</div>
<p>We could also look at the marginal and plot the posterior distribution<br>
of the % loss:</p>
<div>
<pre><code>losses &lt;- 100 * (1 - exp(m1$marginals.lincomb.derived[[1]][,1]))
dens &lt;- m1$marginals.lincomb.derived[[1]][,2]
plot(losses, dens, type = 'l', xlim = c(0, 100),
     xlab = &quot;% loss&quot;, ylab = &quot;density&quot;, las = 1,
      main = &quot;Posterior density for % loss&quot;)
polygon(x = c(losses, rev(losses)),
        y = c(dens, rep(0, length(dens))), col = &quot;turquoise3&quot;)
</code></pre>
</div>
<p><img alt="" src="https://i1.wp.com/www.seascapemodels.org/Images/estimating-popn-decline/abundance-posterior-1.png?w=456" data-recalc-dims="1"></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://www.seascapemodels.org/rstats/2017/06/18/estimating-popn-decline.html"> Bluecology blog</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/GUMbTKljz1Q" height="1"></a></div>