diff --git a/docs/index.html b/docs/index.html
index 75e6ca26..d1a49c33 100755
--- a/docs/index.html
+++ b/docs/index.html
@@ -7,7 +7,7 @@
 <meta name="viewport" content="width=device-width, initial-scale=1.0, user-scalable=yes">
 
 <meta name="author" content="Manny Gimond">
-<meta name="dcterms.date" content="2024-02-03">
+<meta name="dcterms.date" content="2024-02-05">
 
 <title>Exploratory Data Analysis in R</title>
 <style>
@@ -334,7 +334,7 @@
 </li>
           <li class="sidebar-item">
   <div class="sidebar-item-container"> 
-  <a href="./robustness.html" class="sidebar-item-text sidebar-link">
+  <a href="./case_study.html" class="sidebar-item-text sidebar-link">
  <span class="menu-text"><span class="chapter-number">23</span>&nbsp; <span class="chapter-title">A working example: t-tests and re-expression</span></span></a>
   </div>
 </li>
@@ -457,7 +457,7 @@ <h1 class="title">Exploratory Data Analysis in R</h1>
     <div>
     <div class="quarto-title-meta-heading">Published</div>
     <div class="quarto-title-meta-contents">
-      <p class="date">February 3, 2024</p>
+      <p class="date">February 5, 2024</p>
     </div>
   </div>
   
@@ -471,13 +471,13 @@ <h1 class="title">Exploratory Data Analysis in R</h1>
 <h1 class="unnumbered">Preface</h1>
 <p>This book is a compilation of lecture notes used in an <em>Exploratory Data Analysis in R</em> course taught to undergraduates at Colby College. The course assumes little to no background in quantitative analysis nor in computer programming and was first taught in Spring, 2015. The course introduces students to data manipulation in R, data exploration (in the spirit of John Tukey’s EDA) and the R markdown language. Many of the visualization techniques are adopted from William Cleveland’s <em>Data Visualization</em> book.</p>
 <p>The base R plotting environment and the <code>ggplot2</code> ecosystem are used throughout this book. While a chapter is dedicated to the <code>lattice</code> plotting package, its functions are not used outside of that chapter given that <code>ggplot2</code> offers many of <code>lattice</code>’s functionality.</p>
-<p>While great effort is made to adopt a consist plotting environment throughout this book (this being <code>ggplot2</code>, for the most part), a few topics (including the q-q plot and the median polish) will benefit from custom plotting functions available in the <code>tukeyedar</code> package. The <a href="https://github.com/mgimond/tukeyedar/">package</a> can be downloaded from GitHub via the command:</p>
-<div class="cell" data-hash="index_cache/html/unnamed-chunk-1_92abb4a22322b722d1af3b88b82b6e95">
+<p>While great effort is made to adopt a consistent plotting environment throughout this book (this being <code>ggplot2</code>, for the most part), a few topics (including the q-q plot and the median polish) will benefit from custom plotting functions available in the <code>tukeyedar</code> package. The <a href="https://github.com/mgimond/tukeyedar/">package</a> can be downloaded from GitHub via the command:</p>
+<div class="cell">
 <div class="sourceCode cell-code" id="cb1"><pre class="sourceCode r code-with-copy"><code class="sourceCode r"><span id="cb1-1"><a href="#cb1-1" aria-hidden="true" tabindex="-1"></a>devtools<span class="sc">::</span><span class="fu">install_github</span>(<span class="st">"mgimond/tukeyedar"</span>)</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
 </div>
 <p>Note that installing the GitHub package will require that the <code>devtools</code> package be installed first.</p>
 <p>Functions making use of the <code>tukeyedar</code> package will be highlighted in a peach/pink code block as opposed to the default light yellow code block used for all other code blocks. For example, if <code>tukeyedar</code>’s <code>eda_qq</code> function is used, the code block will take on the following appearance:</p>
-<div class="cell" data-hash="index_cache/html/unnamed-chunk-2_7743d95d3556383c2327130315418249">
+<div class="cell">
 <div class="sourceCode cell-code" id="cb2"><pre class="sourceCode r eda code-with-copy"><code class="sourceCode r"><span id="cb2-1"><a href="#cb2-1" aria-hidden="true" tabindex="-1"></a><span class="fu">library</span>(tukeyedar)</span>
 <span id="cb2-2"><a href="#cb2-2" aria-hidden="true" tabindex="-1"></a><span class="fu">eda_qq</span>(Tenor, Bass)</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
 </div>
diff --git a/docs/intro.html b/docs/intro.html
index 429cb45e..edda3dd1 100755
--- a/docs/intro.html
+++ b/docs/intro.html
@@ -333,7 +333,7 @@
 </li>
           <li class="sidebar-item">
   <div class="sidebar-item-container"> 
-  <a href="./robustness.html" class="sidebar-item-text sidebar-link">
+  <a href="./case_study.html" class="sidebar-item-text sidebar-link">
  <span class="menu-text"><span class="chapter-number">23</span>&nbsp; <span class="chapter-title">A working example: t-tests and re-expression</span></span></a>
   </div>
 </li>
@@ -650,7 +650,7 @@ <h3 data-number="1.3.4" class="anchored" data-anchor-id="reproducible-analysis">
 <section id="creating-dynamic-documents-using-r-markdown" class="level2" data-number="1.4">
 <h2 data-number="1.4" class="anchored" data-anchor-id="creating-dynamic-documents-using-r-markdown"><span class="header-section-number">1.4</span> Creating dynamic documents using R Markdown</h2>
 <p>Another source of error in the write-up of a report or publication is the linking of tables, figures and statistical summaries to the write-up. Typically, one saves statistical plots as image files then loads the image into the document. However, the figures may have gone through many different iterations resulting in many different versions of the image file in a working folder. Add to this many other figures, data table files and statistical results from various pieces of software, one quickly realizes the potential for embedding the wrong image files in the document or embedding the wrong statistical summaries in the text. Furthermore, the researcher is then required to properly archive and document the provenance of each figure, data table or statistical summary resulting in a complex structure of files and directories in the project folder thus increasing the odds of an irreproducible analysis.</p>
-<p>Confining all of the analysis to a scripting environment such as R can help, but this still does not alleviate the possibility of loading the wrong figure into the document, or forgetting to update a statistical summary in the text when the original data file was revised. A solution to this potential pitfall is to embed the actual analysis and graphic generation process into the document–such environments are called dynamic documents. In this course, we will use the <a href="http://rmarkdown.rstudio.com/">R Markdown authoring tool</a> which embeds R code into the document. An example of an R Markdown document is this course website which was entirely generated in RMarkdown!. You can view the R Markdown files <a href="https://github.com/mgimond/ES218">on this author’s GitHub repository</a>.</p>
+<p>Confining all of the analysis to a scripting environment such as R can help, but this still does not alleviate the possibility of loading the wrong figure into the document, or forgetting to update a statistical summary in the text when the original data file was revised. A solution to this potential pitfall is to embed the actual analysis and graphic generation process into the document–such environments are called dynamic documents. In this course, we will use the <a href="http://rmarkdown.rstudio.com/">R Markdown authoring tool</a> which embeds R code into the document. An example of an R Markdown document is this course website which was entirely generated in RMarkdown! You can view the R Markdown files <a href="https://github.com/mgimond/ES218">on this author’s GitHub repository</a>.</p>
 
 
 </section>
diff --git a/docs/resistant_lines.html b/docs/resistant_lines.html
index 8b013c4f..26eb1ba9 100644
--- a/docs/resistant_lines.html
+++ b/docs/resistant_lines.html
@@ -563,7 +563,7 @@ <h2 data-number="25.1" class="anchored" data-anchor-id="introduction"><span clas
 </div>
 </div>
 <p>The slope dropped from 2 to 0.86 because of a single point!</p>
-<p>If our goal is to explore what the <em>bulk</em> of the data has to say about the phenomena being investigated, we certainly don’t want a few “maverick” values to hijack the analysis. We therefore need a set of fitting tools that minimize the influence of outliers. There are many options out there; most notable are the <strong>resistant line</strong> (that builds off of Tukey’s <strong>3-point summary line</strong>) and the <strong>bisquare</strong> robust estimation method.</p>
+<p>If our goal is to explore what the <em>bulk</em> of the data has to say about the phenomena being investigated, we certainly don’t want a few “maverick” values to hijack the analysis. We therefore need a set of fitting tools that minimize the influence of outliers. There are many options out there; most notable are <strong>Tukey’s resistant line</strong> and the <strong>bisquare</strong> robust estimation method.</p>
 <p>&nbsp;&nbsp;</p>
 </section>
 <section id="robust-lines" class="level2" data-number="25.2">
@@ -571,20 +571,19 @@ <h2 data-number="25.2" class="anchored" data-anchor-id="robust-lines"><span clas
 <p>&nbsp;&nbsp;</p>
 <section id="tukeys-resistant-line" class="level3" data-number="25.2.1">
 <h3 data-number="25.2.1" class="anchored" data-anchor-id="tukeys-resistant-line"><span class="header-section-number">25.2.1</span> Tukey’s resistant line</h3>
-<p>The idea is simple in concept. It involves dividing the dataset into three approximately equal groups of values and summarizing these groups by computing their respective medians. Two half-slopes are then used to join the three points. Note that points that share the same x value are lumped into the same batch which can lead to unequal group sizes. The motivation behind this plot is to use the three-point summary to provide a robust assessment of the type of relationship between both variables. Such plots are often used to help guide re-expression of the variables x or y or both for the sole purpose of straightening the x-y relationship.</p>
-<p>Let’s look at an example using the last (modified) dataset. First, we divide the plot into three approximately equal batches.</p>
+<p>The idea is simple in concept. It involves dividing the dataset into three approximately equal groups of values and summarizing these groups by computing their respective medians. The batches’ median values are then used to compute the slope and intercept. The motivation behind this plot is to use the three-point summary to provide a robust assessment of the type of relationship between both variables.</p>
+<p>Let’s look at an example of its implementation. We’ll continue with the example used in the previous section. First, we divide the plot into three approximately equal batches.</p>
 <div class="cell" data-small.mar="true" data-hash="resistant_lines_cache/html/unnamed-chunk-7_b8f03c002c624a936661d3fc72e36a2d">
 <div class="cell-output-display">
 <p><img src="resistant_lines_files/figure-html/unnamed-chunk-7-1.png" class="img-fluid" width="288"></p>
 </div>
 </div>
-<p>Next, compute the median x and y values within each section.</p>
+<p>Next, we compute the median x and y values within each section.</p>
 <div class="cell" data-small.mar="true" data-hash="resistant_lines_cache/html/unnamed-chunk-8_3f38d7ceef20e0dfcde207f1ab4a141e">
 <div class="cell-output-display">
 <p><img src="resistant_lines_files/figure-html/unnamed-chunk-8-1.png" class="img-fluid" width="288"></p>
 </div>
 </div>
