-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy path09-Anova.html
799 lines (755 loc) · 64.5 KB
/
09-Anova.html
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
<!DOCTYPE html>
<html >
<head>
<meta charset="UTF-8">
<meta http-equiv="X-UA-Compatible" content="IE=edge">
<title>09-anova</title>
<meta name="description" content="09-anova">
<meta name="generator" content="bookdown and GitBook 2.6.7">
<meta property="og:title" content="09-anova" />
<meta property="og:type" content="book" />
<meta name="twitter:card" content="summary" />
<meta name="twitter:title" content="09-anova" />
<meta name="viewport" content="width=device-width, initial-scale=1">
<meta name="apple-mobile-web-app-capable" content="yes">
<meta name="apple-mobile-web-app-status-bar-style" content="black">
<script src="libs/jquery-2.2.3/jquery.min.js"></script>
<link href="libs/gitbook-2.6.7/css/style.css" rel="stylesheet" />
<link href="libs/gitbook-2.6.7/css/plugin-table.css" rel="stylesheet" />
<link href="libs/gitbook-2.6.7/css/plugin-bookdown.css" rel="stylesheet" />
<link href="libs/gitbook-2.6.7/css/plugin-highlight.css" rel="stylesheet" />
<link href="libs/gitbook-2.6.7/css/plugin-search.css" rel="stylesheet" />
<link href="libs/gitbook-2.6.7/css/plugin-fontsettings.css" rel="stylesheet" />
<link href="libs/pagedtable-1.1/css/pagedtable.css" rel="stylesheet" />
<script src="libs/pagedtable-1.1/js/pagedtable.js"></script>
<style type="text/css">
div.sourceCode { overflow-x: auto; }
table.sourceCode, tr.sourceCode, td.lineNumbers, td.sourceCode {
margin: 0; padding: 0; vertical-align: baseline; border: none; }
table.sourceCode { width: 100%; line-height: 100%; }
td.lineNumbers { text-align: right; padding-right: 4px; padding-left: 4px; color: #aaaaaa; border-right: 1px solid #aaaaaa; }
td.sourceCode { padding-left: 5px; }
code > span.kw { color: #007020; font-weight: bold; } /* Keyword */
code > span.dt { color: #902000; } /* DataType */
code > span.dv { color: #40a070; } /* DecVal */
code > span.bn { color: #40a070; } /* BaseN */
code > span.fl { color: #40a070; } /* Float */
code > span.ch { color: #4070a0; } /* Char */
code > span.st { color: #4070a0; } /* String */
code > span.co { color: #60a0b0; font-style: italic; } /* Comment */
code > span.ot { color: #007020; } /* Other */
code > span.al { color: #ff0000; font-weight: bold; } /* Alert */
code > span.fu { color: #06287e; } /* Function */
code > span.er { color: #ff0000; font-weight: bold; } /* Error */
code > span.wa { color: #60a0b0; font-weight: bold; font-style: italic; } /* Warning */
code > span.cn { color: #880000; } /* Constant */
code > span.sc { color: #4070a0; } /* SpecialChar */
code > span.vs { color: #4070a0; } /* VerbatimString */
code > span.ss { color: #bb6688; } /* SpecialString */
code > span.im { } /* Import */
code > span.va { color: #19177c; } /* Variable */
code > span.cf { color: #007020; font-weight: bold; } /* ControlFlow */
code > span.op { color: #666666; } /* Operator */
code > span.bu { } /* BuiltIn */
code > span.ex { } /* Extension */
code > span.pp { color: #bc7a00; } /* Preprocessor */
code > span.at { color: #7d9029; } /* Attribute */
code > span.do { color: #ba2121; font-style: italic; } /* Documentation */
code > span.an { color: #60a0b0; font-weight: bold; font-style: italic; } /* Annotation */
code > span.cv { color: #60a0b0; font-weight: bold; font-style: italic; } /* CommentVar */
code > span.in { color: #60a0b0; font-weight: bold; font-style: italic; } /* Information */
</style>
</head>
<body>
<div class="book without-animation with-summary font-size-2 font-family-1" data-basepath=".">
<div class="book-summary">
<nav role="navigation">
<ul class="summary">
<li><strong><a href="./">MRDA 2018</a></strong></li>
<li class="divider"></li>
<li class="chapter" data-level="0.1" data-path=""><a href="#comparing-several-means"><i class="fa fa-check"></i><b>0.1</b> Comparing several means</a><ul>
<li class="chapter" data-level="0.1.1" data-path=""><a href="#introduction"><i class="fa fa-check"></i><b>0.1.1</b> Introduction</a></li>
<li class="chapter" data-level="0.1.2" data-path=""><a href="#decomposing-variance"><i class="fa fa-check"></i><b>0.1.2</b> Decomposing variance</a></li>
<li class="chapter" data-level="0.1.3" data-path=""><a href="#one-way-anova"><i class="fa fa-check"></i><b>0.1.3</b> One-way ANOVA</a></li>
</ul></li>
</ul>
</nav>
</div>
<div class="book-body">
<div class="body-inner">
<div class="book-header" role="navigation">
<h1>
<i class="fa fa-circle-o-notch fa-spin"></i><a href="./">09-anova</a>
</h1>
</div>
<div class="page-wrapper" tabindex="-1" role="main">
<div class="page-inner">
<section class="normal" id="section-">
<div id="header">
<h1 class="title">09-anova</h1>
</div>
<div id="comparing-several-means" class="section level2">
<h2><span class="header-section-number">0.1</span> Comparing several means</h2>
<p>This chapter is primarily based on Field, A., Miles J., & Field, Z. (2012): Discovering Statistics Using R. Sage Publications, <strong>chapters 10 & 12</strong>.</p>
<p><a href="./Code/08-anova%20(1).R">You can download the corresponding R-Code here</a></p>
<div id="introduction" class="section level3">
<h3><span class="header-section-number">0.1.1</span> Introduction</h3>
<p>In the previous section we learned how to compare means using a t-test. The t-test has some limitations since it only lets you compare 2 means and you can only use it with one independent variable. However, often we would like to compare means from 3 or more groups. In addition, there may be instances in which you manipulate more than one independent variable. For these applications, ANOVA (ANalysis Of VAriance) can be used. Hence, to conduct ANOVA you need:</p>
<ul>
<li>A metric dependent variable (i.e., measured using an interval or ratio scale)</li>
<li>One or more non-metric (categorical) independent variables (also called factors)</li>
</ul>
<p>A treatment is a particular combination of factor levels, or categories. One-way ANOVA is used when there is only one categorical variable (factor). In this case, a treatment is the same as a factor level. N-way ANOVA is used with two or more factors. Note that we are only going to talk about a single independent variable in the context of ANOVA. If you have multiple independent variables please refere to the chapter on <strong>Regression</strong>.</p>
<p>Let’s use an example to see how ANOVA works. Similar to the previous example it is also imaginable that the music streaming service experiments with a recommendation system for user created playlists. We now have three groups, the control group “A” with the current system, treatment group “B” who have access to playlists created by other users but are not shown recommendations and treatment group “C” who are shown recommendations for user created playlists.<br />
As always, we load and inspect the data first:</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">set.seed</span>(<span class="dv">321</span>)
hours_population_<span class="dv">1</span> <-<span class="st"> </span><span class="kw">rnorm</span>(<span class="dv">25000</span>, <span class="dv">15</span>, <span class="dv">5</span>)
<span class="kw">set.seed</span>(<span class="dv">125671</span>)
sample_<span class="dv">1</span> <-<span class="st"> </span><span class="kw">sample</span>(<span class="dv">1</span><span class="op">:</span><span class="dv">25000</span>, <span class="dt">size =</span> <span class="dv">100</span>, <span class="dt">replace =</span> <span class="ot">FALSE</span>)
sample_<span class="dv">1</span>_hours <-<span class="st"> </span>hours_population_<span class="dv">1</span>[sample_<span class="dv">1</span>]
sample_<span class="dv">1</span>_df <-<span class="st"> </span><span class="kw">data.frame</span>(<span class="dt">hours =</span> <span class="kw">round</span>(sample_<span class="dv">1</span>_hours,<span class="dv">0</span>), <span class="dt">group =</span> <span class="st">"A"</span>)
sample_<span class="dv">1</span>_df<span class="op">$</span>index <-<span class="st"> </span><span class="dv">1</span><span class="op">:</span><span class="dv">100</span>
<span class="kw">set.seed</span>(<span class="dv">321</span>)
hours_population_<span class="dv">2</span> <-<span class="st"> </span><span class="kw">rnorm</span>(<span class="dv">25000</span>, <span class="dv">25</span>, <span class="dv">6</span>)
<span class="kw">set.seed</span>(<span class="dv">125672</span>)
sample_<span class="dv">2</span> <-<span class="st"> </span><span class="kw">sample</span>(<span class="dv">1</span><span class="op">:</span><span class="dv">25000</span>, <span class="dt">size =</span> <span class="dv">100</span>, <span class="dt">replace =</span> <span class="ot">FALSE</span>)
sample_<span class="dv">2</span>_hours <-<span class="st"> </span>hours_population_<span class="dv">2</span>[sample_<span class="dv">2</span>]
sample_<span class="dv">2</span>_df <-<span class="st"> </span><span class="kw">data.frame</span>(<span class="dt">hours =</span> <span class="kw">round</span>(sample_<span class="dv">2</span>_hours,<span class="dv">0</span>), <span class="dt">group =</span> <span class="st">"B"</span>)
sample_<span class="dv">2</span>_df<span class="op">$</span>index <-<span class="st"> </span><span class="dv">1</span><span class="op">:</span><span class="dv">100</span>
<span class="kw">set.seed</span>(<span class="dv">321</span>)
hours_population_<span class="dv">3</span> <-<span class="st"> </span><span class="kw">rnorm</span>(<span class="dv">25000</span>, <span class="dv">35</span>, <span class="dv">6</span>)
<span class="kw">set.seed</span>(<span class="dv">125678</span>)
sample_<span class="dv">3</span> <-<span class="st"> </span><span class="kw">sample</span>(<span class="dv">1</span><span class="op">:</span><span class="dv">25000</span>, <span class="dt">size =</span> <span class="dv">100</span>, <span class="dt">replace =</span> <span class="ot">FALSE</span>)
sample_<span class="dv">3</span>_hours <-<span class="st"> </span>hours_population_<span class="dv">3</span>[sample_<span class="dv">3</span>]
sample_<span class="dv">3</span>_df <-<span class="st"> </span><span class="kw">data.frame</span>(<span class="dt">hours =</span> <span class="kw">round</span>(sample_<span class="dv">3</span>_hours,<span class="dv">0</span>), <span class="dt">group =</span> <span class="st">"C"</span>)
sample_<span class="dv">3</span>_df<span class="op">$</span>index <-<span class="st"> </span><span class="dv">1</span><span class="op">:</span><span class="dv">100</span>
hours_abc <-<span class="st"> </span><span class="kw">rbind</span>(sample_<span class="dv">1</span>_df, sample_<span class="dv">2</span>_df, sample_<span class="dv">3</span>_df)
<span class="kw">head</span>(hours_abc)</code></pre></div>
<div data-pagedtable="false">
<script data-pagedtable-source type="application/json">
{"columns":[{"label":["hours"],"name":[1],"type":["dbl"],"align":["right"]},{"label":["group"],"name":[2],"type":["fctr"],"align":["left"]},{"label":["index"],"name":[3],"type":["int"],"align":["right"]}],"data":[{"1":"13","2":"A","3":"1"},{"1":"14","2":"A","3":"2"},{"1":"16","2":"A","3":"3"},{"1":"17","2":"A","3":"4"},{"1":"9","2":"A","3":"5"},{"1":"13","2":"A","3":"6"}],"options":{"columns":{"min":{},"max":[10]},"rows":{"min":[10],"max":[10]},"pages":{}}}
</script>
</div>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">tail</span>(hours_abc)</code></pre></div>
<div data-pagedtable="false">
<script data-pagedtable-source type="application/json">
{"columns":[{"label":["hours"],"name":[1],"type":["dbl"],"align":["right"]},{"label":["group"],"name":[2],"type":["fctr"],"align":["left"]},{"label":["index"],"name":[3],"type":["int"],"align":["right"]}],"data":[{"1":"32","2":"C","3":"95"},{"1":"37","2":"C","3":"96"},{"1":"34","2":"C","3":"97"},{"1":"37","2":"C","3":"98"},{"1":"30","2":"C","3":"99"},{"1":"28","2":"C","3":"100"}],"options":{"columns":{"min":{},"max":[10]},"rows":{"min":[10],"max":[10]},"pages":{}}}
</script>
</div>
<p>The null hypothesis, typically, is that all means are equal (non-directional hypothesis). Hence, in our case:</p>
<p><span class="math display">\[H_0: \mu_1 = \mu_2 = \mu_3\]</span></p>
<p>The alternative hypothesis is simply that the means are not all equal, i.e.,</p>
<p><span class="math display">\[H_1: \textrm{Means are not all equal}\]</span></p>
<p>If you wanted to put this in mathematical notation, you could also write:</p>
<p><span class="math display">\[H_1: \exists {i,j}: {\mu_i \ne \mu_j} \]</span></p>
<p>To get a first impression if there are any differences in listening times across the experimental groups, we use the <code>describeBy(...)</code> function from the <code>psych</code> package:</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">library</span>(psych)
<span class="kw">describeBy</span>(hours_abc<span class="op">$</span>hours, hours_abc<span class="op">$</span>group) <span class="co">#inspect data</span></code></pre></div>
<pre><code>##
## Descriptive statistics by group
## group: A
## vars n mean sd median trimmed mad min max range skew kurtosis
## X1 1 100 15.47 5.45 15 15.22 4.45 2 36 34 0.64 1.34
## se
## X1 0.55
## --------------------------------------------------------
## group: B
## vars n mean sd median trimmed mad min max range skew kurtosis
## X1 1 100 24.88 6.02 25 24.89 5.93 8 40 32 -0.09 0.04
## se
## X1 0.6
## --------------------------------------------------------
## group: C
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 100 33.66 6 33.5 33.9 5.93 18 44 26 -0.25 -0.63 0.6</code></pre>
<p>In addition, you should visualize the data using appropriate plots:</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="co">#Plot of means</span>
<span class="kw">library</span>(plyr)
<span class="kw">library</span>(ggplot2)
<span class="kw">ggplot</span>(hours_abc, <span class="kw">aes</span>(group, hours)) <span class="op">+</span><span class="st"> </span>
<span class="st"> </span><span class="kw">stat_summary</span>(<span class="dt">fun.y =</span> mean, <span class="dt">geom =</span> <span class="st">"bar"</span>, <span class="dt">fill =</span> <span class="st">"White"</span>, <span class="dt">colour =</span> <span class="st">"Black"</span>) <span class="op">+</span>
<span class="st"> </span><span class="kw">stat_summary</span>(<span class="dt">fun.data =</span> mean_cl_normal, <span class="dt">geom =</span> <span class="st">"pointrange"</span>) <span class="op">+</span><span class="st"> </span>
<span class="st"> </span><span class="kw">labs</span>(<span class="dt">x =</span> <span class="st">"Experimental group "</span>, <span class="dt">y =</span> <span class="st">"Listening time"</span>) <span class="op">+</span><span class="st"> </span>
<span class="st"> </span><span class="kw">theme_bw</span>()</code></pre></div>
<div class="figure" style="text-align: center"><span id="fig:unnamed-chunk-5"></span>
<img src="09-Anova_files/figure-html/unnamed-chunk-5-1.png" alt="Plot of means" width="672" />
<p class="caption">
Figure 1: Plot of means
</p>
</div>
<p>Note that ANOVA is an omnibus test, which means that we test for an overall difference between groups. Hence, the test will only tell you if the group means are different, but it won’t tell you exactly which groups are different from another.</p>
<p>So why don’t we then just conduct a series of t-tests for all combinations of groups (i.e., A vs. B, A vs. C, B vs. C)? The reason is that if we assume each test to be independent, then there is a 5% probability of falsely rejecting the null hypothesis (Type I error) for each test. In our case:</p>
<ul>
<li>A vs. B (α = 0.05)</li>
<li>A vs. C (α = 0.05)</li>
<li>B vs. C (α = 0.05)</li>
</ul>
<p>This means that the overall probability of making a Type I error is 1-(0.95<sup>3</sup>) = 0.143, since the probability of no Type I error is 0.95 for each of the three tests. Consequently, the Type I error probability would be 14.3%, which is above the conventional standard of 5%. This is also known as the family-wise or experiment-wise error.</p>
</div>
<div id="decomposing-variance" class="section level3">
<h3><span class="header-section-number">0.1.2</span> Decomposing variance</h3>
<p>The basic concept underlying ANOVA is the decomposition of the variance in the data. There are three variance components which we need to consider:</p>
<ul>
<li>We calculate how much variability there is between scores: <b>Total sum of squares (SS<sub>T</sub>)</b></li>
<li>We then calculate how much of this variability can be explained by the model we fit to the data (i.e., how much variability is due to the experimental manipulation): <b>Model sum of squares (SS<sub>M</sub>)</b></li>
<li>… and how much cannot be explained (i.e., how much variability is due to individual differences in performance): <b>Residual sum of squares (SS<sub>R</sub>)</b></li>
</ul>
<p>The following figure shows the different variance components using a generalized data matrix:</p>
<div class="figure">
<img src="https://github.com/IMSMWU/Teaching/raw/master/MRDA2017/sum_of_squares.JPG" alt="Decomposing variance" />
<p class="caption">Decomposing variance</p>
</div>
<p>The total variation is determined by the variation between the categories (due to our experimental manipulation) and the within-category variation that is due to extraneous factors (e.g., promotion of artists on a social network):</p>
<p><span class="math display">\[SS_T= SS_M+SS_R\]</span></p>
<p>To get a better feeling how this relates to our data set, we can look at the data in a slightly different way. Specifically, we can use the <code>dcast(...)</code> function from the <code>reshape2</code> package to convert the data to wide format:</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">library</span>(reshape2)
<span class="kw">dcast</span>(hours_abc, index <span class="op">~</span><span class="st"> </span>group, <span class="dt">value.var =</span> <span class="st">"hours"</span>)</code></pre></div>
<div data-pagedtable="false">
<script data-pagedtable-source type="application/json">
{"columns":[{"label":["index"],"name":[1],"type":["int"],"align":["right"]},{"label":["A"],"name":[2],"type":["dbl"],"align":["right"]},{"label":["B"],"name":[3],"type":["dbl"],"align":["right"]},{"label":["C"],"name":[4],"type":["dbl"],"align":["right"]}],"data":[{"1":"1","2":"13","3":"21","4":"30"},{"1":"2","2":"14","3":"18","4":"37"},{"1":"3","2":"16","3":"27","4":"43"},{"1":"4","2":"17","3":"30","4":"31"},{"1":"5","2":"9","3":"27","4":"35"},{"1":"6","2":"13","3":"20","4":"42"},{"1":"7","2":"16","3":"24","4":"37"},{"1":"8","2":"13","3":"28","4":"43"},{"1":"9","2":"12","3":"27","4":"29"},{"1":"10","2":"22","3":"25","4":"38"},{"1":"11","2":"19","3":"15","4":"42"},{"1":"12","2":"26","3":"25","4":"32"},{"1":"13","2":"36","3":"34","4":"26"},{"1":"14","2":"21","3":"19","4":"32"},{"1":"15","2":"16","3":"24","4":"35"},{"1":"16","2":"14","3":"28","4":"33"},{"1":"17","2":"18","3":"27","4":"43"},{"1":"18","2":"18","3":"28","4":"23"},{"1":"19","2":"13","3":"26","4":"44"},{"1":"20","2":"23","3":"23","4":"30"},{"1":"21","2":"18","3":"29","4":"32"},{"1":"22","2":"16","3":"27","4":"29"},{"1":"23","2":"17","3":"25","4":"41"},{"1":"24","2":"13","3":"26","4":"31"},{"1":"25","2":"18","3":"22","4":"41"},{"1":"26","2":"13","3":"11","4":"27"},{"1":"27","2":"17","3":"35","4":"37"},{"1":"28","2":"9","3":"26","4":"37"},{"1":"29","2":"18","3":"20","4":"30"},{"1":"30","2":"20","3":"28","4":"39"},{"1":"31","2":"14","3":"21","4":"27"},{"1":"32","2":"10","3":"25","4":"32"},{"1":"33","2":"11","3":"38","4":"40"},{"1":"34","2":"19","3":"24","4":"35"},{"1":"35","2":"11","3":"29","4":"25"},{"1":"36","2":"13","3":"23","4":"32"},{"1":"37","2":"17","3":"24","4":"39"},{"1":"38","2":"11","3":"36","4":"31"},{"1":"39","2":"18","3":"31","4":"31"},{"1":"40","2":"17","3":"27","4":"23"},{"1":"41","2":"13","3":"33","4":"39"},{"1":"42","2":"23","3":"22","4":"27"},{"1":"43","2":"18","3":"28","4":"36"},{"1":"44","2":"21","3":"23","4":"40"},{"1":"45","2":"8","3":"18","4":"37"},{"1":"46","2":"13","3":"22","4":"28"},{"1":"47","2":"18","3":"40","4":"32"},{"1":"48","2":"12","3":"22","4":"36"},{"1":"49","2":"13","3":"24","4":"25"},{"1":"50","2":"12","3":"30","4":"32"},{"1":"51","2":"16","3":"21","4":"36"},{"1":"52","2":"16","3":"26","4":"33"},{"1":"53","2":"17","3":"20","4":"33"},{"1":"54","2":"21","3":"34","4":"28"},{"1":"55","2":"11","3":"26","4":"28"},{"1":"56","2":"15","3":"30","4":"37"},{"1":"57","2":"23","3":"29","4":"34"},{"1":"58","2":"8","3":"30","4":"22"},{"1":"59","2":"12","3":"16","4":"35"},{"1":"60","2":"10","3":"22","4":"41"},{"1":"61","2":"14","3":"24","4":"29"},{"1":"62","2":"14","3":"18","4":"41"},{"1":"63","2":"31","3":"31","4":"39"},{"1":"64","2":"12","3":"17","4":"24"},{"1":"65","2":"9","3":"27","4":"29"},{"1":"66","2":"22","3":"34","4":"33"},{"1":"67","2":"14","3":"12","4":"34"},{"1":"68","2":"20","3":"30","4":"37"},{"1":"69","2":"14","3":"17","4":"30"},{"1":"70","2":"7","3":"25","4":"35"},{"1":"71","2":"17","3":"36","4":"18"},{"1":"72","2":"19","3":"13","4":"21"},{"1":"73","2":"26","3":"20","4":"39"},{"1":"74","2":"12","3":"22","4":"30"},{"1":"75","2":"24","3":"32","4":"42"},{"1":"76","2":"11","3":"15","4":"33"},{"1":"77","2":"13","3":"28","4":"43"},{"1":"78","2":"9","3":"36","4":"41"},{"1":"79","2":"17","3":"21","4":"36"},{"1":"80","2":"15","3":"20","4":"35"},{"1":"81","2":"21","3":"31","4":"34"},{"1":"82","2":"14","3":"20","4":"43"},{"1":"83","2":"8","3":"31","4":"29"},{"1":"84","2":"15","3":"26","4":"37"},{"1":"85","2":"21","3":"21","4":"27"},{"1":"86","2":"13","3":"17","4":"22"},{"1":"87","2":"11","3":"32","4":"33"},{"1":"88","2":"25","3":"21","4":"43"},{"1":"89","2":"17","3":"8","4":"40"},{"1":"90","2":"8","3":"27","4":"23"},{"1":"91","2":"16","3":"21","4":"42"},{"1":"92","2":"5","3":"29","4":"31"},{"1":"93","2":"6","3":"25","4":"42"},{"1":"94","2":"15","3":"28","4":"30"},{"1":"95","2":"16","3":"20","4":"32"},{"1":"96","2":"21","3":"23","4":"37"},{"1":"97","2":"13","3":"25","4":"34"},{"1":"98","2":"8","3":"27","4":"37"},{"1":"99","2":"23","3":"22","4":"30"},{"1":"100","2":"2","3":"17","4":"28"}],"options":{"columns":{"min":{},"max":[10]},"rows":{"min":[10],"max":[10]},"pages":{}}}
</script>
</div>
<p>In this example, X<sub>1</sub> from the generalized data matrix above would refer to the factor level “A”, X<sub>2</sub> to the level “B”, and X<sub>3</sub> to the level “C”. Y<sub>11</sub> refers to the first data point in the first row (i.e., “13”), Y<sub>12</sub> to the second data point in the first row (i.e., “21”), etc.. The grand mean (<span class="math inline">\(\overline{Y}\)</span>) and the category means (<span class="math inline">\(\overline{Y}_c\)</span>) can be easily computed:</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">mean</span>(hours_abc<span class="op">$</span>hours) <span class="co">#grand mean</span></code></pre></div>
<pre><code>## [1] 24.67</code></pre>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">by</span>(hours_abc<span class="op">$</span>hours, hours_abc<span class="op">$</span>group, mean) <span class="co">#category mean</span></code></pre></div>
<pre><code>## hours_abc$group: A
## [1] 15.47
## --------------------------------------------------------
## hours_abc$group: B
## [1] 24.88
## --------------------------------------------------------
## hours_abc$group: C
## [1] 33.66</code></pre>
<p>To see how each variance component can be derived, let’s look at the data again. The following graph shows the individual observations by experimental group:</p>
<div class="figure" style="text-align: center"><span id="fig:unnamed-chunk-8"></span>
<img src="09-Anova_files/figure-html/unnamed-chunk-8-1.png" alt="Sum of Squares" width="672" />
<p class="caption">
Figure 2: Sum of Squares
</p>
</div>
<div id="total-sum-of-squares" class="section level4">
<h4><span class="header-section-number">0.1.2.1</span> Total sum of squares</h4>
<p>To compute the total variation in the data, we consider the difference between each observation and the grand mean. The grand mean is the mean over all observations in the data set. The vertical lines in the following plot measure how far each observation is away from the grand mean:</p>
<div class="figure" style="text-align: center"><span id="fig:unnamed-chunk-9"></span>
<img src="09-Anova_files/figure-html/unnamed-chunk-9-1.png" alt="Total Sum of Squares" width="672" />
<p class="caption">
Figure 3: Total Sum of Squares
</p>
</div>
<p>The formal representation of the total sum of squares (SS<sub>T</sub>) is:</p>
<p><span class="math display">\[
SS_T= \sum_{i=1}^{N} (Y_i-\bar{Y})^2
\]</span></p>
<p>This means that we need to subtract the grand mean from each individual data point, square the difference, and sum up over all the squared differences. Thus, in our example, the total sum of squares can be calculated as:</p>
<p><span class="math display">\[
\begin{align}
SS_T =&(13−24.67)^2 + (14−24.67)^2 + … + (2−24.67)^2\\
&+(21−24.67)^2 + (18-24.67)^2 + … + (17−24.67)^2\\
&+(30−24.67)^2 + (37−24.67)^2 + … + (28−24.67)^2\\
&=26646.33
\end{align}
\]</span></p>
<p>You could also compute this in R using:</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">SST <-<span class="st"> </span><span class="kw">sum</span>((hours_abc<span class="op">$</span>hours <span class="op">-</span><span class="st"> </span><span class="kw">mean</span>(hours_abc<span class="op">$</span>hours))<span class="op">^</span><span class="dv">2</span>)
SST</code></pre></div>
<pre><code>## [1] 26646.33</code></pre>
<p>For the subsequent analyses, it is important to understand the concept behind the <b>degrees of freedom</b>. Remember that in order to estimate a population value from a sample, we need to hold something in the population constant. In ANOVA, the df are generally one less than the number of values used to calculate the SS. For example, when we estimate the population mean from a sample, we assume that the sample mean is equal to the population mean. Then, in order to estimate the population mean from the sample, all but one scores are free to vary and the remaining score needs to be the value that keeps the population mean constant. In our example, we used all 300 observations to calculate the sum of square, so the total degrees of freedom (df<sub>T</sub>) are:</p>
<span class="math display" id="eq:dfT">\[\begin{equation}
\begin{split}
df_T = N-1=300-1=299
\end{split}
\tag{1}
\end{equation}\]</span>
</div>
<div id="model-sum-of-squares" class="section level4">
<h4><span class="header-section-number">0.1.2.2</span> Model sum of squares</h4>
<p>Now we know that there are 26646.33 units of total variation in our data. Next, we compute how much of the total variation can be explained by the differences between groups (i.e., our experimental manipulation). To compute the explained variation in the data, we consider the difference between the values predicted by our model for each observation (i.e., the group mean) and the grand mean. The group mean refers to the mean value within the experimental group. The vertical lines in the following plot measure how far the predicted value for each observation (i.e., the group mean) is away from the grand mean:</p>
<div class="figure" style="text-align: center"><span id="fig:unnamed-chunk-11"></span>
<img src="09-Anova_files/figure-html/unnamed-chunk-11-1.png" alt="Model Sum of Squares" width="672" />
<p class="caption">
Figure 4: Model Sum of Squares
</p>
</div>
<p>The formal representation of the model sum of squares (SS<sub>M</sub>) is:</p>
<p><span class="math display">\[
SS_M= \sum_{j=1}^{c} n_j(\bar{Y}_j-\bar{Y})^2
\]</span></p>
<p>where c denotes the number of categories (experimental groups). This means that we need to subtract the grand mean from each group mean, square the difference, and sum up over all the squared differences. Thus, in our example, the model sum of squares can be calculated as:</p>
<p><span class="math display">\[
\begin{align}
SS_M &= 100*(15.47−24.67)^2 + 100*(24.88−24.67)^2 + 100*(33.66−24.67)^2 \\
&= 16550.42
\end{align}
\]</span></p>
<p>You could also compute this manually in R using:</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">SSM <-<span class="st"> </span><span class="kw">sum</span>(<span class="dv">100</span><span class="op">*</span>(<span class="kw">by</span>(hours_abc<span class="op">$</span>hours, hours_abc<span class="op">$</span>group, mean) <span class="op">-</span><span class="st"> </span><span class="kw">mean</span>(hours_abc<span class="op">$</span>hours))<span class="op">^</span><span class="dv">2</span>)
SSM</code></pre></div>
<pre><code>## [1] 16550.42</code></pre>
<p>In this case, we used the three group means to calculate the sum of squares, so the model degrees of freedom (df<sub>M</sub>) are:</p>
<p><span class="math display">\[
df_M= c-1=3-1=2
\]</span></p>
</div>
<div id="residual-sum-of-squares" class="section level4">
<h4><span class="header-section-number">0.1.2.3</span> Residual sum of squares</h4>
<p>Lastly, we calculate the amount of variation that cannot be explained by our model. In ANOVA, this is the sum of squared distances between what the model predicts for each data point (i.e., the group means) and the observed values. In other words, this refers to the amount of variation that is caused by extraneous factors, such as differences between product characteristics of the products in the different experimental groups. The vertical lines in the following plot measure how far each observation is away from the group mean:</p>
<div class="figure" style="text-align: center"><span id="fig:unnamed-chunk-13"></span>
<img src="09-Anova_files/figure-html/unnamed-chunk-13-1.png" alt="Residual Sum of Squares" width="672" />
<p class="caption">
Figure 5: Residual Sum of Squares
</p>
</div>
<p>The formal representation of the residual sum of squares (SS<sub>R</sub>) is:</p>
<p><span class="math display">\[
SS_R= \sum_{j=1}^{c} \sum_{i=1}^{n} ({Y}_{ij}-\bar{Y}_{j})^2
\]</span></p>
<p>This means that we need to subtract the group mean from each individual observation, square the difference, and sum up over all the squared differences. Thus, in our example, the model sum of squares can be calculated as:</p>
<p><span class="math display">\[
\begin{align}
SS_R =& (13−15.47)^2 + (14−15.47)^2 + … + (2−15.47)^2 \\
+&(21−24.88)^2 + (18−24.88)^2 + … + (17−24.88)^2 \\
+& (30−33.66)^2 + (37−33.66)^2 + … + (28−33.66)^2 \\
=& 10095.91
\end{align}
\]</span></p>
<p>You could also compute this in R using:</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">SSR <-<span class="st"> </span><span class="kw">sum</span>((hours_abc<span class="op">$</span>hours <span class="op">-</span><span class="st"> </span><span class="kw">rep</span>(<span class="kw">by</span>(hours_abc<span class="op">$</span>hours, hours_abc<span class="op">$</span>group, mean), <span class="dt">each =</span> <span class="dv">100</span>))<span class="op">^</span><span class="dv">2</span>)
SSR</code></pre></div>
<pre><code>## [1] 10095.91</code></pre>
<p>In this case, we used the 10 values for each of the SS for each group, so the residual degrees of freedom (df<sub>R</sub>) are:</p>
<p><span class="math display">\[
\begin{align}
df_R=& (n_1-1)+(n_2-1)+(n_3-1) \\
=&(100-1)+(100-1)+(100-1)=297
\end{align}
\]</span></p>
</div>
<div id="effect-strength" class="section level4">
<h4><span class="header-section-number">0.1.2.4</span> Effect strength</h4>
<p>Once you have computed the different sum of squares, you can investigate the effect strength. <span class="math inline">\(\eta^2\)</span> is a measure of the variation in Y that is explained by X:</p>
<p><span class="math display">\[
\eta^2= \frac{SS_M}{SS_T}=\frac{16550.42}{26646.33}=0.621
\]</span></p>
<p>To compute this in R:</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">eta <-<span class="st"> </span>SSM<span class="op">/</span>SST
eta</code></pre></div>
<pre><code>## [1] 0.6211144</code></pre>
<p>The statistic can only take values between 0 and 1. It is equal to 0 when all the category means are equal, indicating that X has no effect on Y. In contrast, it has a value of 1 when there is no variability within each category of X but there is some variability between categories.</p>
</div>
<div id="test-of-significance" class="section level4">
<h4><span class="header-section-number">0.1.2.5</span> Test of significance</h4>
<p>How can we determine whether the effect of X on Y is significant?</p>
<ul>
<li>First, we calculate the fit of the most basic model (i.e., the grand mean)</li>
<li>Then, we calculate the fit of the “best” model (i.e., the group means)</li>
<li>A good model should fit the data significantly better than the basic model</li>
<li>The F-statistic or F-ratio compares the amount of systematic variance in the data to the amount of unsystematic variance</li>
</ul>
<p>The F-statistic uses the ratio of mean square related to X (explained variation) and the mean square related to the error (unexplained variation):</p>
<p style="text-align:center;">
<span class="math inline">\(\frac{SS_M}{SS_R}\)</span> <br>
</p>
<p>However, since these are summed values, their magnitude is influenced by the number of scores that were summed. For example, to calculate SS<sub>M</sub> we only used the sum of 3 values (the group means), while we used 30 and 27 values to calculate SS<sub>T</sub> and SS<sub>R</sub>, respectively. Thus, we calculate the average sum of squares (“mean square”) to compare the average amount of systematic vs. unsystematic variation by dividing the SS values by the degrees of freedom associated with the respective statistic.</p>
<p>Mean square due to X:</p>
<p><span class="math display">\[
MS_M= \frac{SS_M}{df_M}=\frac{SS_M}{c-1}=\frac{16550.42}{(3-1)}
\]</span></p>
<p>Mean square due to error:</p>
<p><span class="math display">\[
MS_R= \frac{SS_R}{df_R}=\frac{SS_R}{N-c}=\frac{10095.91}{(300-3)}
\]</span></p>
<p>Now, we compare the amount of variability explained by the model (experiment), to the error in the model (variation due to extraneous variables). If the model explains more variability than it can’t explain, then the experimental manipulation has had a significant effect on the outcome (DV). The F-radio can be derived as follows:</p>
<p><span class="math display">\[
F= \frac{MS_M}{MS_R}=\frac{SS_R}{N-c}=\frac{\frac{16550.42}{(3-1)}}{\frac{10095.91}{(300-3)}}=243.44
\]</span></p>
<p>You can easily compute this in R:</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">f_ratio <-<span class="st"> </span>(SSM<span class="op">/</span><span class="dv">2</span>)<span class="op">/</span>(SSR<span class="op">/</span><span class="dv">297</span>)
f_ratio</code></pre></div>
<pre><code>## [1] 243.4389</code></pre>
<p>This statistic follows the F distribution with (m = c – 1) and (n = N – c) degrees of freedom. This means that, like the <span class="math inline">\(\chi^2\)</span> distribution, the shape of the F-distribution depends on the degrees of freedom. In this case, the shape depends on the degrees of freedom associated with the numerator and denominator used to compute the F-ratio. The following figure shows the shape of the F-distribution for different degrees of freedom:</p>
<div class="figure">
<embed src="%22./fdistributions.png%22" />
<p class="caption">The F distribution</p>
</div>
<p>The outcome of the test is one of the following:</p>
<ul>
<li>If the null hypothesis of equal category means is not rejected, then the independent variable does not have a significant effect on the dependent variable</li>
<li>If the null hypothesis is rejected, then the effect of the independent variable is significant</li>
</ul>
<p>For 2 and 297 degrees of freedom, the critical value of F is 3.026 for α=0.05. As usual, you can either look up these values in a table or use the appropriate function in R:</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">f_crit <-<span class="st"> </span><span class="kw">qf</span>(.<span class="dv">95</span>, <span class="dt">df1 =</span> <span class="dv">2</span>, <span class="dt">df2 =</span> <span class="dv">297</span>) <span class="co">#critical value</span>
f_crit </code></pre></div>
<pre><code>## [1] 3.026153</code></pre>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">f_ratio <span class="op">></span><span class="st"> </span>f_crit <span class="co">#test if calculated test statistic is larger than critical value</span></code></pre></div>
<pre><code>## [1] TRUE</code></pre>
<p>The output tells us that the calculated test statistic exceeds the critical value. We can also show the test result visually:</p>
<div class="figure">
<embed src="%22./ftest.png%22" />
<p class="caption">Visual depiction of the test result</p>
</div>
<p>Thus, we conclude that because F<sub>CAL</sub> = 243.44 > F<sub>CR</sub> = 3.03, H<sub>0</sub> is rejected!</p>
<p>Interpretation: one or more of the differences between means are statistically significant.</p>
<p>Reporting: There was a significant effect of promotion on sales levels, F(2,297) = 243.44, p < 0.05, <span class="math inline">\(\eta^2\)</span> = 0.62.</p>
<p>Remember: This doesn’t tell us where the differences between groups lie. To find out which group means exactly differ, we need to use post-hoc procedures (see below).</p>
<p><em>You don’t have to compute these statistics manually!</em> Luckily, there is a function for ANOVA in R, which does the above calculations for you as we will see in the next section.</p>
</div>
</div>
<div id="one-way-anova" class="section level3">
<h3><span class="header-section-number">0.1.3</span> One-way ANOVA</h3>
<div id="basic-anova" class="section level4">
<h4><span class="header-section-number">0.1.3.1</span> Basic ANOVA</h4>
<p>As already indicated, one-way ANOVA is used when there is only one categorical variable (factor). Before conducting ANOVA, you need to check if the assumptions of the test are fulfilled. The assumptions of ANOVA are discussed in the following sections.</p>
<div id="independence-of-observations" class="section level5 unnumbered">
<h5>Independence of observations</h5>
<p>The observations in the groups should be independent. Because we randomly assigned the listeners to the experimental conditions, this assumption can be assumed to be met.</p>
</div>
<div id="distributional-assumptions" class="section level5 unnumbered">
<h5>Distributional assumptions</h5>
<p>ANOVA is relatively immune to violations to the normality assumption when sample sizes are large due to the Central Limit Theorem. However, if your sample is small (i.e., n < 30 per group) you may nevertheless want to check the normality of your data, e.g., by using the Shapiro-Wilk test or QQ-Plot. In our example, we have 100 observations in each group which is plenty but let’s create another example with only 10 observations in each group. In the latter case we cannot rely on the Central Limit Theorem and we should test the normality of our data. This can be done using the Shapiro-Wilk Test, which has the Null Hypothesis that the data is normally distributed. Hence, an insignificant test results means that the data can be assumed to be approximately normally distributed:</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">set.seed</span>(<span class="dv">321</span>)
hours_fewobs <-<span class="st"> </span><span class="kw">data.frame</span>(<span class="dt">hours =</span> <span class="kw">c</span>(<span class="kw">rnorm</span>(<span class="dv">10</span>, <span class="dv">20</span>, <span class="dv">5</span>), <span class="kw">rnorm</span>(<span class="dv">10</span>, <span class="dv">40</span>, <span class="dv">5</span>), <span class="kw">rnorm</span>(<span class="dv">10</span>, <span class="dv">60</span>, <span class="dv">5</span>)),
<span class="dt">group =</span> <span class="kw">c</span>(<span class="kw">rep</span>(<span class="st">'A'</span>, <span class="dv">10</span>), <span class="kw">rep</span>(<span class="st">'B'</span>, <span class="dv">10</span>), <span class="kw">rep</span>(<span class="st">'C'</span>, <span class="dv">10</span>)))
<span class="kw">by</span>(hours_fewobs<span class="op">$</span>hours, hours_fewobs<span class="op">$</span>group, shapiro.test)</code></pre></div>
<pre><code>## hours_fewobs$group: A
##
## Shapiro-Wilk normality test
##
## data: dd[x, ]
## W = 0.91625, p-value = 0.3267
##
## --------------------------------------------------------
## hours_fewobs$group: B
##
## Shapiro-Wilk normality test
##
## data: dd[x, ]
## W = 0.91486, p-value = 0.3161
##
## --------------------------------------------------------
## hours_fewobs$group: C
##
## Shapiro-Wilk normality test
##
## data: dd[x, ]
## W = 0.9595, p-value = 0.7801</code></pre>
<p>Since the test result is insignificant for all groups, we can conclude that the data approximately follow a normal distribution.</p>
<p>We could also test the distributional assumptions visually using a Q-Q plot (i.e., quantile-quantile plot). This plot can be used to assess if a set of data plausibly came from some theoretical distribution such as the Normal distribution. Since this is just a visual check, it is somewhat subjective. But it may help us to judge if our assumption is plausible, and if not, which data points contribute to the violation. A Q-Q plot is a scatterplot created by plotting two sets of quantiles against one another. If both sets of quantiles came from the same distribution, we should see the points forming a line that’s roughly straight. In other words, Q-Q plots take your sample data, sort it in ascending order, and then plot them versus quantiles calculated from a theoretical distribution. Quantiles are often referred to as “percentiles” and refer to the points in your data below which a certain proportion of your data fall. Recall, for example, the standard Normal distribution with a mean of 0 and a standard deviation of 1. Since the 50th percentile (or 0.5 quantile) is 0, half the data lie below 0. The 95th percentile (or 0.95 quantile), is about 1.64, which means that 95 percent of the data lie below 1.64. The 97.5th quantile is about 1.96, which means that 97.5% of the data lie below 1.96. In the Q-Q plot, the number of quantiles is selected to match the size of your sample data.</p>
<p>To create the Q-Q plot for the normal distribution, you may use the <code>qqnorm()</code> function, which takes the data to be tested as an argument. Using the <code>qqline()</code> function subsequently on the data creates the line on which the data points should fall based on the theoretical quantiles. If the individual data points deviate a lot from this line, it means that the data is not likely to follow a normal distribution.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">qqnorm</span>(hours_fewobs[hours_fewobs<span class="op">$</span>group<span class="op">==</span><span class="st">"A"</span>,]<span class="op">$</span>hours)
<span class="kw">qqline</span>(hours_fewobs[hours_fewobs<span class="op">$</span>group<span class="op">==</span><span class="st">"A"</span>,]<span class="op">$</span>hours)</code></pre></div>
<div class="figure" style="text-align: center"><span id="fig:unnamed-chunk-21"></span>
<img src="09-Anova_files/figure-html/unnamed-chunk-21-1.png" alt="Q-Q plot 1" width="672" />
<p class="caption">
Figure 6: Q-Q plot 1
</p>
</div>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">qqnorm</span>(hours_fewobs[hours_fewobs<span class="op">$</span>group<span class="op">==</span><span class="st">"B"</span>,]<span class="op">$</span>hours)
<span class="kw">qqline</span>(hours_fewobs[hours_fewobs<span class="op">$</span>group<span class="op">==</span><span class="st">"B"</span>,]<span class="op">$</span>hours)</code></pre></div>
<div class="figure" style="text-align: center"><span id="fig:unnamed-chunk-21"></span>
<img src="09-Anova_files/figure-html/unnamed-chunk-21-2.png" alt="Q-Q plot 2" width="672" />
<p class="caption">
Figure 6: Q-Q plot 2
</p>
</div>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">qqnorm</span>(hours_fewobs[hours_fewobs<span class="op">$</span>group<span class="op">==</span><span class="st">"C"</span>,]<span class="op">$</span>hours)
<span class="kw">qqline</span>(hours_fewobs[hours_fewobs<span class="op">$</span>group<span class="op">==</span><span class="st">"C"</span>,]<span class="op">$</span>hours)</code></pre></div>
<div class="figure" style="text-align: center"><span id="fig:unnamed-chunk-21"></span>
<img src="09-Anova_files/figure-html/unnamed-chunk-21-3.png" alt="Q-Q plot 3" width="672" />
<p class="caption">
Figure 6: Q-Q plot 3
</p>
</div>
<p>The Q-Q plots suggest an approximately Normal distribution. If the assumption had been violated, you might consider transforming your data or resort to a non-parametric test.</p>
</div>
<div id="homogeneity-of-variance" class="section level5 unnumbered">
<h5>Homogeneity of variance</h5>
<p>Let’s return to our original dataset with 100 observations in each group for the rest of the analysis.</p>
<p>You can test the homogeneity of variances in R using Levene’s test:</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">library</span>(car)
<span class="kw">leveneTest</span>(hours <span class="op">~</span><span class="st"> </span>group, <span class="dt">data =</span> hours_abc, <span class="dt">center =</span> mean)</code></pre></div>
<pre><code>## Levene's Test for Homogeneity of Variance (center = mean)
## Df F value Pr(>F)
## group 2 1.2173 0.2975
## 297</code></pre>
<p>The null hypothesis of the test is that the group variances are equal. Thus, if the test result is significant it means that the variances are not equal. If we cannot reject the null hypothesis (i.e., the group variances are not significantly different), we can proceed with the ANOVA as follows:</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">aov <-<span class="st"> </span><span class="kw">aov</span>(hours <span class="op">~</span><span class="st"> </span>group, <span class="dt">data =</span> hours_abc)
<span class="kw">summary</span>(aov)</code></pre></div>
<pre><code>## Df Sum Sq Mean Sq F value Pr(>F)
## group 2 16550 8275 243.4 <2e-16 ***
## Residuals 297 10096 34
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1</code></pre>
<p>You can see that the p-value is smaller than 0.05. This means that, if there really was no difference between the population means (i.e., the Null hypothesis was true), the probability of the observed differences (or larger differences) is less than 5%.</p>
<p>To compute η<sup>2</sup> from the output, we can extract the relevant sum of squares as follows</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">summary</span>(aov)[[<span class="dv">1</span>]]<span class="op">$</span><span class="st">'Sum Sq'</span>[<span class="dv">1</span>]<span class="op">/</span>(<span class="kw">summary</span>(aov)[[<span class="dv">1</span>]]<span class="op">$</span><span class="st">'Sum Sq'</span>[<span class="dv">1</span>] <span class="op">+</span><span class="st"> </span><span class="kw">summary</span>(aov)[[<span class="dv">1</span>]]<span class="op">$</span><span class="st">'Sum Sq'</span>[<span class="dv">2</span>])</code></pre></div>
<pre><code>## [1] 0.6211144</code></pre>
<p>You can see that the results match the results from our manual computation above (<span class="math inline">\(\eta^2 =\)</span> 0.62).</p>
<p>The <code>aov()</code> function also automatically generates some plots that you can use to judge if the model assumptions are met. We will inspect two of the plots here.</p>
<p>We will use the first plot to inspect if the residual variances are equal across the experimental groups:</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">plot</span>(aov,<span class="dv">1</span>)</code></pre></div>
<p><img src="09-Anova_files/figure-html/unnamed-chunk-25-1.png" width="672" /></p>
<p>Generally, the residual variance (i.e., the range of values on the y-axis) should be the same for different levels of our independent variable. The plot shows, that there are some slight differences. Notably, the range of residuals is higher in group “B” than in group “C”. However, the differences are not that large and since the Levene’s test could not reject the Null of equal variances, we conclude that the variances are similar enough in this case.</p>
<p>The second plot can be used to test the assumption that the residuals are approximately normally distributed. We use a Q-Q plot to test this assumption:</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">plot</span>(aov,<span class="dv">2</span>)</code></pre></div>
<p><img src="09-Anova_files/figure-html/unnamed-chunk-26-1.png" width="672" /></p>
<p>The plot suggests that, the residuals are approximately normally distributed. We could also test this by extracting the residuals from the anova output using the <code>resid()</code> function and using the Shapiro-Wilk test:</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">shapiro.test</span>(<span class="kw">resid</span>(aov))</code></pre></div>
<pre><code>##
## Shapiro-Wilk normality test
##
## data: resid(aov)
## W = 0.99638, p-value = 0.7289</code></pre>
<p>Confirming the impression from the Q-Q plot, we cannot reject the Null that the residuals are approximately normally distributed.</p>
<p>Note that if Levene’s test would have been significant (i.e., variances are not equal), we would have needed to either resort to non-parametric tests (see below), or compute the Welch’s F-ratio instead:</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">oneway.test</span>(hours <span class="op">~</span><span class="st"> </span>group, hours_abc)</code></pre></div>
<pre><code>##
## One-way analysis of means (not assuming equal variances)
##
## data: hours and group
## F = 252.16, num df = 2.00, denom df = 197.56, p-value < 2.2e-16</code></pre>
<p>You can see that the results are fairly similar, since the variances turned out to be fairly equal across groups.</p>
</div>
</div>
<div id="post-hoc-tests" class="section level4">
<h4><span class="header-section-number">0.1.3.2</span> Post-hoc tests</h4>
<p>Provided that significant differences were detected by the overall ANOVA you can find out which group means are different using post hoc procedures. Post hoc procedures are designed to conduct pairwise comparisons of all different combinations of the treatment groups by correcting the level of significance for each test such that the overall Type I error rate (α) across all comparisons remains at 0.05.</p>
<p>In other words, we rejected H<sub>0</sub>: μ<sub>1</sub>= μ<sub>2</sub>= μ<sub>3</sub>, and now we would like to test:</p>
<p>Test1:</p>
<p><span class="math display">\[H_0: \mu_1 = \mu_2\]</span></p>
<p>Test2:</p>
<p><span class="math display">\[H_0: \mu_1 = \mu_3\]</span></p>
<p>Test3:</p>
<p><span class="math display">\[H_0: \mu_2 = \mu_3\]</span></p>
<p>There are several post hoc procedures available to choose from. In this tutorial, we will cover Bonferroni and Tukey’s HSD (“honest significant differences”). Both tests control for family-wise error. Bonferroni tends to have more power when the number of comparisons is small, whereas Tukey’ HSDs is better when testing large numbers of means.</p>
<div id="bonferroni" class="section level5">
<h5><span class="header-section-number">0.1.3.2.1</span> Bonferroni</h5>
<p>One of the most popular (and easiest) methods to correct for the family-wise error rate is to conduct the individual t-tests and divide α by the number of comparisons („k“):</p>
<p><span class="math display">\[
p_{CR}= \frac{\alpha}{k}
\]</span></p>
<p>In our example with three groups:</p>
<p><span class="math display">\[p_{CR}= \frac{0.05}{3}=0.017\]</span></p>
<p>Thus, the “corrected” critical p-value is now 0.017 instead of 0.05 (i.e., the critical t value is higher). You can implement the Bonferroni procedure in R using:</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">bonferroni <-<span class="st"> </span><span class="kw">pairwise.t.test</span>(hours_abc<span class="op">$</span>hours, hours_abc<span class="op">$</span>group, <span class="dt">data =</span> hours_abc, <span class="dt">p.adjust.method =</span> <span class="st">"bonferroni"</span>)
bonferroni</code></pre></div>
<pre><code>##
## Pairwise comparisons using t tests with pooled SD
##
## data: hours_abc$hours and hours_abc$group
##
## A B
## B <2e-16 -
## C <2e-16 <2e-16
##
## P value adjustment method: bonferroni</code></pre>
<p>In the output, you will get the corrected p-values for the individual tests. In our example, we can reject H<sub>0</sub> of equal means for all three tests, since p < 0.05 for all combinations of groups.</p>
<p>Note the difference between the results from the post-hoc test compared to individual t-tests. For example, when we test the “B” vs. “C” groups, the result from a t-test would be:</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">data_subset <-<span class="st"> </span><span class="kw">subset</span>(hours_abc, group <span class="op">!=</span><span class="st"> "A"</span>)
ttest <-<span class="st"> </span><span class="kw">t.test</span>(hours <span class="op">~</span><span class="st"> </span>group, <span class="dt">data =</span> data_subset, <span class="dt">var.equal=</span> <span class="ot">TRUE</span>)
ttest</code></pre></div>
<pre><code>##
## Two Sample t-test
##
## data: hours by group
## t = -10.331, df = 198, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -10.45601 -7.10399
## sample estimates:
## mean in group B mean in group C
## 24.88 33.66</code></pre>
<p>Usually the p-value is lower in the t-test, reflecting the fact that the family-wise error is not corrected (i.e., the test is less conservative). In this case the p-value is extremely small in both cases and thus indistinguishable.</p>
</div>
<div id="tukeys-hsd" class="section level5">
<h5><span class="header-section-number">0.1.3.2.2</span> Tukey’s HSD</h5>
<p>Tukey’s HSD also compares all possible pairs of means (two-by-two combinations; i.e., like a t-test, except that it corrects for family-wise error rate).</p>
<p>Test statistic:</p>
<span class="math display" id="eq:tukey">\[\begin{equation}
\begin{split}
HSD= q\sqrt{\frac{MS_R}{n_c}}
\end{split}
\tag{2}
\end{equation}\]</span>
<p>where:</p>
<ul>
<li>q = value from studentized range table (see e.g., <a href="http://www.real-statistics.com/statistics-tables/studentized-range-q-table/" target="_blank">here</a>)</li>
<li>MS<sub>R</sub> = Mean Square Error from ANOVA</li>
<li>n<sub>c</sub> = number of observations per group</li>
<li>Decision: Reject H<sub>0</sub> if</li>
</ul>
<p><span class="math display">\[|\bar{Y}_i-\bar{Y}_j | > HSD\]</span></p>
<p>The value from the studentized range table can be obtained using the <code>qtukey()</code> function.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">q <-<span class="st"> </span><span class="kw">qtukey</span>(<span class="fl">0.95</span>, <span class="dt">nm =</span> <span class="dv">3</span>, <span class="dt">df =</span> <span class="dv">297</span>)
q</code></pre></div>
<pre><code>## [1] 3.331215</code></pre>
<p>Hence:</p>
<p><span class="math display">\[HSD= 3.33\sqrt{\frac{33.99}{100}}=1.94\]</span></p>
<p>Or, in R:</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">hsd <-<span class="st"> </span>q <span class="op">*</span><span class="st"> </span><span class="kw">sqrt</span>(<span class="kw">summary</span>(aov)[[<span class="dv">1</span>]]<span class="op">$</span><span class="st">'Mean Sq'</span>[<span class="dv">2</span>]<span class="op">/</span><span class="dv">100</span>)
hsd</code></pre></div>
<pre><code>## [1] 1.942214</code></pre>
<p>Since all mean differences between groups are larger than 1.906, we can reject the null hypothesis for all individual tests, confirming the results from the Bonferroni test. To compute Tukey’s HSD, we can use the appropriate function from the <code>multcomp</code> package.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">library</span>(multcomp)
tukeys <-<span class="st"> </span><span class="kw">glht</span>(aov, <span class="dt">linfct =</span> <span class="kw">mcp</span>(<span class="dt">group =</span> <span class="st">"Tukey"</span>))
<span class="kw">summary</span>(tukeys)</code></pre></div>
<pre><code>##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: aov(formula = hours ~ group, data = hours_abc)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## B - A == 0 9.4100 0.8245 11.41 <2e-16 ***
## C - A == 0 18.1900 0.8245 22.06 <2e-16 ***
## C - B == 0 8.7800 0.8245 10.65 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)</code></pre>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">confint</span>(tukeys)</code></pre></div>
<pre><code>##
## Simultaneous Confidence Intervals
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: aov(formula = hours ~ group, data = hours_abc)
##
## Quantile = 2.3558
## 95% family-wise confidence level
##
##
## Linear Hypotheses:
## Estimate lwr upr
## B - A == 0 9.4100 7.4676 11.3524
## C - A == 0 18.1900 16.2476 20.1324
## C - B == 0 8.7800 6.8376 10.7224</code></pre>
<p>We may also plot the result for the mean differences incl. their confidence intervals:</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">plot</span>(tukeys)</code></pre></div>
<div class="figure" style="text-align: center"><span id="fig:unnamed-chunk-34"></span>
<img src="09-Anova_files/figure-html/unnamed-chunk-34-1.png" alt="Tukey's HSD" width="672" />
<p class="caption">
Figure 7: Tukey’s HSD
</p>
</div>
<p>You can see that the CIs do not cross zero, which means that the true difference between group means is unlikely zero.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">mean1 <-<span class="st"> </span><span class="kw">mean</span>(hours_abc[hours_abc<span class="op">$</span>group<span class="op">==</span><span class="st">"A"</span>,<span class="st">"hours"</span>]) <span class="co">#mean group "A"</span>
mean1</code></pre></div>
<pre><code>## [1] 15.47</code></pre>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">mean2 <-<span class="st"> </span><span class="kw">mean</span>(hours_abc[hours_abc<span class="op">$</span>group<span class="op">==</span><span class="st">"B"</span>,<span class="st">"hours"</span>]) <span class="co">#mean group "B"</span>
mean2</code></pre></div>
<pre><code>## [1] 24.88</code></pre>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">mean3 <-<span class="st"> </span><span class="kw">mean</span>(hours_abc[hours_abc<span class="op">$</span>group<span class="op">==</span><span class="st">"C"</span>,<span class="st">"hours"</span>]) <span class="co">#mean group "C"</span>
mean3</code></pre></div>
<pre><code>## [1] 33.66</code></pre>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="co">#CI high vs. medium</span>
mean_diff_high_med <-<span class="st"> </span>mean2<span class="op">-</span>mean1
mean_diff_high_med</code></pre></div>
<pre><code>## [1] 9.41</code></pre>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">ci_med_high_lower <-<span class="st"> </span>mean_diff_high_med<span class="op">-</span>hsd
ci_med_high_upper <-<span class="st"> </span>mean_diff_high_med<span class="op">+</span>hsd
ci_med_high_lower</code></pre></div>
<pre><code>## [1] 7.467786</code></pre>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">ci_med_high_upper</code></pre></div>
<pre><code>## [1] 11.35221</code></pre>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="co">#CI high vs.low</span>
mean_diff_high_low <-<span class="st"> </span>mean3<span class="op">-</span>mean1
mean_diff_high_low</code></pre></div>
<pre><code>## [1] 18.19</code></pre>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">ci_low_high_lower <-<span class="st"> </span>mean_diff_high_low<span class="op">-</span>hsd
ci_low_high_upper <-<span class="st"> </span>mean_diff_high_low<span class="op">+</span>hsd
ci_low_high_lower</code></pre></div>
<pre><code>## [1] 16.24779</code></pre>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">ci_low_high_upper</code></pre></div>
<pre><code>## [1] 20.13221</code></pre>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="co">#CI medium vs.low</span>
mean_diff_med_low <-<span class="st"> </span>mean3<span class="op">-</span>mean2
mean_diff_med_low</code></pre></div>
<pre><code>## [1] 8.78</code></pre>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">ci_low_med_lower <-<span class="st"> </span>mean_diff_med_low<span class="op">-</span>hsd
ci_low_med_upper <-<span class="st"> </span>mean_diff_med_low<span class="op">+</span>hsd
ci_low_med_lower</code></pre></div>
<pre><code>## [1] 6.837786</code></pre>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">ci_low_med_upper</code></pre></div>
<pre><code>## [1] 10.72221</code></pre>
<p>Reporting of post hoc results:</p>
<p>The post hoc tests based on Bonferroni and Tukey’s HSD revealed that people listened to music significantly more when:</p>
<ul>
<li>they had access to user created playlists vs. those who did not,</li>
<li>they got recommendations vs. those who did not. This is true for both the control group “A” as well as treatment “B”.</li>
</ul>
<p><strong>The following video summarizes how to conduct a one-way ANOVA in R</strong></p>
<div align="center">
<iframe width="560" height="315" src="https://www.youtube.com/embed/vl32om-3KpY" frameborder="0" allowfullscreen>
</iframe>
</div>
</div>
</div>
</div>
</div>
</section>
</div>
</div>
</div>
</div>
</div>
<script src="libs/gitbook-2.6.7/js/app.min.js"></script>
<script src="libs/gitbook-2.6.7/js/lunr.js"></script>
<script src="libs/gitbook-2.6.7/js/plugin-search.js"></script>
<script src="libs/gitbook-2.6.7/js/plugin-sharing.js"></script>
<script src="libs/gitbook-2.6.7/js/plugin-fontsettings.js"></script>
<script src="libs/gitbook-2.6.7/js/plugin-bookdown.js"></script>
<script src="libs/gitbook-2.6.7/js/jquery.highlight.js"></script>
<script>
gitbook.require(["gitbook"], function(gitbook) {
gitbook.start({
"sharing": {
"github": false,
"facebook": true,
"twitter": true,
"google": false,
"linkedin": false,
"weibo": false,
"instapaper": false,
"vk": false,
"all": ["facebook", "google", "twitter", "linkedin", "weibo", "instapaper"]
},
"fontsettings": {
"theme": "white",
"family": "sans",
"size": 2
},
"edit": {
"link": null,
"text": null
},
"history": {
"link": null,
"text": null
},
"download": null,
"toc": {
"collapse": "section"
},
"search": false
});
});
</script>
<!-- dynamically load mathjax for compatibility with self-contained -->
<script>
(function () {
var script = document.createElement("script");
script.type = "text/javascript";
var src = "";
if (src === "" || src === "true") src = "https://mathjax.rstudio.com/latest/MathJax.js?config=TeX-MML-AM_CHTML";
if (location.protocol !== "file:" && /^https?:/.test(src))
src = src.replace(/^https?:/, '');
script.src = src;
document.getElementsByTagName("head")[0].appendChild(script);
})();
</script>
</body>
</html>