<p><div><div><div><div data-width="100" data-show-faces="false" data-send="true" data-layout="button_count" data-height="21" data-href="https://www.r-bloggers.com/raccoon-ch-2-3-1-way-anova/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="http://www.quantide.com/raccoon-ch-2-3-1-way-anova/"> R blog | Quantide - R training &amp; consulting</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p><em>This article is part of Quantide’s web book “<a target="_blank" rel="nofollow" href="http://www.quantide.com/raccoon-statistical-models-with-r-web-book/">Raccoon – Statistical Models with R</a>“. Raccoon is Quantide’s third web book after “<a target="_blank" rel="nofollow" href="http://www.quantide.com/rabbit-introduction-to-r/">Rabbit – Introduction to R</a>” and “<a target="_blank" rel="nofollow" href="http://www.quantide.com/ramarro-r-for-developers/">Ramarro – R for Developers</a>“. See the <a target="_blank" rel="nofollow" href="http://www.quantide.com/raccoon-statistical-models-r-free-web-book/">full project here</a>.</em></p>
<p><em>The second chapter of Racooon is focused on T-test and Anova. Through example it shows theory and R code of:</em></p>
<ul>
<li><a target="_blank" rel="nofollow" href="http://www.quantide.com/raccoon-ch-2-1-1-sample-t-test-with-r/">1-sample t-test</a></li>
<li><a target="_blank" rel="nofollow" href="http://www.quantide.com/raccoon-ch-2-2-2-sample-t-test-and-paired-t/">2-sample t-test and Paired t</a></li>
<li>1-way Anova</li>
<li>2-way Anova</li>
<li>Unbalanced Anova with one obs dropped</li>
<li>Fixed effects Nested Anova</li>
</ul>
<p><em>This post is the third section of the chapter, about 1-way Anova.</em></p>
<p><em>Throughout the web-book we will widely use the <strong>package qdata</strong>, containing about 80 datasets. You may find it here: <a target="_blank" rel="nofollow" href="https://github.com/quantide/qdata">https://github.com/quantide/qdata</a>.</em></p>
<pre>require(nortest)
require(car)
require(ggplot2)
require(dplyr)
require(gridExtra)
require(tidyr)
require(qdata)</pre>
<p><h2>Example: Tissues (1-way ANOVA)</h2>
<h3>Data description</h3>
<p>Three quality inspectors of a plant, Henry, Anne, and Andrew (operators) measure the strenght of car seat tissues. The company managers want to test the reproducibility, between operators, of company’s measurement system.<br>
The main goal of study is then to verify if operators’ measures are confrontable.<br>
Each operator measured 25 pieces of car seat tissues.<br>
Globally, 75 samples of tissue randomly chosen from the same production batch were measured.</p>
<h3>Data loading</h3>
<pre>data(tissue)
head(tissue)</pre>
<pre>##   Operator Strength
## 1     Anne    11.30
## 2     Anne    10.62
## 3     Anne    10.36
## 4     Anne    10.23
## 5     Anne    10.42
## 6     Anne    12.64</pre>
<pre>str(tissue)</pre>
<pre>## 'data.frame':    75 obs. of  2 variables:
##  $ Operator: Factor w/ 3 levels &quot;Andrew&quot;,&quot;Anne&quot;,..: 2 2 2 2 2 2 2 2 2 2 ...
##  $ Strength: num  11.3 10.6 10.4 10.2 10.4 ...</pre>
<p>Data is already ordered by operator (<code>Anne</code>, <code>Henry</code>, and <code>Andrew</code>)</p>
<h3>Descriptives</h3>
<p>Plot of <code>Strengt</code> with <code>Operator</code> factor variable in abscissa:</p>
<pre>tissue_mean &lt;- tissue %>% group_by(Operator) %>% summarise(Strength=mean(Strength))

ggp &lt;- ggplot() +
  geom_point(data = tissue, mapping = aes(x=Operator, y=Strength), colour=&quot;darkblue&quot;) + 
  geom_line(data=tissue_mean, mapping = aes(x=Operator,y=Strength), colour=&quot;red&quot;, group=1)

