<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/first-steps-with-mrf-smooths/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="http://www.fromthebottomoftheheap.net/2017/10/19/first-steps-with-mrf-smooths/"> From the Bottom of the Heap - R</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p>
One of the specialist smoother types in the <strong>mgcv</strong> package is the Markov Random Field (MFR) smooth. This smoother essentially allows you to model spatial data with an intrinsic Gaussian Markov random field (GMRF). GRMFs are often used for spatial data measured over discrete spatial regions. MRFs are quite flexible as you can think about them as representing an undirected graph whose nodes are your samples and the connections between the nodes are specified via a neighbourhood structure. I’ve become interested in using these MRF smooths to include information about relationships between species. However, these smooths are not widely documented in the smoothing literature so working out how best to use them to do what we want has been a little tricky once you move beyond the typical spatial examples. As a result I’ve been fiddling with these smooths, fitting them to some spatial data I came across in a tutorial <a rel="nofollow" href="https://pudding.cool/process/regional_smoothing/" target="_blank">Regional Smoothing in R</a> from The Pudding. In this post I take a quick look at how to use the MRF smooth in <strong>mgcv</strong> to model a discrete spatial data set from the US Census Bureau.
</p>
<p>
In that tutorial, the example data are taken from the US Census Bureau via a shapefile prepared by the author. After a little munging — quite a few steps are missing from the tutorial — I managed to get data from the shapefile that matched what was used in the tutorial. The data are on county level percentages of US adults whose highest level of education attainment is a high school diploma. The raw data are shown in the figure below
</p>
<figure>
<img src="https://i2.wp.com/www.fromthebottomoftheheap.net/assets/img/posts/first-steps-with-mrf-smooths-hsd-data-plot-1.png?w=456" data-recalc-dims="1"><br>
</figure>
<p>
To follow along, you’ll need to download the example <a rel="nofollow" href="https://github.com/polygraph-cool/smoothing_tutorial/blob/master/us_county_hs_only.zip" target="_blank">shapefile</a> provided by the author of the post on The Pudding. The shapefile(s) are in a ZIP, which I extracted into the working directory; the code below assumes this.
</p>
<p>
This post will make use of the following set of package; load them now, as shown below, and install any that you may be missing
</p>
<figure>
<pre><code data-lang="r"><span>library</span><span>(</span><span>'rgdal'</span><span>)</span><span>
</span><span>library</span><span>(</span><span>'proj4'</span><span>)</span><span>
</span><span>library</span><span>(</span><span>'spdep'</span><span>)</span><span>
</span><span>library</span><span>(</span><span>'mgcv'</span><span>)</span><span>
</span><span>library</span><span>(</span><span>'ggplot2'</span><span>)</span><span>
</span><span>library</span><span>(</span><span>'dplyr'</span><span>)</span><span>
</span><span>library</span><span>(</span><span>'viridis'</span><span>)</span></code></pre>
</figure>
<p>
Assuming you have extracted the shapefile, we load it into R using <code>readOGR()</code>
</p>
<figure>
<pre><code data-lang="r"><span>shp</span><span> </span><span>&lt;-</span><span> </span><span>readOGR</span><span>(</span><span>'.'</span><span>,</span><span> </span><span>'us_county_hs_only'</span><span>)</span></code></pre>
</figure>
<p>
and do some data munging
</p>
<figure>
<pre><code data-lang="r"><span>## select only mainland US counties
</span><span>states</span><span> </span><span>&lt;-</span><span> </span><span>c</span><span>(</span><span>1</span><span>,</span><span>4</span><span>:</span><span>6</span><span>,</span><span> </span><span>8</span><span>:</span><span>13</span><span>,</span><span> </span><span>16</span><span>:</span><span>42</span><span>,</span><span> </span><span>44</span><span>:</span><span>51</span><span>,</span><span> </span><span>53</span><span>:</span><span>56</span><span>)</span><span>
</span><span>shp</span><span> </span><span>&lt;-</span><span> </span><span>shp</span><span>[</span><span>shp</span><span>$</span><span>STATEFP</span><span> </span><span>%in%</span><span> </span><span>sprintf</span><span>(</span><span>'%02i'</span><span>,</span><span> </span><span>states</span><span>),</span><span> </span><span>]</span><span>
</span><span>df</span><span> </span><span>&lt;-</span><span> </span><span>droplevels</span><span>(</span><span>as</span><span>(</span><span>shp</span><span>,</span><span> </span><span>'data.frame'</span><span>))</span><span>
</span><span>## project data
</span><span>aea.proj</span><span> </span><span>&lt;-</span><span> </span><span>&quot;+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=37.5 +lon_0=-96&quot;</span><span>
</span><span>shp</span><span> </span><span>&lt;-</span><span> </span><span>spTransform</span><span>(</span><span>shp</span><span>,</span><span> </span><span>CRS</span><span>(</span><span>aea.proj</span><span>))</span><span> </span><span># project to Albers
</span><span>shpf</span><span> </span><span>&lt;-</span><span> </span><span>fortify</span><span>(</span><span>shp</span><span>,</span><span> </span><span>region</span><span> </span><span>=</span><span> </span><span>'GEOID'</span><span>)</span><span>
</span><span>## Need a proportion for fitting
</span><span>df</span><span> </span><span>&lt;-</span><span> </span><span>transform</span><span>(</span><span>df</span><span>,</span><span> </span><span>hsd</span><span> </span><span>=</span><span> </span><span>hs_pct</span><span> </span><span>/</span><span> </span><span>100</span><span>)</span></code></pre>
</figure>
<p>
The shapefile contains US Census Bureau data for all US counties, including many that are far from the continental USA. The tutorial from The Pudding doesn’t go into how they removed, or how they drew a map without these additional counties. For our purposes they may cause complications when we try to model them using the MRF smooth. I’m sure the modelling approach can handle data like this, but as I wanted to achieve something that followed the tutorial I’ve removed everything not linked to the continental US landmass, including (I’m sorry!), Alaska and Hawaii — my <strong>ggplot</strong> and mapping skills aren’t yet good enough to move Alaska and Hawaii to the bottom left of such maps.
</p>
<p>
The data were projected using the Albers equal area projection and subsequently passed to the <code>fortify()</code> method from <strong>ggplot2</strong> to get a version of the county polygons suitable for plotting with that package.
</p>
<p>
Finally, I created a new variable <code>hsd</code> which is just the variable <code>hs_pct</code> divided by 100. This creates a proportion that we’ll need for model fitting as you’ll see shortly.
</p>
<p>
Before we can model these data with <code>gam()</code>, we need to create the supporting information that <code>gam()</code> will use to create the MRF smooth penalty. The penalty matrix in an MRF smooth is based on the neighbourhood structure of the observations. There are three ways to pass this information to <code>gam()</code>
</p>
<ol type="1">
<li>
as a list of polygons (not <code>SpatialPolygons</code>, I believe)
</li>
<li>
as a list containing the neighbourhood structure, or
</li>
<li>
the raw penalty matrix itself.
</li>
</ol>
<p>
Options 1 and 3 aren’t easily doable as far I can see — <code>gam()</code> isn’t expecting the sort of object we created when we imported the shapefile and nobody want’s to build a penalty matrix by hand! Thankfully option 2, the neighbourhood structure is relatively easy to create. For that I use the <code>poly2nb()</code> function from the <strong>spdep</strong> package. This function takes a shapefile and works out which regions are neighbours of any other region by virtue of them sharing a border. To make sure everything matches up nicely in the way <code>gam()</code> wants this list, we specify that the region IDs should be the <code>GEOID</code>s from the original data set (the <code>GEOID</code> uniquely identifies each county) and we have to set the <code>names</code> attribute on the neighbouthood list to match these unique IDs
</p>
<figure>
<pre><code data-lang="r"><span>nb</span><span> </span><span>&lt;-</span><span> </span><span>poly2nb</span><span>(</span><span>shp</span><span>,</span><span> </span><span>row.names</span><span> </span><span>=</span><span> </span><span>df</span><span>$</span><span>GEOID</span><span>)</span><span>
</span><span>names</span><span>(</span><span>nb</span><span>)</span><span> </span><span>&lt;-</span><span> </span><span>attr</span><span>(</span><span>nb</span><span>,</span><span> </span><span>&quot;region.id&quot;</span><span>)</span></code></pre>
</figure>
<p>
The result of the previous chunk is a list whose names map on to the levels of the <code>GEOID</code> factor. The values in each element of <code>nb</code> index the elements of <code>nb</code> that are neighbours of the current element
</p>
<figure>
<pre><code data-lang="r"><span>str</span><span>(</span><span>nb</span><span>[</span><span>1</span><span>:</span><span>6</span><span>])</span></code></pre>
</figure>
<figure>
<pre><code data-lang="text">List of 6
 $ 19107: int [1:6] 1417 1464 1632 2277 2278 2851
 $ 19189: int [1:6] 551 1414 2151 2452 2846 2849
 $ 20093: int [1:7] 5 557 1064 1142 1437 1441 2978
 $ 20123: int [1:5] 1469 1565 2648 2966 2977
 $ 20187: int [1:7] 3 554 1142 1441 1620 2142 2238
 $ 21005: int [1:7] 582 583 953 954 1770 1861 2169</code></pre>
