<p><div><div><div><div data-width="100" data-show-faces="false" data-send="true" data-layout="button_count" data-height="21" data-href="https://www.r-bloggers.com/gaussian-predictive-process-models-in-stan/"></div></div><div><a data-url="https://www.r-bloggers.com/gaussian-predictive-process-models-in-stan/" data-counturl="https://www.r-bloggers.com/gaussian-predictive-process-models-in-stan/" data-count="horizontal" data-text="Gaussian predictive process models in Stan" class="twitter-share-button" data-via="rbloggers" href="https://twitter.com/share"></a></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="http://mbjoseph.github.io/2016/08/14/gpp.html"> Maxwell B. Joseph</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p>Gaussian process (GP) models are computationally demanding for large datasets.<br>
Much work has been done to avoid expensive matrix operations that arise in parameter estimation with larger datasets via sparse and/or reduced rank covariance matrices (<a target="_blank" href="http://arxiv.org/pdf/1406.7343.pdf">Datta et al. 2016</a> provide a nice review).<br>
What follows is an implementation of a spatial Gaussian predictive process Poisson GLM in <a target="_blank" href="http://mc-stan.org/">Stan</a>, following <a target="_blank" href="http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2743161/">Finley et al. 2009</a>: <em>Improving the performance of predictive process modeling for large datasets</em>, with a comparison to a full rank GP model in terms of execution time and MCMC efficiency.</p>
<h2>Generative model</h2>
<p>Assume we have $n$ spatially referenced counts $y(s)$ made at spatial locations $s_1, s_2, …, s_n$, which depend on a latent mean zero Gaussian process with an isotropic stationary exponential covariance function:</p>
<p></p>
<p></p>
<p></p>
<p>where $y(s)$ is the response at location $s$, $X$ is an $n \times p$ design matrix, $\beta$ is a length $p$ parameter vector, $\sigma^2$ is a “nugget” parameter, $\eta^2$ is the variance parameter of the Gaussian process, $d_{ij}$ is a spatial distance between locations $s_i$ and $s_j$, and $\phi$ determines how quickly the correlation in $w$ decays as distance increases.<br>
For simplicity, the point locations are uniformly distributed in a 2d unit square spatial region.</p>
<p>To estimate this model in a Bayesian context, we might be faced with taking a Cholesky decomposition of the $n \times n$ matrix $C(d)$ at every iteration in an MCMC algorithm, a costly operation for large $n$.</p>
<h2>Gaussian predictive process representation</h2>
<p>Computational benefits of Gaussian predictive process models arise from the estimation of the latent Gaussian process at $m « n$ locations (knots).<br>
Instead of taking the Cholesky factorization of the $n \times n$ covariance matrix, we instead factorize the $m \times m$ covariance matrix corresponding to the covariance in the latent spatial process among knots.<br>
Knot placement is a non-trivial topic, but for the purpose of illustration let’s place knots on a grid over the spatial region of interest.<br>
Note that here the knots are fixed, but it is possible to model knot locations stochastically as in <a target="_blank" href="http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3268014/">Guhaniyogi et al. 2012</a>.</p>
<p><img alt="" src="https://i2.wp.com/mbjoseph.github.io/images/knots.png?w=456" data-recalc-dims="1"></p>
<p>Knots are shown as stars and the points are observations.</p>
<p>We will replace $w$ above with an estimate of $w$ that is derived from a reduced rank representation of the latent spatial process.<br>
Below, the vector $\tilde{\boldsymbol{\epsilon}}$ corrects for bias (underestimation of $\eta$ and overestimation of $\sigma$) as an extension of Banerjee et al. 2008, and $\mathcal{C}^T(\theta) \mathcal{C}^{*-1}(\theta)$ relates $w$ at desired point locations to the value of the latent GP at the knot locations, where $\mathcal{C}^T(\theta)$ is an $n \times m$ matrix that gets multiplied by the inverse of the $m \times m$ covariance matrix for the latent spatial process at the knots.<br>
For a complete and general description see Finley et al. 2009, but here is the jist of the univariate model:</p>
<p></p>
<p></p>
<p></p>
<p>with priors for $\eta$, $\sigma$, and $\phi$ completing a Bayesian specification.</p>
<p>This approach scales well for larger datasets relative to a full rank GP model.<br>
Comparing the number of effective samples per unit time for the two approaches across a range of sample sizes, the GPP executes more quickly and is more efficient for larger datasets.<br>
Code for this simulation is available on GitHub: <a target="_blank" href="https://github.com/mbjoseph/gpp-speed-test">https://github.com/mbjoseph/gpp-speed-test</a>.</p>
<p><img alt="" src="https://i0.wp.com/mbjoseph.github.io/images/gpp_compare.png?w=456" data-recalc-dims="1"></p>
<h2>Relevant papers</h2>
<p><a target="_blank" href="http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2743161/">Finley, Andrew O., et al. “Improving the performance of predictive process modeling for large datasets.” Computational statistics &amp; data analysis 53.8 (2009): 2873-2884.</a></p>
<p><a target="_blank" href="http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2741335/">Banerjee, Sudipto, et al. “Gaussian predictive process models for large spatial data sets.” Journal of the Royal Statistical Society: Series B (Statistical Methodology) 70.4 (2008): 825-848.</a></p>
<p><a target="_blank" href="http://arxiv.org/pdf/1406.7343.pdf">Datta, Abhirup, et al. “Hierarchical nearest-neighbor Gaussian process models for large geostatistical datasets.” Journal of the American Statistical Association. Accepted (2015).</a></p>
<p><a target="_blank" href="http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3268014/">Guhaniyogi, Rajarshi, et al. “Adaptive Gaussian predictive process models for large spatial datasets.” Environmetrics 22.8 (2011): 997-1007.</a></p>
<p><div><div><div><div data-width="55" data-show-faces="false" data-send="true" data-layout="box_count" data-height="62" data-href="https://www.r-bloggers.com/gaussian-predictive-process-models-in-stan/"></div></div><div><a data-url="https://www.r-bloggers.com/gaussian-predictive-process-models-in-stan/" data-counturl="https://www.r-bloggers.com/gaussian-predictive-process-models-in-stan/" data-count="vertical" data-text="Gaussian predictive process models in Stan" class="twitter-share-button" data-via="rbloggers" href="https://twitter.com/share"></a></div><div></div></div><div></div></div>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="http://mbjoseph.github.io/2016/08/14/gpp.html"> Maxwell B. Joseph</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 title="The R Project for Statistical Computing" rel="nofollow" href="https://www.r-project.org/">R</a> news and <a title="R tutorials" rel="nofollow" href="https://www.r-bloggers.com/?s=tutorial">tutorials</a> on topics such as: <a title="Data science" rel="nofollow" href="https://www.r-bloggers.com/?s=data%20science">Data science</a>, <a title="Big Data" rel="nofollow" href="https://www.r-bloggers.com/?s=Big%20Data">Big Data, <a title="R jobs" rel="nofollow" href="https://www.r-users.com/">R jobs</a>, visualization (<a title="ggplot and ggplot2 tutorials" rel="nofollow" href="https://www.r-bloggers.com/?s=ggplot2">ggplot2</a>, <a title="Boxplots using lattice and ggplot2 tutorials" rel="nofollow" href="https://www.r-bloggers.com/?s=boxplot">Boxplots</a>, <a title="Maps and gis" rel="nofollow" href="https://www.r-bloggers.com/?s=map">maps</a>, <a title="Animation in R" rel="nofollow" href="https://www.r-bloggers.com/?s=animation">animation</a>), programming (<a title="RStudio IDE for R" rel="nofollow" href="https://www.r-bloggers.com/?s=RStudio">RStudio</a>, <a title="Sweave and literate programming" rel="nofollow" href="https://www.r-bloggers.com/?s=sweave">Sweave</a>, <a title="LaTeX in R" rel="nofollow" href="https://www.r-bloggers.com/?s=LaTeX">LaTeX</a>, <a title="SQL and databases" rel="nofollow" href="https://www.r-bloggers.com/?s=SQL">SQL</a>, <a title="Eclipse IDE for R" rel="nofollow" href="https://www.r-bloggers.com/?s=eclipse">Eclipse</a>, <a title="git and github, Version Control System" rel="nofollow" href="https://www.r-bloggers.com/?s=git">git</a>, <a title="Large data in R using Hadoop" rel="nofollow" href="https://www.r-bloggers.com/?s=hadoop">hadoop</a>, <a title="Web Scraping of google, facebook, yahoo, twitter and more using R" rel="nofollow" href="https://www.r-bloggers.com/?s=Web+Scraping">Web Scraping</a>) statistics (<a title="Regressions and ANOVA analysis tutorials" rel="nofollow" href="https://www.r-bloggers.com/?s=regression">regression</a>, <a title="principal component analysis tutorial" rel="nofollow" href="https://www.r-bloggers.com/?s=PCA">PCA</a>, <a title="Time series" rel="nofollow" href="https://www.r-bloggers.com/?s=time+series">time series</a>, <a title="finance trading" rel="nofollow" href="https://www.r-bloggers.com/?s=trading">trading</a>) and more...
<img height="1" alt="" width="1" src="http://feeds.feedburner.com/~r/RBloggers/~4/VhrsuXoJzYM"></a></div>