<p>In the <a href="http://pacha.hk/anova.html">first part</a> of this entry I did show some <code>mtcars</code> examples, where <code>am</code> can be useful to explain ANOVA as its observations are defined as:
</p>
<div>$$
am_i = \begin{cases}1 &amp;\text{ if car } i \text{ is manual} \cr 0 &amp;\text{ if car } i \text{ is automatic}\end{cases}
$$</div>
<p>Now I'll show another example to continue the last example from Part 1 and I'll move to something involved more variables.</p>
<h1>ANOVA with one dummy variable</h1>
<p>Consider a model where the outcome is <code>mpg</code> and the design matrix is <span>\(X = (\vec{1} \: \vec{x}_2)\)</span>.</p>
<p>From the last entry let
</p>
<div>$$
x_2 = \begin{cases}1 &amp;\text{ if car } i \text{ is automatic} \cr 0 &amp;\text{ if car } i \text{ is manual}\end{cases}
$$</div>
<p>This will lead to this estimate:
</p>
<div>$$
\hat{\vec{\beta}} = \begin{bmatrix}\bar{y}_1 \cr \bar{y}_2 - \bar{y}_1\end{bmatrix}
$$</div>
<p>
Fitting the model gives:</p>
<div><pre><span></span>y <span>=</span> mtcars<span>$</span>mpg
x1 <span>=</span> mtcars<span>$</span>am
x2 <span>=</span> <span>ifelse</span><span>(</span>x1 <span>==</span> <span>1</span><span>,</span> <span>0</span><span>,</span> <span>1</span><span>)</span>

fit <span>=</span> lm<span>(</span>y <span>~</span> x2<span>)</span>
<span>summary</span><span>(</span>fit<span>)</span>
</pre></div>
<div><pre><span></span><span>Call</span><span>:</span>
<span>lm</span><span>(</span><span>formula</span> <span>=</span> <span>y</span> <span>~</span> <span>x2</span><span>)</span>

<span>Residuals</span><span>:</span>
    <span>Min</span>      <span>1</span><span>Q</span>  <span>Median</span>      <span>3</span><span>Q</span>     <span>Max</span> 
<span>-</span><span>9.3923</span> <span>-</span><span>3.0923</span> <span>-</span><span>0.2974</span>  <span>3.2439</span>  <span>9.5077</span> 

<span>Coefficients</span><span>:</span>
            <span>Estimate</span> <span>Std</span><span>.</span> <span>Error</span> <span>t</span> <span>value</span> <span>Pr</span><span>(>|</span><span>t</span><span>|)</span>    
<span>(</span><span>Intercept</span><span>)</span>   <span>24.392</span>      <span>1.360</span>  <span>17.941</span>  <span>&lt;</span> <span>2</span><span>e</span><span>-</span><span>16</span> <span>***</span>
<span>x2</span>            <span>-</span><span>7.245</span>      <span>1.764</span>  <span>-</span><span>4.106</span> <span>0.000285</span> <span>***</span>
<span>---</span>
<span>Signif</span><span>.</span> <span>codes</span><span>:</span>  <span>0</span> <span>'***'</span> <span>0.001</span> <span>'**'</span> <span>0.01</span> <span>'*'</span> <span>0.05</span> <span>'.'</span> <span>0.1</span> <span>' '</span> <span>1</span>

