<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/an-introduction-to-greta/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="https://rviews.rstudio.com/2018/04/23/on-first-meeting-greta/"> R Views</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p>I was surprised by <a rel="nofollow" href="https://greta-dev.github.io/greta/" target="_blank"><code>greta</code></a>. I had assumed that the <a rel="nofollow" href="https://tensorflow.rstudio.com/" target="_blank"><code>tensorflow</code></a> and <a rel="nofollow" href="https://rstudio.github.io/reticulate/" target="_blank"><code>reticulate</code></a> packages would eventually enable R developers to look beyond deep learning applications and exploit the <a rel="nofollow" href="https://www.tensorflow.org/" target="_blank">TensorFlow</a> platform to create all manner of production-grade statistical applications. But I wasn’t thinking Bayesian. After all, <a rel="nofollow" href="http://mc-stan.org/" target="_blank">Stan</a> is probably everything a Bayesian modeler could want. Stan is a powerful, production-level probability distribution modeling engine with a slick <a rel="nofollow" href="https://cran.r-project.org/web/packages/rstan/rstan.pdf" target="_blank">R interface</a>, deep documentation, and a dedicated development team.</p>
<p>But <code>greta</code> lets users write TensorFlow-based Bayesian models directly in R! What could be more charming? <code>greta</code> removes the barrier of learning an intermediate modeling language while still promising to deliver high-performance MCMC models that run anywhere TensorFlow can go.</p>
<p>In this post, I’ll introduce you to <code>greta</code> with a simple model used by Richard McElreath in section 8.3 of his iconoclastic book: <a rel="nofollow" href="http://xcelab.net/rm/statistical-rethinking/" target="_blank">Statistical Rethinking: A Bayesian Course with Examples in R and Stan</a>. This model seeks to explain the log of a country’s GDP based on a measure of terrain ruggedness while controlling for whether or not the country is in Africa. I am going to use it just to illustrate MCMC sampling with <code>greta</code>. The extended example in McElreath’s book, however, is a meditation on the subtleties of modeling interactions, and is well worth studying.</p>
<p>First, we load the required packages and fetch the data. <code>DiagrammeR</code> is for plotting the TensorFlow flow diagram of the model, and <code>bayesplot</code> is used to plot trace diagrams of the Markov chains. The rugged data set which provides 52 variables for 234 is fairly interesting, but we will use a trimmed-down data set with only 170 counties and three variables.</p>
<pre><code>library(rethinking)
library(greta)
library(DiagrammeR)
library(bayesplot)
library(ggplot2)

# Example from section 8.3 Statistical Rethinking
data(rugged)
d &lt;- rugged
d$log_gdp &lt;- log(d$rgdppc_2000)
dd &lt;- d[complete.cases(d$rgdppc_2000), ]
dd_trim &lt;- dd[ , c(&quot;log_gdp&quot;,&quot;rugged&quot;,&quot;cont_africa&quot;)]
head(dd_trim)</code></pre>
<pre><code>##     log_gdp rugged cont_africa
## 3  7.492609  0.858           1
## 5  8.216929  3.427           0
## 8  9.933263  0.769           0
## 9  9.407032  0.775           0
## 10 7.792343  2.688           0
## 12 9.212541  0.006           0</code></pre>
<pre><code>set.seed(1234)</code></pre>
<p>In this section of code, we set up the TensorFlow data structures. The first step is to move the data into <code>greta</code> arrays. These data structures behave similarly to R arrays in that they can be manipulated with functions. However, <code>greta</code> doesn’t immediately calculate values for new arrays. It works out the size and shape of the result and creates a place-holder data structure.</p>
<pre><code>#data
g_log_gdp &lt;- as_data(dd_trim$log_gdp)
g_rugged &lt;- as_data(dd_trim$rugged)
g_cont_africa &lt;- as_data(dd_trim$cont_africa)</code></pre>
<p>In this section, we set up the Bayesian model. All parameters need prior probability distributions. Note that the parameters <code>a</code>, <code>bR</code>, <code>bA</code>, <code>bAR</code>, <code>sigma</code>, and <code>mu</code> are all new <code>greta</code> arrays that don’t contain any data. <code>a</code> is 1 x 1 array and <code>mu</code> is a 170 x 1 array with one slot for each observation.</p>
<p>The <code>distribution()</code> function sets up the likelihood function for the model.</p>
<pre><code># Variables and Priors

a &lt;- normal(0, 100)
bR &lt;- normal(0, 10)
bA &lt;- normal(0, 10)
bAR &lt;- normal(0,10)
sigma &lt;- cauchy(0,2,truncation=c(0,Inf))

a  # Look at this greata array</code></pre>
<pre><code>## greta array (variable following a normal distribution)
## 
##      [,1]
## [1,]  ?</code></pre>
<pre><code># operations
mu &lt;- a + bR*g_rugged + bA*g_cont_africa + bAR*g_rugged*g_cont_africa

dim(mu)</code></pre>
<pre><code>## [1] 170   1</code></pre>
<pre><code># likelihood
distribution(g_log_gdp) = normal(mu, sigma)</code></pre>
<p>The <code>model()</code> function does all of the work. It fits the model and produces a fairly complicated object organized as three lists that contain, respectively, the R6 class, TensorFlow structures, and the various <code>greta</code> data arrays.</p>
<pre><code># defining the model
mod &lt;- model(a,bR,bA,bAR,sigma)

str(mod,give.attr=FALSE,max.level=1)</code></pre>
<pre><code>## List of 3
##  $ dag                 :Classes 'dag_class', 'R6' </code></pre>