-<p>Note that we do not include the same point in any two median calculations. This implies that for the mid-third of the data, we do not include the point straddling the left boundary line when computing its median value. Likewise with the right-third of the data.</p>
 <p>The x median values are 2.5, 6, 9.5 and the y median values are 10, 17, 22 respectively.</p>
 <p>The two end medians are used to compute the slope as:</p>
 <p><span class="math display">\[
@@ -602,28 +601,28 @@ <h3 data-number="25.2.1" class="anchored" data-anchor-id="tukeys-resistant-line"
 </div>
 <p>This is not bad for a first iteration. But the outlier’s influence is still present given that the line is not passing through all nine points. Recall that the goal is to <em>not</em> allow one lone observation from influencing a fitted line.</p>
 <p>From this first iteration, we can extract the residuals. A line is then fitted to the residuals following the same procedure outlined above.</p>
-<div class="cell" data-small.mar="true" data-hash="resistant_lines_cache/html/unnamed-chunk-10_f05e07280210339891d5d97d22ddb550">
+<div class="cell" data-small.mar="true" data-hash="resistant_lines_cache/html/unnamed-chunk-10_41f87f969de578f7dc1e7a8ba640d2b3">
 <div class="cell-output-display">
 <p><img src="resistant_lines_files/figure-html/unnamed-chunk-10-1.png" class="img-fluid" width="288"></p>
 </div>
 </div>
-<p>The initial model slope and intercept are 6.43 and 1.71 respectively and the residual’s slope and intercept are -1.22 and 0.24 respectively. The residual slope is then added to the first computed slope and the process is again repeated thus generating the following tweaked slope and updated residuals:</p>
-<div class="cell" data-small.mar="true" data-hash="resistant_lines_cache/html/unnamed-chunk-11_b19fce4230a230eb339f211b686407eb">
+<p>The initial model intercept and slope are 6.429 and 1.714 respectively and the residual’s slope and intercept are -1.224 and 0.245 respectively. The residual slope is then added to the first computed slope giving us an updated slope of 1.959 and the process is again repeated thus generating the following tweaked slope and updated residuals:</p>
+<div class="cell" data-small.mar="true" data-hash="resistant_lines_cache/html/unnamed-chunk-11_7e538dd616e025fe6922c952147b2533">
 <div class="cell-output-display">
 <p><img src="resistant_lines_files/figure-html/unnamed-chunk-11-1.png" class="img-fluid" width="576"></p>
 </div>
 </div>
-<p>The updated slope is now 5.2. The iteration continues until the slope residuals stabilize. Each successive iteration is shown in the following plot with the final line in dark red. One more iteration should be enough to stabilize the residuals. Giving use the following fitted line whose intercept and slope are 5.03 and 1.99 respectively</p>
+<p>The second set of residuals is added to the previously computed slope giving us an updated slope of 1.994. The iteration continues until the slope residuals stabilize. Each successive iteration is shown in the following plot with the final line in dark red. This took three iterations giving us the final fitted line whose intercept and slope are 5.03 and 1.99 respectively–very close to the actual relationship between the two variables.</p>
 <div class="cell" data-small.mar="true" data-hash="resistant_lines_cache/html/unnamed-chunk-12_e85a23ab58acc055ba33325b3854ece6">
 <div class="cell-output-display">
 <p><img src="resistant_lines_files/figure-html/unnamed-chunk-12-1.png" class="img-fluid" width="288"></p>
 </div>
 </div>
-<p>Tukey’s resistant line can be computed using the R built-in <code>line()</code> function. You specify the number of iteration via the <code>iter</code> argument. Usually, five iterations will be enough.</p>
+<p>Tukey’s resistant line can be computed using the R built-in <code>line()</code> function. You specify the maximum number of iterations via the <code>iter</code> argument. Usually, five iterations will be enough. But, the iteration will stop before the maximum number of iterations is reached if the residual slope stabilizes before then.</p>
 <div class="cell" data-hash="resistant_lines_cache/html/unnamed-chunk-13_be2f218c284612df28118901bdba9dea">
 <div class="sourceCode cell-code" id="cb7"><pre class="sourceCode r code-with-copy"><code class="sourceCode r"><span id="cb7-1"><a href="#cb7-1" aria-hidden="true" tabindex="-1"></a>r.fit <span class="ot">&lt;-</span> <span class="fu">line</span>(x,y, <span class="at">iter =</span> <span class="dv">5</span>)</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
 </div>
-<p>You can the resistant to a base plot as follows:</p>
+<p>You can add the resistant line to a base plot as follows:</p>
 <div class="cell" data-hash="resistant_lines_cache/html/unnamed-chunk-14_613f81c207e185507c6d409afd4aae11">
 <div class="sourceCode cell-code" id="cb8"><pre class="sourceCode r code-with-copy"><code class="sourceCode r"><span id="cb8-1"><a href="#cb8-1" aria-hidden="true" tabindex="-1"></a><span class="fu">plot</span>(y <span class="sc">~</span> x, df)</span>
 <span id="cb8-2"><a href="#cb8-2" aria-hidden="true" tabindex="-1"></a><span class="fu">abline</span>(r.fit, <span class="at">col =</span> <span class="st">"red"</span>)</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
@@ -632,11 +631,12 @@ <h3 data-number="25.2.1" class="anchored" data-anchor-id="tukeys-resistant-line"
 </div>
 </div>
 <p>You can also add the line to a <code>ggplot</code>:</p>
-<div class="cell" data-hash="resistant_lines_cache/html/unnamed-chunk-15_7494a1816f080868e5824f192c16c90d">
-<div class="sourceCode cell-code" id="cb9"><pre class="sourceCode r code-with-copy"><code class="sourceCode r"><span id="cb9-1"><a href="#cb9-1" aria-hidden="true" tabindex="-1"></a><span class="fu">ggplot</span>(df, <span class="fu">aes</span>(x,y)) <span class="sc">+</span> <span class="fu">geom_point</span>() <span class="sc">+</span></span>
-<span id="cb9-2"><a href="#cb9-2" aria-hidden="true" tabindex="-1"></a>  <span class="fu">geom_abline</span>(<span class="at">intercept =</span> r.fit<span class="sc">$</span>coefficients[<span class="dv">1</span>], </span>
-<span id="cb9-3"><a href="#cb9-3" aria-hidden="true" tabindex="-1"></a>              <span class="at">slope =</span> r.fit<span class="sc">$</span>coefficients[<span class="dv">2</span>],</span>
-<span id="cb9-4"><a href="#cb9-4" aria-hidden="true" tabindex="-1"></a>              <span class="at">col =</span> <span class="st">"red"</span>)</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
+<div class="cell" data-hash="resistant_lines_cache/html/unnamed-chunk-15_c9fad914b2da127567b8036becee7f07">
+<div class="sourceCode cell-code" id="cb9"><pre class="sourceCode r code-with-copy"><code class="sourceCode r"><span id="cb9-1"><a href="#cb9-1" aria-hidden="true" tabindex="-1"></a><span class="fu">library</span>(ggplot2)</span>
+<span id="cb9-2"><a href="#cb9-2" aria-hidden="true" tabindex="-1"></a><span class="fu">ggplot</span>(df, <span class="fu">aes</span>(x,y)) <span class="sc">+</span> <span class="fu">geom_point</span>() <span class="sc">+</span></span>
+<span id="cb9-3"><a href="#cb9-3" aria-hidden="true" tabindex="-1"></a>  <span class="fu">geom_abline</span>(<span class="at">intercept =</span> r.fit<span class="sc">$</span>coefficients[<span class="dv">1</span>], </span>
+<span id="cb9-4"><a href="#cb9-4" aria-hidden="true" tabindex="-1"></a>              <span class="at">slope =</span> r.fit<span class="sc">$</span>coefficients[<span class="dv">2</span>],</span>
+<span id="cb9-5"><a href="#cb9-5" aria-hidden="true" tabindex="-1"></a>              <span class="at">col =</span> <span class="st">"red"</span>)</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
 <div class="cell-output-display">
 <p><img src="resistant_lines_files/figure-html/unnamed-chunk-15-1.png" class="img-fluid" width="288"></p>
 </div>
@@ -645,7 +645,7 @@ <h3 data-number="25.2.1" class="anchored" data-anchor-id="tukeys-resistant-line"
 <section id="bisquare" class="level3" data-number="25.2.2">
 <h3 data-number="25.2.2" class="anchored" data-anchor-id="bisquare"><span class="header-section-number">25.2.2</span> Bisquare</h3>
 <p>A bisquare fitting strategy makes use of a non-parametric model that assigns <em>weights</em> to observations based on their proximity to the center of the distribution. The closer an observation is to that center, the greater its weight.</p>
-<p>The first step is to run a linear regression model on the data then to extract its residuals. Next, a weighting scheme is fit to the residuals such that the points associated with extreme residuals are assigned the smallest weight and the points associated with the smaller residual values are assigned the largest weight. The regression analysis is then re-run using those same weights thus minimizing the influence of the rogue points. This process is repeated several times until the residuals stabilize. The following figure shows three iterations of the bisquare function whereby the weights (shown as grey text next to each point) start off as <code>1</code> then are modified following the residuals derived from the most recent regression model.</p>
+<p>The first step is to run a linear regression model on the data. The next step is to extract its residuals. Next, a weighting scheme is fitted to the residuals such that the points associated with extreme residuals are assigned the smallest weight and the points associated with the smaller residual values are assigned the largest weight. The regression analysis is then re-run using those same weights thus minimizing the influence of the rogue points. This process is repeated several times until the residuals stabilize. The following figure shows three iterations of the bisquare function whereby the weights (shown as grey text next to each point) start off as <code>1</code> then are modified following the residuals derived from the most recent regression model.</p>
 <div class="cell" data-small.mar="true" data-hash="resistant_lines_cache/html/unnamed-chunk-16_c70a75f3fef1561cb7ed8a2201a63220">
 <div class="cell-output-display">
 <p><img src="resistant_lines_files/figure-html/unnamed-chunk-16-1.png" class="img-fluid" width="768"></p>
@@ -677,9 +677,9 @@ <h3 data-number="25.2.2" class="anchored" data-anchor-id="bisquare"><span class=
 <p>In the above example, the bisquare method does a great job in eliminating the outlier’s influence.</p>
 <section id="built-in-implementation-of-the-bisquare" class="level4" data-number="25.2.2.1">
 <h4 data-number="25.2.2.1" class="anchored" data-anchor-id="built-in-implementation-of-the-bisquare"><span class="header-section-number">25.2.2.1</span> Built-in implementation of the bisquare</h4>
-<p>The <code>MASS</code> package has a robust linear modeling function called <code>rlm</code> that will implement a variation of the aforementioned bisquare estimation technique. Its results may differ slightly from those presented here, but the difference will be insignificant for the most part.</p>
-<p>Note that if you make use of <code>dplyr</code> in a workflow, loading <code>MASS</code> after <code>dplyr</code> will mask dplyr’s <code>select</code> function. This can be problematic. So you either want to load <code>MASS</code> before <code>dplyr</code>, or you can call the function via <code>MASS::rlm</code>. An example of its use follows. The default regression model is added as a dashed line for reference.</p>
-<div class="cell" data-small.mar="true" data-hash="resistant_lines_cache/html/unnamed-chunk-18_98214ea5eb9efd2489658138170116af">
+<p>The <code>MASS</code> package has a robust linear modeling function called <code>rlm</code> that will implement a variation of the aforementioned bisquare estimation technique. Its results may differ slightly from those presented above, but the difference will be insignificant for the most part.</p>
+<p>Note that if you make use of <code>dplyr</code> in a workflow, loading <code>MASS</code> after <code>dplyr</code> will mask dplyr’s <code>select</code> function. This can be problematic. So you either want to load <code>MASS</code> before <code>dplyr</code>, or you can call the function via <code>MASS::rlm</code>. An example of its use follows.</p>
+<div class="cell" data-small.mar="true" data-hash="resistant_lines_cache/html/unnamed-chunk-18_bd557fbeb46fffce48702dba906e7707">
 <div class="sourceCode cell-code" id="cb11"><pre class="sourceCode r code-with-copy"><code class="sourceCode r"><span id="cb11-1"><a href="#cb11-1" aria-hidden="true" tabindex="-1"></a>dat <span class="ot">&lt;-</span> <span class="fu">data.frame</span>(x,y)</span>
 <span id="cb11-2"><a href="#cb11-2" aria-hidden="true" tabindex="-1"></a>M   <span class="ot">&lt;-</span> <span class="fu">lm</span>( y <span class="sc">~</span> x, <span class="at">dat=</span>dat)</span>
 <span id="cb11-3"><a href="#cb11-3" aria-hidden="true" tabindex="-1"></a>Mr  <span class="ot">&lt;-</span> MASS<span class="sc">::</span><span class="fu">rlm</span>( y <span class="sc">~</span> x, <span class="at">dat=</span>dat, <span class="at">psi=</span><span class="st">"psi.bisquare"</span>)</span>
@@ -688,19 +688,19 @@ <h4 data-number="25.2.2.1" class="anchored" data-anchor-id="built-in-implementat
 <span id="cb11-6"><a href="#cb11-6" aria-hidden="true" tabindex="-1"></a><span class="co"># Add the robust line</span></span>
 <span id="cb11-7"><a href="#cb11-7" aria-hidden="true" tabindex="-1"></a><span class="fu">abline</span>( Mr, <span class="at">col=</span><span class="st">"red"</span>)</span>
 <span id="cb11-8"><a href="#cb11-8" aria-hidden="true" tabindex="-1"></a></span>
-<span id="cb11-9"><a href="#cb11-9" aria-hidden="true" tabindex="-1"></a><span class="co"># Add the default regression line</span></span>
+<span id="cb11-9"><a href="#cb11-9" aria-hidden="true" tabindex="-1"></a><span class="co"># Add the default regression line for reference</span></span>
 <span id="cb11-10"><a href="#cb11-10" aria-hidden="true" tabindex="-1"></a><span class="fu">abline</span>( M , <span class="at">col=</span><span class="st">"grey50"</span>, <span class="at">lty=</span><span class="dv">2</span>)</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
 <div class="cell-output-display">
 <p><img src="resistant_lines_files/figure-html/unnamed-chunk-18-1.png" class="img-fluid" width="288"></p>
 </div>
 </div>
+<p>The default regression model is added as a dashed line for reference.</p>
 <p>The function <code>rlm</code> can also be called directly from within <code>ggplot</code>.</p>
-<div class="cell" data-small.mar="true" data-hash="resistant_lines_cache/html/unnamed-chunk-19_5736d0c7f23d2dcbe1a6a29ebaca245d">
+<div class="cell" data-small.mar="true" data-hash="resistant_lines_cache/html/unnamed-chunk-19_48588591ed562775d36db603ab76dc0a">
 <div class="sourceCode cell-code" id="cb12"><pre class="sourceCode r code-with-copy"><code class="sourceCode r"><span id="cb12-1"><a href="#cb12-1" aria-hidden="true" tabindex="-1"></a><span class="fu">library</span>(ggplot2)</span>
 <span id="cb12-2"><a href="#cb12-2" aria-hidden="true" tabindex="-1"></a><span class="fu">ggplot</span>(dat, <span class="fu">aes</span>(<span class="at">x=</span>x, <span class="at">y=</span>y)) <span class="sc">+</span> <span class="fu">geom_point</span>() <span class="sc">+</span> </span>
-<span id="cb12-3"><a href="#cb12-3" aria-hidden="true" tabindex="-1"></a>       <span class="fu">stat_smooth</span>(<span class="at">method =</span> <span class="st">"lm"</span>, <span class="at">se =</span> <span class="cn">FALSE</span>, <span class="at">col =</span> <span class="st">"grey50"</span>, <span class="at">lty =</span> <span class="dv">2</span>) <span class="sc">+</span></span>
-<span id="cb12-4"><a href="#cb12-4" aria-hidden="true" tabindex="-1"></a>       <span class="fu">stat_smooth</span>(<span class="at">method =</span> MASS<span class="sc">::</span>rlm, <span class="at">se =</span> <span class="cn">FALSE</span>, <span class="at">col =</span> <span class="st">"red"</span>,</span>
-<span id="cb12-5"><a href="#cb12-5" aria-hidden="true" tabindex="-1"></a>                   <span class="at">method.args =</span> <span class="fu">list</span>(<span class="at">psi =</span> <span class="st">"psi.bisquare"</span>))</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
+<span id="cb12-3"><a href="#cb12-3" aria-hidden="true" tabindex="-1"></a>       <span class="fu">stat_smooth</span>(<span class="at">method =</span> MASS<span class="sc">::</span>rlm, <span class="at">se =</span> <span class="cn">FALSE</span>, <span class="at">col =</span> <span class="st">"red"</span>,</span>
+<span id="cb12-4"><a href="#cb12-4" aria-hidden="true" tabindex="-1"></a>                   <span class="at">method.args =</span> <span class="fu">list</span>(<span class="at">psi =</span> <span class="st">"psi.bisquare"</span>))</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
 <div class="cell-output-display">
 <p><img src="resistant_lines_files/figure-html/unnamed-chunk-19-1.png" class="img-fluid" width="288"></p>
 </div>
@@ -711,7 +711,7 @@ <h4 data-number="25.2.2.1" class="anchored" data-anchor-id="built-in-implementat
 </section>
 <section id="robust-loess" class="level2" data-number="25.3">
 <h2 data-number="25.3" class="anchored" data-anchor-id="robust-loess"><span class="header-section-number">25.3</span> Robust loess</h2>
-<p>The bisquare estimation method can also be extended to the loess smoothing function by passing the <code>"symmetric"</code> option to the <code>family</code> parameter.</p>
+<p>The loess fit is also not immune to outliers. The bisquare estimation method can be extended to the loess smoothing function by passing the <code>"symmetric"</code> option to the <code>family</code> argument.</p>
 <div class="cell" data-small.mar="true" data-hash="resistant_lines_cache/html/unnamed-chunk-20_a9b9ab144bdaec92a590b3ee3b197474">
 <div class="sourceCode cell-code" id="cb13"><pre class="sourceCode r code-with-copy"><code class="sourceCode r"><span id="cb13-1"><a href="#cb13-1" aria-hidden="true" tabindex="-1"></a><span class="co"># Fit a regular loess model</span></span>
 <span id="cb13-2"><a href="#cb13-2" aria-hidden="true" tabindex="-1"></a>lo <span class="ot">&lt;-</span> <span class="fu">loess</span>(y <span class="sc">~</span> x, dat)</span>
@@ -731,16 +731,14 @@ <h2 data-number="25.3" class="anchored" data-anchor-id="robust-loess"><span clas
 <p><img src="resistant_lines_files/figure-html/unnamed-chunk-20-1.png" class="img-fluid" width="288"></p>
 </div>
 </div>
-<p>The default loess fit is added as a dashed line for reference.</p>
+<p>The robust loess is in red, the default loess fit is added as a dashed line for reference.</p>
 <p>The robust <code>loess</code> can be implemented in <code>ggplot2</code>’s <code>stat_smooth</code> function as follows:</p>
-<div class="cell" data-small.mar="true" data-hash="resistant_lines_cache/html/unnamed-chunk-21_f34299528a71404319578f2d9f21372b">
+<div class="cell" data-small.mar="true" data-hash="resistant_lines_cache/html/unnamed-chunk-21_0a271fc46f5b247bd6ee9b36c233b494">
 <div class="sourceCode cell-code" id="cb14"><pre class="sourceCode r code-with-copy"><code class="sourceCode r"><span id="cb14-1"><a href="#cb14-1" aria-hidden="true" tabindex="-1"></a><span class="fu">library</span>(ggplot2)</span>
 <span id="cb14-2"><a href="#cb14-2" aria-hidden="true" tabindex="-1"></a><span class="fu">ggplot</span>(dat, <span class="fu">aes</span>(<span class="at">x=</span>x, <span class="at">y=</span>y)) <span class="sc">+</span> <span class="fu">geom_point</span>() <span class="sc">+</span> </span>
-<span id="cb14-3"><a href="#cb14-3" aria-hidden="true" tabindex="-1"></a>          <span class="fu">stat_smooth</span>(<span class="at">method =</span> <span class="st">"loess"</span>,                       <span class="co"># Regular loess </span></span>
-<span id="cb14-4"><a href="#cb14-4" aria-hidden="true" tabindex="-1"></a>                      <span class="at">se=</span><span class="cn">FALSE</span>, <span class="at">col=</span><span class="st">"grey50"</span>, <span class="at">lty=</span><span class="dv">2</span>)  <span class="sc">+</span>    </span>
-<span id="cb14-5"><a href="#cb14-5" aria-hidden="true" tabindex="-1"></a>          <span class="fu">stat_smooth</span>(<span class="at">method =</span> <span class="st">"loess"</span>, </span>
-<span id="cb14-6"><a href="#cb14-6" aria-hidden="true" tabindex="-1"></a>                      <span class="at">method.args =</span> <span class="fu">list</span>(<span class="at">family=</span><span class="st">"symmetric"</span>), <span class="co"># Robust loess</span></span>
-<span id="cb14-7"><a href="#cb14-7" aria-hidden="true" tabindex="-1"></a>                      <span class="at">se=</span><span class="cn">FALSE</span>,  <span class="at">col=</span><span class="st">"red"</span>) </span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
+<span id="cb14-3"><a href="#cb14-3" aria-hidden="true" tabindex="-1"></a>          <span class="fu">stat_smooth</span>(<span class="at">method =</span> <span class="st">"loess"</span>, </span>
+<span id="cb14-4"><a href="#cb14-4" aria-hidden="true" tabindex="-1"></a>                      <span class="at">method.args =</span> <span class="fu">list</span>(<span class="at">family=</span><span class="st">"symmetric"</span>), </span>
+<span id="cb14-5"><a href="#cb14-5" aria-hidden="true" tabindex="-1"></a>                      <span class="at">se=</span><span class="cn">FALSE</span>,  <span class="at">col=</span><span class="st">"red"</span>) </span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
 <div class="cell-output-display">
 <p><img src="resistant_lines_files/figure-html/unnamed-chunk-21-1.png" class="img-fluid" width="288"></p>
 </div>
diff --git a/docs/resistant_lines_files/figure-html/unnamed-chunk-10-1.png b/docs/resistant_lines_files/figure-html/unnamed-chunk-10-1.png
index 275a24ff..b6658b22 100644
Binary files a/docs/resistant_lines_files/figure-html/unnamed-chunk-10-1.png and b/docs/resistant_lines_files/figure-html/unnamed-chunk-10-1.png differ
diff --git a/docs/resistant_lines_files/figure-html/unnamed-chunk-11-1.png b/docs/resistant_lines_files/figure-html/unnamed-chunk-11-1.png
index a7a0c067..38d5d42b 100644
Binary files a/docs/resistant_lines_files/figure-html/unnamed-chunk-11-1.png and b/docs/resistant_lines_files/figure-html/unnamed-chunk-11-1.png differ
diff --git a/docs/resistant_lines_files/figure-html/unnamed-chunk-19-1.png b/docs/resistant_lines_files/figure-html/unnamed-chunk-19-1.png
index e31448cb..812b67ea 100644
Binary files a/docs/resistant_lines_files/figure-html/unnamed-chunk-19-1.png and b/docs/resistant_lines_files/figure-html/unnamed-chunk-19-1.png differ
diff --git a/docs/resistant_lines_files/figure-html/unnamed-chunk-21-1.png b/docs/resistant_lines_files/figure-html/unnamed-chunk-21-1.png
index 9547802f..812b67ea 100644
Binary files a/docs/resistant_lines_files/figure-html/unnamed-chunk-21-1.png and b/docs/resistant_lines_files/figure-html/unnamed-chunk-21-1.png differ
diff --git a/docs/search.json b/docs/search.json
index 32ad6397..11320cb4 100755
--- a/docs/search.json
+++ b/docs/search.json
@@ -4,7 +4,7 @@
     "href": "index.html",
     "title": "Exploratory Data Analysis in R",
     "section": "",
-    "text": "Preface\nThis book is a compilation of lecture notes used in an Exploratory Data Analysis in R course taught to undergraduates at Colby College. The course assumes little to no background in quantitative analysis nor in computer programming and was first taught in Spring, 2015. The course introduces students to data manipulation in R, data exploration (in the spirit of John Tukey’s EDA) and the R markdown language. Many of the visualization techniques are adopted from William Cleveland’s Data Visualization book.\nThe base R plotting environment and the ggplot2 ecosystem are used throughout this book. While a chapter is dedicated to the lattice plotting package, its functions are not used outside of that chapter given that ggplot2 offers many of lattice’s functionality.\nWhile great effort is made to adopt a consist plotting environment throughout this book (this being ggplot2, for the most part), a few topics (including the q-q plot and the median polish) will benefit from custom plotting functions available in the tukeyedar package. The package can be downloaded from GitHub via the command:\n\ndevtools::install_github(\"mgimond/tukeyedar\")\n\nNote that installing the GitHub package will require that the devtools package be installed first.\nFunctions making use of the tukeyedar package will be highlighted in a peach/pink code block as opposed to the default light yellow code block used for all other code blocks. For example, if tukeyedar’s eda_qq function is used, the code block will take on the following appearance:\n\nlibrary(tukeyedar)\neda_qq(Tenor, Bass)\n\nThe tukeyedar functions are built off of base R graphics and require R vesion 4.1 or greater.\n\nManuel “Manny” Gimond"
+    "text": "Preface\nThis book is a compilation of lecture notes used in an Exploratory Data Analysis in R course taught to undergraduates at Colby College. The course assumes little to no background in quantitative analysis nor in computer programming and was first taught in Spring, 2015. The course introduces students to data manipulation in R, data exploration (in the spirit of John Tukey’s EDA) and the R markdown language. Many of the visualization techniques are adopted from William Cleveland’s Data Visualization book.\nThe base R plotting environment and the ggplot2 ecosystem are used throughout this book. While a chapter is dedicated to the lattice plotting package, its functions are not used outside of that chapter given that ggplot2 offers many of lattice’s functionality.\nWhile great effort is made to adopt a consistent plotting environment throughout this book (this being ggplot2, for the most part), a few topics (including the q-q plot and the median polish) will benefit from custom plotting functions available in the tukeyedar package. The package can be downloaded from GitHub via the command:\n\ndevtools::install_github(\"mgimond/tukeyedar\")\n\nNote that installing the GitHub package will require that the devtools package be installed first.\nFunctions making use of the tukeyedar package will be highlighted in a peach/pink code block as opposed to the default light yellow code block used for all other code blocks. For example, if tukeyedar’s eda_qq function is used, the code block will take on the following appearance:\n\nlibrary(tukeyedar)\neda_qq(Tenor, Bass)\n\nThe tukeyedar functions are built off of base R graphics and require R vesion 4.1 or greater.\n\nManuel “Manny” Gimond"
   },
   {
     "objectID": "intro.html#what-is-exploratory-data-analysis-eda",
@@ -32,7 +32,7 @@
     "href": "intro.html#creating-dynamic-documents-using-r-markdown",
     "title": "1  EDA",
     "section": "1.4 Creating dynamic documents using R Markdown",
-    "text": "1.4 Creating dynamic documents using R Markdown\nAnother source of error in the write-up of a report or publication is the linking of tables, figures and statistical summaries to the write-up. Typically, one saves statistical plots as image files then loads the image into the document. However, the figures may have gone through many different iterations resulting in many different versions of the image file in a working folder. Add to this many other figures, data table files and statistical results from various pieces of software, one quickly realizes the potential for embedding the wrong image files in the document or embedding the wrong statistical summaries in the text. Furthermore, the researcher is then required to properly archive and document the provenance of each figure, data table or statistical summary resulting in a complex structure of files and directories in the project folder thus increasing the odds of an irreproducible analysis.\nConfining all of the analysis to a scripting environment such as R can help, but this still does not alleviate the possibility of loading the wrong figure into the document, or forgetting to update a statistical summary in the text when the original data file was revised. A solution to this potential pitfall is to embed the actual analysis and graphic generation process into the document–such environments are called dynamic documents. In this course, we will use the R Markdown authoring tool which embeds R code into the document. An example of an R Markdown document is this course website which was entirely generated in RMarkdown!. You can view the R Markdown files on this author’s GitHub repository."
+    "text": "1.4 Creating dynamic documents using R Markdown\nAnother source of error in the write-up of a report or publication is the linking of tables, figures and statistical summaries to the write-up. Typically, one saves statistical plots as image files then loads the image into the document. However, the figures may have gone through many different iterations resulting in many different versions of the image file in a working folder. Add to this many other figures, data table files and statistical results from various pieces of software, one quickly realizes the potential for embedding the wrong image files in the document or embedding the wrong statistical summaries in the text. Furthermore, the researcher is then required to properly archive and document the provenance of each figure, data table or statistical summary resulting in a complex structure of files and directories in the project folder thus increasing the odds of an irreproducible analysis.\nConfining all of the analysis to a scripting environment such as R can help, but this still does not alleviate the possibility of loading the wrong figure into the document, or forgetting to update a statistical summary in the text when the original data file was revised. A solution to this potential pitfall is to embed the actual analysis and graphic generation process into the document–such environments are called dynamic documents. In this course, we will use the R Markdown authoring tool which embeds R code into the document. An example of an R Markdown document is this course website which was entirely generated in RMarkdown! You can view the R Markdown files on this author’s GitHub repository."
   },
   {
     "objectID": "intro.html#footnotes",
@@ -886,21 +886,21 @@
     "href": "resistant_lines.html#introduction",
     "title": "25  Resistant lines",
     "section": "25.1 Introduction",
-    "text": "25.1 Introduction\nOrdinary least squares regression lines (those created with the lm() function) suffer from sensitivity to outliers. Because lm’s best fit line makes use of the mean (which is not a robust measure of location), its breakdown point is \\(1/n\\), meaning that all it takes is for one data point to behave differently from the rest of the points to significantly alter the slope and intercept of the best fit line. For example, let’s start with a well behaved dataset where all points are perfectly aligned, and fit this batch with a regression line:\n\nx &lt;- seq(1:11)\ny &lt;- 5 + 2 * x\nplot(y ~ x, pch = 20)\nM &lt;- lm(y ~ x)\nabline( M , col = \"red\")\n\n\n\n\nAs expected, we have a perfect fit. And the regression model’s coefficients match those used to create the data.\n\ncoef(M)\n\n(Intercept)           x \n          5           2 \n\n\nNow, what if one of the points is re-positioned in the plot, what happens to the regression line?\n\ny[11] &lt;- 2\nplot(y ~ x, pch = 20)\nM &lt;- lm(y ~ x)\nabline( M , col = \"red\")\n\n\n\n\nNote the significant change in the line’s characteristics, its intercept and slope are now:\n\ncoef(M)\n\n(Intercept)           x \n  9.5454545   0.8636364 \n\n\nThe slope dropped from 2 to 0.86 because of a single point!\nIf our goal is to explore what the bulk of the data has to say about the phenomena being investigated, we certainly don’t want a few “maverick” values to hijack the analysis. We therefore need a set of fitting tools that minimize the influence of outliers. There are many options out there; most notable are the resistant line (that builds off of Tukey’s 3-point summary line) and the bisquare robust estimation method."
+    "text": "25.1 Introduction\nOrdinary least squares regression lines (those created with the lm() function) suffer from sensitivity to outliers. Because lm’s best fit line makes use of the mean (which is not a robust measure of location), its breakdown point is \\(1/n\\), meaning that all it takes is for one data point to behave differently from the rest of the points to significantly alter the slope and intercept of the best fit line. For example, let’s start with a well behaved dataset where all points are perfectly aligned, and fit this batch with a regression line:\n\nx &lt;- seq(1:11)\ny &lt;- 5 + 2 * x\nplot(y ~ x, pch = 20)\nM &lt;- lm(y ~ x)\nabline( M , col = \"red\")\n\n\n\n\nAs expected, we have a perfect fit. And the regression model’s coefficients match those used to create the data.\n\ncoef(M)\n\n(Intercept)           x \n          5           2 \n\n\nNow, what if one of the points is re-positioned in the plot, what happens to the regression line?\n\ny[11] &lt;- 2\nplot(y ~ x, pch = 20)\nM &lt;- lm(y ~ x)\nabline( M , col = \"red\")\n\n\n\n\nNote the significant change in the line’s characteristics, its intercept and slope are now:\n\ncoef(M)\n\n(Intercept)           x \n  9.5454545   0.8636364 \n\n\nThe slope dropped from 2 to 0.86 because of a single point!\nIf our goal is to explore what the bulk of the data has to say about the phenomena being investigated, we certainly don’t want a few “maverick” values to hijack the analysis. We therefore need a set of fitting tools that minimize the influence of outliers. There are many options out there; most notable are Tukey’s resistant line and the bisquare robust estimation method."
   },
   {
     "objectID": "resistant_lines.html#robust-lines",
     "href": "resistant_lines.html#robust-lines",
     "title": "25  Resistant lines",
     "section": "25.2 Robust lines",
-    "text": "25.2 Robust lines\n  \n\n25.2.1 Tukey’s resistant line\nThe idea is simple in concept. It involves dividing the dataset into three approximately equal groups of values and summarizing these groups by computing their respective medians. Two half-slopes are then used to join the three points. Note that points that share the same x value are lumped into the same batch which can lead to unequal group sizes. The motivation behind this plot is to use the three-point summary to provide a robust assessment of the type of relationship between both variables. Such plots are often used to help guide re-expression of the variables x or y or both for the sole purpose of straightening the x-y relationship.\nLet’s look at an example using the last (modified) dataset. First, we divide the plot into three approximately equal batches.\n\n\n\n\n\nNext, compute the median x and y values within each section.\n\n\n\n\n\nNote that we do not include the same point in any two median calculations. This implies that for the mid-third of the data, we do not include the point straddling the left boundary line when computing its median value. Likewise with the right-third of the data.\nThe x median values are 2.5, 6, 9.5 and the y median values are 10, 17, 22 respectively.\nThe two end medians are used to compute the slope as:\n\\[\nb = \\frac{y_r - y_l}{x_r-x_l}\n\\] where the subscripts \\(r\\) and \\(l\\) reference the median values for the right-most and left-most batches.\nOnce the slope is computed, the intercept can be computed as follows:\n\\[\nmedian(y_{l,m,r} - b * x_{l,m,r})\n\\]\nwhere \\((x,y)_{l,m,r}\\) are the median x and y values for each batch.\n\n\n\n\n\nThis is not bad for a first iteration. But the outlier’s influence is still present given that the line is not passing through all nine points. Recall that the goal is to not allow one lone observation from influencing a fitted line.\nFrom this first iteration, we can extract the residuals. A line is then fitted to the residuals following the same procedure outlined above.\n\n\n\n\n\nThe initial model slope and intercept are 6.43 and 1.71 respectively and the residual’s slope and intercept are -1.22 and 0.24 respectively. The residual slope is then added to the first computed slope and the process is again repeated thus generating the following tweaked slope and updated residuals:\n\n\n\n\n\nThe updated slope is now 5.2. The iteration continues until the slope residuals stabilize. Each successive iteration is shown in the following plot with the final line in dark red. One more iteration should be enough to stabilize the residuals. Giving use the following fitted line whose intercept and slope are 5.03 and 1.99 respectively\n\n\n\n\n\nTukey’s resistant line can be computed using the R built-in line() function. You specify the number of iteration via the iter argument. Usually, five iterations will be enough.\n\nr.fit &lt;- line(x,y, iter = 5)\n\nYou can the resistant to a base plot as follows:\n\nplot(y ~ x, df)\nabline(r.fit, col = \"red\")\n\n\n\n\nYou can also add the line to a ggplot:\n\nggplot(df, aes(x,y)) + geom_point() +\n  geom_abline(intercept = r.fit$coefficients[1], \n              slope = r.fit$coefficients[2],\n              col = \"red\")\n\n\n\n\n\n\n25.2.2 Bisquare\nA bisquare fitting strategy makes use of a non-parametric model that assigns weights to observations based on their proximity to the center of the distribution. The closer an observation is to that center, the greater its weight.\nThe first step is to run a linear regression model on the data then to extract its residuals. Next, a weighting scheme is fit to the residuals such that the points associated with extreme residuals are assigned the smallest weight and the points associated with the smaller residual values are assigned the largest weight. The regression analysis is then re-run using those same weights thus minimizing the influence of the rogue points. This process is repeated several times until the residuals stabilize. The following figure shows three iterations of the bisquare function whereby the weights (shown as grey text next to each point) start off as 1 then are modified following the residuals derived from the most recent regression model.\n\n\n\n\n\nThere are different weighing strategies that can be implemented. One such implementation is presented next.\n\n# Create the bisquare function\nwt.bisquare &lt;- function(u, c = 6) {\n   ifelse( abs(u/c) &lt; 1, (1-(u/c)^2)^2, 0)\n}\n\n# Assign an equal weight to all points\nwt &lt;- rep(1, length(x))\n\n# Compute the regression, then assign weights based on residual values\nfor(i in 1:10){\n  dat.lm &lt;- lm(y ~ x ,weights=wt)\n  wt &lt;- wt.bisquare( dat.lm$res/ median(abs(dat.lm$res)), c = 6 )\n}\n\n# Plot the data and the resulting line\nplot(x, y, pch = 20)\nabline(dat.lm, col = rgb(1,0,0,0.3))\n\n\n\n\nIn the above example, the bisquare method does a great job in eliminating the outlier’s influence.\n\n25.2.2.1 Built-in implementation of the bisquare\nThe MASS package has a robust linear modeling function called rlm that will implement a variation of the aforementioned bisquare estimation technique. Its results may differ slightly from those presented here, but the difference will be insignificant for the most part.\nNote that if you make use of dplyr in a workflow, loading MASS after dplyr will mask dplyr’s select function. This can be problematic. So you either want to load MASS before dplyr, or you can call the function via MASS::rlm. An example of its use follows. The default regression model is added as a dashed line for reference.\n\ndat &lt;- data.frame(x,y)\nM   &lt;- lm( y ~ x, dat=dat)\nMr  &lt;- MASS::rlm( y ~ x, dat=dat, psi=\"psi.bisquare\")\nplot(y ~ x,dat, pch=16, col=rgb(0,0,0,0.2))\n\n# Add the robust line\nabline( Mr, col=\"red\")\n\n# Add the default regression line\nabline( M , col=\"grey50\", lty=2)\n\n\n\n\nThe function rlm can also be called directly from within ggplot.\n\nlibrary(ggplot2)\nggplot(dat, aes(x=x, y=y)) + geom_point() + \n       stat_smooth(method = \"lm\", se = FALSE, col = \"grey50\", lty = 2) +\n       stat_smooth(method = MASS::rlm, se = FALSE, col = \"red\",\n                   method.args = list(psi = \"psi.bisquare\"))"
+    "text": "25.2 Robust lines\n  \n\n25.2.1 Tukey’s resistant line\nThe idea is simple in concept. It involves dividing the dataset into three approximately equal groups of values and summarizing these groups by computing their respective medians. The batches’ median values are then used to compute the slope and intercept. The motivation behind this plot is to use the three-point summary to provide a robust assessment of the type of relationship between both variables.\nLet’s look at an example of its implementation. We’ll continue with the example used in the previous section. First, we divide the plot into three approximately equal batches.\n\n\n\n\n\nNext, we compute the median x and y values within each section.\n\n\n\n\n\nThe x median values are 2.5, 6, 9.5 and the y median values are 10, 17, 22 respectively.\nThe two end medians are used to compute the slope as:\n\\[\nb = \\frac{y_r - y_l}{x_r-x_l}\n\\] where the subscripts \\(r\\) and \\(l\\) reference the median values for the right-most and left-most batches.\nOnce the slope is computed, the intercept can be computed as follows:\n\\[\nmedian(y_{l,m,r} - b * x_{l,m,r})\n\\]\nwhere \\((x,y)_{l,m,r}\\) are the median x and y values for each batch.\n\n\n\n\n\nThis is not bad for a first iteration. But the outlier’s influence is still present given that the line is not passing through all nine points. Recall that the goal is to not allow one lone observation from influencing a fitted line.\nFrom this first iteration, we can extract the residuals. A line is then fitted to the residuals following the same procedure outlined above.\n\n\n\n\n\nThe initial model intercept and slope are 6.429 and 1.714 respectively and the residual’s slope and intercept are -1.224 and 0.245 respectively. The residual slope is then added to the first computed slope giving us an updated slope of 1.959 and the process is again repeated thus generating the following tweaked slope and updated residuals:\n\n\n\n\n\nThe second set of residuals is added to the previously computed slope giving us an updated slope of 1.994. The iteration continues until the slope residuals stabilize. Each successive iteration is shown in the following plot with the final line in dark red. This took three iterations giving us the final fitted line whose intercept and slope are 5.03 and 1.99 respectively–very close to the actual relationship between the two variables.\n\n\n\n\n\nTukey’s resistant line can be computed using the R built-in line() function. You specify the maximum number of iterations via the iter argument. Usually, five iterations will be enough. But, the iteration will stop before the maximum number of iterations is reached if the residual slope stabilizes before then.\n\nr.fit &lt;- line(x,y, iter = 5)\n\nYou can add the resistant line to a base plot as follows:\n\nplot(y ~ x, df)\nabline(r.fit, col = \"red\")\n\n\n\n\nYou can also add the line to a ggplot:\n\nlibrary(ggplot2)\nggplot(df, aes(x,y)) + geom_point() +\n  geom_abline(intercept = r.fit$coefficients[1], \n              slope = r.fit$coefficients[2],\n              col = \"red\")\n\n\n\n\n\n\n25.2.2 Bisquare\nA bisquare fitting strategy makes use of a non-parametric model that assigns weights to observations based on their proximity to the center of the distribution. The closer an observation is to that center, the greater its weight.\nThe first step is to run a linear regression model on the data. The next step is to extract its residuals. Next, a weighting scheme is fitted to the residuals such that the points associated with extreme residuals are assigned the smallest weight and the points associated with the smaller residual values are assigned the largest weight. The regression analysis is then re-run using those same weights thus minimizing the influence of the rogue points. This process is repeated several times until the residuals stabilize. The following figure shows three iterations of the bisquare function whereby the weights (shown as grey text next to each point) start off as 1 then are modified following the residuals derived from the most recent regression model.\n\n\n\n\n\nThere are different weighing strategies that can be implemented. One such implementation is presented next.\n\n# Create the bisquare function\nwt.bisquare &lt;- function(u, c = 6) {\n   ifelse( abs(u/c) &lt; 1, (1-(u/c)^2)^2, 0)\n}\n\n# Assign an equal weight to all points\nwt &lt;- rep(1, length(x))\n\n# Compute the regression, then assign weights based on residual values\nfor(i in 1:10){\n  dat.lm &lt;- lm(y ~ x ,weights=wt)\n  wt &lt;- wt.bisquare( dat.lm$res/ median(abs(dat.lm$res)), c = 6 )\n}\n\n# Plot the data and the resulting line\nplot(x, y, pch = 20)\nabline(dat.lm, col = rgb(1,0,0,0.3))\n\n\n\n\nIn the above example, the bisquare method does a great job in eliminating the outlier’s influence.\n\n25.2.2.1 Built-in implementation of the bisquare\nThe MASS package has a robust linear modeling function called rlm that will implement a variation of the aforementioned bisquare estimation technique. Its results may differ slightly from those presented above, but the difference will be insignificant for the most part.\nNote that if you make use of dplyr in a workflow, loading MASS after dplyr will mask dplyr’s select function. This can be problematic. So you either want to load MASS before dplyr, or you can call the function via MASS::rlm. An example of its use follows.\n\ndat &lt;- data.frame(x,y)\nM   &lt;- lm( y ~ x, dat=dat)\nMr  &lt;- MASS::rlm( y ~ x, dat=dat, psi=\"psi.bisquare\")\nplot(y ~ x,dat, pch=16, col=rgb(0,0,0,0.2))\n\n# Add the robust line\nabline( Mr, col=\"red\")\n\n# Add the default regression line for reference\nabline( M , col=\"grey50\", lty=2)\n\n\n\n\nThe default regression model is added as a dashed line for reference.\nThe function rlm can also be called directly from within ggplot.\n\nlibrary(ggplot2)\nggplot(dat, aes(x=x, y=y)) + geom_point() + \n       stat_smooth(method = MASS::rlm, se = FALSE, col = \"red\",\n                   method.args = list(psi = \"psi.bisquare\"))"
   },
   {
     "objectID": "resistant_lines.html#robust-loess",
     "href": "resistant_lines.html#robust-loess",
     "title": "25  Resistant lines",
     "section": "25.3 Robust loess",
-    "text": "25.3 Robust loess\nThe bisquare estimation method can also be extended to the loess smoothing function by passing the \"symmetric\" option to the family parameter.\n\n# Fit a regular loess model\nlo &lt;- loess(y ~ x, dat)\n\n# Fit a robust loess model\nlor &lt;- loess(y ~ x, dat, family = \"symmetric\")\n\n# Plot the data\nplot(y ~ x, dat, pch = 16, col = rgb(0,0,0,0.2))\n\n# Add the default loess\nlines(dat$x, predict(lo), col = \"grey50\", lty = 2)\n\n# Add the robust loess\nlines(dat$x, predict(lor), col = \"red\")\n\n\n\n\nThe default loess fit is added as a dashed line for reference.\nThe robust loess can be implemented in ggplot2’s stat_smooth function as follows:\n\nlibrary(ggplot2)\nggplot(dat, aes(x=x, y=y)) + geom_point() + \n          stat_smooth(method = \"loess\",                       # Regular loess \n                      se=FALSE, col=\"grey50\", lty=2)  +    \n          stat_smooth(method = \"loess\", \n                      method.args = list(family=\"symmetric\"), # Robust loess\n                      se=FALSE,  col=\"red\") \n\n\n\n\nThe loess parameters are passed to the loess function via the call to method.args. The red curve is that of the robust loess and the grey dashed curve is that of the standard loess."
+    "text": "25.3 Robust loess\nThe loess fit is also not immune to outliers. The bisquare estimation method can be extended to the loess smoothing function by passing the \"symmetric\" option to the family argument.\n\n# Fit a regular loess model\nlo &lt;- loess(y ~ x, dat)\n\n# Fit a robust loess model\nlor &lt;- loess(y ~ x, dat, family = \"symmetric\")\n\n# Plot the data\nplot(y ~ x, dat, pch = 16, col = rgb(0,0,0,0.2))\n\n# Add the default loess\nlines(dat$x, predict(lo), col = \"grey50\", lty = 2)\n\n# Add the robust loess\nlines(dat$x, predict(lor), col = \"red\")\n\n\n\n\nThe robust loess is in red, the default loess fit is added as a dashed line for reference.\nThe robust loess can be implemented in ggplot2’s stat_smooth function as follows:\n\nlibrary(ggplot2)\nggplot(dat, aes(x=x, y=y)) + geom_point() + \n          stat_smooth(method = \"loess\", \n                      method.args = list(family=\"symmetric\"), \n                      se=FALSE,  col=\"red\") \n\n\n\n\nThe loess parameters are passed to the loess function via the call to method.args. The red curve is that of the robust loess and the grey dashed curve is that of the standard loess."
   },
   {
     "objectID": "example.html#an-eda-example-co2-analysis",
@@ -956,7 +956,7 @@
     "href": "two_way.html#introduction",
     "title": "28  Analyzing two-way tables: The median polish",
     "section": "28.1 Introduction",
-    "text": "28.1 Introduction\nBefore tackling two-way tables, let’s explore a simpler one-way table. The following dataset shows mean annual income by occupation for the year 2014 (src: Bureau of Labor Statistics)\n\n\n\n\n\n\n\n\n\n\nIncome\n\n\n\n\nArchitecture and Engineering\n\n\n112490\n\n\n\n\nArts, Design, Entertainment, Sports, and Media\n\n\n72410\n\n\n\n\nBuilding and Grounds Cleaning and Maintenance\n\n\n83970\n\n\n\n\nBusiness and Financial Operations\n\n\n81520\n\n\n\n\nCommunity and Social Service\n\n\n70070\n\n\n\n\nComputer and Mathematical\n\n\n45310\n\n\n\n\nConstruction and Extraction\n\n\n101110\n\n\n\n\nEducation, Training, and Library\n\n\n52210\n\n\n\n\nFarming, Fishing, and Forestry\n\n\n55790\n\n\n\n\nFood Preparation and Serving Related\n\n\n76010\n\n\n\n\nHealthcare Practitioners and Technical\n\n\n28820\n\n\n\n\nHealthcare Support\n\n\n43980\n\n\n\n\nInstallation, Maintenance, and Repair\n\n\n21980\n\n\n\n\nLegal\n\n\n26370\n\n\n\n\nLife, Physical, and Social Science\n\n\n24980\n\n\n\n\nManagement\n\n\n38660\n\n\n\n\nOffice and Administrative Support\n\n\n35530\n\n\n\n\nPersonal Care and Service\n\n\n25160\n\n\n\n\nProduction\n\n\n46600\n\n\n\n\nProtective Service\n\n\n45220\n\n\n\n\nSales and Related\n\n\n35490\n\n\n\n\nTransportation and Material Moving\n\n\n34460\n\n\n\n\n\nSo, what can we glean from this dataset other than the raw values associated with each occupation type? For starters, we may be interested in coming up with a summary of income values such as the median, or $45265 in our working example. This is a single value that can be used to characterize the overall income value. However, the income values in our batch are not all equal to $45265–they vary above or below this median value. A good way to measure the variability about the median is to subtract the median income value from each occupation income value. We can do this in a table, but will choose a dotplot chart instead to facilitate comparison.\n\n\n\n\n\nWe have separated the income values into two parts: an overall value and the difference between each value and the overall–the residuals. This can be written in the following algebraic form:\n\\[\nIncome_{Occupation} = Income_{median} + Income_{residual}\n\\]\nFor example, ‘Healthcare Support’ average income of $43980 can be separated into the the median, $45265, and the residual, $-1285.\nSo in addition to gleaning a median income value from this batch of values, we have also devised a way to facilitate the comparison of the income values to one another. Alternatively, we could have assigned colors to each cell of the table reflecting their values relative to the overall median.\n\n\n\n\n\nThe income value can be thought of as the response to the occupation type (also known as a factor). In other words, the factor determines the value in our table. We can re-write the last equation in a more formal way:\n\\[\ny_{i} = \\mu + \\epsilon_i\n\\]\nwhere \\(y_{i}\\) is the response variable (e.g. income) for factor \\(i\\) (e.g. occupation), \\(\\mu\\) is the overall typical value (hereafter referred to as the common value) and is usually the mean or median, and \\(\\epsilon_i\\) is the residual. Next, we will see how the idea of separating the values into components can be expanded to two-way (aka two-factor) tables."
+    "text": "28.1 Introduction\nBefore tackling two-way tables, let’s explore a simpler one-way table. The following dataset shows mean annual income by occupation for the year 2014 (src: Bureau of Labor Statistics)\n\n\n\n\n\n\n\n\n\n\nIncome\n\n\n\n\nArchitecture and Engineering\n\n\n112490\n\n\n\n\nArts, Design, Entertainment, Sports, and Media\n\n\n72410\n\n\n\n\nBuilding and Grounds Cleaning and Maintenance\n\n\n83970\n\n\n\n\nBusiness and Financial Operations\n\n\n81520\n\n\n\n\nCommunity and Social Service\n\n\n70070\n\n\n\n\nComputer and Mathematical\n\n\n45310\n\n\n\n\nConstruction and Extraction\n\n\n101110\n\n\n\n\nEducation, Training, and Library\n\n\n52210\n\n\n\n\nFarming, Fishing, and Forestry\n\n\n55790\n\n\n\n\nFood Preparation and Serving Related\n\n\n76010\n\n\n\n\nHealthcare Practitioners and Technical\n\n\n28820\n\n\n\n\nHealthcare Support\n\n\n43980\n\n\n\n\nInstallation, Maintenance, and Repair\n\n\n21980\n\n\n\n\nLegal\n\n\n26370\n\n\n\n\nLife, Physical, and Social Science\n\n\n24980\n\n\n\n\nManagement\n\n\n38660\n\n\n\n\nOffice and Administrative Support\n\n\n35530\n\n\n\n\nPersonal Care and Service\n\n\n25160\n\n\n\n\nProduction\n\n\n46600\n\n\n\n\nProtective Service\n\n\n45220\n\n\n\n\nSales and Related\n\n\n35490\n\n\n\n\nTransportation and Material Moving\n\n\n34460\n\n\n\n\n\nSo, what can we glean from this dataset other than the raw values associated with each occupation type? For starters, we may be interested in coming up with a summary of income values such as the median, or $45265 in our working example. This is a single value that can be used to characterize the overall income value. However, the income values in our batch are not all equal to $45265–they vary above or below this median value. A good way to measure the variability about the median is to subtract the median income value from each occupation income value. We can do this in a table, but will choose a dotplot chart instead to facilitate comparison.\n\n\n\n\n\nWe have separated the income values into two parts: an overall value and the difference between each value and the overall–i.e. the residuals. This can be written in the following algebraic form:\n\\[\nIncome_{Occupation} = Income_{median} + Income_{residual}\n\\]\nFor example, ‘Healthcare Support’ average income of $43980 can be separated into the median, $45265, and its residual, $-1285.\nSo, in addition to gleaning a median income value from this batch of values, we have also devised a way to facilitate the comparison of the income values to one another. Alternatively, we could have assigned colors to each cell of the table reflecting their values relative to the overall median.\n\n\n\n\n\nThe income value can be thought of as the response to the occupation type (also known as a factor). In other words, the factor determines the value in our table. We can re-write the last equation in a more formal way:\n\\[\ny_{i} = \\mu + \\epsilon_i\n\\]\nwhere \\(y_{i}\\) is the response variable (e.g. income) for factor \\(i\\) (e.g. occupation), \\(\\mu\\) is the overall typical value (hereafter referred to as the common value) and is usually the mean or median, and \\(\\epsilon_i\\) is the residual. Next, we will see how the idea of separating the values into components can be expanded to two-way (aka two-factor) tables."
   },
   {
     "objectID": "two_way.html#anatomy-of-a-two-way-table",
@@ -970,21 +970,21 @@
     "href": "two_way.html#analysis-workflow",
     "title": "28  Analyzing two-way tables: The median polish",
     "section": "28.3 Analysis workflow",
-    "text": "28.3 Analysis workflow\nLet’s first create the dataframe.\n\ndf &lt;- data.frame(row.names = c(\"NE\",\"NC\",\"S\",\"W\"),\n                 ed8       = c(25.3,32.1,38.8,25.4), \n                 ed9to11   = c(25.3,29,31,21.1),\n                 ed12      = c(18.2,18.8,19.3,20.3),\n                 ed13to15  = c(18.3,24.3,15.7,24),\n                 ed16      = c(16.3,19,16.8,17.5)\n                 )\n\n\n28.3.1 Visualizing the data\nIt’s often easier to look at a graphical representation of the data than a tabular one. Even a table as small as this can benefit from a plot.\nWe will adopt Bill Cleveland’s dotplot for this purpose. R has a built-in dotplot function called dotchart. It requires that the table be stored as a matrix instead of a dataframe; we will therefore convert df to a matrix by wrapping it with the as.matrix() function.\n\ndotchart( as.matrix(df), cex=0.7)\n\n\n\n\nThe plot helps visualize any differences in mortality rates across different father educational attainment levels. There seems to be a gradual decrease in child mortality with increasing father educational attainment.\nBut the plot does not help spot differences across regions (except for the ed12 group). We can generate another plot where region becomes the main grouping factor. We do this by wrapping the matrix with the transpose function t().\n\ndotchart( t(as.matrix(df)), cex=0.7)\n\n\n\n\nAt first glance, there seems to be higher death rates for the NC and S regions and relatively lower rates for the W and NE regions. But our eyes may be fooled by outliers in the data.\nNext, we’ll generate side-by-side boxplots to compare the effects between both categories. Note that we’ll need to create a long version of the table using the tidyr package.\n\nlibrary(tidyr)\nlibrary(dplyr)\n\ndf.l &lt;- df %&gt;%\n  mutate(Region = as.factor(row.names(.)) )  %&gt;%\n  pivot_longer(names_to = \"Edu\", values_to = \"Deaths\", -Region) %&gt;% \n  mutate(Edu = factor(Edu, levels=names(df)))\n\n# side-by-side plot\nOP &lt;- par(mfrow=c(1,2))\nplot(Deaths ~ Region + Edu, df.l)\npar(OP)\n\n\n\n\nSo at this point the plots suggest that there may be a father’s educational attainment effect as well as a regional effect on infant mortality with the former effect being possibly more important.\nNext we will attempt to quantify these effects using Tukey’s median polish.\n\n\n28.3.2 Median polish\nA median polish is an iterative process whereby we seek to extract the effects from the data into their row and column components. The steps involved in each iteration are described next.\n\n28.3.2.1 Step 1: Compute overall median and residual table\nFirst, we will compute the median value for all values in the dataset; this will be our first estimate of the common value. The resulting median value is placed in the upper left-hand margin of the table. A residual table is created by taking the difference between the original value and the overall median.\n\nYou’ll also note that the row and column effect values are initially populated with 0 values.\n\n\n28.3.2.2 Step 2: Compute the row medians\nNext, row medians are computed (shown in red in the following figure) from the residual table.\n\nNote that in addition to computing the row median values for the response variable, we are also computing the median value for the column effects in the upper margin (which happens to be zero since we have no column effect values yet).\n\n\n28.3.2.3 Step 3: Create a new residual table from the row medians\nThe row medians are added to the left hand margin of the new residual table (shown in green). These values represent the row effects in this first iteration. A new set of residual values is created from the row medians where each cell takes on the value of the subtraction of the row median from each response variable in that row.\n\nThe subtraction is also done for the column effect values (even though all values remain zero) and the global median (common effect value).\n\n\n28.3.2.4 Step 4: Compute the column medians\nNext, column medians are computed (shown in red in the following figure) from the residual table. Note that you are NOT computing the column medians from the original response variables. Note too that we are also computing the median value for the row effect values (left margin).\n\n\n\n28.3.2.5 Step 5: Create a new residual table from the column medians\nA new residual table is created from the column medians where each new cell takes on the value of the subtraction of the column median from each initial residual value in that column. For example, the first upper-left cell’s residual is \\(7.0 - 7.4 = -0.4\\).\nThe column medians are also added to the column effect margin. Note that we are also adding the row effect median value, -0.5, to the common effect cell in the upper left-hand corner.\n\nWe have now completed our first row and column smoothing iteration. However, there may be more row and column effects that can be pulled from the residuals. We therefore move on to a second iteration.\n\n\n28.3.2.6 Step 6: Second iteration – row effects\nNext, we compute the row medians from the residuals, then add the column effect median to the top margin (and the common value) and subtract the row medians from the residuals. as in step 3.\n\n\n\n28.3.2.7 Step 7: Second iteration – column effects\nTo wrap-up the second iteration, we compute the column median values from the residuals then subtract the medians from those residuals. We also add the medians to the row effect values and the common value.\n\n\n\n28.3.2.8 When do we stop iterating?\nThe goal is to iterate through the row and column smoothing operations until the row and column effect medians are close to 0. However, Hoaglin et al. (1983) warn against “using unrestrained iteration” and suggest that a few steps should be more than adequate in most instances. In our working example, a third iteration may be warranted. The complete workflow for this third iteration follows.\n\nThe final version of the table (along with the column and row values) is shown below:\n\nA progression of each iteration is shown side-by-side in the next figure. Note how, at the end of the iterations, the common factor ends up contributing the largest effect. Also note the one abnormally high residual in the southern region associated with 8th grade education.\n\n\n\n\n\nThe left-most table is the original data showing death rates. The second table shows the outcome of the first round of polishing (including the initial overall median value of 20.2). The third and forth table show the second and third iterations of the smoothing operations. Additional iterations are not deemed necessary given that little more can be extracted from the residuals.\n\n\n\n28.3.3 Interpreting the median polish\nAs noted earlier, a two-way table represents the relationship between the response variable, \\(y\\), and the two categories as:\n\\[\ny_{ij} = \\mu + \\alpha_i + \\beta_j + \\epsilon_{ij}\n\\]\nIn our working example, \\(\\mu\\) = 20.6; \\(\\alpha_i\\) = -1.5, 2.6 -0.4, and 0.4 for \\(i\\) = NE, Nc, S and W respectively; \\(\\beta_j\\) = 7.6, 6.0, -0.9, 0.2 and -3.5 for \\(j\\) = ed8, ed9to11, ed12, ed13to15 and ed16 respectively.\nThe residuals represent the portion of the mortality rates that can’t be explained by either factors.\nSo the mortality rate in the upper left-hand cell from the original table can be deconstructed as:\n\nThe examination of the table suggests that the infant mortality rate is greatest for fathers who did not attain more than 8 years of school (i.e. who has not completed high school) as noted by the high column effect value of 7.6. This is the rate of infant mortality relative to the overall median (i.e. on average, 20.6 infants per thousand die every year and the rate goes up to 7.6 + 20.6 for infants whose father has not passed the 8th grade). Infants whose father has completed more than 16 years of school (i.e. who has completed college) have a lower rate of mortality as indicated by the low effect value of -3.5 (i.e. 3.5 fewer depths than average). The effects from regions also show higher infant mortality rates for North Central and Western regions (with effect values of 2.6 and 0.4 respectively) and lower rates for the northeastern and southern regions; however the regional effect does not appear to be as dominant as that of the father’s educational attainment.\nIt’s also important to look at the distribution of the residual numbers across the two-way table. One should identify unusually high or low residuals in the table. Such residuals may warrant further investigation (e.g. the high southern region residual value of 10.8 may need further exploration).\n\n\n28.3.4 Test for non-additivity\nThus far, we have assumed an additive relationship between the effects (factors). But this additive model may not be the best fit for our data. A good way to test this is by generating a Tukey Additivity Plot where we plot residuals vs. the comparison value, \\(cv_{ij}\\), defined as \\(\\alpha_i \\beta_j / \\mu\\). If the plot is devoid of any obvious trend or pattern we can conclude that our dataset is consistent with an additive model. Such seems to be the case with our working example as shown in the following plot.\n\n\n\n\n\nIf the diagnostic plot revealed a trend, its shape–or more specifically, its slope–could be used in helping define an appropriate transformation for the data. A rule of thumb is to apply a \\((1 – slope)\\) power transformation to the original data. If the resulting power is not appropriate for the dataset then the \\(cv_{ij}\\) can be added to the additive model as follows:\n\\[\ny_{ij} = \\mu + \\alpha_i + \\beta_j + k(\\alpha_i \\beta_j / \\mu)  + \\epsilon_{ij}\n\\] where the constant \\(k\\) is the slope estimated from the Tukey Additivity plot."
+    "text": "28.3 Analysis workflow\nLet’s first create the dataframe.\n\ndf &lt;- data.frame(row.names = c(\"NE\",\"NC\",\"S\",\"W\"),\n                 ed8       = c(25.3,32.1,38.8,25.4), \n                 ed9to11   = c(25.3,29,31,21.1),\n                 ed12      = c(18.2,18.8,19.3,20.3),\n                 ed13to15  = c(18.3,24.3,15.7,24),\n                 ed16      = c(16.3,19,16.8,17.5)\n                 )\n\n\n28.3.1 Visualizing the data\nIt’s often easier to look at a graphical representation of the data than a tabular one. Even a table as small as this can benefit from a plot.\nWe will adopt Bill Cleveland’s dotplot for this purpose. R has a built-in dotplot function called dotchart. It requires that the table be stored as a matrix instead of a dataframe; we will therefore convert df to a matrix by wrapping it with the as.matrix() function.\n\ndotchart( as.matrix(df), cex=0.7)\n\n\n\n\nThe plot helps visualize any differences in mortality rates across different father educational attainment levels. There seems to be a gradual decrease in child mortality with increasing father educational attainment.\nBut the plot does not help spot differences across regions (except for the ed12 group). We can generate another plot where region becomes the main grouping factor. We do this by wrapping the matrix with the transpose function t().\n\ndotchart( t(as.matrix(df)), cex=0.7)\n\n\n\n\nAt first glance, there seems to be higher death rates for the NC and S regions and relatively lower rates for the W and NE regions. But our eyes may be fooled by outliers in the data.\nNext, we’ll generate side-by-side boxplots to compare the effects between both categories. Note that we’ll need to create a long version of the table using the tidyr package.\n\nlibrary(tidyr)\nlibrary(dplyr)\n\ndf.l &lt;- df %&gt;%\n  mutate(Region = as.factor(row.names(.)) )  %&gt;%\n  pivot_longer(names_to = \"Edu\", values_to = \"Deaths\", -Region) %&gt;% \n  mutate(Edu = factor(Edu, levels=names(df)))\n\n# side-by-side plot\nOP &lt;- par(mfrow=c(1,2))\nplot(Deaths ~ Region + Edu, df.l)\npar(OP)\n\n\n\n\nSo at this point the plots suggest that there may be a father’s educational attainment effect as well as a regional effect on infant mortality with the former effect being possibly more important.\nNext we will attempt to quantify these effects using Tukey’s median polish.\n\n\n28.3.2 Median polish\nA median polish is an iterative process whereby we seek to extract the effects from the data into their row and column components. The steps involved in each iteration are described next.\n\n28.3.2.1 Step 1: Compute overall median and residual table\nFirst, we will compute the median value for all values in the dataset; this will be our first estimate of the common value. The resulting median value is placed in the upper left-hand margin of the table. A residual table is created by taking the difference between the original value and the overall median.\n\nYou’ll also note that the row and column effect values are initially populated with 0 values.\n\n\n28.3.2.2 Step 2: Compute the row medians\nNext, row medians are computed (shown in red in the following figure) from the residual table.\n\nNote that in addition to computing the row median values for the response variable, we are also computing the median value for the column effects in the upper margin (which happens to be zero since we have no column effect values yet).\n\n\n28.3.2.3 Step 3: Create a new residual table from the row medians\nThe row medians are added to the left hand margin of the new residual table (shown in green). These values represent the row effects in this first iteration. A new set of residual values is created from the row medians where each cell takes on the value of the subtraction of the row median from each response variable in that row.\n\nThe subtraction is also done for the column effect values (even though all values remain zero) and the global median (common effect value).\n\n\n28.3.2.4 Step 4: Compute the column medians\nNext, column medians are computed (shown in red in the following figure) from the residual table. Note that you are NOT computing the column medians from the original response variables. Note too that we are also computing the median value for the row effect values (left margin).\n\n\n\n28.3.2.5 Step 5: Create a new residual table from the column medians\nA new residual table is created from the column medians where each new cell takes on the value of the subtraction of the column median from each initial residual value in that column. For example, the first upper-left cell’s residual is \\(7.0 - 7.4 = -0.4\\).\nThe column medians are also added to the column effect margin. Note that we are also adding the row effect median value, -0.5, to the common effect cell in the upper left-hand corner.\n\nWe have now completed our first row and column smoothing iteration. However, there may be more row and column effects that can be pulled from the residuals. We therefore move on to a second iteration.\n\n\n28.3.2.6 Step 6: Second iteration – row effects\nNext, we compute the row medians from the residuals, then add the column effect median to the top margin (and the common value) and subtract the row medians from the residuals. as in step 3.\n\n\n\n28.3.2.7 Step 7: Second iteration – column effects\nTo wrap-up the second iteration, we compute the column median values from the residuals then subtract the medians from those residuals. We also add the medians to the row effect values and the common value.\n\n\n\n28.3.2.8 When do we stop iterating?\nThe goal is to iterate through the row and column smoothing operations until the row and column effect medians are close to 0. However, Hoaglin et al. (1983) warn against “using unrestrained iteration” and suggest that a few steps should be more than adequate in most instances. In our working example, a third iteration may be warranted. The complete workflow for this third iteration follows.\n\nThe final version of the table (along with the column and row values) is shown below:\n\nA progression of each iteration is shown side-by-side in the next figure. The left-most table is the original data showing death rates. The second table shows the outcome of the first round of polishing (including the initial overall median value of 20.2). The third and forth table show the second and third iterations of the smoothing operations. Additional iterations are not deemed necessary given that little more can be extracted from the residuals.\n\n\n\n\n\nA few takeway observations from the analysis follow:\n\nThe common factor ends up contributing the largest effect.\nBetween Region and Education effects, the latter explains more of the variability in deaths with a range of 7.43 - (-3.7) = 11.13 deaths per 1000 births. The Region effect explains about 2.3 - (-1.46) = 3.76 deaths per 1000 births of the variability in the data.\nOne abnormally high residual–10.97 deaths above the overall median deaths–is associated with the southern region and an 8th grade education.\n\n\n\n\n28.3.3 Interpreting the median polish\nAs noted earlier, a two-way table represents the relationship between the response variable, \\(y\\), and the two categories as:\n\\[\ny_{ij} = \\mu + \\alpha_i + \\beta_j + \\epsilon_{ij}\n\\]\nIn our working example, \\(\\mu\\) = 20.6; \\(\\alpha_i\\) = -1.5, 2.6 -0.4, and 0.4 for \\(i\\) = NE, Nc, S and W respectively; \\(\\beta_j\\) = 7.6, 6.0, -0.9, 0.2 and -3.5 for \\(j\\) = ed8, ed9to11, ed12, ed13to15 and ed16 respectively.\nThe residuals represent the portion of the mortality rates that can’t be explained by either factors.\nSo the mortality rate in the upper left-hand cell from the original table can be deconstructed as:\n\nThe examination of the table suggests that the infant mortality rate is greatest for fathers who did not attain more than 8 years of school (i.e. who has not completed high school) as noted by the high column effect value of 7.6. This is the rate of infant mortality relative to the overall median (i.e. on average, 20.6 infants per thousand die every year and the rate goes up to 7.6 + 20.6 for infants whose father has not passed the 8th grade). Infants whose father has completed more than 16 years of school (i.e. who has completed college) have a lower rate of mortality as indicated by the low effect value of -3.5 (i.e. 3.5 fewer depths than average). The effects from regions also show higher infant mortality rates for North Central and Western regions (with effect values of 2.6 and 0.4 respectively) and lower rates for the northeastern and southern regions; however the regional effect does not appear to be as dominant as that of the father’s educational attainment.\nIt’s also important to look at the distribution of the residual numbers across the two-way table. One should identify unusually high or low residuals in the table. Such residuals may warrant further investigation (e.g. the high southern region residual value of 10.8 may need further exploration).\n\n\n28.3.4 Test for non-additivity\nThus far, we have assumed an additive relationship between the effects (factors). But this additive model may not be the best fit for our data. A good way to test this is by generating a Tukey Additivity Plot where we plot residuals vs. the comparison value, \\(cv_{ij}\\), defined as \\(\\alpha_i \\beta_j / \\mu\\). If the plot is devoid of any obvious trend or pattern we can conclude that our dataset is consistent with an additive model. Such seems to be the case with our working example as shown in the following plot.\n\n\n\n\n\nIf the diagnostic plot revealed a trend, its shape–or more specifically, its slope–could be used in helping define an appropriate transformation for the data. A rule of thumb is to apply a \\((1 – slope)\\) power transformation to the original data. If the resulting power is not appropriate for the dataset then the \\(cv_{ij}\\) can be added to the additive model as follows:\n\\[\ny_{ij} = \\mu + \\alpha_i + \\beta_j + k(\\alpha_i \\beta_j / \\mu)  + \\epsilon_{ij}\n\\] where the constant \\(k\\) is the slope estimated from the Tukey Additivity plot."
   },
   {
     "objectID": "two_way.html#implementing-the-median-polish-in-r",
     "href": "two_way.html#implementing-the-median-polish-in-r",
     "title": "28  Analyzing two-way tables: The median polish",
     "section": "28.4 Implementing the median polish in R",
-    "text": "28.4 Implementing the median polish in R\n\n28.4.1 Using base R’s medpolish\nThe steps outlined in the previous section can be easily implemented using pen and paper or a spreadsheet environment for larger datasets. R has a built-in function called medpolish() which does this for us. We can define the maximum number of iterations by setting the maxiter= parameter but note that medpolish will, by default, automatically estimate the best number of iterations for us.\n\ndf.med &lt;- medpolish( df)\n\n1: 47.1\n2: 42.9\n3: 42.45\nFinal: 42.25\n\n\nThe median polish ran through 4 iterations–one more than was run in our step-by-step example.\nThe four values printed to the console gives us the sum of absolute residuals at each iteration. The output from the model is stored in the df.med object. To see its contents , simply type it at the command line.\n\ndf.med\n\n\nMedian Polish Results (Dataset: \"df\")\n\nOverall: 20.6\n\nRow Effects:\n   NE    NC     S     W \n-1.50  2.55 -0.35  0.35 \n\nColumn Effects:\n     ed8  ed9to11     ed12 ed13to15     ed16 \n   7.575    6.025   -0.925    0.175   -3.450 \n\nResiduals:\n      ed8 ed9to11   ed12 ed13to15  ed16\nNE -1.375   0.175  0.025   -0.975  0.65\nNC  1.375  -0.175 -3.425    0.975 -0.70\nS  10.975   4.725 -0.025   -4.725  0.00\nW  -3.125  -5.875  0.275    2.875  0.00\n\n\nAll three effects are displayed as well as the residuals (note that the precision returned is greater than that used in our earlier analysis).\nTo generate the Tukey additivity plot, simply wrap the median polish object with the plot command as in:\n\nplot( df.med )\n\n\n\n\n\n\n28.4.2 Using eda_pol\nThe tukeyedar package has a custom function called eda_pol that will output the polished table as a graphic element.\nUnlike the base medpolish function, eda_pol requires that the data be in long form. In creating the long form version of the data, we’ll explicitly set the levels for educational attainment and region.\n\ngrd &lt;- c(\"ed8\", \"ed9-11\", \"ed12\", \"ed13-15\", \"ed16\")\nreg &lt;- c(\"NE\", \"NC\", \"S\", \"W\")\ndfl &lt;- data.frame(Region =  factor(rep( c(\"NE\", \"NC\", \"S\", \"W\"), each = 5),\n                                   levels = reg),\n                 Education = factor(rep( grd , 4), levels = grd),\n                 Deaths = c(25.3, 25.3, 18.2, 18.3, 16.3, 32.1, 29, 18.8,\n                          24.3, 19, 38.8, 31, 19.3, 15.7, 16.8, 25.4, \n                          21.1, 20.3, 24, 17.5))\n\nThe eda_pol requires that you specify which variables will be mapped to the row and column effects.\n\nlibrary(tukeyedar)\n\ndf.pol &lt;- eda_pol(dfl, row = Region, col = Education , val = Deaths )\n\n\n\n\n\n\nHere, the eda_pol function required five iterations which explains the slight difference in values from those generated with the medpolish function. You can control the maximum number of iterations in eda_pol with the maxiter argument.\nBy default, the column and row values are listed in the same order as they are presented in the input data frame (either defined alphabetically or defined by their levels). You can have the function sort the values by their effect size using the sort argument.\n\ndf.pol &lt;- eda_pol(dfl, row = Region, col = Education, val = Deaths, sort = TRUE )\n\n\n\n\n\n\nThe df.pol object is a list that stores output values such as the common, row and column effects.\n\ndf.pol$global\n\n[1] 20.85\n\ndf.pol$row\n\n  Region  effect\n1     NE -1.4625\n2     NC  2.3000\n3      S -0.3500\n4      W  0.3500\n\ndf.pol$col\n\n  Education   effect\n1       ed8  7.43125\n2    ed9-11  5.88125\n3      ed12 -1.19375\n4   ed13-15  0.03125\n5      ed16 -3.70000\n\n\nThe eda_pol will also output a diagnostic plot.\n\nplot(df.pol, type = \"diagnostic\")\n\n$slope\n    cv \n1.3688 \n\n\n\n\n\nIn addition to generating the plot, the function will also return the slope of the fitted line.\nFor a comprehensive review of the median polish and the eda_pol function, see the eda_pol vignette."
+    "text": "28.4 Implementing the median polish in R\n\n28.4.1 Using base R’s medpolish\nThe steps outlined in the previous section can be easily implemented using pen and paper or a spreadsheet environment for larger datasets. R has a built-in function called medpolish() which does this for us. We can define the maximum number of iterations by setting the maxiter= parameter but note that medpolish will, by default, automatically estimate the best number of iterations for us.\n\ndf.med &lt;- medpolish( df)\n\n1: 47.1\n2: 42.9\n3: 42.45\nFinal: 42.25\n\n\nThe median polish ran through 4 iterations–one more than was run in our step-by-step example.\nThe four values printed to the console gives us the sum of absolute residuals at each iteration. The output from the model is stored in the df.med object. To see its contents , simply type it at the command line.\n\ndf.med\n\n\nMedian Polish Results (Dataset: \"df\")\n\nOverall: 20.775\n\nRow Effects:\n    NE     NC      S      W \n-1.475  2.375 -0.350  0.350 \n\nColumn Effects:\n     ed8  ed9to11     ed12 ed13to15     ed16 \n  7.4750   5.9250  -1.1125   0.0750  -3.6250 \n\nResiduals:\n      ed8 ed9to11    ed12 ed13to15   ed16\nNE -1.475   0.075  0.0125   -1.075  0.625\nNC  1.475  -0.075 -3.2375    1.075 -0.525\nS  10.900   4.650 -0.0125   -4.800  0.000\nW  -3.200  -5.950  0.2875    2.800  0.000\n\n\nAll three effects are displayed as well as the residuals (note that the precision returned is greater than that used in our earlier analysis).\nTo generate the Tukey additivity plot, simply wrap the median polish object with the plot command as in:\n\nplot( df.med )\n\n\n\n\n\n\n28.4.2 Using eda_pol\nThe tukeyedar package has a custom function called eda_pol that will output the polished table as a graphic element.\nUnlike the base medpolish function, eda_pol requires that the data be in long form. In creating the long form version of the data, we’ll explicitly set the levels for educational attainment and region.\n\ngrd &lt;- c(\"ed8\", \"ed9-11\", \"ed12\", \"ed13-15\", \"ed16\")\nreg &lt;- c(\"NE\", \"NC\", \"S\", \"W\")\ndfl &lt;- data.frame(Region =  factor(rep( c(\"NE\", \"NC\", \"S\", \"W\"), each = 5),\n                                   levels = reg),\n                 Education = factor(rep( grd , 4), levels = grd),\n                 Deaths = c(25.3, 25.3, 18.2, 18.3, 16.3, 32.1, 29, 18.8,\n                          24.3, 19, 38.8, 31, 19.3, 15.7, 16.8, 25.4, \n                          21.1, 20.3, 24, 17.5))\n\nThe eda_pol requires that you specify which variables will be mapped to the row and column effects.\n\nlibrary(tukeyedar)\n\ndf.pol &lt;- eda_pol(dfl, row = Region, col = Education , val = Deaths )\n\n\n\n\n\n\nHere, the eda_pol function required five iterations which explains the slight difference in values from those generated with the medpolish function. You can control the maximum number of iterations in eda_pol with the maxiter argument.\nBy default, the column and row values are listed in the same order as they are presented in the input data frame (either defined alphabetically or defined by their levels). You can have the function sort the values by their effect size using the sort argument.\n\ndf.pol &lt;- eda_pol(dfl, row = Region, col = Education, val = Deaths, sort = TRUE )\n\n\n\n\n\n\nThe df.pol object is a list that stores output values such as the common, row and column effects.\n\ndf.pol$global\n\n[1] 20.85\n\ndf.pol$row\n\n  Region  effect\n1     NE -1.4625\n2     NC  2.3000\n3      S -0.3500\n4      W  0.3500\n\ndf.pol$col\n\n  Education   effect\n1       ed8  7.43125\n2    ed9-11  5.88125\n3      ed12 -1.19375\n4   ed13-15  0.03125\n5      ed16 -3.70000\n\n\nThe eda_pol will also output a diagnostic plot.\n\nplot(df.pol, type = \"diagnostic\")\n\n$slope\n    cv \n1.3688 \n\n\n\n\n\nIn addition to generating the plot, the function will also return the slope of the fitted line in the console.\nFor a comprehensive review of the median polish and the eda_pol function, see the eda_pol vignette."
   },
   {
     "objectID": "two_way.html#what-if-we-use-the-mean-instead-of-the-median",
     "href": "two_way.html#what-if-we-use-the-mean-instead-of-the-median",
     "title": "28  Analyzing two-way tables: The median polish",
     "section": "28.5 What if we use the mean instead of the median?",
-    "text": "28.5 What if we use the mean instead of the median?\nThe procedure is similar with some notable differences. First, we compute the global mean instead of the medians (the common effect) then subtract it from all values in the table. Next, we compute the row means (the row effect) then subtract each mean from all values in its associated row. We finally compute the column means (from the residuals) to give us the column effect. That’s it, unlike the median polish, we do not iterate the smoothing sequences. An example of a “mean” polish applied to our data follows:\n\nThe results differ slightly from those produced using the median polish. Recall that the mean is sensitive to outliers and medians are not. If a robust estimate of the effects is sought, stick with the median polish.\nSo what can we gain from the “mean” polish? Well, as it turns out, it serves as the basis of the two-way ANOVA.\n\n28.5.1 Two-way ANOVA\nA two-way ANOVA assesses whether the factors (categories) have a significant effect on the outcome. Its implementation can be conceptualized as a regression analysis where each row and column level is treated as a dummy variable. The computed regression coefficients are the same as the levels computed using the “mean” polish outlined in the previous step except for one notable difference: the ANOVA adds one of the column level values and one of the row level values to the grand mean. It then subtracts those values from their respective row/column effects.\n\nThe resulting row/column levels are the regression coefficients:\n\nM &lt;- lm(Deaths ~ Education + Region, dfl)\ncoef(M)\n\n     (Intercept)  Educationed9-11    Educationed12 Educationed13-15    Educationed16         RegionNE          RegionS \n          32.215           -3.800          -11.250           -9.825          -13.000           -3.960           -0.320 \n         RegionW \n          -2.980 \n\n\nThe regression formula takes on the form:\n\\[\n\\begin{split}\nDeaths = 32.21 - 3.8(ed9to11) -11.25(ed12) -9.83(ed13to15) -13(ed16) \\\\\n         -3.96(RegionNE) -0.320(RegionS)  -2.980(RegionW) + \\epsilon_{ij}\n\\end{split}\n\\]\nSo, for example, to compute the first value in the raw table (death rate = 25.3) from the formula, substitute the variables as follows:\n\\[\n\\begin{split}\nDeaths &= 32.21 - 3.8(0) -11.25(0) -9.83(0) -13(0)&& \\\\\n         &-3.96(1) -0.320(0)  -2.980(0) - 3.0& \\\\\n       &= 32.21 -3.96 -3.0&& \\\\\n       &= 25.3&&\n\\end{split}\n\\]\nTo assess if any of the factors have a significant effect on the Death variables, simply wrap the regression with an anova() function.\n\nanova(M)\n\nAnalysis of Variance Table\n\nResponse: Deaths\n          Df Sum Sq Mean Sq F value  Pr(&gt;F)   \nEdu        4 478.52 119.630  7.7539 0.00251 **\nRegion     3  57.44  19.146  1.2410 0.33801   \nResiduals 12 185.14  15.428                   \n---\nSignif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1\n\n\nThe results suggest that the father’s educational attainment has a significant effect on infant mortality (p = 0.0025) whereas the region does not (p = 0.338).\n\n\n28.5.2 Visualizing a two-way ANOVA as a table\nYou can view the output of a two-way ANOVA as a “mean polish” table using the eda_pol function. Simply set the stat argument to mean.\n\ndf.ano &lt;- eda_pol(dfl, row = Region, col = Education, val = Deaths, stat = mean)"
+    "text": "28.5 What if we use the mean instead of the median?\nThe procedure is similar with some notable differences. First, we compute the global mean instead of the medians (the common effect) then subtract it from all values in the table. Next, we compute the row means (the row effect) then subtract each mean from all values in its associated row. We finally compute the column means (from the residuals) to give us the column effect. That’s it, unlike the median polish, we do not iterate the smoothing sequences. An example of a “mean” polish applied to our data follows:\n\nThe results differ slightly from those produced using the median polish. Recall that the mean is sensitive to outliers and medians are not. If a robust estimate of the effects is sought, stick with the median polish.\nSo what can we gain from the “mean” polish? Well, as it turns out, it serves as the basis of the two-way ANOVA.\n\n28.5.1 Two-way ANOVA\nA two-way ANOVA assesses whether the factors (categories) have a significant effect on the outcome. Its implementation can be conceptualized as a regression analysis where each row and column level is treated as a dummy variable. The computed regression coefficients are the same as the levels computed using the “mean” polish outlined in the previous step except for one notable difference: the ANOVA adds one of the column level values and one of the row level values to the grand mean. It then subtracts those values from their respective row/column effects.\n\nThe resulting row/column levels are the regression coefficients:\n\nM &lt;- lm(Deaths ~ Education + Region, dfl)\ncoef(M)\n\n     (Intercept)  Educationed9-11    Educationed12 Educationed13-15    Educationed16         RegionNC          RegionS \n          28.255           -3.800          -11.250           -9.825          -13.000            3.960            3.640 \n         RegionW \n           0.980 \n\n\nThe regression formula takes on the form:\n\\[\n\\begin{split}\nDeaths = 32.21 - 3.8(ed9to11) -11.25(ed12) -9.83(ed13to15) -13(ed16) \\\\\n         -3.96(RegionNE) -0.320(RegionS)  -2.980(RegionW) + \\epsilon_{ij}\n\\end{split}\n\\]\nSo, for example, to compute the first value in the raw table (death rate = 25.3) from the formula, substitute the variables as follows:\n\\[\n\\begin{split}\nDeaths &= 32.21 - 3.8(0) -11.25(0) -9.83(0) -13(0)&& \\\\\n         &-3.96(1) -0.320(0)  -2.980(0) - 3.0& \\\\\n       &= 32.21 -3.96 -3.0&& \\\\\n       &= 25.3&&\n\\end{split}\n\\]\nTo assess if any of the factors have a significant effect on the Death variables, simply wrap the regression with an anova() function.\n\nanova(M)\n\nAnalysis of Variance Table\n\nResponse: Deaths\n          Df Sum Sq Mean Sq F value  Pr(&gt;F)   \nEducation  4 478.52 119.630  7.7539 0.00251 **\nRegion     3  57.44  19.146  1.2410 0.33801   \nResiduals 12 185.14  15.428                   \n---\nSignif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1\n\n\nThe results suggest that the father’s educational attainment has a significant effect on infant mortality (p = 0.0025) whereas the region does not (p = 0.338).\n\n\n28.5.2 Visualizing a two-way ANOVA as a table\nYou can view the output of a two-way ANOVA as a “mean polish” table using the eda_pol function. Simply set the stat argument to mean.\n\ndf.ano &lt;- eda_pol(dfl, row = Region, col = Education, val = Deaths, stat = mean)"
   },
   {
     "objectID": "two_way.html#references",
@@ -1181,5 +1181,12 @@
     "title": "27  Slicing data",
     "section": "27.2 Case study: Trends in daily temperature extremes",
     "text": "27.2 Case study: Trends in daily temperature extremes\n\nDisclaimer: the analysis presented here is only exploratory and does not mirror the complete analysis conducted by Vincent et al. nor the one conducted by Stone.\n\n\n27.2.1 Original analysis\nThe following data are pulled from the paper titled “Observed Trends in Indices of Daily Temperature Extremes in South America 1960-2000” (Vincent et al., 2005) and represent the percentage of nights with temperatures greater than or colder than the 90th and 10th percentiles respectively within each year. The percentiles were calculated for the 1961 to 2000 period.\n\nlibrary(tidyr)\n\nYear &lt;- 1960:2000\nPerC &lt;- c(11.69,9.33,14.35,10.73,14.15,11.16,13,12.13,14.25,10.01,11.94,14.35,\n          10.83,9.38,11.5,10.44,12.66,7.55,9.77,9.81,8.9,8.51,7.02,6.83,9.67,\n          7.84,7.11,8.56,10.59,7.93,8.85,8.8,8.75,8.18,7.16,9.91,10.15,6.58,\n          6.44,9.43,8.03)\nPerH &lt;- c(8.62,10.1,6.67,11.13,5.71,9.48,7.63,8.12,7.2,9.64,8.42,5.71,11.72,\n          11.32,7.2,7.17,7.46,13.17,9.28,8.75,12.38,10,13.83,17.59,10.14,\n          9.84,11.23,14.39,9.44,8.26,12.15,12.45,13.14,13.67,15.22,11.79,11.16,\n          20.37,17.56,11.13,11.49)\ndf2 &lt;- data.frame(Year, PerC, PerH)\n\ndf2.l &lt;- pivot_longer(df2, names_to = \"Temp\", values_to = \"Percent\", -Year)\n\nLet’s plot the data and fit a straight line to the points.\n\nggplot(df2.l, aes(x = Year, y = Percent)) + geom_point() +\n  stat_smooth(method = \"lm\", se = FALSE) + facet_wrap(~ Temp, nrow = 1)\n\n\n\n\nThe plot on the left shows percent cold nights and the one on the right shows percent hot nights. At first glance, the trends seem real and monotonic.\nNext we’ll fit a loess to see if the trends are indeed monotonic. To minimize the undue influence of end values in the plot, we’ll implement loess’ bisquare estimation method via the family=symmetric option. We’ll also use a small span to help identify any “kinks” in the patterns.\n\nggplot(df2.l, aes(x = Year, y = Percent)) + geom_point() +\n  stat_smooth(method = \"loess\", se = FALSE, span = 0.5, \n              method.args = list(degree = 1, family = \"symmetric\"))  + \n  facet_wrap(~ Temp, nrow = 1)\n\n\n\n\nThe patterns seem to be segmented around the 1975-1980 period for both plots suggesting that the observed trends may not be monotonic. In fact, there appears to be a prominent kink in the percent cold data around the mid to late 1970`s. A similar, but not as prominent kink can also be observed in the percent hot data at around the same time period.\n\n\n27.2.2 Changepoint\nIn a comment to Vincent et al.’s paper, R.J. Stone argues that the trend observed in the percent hot and cold dataset is not monotonic but segmented instead. In other words, there is an abrupt change in patterns for both datasets that make it seem as though a monotonic trend exists when in fact the data may follow relatively flat patterns for two different segments of time. He notes that the abrupt change (which he refers to as a changepoint) occurs around the 1976 and 1977 period. He suggests that this time period coincides with a change in the Pacific Decadal Oscillation (PDO) pattern. PDO refers to an ocean/atmospheric weather pattern that spans multiple decades and that is believed to impact global climate.\nThe following chunk of code loads the PDO data, then summarizes the data by year before plotting the resulting dataset.\n\ndf3  &lt;- read.table(\"http://mgimond.github.io/ES218/Data/PDO.dat\", \n                   header = TRUE, na.strings = \"-9999\")\npdo &lt;- df3 %&gt;%\n  pivot_longer(names_to = \"Month\", values_to = \"PDO\", -YEAR) %&gt;%\n  group_by(YEAR) %&gt;%\n  summarise(PDO = median(PDO) )\n\nggplot(pdo, aes(x = YEAR, y = PDO)) + geom_line() + \n  stat_smooth(se = FALSE, span = 0.25) +\n  geom_vline(xintercept = c(1960, 1976, 2000), lty = 3)\n\n\n\n\nThe contrast in PDO indexes between the 1960-1976 period and the 1976-2000 period is obvious with the pre-1977 index values appearing to remain relatively flat over a 15 year period and with the post-1977 index appearing to show a gradual increase towards a peak around the early 1990’s.\nTo see if distinct patterns emerge from the percent hot and cold data before and after 1976, we’ll split the data into two segments using a cutoff year of 1976-1977. Values associated with a period prior to 1977 will be assigned a seg value of Before and those associated with a post-1977 period will be assigned a seg value of After.\n\ndf2.l$seg &lt;- ifelse(df2.l$Year &lt; 1977, \"Before\", \"After\")\n\nNext, we’ll plot the data across four facets:\n\nggplot(df2.l, aes(x = Year, y = Percent)) + geom_point() +\n  stat_smooth(method = \"lm\", se = FALSE) + facet_grid(seg ~ Temp)\n\n\n\n\nWe can also choose to map seg to the color aesthetics which will split the points by color with the added benefit of fitting two separate models to each batch.\n\nggplot(df2.l, aes(x = Year, y = Percent, col = seg)) + geom_point() +\n  stat_smooth(method = \"lm\", se=FALSE) + facet_wrap( ~ Temp, nrow = 1)\n\n\n\n\nTo check for “straightness” in the fits, we’ll fit a loess to the points.\n\nggplot(df2.l, aes(x = Year, y = Percent, col = seg)) + geom_point() +\n  stat_smooth(method = \"loess\", se = FALSE, method.args = list(degree = 1)) + \n  facet_wrap( ~ Temp, nrow = 1)\n\n\n\n\nThere is a clear “stair-step” pattern for the percent cold nights. However, there seems to be an upward trend in the percent of hot nights for the post-1977 period which could imply that in addition to the PDO effect, another process could be at play.\nThis little exercise highlights the ease in which an analysis can follow different (and seemingly sound) paths. It also serves as a reminder of the importance of domain knowledge when exploring data."
+  },
+  {
+    "objectID": "sl_plot.html#example-the-iris-dataset",
+    "href": "sl_plot.html#example-the-iris-dataset",
+    "title": "20  Spread-level plots",
+    "section": "20.3 Example: the iris dataset",
+    "text": "20.3 Example: the iris dataset\nR has a built-in dataset called iris that provide measurements of sepal and petal dimensions for three different species of the iris family. In this next example, we will plot the spreads of the Petal.Length residuals (after removing their group median values) to their group medians.\n\n# Create two new columns: group median and group residuals\n df1 &lt;- iris %&gt;%\n   group_by(Species)  %&gt;%\n   mutate( Median = median(Petal.Length),\n           Residuals = sqrt(abs( Petal.Length - Median)))  \n\n# Generate the s-l plot \n ggplot(df1, aes(x = Median, y = Residuals)) + \n   geom_jitter(alpha = 0.4, width = 0.05, height = 0) +\n   stat_summary(fun = median, geom = \"line\", col = \"red\") +\n   ylab(expression(sqrt( abs(\" Residuals \")))) +\n   geom_text(aes(x = Median, y = 1.3, label = Species))\n\n\n\n\nA monotonic spread is apparent in this dataset too, i.e. as the median length of the Petal increases, so does the spread."
   }
 ]
