<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/smoothing-a-time-series-with-a-bayesian-model/"></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/21/bayesian-smoothing.html"> Bluecology blog</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<h2>Smoothing a time-series with a Bayesian model</h2>
<p>Recently I looked at <a rel="nofollow" href="http://www.seascapemodels.org/rstats/2017/06/18/estimating-popn-decline.html" target="_blank">fitting a smoother to a<br>
time-series</a> using<br>
Bayesian modelling.<br>
Now I will look at how you can control the smoothness by using more or<br>
less informative priors on the precision (1/variance) of the random<br>
effect.</p>
<p>We will use the <a rel="nofollow" href="http://www.seascapemodels.org/rstats/2017/06/18/estimating-popn-decline.html" target="_blank">same dataset as the last<br>
post</a>.</p>
<p>To control the priors for an <a rel="nofollow" href="http://www.r-inla.org/" target="_blank">R-INLA</a> model, we<br>
use the <code>hyper</code> argument (not hyperactive, but hyper-parameters):</p>
<div>
<pre><code>library(INLA)
f3 &lt;- y ~ 1 + f(z, model = &quot;rw1&quot;, scale.model = TRUE,
hyper = list(theta = list(prior=&quot;pc.prec&quot;, param=c(1,0.01))))
</code></pre>
</div>
<p>We can control the level of smoothing through <code>param=c(theta1,0.01)</code>. A<br>
value of 1 (theta1) is a reasonable starting point (<a rel="nofollow" href="http://www.math.ntnu.no/inla/r-inla.org/doc/latent/rw1.pdf" target="_blank">based on the INLA<br>
documentation</a>).<br>
Lower values will result in a smoother fit.</p>
<p>The <code>pc.param</code> stands for Penalized complexity parameters (you could<br>
also use a <code>loggamma</code> prior here). My understanding of <a rel="nofollow" href="https://arxiv.org/abs/1403.4630" target="_blank">penalized<br>
complexity priors</a> is that they shrink<br>
the parameter estimate towards a ‘base-model’ that is less flexible. In<br>
this case, we are shrinking the standard deviation (AKA the flexibility)<br>
of the random walk (ie how much sequential data points deviate from each<br>
other) towards zero. Ultimately if we set theta1 near zero the smoother<br>
will be a straight line.</p>
<p>So let’s fit three models with theta1 varying and see how it affects the<br>
smoothness of the fit:</p>
<div>
<pre><code>f1 &lt;- y ~ 1 + f(z, model = &quot;rw1&quot;, scale.model = TRUE,
hyper = list(theta = list(prior=&quot;pc.prec&quot;, param=c(0.01,0.01))))

f2 &lt;- y ~ 1 + f(z, model = &quot;rw1&quot;, scale.model = TRUE,
hyper = list(theta = list(prior=&quot;pc.prec&quot;, param=c(0.05,0.01))))

f3 &lt;- y ~ 1 + f(z, model = &quot;rw1&quot;, scale.model = TRUE,
hyper = list(theta = list(prior=&quot;pc.prec&quot;, param=c(1,0.01))))

m1 &lt;- inla(f1,family = &quot;nbinomial&quot;, data  =  dat,
    control.predictor = list(compute = TRUE, link = 1)
)

m2 &lt;- inla(f2,family = &quot;nbinomial&quot;, data  =  dat,
    control.predictor = list(compute = TRUE, link = 1)
)

m3 &lt;- inla(f3,family = &quot;nbinomial&quot;, data  =  dat,
    control.predictor = list(compute = TRUE, link = 1)
)
</code></pre>
</div>
<p>Here are the resulting fits:</p>
<div>
<pre><code>plot(dat$z, dat$y, col = 'grey80', type = 'l', lwd = 1, xlab = &quot;years&quot;, ylab = &quot;Abundance&quot;, las = 1)

lines(dat$z, m1$summary.fitted.values$mean, col = &quot;skyblue&quot;, lwd = 2)
lines(dat$z, m2$summary.fitted.values$mean, col = &quot;steelblue&quot;, lwd = 2)
lines(dat$z, m3$summary.fitted.values$mean, col = &quot;darkblue&quot;, lwd = 2)