<span>Residual</span> <span>standard</span> <span>error</span><span>:</span> <span>4.902</span> <span>on</span> <span>30</span> <span>degrees</span> <span>of</span> <span>freedom</span>
<span>Multiple</span> <span>R</span><span>-</span><span>squared</span><span>:</span>  <span>0.3598</span><span>,</span>    <span>Adjusted</span> <span>R</span><span>-</span><span>squared</span><span>:</span>  <span>0.3385</span> 
<span>F</span><span>-</span><span>statistic</span><span>:</span> <span>16.86</span> <span>on</span> <span>1</span> <span>and</span> <span>30</span> <span>DF</span><span>,</span>  <span>p</span><span>-</span><span>value</span><span>:</span> <span>0.000285</span>
</pre></div>
<p>So to see the relationship between the estimates and the group means I need additional steps:</p>
<div><pre><span></span>x0 <span>=</span> <span>rep</span><span>(</span><span>1</span><span>,</span><span>length</span><span>(</span>y<span>))</span>
X <span>=</span> <span>cbind</span><span>(</span>x0<span>,</span>x2<span>)</span>
beta <span>=</span> <span>solve</span><span>(</span><span>t</span><span>(</span>X<span>)</span><span>%*%</span>X<span>)</span> <span>%*%</span> <span>(</span><span>t</span><span>(</span>X<span>)</span><span>%*%</span>y<span>)</span>
<span>beta</span>
</pre></div>
<div><pre><span></span>        [,1]
x0 24.392308
x2 -7.244939
</pre></div>
<p>I did obtain the same estimates with <code>lm</code> command so now I calculate the group means:</p>
<div><pre><span></span>x1 <span>=</span> <span>ifelse</span><span>(</span>x1 <span>==</span> <span>0</span><span>,</span> <span>NA</span><span>,</span> x1<span>)</span>
x2 <span>=</span> <span>ifelse</span><span>(</span>x2 <span>==</span> <span>0</span><span>,</span> <span>NA</span><span>,</span> x2<span>)</span>

m1 <span>=</span> <span>mean</span><span>(</span>y<span>*</span>x1<span>,</span> na.rm <span>=</span> <span>TRUE</span><span>)</span>
m2 <span>=</span> <span>mean</span><span>(</span>y<span>*</span>x2<span>,</span> na.rm <span>=</span> <span>TRUE</span><span>)</span>

beta0 <span>=</span> m1
beta2 <span>=</span> m2<span>-</span>m1

beta0<span>;</span>beta2
</pre></div>
<div><pre><span></span>[1] 24.39231
</pre></div>
<div><pre><span></span>[1] -7.244939
</pre></div>
<p>In this case this means that the slope for the two groups is the same but the intercept is different, and therefore exists a negative effect of automatic transmission on miles per gallon in average terms.</p>
<p>Again I'll verify the equivalency between <code>lm</code> and <code>aov</code> in this particular case:</p>
<div><pre><span></span>y <span>=</span> mtcars<span>$</span>mpg
x1 <span>=</span> mtcars<span>$</span>am
x2 <span>=</span> <span>ifelse</span><span>(</span>x1 <span>==</span> <span>1</span><span>,</span> <span>0</span><span>,</span> <span>1</span><span>)</span>

fit2 <span>=</span> aov<span>(</span>y <span>~</span> x2<span>)</span>
fit2<span>$</span>coefficients
</pre></div>
<div><pre><span></span>(Intercept)          x2 
  24.392308   -7.244939 
</pre></div>
<p>I can calculate the residuals by hand:</p>
<div><pre><span></span>mean_mpg <span>=</span> <span>mean</span><span>(</span>mtcars<span>$</span>mpg<span>)</span>
fitted_mpg <span>=</span> fit3<span>$</span>coefficients<span>[</span><span>1</span><span>]</span> <span>+</span> fit3<span>$</span>coefficients<span>[</span><span>2</span><span>]</span><span>*</span>mtcars<span>$</span>am
observed_mpg <span>=</span> mtcars<span>$</span>mpg

TSS <span>=</span> <span>sum</span><span>((</span>observed_mpg <span>-</span> mean_mpg<span>)</span><span>^</span><span>2</span><span>)</span> 
ESS <span>=</span> <span>sum</span><span>((</span>fitted_mpg <span>-</span> mean_mpg<span>)</span><span>^</span><span>2</span><span>)</span>
RSS <span>=</span> <span>sum</span><span>((</span>observed_mpg <span>-</span> fitted_mpg<span>)</span><span>^</span><span>2</span><span>)</span>