\ No newline at end of file
diff --git a/docs/sl_plot.html b/docs/sl_plot.html
index 34f694ef..63f20289 100755
--- a/docs/sl_plot.html
+++ b/docs/sl_plot.html
@@ -336,7 +336,7 @@
 </li>
           <li class="sidebar-item">
   <div class="sidebar-item-container"> 
-  <a href="./robustness.html" class="sidebar-item-text sidebar-link">
+  <a href="./case_study.html" class="sidebar-item-text sidebar-link">
  <span class="menu-text"><span class="chapter-number">23</span>&nbsp; <span class="chapter-title">A working example: t-tests and re-expression</span></span></a>
   </div>
 </li>
@@ -435,7 +435,7 @@ <h2 id="toc-title">Table of contents</h2>
   <ul>
   <li><a href="#introduction" id="toc-introduction" class="nav-link active" data-scroll-target="#introduction"><span class="header-section-number">20.1</span> Introduction</a></li>
   <li><a href="#constructing-the-s-l-plot" id="toc-constructing-the-s-l-plot" class="nav-link" data-scroll-target="#constructing-the-s-l-plot"><span class="header-section-number">20.2</span> Constructing the s-l plot</a></li>
-  <li><a href="#rs-built-in-iris-dataset" id="toc-rs-built-in-iris-dataset" class="nav-link" data-scroll-target="#rs-built-in-iris-dataset"><span class="header-section-number">20.3</span> R’s built in <code>iris</code> dataset</a></li>
+  <li><a href="#example-the-iris-dataset" id="toc-example-the-iris-dataset" class="nav-link" data-scroll-target="#example-the-iris-dataset"><span class="header-section-number">20.3</span> Example: the <code>iris</code> dataset</a></li>
   <li><a href="#how-can-we-stabilize-spreads-in-a-dataset" id="toc-how-can-we-stabilize-spreads-in-a-dataset" class="nav-link" data-scroll-target="#how-can-we-stabilize-spreads-in-a-dataset"><span class="header-section-number">20.4</span> How can we stabilize spreads in a dataset?</a>
   <ul class="collapse">
   <li><a href="#variations-of-the-s-l-plot" id="toc-variations-of-the-s-l-plot" class="nav-link" data-scroll-target="#variations-of-the-s-l-plot"><span class="header-section-number">20.4.1</span> Variations of the S-L plot</a></li>