legend('topright', legend = c(&quot;data&quot;, &quot;theta1 = 0.01&quot;, &quot;theta1 = 0.05&quot;, &quot;theta1 = 1&quot;),
      lty = 1, col = c(&quot;grey80&quot;, &quot;skyblue&quot;, &quot;steelblue&quot;, &quot;darkblue&quot;), lwd = 2)
</code></pre>
</div>
<p><img alt="" src="https://i1.wp.com/www.seascapemodels.org/Images/bayesian-smoothing_files/fitted-models-vary-theta-1.png?w=456" data-recalc-dims="1"></p>
<p>The exact value to use for theta1 will vary depending on our data. If<br>
the data are more informative (e.g. a longer time-series) we will have<br>
to use a smaller value to create a smoother fit.</p>
<div>
<pre><code>dat2 &lt;- list(y = dat$y[1:40], z = 1:40)
m4 &lt;- inla(f2,family = &quot;nbinomial&quot;, data  =  dat2,
    control.predictor = list(compute = TRUE, link = 1)
)
plot(dat$z, dat$y, col = 'grey80', type = 'l', lwd = 1, xlab = &quot;years&quot;, ylab = &quot;Abundance&quot;, las = 1)

lines(dat$z, m2$summary.fitted.values$mean, col = &quot;skyblue&quot;, lwd = 2)
lines(dat2$z, m4$summary.fitted.values$mean, col = &quot;steelblue&quot;, lwd = 2)

legend('topright', legend = c(&quot;data&quot;, &quot;theta1 = 0.05, 50 data points&quot;, &quot;theta1 = 0.05, 40 data points&quot;),
      lty = 1, col = c(&quot;grey80&quot;, &quot;skyblue&quot;, &quot;steelblue&quot;), lwd = 2)
</code></pre>
</div>
<p><img alt="" src="https://i1.wp.com/www.seascapemodels.org/Images/bayesian-smoothing_files/fitted-models-vary-n-1.png?w=456" data-recalc-dims="1"></p>
<p>And that is pretty much it. I haven’t read how to choose the ‘optimal’<br>
value of theta, but presumably one could do it with cross validation or<br>
perhaps a model selection measure. Message me on<br>
<a rel="nofollow" href="https://twitter.com/bluecology/" target="_blank">Twitter</a> if you have seen an example<br>
of doing this.</p>
<p>It is amazing to see how the use of priors in statistics has changed<br>
since I first learned about Bayesian statistics. It used to be that you<br>
would use informative priors only if you had strong prior beliefs about<br>
a parameter estimate, such as from expert elicitation, or repeat<br>
experiments. If you didn’t have strong prior beleifs, then the view (at<br>
least amongst many ecologists) was that it was most cautious to use a<br>
highly uninformative prior like the good old Winbugs gamma(0.01, 0.01)<br>
prior that was used for precision parameters.</p>
<p>Now the experts are encouraging us to use weakly informative priors,<br>
even when little prior evidence is available. The case being that too<br>
broad a prior can slow computations and result in ridiculous results.<br>
Consider the broad gamma(0.01, 0.01) prior: it amounts to giving equal<br>
weight to a standard deviation of 1 as it does to an SD of 10000. The<br>
end result is that this ‘uninformative prior’ can bias your estimates of<br>
the SD to be too high.</p>
<p>As demonstrated here, another nice feature of informative priors is they<br>
can be used to control ‘shrinkage’. Here by varying theta1 we could<br>
shrink our model fit towards a base case model (e.g. a null hypothesis)<br>
that had no temporal variation.</p>
<p>If you are interested to learn more, it is worth reading at least the<br>
Introduction of <a rel="nofollow" href="https://arxiv.org/abs/1403.4630" target="_blank">Simpson et al. Penalized Complexity Priors<br>
pub</a>. Other good general guidance can<br>
be found on the <a rel="nofollow" href="https://github.com/stan-dev/stan/wiki/Prior-Choice-Recommendations" target="_blank">STAN<br>
wiki</a>.</p>
<p>There is clearly a lot more fun we can have by playing around with<br>
priors. I anticipate that applied scientists, like my ecological<br>
colleagues, will soon start paying much more attention to prior choice.</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/21/bayesian-smoothing.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/_7zJQVHjKAA" height="1"></a></div>