print(ggp)</pre>
<p>
<span><img sizes="(max-width: 450px) 100vw, 450px" alt="anova14" class="wp-image-2900" src="https://i2.wp.com/www.quantide.com/wp-content/uploads/2016/09/Anova14-300x300.png?resize=450%2C450" data-recalc-dims="1"></span></p>
<p>Stripchart of Strenght by Operator with connect line</p>
<p>A different way to show similar information is through a plot of univariate effects:</p>
<pre>plot.design(Strength ~ Operator, data = tissue)</pre>
<p><p><img sizes="(max-width: 450px) 100vw, 450px" alt="anova15" class="wp-image-2901" src="https://i1.wp.com/www.quantide.com/wp-content/uploads/2016/09/Anova15-300x300.png?resize=450%2C450" data-recalc-dims="1">Plot of univariate effects of Operator on Strenght</p>
<p>This plot shows the mean for each level of grouping factor compared to the grand mean (the center wider line).</p>
<h3>Inference and models</h3>
<p>Now we will try models with different contrast types.<br>
Initially the model with default contrasts will be fitted to data, then the <code>contr.sum</code> and the <code>contr.helmert</code> contrast models will be shown.</p>
<p>Default contrasts:</p>
<pre>options(&quot;contrasts&quot;)</pre>

<pre>## $contrasts
## [1] &quot;contr.treatment&quot; &quot;contr.poly&quot;</pre>

<pre>options(contrasts = c(&quot;contr.treatment&quot;, &quot;contr.poly&quot;))
fm_treatment &lt;- aov(Strength~Operator, data = tissue)
summary(fm_treatment)</pre>

<pre>##             Df Sum Sq Mean Sq F value Pr(>F)  
## Operator     2   6.62    3.31   3.851 0.0258 *
## Residuals   72  61.90    0.86                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1</pre>

<pre>summary.lm(fm_treatment)</pre>

<pre>## 
## Call:
## aov(formula = Strength ~ Operator, data = tissue)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.9064 -0.5482 -0.0388  0.4656  2.4100 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    10.4364     0.1854  56.280  &lt; 2e-16 ***
## OperatorAnne   -0.2064     0.2622  -0.787  0.43384    
## OperatorHenry  -0.7076     0.2622  -2.698  0.00868 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9272 on 72 degrees of freedom
## Multiple R-squared:  0.09663,    Adjusted R-squared:  0.07154 
## F-statistic: 3.851 on 2 and 72 DF,  p-value: 0.02577</pre>
<p>
<code>aov()</code>’s predefined output gives a global p-value of 0.0258 on <code>Operator</code> effect.</p>
<p>Note that <code>summary.lm()</code> shows the <code>aov</code> resulting object <code>fm_treatment</code> as if it simply were an <code>lm</code>-type object.</p>
<p>From its output, Andrew’s mean appears significantly different from 0 (\(beta_0\) = 10.4364, p = 2e-16),</p>
<p><span>Henry seems significantly different from Andrew (0.7076 lesser, p-value=0.00868), whereas Anne’s seem not significantly different from Andrew (0.2064 lesser, p-value=0.43384).  </span></p>
<p>Instead of <code>summary()</code> and <code>summary.lm()</code>, <code>summary.aov()</code> and <code>summary()</code> may be used, obtaining the same outputs:</p>
<pre>fm_treatment &lt;- lm(Strength~Operator, data = tissue)
summary.aov(fm_treatment)</pre>

<pre>##             Df Sum Sq Mean Sq F value Pr(>F)  
## Operator     2   6.62    3.31   3.851 0.0258 *
## Residuals   72  61.90    0.86                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1</pre>

<pre>summary(fm_treatment)</pre>

<pre>## 
## Call:
## lm(formula = Strength ~ Operator, data = tissue)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.9064 -0.5482 -0.0388  0.4656  2.4100 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    10.4364     0.1854  56.280  &lt; 2e-16 ***
## OperatorAnne   -0.2064     0.2622  -0.787  0.43384    
## OperatorHenry  -0.7076     0.2622  -2.698  0.00868 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9272 on 72 degrees of freedom
## Multiple R-squared:  0.09663,    Adjusted R-squared:  0.07154 
## F-statistic: 3.851 on 2 and 72 DF,  p-value: 0.02577</pre>
<p>
In this sense, <code>aov()</code>,<code>lm()</code> and their resulting objects are substantially interchangeable.</p>
<p>Now let’s see what changes if constrasts settings are modified. The first contrast setting is <code>contr.sum</code>:</p>
<pre>options(contrasts = c(&quot;contr.sum&quot;, &quot;contr.poly&quot;))
fm_sum &lt;- aov(Strength~Operator, data = tissue)
summary(fm_sum)</pre>