@@ -463,7 +463,7 @@ <h1 class="title"><span class="chapter-number">20</span>&nbsp; <span class="chap
 
 </header>
 
-<div class="cell" data-hash="sl_plot_cache/html/unnamed-chunk-2_ddef71e58918eb94c9077c76bf467040">
+<div class="cell" data-hash="sl_plot_cache/html/unnamed-chunk-2_a9e0c5c25fa201d342586fad77e0b216">
 <div class="cell-output-display">
 
 <table data-quarto-disable-processing="true" class="table" style="width: auto !important; ">
@@ -476,6 +476,13 @@ <h1 class="title"><span class="chapter-number">20</span>&nbsp; <span class="chap
              padding-right:4px !important;
              font-size: 0.8em;
              border-radius: 5px;"> dplyr </th>
+   <th style="text-align:left;color: rgba(85, 85, 85, 255) !important;background-color: rgba(221, 221, 221, 255) !important;text-align: center;border: 1px solid white !important;
+             font-family: 'Source Code Pro', 'Open Sans';
+             padding:1px !important;
+             padding-left:4px !important;
+             padding-right:4px !important;
+             font-size: 0.8em;
+             border-radius: 5px;"> ggplot2 </th>
   </tr>
  </thead>
 <tbody>
@@ -487,6 +494,13 @@ <h1 class="title"><span class="chapter-number">20</span>&nbsp; <span class="chap
              padding-right:4px !important;
              font-size: 0.8em;
              border-radius: 5px;"> 1.1.4 </td>
+   <td style="text-align:left;color: darkred !important;background-color: rgba(250, 232, 232, 255) !important;text-align: center;border: 1px solid white;
+             font-family: 'Open Sans', Arial;
+             padding:1px !important;
+             padding-left:4px !important;
+             padding-right:4px !important;
+             font-size: 0.8em;
+             border-radius: 5px;"> 3.4.4 </td>
   </tr>
 </tbody>
 </table>
@@ -525,8 +539,8 @@ <h2 data-number="20.2" class="anchored" data-anchor-id="constructing-the-s-l-plo
 <p>The singer dataset does not seem to exhibit any dependence between a voice part’s spread and its median value.</p>
 <p>Next, we’ll look at an example of a dataset that <em>does</em> exhibit a dependence between spread and fitted values.</p>
 </section>
-<section id="rs-built-in-iris-dataset" class="level2" data-number="20.3">
-<h2 data-number="20.3" class="anchored" data-anchor-id="rs-built-in-iris-dataset"><span class="header-section-number">20.3</span> R’s built in <code>iris</code> dataset</h2>
+<section id="example-the-iris-dataset" class="level2" data-number="20.3">
+<h2 data-number="20.3" class="anchored" data-anchor-id="example-the-iris-dataset"><span class="header-section-number">20.3</span> Example: the <code>iris</code> dataset</h2>
 <p>R has a built-in dataset called <code>iris</code> that provide measurements of sepal and petal dimensions for three different species of the iris family. In this next example, we will plot the spreads of the <code>Petal.Length</code> residuals (after removing their group median values) to their group medians.</p>
 <div class="cell" data-hash="sl_plot_cache/html/unnamed-chunk-4_c31935faa3125f3844ba51b516baafae">
 <div class="sourceCode cell-code" id="cb2"><pre class="sourceCode r code-with-copy"><code class="sourceCode r"><span id="cb2-1"><a href="#cb2-1" aria-hidden="true" tabindex="-1"></a><span class="co"># Create two new columns: group median and group residuals</span></span>
diff --git a/docs/sl_plot_files/figure-html/unnamed-chunk-3-1.png b/docs/sl_plot_files/figure-html/unnamed-chunk-3-1.png
index e9292ca6..d1ef501b 100755
Binary files a/docs/sl_plot_files/figure-html/unnamed-chunk-3-1.png and b/docs/sl_plot_files/figure-html/unnamed-chunk-3-1.png differ
diff --git a/docs/sl_plot_files/figure-html/unnamed-chunk-4-1.png b/docs/sl_plot_files/figure-html/unnamed-chunk-4-1.png
index 8639b65a..b54cbc3e 100755
Binary files a/docs/sl_plot_files/figure-html/unnamed-chunk-4-1.png and b/docs/sl_plot_files/figure-html/unnamed-chunk-4-1.png differ
diff --git a/docs/sl_plot_files/figure-html/unnamed-chunk-7-1.png b/docs/sl_plot_files/figure-html/unnamed-chunk-7-1.png
index 626d5267..74d9fea7 100644
Binary files a/docs/sl_plot_files/figure-html/unnamed-chunk-7-1.png and b/docs/sl_plot_files/figure-html/unnamed-chunk-7-1.png differ
diff --git a/docs/two_way.html b/docs/two_way.html
index dbfeaccb..1ec8dfdf 100644
--- a/docs/two_way.html
+++ b/docs/two_way.html
@@ -478,15 +478,6 @@ <h1 class="title"><span class="chapter-number">28</span>&nbsp; <span class="chap
 
 </header>
 
-<div class="cell">
-<div class="cell-output cell-output-stderr">
-<pre><code>Clean these cache files?
-
-two_way_cache/html/clock_a80bfd981f796072ed50d2d6f8d5d27e.RData
-two_way_cache/html/clock_a80bfd981f796072ed50d2d6f8d5d27e.rdb
-two_way_cache/html/clock_a80bfd981f796072ed50d2d6f8d5d27e.rdx</code></pre>
-</div>
-</div>
 <div class="cell" data-hash="two_way_cache/html/unnamed-chunk-2_d6b817bd2dc7f751b645b175972dbeb2">
 <div class="cell-output-display">
 
@@ -551,7 +542,7 @@ <h2 data-number="28.1" class="anchored" data-anchor-id="introduction"><span clas
 <div class="cell" data-hash="two_way_cache/html/unnamed-chunk-3_bb5971ce07d8d8c0b8943589b38c4d73">
 
 <!-- html table generated in R 4.3.2 by xtable 1.8-4 package -->
-<!-- Sun Feb  4 16:25:02 2024 -->
+<!-- Mon Feb  5 09:05:50 2024 -->
 
 <table border="1">
 <tbody><tr>
@@ -746,12 +737,12 @@ <h2 data-number="28.1" class="anchored" data-anchor-id="introduction"><span clas
 <p><img src="two_way_files/figure-html/unnamed-chunk-4-1.png" class="img-fluid" width="576"></p>
 </div>
 </div>
-<p>We have separated the income values into two parts: an <strong>overall</strong> value and the difference between each value and the overall–the <strong>residuals</strong>. This can be written in the following algebraic form:</p>
+<p>We have separated the income values into two parts: an <strong>overall</strong> value and the difference between each value and the overall–i.e.&nbsp;the <strong>residuals</strong>. This can be written in the following algebraic form:</p>
 <p><span class="math display">\[
 Income_{Occupation} = Income_{median} + Income_{residual}
 \]</span></p>
-<p>For example, ‘Healthcare Support’ average income of $43980 can be separated into the the median, $45265, and the residual, $-1285.</p>
-<p>So in addition to gleaning a median income value from this batch of values, we have also devised a way to facilitate the comparison of the income values to one another. Alternatively, we could have assigned colors to each cell of the table reflecting their values <em>relative</em> to the overall median.</p>
+<p>For example, ‘Healthcare Support’ average income of $43980 can be separated into the median, $45265, and its residual, $-1285.</p>
+<p>So, in addition to gleaning a median income value from this batch of values, we have also devised a way to facilitate the comparison of the income values to one another. Alternatively, we could have assigned colors to each cell of the table reflecting their values <em>relative</em> to the overall median.</p>
 <div class="cell" data-small.mar="true" data-hash="two_way_cache/html/unnamed-chunk-5_4a4743b3e9246aea859a9b3f2be216dc">
 <div class="cell-output-display">
 <p><img src="two_way_files/figure-html/unnamed-chunk-5-1.png" class="img-fluid" width="384"></p>
@@ -785,20 +776,20 @@ <h2 data-number="28.2" class="anchored" data-anchor-id="anatomy-of-a-two-way-tab
 <h2 data-number="28.3" class="anchored" data-anchor-id="analysis-workflow"><span class="header-section-number">28.3</span> Analysis workflow</h2>
 <p>Let’s first create the dataframe.</p>
 <div class="cell" data-hash="two_way_cache/html/unnamed-chunk-6_f5a826f622c2f91a29ff02f720232988">
-<div class="sourceCode cell-code" id="cb2"><pre class="sourceCode r code-with-copy"><code class="sourceCode r"><span id="cb2-1"><a href="#cb2-1" aria-hidden="true" tabindex="-1"></a>df <span class="ot">&lt;-</span> <span class="fu">data.frame</span>(<span class="at">row.names =</span> <span class="fu">c</span>(<span class="st">"NE"</span>,<span class="st">"NC"</span>,<span class="st">"S"</span>,<span class="st">"W"</span>),</span>
-<span id="cb2-2"><a href="#cb2-2" aria-hidden="true" tabindex="-1"></a>                 <span class="at">ed8       =</span> <span class="fu">c</span>(<span class="fl">25.3</span>,<span class="fl">32.1</span>,<span class="fl">38.8</span>,<span class="fl">25.4</span>), </span>
-<span id="cb2-3"><a href="#cb2-3" aria-hidden="true" tabindex="-1"></a>                 <span class="at">ed9to11   =</span> <span class="fu">c</span>(<span class="fl">25.3</span>,<span class="dv">29</span>,<span class="dv">31</span>,<span class="fl">21.1</span>),</span>
-<span id="cb2-4"><a href="#cb2-4" aria-hidden="true" tabindex="-1"></a>                 <span class="at">ed12      =</span> <span class="fu">c</span>(<span class="fl">18.2</span>,<span class="fl">18.8</span>,<span class="fl">19.3</span>,<span class="fl">20.3</span>),</span>
-<span id="cb2-5"><a href="#cb2-5" aria-hidden="true" tabindex="-1"></a>                 <span class="at">ed13to15  =</span> <span class="fu">c</span>(<span class="fl">18.3</span>,<span class="fl">24.3</span>,<span class="fl">15.7</span>,<span class="dv">24</span>),</span>
-<span id="cb2-6"><a href="#cb2-6" aria-hidden="true" tabindex="-1"></a>                 <span class="at">ed16      =</span> <span class="fu">c</span>(<span class="fl">16.3</span>,<span class="dv">19</span>,<span class="fl">16.8</span>,<span class="fl">17.5</span>)</span>
-<span id="cb2-7"><a href="#cb2-7" aria-hidden="true" tabindex="-1"></a>                 )</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
+<div class="sourceCode cell-code" id="cb1"><pre class="sourceCode r code-with-copy"><code class="sourceCode r"><span id="cb1-1"><a href="#cb1-1" aria-hidden="true" tabindex="-1"></a>df <span class="ot">&lt;-</span> <span class="fu">data.frame</span>(<span class="at">row.names =</span> <span class="fu">c</span>(<span class="st">"NE"</span>,<span class="st">"NC"</span>,<span class="st">"S"</span>,<span class="st">"W"</span>),</span>
+<span id="cb1-2"><a href="#cb1-2" aria-hidden="true" tabindex="-1"></a>                 <span class="at">ed8       =</span> <span class="fu">c</span>(<span class="fl">25.3</span>,<span class="fl">32.1</span>,<span class="fl">38.8</span>,<span class="fl">25.4</span>), </span>
+<span id="cb1-3"><a href="#cb1-3" aria-hidden="true" tabindex="-1"></a>                 <span class="at">ed9to11   =</span> <span class="fu">c</span>(<span class="fl">25.3</span>,<span class="dv">29</span>,<span class="dv">31</span>,<span class="fl">21.1</span>),</span>
+<span id="cb1-4"><a href="#cb1-4" aria-hidden="true" tabindex="-1"></a>                 <span class="at">ed12      =</span> <span class="fu">c</span>(<span class="fl">18.2</span>,<span class="fl">18.8</span>,<span class="fl">19.3</span>,<span class="fl">20.3</span>),</span>
+<span id="cb1-5"><a href="#cb1-5" aria-hidden="true" tabindex="-1"></a>                 <span class="at">ed13to15  =</span> <span class="fu">c</span>(<span class="fl">18.3</span>,<span class="fl">24.3</span>,<span class="fl">15.7</span>,<span class="dv">24</span>),</span>
+<span id="cb1-6"><a href="#cb1-6" aria-hidden="true" tabindex="-1"></a>                 <span class="at">ed16      =</span> <span class="fu">c</span>(<span class="fl">16.3</span>,<span class="dv">19</span>,<span class="fl">16.8</span>,<span class="fl">17.5</span>)</span>
+<span id="cb1-7"><a href="#cb1-7" aria-hidden="true" tabindex="-1"></a>                 )</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
 </div>
 <section id="visualizing-the-data" class="level3" data-number="28.3.1">
 <h3 data-number="28.3.1" class="anchored" data-anchor-id="visualizing-the-data"><span class="header-section-number">28.3.1</span> Visualizing the data</h3>
 <p>It’s often easier to look at a graphical representation of the data than a tabular one. Even a table as small as this can benefit from a plot.</p>
 <p>We will adopt Bill Cleveland’s dotplot for this purpose. R has a built-in dotplot function called <code>dotchart</code>. It requires that the table be stored as a matrix instead of a dataframe; we will therefore convert <code>df</code> to a matrix by wrapping it with the <code>as.matrix()</code> function.</p>
 <div class="cell" data-small.mar="true" data-hash="two_way_cache/html/unnamed-chunk-7_14d1690872167f9cd651a876a342ea16">
-<div class="sourceCode cell-code" id="cb3"><pre class="sourceCode r code-with-copy"><code class="sourceCode r"><span id="cb3-1"><a href="#cb3-1" aria-hidden="true" tabindex="-1"></a><span class="fu">dotchart</span>( <span class="fu">as.matrix</span>(df), <span class="at">cex=</span><span class="fl">0.7</span>)</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
+<div class="sourceCode cell-code" id="cb2"><pre class="sourceCode r code-with-copy"><code class="sourceCode r"><span id="cb2-1"><a href="#cb2-1" aria-hidden="true" tabindex="-1"></a><span class="fu">dotchart</span>( <span class="fu">as.matrix</span>(df), <span class="at">cex=</span><span class="fl">0.7</span>)</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
 <div class="cell-output-display">
 <p><img src="two_way_files/figure-html/unnamed-chunk-7-1.png" class="img-fluid" width="288"></p>
 </div>
@@ -806,7 +797,7 @@ <h3 data-number="28.3.1" class="anchored" data-anchor-id="visualizing-the-data">
 <p>The plot helps visualize any differences in mortality rates across different father educational attainment levels. There seems to be a gradual decrease in child mortality with increasing father educational attainment.</p>
 <p>But the plot does not help spot differences across regions (except for the <code>ed12</code> group). We can generate another plot where region becomes the main grouping factor. We do this by wrapping the matrix with the transpose function <code>t()</code>.</p>
 <div class="cell" data-small.mar="true" data-hash="two_way_cache/html/unnamed-chunk-8_5a724aeed8ecace6e93a367317113d23">
-<div class="sourceCode cell-code" id="cb4"><pre class="sourceCode r code-with-copy"><code class="sourceCode r"><span id="cb4-1"><a href="#cb4-1" aria-hidden="true" tabindex="-1"></a><span class="fu">dotchart</span>( <span class="fu">t</span>(<span class="fu">as.matrix</span>(df)), <span class="at">cex=</span><span class="fl">0.7</span>)</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
+<div class="sourceCode cell-code" id="cb3"><pre class="sourceCode r code-with-copy"><code class="sourceCode r"><span id="cb3-1"><a href="#cb3-1" aria-hidden="true" tabindex="-1"></a><span class="fu">dotchart</span>( <span class="fu">t</span>(<span class="fu">as.matrix</span>(df)), <span class="at">cex=</span><span class="fl">0.7</span>)</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
 <div class="cell-output-display">
 <p><img src="two_way_files/figure-html/unnamed-chunk-8-1.png" class="img-fluid" width="288"></p>
 </div>
@@ -814,18 +805,18 @@ <h3 data-number="28.3.1" class="anchored" data-anchor-id="visualizing-the-data">
 <p>At first glance, there seems to be higher death rates for the <code>NC</code> and <code>S</code> regions and relatively lower rates for the <code>W</code> and <code>NE</code> regions. But our eyes may be fooled by outliers in the data.</p>
 <p>Next, we’ll generate side-by-side boxplots to compare the effects between both categories. Note that we’ll need to create a <em>long</em> version of the table using the <code>tidyr</code> package.</p>
 <div class="cell" data-small.mar="true" data-hash="two_way_cache/html/unnamed-chunk-9_42f9c831a6011131f3d0a7e6519ebf81">
-<div class="sourceCode cell-code" id="cb5"><pre class="sourceCode r code-with-copy"><code class="sourceCode r"><span id="cb5-1"><a href="#cb5-1" aria-hidden="true" tabindex="-1"></a><span class="fu">library</span>(tidyr)</span>
-<span id="cb5-2"><a href="#cb5-2" aria-hidden="true" tabindex="-1"></a><span class="fu">library</span>(dplyr)</span>
-<span id="cb5-3"><a href="#cb5-3" aria-hidden="true" tabindex="-1"></a></span>
-<span id="cb5-4"><a href="#cb5-4" aria-hidden="true" tabindex="-1"></a>df.l <span class="ot">&lt;-</span> df <span class="sc">%&gt;%</span></span>
-<span id="cb5-5"><a href="#cb5-5" aria-hidden="true" tabindex="-1"></a>  <span class="fu">mutate</span>(<span class="at">Region =</span> <span class="fu">as.factor</span>(<span class="fu">row.names</span>(.)) )  <span class="sc">%&gt;%</span></span>
-<span id="cb5-6"><a href="#cb5-6" aria-hidden="true" tabindex="-1"></a>  <span class="fu">pivot_longer</span>(<span class="at">names_to =</span> <span class="st">"Edu"</span>, <span class="at">values_to =</span> <span class="st">"Deaths"</span>, <span class="sc">-</span>Region) <span class="sc">%&gt;%</span> </span>
-<span id="cb5-7"><a href="#cb5-7" aria-hidden="true" tabindex="-1"></a>  <span class="fu">mutate</span>(<span class="at">Edu =</span> <span class="fu">factor</span>(Edu, <span class="at">levels=</span><span class="fu">names</span>(df)))</span>
-<span id="cb5-8"><a href="#cb5-8" aria-hidden="true" tabindex="-1"></a></span>
-<span id="cb5-9"><a href="#cb5-9" aria-hidden="true" tabindex="-1"></a><span class="co"># side-by-side plot</span></span>
-<span id="cb5-10"><a href="#cb5-10" aria-hidden="true" tabindex="-1"></a>OP <span class="ot">&lt;-</span> <span class="fu">par</span>(<span class="at">mfrow=</span><span class="fu">c</span>(<span class="dv">1</span>,<span class="dv">2</span>))</span>
-<span id="cb5-11"><a href="#cb5-11" aria-hidden="true" tabindex="-1"></a><span class="fu">plot</span>(Deaths <span class="sc">~</span> Region <span class="sc">+</span> Edu, df.l)</span>
-<span id="cb5-12"><a href="#cb5-12" aria-hidden="true" tabindex="-1"></a><span class="fu">par</span>(OP)</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
+<div class="sourceCode cell-code" id="cb4"><pre class="sourceCode r code-with-copy"><code class="sourceCode r"><span id="cb4-1"><a href="#cb4-1" aria-hidden="true" tabindex="-1"></a><span class="fu">library</span>(tidyr)</span>
+<span id="cb4-2"><a href="#cb4-2" aria-hidden="true" tabindex="-1"></a><span class="fu">library</span>(dplyr)</span>
+<span id="cb4-3"><a href="#cb4-3" aria-hidden="true" tabindex="-1"></a></span>
+<span id="cb4-4"><a href="#cb4-4" aria-hidden="true" tabindex="-1"></a>df.l <span class="ot">&lt;-</span> df <span class="sc">%&gt;%</span></span>
+<span id="cb4-5"><a href="#cb4-5" aria-hidden="true" tabindex="-1"></a>  <span class="fu">mutate</span>(<span class="at">Region =</span> <span class="fu">as.factor</span>(<span class="fu">row.names</span>(.)) )  <span class="sc">%&gt;%</span></span>
+<span id="cb4-6"><a href="#cb4-6" aria-hidden="true" tabindex="-1"></a>  <span class="fu">pivot_longer</span>(<span class="at">names_to =</span> <span class="st">"Edu"</span>, <span class="at">values_to =</span> <span class="st">"Deaths"</span>, <span class="sc">-</span>Region) <span class="sc">%&gt;%</span> </span>
+<span id="cb4-7"><a href="#cb4-7" aria-hidden="true" tabindex="-1"></a>  <span class="fu">mutate</span>(<span class="at">Edu =</span> <span class="fu">factor</span>(Edu, <span class="at">levels=</span><span class="fu">names</span>(df)))</span>
+<span id="cb4-8"><a href="#cb4-8" aria-hidden="true" tabindex="-1"></a></span>
+<span id="cb4-9"><a href="#cb4-9" aria-hidden="true" tabindex="-1"></a><span class="co"># side-by-side plot</span></span>
+<span id="cb4-10"><a href="#cb4-10" aria-hidden="true" tabindex="-1"></a>OP <span class="ot">&lt;-</span> <span class="fu">par</span>(<span class="at">mfrow=</span><span class="fu">c</span>(<span class="dv">1</span>,<span class="dv">2</span>))</span>
+<span id="cb4-11"><a href="#cb4-11" aria-hidden="true" tabindex="-1"></a><span class="fu">plot</span>(Deaths <span class="sc">~</span> Region <span class="sc">+</span> Edu, df.l)</span>
+<span id="cb4-12"><a href="#cb4-12" aria-hidden="true" tabindex="-1"></a><span class="fu">par</span>(OP)</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
 <div class="cell-output-display">
 <p><img src="two_way_files/figure-html/unnamed-chunk-9-1.png" class="img-fluid" width="576"></p>
 </div>