TSS<span>;</span>ESS<span>;</span>RSS
</pre></div>
<div><pre><span></span>[1] 1126.047
</pre></div>
<div><pre><span></span>[1] 405.1506
</pre></div>
<div><pre><span></span>[1] 720.8966
</pre></div>
<p>Here its verified that <span>\(TSS = ESS + RSS\)</span> but aside from that I can extract information from <code>aov</code>:</p>
<div><pre><span></span><span>summary</span><span>(</span>fit2<span>)</span>
</pre></div>
<div><pre><span></span>            Df Sum Sq Mean Sq F value   Pr(>F)    
x2           1  405.2   405.2   16.86 0.000285 ***
Residuals   30  720.9    24.0                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
</pre></div>
<p>And check that, as expected, <span>\(ESS\)</span> is the variance explained by <code>x2</code>. </p>
<p>I also can run ANOVA over <code>lm</code> with:</p>
<div><pre><span></span>anova<span>(</span>fit<span>)</span>
</pre></div>
<div><pre><span></span>Analysis of Variance Table

Response: y
          Df Sum Sq Mean Sq F value   Pr(>F)    
x2         1 405.15  405.15   16.86 0.000285 ***
Residuals 30 720.90   24.03                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
</pre></div>
<p>The table provides information on the effect of <code>x2</code> over <code>y</code>. In this case the null hypothesis is rejected because of the large F-value and the associated p-values. </p>
<p>Considering a 0.05 significance threshold I can say, with 95% of confidence, that the regression slope is statistically different from zero or that there is a difference in group means between automatic and manual transmission.</p>
<h1>ANOVA with three dummy variables</h1>
<p>Now let's explore something more complex than <code>am</code>. Reading the documentation I wonder if <code>cyl</code> has an impact on <code>mpg</code> so I explore that variable:</p>
<div><pre><span></span>str<span>(</span>mtcars<span>)</span>
</pre></div>
<div><pre><span></span>'data.frame':   32 obs. of  11 variables:
 $ mpg : num  21 21 22.8 21.4 18.7 18.1 14.3 24.4 22.8 19.2 ...
 $ cyl : num  6 6 4 6 8 6 8 4 4 6 ...
 $ disp: num  160 160 108 258 360 ...
 $ hp  : num  110 110 93 110 175 105 245 62 95 123 ...
 $ drat: num  3.9 3.9 3.85 3.08 3.15 2.76 3.21 3.69 3.92 3.92 ...
 $ wt  : num  2.62 2.88 2.32 3.21 3.44 ...
 $ qsec: num  16.5 17 18.6 19.4 17 ...
 $ vs  : num  0 0 1 1 0 1 0 1 1 1 ...
 $ am  : num  1 1 1 0 0 0 0 0 0 0 ...
 $ gear: num  4 4 4 3 3 3 3 4 4 4 ...
 $ carb: num  4 4 1 1 2 1 4 2 2 4 ...
</pre></div>
<div><pre><span></span><span>unique</span><span>(</span>mtcars<span>$</span>cyl<span>)</span>
</pre></div>
<div><pre><span></span>[1] 6 4 8
</pre></div>
<p>One (wrong) possibility is to write:</p>
<div><pre><span></span>y <span>=</span> mtcars<span>$</span>mpg
x1 <span>=</span> mtcars<span>$</span>cyl<span>;</span> x1 <span>=</span> <span>ifelse</span><span>(</span>x1 <span>==</span> <span>4</span><span>,</span> <span>1</span><span>,</span> <span>0</span><span>)</span>
x2 <span>=</span> mtcars<span>$</span>cyl<span>;</span> x2 <span>=</span> <span>ifelse</span><span>(</span>x1 <span>==</span> <span>6</span><span>,</span> <span>1</span><span>,</span> <span>0</span><span>)</span>
x3 <span>=</span> mtcars<span>$</span>cyl<span>;</span> x3 <span>=</span> <span>ifelse</span><span>(</span>x1 <span>==</span> <span>8</span><span>,</span> <span>1</span><span>,</span> <span>0</span><span>)</span>