</figure>
<p>
With that done we can now fit the GAM. Fitting this is going to take a wee while (over 3 hours for the full rank MRF, using 6 threads, on a reasonably powerful 3-year old workstation with dual 4-core Xeon processors). To specify an MRF smooth we use the <code>bs</code> argument to the <code>s()</code> function, setting it to <code>bs = ‘mrf’</code>. The neighbourhood list is passed via the <code>xt</code> argument, which takes a list as a value; here we specify a component <code>nb</code> which takes our neighbourhood list <code>nb</code>. The final set-up variable to consider is whether to fit a full rank MRF, where a coefficient for each county will be estimated, or a reduced rank MRF, wherein the MRF is represented using fewer coefficients and counties are mapped to the smaller set of coefficients. The rank of the MRF smooth is set using the <code>k</code> argument. The default is to fit a full rank MRF, whilst setting <code>k &lt; NROW(data)</code> will result ins a reduced-rank MRF being etimated.
</p>
<p>
The full rank MRF model is estimated using
</p>
<figure>
<pre><code data-lang="r"><span>ctrl</span><span> </span><span>&lt;-</span><span> </span><span>gam.control</span><span>(</span><span>nthreads</span><span> </span><span>=</span><span> </span><span>6</span><span>)</span><span> </span><span># use 6 parallel threads, reduce if fewer physical CPU cores
</span><span>m</span><span>1</span><span> </span><span>&lt;-</span><span> </span><span>gam</span><span>(</span><span>hsd</span><span> </span><span>~</span><span> </span><span>s</span><span>(</span><span>GEOID</span><span>,</span><span> </span><span>bs</span><span> </span><span>=</span><span> </span><span>'mrf'</span><span>,</span><span> </span><span>xt</span><span> </span><span>=</span><span> </span><span>list</span><span>(</span><span>nb</span><span> </span><span>=</span><span> </span><span>nb</span><span>)),</span><span> </span><span># define MRF smooth
</span><span> </span><span>data</span><span> </span><span>=</span><span> </span><span>df</span><span>,</span><span>
</span><span>method</span><span> </span><span>=</span><span> </span><span>'REML'</span><span>,</span><span> </span><span># fast version of REML smoothness selection
</span><span> </span><span>control</span><span> </span><span>=</span><span> </span><span>ctrl</span><span>,</span><span>
</span><span>family</span><span> </span><span>=</span><span> </span><span>betar</span><span>())</span><span> </span><span>#</span><span> </span><span>fit</span><span> </span><span>a</span><span> </span><span>beta</span><span> </span><span>regression</span></code></pre>
</figure>
<p>
As the response is a proportion, the fitted GAM uses the beta distribution as the conditional distribution of the response. The default link in the logit, just as it is in for the binomial distribution, and insures that fitted values on the scale of the linear predictor are mapped onto the allowed range for proportions of 0–1.
</p>
<p>
The final model uses in the region of 1700 effective degrees of freedom. This is the smoothness penalty at work; rather than 3108 individual coefficients, the smoothness invoked to try to arrange for neighbouring counties to have similar coefficients has shrunk away almost half of the complexity implied by the full rank MRF.
</p>
<figure>
<pre><code data-lang="r"><span>summary</span><span>(</span><span>m</span><span>1</span><span>)</span></code></pre>
</figure>
<figure>
<pre><code data-lang="text">Family: Beta regression(179.532) 
Link function: logit 