@@ -882,13 +873,18 @@ <h4 data-number="28.3.2.8" class="anchored" data-anchor-id="when-do-we-stop-iter
 <p><img src="img/twoway_08.jpg" class="img-fluid"></p>
 <p>The final version of the table (along with the column and row values) is shown below:</p>
 <p><img src="img/twoway_09.jpg" class="img-fluid"></p>
-<p>A progression of each iteration is shown side-by-side in the next figure. Note how, at the end of the iterations, the <em>common factor</em> ends up contributing the largest effect. Also note the one abnormally high residual in the southern region associated with 8<sup>th</sup> grade education.</p>
+<p>A progression of each iteration is shown side-by-side in the next figure. The left-most table is the original data showing death rates. The second table shows the outcome of the first round of polishing (including the initial overall median value of 20.2). The third and forth table show the second and third iterations of the smoothing operations. Additional iterations are not deemed necessary given that little more can be extracted from the residuals.</p>
 <div class="cell" data-hash="two_way_cache/html/unnamed-chunk-10_90374ca398847c87af3318ada6e07f19">
 <div class="cell-output-display">
 <p><img src="two_way_files/figure-html/unnamed-chunk-10-1.png" class="img-fluid" width="1056"></p>
 </div>
 </div>
-<p>The left-most table is the original data showing death rates. The second table shows the outcome of the first round of polishing (including the initial overall median value of 20.2). The third and forth table show the second and third iterations of the smoothing operations. Additional iterations are not deemed necessary given that little more can be extracted from the residuals.</p>
+<p>A few takeway observations from the analysis follow:</p>
+<ul>
+<li>The <em>common factor</em> ends up contributing the largest effect.</li>
+<li>Between <code>Region</code> and <code>Education</code> effects, the latter explains more of the variability in deaths with a range of 7.43 - (-3.7) = 11.13 deaths per 1000 births. The <code>Region</code> effect explains about 2.3 - (-1.46) = 3.76 deaths per 1000 births of the variability in the data.</li>
+<li>One abnormally high residual–10.97 deaths above the overall median deaths–is associated with the southern region and an 8<sup>th</sup> grade education.</li>
+</ul>
 </section>
 </section>
 <section id="interpreting-the-median-polish" class="level3" data-number="28.3.3">
@@ -924,7 +920,7 @@ <h2 data-number="28.4" class="anchored" data-anchor-id="implementing-the-median-
 <h3 data-number="28.4.1" class="anchored" data-anchor-id="using-base-rs-medpolish"><span class="header-section-number">28.4.1</span> Using base R’s <code>medpolish</code></h3>
 <p>The steps outlined in the previous section can be easily implemented using pen and paper or a spreadsheet environment for larger datasets. R has a built-in function called <code>medpolish()</code> which does this for us. We can define the maximum number of iterations by setting the <code>maxiter=</code> parameter but note that <code>medpolish</code> will, by default, automatically estimate the best number of iterations for us.</p>
 <div class="cell" data-hash="two_way_cache/html/unnamed-chunk-12_b33d0905705eb525559f39c2bdc77e3f">
-<div class="sourceCode cell-code" id="cb6"><pre class="sourceCode r code-with-copy"><code class="sourceCode r"><span id="cb6-1"><a href="#cb6-1" aria-hidden="true" tabindex="-1"></a>df.med <span class="ot">&lt;-</span> <span class="fu">medpolish</span>( df)</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
+<div class="sourceCode cell-code" id="cb5"><pre class="sourceCode r code-with-copy"><code class="sourceCode r"><span id="cb5-1"><a href="#cb5-1" aria-hidden="true" tabindex="-1"></a>df.med <span class="ot">&lt;-</span> <span class="fu">medpolish</span>( df)</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
 <div class="cell-output cell-output-stdout">
 <pre><code>1: 47.1
 2: 42.9
@@ -935,33 +931,33 @@ <h3 data-number="28.4.1" class="anchored" data-anchor-id="using-base-rs-medpolis
 <p>The median polish ran through 4 iterations–one more than was run in our step-by-step example.</p>
 <p>The four values printed to the console gives us the sum of absolute residuals at each iteration. The output from the model is stored in the <code>df.med</code> object. To see its contents , simply type it at the command line.</p>
 <div class="cell" data-hash="two_way_cache/html/unnamed-chunk-13_32619a5000872bc3f93103779d81e8c4">
-<div class="sourceCode cell-code" id="cb8"><pre class="sourceCode r code-with-copy"><code class="sourceCode r"><span id="cb8-1"><a href="#cb8-1" aria-hidden="true" tabindex="-1"></a>df.med</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
+<div class="sourceCode cell-code" id="cb7"><pre class="sourceCode r code-with-copy"><code class="sourceCode r"><span id="cb7-1"><a href="#cb7-1" aria-hidden="true" tabindex="-1"></a>df.med</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
 <div class="cell-output cell-output-stdout">
 <pre><code>
 Median Polish Results (Dataset: "df")
 
-Overall: 20.6
+Overall: 20.775
 
 Row Effects:
-   NE    NC     S     W 
--1.50  2.55 -0.35  0.35 
+    NE     NC      S      W 
+-1.475  2.375 -0.350  0.350 
 
 Column Effects:
      ed8  ed9to11     ed12 ed13to15     ed16 
-   7.575    6.025   -0.925    0.175   -3.450 
+  7.4750   5.9250  -1.1125   0.0750  -3.6250 
 
 Residuals:
-      ed8 ed9to11   ed12 ed13to15  ed16
-NE -1.375   0.175  0.025   -0.975  0.65
-NC  1.375  -0.175 -3.425    0.975 -0.70
-S  10.975   4.725 -0.025   -4.725  0.00
-W  -3.125  -5.875  0.275    2.875  0.00</code></pre>
+      ed8 ed9to11    ed12 ed13to15   ed16
+NE -1.475   0.075  0.0125   -1.075  0.625
+NC  1.475  -0.075 -3.2375    1.075 -0.525
+S  10.900   4.650 -0.0125   -4.800  0.000
+W  -3.200  -5.950  0.2875    2.800  0.000</code></pre>
 </div>
 </div>
 <p>All three effects are displayed as well as the residuals (note that the precision returned is greater than that used in our earlier analysis).</p>
 <p>To generate the Tukey additivity plot, simply wrap the median polish object with the <code>plot</code> command as in:</p>
 <div class="cell" data-small.mar="true" data-hash="two_way_cache/html/unnamed-chunk-14_09c5df91ce09cbcd4b168d3429c6f89a">
-<div class="sourceCode cell-code" id="cb10"><pre class="sourceCode r code-with-copy"><code class="sourceCode r"><span id="cb10-1"><a href="#cb10-1" aria-hidden="true" tabindex="-1"></a><span class="fu">plot</span>( df.med )</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
+<div class="sourceCode cell-code" id="cb9"><pre class="sourceCode r code-with-copy"><code class="sourceCode r"><span id="cb9-1"><a href="#cb9-1" aria-hidden="true" tabindex="-1"></a><span class="fu">plot</span>( df.med )</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
 <div class="cell-output-display">
 <p><img src="two_way_files/figure-html/unnamed-chunk-14-1.png" class="img-fluid" width="288"></p>
 </div>
@@ -972,30 +968,30 @@ <h3 data-number="28.4.2" class="anchored" data-anchor-id="using-eda_pol"><span c
 <p>The <code>tukeyedar</code> package has a custom function called <code>eda_pol</code> that will output the polished table as a graphic element.</p>
 <p>Unlike the base <code>medpolish</code> function, <code>eda_pol</code> requires that the data be in long form. In creating the long form version of the data, we’ll explicitly set the levels for educational attainment and region.</p>
 <div class="cell" data-hash="two_way_cache/html/unnamed-chunk-15_d311464cbe1865576659f04c00e2fe0b">
-<div class="sourceCode cell-code" id="cb11"><pre class="sourceCode r code-with-copy"><code class="sourceCode r"><span id="cb11-1"><a href="#cb11-1" aria-hidden="true" tabindex="-1"></a>grd <span class="ot">&lt;-</span> <span class="fu">c</span>(<span class="st">"ed8"</span>, <span class="st">"ed9-11"</span>, <span class="st">"ed12"</span>, <span class="st">"ed13-15"</span>, <span class="st">"ed16"</span>)</span>
-<span id="cb11-2"><a href="#cb11-2" aria-hidden="true" tabindex="-1"></a>reg <span class="ot">&lt;-</span> <span class="fu">c</span>(<span class="st">"NE"</span>, <span class="st">"NC"</span>, <span class="st">"S"</span>, <span class="st">"W"</span>)</span>
-<span id="cb11-3"><a href="#cb11-3" aria-hidden="true" tabindex="-1"></a>dfl <span class="ot">&lt;-</span> <span class="fu">data.frame</span>(<span class="at">Region =</span>  <span class="fu">factor</span>(<span class="fu">rep</span>( <span class="fu">c</span>(<span class="st">"NE"</span>, <span class="st">"NC"</span>, <span class="st">"S"</span>, <span class="st">"W"</span>), <span class="at">each =</span> <span class="dv">5</span>),</span>
-<span id="cb11-4"><a href="#cb11-4" aria-hidden="true" tabindex="-1"></a>                                   <span class="at">levels =</span> reg),</span>
-<span id="cb11-5"><a href="#cb11-5" aria-hidden="true" tabindex="-1"></a>                 <span class="at">Education =</span> <span class="fu">factor</span>(<span class="fu">rep</span>( grd , <span class="dv">4</span>), <span class="at">levels =</span> grd),</span>
-<span id="cb11-6"><a href="#cb11-6" aria-hidden="true" tabindex="-1"></a>                 <span class="at">Deaths =</span> <span class="fu">c</span>(<span class="fl">25.3</span>, <span class="fl">25.3</span>, <span class="fl">18.2</span>, <span class="fl">18.3</span>, <span class="fl">16.3</span>, <span class="fl">32.1</span>, <span class="dv">29</span>, <span class="fl">18.8</span>,</span>
-<span id="cb11-7"><a href="#cb11-7" aria-hidden="true" tabindex="-1"></a>                          <span class="fl">24.3</span>, <span class="dv">19</span>, <span class="fl">38.8</span>, <span class="dv">31</span>, <span class="fl">19.3</span>, <span class="fl">15.7</span>, <span class="fl">16.8</span>, <span class="fl">25.4</span>, </span>
-<span id="cb11-8"><a href="#cb11-8" aria-hidden="true" tabindex="-1"></a>                          <span class="fl">21.1</span>, <span class="fl">20.3</span>, <span class="dv">24</span>, <span class="fl">17.5</span>))</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
+<div class="sourceCode cell-code" id="cb10"><pre class="sourceCode r code-with-copy"><code class="sourceCode r"><span id="cb10-1"><a href="#cb10-1" aria-hidden="true" tabindex="-1"></a>grd <span class="ot">&lt;-</span> <span class="fu">c</span>(<span class="st">"ed8"</span>, <span class="st">"ed9-11"</span>, <span class="st">"ed12"</span>, <span class="st">"ed13-15"</span>, <span class="st">"ed16"</span>)</span>
+<span id="cb10-2"><a href="#cb10-2" aria-hidden="true" tabindex="-1"></a>reg <span class="ot">&lt;-</span> <span class="fu">c</span>(<span class="st">"NE"</span>, <span class="st">"NC"</span>, <span class="st">"S"</span>, <span class="st">"W"</span>)</span>
+<span id="cb10-3"><a href="#cb10-3" aria-hidden="true" tabindex="-1"></a>dfl <span class="ot">&lt;-</span> <span class="fu">data.frame</span>(<span class="at">Region =</span>  <span class="fu">factor</span>(<span class="fu">rep</span>( <span class="fu">c</span>(<span class="st">"NE"</span>, <span class="st">"NC"</span>, <span class="st">"S"</span>, <span class="st">"W"</span>), <span class="at">each =</span> <span class="dv">5</span>),</span>
+<span id="cb10-4"><a href="#cb10-4" aria-hidden="true" tabindex="-1"></a>                                   <span class="at">levels =</span> reg),</span>
+<span id="cb10-5"><a href="#cb10-5" aria-hidden="true" tabindex="-1"></a>                 <span class="at">Education =</span> <span class="fu">factor</span>(<span class="fu">rep</span>( grd , <span class="dv">4</span>), <span class="at">levels =</span> grd),</span>
+<span id="cb10-6"><a href="#cb10-6" aria-hidden="true" tabindex="-1"></a>                 <span class="at">Deaths =</span> <span class="fu">c</span>(<span class="fl">25.3</span>, <span class="fl">25.3</span>, <span class="fl">18.2</span>, <span class="fl">18.3</span>, <span class="fl">16.3</span>, <span class="fl">32.1</span>, <span class="dv">29</span>, <span class="fl">18.8</span>,</span>
+<span id="cb10-7"><a href="#cb10-7" aria-hidden="true" tabindex="-1"></a>                          <span class="fl">24.3</span>, <span class="dv">19</span>, <span class="fl">38.8</span>, <span class="dv">31</span>, <span class="fl">19.3</span>, <span class="fl">15.7</span>, <span class="fl">16.8</span>, <span class="fl">25.4</span>, </span>
+<span id="cb10-8"><a href="#cb10-8" aria-hidden="true" tabindex="-1"></a>                          <span class="fl">21.1</span>, <span class="fl">20.3</span>, <span class="dv">24</span>, <span class="fl">17.5</span>))</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
 </div>
 <p>The <code>eda_pol</code> requires that you specify which variables will be mapped to the row and column effects.</p>
 <div class="cell" data-hash="two_way_cache/html/unnamed-chunk-16_e9e1b42d5c1c6d97e0de8cc86a753169">
-<div class="sourceCode cell-code" id="cb12"><pre class="sourceCode r eda code-with-copy"><code class="sourceCode r"><span id="cb12-1"><a href="#cb12-1" aria-hidden="true" tabindex="-1"></a><span class="fu">library</span>(tukeyedar)</span>
-<span id="cb12-2"><a href="#cb12-2" aria-hidden="true" tabindex="-1"></a></span>
-<span id="cb12-3"><a href="#cb12-3" aria-hidden="true" tabindex="-1"></a>df.pol <span class="ot">&lt;-</span> <span class="fu">eda_pol</span>(dfl, <span class="at">row =</span> Region, <span class="at">col =</span> Education , <span class="at">val =</span> Deaths )</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
+<div class="sourceCode cell-code" id="cb11"><pre class="sourceCode r eda code-with-copy"><code class="sourceCode r"><span id="cb11-1"><a href="#cb11-1" aria-hidden="true" tabindex="-1"></a><span class="fu">library</span>(tukeyedar)</span>
+<span id="cb11-2"><a href="#cb11-2" aria-hidden="true" tabindex="-1"></a></span>
+<span id="cb11-3"><a href="#cb11-3" aria-hidden="true" tabindex="-1"></a>df.pol <span class="ot">&lt;-</span> <span class="fu">eda_pol</span>(dfl, <span class="at">row =</span> Region, <span class="at">col =</span> Education , <span class="at">val =</span> Deaths )</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
 </div>
 <div class="cell" data-hash="two_way_cache/html/unnamed-chunk-17_0485061db5138b5f8809db9973d13722">
 <div class="cell-output-display">
 <p><img src="two_way_files/figure-html/unnamed-chunk-17-1.png" class="img-fluid" width="412"></p>
 </div>
 </div>
-<p>Here, the <code>eda_pol</code> function required five iterations which explains the slight difference in values from those generated with the medpolish function. You can control the maximum number of iterations in <code>eda_pol</code> with the <code>maxiter</code> argument.</p>
+<p>Here, the <code>eda_pol</code> function required five iterations which explains the slight difference in values from those generated with the <code>medpolish</code> function. You can control the maximum number of iterations in <code>eda_pol</code> with the <code>maxiter</code> argument.</p>
 <p>By default, the column and row values are listed in the same order as they are presented in the input data frame (either defined alphabetically or defined by their levels). You can have the function sort the values by their effect size using the <code>sort</code> argument.</p>
 <div class="cell" data-hash="two_way_cache/html/unnamed-chunk-18_a75fa6025b012257c3f986da05c23b65">
-<div class="sourceCode cell-code" id="cb13"><pre class="sourceCode r eda code-with-copy"><code class="sourceCode r"><span id="cb13-1"><a href="#cb13-1" aria-hidden="true" tabindex="-1"></a>df.pol <span class="ot">&lt;-</span> <span class="fu">eda_pol</span>(dfl, <span class="at">row =</span> Region, <span class="at">col =</span> Education, <span class="at">val =</span> Deaths, <span class="at">sort =</span> <span class="cn">TRUE</span> )</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
+<div class="sourceCode cell-code" id="cb12"><pre class="sourceCode r eda code-with-copy"><code class="sourceCode r"><span id="cb12-1"><a href="#cb12-1" aria-hidden="true" tabindex="-1"></a>df.pol <span class="ot">&lt;-</span> <span class="fu">eda_pol</span>(dfl, <span class="at">row =</span> Region, <span class="at">col =</span> Education, <span class="at">val =</span> Deaths, <span class="at">sort =</span> <span class="cn">TRUE</span> )</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
 </div>
 <div class="cell" data-hash="two_way_cache/html/unnamed-chunk-19_a4c8ef0614833d803082b2ce2a30b5b3">
 <div class="cell-output-display">
@@ -1004,11 +1000,11 @@ <h3 data-number="28.4.2" class="anchored" data-anchor-id="using-eda_pol"><span c
 </div>
 <p>The <code>df.pol</code> object is a list that stores output values such as the common, row and column effects.</p>
 <div class="cell" data-hash="two_way_cache/html/unnamed-chunk-20_c867fa9dd2b6b52aafdcfef896e2312d">
-<div class="sourceCode cell-code" id="cb14"><pre class="sourceCode r code-with-copy"><code class="sourceCode r"><span id="cb14-1"><a href="#cb14-1" aria-hidden="true" tabindex="-1"></a>df.pol<span class="sc">$</span>global</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
+<div class="sourceCode cell-code" id="cb13"><pre class="sourceCode r code-with-copy"><code class="sourceCode r"><span id="cb13-1"><a href="#cb13-1" aria-hidden="true" tabindex="-1"></a>df.pol<span class="sc">$</span>global</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
 <div class="cell-output cell-output-stdout">
 <pre><code>[1] 20.85</code></pre>
 </div>
-<div class="sourceCode cell-code" id="cb16"><pre class="sourceCode r code-with-copy"><code class="sourceCode r"><span id="cb16-1"><a href="#cb16-1" aria-hidden="true" tabindex="-1"></a>df.pol<span class="sc">$</span>row</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
+<div class="sourceCode cell-code" id="cb15"><pre class="sourceCode r code-with-copy"><code class="sourceCode r"><span id="cb15-1"><a href="#cb15-1" aria-hidden="true" tabindex="-1"></a>df.pol<span class="sc">$</span>row</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
 <div class="cell-output cell-output-stdout">
 <pre><code>  Region  effect
 1     NE -1.4625
@@ -1016,7 +1012,7 @@ <h3 data-number="28.4.2" class="anchored" data-anchor-id="using-eda_pol"><span c
 3      S -0.3500
 4      W  0.3500</code></pre>
 </div>
-<div class="sourceCode cell-code" id="cb18"><pre class="sourceCode r code-with-copy"><code class="sourceCode r"><span id="cb18-1"><a href="#cb18-1" aria-hidden="true" tabindex="-1"></a>df.pol<span class="sc">$</span>col</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
+<div class="sourceCode cell-code" id="cb17"><pre class="sourceCode r code-with-copy"><code class="sourceCode r"><span id="cb17-1"><a href="#cb17-1" aria-hidden="true" tabindex="-1"></a>df.pol<span class="sc">$</span>col</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
 <div class="cell-output cell-output-stdout">
 <pre><code>  Education   effect
 1       ed8  7.43125
@@ -1028,7 +1024,7 @@ <h3 data-number="28.4.2" class="anchored" data-anchor-id="using-eda_pol"><span c
 </div>
 <p>The <code>eda_pol</code> will also output a diagnostic plot.</p>
 <div class="cell" data-hash="two_way_cache/html/unnamed-chunk-21_da3b4d2d15c439468dd98d9b3e10f4db">
-<div class="sourceCode cell-code" id="cb20"><pre class="sourceCode r eda code-with-copy"><code class="sourceCode r"><span id="cb20-1"><a href="#cb20-1" aria-hidden="true" tabindex="-1"></a><span class="fu">plot</span>(df.pol, <span class="at">type =</span> <span class="st">"diagnostic"</span>)</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
+<div class="sourceCode cell-code" id="cb19"><pre class="sourceCode r eda code-with-copy"><code class="sourceCode r"><span id="cb19-1"><a href="#cb19-1" aria-hidden="true" tabindex="-1"></a><span class="fu">plot</span>(df.pol, <span class="at">type =</span> <span class="st">"diagnostic"</span>)</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
 <div class="cell-output cell-output-stdout">
 <pre><code>$slope
     cv 
@@ -1038,7 +1034,7 @@ <h3 data-number="28.4.2" class="anchored" data-anchor-id="using-eda_pol"><span c
 <p><img src="two_way_files/figure-html/unnamed-chunk-21-1.png" class="img-fluid" width="288"></p>
 </div>
 </div>
-<p>In addition to generating the plot, the function will also return the slope of the fitted line.</p>
+<p>In addition to generating the plot, the function will also return the slope of the fitted line in the console.</p>
 <p>For a comprehensive review of the median polish and the <code>eda_pol</code> function, see the <code>eda_pol</code> <a href="https://mgimond.github.io/tukeyedar/articles/polish.html">vignette</a>.</p>
 </section>
 </section>
@@ -1054,13 +1050,13 @@ <h3 data-number="28.5.1" class="anchored" data-anchor-id="two-way-anova"><span c
 <p><img src="img/twoway_12.png" class="img-fluid"></p>
 <p>The resulting row/column levels are the regression coefficients:</p>
 <div class="cell" data-hash="two_way_cache/html/unnamed-chunk-22_9857b00948e76646f2b42871c3c01715">
-<div class="sourceCode cell-code" id="cb22"><pre class="sourceCode r code-with-copy"><code class="sourceCode r"><span id="cb22-1"><a href="#cb22-1" aria-hidden="true" tabindex="-1"></a>M <span class="ot">&lt;-</span> <span class="fu">lm</span>(Deaths <span class="sc">~</span> Education <span class="sc">+</span> Region, dfl)</span>
-<span id="cb22-2"><a href="#cb22-2" aria-hidden="true" tabindex="-1"></a><span class="fu">coef</span>(M)</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
+<div class="sourceCode cell-code" id="cb21"><pre class="sourceCode r code-with-copy"><code class="sourceCode r"><span id="cb21-1"><a href="#cb21-1" aria-hidden="true" tabindex="-1"></a>M <span class="ot">&lt;-</span> <span class="fu">lm</span>(Deaths <span class="sc">~</span> Education <span class="sc">+</span> Region, dfl)</span>
+<span id="cb21-2"><a href="#cb21-2" aria-hidden="true" tabindex="-1"></a><span class="fu">coef</span>(M)</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
 <div class="cell-output cell-output-stdout">
-<pre><code>     (Intercept)  Educationed9-11    Educationed12 Educationed13-15    Educationed16         RegionNE          RegionS 
-          32.215           -3.800          -11.250           -9.825          -13.000           -3.960           -0.320 
+<pre><code>     (Intercept)  Educationed9-11    Educationed12 Educationed13-15    Educationed16         RegionNC          RegionS 
+          28.255           -3.800          -11.250           -9.825          -13.000            3.960            3.640 
          RegionW 
-          -2.980 </code></pre>
+           0.980 </code></pre>
 </div>
 </div>
 <p>The regression formula takes on the form:</p>
@@ -1081,13 +1077,13 @@ <h3 data-number="28.5.1" class="anchored" data-anchor-id="two-way-anova"><span c
 \]</span></p>
 <p>To assess if any of the factors have a significant effect on the <code>Death</code> variables, simply wrap the regression with an <code>anova()</code> function.</p>
 <div class="cell" data-hash="two_way_cache/html/unnamed-chunk-23_ce64700903bec89d35ce4c76b73761d0">
-<div class="sourceCode cell-code" id="cb24"><pre class="sourceCode r code-with-copy"><code class="sourceCode r"><span id="cb24-1"><a href="#cb24-1" aria-hidden="true" tabindex="-1"></a><span class="fu">anova</span>(M)</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
+<div class="sourceCode cell-code" id="cb23"><pre class="sourceCode r code-with-copy"><code class="sourceCode r"><span id="cb23-1"><a href="#cb23-1" aria-hidden="true" tabindex="-1"></a><span class="fu">anova</span>(M)</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
 <div class="cell-output cell-output-stdout">
 <pre><code>Analysis of Variance Table
 
 Response: Deaths
           Df Sum Sq Mean Sq F value  Pr(&gt;F)   
-Edu        4 478.52 119.630  7.7539 0.00251 **
+Education  4 478.52 119.630  7.7539 0.00251 **
 Region     3  57.44  19.146  1.2410 0.33801   
 Residuals 12 185.14  15.428                   
 ---
@@ -1100,7 +1096,7 @@ <h3 data-number="28.5.1" class="anchored" data-anchor-id="two-way-anova"><span c
 <h3 data-number="28.5.2" class="anchored" data-anchor-id="visualizing-a-two-way-anova-as-a-table"><span class="header-section-number">28.5.2</span> Visualizing a two-way ANOVA as a table</h3>
 <p>You can view the output of a two-way ANOVA as a “mean polish” table using the <code>eda_pol</code> function. Simply set the <code>stat</code> argument to <code>mean</code>.</p>
 <div class="cell" data-hash="two_way_cache/html/unnamed-chunk-24_53efb0f64413d89e047d527343fced23">
-<div class="sourceCode cell-code" id="cb26"><pre class="sourceCode r eda code-with-copy"><code class="sourceCode r"><span id="cb26-1"><a href="#cb26-1" aria-hidden="true" tabindex="-1"></a>df.ano <span class="ot">&lt;-</span> <span class="fu">eda_pol</span>(dfl, <span class="at">row =</span> Region, <span class="at">col =</span> Education, <span class="at">val =</span> Deaths, <span class="at">stat =</span> mean)</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
+<div class="sourceCode cell-code" id="cb25"><pre class="sourceCode r eda code-with-copy"><code class="sourceCode r"><span id="cb25-1"><a href="#cb25-1" aria-hidden="true" tabindex="-1"></a>df.ano <span class="ot">&lt;-</span> <span class="fu">eda_pol</span>(dfl, <span class="at">row =</span> Region, <span class="at">col =</span> Education, <span class="at">val =</span> Deaths, <span class="at">stat =</span> mean)</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
 </div>
 <div class="cell" data-hash="two_way_cache/html/unnamed-chunk-25_b9389ebf43a5cb11a94361e7db6b0a52">
 <div class="cell-output-display">