<pre>##             Df Sum Sq Mean Sq F value Pr(>F)  
## Operator     2   6.62    3.31   3.851 0.0258 *
## Residuals   72  61.90    0.86                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1</pre>

<pre>summary.lm(fm_sum)</pre>

<pre>## 
## Call:
## aov(formula = Strength ~ Operator, data = tissue)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.9064 -0.5482 -0.0388  0.4656  2.4100 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 10.13173    0.10706  94.635   &lt;2e-16 ***
## Operator1    0.30467    0.15141   2.012   0.0479 *  
## Operator2    0.09827    0.15141   0.649   0.5184    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9272 on 72 degrees of freedom
## Multiple R-squared:  0.09663,    Adjusted R-squared:  0.07154 
## F-statistic: 3.851 on 2 and 72 DF,  p-value: 0.02577</pre>
<p>
<span>The interpretation of the <code>contr.sum</code> coefficients is not very simple: contrasts are the \(k-1\) differences (where \(k\) is the number of factor levels) between each level and a reference level; the resulting estimated coefficients contain the differences of each level mean (excluded the reference level) with respect to the grand mean; in fact, the model intercept is equal to grand mean.  </span></p>
<p><span><code>aov()</code>, however, gives the same results obtained with default contrasts.</span></p>
<p>Now let’s change contrast settings by selecting Helmert contrasts</p>
<pre>options(contrasts = c(&quot;contr.helmert&quot;, &quot;contr.poly&quot;))
fm_helm &lt;- aov(Strength~Operator, data = tissue)
summary(fm_helm)</pre>

<pre>##             Df Sum Sq Mean Sq F value Pr(>F)  
## Operator     2   6.62    3.31   3.851 0.0258 *
## Residuals   72  61.90    0.86                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1</pre>

<pre>summary.lm(fm_helm)</pre>