Formula:
hsd ~ s(GEOID, bs = &quot;mrf&quot;, xt = list(nb = nb))

Parametric coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -0.63806    0.00283  -225.5   &lt;2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Approximate significance of smooth terms:
          edf Ref.df Chi.sq p-value    
s(GEOID) 1732   3107   9382  &lt;2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

R-sq.(adj) =  0.769   Deviance explained = 89.7%
-REML =  -4544  Scale est. = 1         n = 3108</code></pre>
</figure>
<p>
Whilst the penalty enforces smoothness, further smoothness can be enforced by fitting a reduced rank MARF. In the next code block I fit models with <code>k = 300</code> and <code>k = 30</code> respectively, which imply considerable smoothing relative to the full rank model.
</p>
<figure>
<pre><code data-lang="r"><span>## rank 300 MRF
</span><span>m</span><span>2</span><span> </span><span>&lt;-</span><span> </span><span>gam</span><span>(</span><span>hsd</span><span> </span><span>~</span><span> </span><span>s</span><span>(</span><span>GEOID</span><span>,</span><span> </span><span>bs</span><span> </span><span>=</span><span> </span><span>'mrf'</span><span>,</span><span> </span><span>k</span><span> </span><span>=</span><span> </span><span>300</span><span>,</span><span> </span><span>xt</span><span> </span><span>=</span><span> </span><span>list</span><span>(</span><span>nb</span><span> </span><span>=</span><span> </span><span>nb</span><span>)),</span><span>
</span><span>data</span><span> </span><span>=</span><span> </span><span>df</span><span>,</span><span> </span><span>method</span><span> </span><span>=</span><span> </span><span>'REML'</span><span>,</span><span> </span><span>control</span><span> </span><span>=</span><span> </span><span>ctrl</span><span>,</span><span>
</span><span>family</span><span> </span><span>=</span><span> </span><span>betar</span><span>())</span><span>
</span><span>## rank 30 MRF
</span><span>m</span><span>3</span><span> </span><span>&lt;-</span><span> </span><span>gam</span><span>(</span><span>hsd</span><span> </span><span>~</span><span> </span><span>s</span><span>(</span><span>GEOID</span><span>,</span><span> </span><span>bs</span><span> </span><span>=</span><span> </span><span>'mrf'</span><span>,</span><span> </span><span>k</span><span> </span><span>=</span><span> </span><span>30</span><span>,</span><span> </span><span>xt</span><span> </span><span>=</span><span> </span><span>list</span><span>(</span><span>nb</span><span> </span><span>=</span><span> </span><span>nb</span><span>)),</span><span>
</span><span>data</span><span> </span><span>=</span><span> </span><span>df</span><span>,</span><span> </span><span>method</span><span> </span><span>=</span><span> </span><span>'REML'</span><span>,</span><span> </span><span>control</span><span> </span><span>=</span><span> </span><span>ctrl</span><span>,</span><span>
</span><span>family</span><span> </span><span>=</span><span> </span><span>betar</span><span>())</span></code></pre>
</figure>
<p>
To visualise the different fits we need to generate predicted values on the response scale for each county and add this data to the county data <code>df</code>
</p>
<figure>
<pre><code data-lang="r"><span>df</span><span> </span><span>&lt;-</span><span> </span><span>transform</span><span>(</span><span>df</span><span>,</span><span>
</span><span>mrfFull</span><span> </span><span>=</span><span> </span><span>predict</span><span>(</span><span>m</span><span>1</span><span>,</span><span> </span><span>type</span><span> </span><span>=</span><span> </span><span>'response'</span><span>),</span><span>
</span><span>mrfRrank300</span><span> </span><span>=</span><span> </span><span>predict</span><span>(</span><span>m</span><span>2</span><span>,</span><span> </span><span>type</span><span> </span><span>=</span><span> </span><span>'response'</span><span>),</span><span>
</span><span>mrfRrank30</span><span> </span><span>=</span><span> </span><span>predict</span><span>(</span><span>m</span><span>3</span><span>,</span><span> </span><span>type</span><span> </span><span>=</span><span> </span><span>'response'</span><span>))</span></code></pre>
</figure>
<p>
Before we can plot these fitted values we need to merge <code>df</code> with the fortified shapefile
</p>
<figure>
<pre><code data-lang="r"><span>## merge data with fortified shapefile
</span><span>mdata</span><span> </span><span>&lt;-</span><span> </span><span>left_join</span><span>(</span><span>shpf</span><span>,</span><span> </span><span>df</span><span>,</span><span> </span><span>by</span><span> </span><span>=</span><span> </span><span>c</span><span>(</span><span>'id'</span><span> </span><span>=</span><span> </span><span>'GEOID'</span><span>))</span></code></pre>
</figure>
<figure>
<pre><code data-lang="text">Warning: Column `id`/`GEOID` joining character vector and factor, coercing
into character vector</code></pre>
</figure>
<p>
To facilitate plotting with <strong>ggplot2</strong> I begin by creating some fixed plot components, like the theme, scale, and labels
</p>
<figure>
<pre><code data-lang="r"><span>theme_map</span><span> </span><span>&lt;-</span><span> </span><span>function</span><span>(</span><span>...</span><span>)</span><span> </span><span>{</span><span>
</span><span>theme_minimal</span><span>()</span><span> </span><span>+</span><span>
</span><span>theme</span><span>(</span><span>...</span><span>,</span><span>
</span><span>axis.line</span><span> </span><span>=</span><span> </span><span>element_blank</span><span>(),</span><span>
</span><span>axis.text.x</span><span> </span><span>=</span><span> </span><span>element_blank</span><span>(),</span><span>
</span><span>axis.text.y</span><span> </span><span>=</span><span> </span><span>element_blank</span><span>(),</span><span>
</span><span>axis.ticks</span><span> </span><span>=</span><span> </span><span>element_blank</span><span>(),</span><span>
</span><span>axis.title.x</span><span> </span><span>=</span><span> </span><span>element_blank</span><span>(),</span><span>
</span><span>axis.title.y</span><span> </span><span>=</span><span> </span><span>element_blank</span><span>(),</span><span>
</span><span>panel.border</span><span> </span><span>=</span><span> </span><span>element_blank</span><span>())</span><span>
</span><span>}</span><span>
</span><span>myTheme</span><span> </span><span>&lt;-</span><span> </span><span>theme_map</span><span>(</span><span>legend.position</span><span> </span><span>=</span><span> </span><span>'bottom'</span><span>)</span><span>
</span><span>myScale</span><span> </span><span>&lt;-</span><span> </span><span>scale_fill_viridis</span><span>(</span><span>name</span><span> </span><span>=</span><span> </span><span>'%'</span><span>,</span><span> </span><span>option</span><span> </span><span>=</span><span> </span><span>'plasma'</span><span>,</span><span>
</span><span>limits</span><span> </span><span>=</span><span> </span><span>c</span><span>(</span><span>0.1</span><span>,</span><span> </span><span>0.55</span><span>),</span><span>
</span><span>labels</span><span> </span><span>=</span><span> </span><span>function</span><span>(</span><span>x</span><span>)</span><span> </span><span>x</span><span> </span><span>*</span><span> </span><span>100</span><span>,</span><span>
</span><span>guide</span><span> </span><span>=</span><span> </span><span>guide_colorbar</span><span>(</span><span>direction</span><span> </span><span>=</span><span> </span><span>&quot;horizontal&quot;</span><span>,</span><span>
</span><span>barheight</span><span> </span><span>=</span><span> </span><span>unit</span><span>(</span><span>2</span><span>,</span><span> </span><span>units</span><span> </span><span>=</span><span> </span><span>&quot;mm&quot;</span><span>),</span><span>
</span><span>barwidth</span><span> </span><span>=</span><span> </span><span>unit</span><span>(</span><span>75</span><span>,</span><span> </span><span>units</span><span> </span><span>=</span><span> </span><span>&quot;mm&quot;</span><span>),</span><span>
</span><span>title.position</span><span> </span><span>=</span><span> </span><span>'left'</span><span>,</span><span>
</span><span>title.hjust</span><span> </span><span>=</span><span> </span><span>0.5</span><span>,</span><span>
</span><span>label.hjust</span><span> </span><span>=</span><span> </span><span>0.5</span><span>))</span><span>
</span><span>myLabs</span><span> </span><span>&lt;-</span><span> </span><span>labs</span><span>(</span><span>x</span><span> </span><span>=</span><span> </span><span>NULL</span><span>,</span><span> </span><span>y</span><span> </span><span>=</span><span> </span><span>NULL</span><span>,</span><span> </span><span>title</span><span> </span><span>=</span><span> </span><span>'US Adult Education'</span><span>,</span><span>
</span><span>subtitle</span><span> </span><span>=</span><span> </span><span>'% of adults where high school diploma is highest level education'</span><span>,</span><span>
</span><span>caption</span><span> </span><span>=</span><span> </span><span>'Source: US Census Bureau'</span><span>)</span></code></pre>
</figure>
<p>
I took many of these settings from Timo Grossenbacher’s excellent <a rel="nofollow" href="https://timogrossenbacher.ch/2016/12/beautiful-thematic-maps-with-ggplot2-only/" target="_blank">post on mapping regional demographic data in Switzerland</a>.
</p>
<p>
Now we can plot the fitted proportions. Note that whilst we plot proportions, the colour bar labels are in percentages in keeping with the original data (see the definition for <code>my_scale</code> to see how this was achieved).
</p>
<p>
Fitted values from the full rank MRF are shown below
</p>
<figure>
<pre><code data-lang="r"><span>ggplot</span><span>(</span><span>mdata</span><span>,</span><span> </span><span>aes</span><span>(</span><span>x</span><span> </span><span>=</span><span> </span><span>long</span><span>,</span><span> </span><span>y</span><span> </span><span>=</span><span> </span><span>lat</span><span>,</span><span> </span><span>group</span><span> </span><span>=</span><span> </span><span>group</span><span>))</span><span> </span><span>+</span><span>
</span><span>geom_polygon</span><span>(</span><span>aes</span><span>(</span><span>fill</span><span> </span><span>=</span><span> </span><span>mrfFull</span><span>))</span><span> </span><span>+</span><span>
</span><span>geom_path</span><span>(</span><span>col</span><span> </span><span>=</span><span> </span><span>'black'</span><span>,</span><span> </span><span>alpha</span><span> </span><span>=</span><span> </span><span>0.5</span><span>,</span><span> </span><span>size</span><span> </span><span>=</span><span> </span><span>0.1</span><span>)</span><span> </span><span>+</span><span>
</span><span>coord_equal</span><span>()</span><span> </span><span>+</span><span>
</span><span>myTheme</span><span> </span><span>+</span><span> </span><span>myScale</span><span> </span><span>+</span><span> </span><span>myLabs</span></code></pre>
</figure>
<figure>
<img src="https://i2.wp.com/www.fromthebottomoftheheap.net/assets/img/posts/first-steps-with-mrf-smooths-plot-full-rank-mrf-1.png?w=456" data-recalc-dims="1"><br>
</figure>
<p>
This model expains about 90% of the deviance in the original data. Whilst some smoothing is evident, the fitted values show a considerable about of non-spatial variation. This is most likely due to not including important covariates, such as country average income, which might explain some of the finer scale structure; neighbouring counties with quite different proportions. A more considered analysis would include these and other relevant predictors alongside the MRF.
</p>
<p>
Smoother surfaces can be achieved via the reduced rank MRFs. First the rank 300 MRF
</p>
<figure>
<pre><code data-lang="r"><span>ggplot</span><span>(</span><span>mdata</span><span>,</span><span> </span><span>aes</span><span>(</span><span>x</span><span> </span><span>=</span><span> </span><span>long</span><span>,</span><span> </span><span>y</span><span> </span><span>=</span><span> </span><span>lat</span><span>,</span><span> </span><span>group</span><span> </span><span>=</span><span> </span><span>group</span><span>))</span><span> </span><span>+</span><span>
</span><span>geom_polygon</span><span>(</span><span>aes</span><span>(</span><span>fill</span><span> </span><span>=</span><span> </span><span>mrfRrank300</span><span>))</span><span> </span><span>+</span><span>
</span><span>geom_path</span><span>(</span><span>col</span><span> </span><span>=</span><span> </span><span>'black'</span><span>,</span><span> </span><span>alpha</span><span> </span><span>=</span><span> </span><span>0.5</span><span>,</span><span> </span><span>size</span><span> </span><span>=</span><span> </span><span>0.1</span><span>)</span><span> </span><span>+</span><span>
</span><span>coord_equal</span><span>()</span><span> </span><span>+</span><span>
</span><span>myTheme</span><span> </span><span>+</span><span> </span><span>myScale</span><span> </span><span>+</span><span> </span><span>myLabs</span></code></pre>
</figure>
<figure>
<img src="https://i0.wp.com/www.fromthebottomoftheheap.net/assets/img/posts/first-steps-with-mrf-smooths-plot-rank-300-mrf-1.png?w=456" data-recalc-dims="1"><br>
</figure>
<p>
and next the rank 30 MRF
</p>
<figure>
<pre><code data-lang="r"><span>ggplot</span><span>(</span><span>mdata</span><span>,</span><span> </span><span>aes</span><span>(</span><span>x</span><span> </span><span>=</span><span> </span><span>long</span><span>,</span><span> </span><span>y</span><span> </span><span>=</span><span> </span><span>lat</span><span>,</span><span> </span><span>group</span><span> </span><span>=</span><span> </span><span>group</span><span>))</span><span> </span><span>+</span><span>
</span><span>geom_polygon</span><span>(</span><span>aes</span><span>(</span><span>fill</span><span> </span><span>=</span><span> </span><span>mrfRrank30</span><span>))</span><span> </span><span>+</span><span>
</span><span>geom_path</span><span>(</span><span>col</span><span> </span><span>=</span><span> </span><span>'black'</span><span>,</span><span> </span><span>alpha</span><span> </span><span>=</span><span> </span><span>0.5</span><span>,</span><span> </span><span>size</span><span> </span><span>=</span><span> </span><span>0.1</span><span>)</span><span> </span><span>+</span><span>
</span><span>coord_equal</span><span>()</span><span> </span><span>+</span><span>
</span><span>myTheme</span><span> </span><span>+</span><span> </span><span>myScale</span><span> </span><span>+</span><span> </span><span>myLabs</span></code></pre>
</figure>
<figure>
<img src="https://i1.wp.com/www.fromthebottomoftheheap.net/assets/img/posts/first-steps-with-mrf-smooths-plot-rank-30-mrf-1.png?w=456" data-recalc-dims="1"><br>
</figure>
<p>
As can be clearly seen from the plots, the degree of smoothness can be controlled effectively via the <code>k</code> argument.
</p>
<p>
In a future post I’ll take a closer look at using MRFs alongside other covariates as part of a model complex spatial modeling exercise.</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.fromthebottomoftheheap.net/2017/10/19/first-steps-with-mrf-smooths/"> From the Bottom of the Heap - 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/YVEIGNNRCWU" height="1"></a></div>