fit <span>=</span> lm<span>(</span>y <span>~</span> x1 <span>+</span> x2 <span>+</span> x3<span>)</span>
<span>summary</span><span>(</span>fit<span>)</span>
</pre></div>
<div><pre><span></span><span>Call</span><span>:</span>
<span>lm</span><span>(</span><span>formula</span> <span>=</span> <span>y</span> <span>~</span> <span>x1</span> <span>+</span> <span>x2</span> <span>+</span> <span>x3</span><span>)</span>

<span>Residuals</span><span>:</span>
    <span>Min</span>      <span>1</span><span>Q</span>  <span>Median</span>      <span>3</span><span>Q</span>     <span>Max</span> 
<span>-</span><span>6.2476</span> <span>-</span><span>2.2846</span> <span>-</span><span>0.4556</span>  <span>2.6774</span>  <span>7.2364</span> 

<span>Coefficients</span><span>:</span> <span>(</span><span>2</span> <span>not</span> <span>defined</span> <span>because</span> <span>of</span> <span>singularities</span><span>)</span>
            <span>Estimate</span> <span>Std</span><span>.</span> <span>Error</span> <span>t</span> <span>value</span> <span>Pr</span><span>(>|</span><span>t</span><span>|)</span>    
<span>(</span><span>Intercept</span><span>)</span>  <span>16.6476</span>     <span>0.7987</span>  <span>20.844</span>  <span>&lt;</span> <span>2</span><span>e</span><span>-</span><span>16</span> <span>***</span>
<span>x1</span>           <span>10.0160</span>     <span>1.3622</span>   <span>7.353</span> <span>3.44</span><span>e</span><span>-</span><span>08</span> <span>***</span>
<span>x2</span>                <span>NA</span>         <span>NA</span>      <span>NA</span>       <span>NA</span>    
<span>x3</span>                <span>NA</span>         <span>NA</span>      <span>NA</span>       <span>NA</span>    
<span>---</span>
<span>Signif</span><span>.</span> <span>codes</span><span>:</span>  <span>0</span> <span>'***'</span> <span>0.001</span> <span>'**'</span> <span>0.01</span> <span>'*'</span> <span>0.05</span> <span>'.'</span> <span>0.1</span> <span>' '</span> <span>1</span>

<span>Residual</span> <span>standard</span> <span>error</span><span>:</span> <span>3.66</span> <span>on</span> <span>30</span> <span>degrees</span> <span>of</span> <span>freedom</span>
<span>Multiple</span> <span>R</span><span>-</span><span>squared</span><span>:</span>  <span>0.6431</span><span>,</span>    <span>Adjusted</span> <span>R</span><span>-</span><span>squared</span><span>:</span>  <span>0.6312</span> 
<span>F</span><span>-</span><span>statistic</span><span>:</span> <span>54.06</span> <span>on</span> <span>1</span> <span>and</span> <span>30</span> <span>DF</span><span>,</span>  <span>p</span><span>-</span><span>value</span><span>:</span> <span>3.436</span><span>e</span><span>-</span><span>08</span>
</pre></div>
<p>Here the <code>NAs</code> mean there are variables that are linearly related to the other variables (e.g. the variable pointed with <code>NA</code> is an average of one or more of the rest of the variables, like <span>\(x_2 = 2x_1 + x_3\)</span> or another linear combination), then there's no unique solution to the regression without dropping variables. </p>
<p>But R has a command called <code>as.factor()</code> that is useful in these cases and also can save you some lines of code in other cases:</p>
<div><pre><span></span>fit2 <span>=</span> lm<span>(</span>mpg <span>~</span> <span>as.factor</span><span>(</span>cyl<span>),</span> data <span>=</span> mtcars<span>)</span>
<span>summary</span><span>(</span>fit2<span>)</span>
</pre></div>
<div><pre><span></span><span>Call</span><span>:</span>
<span>lm</span><span>(</span><span>formula</span> <span>=</span> <span>mpg</span> <span>~</span> <span>as</span><span>.</span><span>factor</span><span>(</span><span>cyl</span><span>),</span> <span>data</span> <span>=</span> <span>mtcars</span><span>)</span>