diff --git a/docs/two_way_files/figure-html/unnamed-chunk-14-1.png b/docs/two_way_files/figure-html/unnamed-chunk-14-1.png
index 2691a1ce..2e33439e 100644
Binary files a/docs/two_way_files/figure-html/unnamed-chunk-14-1.png and b/docs/two_way_files/figure-html/unnamed-chunk-14-1.png differ
diff --git a/example.knit.md b/example.knit.md
new file mode 100644
index 00000000..525282fc
--- /dev/null
+++ b/example.knit.md
@@ -0,0 +1,444 @@
+---
+output: html_document
+editor_options: 
+  chunk_output_type: console
+---
+
+# The third R of EDA: Residuals
+
+
+::: {.cell}
+
+:::
+
+::: {.cell hash='example_cache/html/unnamed-chunk-2_fd93e819e7915059036d51d67790e90a'}
+::: {.cell-output-display}
+`````{=html}
+<table data-quarto-disable-processing="true" class="table" style="width: auto !important; ">
+ <thead>
+  <tr>
+   <th style="text-align:left;color: rgba(85, 85, 85, 255) !important;background-color: rgba(221, 221, 221, 255) !important;text-align: center;border: 1px solid white !important;
+             font-family: 'Source Code Pro', 'Open Sans';
+             padding:1px !important;
+             padding-left:4px !important;
+             padding-right:4px !important;
+             font-size: 0.8em;
+             border-radius: 5px;"> dplyr </th>
+   <th style="text-align:left;color: rgba(85, 85, 85, 255) !important;background-color: rgba(221, 221, 221, 255) !important;text-align: center;border: 1px solid white !important;
+             font-family: 'Source Code Pro', 'Open Sans';
+             padding:1px !important;
+             padding-left:4px !important;
+             padding-right:4px !important;
+             font-size: 0.8em;
+             border-radius: 5px;"> ggplot2 </th>
+   <th style="text-align:left;color: rgba(85, 85, 85, 255) !important;background-color: rgba(221, 221, 221, 255) !important;text-align: center;border: 1px solid white !important;
+             font-family: 'Source Code Pro', 'Open Sans';
+             padding:1px !important;
+             padding-left:4px !important;
+             padding-right:4px !important;
+             font-size: 0.8em;
+             border-radius: 5px;"> tidyr </th>
+   <th style="text-align:left;color: rgba(85, 85, 85, 255) !important;background-color: rgba(221, 221, 221, 255) !important;text-align: center;border: 1px solid white !important;
+             font-family: 'Source Code Pro', 'Open Sans';
+             padding:1px !important;
+             padding-left:4px !important;
+             padding-right:4px !important;
+             font-size: 0.8em;
+             border-radius: 5px;"> lubridate </th>
+  </tr>
+ </thead>
+<tbody>
+  <tr>
+   <td style="text-align:left;color: darkred !important;background-color: rgba(250, 232, 232, 255) !important;text-align: center;border: 1px solid white;
+             font-family: 'Open Sans', Arial;
+             padding:1px !important;
+             padding-left:4px !important;
+             padding-right:4px !important;
+             font-size: 0.8em;
+             border-radius: 5px;"> 1.1.4 </td>
+   <td style="text-align:left;color: darkred !important;background-color: rgba(250, 232, 232, 255) !important;text-align: center;border: 1px solid white;
+             font-family: 'Open Sans', Arial;
+             padding:1px !important;
+             padding-left:4px !important;
+             padding-right:4px !important;
+             font-size: 0.8em;
+             border-radius: 5px;"> 3.4.4 </td>
+   <td style="text-align:left;color: darkred !important;background-color: rgba(250, 232, 232, 255) !important;text-align: center;border: 1px solid white;
+             font-family: 'Open Sans', Arial;
+             padding:1px !important;
+             padding-left:4px !important;
+             padding-right:4px !important;
+             font-size: 0.8em;
+             border-radius: 5px;"> 1.3.1 </td>
+   <td style="text-align:left;color: darkred !important;background-color: rgba(250, 232, 232, 255) !important;text-align: center;border: 1px solid white;
+             font-family: 'Open Sans', Arial;
+             padding:1px !important;
+             padding-left:4px !important;
+             padding-right:4px !important;
+             font-size: 0.8em;
+             border-radius: 5px;"> 1.9.3 </td>
+  </tr>
+</tbody>
+</table>
+
+`````
+:::
+:::
+
+
+So far, we've used residuals to help diagnose model fits in both univariate and bivariate analyses. We've also learned how to leverage re-expression and robust methods (the two other R's of EDA) in helping improve model fits.
+
+In this chapter, we work through an example that highlights the importance of exploring residuals to reveal additional patterns in the data--patterns that may not necessarily be observed in an initial model fit.
+
+## An EDA example: CO2 analysis
+
+Let's start off by plotting the atmospheric CO2 concentrations (in ppm) over time. The original dataset was pulled from [NOAA's website](http://www.esrl.noaa.gov/gmd/ccgg/trends/).
+
+
+::: {.cell hash='example_cache/html/unnamed-chunk-3_8ff541de663eb986fb7aa84d57d9e4a8'}
+
+```{.r .cell-code}
+library(lubridate)
+library(dplyr)
+library(ggplot2)
+
+dat <- read.csv("http://mgimond.github.io/ES218/Data/co2_2020.csv")
+
+dat2 <- dat %>%
+      mutate( Date = ymd(paste(Year,"/", Month,"/15")),
+              Year = decimal_date(Date),
+              CO2 = Average)    %>%
+     select(Year, Date, CO2)
+```
+:::
+
+
+
+Next, let's look at the data:
+      
+
+::: {.cell hash='example_cache/html/unnamed-chunk-4_18cffa525f59a03421c6392b95970664'}
+
+```{.r .cell-code}
+ggplot(dat2, aes(x = Year , y = CO2)) + geom_point(size = 0.6)
+```
+
+::: {.cell-output-display}
+![](example_files/figure-html/unnamed-chunk-4-1.png){width=576}
+:::
+:::
+
+
+Let's see if plotting the data using lines instead of points helps identify any underlying patterns.
+      
+
+::: {.cell hash='example_cache/html/unnamed-chunk-5_55c2b633da90110e2e7c31caf6a481d5'}
+
+```{.r .cell-code}
+ggplot(dat2, aes(x = Year , y = CO2)) + geom_line()
+```
+
+::: {.cell-output-display}
+![](example_files/figure-html/unnamed-chunk-5-1.png){width=576}
+:::
+:::
+
+
+We note two patterns of interest: an overall upward trend, and a cyclical trend. We will first *smooth* out the overall trend by finding a model that best explains the overall variability, then we'll remove that trend from the data to explore the cyclical component.
+
+## Exploring the overall trend
+
+We can attempt to model the overall trend by fitting a straight line to the data using a 1^st^ order polynomial. The fitted line is displayed in red in the following plot. We will use `ggplot2`'s built-in smoothing function to plot the regression line.
+
+
+
+::: {.cell hash='example_cache/html/unnamed-chunk-6_2161d4afd86c476bae6977716b12b4a3'}
+
+```{.r .cell-code}
+ggplot(dat2, aes(x = Year, y = CO2)) + geom_line() + 
+             stat_smooth(method = "lm", se = FALSE, col = "red") 
+```
+
+::: {.cell-output-display}
+![](example_files/figure-html/unnamed-chunk-6-1.png){width=576}
+:::
+:::
+
+
+Next, we subtract the modeled line from the CO2 data and plot the result--this gives us the **residuals** and can be thought of as representing what the linear model does *not* explain about the data. But first, we will need to create a regression model object since it will provide us with the residuals needed to generate the residual plot.
+
+
+::: {.cell hash='example_cache/html/unnamed-chunk-7_a35d55aa4f50ac91d7535d1ce418686b'}
+
+```{.r .cell-code}
+# Generate a model object from the regression fit
+M1 <- lm(CO2 ~ Year, dat2) 
+dat2$res1 <- M1$residuals
+
+# Plot the residuals from the model object
+ggplot(dat2, aes(x = Year, y = res1)) + geom_point()
+```
+
+::: {.cell-output-display}
+![](example_files/figure-html/unnamed-chunk-7-1.png){width=576}
+:::
+:::
+
+
+An overall trend is still present, despite having attempted to control for it. This suggests that our simple linear model does not do a good job in *capturing* the shape of the overall trend. 
+
+It appears that the overall trend is slightly convex and has a small peak around the 1990's. This suggests that we should move on to a 2^nd^ order polynomial of the form:
+
+$$
+CO2_{trend} = a + b(Year) + c(Year^2) 
+$$
+
+The fitted line looks like this: 
+
+
+::: {.cell hash='example_cache/html/unnamed-chunk-8_3a3326053a34672abde81f4ec6b7a77f'}
+
+```{.r .cell-code}
+ggplot(dat2, aes(x = Year, y = CO2)) + geom_line() + 
+           stat_smooth(method = "lm", formula = y ~ x + I(x^2), 
+                       se = FALSE, col = "red")
+```
+
+::: {.cell-output-display}
+![](example_files/figure-html/unnamed-chunk-8-1.png){width=576}
+:::
+:::
+
+
+Now, let's look at the residuals. In helping discern any pattern in our residuals, we add a loess smooth to the plot.
+
+
+::: {.cell hash='example_cache/html/unnamed-chunk-9_737f2732efa31eb98a846acf9540f805'}
+
+```{.r .cell-code}
+M2 <- lm(CO2 ~ Year + I(Year^2), dat2)
+dat2$res2 <- M2$residuals
+
+ggplot(dat2 , aes(x = Year, y = res2)) + geom_point() +
+           stat_smooth(method = "loess", se = FALSE)
+```
+
+::: {.cell-output-display}
+![](example_files/figure-html/unnamed-chunk-9-1.png){width=576}
+:::
+:::
+
+
+This is an improvement over the simple line model (which was a 1^st^ order polynomial fit).  So what we have learned so far is that the overall trend is not perfectly linear but instead follows a parabolic like trajectory with a small hump halfway across the time span. However, we can still make out a "W" shaped trend in the residuals which can hamper our analysis of the oscillatory patterns in the data. We can try a 3^rd^ order polynomial and see if it can capture the 1980 to 1995 hump.
+
+
+
+::: {.cell hash='example_cache/html/unnamed-chunk-10_b5a0f63a1435b46624cbc8528f9f1a7e'}
+
+```{.r .cell-code}
+ggplot(dat2, aes(x = Year, y = CO2)) + geom_line() + 
+           stat_smooth(method = "lm", formula = y ~ x + I(x^2) + I(x^3), 
+                       se = FALSE, col = "red")
+```
+
+::: {.cell-output-display}
+![](example_files/figure-html/unnamed-chunk-10-1.png){width=576}
+:::
+:::
+
+
+Now, let's look at the residuals. 
+
+
+::: {.cell hash='example_cache/html/unnamed-chunk-11_b72a0542c1f88351317336ac149a8bfd'}
+
+```{.r .cell-code}
+M2 <- lm(CO2 ~ Year + I(Year^2) + I(Year^3), dat2)
+dat2$res2 <- M2$residuals
+
+ggplot(dat2 , aes(x = Year, y = res2)) + geom_point() +
+           stat_smooth(method = "loess", se = FALSE)
+```
+
+::: {.cell-output-display}
+![](example_files/figure-html/unnamed-chunk-11-1.png){width=576}
+:::
+:::
+
+
+This does not seem to be an improvement over the 2^nd^ order polynomial. At this point, we could play with different order polynomials in an attempt to smooth the trend or, we may opt for a non-parametric fit. When the goal is to *peel* off one pattern to explore any underlying pattern we should not limit ourselves to parametric fits (which impose a mathematical model on our data) and instead explore non-parametric smoothing techniques that do not impose any structure on the data whatsoever. We'll make use of the loess fit which is shown in the following figure.
+
+
+::: {.cell hash='example_cache/html/unnamed-chunk-12_83a89c9aa33856f9a0d1116388596410'}
+
+```{.r .cell-code}
+ggplot(dat2, aes(x = Year, y = CO2)) + geom_line() + 
+             stat_smooth(method = "loess", span = 1/4, se = FALSE, col = "red")
+```
+
+::: {.cell-output-display}
+![](example_files/figure-html/unnamed-chunk-12-1.png){width=576}
+:::
+:::
+
+
+At first glance, this may not look any different from our 2^nd^ or 3^rd^ order polynomial model. But the resulting residuals suggest that the loess smooth did a better job in removing any decadal patterns in our batch of CO2 values.
+
+
+::: {.cell hash='example_cache/html/unnamed-chunk-13_6b9fcf514f270f4e9f3bef4be83df413'}
+
+```{.r .cell-code}
+lo <- loess(CO2 ~ Year, dat2, span=1/4)
+dat2$res3 <- lo$residuals
+
+ggplot(dat2,  aes(x = Year, y = res3)) + geom_point() +
+           stat_smooth(method = "loess", se = FALSE, span = 0.2)
+```
+
+::: {.cell-output-display}
+![](example_files/figure-html/unnamed-chunk-13-1.png){width=576}
+:::
+:::
+
+
+## Exploring the seasonal component
+
+Let's now explore the cyclical pattern in the residuals. Note the residuals' y-axis values: they are three orders of magnitude less than the overall range of CO2 values. This indicates that the oscillating component of the data is relatively small (the CO2 values have a range of [313.33, 417.31] ppm whereas the residuals have a range of [-4.34, 4.34] ppm).
+
+Now, we may be tempted to fit a *smooth* to the residuals but that may not prove to be fruitful. Instead, let's see if the oscillation follows a 12 month cycle. We'll group all the residual values by month of the year. In other words, we will remove the year *tag* associated with each value and explore those values as a function of month alone. Each month's batch of values is distilled into a boxplot.
+
+
+::: {.cell hash='example_cache/html/unnamed-chunk-14_34a9ab0e1daca67c09ee77a017a7e443'}
+
+```{.r .cell-code}
+# Aggregate residuals by year
+dat2$Month <- month(dat2$Date, label=TRUE) 
+
+ggplot(dat2, aes(x = Month, y = res3)) + geom_boxplot()         
+```
+
+::: {.cell-output-display}
+![](example_files/figure-html/unnamed-chunk-14-1.png){width=576}
+:::
+:::
+
+
+It's clear from the plot that the oscillation follows a yearly cycle: a peak in the spring and a dip in the fall. This cycle is explained in part by the increased land mass in the northern hemisphere relative to the southern hemisphere. Because plants (and by extension photosynthesis) go dormant during the winter months in the northern hemisphere, atmospheric CO2 is no longer being photosynthesized; this despite the southern hemisphere's photosynthesis peak during the October-March period (a result of the southern hemisphere's smaller land mass). Other factors such as increased CO2 emissions during the winter months  may also contribute to the oscillatory nature of atmospheric CO2 concentrations.
+
+Thus far, we have uncovered three patterns of interest: an overall trend, a "hump" around the 1980-1995 time period, and a cyclical component. Note that to effectively explore the cyclical pattern we had to **de-trend** (or **smooth**) the data. Next, we should *smooth* the seasonal component of the data to see if another pattern emerges. We may, for example, smooth the data by subtracting the monthly median from each residual batch leaving us with the next batch of residual values to explore:
+
+
+::: {.cell hash='example_cache/html/unnamed-chunk-15_eda46418950b965904f820dc8f0de69d'}
+
+```{.r .cell-code}
+d3 <- dat2 %>% 
+  group_by(month(Date)) %>%
+  mutate(Resid.lev = res3 - median(res3)) 
+
+ggplot(d3, aes(x = Month, y = Resid.lev)) + geom_boxplot()     
+```
+
+::: {.cell-output-display}
+![](example_files/figure-html/unnamed-chunk-15-1.png){width=768}
+:::
+:::
+
+
+All the boxplots are lined up along their median values. We are now exploring the data after having accounted for both overall trend and seasonality. What can be gleaned from this dataset? We may want to explore the skewed nature of the residuals in March, or the narrower range of CO2 values for the fall months, for example. 
+
+We may also be interested in seeing if any trend in CO2 concentrations exists within each month. Let's explore this by fitting a first order polynomial to the monthly data.
+
+
+::: {.cell hash='example_cache/html/unnamed-chunk-16_8696b7ec21d0b5f901e1b7fc5a1fa3a0'}
+
+```{.r .cell-code}
+ggplot(d3, aes(x = Year, y = Resid.lev)) + geom_point(col = "grey", cex = .6)  + 
+     stat_smooth(se = FALSE, method = "lm") + facet_wrap(~ Month, nrow = 1) +
+     theme(axis.text.x=element_text(angle=45, hjust = 1.5, vjust = 1.6,  size = 7))
+```
+
+::: {.cell-output-display}
+![](example_files/figure-html/unnamed-chunk-16-1.png){width=768}
+:::
+:::
+
+
+Note the overall increase in CO2 for the winter and spring months followed by a similar (but not as persistent) decrease in CO2 loss in the summer and fall months. This suggests either an increase in cyclical amplitude with the deviation from a central value being possibly greater for the winter/spring months. It could also suggest a gradual offset in the periodicity over the 60 year period. More specifically, a leftward shift in periodicity.
+
+Let's fit a loess to see if the trends within each month are monotonic. We'll adopt a robust loess (`family = "symmetric"`) to control for possible outliers in the data.
+
+
+::: {.cell hash='example_cache/html/unnamed-chunk-17_877c39917518a898b62a09cb5cdb7629'}
+
+```{.r .cell-code}
+ggplot(d3, aes(x = Year, y = Resid.lev)) + geom_point(col = "grey", cex = .6)  + 
+     stat_smooth(method = "loess", se = FALSE, 
+                 method.args = list(family = "symmetric")) + 
+     facet_wrap(~ Month, nrow = 1) +
+     theme(axis.text.x=element_text(angle=45, hjust = 1.5, vjust = 1.6,  size = 7))
+```
+
+::: {.cell-output-display}
+![](example_files/figure-html/unnamed-chunk-17-1.png){width=768}
+:::
+:::
+
+
+The trends are clearly not monotonic suggesting that the processes driving the cyclical concentrations of CO2 are not systematic. However, despite the curved nature of the fits, the overall trends observed with the straight lines still remain.
+
+## Exploring what remains
+
+Now that we've fitted the monthly medians to the seasonal component of the data, let's plot the residuals across the years. We'll make use of the robust version of the loess given that the data may be quite noisy.
+
+
+::: {.cell hash='example_cache/html/unnamed-chunk-18_4e96c68db86bc64e05e956e9304ef83e'}
+
+```{.r .cell-code}
+ggplot(d3, aes(x = Year, y = Resid.lev)) +  geom_line(col = "grey") + 
+     stat_smooth(method = "loess", se = FALSE, span = 0.1 ,
+                 method.args = list(family = "symmetric", degree = 2))    
+```
+
+::: {.cell-output-display}
+![](example_files/figure-html/unnamed-chunk-18-1.png){width=768}
+:::
+:::
+
+
+Recall that at this points, we have accounted for the overall trend (which has by far the greatest effect with a range of over 90 ppm) and the seasonal component (which has a smaller effect with a range of about 6 ppm). The 3^rd^ set of residuals account for about 2 ppm of the total variability in the dataset. Note that we can still make out the seasonal fluctuations if we look closely at the data. But we can also make out broader patterns of interest not consistent with random noise and not explained by seasonal variability. Some of the patterns could be explained by the southern oscillation index (SOI), for example,  which is a measure of the difference in atmospheric pressure across the pacific basin. When the index is positive the resulting winds drive upwelling currents that release additional CO2 into the atmosphere. 
+
+We'll create a simple overlay the SOI data on top of our last plot, but note that a scatter plot of the two variables with a lagging component would be best suited for this analysis.
+
+
+::: {.cell hash='example_cache/html/unnamed-chunk-19_af66d480a9a97d8f1f14f998181dd648'}
+
+```{.r .cell-code}
+library(tidyr)
+
+soi <- read.csv("http://mgimond.github.io/ES218/Data/Southern_oscillation.csv")
+
+soi <- gather(soi, key = Month, value = Index, -1)
+Date <- ymd(paste(soi$YEAR, soi$Month, 15) )
+soi$Yr <- decimal_date(Date)
+
+ggplot(d3, aes(x = Year, y = Resid.lev)) +  
+             geom_line(col = "blue", alpha = 0.2) + 
+             stat_smooth(method = "loess", se = FALSE, span = 0.1, lwd = 0.5,
+                         method.args = list(family = "symmetric", degree = 2)) +
+             geom_line(dat = soi, aes(Yr, Index/2, colour = "red"), alpha = 0.3) +
+             stat_smooth(dat = soi, aes(Yr, Index/2, colour = "red"), lwd = 0.5,
+                         method = "loess", se = FALSE, span = 0.2,
+                         method.args = list(family = "symmetric", degree = 2)) +
+             guides(color = FALSE)
+```
+
+::: {.cell-output-display}
+![](example_files/figure-html/unnamed-chunk-19-1.png){width=768}
+:::
+:::
+
+
+The blue lines represent the residuals (with the bold blue line showing the loess fit) and the red lines represent the SOI trend (with the bold line showing its loess fit). Note that we halved the SOI values to match the range along the y-scale; the intent is to facilitate the comparison of SOI trend to CO2 trend. There appears to be matching maxima (notably around 1975 and 1988) but we also note some opposing trends such as in the period after 2010. When the index drops below 0, "El Nino" events become more prominent and upwelling flows are suppressed. Though this may prevent CO2 trapped in the deeper Pacific oceans from reaching the surface, "El Nino" events can also have a measurable impact on land climate thus influencing the rate of vegetation growth and, by extension, CO2 uptake. El nino is also believed to increase the rate of forest and brush fires around the world which may also contribute to atmospheric CO2 concentrations.  
+
+It's important to realize that this is one of many ways in which we could have proceeded with the analysis. For example, we could have started off the analysis by removing the seasonal component from the data, then analyzed the long term trend. 
+
diff --git a/example_files/execute-results/html.json b/example_files/execute-results/html.json
new file mode 100644
index 00000000..cdeade25
--- /dev/null
+++ b/example_files/execute-results/html.json
@@ -0,0 +1,18 @@
+{
+  "hash": "e0ea20c732c9a2b53fab033d2b69fa12",
+  "result": {
+    "markdown": "# The third R of EDA: Residuals\n\n\n::: {.cell}\n\n:::\n\n::: {.cell hash='example_cache/html/unnamed-chunk-2_fd93e819e7915059036d51d67790e90a'}\n::: {.cell-output-display}\n`````{=html}\n<table class=\"table\" style=\"width: auto !important; \">\n <thead>\n  <tr>\n   <th style=\"text-align:left;color: #555555 !important;background-color: #dddddd !important;text-align: center;border: 1px solid white !important;\n             font-family: 'Source Code Pro', 'Open Sans';\n             padding:1px !important;\n             padding-left:4px !important;\n             padding-right:4px !important;\n             font-size: 0.8em;\n             border-radius: 5px;\"> dplyr </th>\n   <th style=\"text-align:left;color: #555555 !important;background-color: #dddddd !important;text-align: center;border: 1px solid white !important;\n             font-family: 'Source Code Pro', 'Open Sans';\n             padding:1px !important;\n             padding-left:4px !important;\n             padding-right:4px !important;\n             font-size: 0.8em;\n             border-radius: 5px;\"> ggplot2 </th>\n   <th style=\"text-align:left;color: #555555 !important;background-color: #dddddd !important;text-align: center;border: 1px solid white !important;\n             font-family: 'Source Code Pro', 'Open Sans';\n             padding:1px !important;\n             padding-left:4px !important;\n             padding-right:4px !important;\n             font-size: 0.8em;\n             border-radius: 5px;\"> tidyr </th>\n   <th style=\"text-align:left;color: #555555 !important;background-color: #dddddd !important;text-align: center;border: 1px solid white !important;\n             font-family: 'Source Code Pro', 'Open Sans';\n             padding:1px !important;\n             padding-left:4px !important;\n             padding-right:4px !important;\n             font-size: 0.8em;\n             border-radius: 5px;\"> lubridate </th>\n  </tr>\n </thead>\n<tbody>\n  <tr>\n   <td style=\"text-align:left;color: darkred !important;background-color: #FAE8E8 !important;text-align: center;border: 1px solid white;\n             font-family: 'Open Sans', Arial;\n             padding:1px !important;\n             padding-left:4px !important;\n             padding-right:4px !important;\n             font-size: 0.8em;\n             border-radius: 5px;\"> 1.1.4 </td>\n   <td style=\"text-align:left;color: darkred !important;background-color: #FAE8E8 !important;text-align: center;border: 1px solid white;\n             font-family: 'Open Sans', Arial;\n             padding:1px !important;\n             padding-left:4px !important;\n             padding-right:4px !important;\n             font-size: 0.8em;\n             border-radius: 5px;\"> 3.4.4 </td>\n   <td style=\"text-align:left;color: darkred !important;background-color: #FAE8E8 !important;text-align: center;border: 1px solid white;\n             font-family: 'Open Sans', Arial;\n             padding:1px !important;\n             padding-left:4px !important;\n             padding-right:4px !important;\n             font-size: 0.8em;\n             border-radius: 5px;\"> 1.3.0 </td>\n   <td style=\"text-align:left;color: darkred !important;background-color: #FAE8E8 !important;text-align: center;border: 1px solid white;\n             font-family: 'Open Sans', Arial;\n             padding:1px !important;\n             padding-left:4px !important;\n             padding-right:4px !important;\n             font-size: 0.8em;\n             border-radius: 5px;\"> 1.9.3 </td>\n  </tr>\n</tbody>\n</table>\n\n`````\n:::\n:::\n\n\n## An EDA example: CO2 analysis\n\nLet's start off by plotting the atmospheric CO2 concentrations (in ppm). The original dataset was pulled from [NOAA's website](http://www.esrl.noaa.gov/gmd/ccgg/trends/).\n\n\n::: {.cell hash='example_cache/html/unnamed-chunk-3_8ff541de663eb986fb7aa84d57d9e4a8'}\n\n```{.r .cell-code}\nlibrary(lubridate)\nlibrary(dplyr)\nlibrary(ggplot2)\n\ndat <- read.csv(\"http://mgimond.github.io/ES218/Data/co2_2020.csv\")\n\ndat2 <- dat %>%\n      mutate( Date = ymd(paste(Year,\"/\", Month,\"/15\")),\n              Year = decimal_date(Date),\n              CO2 = Average)    %>%\n     select(Year, Date, CO2)\n```\n:::\n\n\n\nNext, let's look at the data:\n      \n\n::: {.cell hash='example_cache/html/unnamed-chunk-4_18cffa525f59a03421c6392b95970664'}\n\n```{.r .cell-code}\nggplot(dat2, aes(x = Year , y = CO2)) + geom_point(size = 0.6)\n```\n\n::: {.cell-output-display}\n![](example_files/figure-html/unnamed-chunk-4-1.png){width=576}\n:::\n:::\n\n\nLet's see if plotting the data using lines instead of points helps identify any underlying patterns.\n      \n\n::: {.cell hash='example_cache/html/unnamed-chunk-5_55c2b633da90110e2e7c31caf6a481d5'}\n\n```{.r .cell-code}\nggplot(dat2, aes(x = Year , y = CO2)) + geom_line()\n```\n\n::: {.cell-output-display}\n![](example_files/figure-html/unnamed-chunk-5-1.png){width=576}\n:::\n:::\n\n\nWe note two patterns of interest: an overall upward trend, and a cyclical trend. We will first *smooth* out the overall trend by finding a model that best explains the overall variability, then we'll remove that trend from the data to explore the cyclical component.\n\n## Exploring the overall trend\n\nWe can attempt to model the overall trend by fitting a straight line to the data using a standard regression analysis procedure. The fitted line is displayed in red in the following plot. We will use `ggplot2`'s built-in smoothing function to plot the regression line.\n\n\n\n::: {.cell hash='example_cache/html/unnamed-chunk-6_2161d4afd86c476bae6977716b12b4a3'}\n\n```{.r .cell-code}\nggplot(dat2, aes(x = Year, y = CO2)) + geom_line() + \n             stat_smooth(method = \"lm\", se = FALSE, col = \"red\") \n```\n\n::: {.cell-output-display}\n![](example_files/figure-html/unnamed-chunk-6-1.png){width=576}\n:::\n:::\n\n\nNext, we subtract the modeled line from the CO2 data and plot the result--this gives us the **residuals** and can be thought of as representing what the linear model does *not* explain about the data. But first, we will need to create a regression model object since it will provide us with the residuals needed to generate the residual plot.\n\n\n::: {.cell hash='example_cache/html/unnamed-chunk-7_a35d55aa4f50ac91d7535d1ce418686b'}\n\n```{.r .cell-code}\n# Generate a model object from the regression fit\nM1 <- lm(CO2 ~ Year, dat2) \ndat2$res1 <- M1$residuals\n\n# Plot the residuals from the model object\nggplot(dat2, aes(x = Year, y = res1)) + geom_point()\n```\n\n::: {.cell-output-display}\n![](example_files/figure-html/unnamed-chunk-7-1.png){width=576}\n:::\n:::\n\n\nAn overall trend is still present, despite having attempted to control for it. This suggests that our simple linear model does not do a good job in *capturing* the shape of the overall trend. \n\nIt appears that the overall trend is slightly convex and has a small peak around the 1990's. We can first try to fit the trend using a 2^nd^ order polynomial of the form:\n\n$$\nCO2_{trend} = a + b(Year) + c(Year^2) \n$$\n\nThe fitted line looks like this: \n\n\n::: {.cell hash='example_cache/html/unnamed-chunk-8_3a3326053a34672abde81f4ec6b7a77f'}\n\n```{.r .cell-code}\nggplot(dat2, aes(x = Year, y = CO2)) + geom_line() + \n           stat_smooth(method = \"lm\", formula = y ~ x + I(x^2), \n                       se = FALSE, col = \"red\")\n```\n\n::: {.cell-output-display}\n![](example_files/figure-html/unnamed-chunk-8-1.png){width=576}\n:::\n:::\n\n\nNow, let's look at the residuals. In helping discern any pattern in our residuals, we add a loess smooth to the plot.\n\n\n::: {.cell hash='example_cache/html/unnamed-chunk-9_737f2732efa31eb98a846acf9540f805'}\n\n```{.r .cell-code}\nM2 <- lm(CO2 ~ Year + I(Year^2), dat2)\ndat2$res2 <- M2$residuals\n\nggplot(dat2 , aes(x = Year, y = res2)) + geom_point() +\n           stat_smooth(method = \"loess\", se = FALSE)\n```\n\n::: {.cell-output-display}\n![](example_files/figure-html/unnamed-chunk-9-1.png){width=576}\n:::\n:::\n\n\nThis is an improvement over the simple line model (which was a 1^st^ order polynomial fit).  So what we have learned so far is that the overall trend is not perfectly linear but instead follows a parabolic like trajectory with a small hump halfway across the time span. However, we can still make out a \"W\" shaped trend in the residuals which can hamper our analysis of the oscillatory patterns in the data. We can try a 3^rd^ order polynomial and see if it can capture the 1980 to 1995 hump.\n\n\n\n::: {.cell hash='example_cache/html/unnamed-chunk-10_b5a0f63a1435b46624cbc8528f9f1a7e'}\n\n```{.r .cell-code}\nggplot(dat2, aes(x = Year, y = CO2)) + geom_line() + \n           stat_smooth(method = \"lm\", formula = y ~ x + I(x^2) + I(x^3), \n                       se = FALSE, col = \"red\")\n```\n\n::: {.cell-output-display}\n![](example_files/figure-html/unnamed-chunk-10-1.png){width=576}\n:::\n:::\n\n\nNow, let's look at the residuals. \n\n\n::: {.cell hash='example_cache/html/unnamed-chunk-11_b72a0542c1f88351317336ac149a8bfd'}\n\n```{.r .cell-code}\nM2 <- lm(CO2 ~ Year + I(Year^2) + I(Year^3), dat2)\ndat2$res2 <- M2$residuals\n\nggplot(dat2 , aes(x = Year, y = res2)) + geom_point() +\n           stat_smooth(method = \"loess\", se = FALSE)\n```\n\n::: {.cell-output-display}\n![](example_files/figure-html/unnamed-chunk-11-1.png){width=576}\n:::\n:::\n\n\nThis does not seem to be an improvement over the 2^nd^ order polynomial. At this point, we could play with different order polynomials in an attempt to smooth the trend or, we may opt for a non-parametric fit. When the goal is to *peel* off one pattern to explore any underlying pattern we should not limit ourselves to parametric fits (which impose a mathematical model on our data) and instead explore non-parametric smoothing techniques that do not impose any structure on the data whatsoever. An example of such a smoothing technique is the loess fit which is shown in the following figure.\n\n\n::: {.cell hash='example_cache/html/unnamed-chunk-12_83a89c9aa33856f9a0d1116388596410'}\n\n```{.r .cell-code}\nggplot(dat2, aes(x = Year, y = CO2)) + geom_line() + \n             stat_smooth(method = \"loess\", span = 1/4, se = FALSE, col = \"red\")\n```\n\n::: {.cell-output-display}\n![](example_files/figure-html/unnamed-chunk-12-1.png){width=576}\n:::\n:::\n\n\nAt first glance, this may not look any different from our 2^nd^ or 3^rd^ order polynomial model. But the resulting residuals suggest that the loess smooth did a better job in removing any decadal patterns in our batch of CO2 values.\n\n\n::: {.cell hash='example_cache/html/unnamed-chunk-13_6b9fcf514f270f4e9f3bef4be83df413'}\n\n```{.r .cell-code}\nlo <- loess(CO2 ~ Year, dat2, span=1/4)\ndat2$res3 <- lo$residuals\n\nggplot(dat2,  aes(x = Year, y = res3)) + geom_point() +\n           stat_smooth(method = \"loess\", se = FALSE, span = 0.2)\n```\n\n::: {.cell-output-display}\n![](example_files/figure-html/unnamed-chunk-13-1.png){width=576}\n:::\n:::\n\n\n## Exploring the seasonal component\n\nLet's now explore the cyclical pattern in the residuals. Note the y-axis values: they are three orders of magnitude less than the overall range of CO2 values. This indicates that the oscillating nature along the overall trend is relatively small (the CO2 values have a range of [313.33, 417.31] ppm whereas the residuals have a range of [-4.34, 4.34] ppm).\n\nNow, we may be tempted to fit a *smooth* (as Tukey would say) to the residuals but that may not prove to be fruitful. Instead let's see if the oscillation follows a 12 month cycle. We'll group all the residual values by month of the year. In other words, we will remove the year *tag* associated with each value and explore those values as a function of month alone. Each month's batch of values is distilled into a boxplot.\n\n\n::: {.cell hash='example_cache/html/unnamed-chunk-14_34a9ab0e1daca67c09ee77a017a7e443'}\n\n```{.r .cell-code}\n# Aggregate residuals by year\ndat2$Month <- month(dat2$Date, label=TRUE) \n\nggplot(dat2, aes(x = Month, y = res3)) + geom_boxplot()         \n```\n\n::: {.cell-output-display}\n![](example_files/figure-html/unnamed-chunk-14-1.png){width=576}\n:::\n:::\n\n\nIt's clear from the plot that the oscillation follows a yearly cycle: a peak in the spring and a dip in the fall. This cycle is explained in part by the increased land mass in the northern hemisphere relative to the southern hemisphere. Because plants (and by extension photosynthesis) go dormant during the winter months in the northern hemisphere, atmospheric CO2 is no longer being photosynthesized; this despite the southern hemisphere's photosynthesis peak during the October-March period (a result of the southern hemisphere's smaller land mass). Other factors such as increased CO2 emissions during the winter months  may also contribute to the oscillatory nature of atmospheric CO2 concentrations.\n\nThus far, we have uncovered three patterns of interest: an overall trend, a \"hump\" around the 1980-1995 time period, and a cyclical component. Note that to effectively explore the cyclical pattern we had to **de-trend** (or **smooth**) the data. Next, we should *smooth* the seasonal component of the data to see if another pattern emerges. We may, for example, smooth the data by subtracting the monthly median from each residual batch leaving us with the next batch of residual values to explore:\n\n\n::: {.cell hash='example_cache/html/unnamed-chunk-15_eda46418950b965904f820dc8f0de69d'}\n\n```{.r .cell-code}\nd3 <- dat2 %>% \n  group_by(month(Date)) %>%\n  mutate(Resid.lev = res3 - median(res3)) \n\nggplot(d3, aes(x = Month, y = Resid.lev)) + geom_boxplot()     \n```\n\n::: {.cell-output-display}\n![](example_files/figure-html/unnamed-chunk-15-1.png){width=768}\n:::\n:::\n\n\nAll the boxplots are lined up along their median values. We are now exploring the data after having accounted for both overall trend and seasonality. What can be gleaned from this dataset? We may want to explore the skewed nature of the residuals in March, or the narrower range of CO2 values for the fall months, for example. \n\nWe may also be interested in seeing if any trend in CO2 concentrations exists within each month. Let's explore this by fitting a first order polynomial to the monthly data.\n\n\n::: {.cell hash='example_cache/html/unnamed-chunk-16_8696b7ec21d0b5f901e1b7fc5a1fa3a0'}\n\n```{.r .cell-code}\nggplot(d3, aes(x = Year, y = Resid.lev)) + geom_point(col = \"grey\", cex = .6)  + \n     stat_smooth(se = FALSE, method = \"lm\") + facet_wrap(~ Month, nrow = 1) +\n     theme(axis.text.x=element_text(angle=45, hjust = 1.5, vjust = 1.6,  size = 7))\n```\n\n::: {.cell-output-display}\n![](example_files/figure-html/unnamed-chunk-16-1.png){width=768}\n:::\n:::\n\n\nNote the overall increase in CO2 for the winter and spring months followed by a similar (but not as persistent) decrease in CO2 loss in the summer and fall months. This suggests either an increase in cyclical amplitude with the deviation from a central value being possibly greater for the winter/spring months. It could also suggest a gradual offset in the periodicity over the 60 year period. More specifically, a leftward shift in periodicity.\n\nLet's fit a loess to see if the trends within each month are monotonic. We'll adopt a robust loess (`family = \"symmetric\"`) to control for possible outliers in the data.\n\n\n::: {.cell hash='example_cache/html/unnamed-chunk-17_877c39917518a898b62a09cb5cdb7629'}\n\n```{.r .cell-code}\nggplot(d3, aes(x = Year, y = Resid.lev)) + geom_point(col = \"grey\", cex = .6)  + \n     stat_smooth(method = \"loess\", se = FALSE, \n                 method.args = list(family = \"symmetric\")) + \n     facet_wrap(~ Month, nrow = 1) +\n     theme(axis.text.x=element_text(angle=45, hjust = 1.5, vjust = 1.6,  size = 7))\n```\n\n::: {.cell-output-display}\n![](example_files/figure-html/unnamed-chunk-17-1.png){width=768}\n:::\n:::\n\n\nThe trends are clearly not monotonic suggesting that the processes driving the cyclical concentrations of CO2 are not systematic. However, despite the curved nature of the fits, the overall trends observed with the straight lines still remain.\n\n## Exploring what remains\n\nNow that we've fitted the monthly medians to the seasonal component of the data, let's plot the residuals across the years. We'll make use of the robust version of the loess given that the data may be quite noisy.\n\n\n::: {.cell hash='example_cache/html/unnamed-chunk-18_4e96c68db86bc64e05e956e9304ef83e'}\n\n```{.r .cell-code}\nggplot(d3, aes(x = Year, y = Resid.lev)) +  geom_line(col = \"grey\") + \n     stat_smooth(method = \"loess\", se = FALSE, span = 0.1 ,\n                 method.args = list(family = \"symmetric\", degree = 2))    \n```\n\n::: {.cell-output-display}\n![](example_files/figure-html/unnamed-chunk-18-1.png){width=768}\n:::\n:::\n\n\nRecall that at this points, we have accounted for the overall trend (which has by far the greatest effect with a range of over 90 ppm) and the seasonal component (which has a smaller effect with a range of about 6 ppm). The 3^rd^ set of residuals account for about 2 ppm of the total variability in the dataset. Note that we can still make out the seasonal fluctuations if we look closely at the data. But we can also make out broader patterns of interest not consistent with random noise and not explained by seasonal variability. Some of the patterns could be explained by the southern oscillation index (SOI), for example,  which is a measure of the difference in atmospheric pressure across the pacific basin. When the index is positive the resulting winds drive upwelling currents that release additional CO2 into the atmosphere. \n\nWe'll create a simple overlay the SOI data on top of our last plot, but note that a scatter plot of the two variables with a lagging component would be best suited for this analysis.\n\n\n::: {.cell hash='example_cache/html/unnamed-chunk-19_af66d480a9a97d8f1f14f998181dd648'}\n\n```{.r .cell-code}\nlibrary(tidyr)\n\nsoi <- read.csv(\"http://mgimond.github.io/ES218/Data/Southern_oscillation.csv\")\n\nsoi <- gather(soi, key = Month, value = Index, -1)\nDate <- ymd(paste(soi$YEAR, soi$Month, 15) )\nsoi$Yr <- decimal_date(Date)\n\nggplot(d3, aes(x = Year, y = Resid.lev)) +  \n             geom_line(col = \"blue\", alpha = 0.2) + \n             stat_smooth(method = \"loess\", se = FALSE, span = 0.1, lwd = 0.5,\n                         method.args = list(family = \"symmetric\", degree = 2)) +\n             geom_line(dat = soi, aes(Yr, Index/2, colour = \"red\"), alpha = 0.3) +\n             stat_smooth(dat = soi, aes(Yr, Index/2, colour = \"red\"), lwd = 0.5,\n                         method = \"loess\", se = FALSE, span = 0.2,\n                         method.args = list(family = \"symmetric\", degree = 2)) +\n             guides(color = FALSE)\n```\n\n::: {.cell-output-display}\n![](example_files/figure-html/unnamed-chunk-19-1.png){width=768}\n:::\n:::\n\n\nThe blue lines represent the residuals (with the bold blue line showing the loess fit) and the red lines represent the SOI trend (with the bold line showing its loess fit). Note that we halved the SOI values to match the range along the y-scale; the intent is to facilitate the comparison of SOI trend to CO2 trend. There appears to be matching maxima (notably around 1975 and 1988) but we also note some opposing trends such as in the period after 2010. When the index drops below 0, \"El Nino\" events become more prominent and upwelling flows are suppressed. Though this may prevent CO2 trapped in the deeper Pacific oceans from reaching the surface, \"El Nino\" events can also have a measurable impact on land climate thus influencing the rate of vegetation growth and, by extension, CO2 uptake. El nino is also believed to increase the rate of forest and brush fires around the world which may also contribute to atmospheric CO2 concentrations.  \n\nIt's important to realize that this is one of many ways in which we could have proceeded with the analysis. For example, we could have started off the analysis by removing the seasonal component from the data, then analyzed the long term trend. This is the approach taken by William Cleveland on page 165 of the course's text.\n",
+    "supporting": [],
+    "filters": [
+      "rmarkdown/pagebreak.lua"
+    ],
+    "includes": {
+      "include-in-header": [
+        "<script src=\"site_libs/kePrint-0.0.1/kePrint.js\"></script>\r\n<link href=\"site_libs/lightable-0.0.1/lightable.css\" rel=\"stylesheet\" />\r\n"
+      ]
+    },
+    "engineDependencies": {},
+    "preserve": {},
+    "postProcess": true
+  }
+}
\ No newline at end of file
diff --git a/index.qmd b/index.qmd
index dbd1588f..364776dc 100644
--- a/index.qmd
+++ b/index.qmd
@@ -5,7 +5,7 @@ This book is a compilation of lecture notes used in an *Exploratory Data Analysi
 
 The base R plotting environment and the `ggplot2` ecosystem are used throughout this book. While a chapter is dedicated to the `lattice` plotting package, its functions are not used outside of that chapter given that `ggplot2` offers many of `lattice`'s functionality.
 
-While great effort is made to adopt a consist plotting environment throughout this book (this being `ggplot2`, for the most part), a few topics (including the q-q plot and the median polish) will benefit from custom plotting functions available in the `tukeyedar` package. The [package](https://github.com/mgimond/tukeyedar/) can be downloaded from GitHub via the command:
+While great effort is made to adopt a consistent plotting environment throughout this book (this being `ggplot2`, for the most part), a few topics (including the q-q plot and the median polish) will benefit from custom plotting functions available in the `tukeyedar` package. The [package](https://github.com/mgimond/tukeyedar/) can be downloaded from GitHub via the command:
 
 ```{r eval = FALSE}
 devtools::install_github("mgimond/tukeyedar")
diff --git a/intro.qmd b/intro.qmd
index 6f7af06c..fda20cea 100644
--- a/intro.qmd
+++ b/intro.qmd
@@ -126,7 +126,7 @@ Many journals now require or *strongly encourage* authors to *"make materials, d
 
 Another source of error in the write-up of a  report or publication is the linking of tables, figures and statistical summaries to the write-up. Typically, one saves statistical plots as image files then loads the image into the document. However, the figures may have gone through many different iterations resulting in many different versions of the image file in a working folder. Add to this many other figures, data table files and statistical results from various pieces of software, one quickly realizes the potential for embedding the wrong image files in the document or embedding the wrong statistical summaries in the text. Furthermore, the researcher is then required to properly archive and document the provenance of each figure, data table or statistical summary resulting in a complex structure of files and directories in the project folder thus increasing the odds of an irreproducible analysis.
 
-Confining all of the analysis to a scripting environment such as R can help, but this still does not alleviate the possibility of loading the wrong figure into the document, or forgetting to update a statistical summary in the text when the original data file was revised. A solution to this potential pitfall is to embed the actual analysis and graphic generation process into the document--such environments are called dynamic documents. In this course, we will use the [R Markdown authoring tool](http://rmarkdown.rstudio.com/) which embeds R code into the document. An example of an R Markdown document is this course website which was entirely generated in RMarkdown!. You can view the  R Markdown files [on this author's GitHub repository](https://github.com/mgimond/ES218).
+Confining all of the analysis to a scripting environment such as R can help, but this still does not alleviate the possibility of loading the wrong figure into the document, or forgetting to update a statistical summary in the text when the original data file was revised. A solution to this potential pitfall is to embed the actual analysis and graphic generation process into the document--such environments are called dynamic documents. In this course, we will use the [R Markdown authoring tool](http://rmarkdown.rstudio.com/) which embeds R code into the document. An example of an R Markdown document is this course website which was entirely generated in RMarkdown! You can view the  R Markdown files [on this author's GitHub repository](https://github.com/mgimond/ES218).
 
 
 [^1]: Tukey, John W. *Exploratory Data Analysis*. 1977. Addison-Wesley.
diff --git a/resistant_lines.qmd b/resistant_lines.qmd
index a26b0fb6..b45ef4e9 100644
--- a/resistant_lines.qmd
+++ b/resistant_lines.qmd
@@ -50,7 +50,7 @@ coef(M)
 
 The slope dropped from 2 to 0.86 because of a single point!
 
-If our goal is to explore what the *bulk* of the data has to say about the phenomena being investigated, we certainly don't want a few "maverick" values to hijack the analysis. We therefore need a set of fitting tools that minimize the influence of outliers. There are many options out there; most notable are the **resistant line** (that builds off of Tukey's **3-point summary line**) and the **bisquare** robust estimation method.
+If our goal is to explore what the *bulk* of the data has to say about the phenomena being investigated, we certainly don't want a few "maverick" values to hijack the analysis. We therefore need a set of fitting tools that minimize the influence of outliers. There are many options out there; most notable are  **Tukey's resistant line** and the **bisquare** robust estimation method.
 
 &nbsp;&nbsp;
 
@@ -59,9 +59,9 @@ If our goal is to explore what the *bulk* of the data has to say about the pheno
 
 ### Tukey's resistant line
 
-The idea is simple in concept. It involves dividing the dataset into three approximately equal groups of values and summarizing these groups by computing their respective medians. Two half-slopes are then used to join the three points. Note that points that share the same x value are lumped into the same batch which can lead to unequal group sizes. The motivation behind this plot is to use the three-point summary to provide a robust assessment of the type of relationship between both variables. Such plots are often used to help guide re-expression of the variables x or y or both for the sole purpose of straightening the x-y relationship.
+The idea is simple in concept. It involves dividing the dataset into three approximately equal groups of values and summarizing these groups by computing their respective medians. The batches’ median values are then used to compute the slope and intercept. The motivation behind this plot is to use the three-point summary to provide a robust assessment of the type of relationship between both variables. 
 
-Let's look at an example using the last (modified) dataset. First, we divide the plot into three approximately equal batches.
+Let's look at an example of its implementation. We'll continue with the example used in the previous section. First, we divide the plot into three approximately equal batches.
 
 
 ```{r fig.height=3, fig.width=3,  small.mar=TRUE,echo=FALSE}
@@ -71,15 +71,12 @@ Mr <- eda_rline(df, x, y, iter=1)
 plot(Mr, model=F, pt3=F, fit = F)
 ```
 
-
-Next, compute the median x and y values within each section. 
+Next, we compute the median x and y values within each section. 
 
 ```{r fig.height=3, fig.width=3,  small.mar=TRUE,echo=FALSE}
 plot(Mr, model=F, pt3=T, fit = F)
 ```
 
-Note that we do not include the same point in any two median calculations. This implies that for the mid-third of the data, we do not include the point straddling the left boundary line when computing its median value. Likewise with the right-third of the data. 
-
 The x median values are `r Mr$xmed` and the y median values are `r Mr$ymed` respectively.
 
 The two end medians are used to compute the slope as:
@@ -111,13 +108,13 @@ From this first iteration, we can extract the residuals. A line is then fitted t
 
 ```{r fig.height=3, fig.width=3,  small.mar=TRUE,echo=FALSE}
 df$r1 <- residuals(rl1)
-plot(eda_rline(df, x, r1), model=F, pt3=T, fit = F)
+plot(eda_rline(df, x, r1), model=F, pt3=T, fit = F, ylab = "Residuals")
 rl1r <- line(x,df$r1,iter=1)
 abline(rl1r,col="red")
 ```
 
-The initial model slope and intercept are `r round(rl1$coefficients[1],2)` and `r round(rl1$coefficients[2],2)` respectively and the residual's slope and intercept are `r round(rl1r$coefficients[1],2)` and `r round(rl1r$coefficients[2],2)` respectively.
-The residual slope is then added to the first computed slope and the process is again repeated thus generating the following tweaked slope and updated residuals:
+The initial model intercept and slope are `r round(rl1$coefficients[1],3)` and `r round(rl1$coefficients[2],3)` respectively and the residual's slope and intercept are `r round(rl1r$coefficients[1],3)` and `r round(rl1r$coefficients[2],3)` respectively.
+The residual slope is then added to the first computed slope giving us an updated slope of `r round(rl1$coefficients[2] + rl1r$coefficients[2] ,3)` and the process is again repeated thus generating the following tweaked slope and updated residuals:
 
 ```{r fig.height=3, fig.width=6,  small.mar=TRUE,echo=FALSE}
 OP <- par(mfrow = c(1, 2))
@@ -126,13 +123,13 @@ rl2 <- line(x,y,iter=2)
 abline(rl2, col="red")
 
 df$r2 <- residuals(rl2)
-plot(eda_rline(df, x, r2), model=F, pt3=T, fit = F)
+plot(eda_rline(df, x, r2), model=F, pt3=T, fit = F, ylab = "Residuals")
 rl2r <- line(x,df$r2,iter=1)
 abline(rl2r,col="red")
 par(OP)
 ```
 
-The updated slope is now `r round(rl2$coefficients[1],2)`.  The iteration continues until the slope residuals stabilize. Each successive iteration is shown in the following plot with the final line in dark red. One more iteration should be enough to stabilize the residuals. Giving use the following fitted line whose intercept and slope are `r round(line(x,y,iter=3)$coefficients[1],2)` and `r round(line(x,y,iter=3)$coefficients[2],2)` respectively
+The second set of residuals is added to the previously computed slope giving us an updated slope of `r round(rl2$coefficients[2] +  rl2r$coefficients[2],3)`. The iteration continues until the slope residuals stabilize. Each successive iteration is shown in the following plot with the final line in dark red. This took three iterations giving us the final fitted line whose intercept and slope are `r round(line(x,y,iter=3)$coefficients[1],2)` and `r round(line(x,y,iter=3)$coefficients[2],2)` respectively--very close to the actual relationship between the two variables.
 
 ```{r fig.height=3, fig.width=3,  small.mar=TRUE,echo=FALSE}
 df$r1 <- residuals(rl1)
@@ -142,13 +139,13 @@ abline(line(x,y,iter=2),col="#FF8888")
 abline(line(x,y,iter=3),col="#FF0000")
 ```
 
-Tukey's resistant line can be computed using the R built-in `line()` function. You specify the number of iteration via the `iter` argument. Usually, five iterations will be enough.
+Tukey's resistant line can be computed using the R built-in `line()` function. You specify the maximum number of iterations via the `iter` argument. Usually, five iterations will be enough. But, the iteration will stop before the maximum number of iterations is reached if the residual slope stabilizes before then.
 
 ```{r}
 r.fit <- line(x,y, iter = 5)
 ```
 
-You can the resistant to a base plot as follows:
+You can add the resistant line to a base plot as follows:
 
 ```{r fig.height=3, fig.width=3, echo = 2:3}
 OP <- par(mar = c(4,4,1,0))
@@ -160,6 +157,7 @@ par(OP)
 You can also add the line to a `ggplot`:
 
 ```{r fig.height=3, fig.width=3}
+library(ggplot2)
 ggplot(df, aes(x,y)) + geom_point() +
   geom_abline(intercept = r.fit$coefficients[1], 
               slope = r.fit$coefficients[2],
@@ -171,7 +169,7 @@ ggplot(df, aes(x,y)) + geom_point() +
 
 A bisquare fitting strategy makes use of a non-parametric model that assigns *weights* to observations based on their proximity to the center of the distribution. The closer an observation is to that center, the greater its weight. 
 
-The first step is to run a linear regression model on the data then to extract its residuals. Next, a weighting scheme is fit to the residuals such that the points associated with extreme residuals are assigned the smallest weight and the points associated with the smaller residual values are assigned the largest weight. The regression analysis is then re-run using those same weights thus minimizing the influence of the rogue points. This process is repeated several times until the residuals stabilize. The following figure shows three iterations of the bisquare function whereby the weights (shown as grey text next to each point) start off as `1` then are modified following the residuals derived from the most recent regression model.
+The first step is to run a linear regression model on the data. The next step is to extract its residuals. Next, a weighting scheme is fitted to the residuals such that the points associated with extreme residuals are assigned the smallest weight and the points associated with the smaller residual values are assigned the largest weight. The regression analysis is then re-run using those same weights thus minimizing the influence of the rogue points. This process is repeated several times until the residuals stabilize. The following figure shows three iterations of the bisquare function whereby the weights (shown as grey text next to each point) start off as `1` then are modified following the residuals derived from the most recent regression model.
 
 ```{r fig.height=2.5, fig.width=8, small.mar=TRUE, echo = FALSE}
 wt.bisquare <- function(u, c = 6) {
@@ -233,9 +231,9 @@ In the above example, the bisquare method does a great job in eliminating the ou
 
 #### Built-in implementation of the bisquare
 
-The `MASS` package has a robust linear modeling function called `rlm` that will implement a variation of the aforementioned bisquare estimation technique. Its results may differ slightly from those presented here, but the difference will be insignificant for the most part. 
+The `MASS` package has a robust linear modeling function called `rlm` that will implement a variation of the aforementioned bisquare estimation technique. Its results may differ slightly from those presented above, but the difference will be insignificant for the most part. 
 
-Note that if you make use of `dplyr` in a workflow, loading `MASS` after `dplyr` will mask dplyr's `select` function. This can be problematic. So you either want to load `MASS` before `dplyr`, or you can call the function via `MASS::rlm`. An example of its use follows. The default regression model is added as a dashed line for reference.
+Note that if you make use of `dplyr` in a workflow, loading `MASS` after `dplyr` will mask dplyr's `select` function. This can be problematic. So you either want to load `MASS` before `dplyr`, or you can call the function via `MASS::rlm`. An example of its use follows. 
 
 ```{r fig.height=2.5, fig.width=3,  small.mar=TRUE, fig.show='hold'}
 dat <- data.frame(x,y)
@@ -246,16 +244,17 @@ plot(y ~ x,dat, pch=16, col=rgb(0,0,0,0.2))
 # Add the robust line
 abline( Mr, col="red")
 
-# Add the default regression line
+# Add the default regression line for reference
 abline( M , col="grey50", lty=2)
 ```
 
+The default regression model is added as a dashed line for reference.
+
 The function `rlm` can also be called directly from within `ggplot`.
 
 ```{r fig.height=2.5, fig.width=3,  small.mar=TRUE, fig.show='hold'}
 library(ggplot2)
 ggplot(dat, aes(x=x, y=y)) + geom_point() + 
-       stat_smooth(method = "lm", se = FALSE, col = "grey50", lty = 2) +
        stat_smooth(method = MASS::rlm, se = FALSE, col = "red",
                    method.args = list(psi = "psi.bisquare"))
 ``` 
@@ -264,7 +263,7 @@ ggplot(dat, aes(x=x, y=y)) + geom_point() +
 
 ## Robust loess
 
-The bisquare estimation method can also be extended to the loess smoothing function by passing the `"symmetric"` option to the `family` parameter.
+The loess fit is also not immune to outliers. The bisquare estimation method can be extended to the loess smoothing function by passing the `"symmetric"` option to the `family` argument.
 
 
 ```{r fig.height=2.5, fig.width=3,  small.mar=TRUE,  fig.show='hold'}
@@ -284,19 +283,16 @@ lines(dat$x, predict(lo), col = "grey50", lty = 2)
 lines(dat$x, predict(lor), col = "red")
 ```
 
-The default loess fit is added as a dashed line for reference.
+The robust loess is in red, the default loess fit is added as a dashed line for reference.
 
 The robust `loess` can be implemented in `ggplot2`'s `stat_smooth` function as follows:
 
 ```{r fig.height=2.5, fig.width=3,  small.mar=TRUE, fig.show='hold'}
 library(ggplot2)
 ggplot(dat, aes(x=x, y=y)) + geom_point() + 
-          stat_smooth(method = "loess",                       # Regular loess 
-                      se=FALSE, col="grey50", lty=2)  +    
           stat_smooth(method = "loess", 
-                      method.args = list(family="symmetric"), # Robust loess
+                      method.args = list(family="symmetric"), 
                       se=FALSE,  col="red") 
-      
 ```
 
 The loess parameters are passed to the `loess` function via the call to `method.args`. The red curve is that of the robust loess and the grey dashed curve is that of the standard loess.
diff --git a/resistant_lines_cache/html/unnamed-chunk-10_41f87f969de578f7dc1e7a8ba640d2b3.RData b/resistant_lines_cache/html/unnamed-chunk-10_41f87f969de578f7dc1e7a8ba640d2b3.RData
new file mode 100644
index 00000000..dfbc38a4
Binary files /dev/null and b/resistant_lines_cache/html/unnamed-chunk-10_41f87f969de578f7dc1e7a8ba640d2b3.RData differ
diff --git a/resistant_lines_cache/html/unnamed-chunk-10_f05e07280210339891d5d97d22ddb550.rdb b/resistant_lines_cache/html/unnamed-chunk-10_41f87f969de578f7dc1e7a8ba640d2b3.rdb
similarity index 100%
rename from resistant_lines_cache/html/unnamed-chunk-10_f05e07280210339891d5d97d22ddb550.rdb
rename to resistant_lines_cache/html/unnamed-chunk-10_41f87f969de578f7dc1e7a8ba640d2b3.rdb
diff --git a/resistant_lines_cache/html/unnamed-chunk-10_f05e07280210339891d5d97d22ddb550.rdx b/resistant_lines_cache/html/unnamed-chunk-10_41f87f969de578f7dc1e7a8ba640d2b3.rdx
similarity index 100%
rename from resistant_lines_cache/html/unnamed-chunk-10_f05e07280210339891d5d97d22ddb550.rdx
rename to resistant_lines_cache/html/unnamed-chunk-10_41f87f969de578f7dc1e7a8ba640d2b3.rdx
diff --git a/resistant_lines_cache/html/unnamed-chunk-10_f05e07280210339891d5d97d22ddb550.RData b/resistant_lines_cache/html/unnamed-chunk-10_f05e07280210339891d5d97d22ddb550.RData
deleted file mode 100644
index 964e503a..00000000
Binary files a/resistant_lines_cache/html/unnamed-chunk-10_f05e07280210339891d5d97d22ddb550.RData and /dev/null differ
diff --git a/resistant_lines_cache/html/unnamed-chunk-11_7e538dd616e025fe6922c952147b2533.RData b/resistant_lines_cache/html/unnamed-chunk-11_7e538dd616e025fe6922c952147b2533.RData
new file mode 100644
index 00000000..963fdee4
Binary files /dev/null and b/resistant_lines_cache/html/unnamed-chunk-11_7e538dd616e025fe6922c952147b2533.RData differ
diff --git a/resistant_lines_cache/html/unnamed-chunk-11_b19fce4230a230eb339f211b686407eb.rdb b/resistant_lines_cache/html/unnamed-chunk-11_7e538dd616e025fe6922c952147b2533.rdb
similarity index 100%
rename from resistant_lines_cache/html/unnamed-chunk-11_b19fce4230a230eb339f211b686407eb.rdb
rename to resistant_lines_cache/html/unnamed-chunk-11_7e538dd616e025fe6922c952147b2533.rdb
diff --git a/resistant_lines_cache/html/unnamed-chunk-11_b19fce4230a230eb339f211b686407eb.rdx b/resistant_lines_cache/html/unnamed-chunk-11_7e538dd616e025fe6922c952147b2533.rdx
similarity index 100%
rename from resistant_lines_cache/html/unnamed-chunk-11_b19fce4230a230eb339f211b686407eb.rdx
rename to resistant_lines_cache/html/unnamed-chunk-11_7e538dd616e025fe6922c952147b2533.rdx
diff --git a/resistant_lines_cache/html/unnamed-chunk-11_b19fce4230a230eb339f211b686407eb.RData b/resistant_lines_cache/html/unnamed-chunk-11_b19fce4230a230eb339f211b686407eb.RData
deleted file mode 100644
index 0b2c9bf6..00000000
Binary files a/resistant_lines_cache/html/unnamed-chunk-11_b19fce4230a230eb339f211b686407eb.RData and /dev/null differ
diff --git a/resistant_lines_cache/html/unnamed-chunk-15_7494a1816f080868e5824f192c16c90d.RData b/resistant_lines_cache/html/unnamed-chunk-15_7494a1816f080868e5824f192c16c90d.RData
deleted file mode 100644
index 76e79273..00000000
Binary files a/resistant_lines_cache/html/unnamed-chunk-15_7494a1816f080868e5824f192c16c90d.RData and /dev/null differ
diff --git a/resistant_lines_cache/html/unnamed-chunk-15_c9fad914b2da127567b8036becee7f07.RData b/resistant_lines_cache/html/unnamed-chunk-15_c9fad914b2da127567b8036becee7f07.RData
new file mode 100644
index 00000000..78c31424
Binary files /dev/null and b/resistant_lines_cache/html/unnamed-chunk-15_c9fad914b2da127567b8036becee7f07.RData differ
diff --git a/resistant_lines_cache/html/unnamed-chunk-15_7494a1816f080868e5824f192c16c90d.rdb b/resistant_lines_cache/html/unnamed-chunk-15_c9fad914b2da127567b8036becee7f07.rdb
similarity index 100%
rename from resistant_lines_cache/html/unnamed-chunk-15_7494a1816f080868e5824f192c16c90d.rdb
rename to resistant_lines_cache/html/unnamed-chunk-15_c9fad914b2da127567b8036becee7f07.rdb
diff --git a/resistant_lines_cache/html/unnamed-chunk-15_7494a1816f080868e5824f192c16c90d.rdx b/resistant_lines_cache/html/unnamed-chunk-15_c9fad914b2da127567b8036becee7f07.rdx
similarity index 100%
rename from resistant_lines_cache/html/unnamed-chunk-15_7494a1816f080868e5824f192c16c90d.rdx
rename to resistant_lines_cache/html/unnamed-chunk-15_c9fad914b2da127567b8036becee7f07.rdx
diff --git a/resistant_lines_cache/html/unnamed-chunk-16_c70a75f3fef1561cb7ed8a2201a63220.rdb b/resistant_lines_cache/html/unnamed-chunk-16_c70a75f3fef1561cb7ed8a2201a63220.rdb
index 9080c93f..226c2675 100644
Binary files a/resistant_lines_cache/html/unnamed-chunk-16_c70a75f3fef1561cb7ed8a2201a63220.rdb and b/resistant_lines_cache/html/unnamed-chunk-16_c70a75f3fef1561cb7ed8a2201a63220.rdb differ
diff --git a/resistant_lines_cache/html/unnamed-chunk-16_c70a75f3fef1561cb7ed8a2201a63220.rdx b/resistant_lines_cache/html/unnamed-chunk-16_c70a75f3fef1561cb7ed8a2201a63220.rdx
index 8c5cd586..e0d4c27a 100644
Binary files a/resistant_lines_cache/html/unnamed-chunk-16_c70a75f3fef1561cb7ed8a2201a63220.rdx and b/resistant_lines_cache/html/unnamed-chunk-16_c70a75f3fef1561cb7ed8a2201a63220.rdx differ
diff --git a/resistant_lines_cache/html/unnamed-chunk-17_2a05749d67b8200d4ea3d0f35a10d198.rdb b/resistant_lines_cache/html/unnamed-chunk-17_2a05749d67b8200d4ea3d0f35a10d198.rdb
index ba27f502..5d2ea414 100644
Binary files a/resistant_lines_cache/html/unnamed-chunk-17_2a05749d67b8200d4ea3d0f35a10d198.rdb and b/resistant_lines_cache/html/unnamed-chunk-17_2a05749d67b8200d4ea3d0f35a10d198.rdb differ
diff --git a/resistant_lines_cache/html/unnamed-chunk-18_98214ea5eb9efd2489658138170116af.RData b/resistant_lines_cache/html/unnamed-chunk-18_98214ea5eb9efd2489658138170116af.RData
deleted file mode 100644
index 59bd1767..00000000
Binary files a/resistant_lines_cache/html/unnamed-chunk-18_98214ea5eb9efd2489658138170116af.RData and /dev/null differ
diff --git a/resistant_lines_cache/html/unnamed-chunk-18_bd557fbeb46fffce48702dba906e7707.RData b/resistant_lines_cache/html/unnamed-chunk-18_bd557fbeb46fffce48702dba906e7707.RData
new file mode 100644
index 00000000..fa7e6811
Binary files /dev/null and b/resistant_lines_cache/html/unnamed-chunk-18_bd557fbeb46fffce48702dba906e7707.RData differ
diff --git a/resistant_lines_cache/html/unnamed-chunk-18_98214ea5eb9efd2489658138170116af.rdb b/resistant_lines_cache/html/unnamed-chunk-18_bd557fbeb46fffce48702dba906e7707.rdb
similarity index 100%
rename from resistant_lines_cache/html/unnamed-chunk-18_98214ea5eb9efd2489658138170116af.rdb
rename to resistant_lines_cache/html/unnamed-chunk-18_bd557fbeb46fffce48702dba906e7707.rdb
diff --git a/resistant_lines_cache/html/unnamed-chunk-18_98214ea5eb9efd2489658138170116af.rdx b/resistant_lines_cache/html/unnamed-chunk-18_bd557fbeb46fffce48702dba906e7707.rdx
similarity index 100%
rename from resistant_lines_cache/html/unnamed-chunk-18_98214ea5eb9efd2489658138170116af.rdx
rename to resistant_lines_cache/html/unnamed-chunk-18_bd557fbeb46fffce48702dba906e7707.rdx
diff --git a/resistant_lines_cache/html/unnamed-chunk-19_48588591ed562775d36db603ab76dc0a.RData b/resistant_lines_cache/html/unnamed-chunk-19_48588591ed562775d36db603ab76dc0a.RData
new file mode 100644
index 00000000..27a4ed94
Binary files /dev/null and b/resistant_lines_cache/html/unnamed-chunk-19_48588591ed562775d36db603ab76dc0a.RData differ
diff --git a/resistant_lines_cache/html/unnamed-chunk-19_5736d0c7f23d2dcbe1a6a29ebaca245d.rdb b/resistant_lines_cache/html/unnamed-chunk-19_48588591ed562775d36db603ab76dc0a.rdb
similarity index 100%
rename from resistant_lines_cache/html/unnamed-chunk-19_5736d0c7f23d2dcbe1a6a29ebaca245d.rdb
rename to resistant_lines_cache/html/unnamed-chunk-19_48588591ed562775d36db603ab76dc0a.rdb
diff --git a/resistant_lines_cache/html/unnamed-chunk-19_5736d0c7f23d2dcbe1a6a29ebaca245d.rdx b/resistant_lines_cache/html/unnamed-chunk-19_48588591ed562775d36db603ab76dc0a.rdx
similarity index 100%
rename from resistant_lines_cache/html/unnamed-chunk-19_5736d0c7f23d2dcbe1a6a29ebaca245d.rdx
rename to resistant_lines_cache/html/unnamed-chunk-19_48588591ed562775d36db603ab76dc0a.rdx
diff --git a/resistant_lines_cache/html/unnamed-chunk-19_5736d0c7f23d2dcbe1a6a29ebaca245d.RData b/resistant_lines_cache/html/unnamed-chunk-19_5736d0c7f23d2dcbe1a6a29ebaca245d.RData
deleted file mode 100644
index c6aa71e3..00000000
Binary files a/resistant_lines_cache/html/unnamed-chunk-19_5736d0c7f23d2dcbe1a6a29ebaca245d.RData and /dev/null differ
diff --git a/resistant_lines_cache/html/unnamed-chunk-1_0e263fe5f98ae08dcce973f6ebb78bcf.RData b/resistant_lines_cache/html/unnamed-chunk-1_0e263fe5f98ae08dcce973f6ebb78bcf.RData
deleted file mode 100644
index be0382cc..00000000
Binary files a/resistant_lines_cache/html/unnamed-chunk-1_0e263fe5f98ae08dcce973f6ebb78bcf.RData and /dev/null differ
diff --git a/resistant_lines_cache/html/unnamed-chunk-1_0e263fe5f98ae08dcce973f6ebb78bcf.rdb b/resistant_lines_cache/html/unnamed-chunk-1_0e263fe5f98ae08dcce973f6ebb78bcf.rdb
deleted file mode 100644
index 9fd2f85d..00000000
Binary files a/resistant_lines_cache/html/unnamed-chunk-1_0e263fe5f98ae08dcce973f6ebb78bcf.rdb and /dev/null differ
diff --git a/resistant_lines_cache/html/unnamed-chunk-1_0e263fe5f98ae08dcce973f6ebb78bcf.rdx b/resistant_lines_cache/html/unnamed-chunk-1_0e263fe5f98ae08dcce973f6ebb78bcf.rdx
deleted file mode 100644
index 184258b4..00000000
Binary files a/resistant_lines_cache/html/unnamed-chunk-1_0e263fe5f98ae08dcce973f6ebb78bcf.rdx and /dev/null differ
diff --git a/resistant_lines_cache/html/unnamed-chunk-21_0a271fc46f5b247bd6ee9b36c233b494.RData b/resistant_lines_cache/html/unnamed-chunk-21_0a271fc46f5b247bd6ee9b36c233b494.RData
new file mode 100644
index 00000000..8d8e76a2
Binary files /dev/null and b/resistant_lines_cache/html/unnamed-chunk-21_0a271fc46f5b247bd6ee9b36c233b494.RData differ
diff --git a/resistant_lines_cache/html/unnamed-chunk-21_f34299528a71404319578f2d9f21372b.rdb b/resistant_lines_cache/html/unnamed-chunk-21_0a271fc46f5b247bd6ee9b36c233b494.rdb
similarity index 100%
rename from resistant_lines_cache/html/unnamed-chunk-21_f34299528a71404319578f2d9f21372b.rdb
rename to resistant_lines_cache/html/unnamed-chunk-21_0a271fc46f5b247bd6ee9b36c233b494.rdb
diff --git a/resistant_lines_cache/html/unnamed-chunk-21_f34299528a71404319578f2d9f21372b.rdx b/resistant_lines_cache/html/unnamed-chunk-21_0a271fc46f5b247bd6ee9b36c233b494.rdx
similarity index 100%
rename from resistant_lines_cache/html/unnamed-chunk-21_f34299528a71404319578f2d9f21372b.rdx
rename to resistant_lines_cache/html/unnamed-chunk-21_0a271fc46f5b247bd6ee9b36c233b494.rdx
diff --git a/resistant_lines_cache/html/unnamed-chunk-21_f34299528a71404319578f2d9f21372b.RData b/resistant_lines_cache/html/unnamed-chunk-21_f34299528a71404319578f2d9f21372b.RData
deleted file mode 100644
index 2dfdc6cf..00000000
Binary files a/resistant_lines_cache/html/unnamed-chunk-21_f34299528a71404319578f2d9f21372b.RData and /dev/null differ
diff --git a/resistant_lines_cache/html/unnamed-chunk-2_fe5120f1633224257ee4d0dcff357d33.RData b/resistant_lines_cache/html/unnamed-chunk-2_fe5120f1633224257ee4d0dcff357d33.RData
index 487eb42b..aa998dc2 100644
Binary files a/resistant_lines_cache/html/unnamed-chunk-2_fe5120f1633224257ee4d0dcff357d33.RData and b/resistant_lines_cache/html/unnamed-chunk-2_fe5120f1633224257ee4d0dcff357d33.RData differ
diff --git a/resistant_lines_cache/html/unnamed-chunk-2_fe5120f1633224257ee4d0dcff357d33.rdb b/resistant_lines_cache/html/unnamed-chunk-2_fe5120f1633224257ee4d0dcff357d33.rdb
index 2bdfa73a..ad115d6c 100644
Binary files a/resistant_lines_cache/html/unnamed-chunk-2_fe5120f1633224257ee4d0dcff357d33.rdb and b/resistant_lines_cache/html/unnamed-chunk-2_fe5120f1633224257ee4d0dcff357d33.rdb differ
diff --git a/resistant_lines_cache/html/unnamed-chunk-2_fe5120f1633224257ee4d0dcff357d33.rdx b/resistant_lines_cache/html/unnamed-chunk-2_fe5120f1633224257ee4d0dcff357d33.rdx
index a3f17e1f..2e9e48cb 100644
Binary files a/resistant_lines_cache/html/unnamed-chunk-2_fe5120f1633224257ee4d0dcff357d33.rdx and b/resistant_lines_cache/html/unnamed-chunk-2_fe5120f1633224257ee4d0dcff357d33.rdx differ
diff --git a/resistant_lines_files/figure-html/unnamed-chunk-10-1.png b/resistant_lines_files/figure-html/unnamed-chunk-10-1.png
index 275a24ff..b6658b22 100644
Binary files a/resistant_lines_files/figure-html/unnamed-chunk-10-1.png and b/resistant_lines_files/figure-html/unnamed-chunk-10-1.png differ
diff --git a/resistant_lines_files/figure-html/unnamed-chunk-11-1.png b/resistant_lines_files/figure-html/unnamed-chunk-11-1.png
index a7a0c067..38d5d42b 100644
Binary files a/resistant_lines_files/figure-html/unnamed-chunk-11-1.png and b/resistant_lines_files/figure-html/unnamed-chunk-11-1.png differ
diff --git a/resistant_lines_files/figure-html/unnamed-chunk-13-1.png b/resistant_lines_files/figure-html/unnamed-chunk-13-1.png
deleted file mode 100644
index 8cba035d..00000000
Binary files a/resistant_lines_files/figure-html/unnamed-chunk-13-1.png and /dev/null differ
diff --git a/resistant_lines_files/figure-html/unnamed-chunk-19-1.png b/resistant_lines_files/figure-html/unnamed-chunk-19-1.png
index e31448cb..812b67ea 100644
Binary files a/resistant_lines_files/figure-html/unnamed-chunk-19-1.png and b/resistant_lines_files/figure-html/unnamed-chunk-19-1.png differ
diff --git a/resistant_lines_files/figure-html/unnamed-chunk-21-1.png b/resistant_lines_files/figure-html/unnamed-chunk-21-1.png
index 9547802f..812b67ea 100644
Binary files a/resistant_lines_files/figure-html/unnamed-chunk-21-1.png and b/resistant_lines_files/figure-html/unnamed-chunk-21-1.png differ
diff --git a/sl_plot.qmd b/sl_plot.qmd
index 8667f68c..87314140 100755
--- a/sl_plot.qmd
+++ b/sl_plot.qmd
@@ -12,7 +12,7 @@ source("libs/Common.R")
 
 
 ```{r echo = FALSE}
-pkg_ver(c("dplyr"))
+pkg_ver(c("dplyr", "ggplot2"))
 ```
 
 
@@ -48,7 +48,7 @@ The singer dataset does not seem to exhibit any dependence between a voice part'
 
 Next, we'll look at an example of a dataset that *does* exhibit a dependence between spread and fitted values.
 
-## R's built in `iris` dataset
+## Example: the `iris` dataset
 
 R has a built-in dataset called `iris` that provide measurements of sepal and petal dimensions for three different species of the iris family. In this next example, we will plot the spreads of the `Petal.Length` residuals (after removing their group median values) to their group medians.
 
diff --git a/sl_plot_cache/html/unnamed-chunk-1_9d27599af8b18542bf589d5de9a52edf.RData b/sl_plot_cache/html/unnamed-chunk-1_9d27599af8b18542bf589d5de9a52edf.RData
deleted file mode 100644
index 1c228d73..00000000
Binary files a/sl_plot_cache/html/unnamed-chunk-1_9d27599af8b18542bf589d5de9a52edf.RData and /dev/null differ
diff --git a/sl_plot_cache/html/unnamed-chunk-1_9d27599af8b18542bf589d5de9a52edf.rdb b/sl_plot_cache/html/unnamed-chunk-1_9d27599af8b18542bf589d5de9a52edf.rdb
deleted file mode 100644
index 9fd2f85d..00000000
Binary files a/sl_plot_cache/html/unnamed-chunk-1_9d27599af8b18542bf589d5de9a52edf.rdb and /dev/null differ
diff --git a/sl_plot_cache/html/unnamed-chunk-1_9d27599af8b18542bf589d5de9a52edf.rdx b/sl_plot_cache/html/unnamed-chunk-1_9d27599af8b18542bf589d5de9a52edf.rdx
deleted file mode 100644
index 184258b4..00000000
Binary files a/sl_plot_cache/html/unnamed-chunk-1_9d27599af8b18542bf589d5de9a52edf.rdx and /dev/null differ
diff --git a/sl_plot_cache/html/unnamed-chunk-2_a9e0c5c25fa201d342586fad77e0b216.RData b/sl_plot_cache/html/unnamed-chunk-2_a9e0c5c25fa201d342586fad77e0b216.RData
new file mode 100644
index 00000000..3734f719
Binary files /dev/null and b/sl_plot_cache/html/unnamed-chunk-2_a9e0c5c25fa201d342586fad77e0b216.RData differ
diff --git a/sl_plot_cache/html/unnamed-chunk-2_a9e0c5c25fa201d342586fad77e0b216.rdb b/sl_plot_cache/html/unnamed-chunk-2_a9e0c5c25fa201d342586fad77e0b216.rdb
new file mode 100644
index 00000000..ad115d6c
Binary files /dev/null and b/sl_plot_cache/html/unnamed-chunk-2_a9e0c5c25fa201d342586fad77e0b216.rdb differ
diff --git a/sl_plot_cache/html/unnamed-chunk-2_a9e0c5c25fa201d342586fad77e0b216.rdx b/sl_plot_cache/html/unnamed-chunk-2_a9e0c5c25fa201d342586fad77e0b216.rdx
new file mode 100644
index 00000000..03351299
Binary files /dev/null and b/sl_plot_cache/html/unnamed-chunk-2_a9e0c5c25fa201d342586fad77e0b216.rdx differ
diff --git a/sl_plot_cache/html/unnamed-chunk-2_ddef71e58918eb94c9077c76bf467040.RData b/sl_plot_cache/html/unnamed-chunk-2_ddef71e58918eb94c9077c76bf467040.RData
deleted file mode 100755
index 1e5b130c..00000000
Binary files a/sl_plot_cache/html/unnamed-chunk-2_ddef71e58918eb94c9077c76bf467040.RData and /dev/null differ
diff --git a/sl_plot_cache/html/unnamed-chunk-2_ddef71e58918eb94c9077c76bf467040.rdb b/sl_plot_cache/html/unnamed-chunk-2_ddef71e58918eb94c9077c76bf467040.rdb
deleted file mode 100755
index 2bdfa73a..00000000
Binary files a/sl_plot_cache/html/unnamed-chunk-2_ddef71e58918eb94c9077c76bf467040.rdb and /dev/null differ
diff --git a/sl_plot_cache/html/unnamed-chunk-2_ddef71e58918eb94c9077c76bf467040.rdx b/sl_plot_cache/html/unnamed-chunk-2_ddef71e58918eb94c9077c76bf467040.rdx
deleted file mode 100755
index b57d653d..00000000
Binary files a/sl_plot_cache/html/unnamed-chunk-2_ddef71e58918eb94c9077c76bf467040.rdx and /dev/null differ
diff --git a/sl_plot_cache/html/unnamed-chunk-3_9187cec471b10956d63ffc9a30e6f07a.RData b/sl_plot_cache/html/unnamed-chunk-3_9187cec471b10956d63ffc9a30e6f07a.RData
index fab7d80f..f8bca39c 100755
Binary files a/sl_plot_cache/html/unnamed-chunk-3_9187cec471b10956d63ffc9a30e6f07a.RData and b/sl_plot_cache/html/unnamed-chunk-3_9187cec471b10956d63ffc9a30e6f07a.RData differ
diff --git a/sl_plot_cache/html/unnamed-chunk-4_c31935faa3125f3844ba51b516baafae.RData b/sl_plot_cache/html/unnamed-chunk-4_c31935faa3125f3844ba51b516baafae.RData
index d042ec9d..1ff0c318 100644
Binary files a/sl_plot_cache/html/unnamed-chunk-4_c31935faa3125f3844ba51b516baafae.RData and b/sl_plot_cache/html/unnamed-chunk-4_c31935faa3125f3844ba51b516baafae.RData differ
diff --git a/sl_plot_cache/html/unnamed-chunk-5_e930bbfdb7e1af934707bcd46bff54b4.RData b/sl_plot_cache/html/unnamed-chunk-5_e930bbfdb7e1af934707bcd46bff54b4.RData
index 363bdb0a..2b7fe0a9 100644
Binary files a/sl_plot_cache/html/unnamed-chunk-5_e930bbfdb7e1af934707bcd46bff54b4.RData and b/sl_plot_cache/html/unnamed-chunk-5_e930bbfdb7e1af934707bcd46bff54b4.RData differ
diff --git a/sl_plot_cache/html/unnamed-chunk-6_6d3fd12613d9bfbd45d59b658e228128.RData b/sl_plot_cache/html/unnamed-chunk-6_6d3fd12613d9bfbd45d59b658e228128.RData
index 3fd75a18..7e626887 100644
Binary files a/sl_plot_cache/html/unnamed-chunk-6_6d3fd12613d9bfbd45d59b658e228128.RData and b/sl_plot_cache/html/unnamed-chunk-6_6d3fd12613d9bfbd45d59b658e228128.RData differ
diff --git a/sl_plot_cache/html/unnamed-chunk-7_2eb9368a33f85e5ca3bb37b49231d052.RData b/sl_plot_cache/html/unnamed-chunk-7_2eb9368a33f85e5ca3bb37b49231d052.RData
index a89561e4..ec5d9767 100644
Binary files a/sl_plot_cache/html/unnamed-chunk-7_2eb9368a33f85e5ca3bb37b49231d052.RData and b/sl_plot_cache/html/unnamed-chunk-7_2eb9368a33f85e5ca3bb37b49231d052.RData differ
diff --git a/sl_plot_files/figure-html/unnamed-chunk-3-1.png b/sl_plot_files/figure-html/unnamed-chunk-3-1.png
index e9292ca6..d1ef501b 100755
Binary files a/sl_plot_files/figure-html/unnamed-chunk-3-1.png and b/sl_plot_files/figure-html/unnamed-chunk-3-1.png differ
diff --git a/sl_plot_files/figure-html/unnamed-chunk-4-1.png b/sl_plot_files/figure-html/unnamed-chunk-4-1.png
index 8639b65a..b54cbc3e 100755
Binary files a/sl_plot_files/figure-html/unnamed-chunk-4-1.png and b/sl_plot_files/figure-html/unnamed-chunk-4-1.png differ
diff --git a/sl_plot_files/figure-html/unnamed-chunk-7-1.png b/sl_plot_files/figure-html/unnamed-chunk-7-1.png
index 626d5267..74d9fea7 100644
Binary files a/sl_plot_files/figure-html/unnamed-chunk-7-1.png and b/sl_plot_files/figure-html/unnamed-chunk-7-1.png differ
diff --git a/two_way.qmd b/two_way.qmd
index 35c66c81..e36bdabd 100644
--- a/two_way.qmd
+++ b/two_way.qmd
@@ -53,15 +53,15 @@ dotchart( as.matrix(df)- median(df$Income), pch=16,cex=0.6)
 abline(v = 0, col="red")
 ```
 
-We have separated the income values into two parts: an **overall** value and the difference between each value and the overall--the **residuals**. This can be written in the following algebraic form:
+We have separated the income values into two parts: an **overall** value and the difference between each value and the overall--i.e. the **residuals**. This can be written in the following algebraic form:
 
 $$
 Income_{Occupation} = Income_{median} + Income_{residual}
 $$
 
-For example, 'Healthcare Support' average income of  $`r sprintf("%i",df$Income[12])` can be separated into the the median, $`r sprintf("%i",median(df$Income))`, and the residual, $`r sprintf("%i",df$Income[12]-median(df$Income))`. 
+For example, 'Healthcare Support' average income of  $`r sprintf("%i",df$Income[12])` can be separated into the median, $`r sprintf("%i",median(df$Income))`, and its residual, $`r sprintf("%i",df$Income[12]-median(df$Income))`. 
 
-So in addition to gleaning a median income value from this batch of values, we have also devised a way to facilitate the comparison of the income values to one another. Alternatively, we could have assigned colors to each cell of the table reflecting their values *relative* to the overall median.
+So, in addition to gleaning a median income value from this batch of values, we have also devised a way to facilitate the comparison of the income values to one another. Alternatively, we could have assigned colors to each cell of the table reflecting their values *relative* to the overall median.
 
 ```{r small.mar=TRUE, fig.width=4,fig.height=5, echo=FALSE}
 library(RColorBrewer)
@@ -233,7 +233,7 @@ The final version of the table (along with the column and row values) is shown b
 
 ![](img/twoway_09.jpg)
 
-A progression of each iteration is shown side-by-side in the next figure. Note how, at the end of the iterations, the *common factor* ends up contributing the largest effect. Also note the one abnormally high residual in the southern region associated with 8^th^ grade education.
+A progression of each iteration is shown side-by-side in the next figure. The left-most table is the original data showing death rates. The second table shows the outcome of the first round of polishing (including the initial overall median value of 20.2). The third and forth table show the second and third iterations of the smoothing operations. Additional iterations are not deemed necessary given that little more can be extracted from the residuals.
 
 ```{r echo = FALSE, fig.width = 11, fig.height=2.5}
 library(tukeyedar)
@@ -259,7 +259,11 @@ OP <- par(mfrow = c(1,4))
 par(OP)
 ```
 
-The left-most table is the original data showing death rates. The second table shows the outcome of the first round of polishing (including the initial overall median value of 20.2). The third and forth table show the second and third iterations of the smoothing operations. Additional iterations are not deemed necessary given that little more can be extracted from the residuals.
+A few takeway observations from the analysis follow:
+
+ * The *common factor* ends up contributing the largest effect. 
+ * Between `Region` and `Education` effects, the latter explains more of the variability in deaths with a range of 7.43 - (-3.7) = 11.13 deaths per 1000 births. The `Region` effect explains about 2.3 - (-1.46) = 3.76 deaths per 1000 births of the variability in the data.
+ * One abnormally high residual--10.97 deaths above the overall median deaths--is associated with the southern region and an 8^th^ grade education.
 
 ### Interpreting the median polish
 
@@ -356,7 +360,7 @@ df.pol <- eda_pol(dfl, row = Region, col = Education , val = Deaths,
                    res.size = 0.8, row.size = 0.8, col.size = 0.8)
 ```
 
-Here, the `eda_pol` function required five iterations which explains the slight difference in values from those generated with the medpolish function. You can control the maximum number of iterations in `eda_pol` with the `maxiter` argument.
+Here, the `eda_pol` function required five iterations which explains the slight difference in values from those generated with the `medpolish` function. You can control the maximum number of iterations in `eda_pol` with the `maxiter` argument.
 
 By default, the column and row values are listed in the same order as they are presented in the input data frame (either defined alphabetically or defined by their levels). You can have the function sort the values by their effect size using the `sort` argument.
 
@@ -383,7 +387,7 @@ The `eda_pol` will also output a diagnostic plot.
 plot(df.pol, type = "diagnostic")
 ```
 
-In addition to generating the plot, the function will also return the slope of the fitted line.
+In addition to generating the plot, the function will also return the slope of the fitted line in the console.
 
 For a comprehensive review of the median polish and the `eda_pol` function, see the `eda_pol` [vignette](https://mgimond.github.io/tukeyedar/articles/polish.html).
 
diff --git a/two_way_cache/html/__packages b/two_way_cache/html/__packages
index 5b3351c0..e90ec1ea 100644
--- a/two_way_cache/html/__packages
+++ b/two_way_cache/html/__packages
@@ -3,5 +3,4 @@ RColorBrewer
 ggplot2
 tidyr
 dplyr
-kableExtra
 tukeyedar
diff --git a/two_way_cache/html/clock_a80bfd981f796072ed50d2d6f8d5d27e.RData b/two_way_cache/html/clock_a80bfd981f796072ed50d2d6f8d5d27e.RData
deleted file mode 100644
index de005e13..00000000
Binary files a/two_way_cache/html/clock_a80bfd981f796072ed50d2d6f8d5d27e.RData and /dev/null differ
diff --git a/two_way_cache/html/clock_a80bfd981f796072ed50d2d6f8d5d27e.rdb b/two_way_cache/html/clock_a80bfd981f796072ed50d2d6f8d5d27e.rdb
deleted file mode 100644
index e69de29b..00000000
diff --git a/two_way_cache/html/clock_a80bfd981f796072ed50d2d6f8d5d27e.rdx b/two_way_cache/html/clock_a80bfd981f796072ed50d2d6f8d5d27e.rdx
deleted file mode 100644
index 77f93c8d..00000000
Binary files a/two_way_cache/html/clock_a80bfd981f796072ed50d2d6f8d5d27e.rdx and /dev/null differ
diff --git a/two_way_cache/html/unnamed-chunk-13_32619a5000872bc3f93103779d81e8c4.RData b/two_way_cache/html/unnamed-chunk-13_32619a5000872bc3f93103779d81e8c4.RData
index 165cdbb3..b8f819c8 100644
Binary files a/two_way_cache/html/unnamed-chunk-13_32619a5000872bc3f93103779d81e8c4.RData and b/two_way_cache/html/unnamed-chunk-13_32619a5000872bc3f93103779d81e8c4.RData differ
diff --git a/two_way_cache/html/unnamed-chunk-1_6737535602c30bc1c6e50ac19f5d2b0c.RData b/two_way_cache/html/unnamed-chunk-1_6737535602c30bc1c6e50ac19f5d2b0c.RData
deleted file mode 100644
index 78aa8874..00000000
Binary files a/two_way_cache/html/unnamed-chunk-1_6737535602c30bc1c6e50ac19f5d2b0c.RData and /dev/null differ
diff --git a/two_way_cache/html/unnamed-chunk-1_6737535602c30bc1c6e50ac19f5d2b0c.rdb b/two_way_cache/html/unnamed-chunk-1_6737535602c30bc1c6e50ac19f5d2b0c.rdb
deleted file mode 100644
index 9fd2f85d..00000000
Binary files a/two_way_cache/html/unnamed-chunk-1_6737535602c30bc1c6e50ac19f5d2b0c.rdb and /dev/null differ
diff --git a/two_way_cache/html/unnamed-chunk-1_6737535602c30bc1c6e50ac19f5d2b0c.rdx b/two_way_cache/html/unnamed-chunk-1_6737535602c30bc1c6e50ac19f5d2b0c.rdx
deleted file mode 100644
index 184258b4..00000000
Binary files a/two_way_cache/html/unnamed-chunk-1_6737535602c30bc1c6e50ac19f5d2b0c.rdx and /dev/null differ
diff --git a/two_way_cache/html/unnamed-chunk-22_9857b00948e76646f2b42871c3c01715.RData b/two_way_cache/html/unnamed-chunk-22_9857b00948e76646f2b42871c3c01715.RData
index d0b5b0f5..65ecec06 100644
Binary files a/two_way_cache/html/unnamed-chunk-22_9857b00948e76646f2b42871c3c01715.RData and b/two_way_cache/html/unnamed-chunk-22_9857b00948e76646f2b42871c3c01715.RData differ
diff --git a/two_way_cache/html/unnamed-chunk-22_9857b00948e76646f2b42871c3c01715.rdb b/two_way_cache/html/unnamed-chunk-22_9857b00948e76646f2b42871c3c01715.rdb
index 8257839b..1c85c378 100644
Binary files a/two_way_cache/html/unnamed-chunk-22_9857b00948e76646f2b42871c3c01715.rdb and b/two_way_cache/html/unnamed-chunk-22_9857b00948e76646f2b42871c3c01715.rdb differ
diff --git a/two_way_cache/html/unnamed-chunk-22_9857b00948e76646f2b42871c3c01715.rdx b/two_way_cache/html/unnamed-chunk-22_9857b00948e76646f2b42871c3c01715.rdx
index 04e95d5d..7e025770 100644
Binary files a/two_way_cache/html/unnamed-chunk-22_9857b00948e76646f2b42871c3c01715.rdx and b/two_way_cache/html/unnamed-chunk-22_9857b00948e76646f2b42871c3c01715.rdx differ
diff --git a/two_way_cache/html/unnamed-chunk-23_ce64700903bec89d35ce4c76b73761d0.RData b/two_way_cache/html/unnamed-chunk-23_ce64700903bec89d35ce4c76b73761d0.RData
index 7949a8fb..0132b436 100644
Binary files a/two_way_cache/html/unnamed-chunk-23_ce64700903bec89d35ce4c76b73761d0.RData and b/two_way_cache/html/unnamed-chunk-23_ce64700903bec89d35ce4c76b73761d0.RData differ
diff --git a/two_way_cache/html/unnamed-chunk-2_d6b817bd2dc7f751b645b175972dbeb2.RData b/two_way_cache/html/unnamed-chunk-2_d6b817bd2dc7f751b645b175972dbeb2.RData
index a601f4fd..7559b8bb 100644
Binary files a/two_way_cache/html/unnamed-chunk-2_d6b817bd2dc7f751b645b175972dbeb2.RData and b/two_way_cache/html/unnamed-chunk-2_d6b817bd2dc7f751b645b175972dbeb2.RData differ
diff --git a/two_way_cache/html/unnamed-chunk-2_d6b817bd2dc7f751b645b175972dbeb2.rdb b/two_way_cache/html/unnamed-chunk-2_d6b817bd2dc7f751b645b175972dbeb2.rdb
index 2bdfa73a..ad115d6c 100644
Binary files a/two_way_cache/html/unnamed-chunk-2_d6b817bd2dc7f751b645b175972dbeb2.rdb and b/two_way_cache/html/unnamed-chunk-2_d6b817bd2dc7f751b645b175972dbeb2.rdb differ
diff --git a/two_way_cache/html/unnamed-chunk-2_d6b817bd2dc7f751b645b175972dbeb2.rdx b/two_way_cache/html/unnamed-chunk-2_d6b817bd2dc7f751b645b175972dbeb2.rdx
index 40ea98e1..81091615 100644
Binary files a/two_way_cache/html/unnamed-chunk-2_d6b817bd2dc7f751b645b175972dbeb2.rdx and b/two_way_cache/html/unnamed-chunk-2_d6b817bd2dc7f751b645b175972dbeb2.rdx differ
diff --git a/two_way_cache/html/unnamed-chunk-3_bb5971ce07d8d8c0b8943589b38c4d73.RData b/two_way_cache/html/unnamed-chunk-3_bb5971ce07d8d8c0b8943589b38c4d73.RData
index 47069965..4dfb420d 100644
Binary files a/two_way_cache/html/unnamed-chunk-3_bb5971ce07d8d8c0b8943589b38c4d73.RData and b/two_way_cache/html/unnamed-chunk-3_bb5971ce07d8d8c0b8943589b38c4d73.RData differ
diff --git a/two_way_files/figure-html/clock-1.png b/two_way_files/figure-html/clock-1.png
deleted file mode 100644
index e4bbe66e..00000000
Binary files a/two_way_files/figure-html/clock-1.png and /dev/null differ
diff --git a/two_way_files/figure-html/unnamed-chunk-13-1.png b/two_way_files/figure-html/unnamed-chunk-13-1.png
deleted file mode 100644
index 2691a1ce..00000000
Binary files a/two_way_files/figure-html/unnamed-chunk-13-1.png and /dev/null differ
diff --git a/two_way_files/figure-html/unnamed-chunk-14-1.png b/two_way_files/figure-html/unnamed-chunk-14-1.png
index 2691a1ce..2e33439e 100644
Binary files a/two_way_files/figure-html/unnamed-chunk-14-1.png and b/two_way_files/figure-html/unnamed-chunk-14-1.png differ
diff --git a/two_way_files/figure-html/unnamed-chunk-15-1.png b/two_way_files/figure-html/unnamed-chunk-15-1.png
deleted file mode 100644
index 8e220fa6..00000000
Binary files a/two_way_files/figure-html/unnamed-chunk-15-1.png and /dev/null differ
diff --git a/two_way_files/figure-html/unnamed-chunk-16-1.png b/two_way_files/figure-html/unnamed-chunk-16-1.png
deleted file mode 100644
index 8e220fa6..00000000
Binary files a/two_way_files/figure-html/unnamed-chunk-16-1.png and /dev/null differ
diff --git a/two_way_files/figure-html/unnamed-chunk-18-1.png b/two_way_files/figure-html/unnamed-chunk-18-1.png
deleted file mode 100644
index e03bef5a..00000000
Binary files a/two_way_files/figure-html/unnamed-chunk-18-1.png and /dev/null differ
diff --git a/two_way_files/figure-html/unnamed-chunk-20-1.png b/two_way_files/figure-html/unnamed-chunk-20-1.png
deleted file mode 100644
index 2575dc35..00000000
Binary files a/two_way_files/figure-html/unnamed-chunk-20-1.png and /dev/null differ