<pre>## 
## Call:
## aov(formula = Strength ~ Operator, data = tissue)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.9064 -0.5482 -0.0388  0.4656  2.4100 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  10.1317     0.1071  94.635  &lt; 2e-16 ***
## Operator1    -0.1032     0.1311  -0.787  0.43384    
## Operator2    -0.2015     0.0757  -2.661  0.00959 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9272 on 72 degrees of freedom
## Multiple R-squared:  0.09663,    Adjusted R-squared:  0.07154 
## F-statistic: 3.851 on 2 and 72 DF,  p-value: 0.02577</pre>
<p>
<span>Helmert’s contrasts are even more complex than the others; they analyze:  </span></p>
<ol>
<li><span>Initially, the difference between first level mean and the second level mean</span></li>
<li><span>Then, the difference between the mean of first two (pooled) levels and the third level mean</span></li>
<li><span>Then, the difference between the mean of first three (pooled) levels and the fourth level mean</span></li>
<li><span>And so on for \(k-1\) contrasts</span></li>
</ol>
<p>The main advantage of using Helmert contrasts is that they are orthogonal; the tests on linear model coefficients are then stocastically independent: the Type I and Type II error rate for the tests on Helmert coefficients are unrelated. The main disavantage, of course, is their complexity in coefficients explanations.</p>
<p>Now, let’s restore predefined contrasts:</p>
<pre>options(contrasts = c(&quot;contr.treatment&quot;, &quot;contr.poly&quot;))</pre>
<p>
Summarizing the results just seen about contrasts:</p>
<ul>
<li><code>contr.treatment</code> makes easier to read results, but the tests on <em>coefficients</em> are NOT independent</li>
<li><code>contr.sum</code> makes NO independent tests on <em>coefficients</em>, but compares pairs of means; also, the coefficients represent the distance of level’s mean with respect to grand mean, however they are not always easily interpretable</li>
<li><code>contr.helmert</code> makes independent tests on <em>coefficients</em>, but the coefficients are more difficult to interpret</li>
<li>anyway, until now, the global ANOVA (<code>aov()</code>)results are always the same: <code>aov()</code> analyzes effects globally, not single means</li>
</ul>
<h3>Residuals analysis</h3>
<p>By plotting a model object (<code>lm</code> or <code>aov</code>) up to six residuals diagnostic plots may be shown.</p>
<p>By default, the <code>plot()</code> method applied on a model object shows four graphs:</p>
<ol>
<li><span>A plot of raw residuals against fitted values</span></li>
<li><span>A Normal Q-Q plot on standardized residuals</span></li>
<li><span>A scale-Location plot of \(sqrt{vert std. residuals vert}\) against fitted values. Here, \(std. residuals\) means raw residuals divided by \(hatsigma cdot sqrt{1-h_{ii}}\), where \(h_{ii}\) are the diagonal entries of the “hat” matrix (see Appendices chapter).</span></li>
<li><span>A plot of standardized Pearson residuals against leverages. If the leverages are constant (as is the case in the balanced design of this example) the plot uses factor level combinations instead of the leverages for the x-axis.</span></li>
</ol>
<p>The next lines of code split the graphical device in 4 areas (2 x 2) and plot diagnostic graphs in only one display</p>
<pre>op &lt;- par(mfrow = c(2, 2))
plot(fm_treatment)
par(op)</pre>
<p>
<img sizes="(max-width: 450px) 100vw, 450px" alt="anova16" class="wp-image-2902" src="https://i0.wp.com/www.quantide.com/wp-content/uploads/2016/09/Anova16-300x300.png?resize=450%2C450" data-recalc-dims="1"></p>
<p>Compound diagnostic residual plot of Strenght Vs. Operator ANOVA model</p>
<p>The last line of above code restores the graphical device.</p>
<p>The residual plots confirm that normality and homoscedasticity assumptions are met, and no outliers appear. This confirms the model results.</p>
<p> </p>
<h2></h2>
<p>The post <a target="_blank" rel="nofollow" href="http://www.quantide.com/raccoon-ch-2-3-1-way-anova/">Raccoon | Ch 2.3 – 1-way Anova</a> appeared first on <a target="_blank" rel="nofollow" href="http://www.quantide.com/">Quantide – R training &amp; consulting</a>.</p>
<p>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="http://www.quantide.com/raccoon-ch-2-3-1-way-anova/"> R blog | Quantide - R training &amp; consulting</a></strong>.</div>
<hr>
<a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="https://feedburner.google.com/fb/a/mailverify?uri=RBloggers">daily e-mail updates</a></strong> about <a title="The R Project for Statistical Computing" rel="nofollow" href="https://www.r-project.org/">R</a> news and <a title="R tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/tutorial">tutorials</a> on topics such as: <a title="Data science" rel="nofollow" href="https://www.r-bloggers.com/search/data%20science">Data science</a>, <a title="Big Data" rel="nofollow" href="https://www.r-bloggers.com/search/Big%20Data">Big Data, <a title="R jobs" rel="nofollow" href="https://www.r-users.com/">R jobs</a>, visualization (<a title="ggplot and ggplot2 tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/ggplot2">ggplot2</a>, <a title="Boxplots using lattice and ggplot2 tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/boxplot">Boxplots</a>, <a title="Maps and gis" rel="nofollow" href="https://www.r-bloggers.com/search/map">maps</a>, <a title="Animation in R" rel="nofollow" href="https://www.r-bloggers.com/search/animation">animation</a>), programming (<a title="RStudio IDE for R" rel="nofollow" href="https://www.r-bloggers.com/search/RStudio">RStudio</a>, <a title="Sweave and literate programming" rel="nofollow" href="https://www.r-bloggers.com/search/sweave">Sweave</a>, <a title="LaTeX in R" rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX">LaTeX</a>, <a title="SQL and databases" rel="nofollow" href="https://www.r-bloggers.com/search/SQL">SQL</a>, <a title="Eclipse IDE for R" rel="nofollow" href="https://www.r-bloggers.com/search/eclipse">Eclipse</a>, <a title="git and github, Version Control System" rel="nofollow" href="https://www.r-bloggers.com/search/git">git</a>, <a title="Large data in R using Hadoop" rel="nofollow" href="https://www.r-bloggers.com/search/hadoop">hadoop</a>, <a title="Web Scraping of google, facebook, yahoo, twitter and more using R" rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping">Web Scraping</a>) statistics (<a title="Regressions and ANOVA analysis tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/regression">regression</a>, <a title="principal component analysis tutorial" rel="nofollow" href="https://www.r-bloggers.com/search/PCA">PCA</a>, <a title="Time series" rel="nofollow" href="https://www.r-bloggers.com/search/time+series">time series</a>, <a title="finance trading" rel="nofollow" href="https://www.r-bloggers.com/search/trading">trading</a>) and more...
<img height="1" alt="" width="1" src="http://feeds.feedburner.com/~r/RBloggers/~4/cufmc-rAw_g"></a></div>