<span>Residuals</span><span>:</span>
    <span>Min</span>      <span>1</span><span>Q</span>  <span>Median</span>      <span>3</span><span>Q</span>     <span>Max</span> 
<span>-</span><span>5.2636</span> <span>-</span><span>1.8357</span>  <span>0.0286</span>  <span>1.3893</span>  <span>7.2364</span> 

<span>Coefficients</span><span>:</span>
                <span>Estimate</span> <span>Std</span><span>.</span> <span>Error</span> <span>t</span> <span>value</span> <span>Pr</span><span>(>|</span><span>t</span><span>|)</span>    
<span>(</span><span>Intercept</span><span>)</span>      <span>26.6636</span>     <span>0.9718</span>  <span>27.437</span>  <span>&lt;</span> <span>2</span><span>e</span><span>-</span><span>16</span> <span>***</span>
<span>as</span><span>.</span><span>factor</span><span>(</span><span>cyl</span><span>)</span><span>6</span>  <span>-</span><span>6.9208</span>     <span>1.5583</span>  <span>-</span><span>4.441</span> <span>0.000119</span> <span>***</span>
<span>as</span><span>.</span><span>factor</span><span>(</span><span>cyl</span><span>)</span><span>8</span> <span>-</span><span>11.5636</span>     <span>1.2986</span>  <span>-</span><span>8.905</span> <span>8.57</span><span>e</span><span>-</span><span>10</span> <span>***</span>
<span>---</span>
<span>Signif</span><span>.</span> <span>codes</span><span>:</span>  <span>0</span> <span>'***'</span> <span>0.001</span> <span>'**'</span> <span>0.01</span> <span>'*'</span> <span>0.05</span> <span>'.'</span> <span>0.1</span> <span>' '</span> <span>1</span>

<span>Residual</span> <span>standard</span> <span>error</span><span>:</span> <span>3.223</span> <span>on</span> <span>29</span> <span>degrees</span> <span>of</span> <span>freedom</span>
<span>Multiple</span> <span>R</span><span>-</span><span>squared</span><span>:</span>  <span>0.7325</span><span>,</span>    <span>Adjusted</span> <span>R</span><span>-</span><span>squared</span><span>:</span>  <span>0.714</span> 
<span>F</span><span>-</span><span>statistic</span><span>:</span>  <span>39.7</span> <span>on</span> <span>2</span> <span>and</span> <span>29</span> <span>DF</span><span>,</span>  <span>p</span><span>-</span><span>value</span><span>:</span> <span>4.979</span><span>e</span><span>-</span><span>09</span>
</pre></div>
<p>The <code>aov</code> version of this is:</p>
<div><pre><span></span>fit3 <span>=</span> aov<span>(</span>mpg <span>~</span> <span>as.factor</span><span>(</span>cyl<span>),</span> data <span>=</span> mtcars<span>)</span>
TukeyHSD<span>(</span>fit3<span>)</span>
</pre></div>
<div><pre><span></span>  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = mpg ~ as.factor(cyl), data = mtcars)

$`as.factor(cyl)`
          diff        lwr        upr     p adj
6-4  -6.920779 -10.769350 -3.0722086 0.0003424
8-4 -11.563636 -14.770779 -8.3564942 0.0000000
8-6  -4.642857  -8.327583 -0.9581313 0.0112287
</pre></div>
<p>As I said many times in this entry, ANOVA is linear regression. Interpreting the coefficients is up to you.</p>