forked from lindeloev/tests-as-linear
-
Notifications
You must be signed in to change notification settings - Fork 0
/
index.html
1703 lines (1556 loc) · 109 KB
/
index.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
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
<!DOCTYPE html>
<html xmlns="http://www.w3.org/1999/xhtml">
<head>
<meta charset="utf-8" />
<meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
<meta name="generator" content="pandoc" />
<title>Common statistical tests are linear models (or: how to teach stats)</title>
<script src="index_files/jquery-1.12.4/jquery.min.js"></script>
<meta name="viewport" content="width=device-width, initial-scale=1" />
<link href="index_files/bootstrap-3.3.5/css/bootstrap.min.css" rel="stylesheet" />
<script src="index_files/bootstrap-3.3.5/js/bootstrap.min.js"></script>
<script src="index_files/bootstrap-3.3.5/shim/html5shiv.min.js"></script>
<script src="index_files/bootstrap-3.3.5/shim/respond.min.js"></script>
<script src="index_files/jqueryui-1.11.4/jquery-ui.min.js"></script>
<link href="index_files/tocify-1.9.1/jquery.tocify.css" rel="stylesheet" />
<script src="index_files/tocify-1.9.1/jquery.tocify.js"></script>
<script src="index_files/navigation-1.1/tabsets.js"></script>
<link href="index_files/highlightjs-9.12.0/default.css" rel="stylesheet" />
<script src="index_files/highlightjs-9.12.0/highlight.js"></script>
<script src="index_files/htmlwidgets-1.3/htmlwidgets.js"></script>
<link href="index_files/datatables-css-0.0.0/datatables-crosstalk.css" rel="stylesheet" />
<script src="index_files/datatables-binding-0.5/datatables.js"></script>
<link href="index_files/dt-core-1.10.16/css/jquery.dataTables.min.css" rel="stylesheet" />
<link href="index_files/dt-core-1.10.16/css/jquery.dataTables.extra.css" rel="stylesheet" />
<script src="index_files/dt-core-1.10.16/js/jquery.dataTables.min.js"></script>
<link href="index_files/crosstalk-1.0.0/css/crosstalk.css" rel="stylesheet" />
<script src="index_files/crosstalk-1.0.0/js/crosstalk.min.js"></script>
<meta name="twitter:card" content="summary_large_image">
<meta name="twitter:site" content="@jonaslindeloev">
<meta name="twitter:title" content="Common statistical tests are linear models (or: how to teach stats)">
<meta name="twitter:image" content="https://lindeloev.github.io/tests-as-linear/linear_tests_cheat_sheet.png">
<!-- Global site tag (gtag.js) - Google Analytics -->
<script async src="https://www.googletagmanager.com/gtag/js?id=UA-1026978-2"></script>
<script>
window.dataLayer = window.dataLayer || [];
function gtag(){dataLayer.push(arguments);}
gtag('js', new Date());
gtag('config', 'UA-1026978-2');
</script>
<style type="text/css">code{white-space: pre;}</style>
<style type="text/css">
pre:not([class]) {
background-color: white;
}
</style>
<script type="text/javascript">
if (window.hljs) {
hljs.configure({languages: []});
hljs.initHighlightingOnLoad();
if (document.readyState && document.readyState === "complete") {
window.setTimeout(function() { hljs.initHighlighting(); }, 0);
}
}
</script>
<style type="text/css">
h1 {
font-size: 34px;
}
h1.title {
font-size: 38px;
}
h2 {
font-size: 30px;
}
h3 {
font-size: 24px;
}
h4 {
font-size: 18px;
}
h5 {
font-size: 16px;
}
h6 {
font-size: 12px;
}
.table th:not([align]) {
text-align: left;
}
</style>
<style type = "text/css">
.main-container {
max-width: 940px;
margin-left: auto;
margin-right: auto;
}
code {
color: inherit;
background-color: rgba(0, 0, 0, 0.04);
}
img {
max-width:100%;
height: auto;
}
.tabbed-pane {
padding-top: 12px;
}
.html-widget {
margin-bottom: 20px;
}
button.code-folding-btn:focus {
outline: none;
}
summary {
display: list-item;
}
</style>
<!-- tabsets -->
<style type="text/css">
.tabset-dropdown > .nav-tabs {
display: inline-table;
max-height: 500px;
min-height: 44px;
overflow-y: auto;
background: white;
border: 1px solid #ddd;
border-radius: 4px;
}
.tabset-dropdown > .nav-tabs > li.active:before {
content: "";
font-family: 'Glyphicons Halflings';
display: inline-block;
padding: 10px;
border-right: 1px solid #ddd;
}
.tabset-dropdown > .nav-tabs.nav-tabs-open > li.active:before {
content: "";
border: none;
}
.tabset-dropdown > .nav-tabs.nav-tabs-open:before {
content: "";
font-family: 'Glyphicons Halflings';
display: inline-block;
padding: 10px;
border-right: 1px solid #ddd;
}
.tabset-dropdown > .nav-tabs > li.active {
display: block;
}
.tabset-dropdown > .nav-tabs > li > a,
.tabset-dropdown > .nav-tabs > li > a:focus,
.tabset-dropdown > .nav-tabs > li > a:hover {
border: none;
display: inline-block;
border-radius: 4px;
}
.tabset-dropdown > .nav-tabs.nav-tabs-open > li {
display: block;
float: none;
}
.tabset-dropdown > .nav-tabs > li {
display: none;
}
</style>
<script>
$(document).ready(function () {
window.buildTabsets("TOC");
});
$(document).ready(function () {
$('.tabset-dropdown > .nav-tabs > li').click(function () {
$(this).parent().toggleClass('nav-tabs-open')
});
});
</script>
<!-- code folding -->
<script>
$(document).ready(function () {
// move toc-ignore selectors from section div to header
$('div.section.toc-ignore')
.removeClass('toc-ignore')
.children('h1,h2,h3,h4,h5').addClass('toc-ignore');
// establish options
var options = {
selectors: "h1,h2",
theme: "bootstrap3",
context: '.toc-content',
hashGenerator: function (text) {
return text.replace(/[.\\/?&!#<>]/g, '').replace(/\s/g, '_').toLowerCase();
},
ignoreSelector: ".toc-ignore",
scrollTo: 0
};
options.showAndHide = false;
options.smoothScroll = true;
// tocify
var toc = $("#TOC").tocify(options).data("toc-tocify");
});
</script>
<style type="text/css">
#TOC {
margin: 25px 0px 20px 0px;
}
@media (max-width: 768px) {
#TOC {
position: relative;
width: 100%;
}
}
.toc-content {
padding-left: 30px;
padding-right: 40px;
}
div.main-container {
max-width: 1200px;
}
div.tocify {
width: 20%;
max-width: 260px;
max-height: 85%;
}
@media (min-width: 768px) and (max-width: 991px) {
div.tocify {
width: 25%;
}
}
@media (max-width: 767px) {
div.tocify {
width: 100%;
max-width: none;
}
}
.tocify ul, .tocify li {
line-height: 20px;
}
.tocify-subheader .tocify-item {
font-size: 0.90em;
padding-left: 25px;
text-indent: 0;
}
.tocify .list-group-item {
border-radius: 0px;
}
.tocify-subheader {
display: inline;
}
.tocify-subheader .tocify-item {
font-size: 0.95em;
}
</style>
</head>
<body>
<div class="container-fluid main-container">
<!-- setup 3col/9col grid for toc_float and main content -->
<div class="row-fluid">
<div class="col-xs-12 col-sm-4 col-md-3">
<div id="TOC" class="tocify">
</div>
</div>
<div class="toc-content col-xs-12 col-sm-8 col-md-9">
<div class="fluid-row" id="header">
<h1 class="title toc-ignore">Common statistical tests are linear models (or: how to teach stats)</h1>
</div>
<p><link rel="stylesheet" type="text/css" href="include/style.css"></p>
<!-- From https://stackoverflow.com/a/37839683/1297830 -->
<link rel="stylesheet" type="text/css" href="include/hideOutput.css">
<script src="include/hideOutput.js"></script>
<p>By Jonas Kristoffer Lindeløv (<a href="https://lindeloev.net">blog</a>, <a href="http://personprofil.aau.dk/117060">profile</a>).<br />Last updated: 02 April, 2019 (See <a href="https://github.com/lindeloev/tests-as-linear/commits/master">changelog</a>).</p>
<!-- Social sharing. From simplesharebuttons.com -->
<style type="text/css">
#share-buttons img {
width: 40px;
padding-right: 15px;
border: 0;
box-shadow: 0;
display: inline;
vertical-align: top;
}
</style>
<div id="share-buttons">
<!-- Twitter -->
<script async src="https://platform.twitter.com/widgets.js" charset="utf-8"></script>
<p><a href="https://twitter.com/intent/tweet?text=Common%20statistical%20tests%20are%20linear%20models%20(or:%20how%20to%20teach%20stats)%20https%3A%2F%2Flindeloev.github.io%2Ftests-as-linear%20via%20%40jonaslindeloev" class="twitter-hashtag-button" data-size="large" data-related="jonaslindeloev" data-show-count="false">Share on Twitter</a> <!-- Facebook --><a href="http://www.facebook.com/sharer.php?u=https://lindeloev.github.io/tests-as-linear/" target="_blank"><img src="https://simplesharebuttons.com/images/somacro/facebook.png" alt="Facebook" /></a><!-- LinkedIn --><a href="http://www.linkedin.com/shareArticle?mini=true&url=https://lindeloev.github.io/tests-as-linear/" target="_blank"><img src="https://simplesharebuttons.com/images/somacro/linkedin.png" alt="LinkedIn" /></a><!-- Digg --><a href="http://www.digg.com/submit?url=https://lindeloev.github.io/tests-as-linear/" target="_blank"><img src="https://simplesharebuttons.com/images/somacro/diggit.png" alt="Digg" /></a><!-- Reddit --><a href="http://reddit.com/submit?url=https://lindeloev.github.io/tests-as-linear/&title=Common statistical tests are linear models (or: how to teach stats)" target="_blank"><img src="https://simplesharebuttons.com/images/somacro/reddit.png" alt="Reddit" /></a><!-- Email --><a href="mailto:?Subject=Common statistical tests are linear models (or: how to teach stats)&Body=https://lindeloev.github.io/tests-as-linear/"><img src="https://simplesharebuttons.com/images/somacro/email.png" alt="Email" /></a></p>
</div>
<p><br /></p>
<p>This document is summarised in the table below. It shows the linear models underlying common parametric and “non-parametric”" tests. Formulating all the tests in the same language highlights the many similarities between them. Get it <a href="linear_tests_cheat_sheet.png">as an image</a> or <a href="linear_tests_cheat_sheet.pdf">as a PDF</a>.</p>
<hr />
<p><a href="linear_tests_cheat_sheet.pdf"><img src="linear_tests_cheat_sheet.png" /></a></p>
<hr />
<div id="the-simplicity-underlying-common-tests" class="section level1">
<h1><span class="header-section-number">1</span> The simplicity underlying common tests</h1>
<p>Most of the common statistical models (t-test, correlation, ANOVA; chi-square, etc.) are special cases of linear models or a very close approximation. This beautiful simplicity means that there is less to learn. In particular, it all comes down to <span class="math inline">\(y = a \cdot x + b\)</span> which most students know from highschool. Unfortunately, stats intro courses are usually taught as if each test is an independent tool, needlessly making life more complicated for students and teachers alike.</p>
<p>This needless complexity multiplies when students try to rote learn the parametric assumptions underlying each test separately rather than deducing them from the linear model.</p>
<p>For this reason, I think that teaching linear models first and foremost and <em>then</em> name-dropping the special cases along the way makes for an excellent teaching strategy, emphasizing <em>understanding</em> over rote learning. Since linear models are the same across frequentist, Bayesian, and permutation-based inferences, I’d argue that it’s better to start with modeling than p-values, type-1 errors, Bayes factors, or other inferences.</p>
<p>Concerning the teaching of <em>“non-parametric”</em> tests in intro-courses, I think that we can justify <a href="https://en.wikipedia.org/wiki/Lie-to-children">lying-to-children</a> and teach “non-parametric”" tests as if they are merely ranked versions of the corresponding parametric tests. It is much better for students to think “ranks!” than to believe that you can magically throw away assumptions. Indeed, the Bayesian equivalents of “non-parametric”" tests implemented in <a href="https://jasp-stats.org">JASP</a> <a href="https://arxiv.org/abs/1712.06941">literally just do (latent) ranking</a> and that’s it. For the frequentist “non-parametric”" tests considered here, this approach is highly accurate for N > 15.</p>
<p>Use the menu to jump to your favourite section. There are links to lots of similar (though more scattered) stuff under <a href="#links">sources</a> and <a href="#course">teaching materials</a>. I hope that you will join in suggesting improvements or submitting improvements yourself in <a href="https://github.com/lindeloev/tests-as-linear">the Github repo to this page</a>. Let’s make it awesome!</p>
</div>
<div id="settings-and-toy-data" class="section level1">
<h1><span class="header-section-number">2</span> Settings and toy data</h1>
Unfold this if you want to see functions and other settings for this notebook:
<div class="fold s">
<pre class="r"><code># Load packages for data handling and plotting
library(tidyverse)
library(patchwork)
library(broom)
# Reproducible "random" results
set.seed(40)
# To show tables. Rounds
print_df = function(D, decimals=4, navigate=FALSE) {
DT::datatable(mutate_if(D, is.numeric, round, decimals),
rownames = FALSE,
options = list(
searching=FALSE,
lengthChange=FALSE,
ordering=FALSE,
autoWidth=TRUE,
bPaginate=navigate,
bInfo=navigate,
paging=navigate
)
)
}
# Generate normal data with known parameters
rnorm_fixed = function(N, mu=0, sd=1) scale(rnorm(N))*sd + mu
# Plot style.
theme_axis = function(P, jitter=FALSE, xlim=c(-0.5, 2), ylim=c(-0.5, 2), legend.position=NULL) {
P = P + theme_bw(15) +
geom_segment(x=-1000, xend=1000, y=0, yend=0, lty=2, color='dark gray', lwd=0.5) +
geom_segment(x=0, xend=0, y=-1000, yend=1000, lty=2, color='dark gray', lwd=0.5) +
coord_cartesian(xlim=xlim, ylim=ylim) +
theme(axis.title = element_blank(),
axis.text = element_blank(),
axis.ticks = element_blank(),
panel.border = element_blank(),
panel.grid = element_blank(),
legend.position = legend.position)
# Return jittered or non-jittered plot?
if(jitter) {
P + geom_jitter(width=0.1, size=2)
}
else {
P + geom_point(size=2)
}
}</code></pre>
</div>
<p>For a start, we’ll keep it simple and play with three standard normals in wide (<code>a</code>, <code>b</code>, <code>c</code>) and long format (<code>value</code>, <code>group</code>):</p>
<pre class="r"><code># Wide format (sort of)
y = rnorm_fixed(50, mu=0.3, sd=2) # Almost zero mean
x = rnorm_fixed(50, mu=0, sd=1) # Used in correlation where this is on x-axis
y2 = rnorm_fixed(50, mu=0.5, sd=1.5) # Used in two means
# Long format data with indicator
value = c(y, y2)
group = rep(c('y1', 'y2'), each = 50)</code></pre>
</div>
<div id="correlation" class="section level1">
<h1><span class="header-section-number">3</span> Pearson and Spearman correlation</h1>
<div id="theory-as-linear-models" class="section level3">
<h3><span class="header-section-number">3.0.1</span> Theory: As linear models</h3>
<p>Model: the recipe for <span class="math inline">\(y\)</span> is a slope (<span class="math inline">\(\beta_1\)</span>) times <span class="math inline">\(x\)</span> plus an intercept (<span class="math inline">\(\beta_0\)</span>, aka a straight line).</p>
<p><span class="math inline">\(y = \beta_0 + \beta_1 x \qquad \mathcal{H}_0: \beta_1 = 0\)</span></p>
<p>… which is a math-y way of writing the good old <span class="math inline">\(y = ax + b\)</span> (here ordered as <span class="math inline">\(y = b + ax\)</span>). In R we are lazy and write <code>y ~ 1 + x</code> which R reads like <code>y = 1*number + x*othernumber</code> and the task of t-tests, lm, etc., is simply to find the numbers that best predict <span class="math inline">\(y\)</span>.</p>
<p>Either way you write it, it’s an intercept (<span class="math inline">\(\beta_0\)</span>) and a slope (<span class="math inline">\(\beta_1\)</span>) yielding a straight line:</p>
<div class="fold s">
<pre class="r"><code># Fixed correlation
D_correlation = data.frame(MASS::mvrnorm(30, mu=c(0.9, 0.9), Sigma=matrix(c(1, 0.8, 1, 0.8), ncol=2), empirical=TRUE)) # Correlated data
# Add labels (for next plot)
D_correlation$label_num = sprintf('(%.1f,%.1f)', D_correlation$X1, D_correlation$X2)
D_correlation$label_rank = sprintf('(%i,%i)', rank(D_correlation$X1), rank(D_correlation$X2))
# Plot it
fit = lm(I(X2*0.5+0.4) ~ I(X1*0.5+0.2), D_correlation)
intercept_pearson = coefficients(fit)[1]
P_pearson = ggplot(D_correlation, aes(x=X1*0.5+0.2, y=X2*0.5+0.4)) +
geom_smooth(method=lm, se=FALSE, lwd=2, aes(colour='beta_1')) +
geom_segment(x=-100, xend=100,
y=intercept_pearson, yend=intercept_pearson,
lwd=2, aes(color="beta_0")) +
scale_color_manual(name=NULL, values=c("blue", "red"), labels=c(bquote(beta[0]*" (intercept)"), bquote(beta[1]*" (slope)")))
theme_axis(P_pearson, legend.position=c(0.4, 0.9))</code></pre>
<p><img src="index_files/figure-html/unnamed-chunk-4-1.png" width="576" style="display: block; margin: auto;" /></p>
</div>
<p>This is often simply called a <strong>regression</strong> model which can be extended to <strong>multiple regression</strong> where there are several <span class="math inline">\(\beta\)</span>s and on the right-hand side multiplied with the predictors. Everything below, from <a href="#t1">one-sample t-test</a> to <a href="#anova2">two-way ANOVA</a> are just special cases of this system. Nothing more, nothing less.</p>
<p>As the name implies, the <strong>Spearman rank correlation</strong> is a <strong>Pearson correlation</strong> on rank-transformed <span class="math inline">\(x\)</span> and <span class="math inline">\(y\)</span>:</p>
<p><span class="math inline">\(rank(y) = \beta_0 + \beta_1 \cdot rank(x) \qquad \mathcal{H}_0: \beta_1 = 0\)</span></p>
<p>I’ll introduce <a href="#rank">ranks</a> in a minute. For now, notice that the correlation coefficient of the linear model is identical to a “real” Pearson correlation, but p-values are an approximation which is is <a href="simulations/simulate_spearman.html">appropriate for samples greater than N=10 and almost perfect when N > 20</a>.</p>
<p>Such a nice and non-mysterious equivalence that many students are left unaware of! Visualizing them side by side including data labels, we see this rank-transformation in action:</p>
<div class="fold s">
<pre class="r"><code># Spearman intercept
intercept_spearman = coefficients(lm(rank(X2) ~ rank(X1), D_correlation))[1]
# Spearman plot
P_spearman = ggplot(D_correlation, aes(x=rank(X1), y=rank(X2))) +
geom_smooth(method=lm, se=FALSE, lwd=2, aes(color='beta_1')) +
geom_text(aes(label=label_rank), nudge_y=1, size=3, color='dark gray') +
geom_segment(x=-100, xend=100,
y=intercept_spearman, yend=intercept_spearman,
lwd=2, aes(color='beta_0')) +
scale_color_manual(name=NULL, values=c("blue", "red"), labels=c(bquote(beta[0]*" (intercept)"), bquote(beta[1]*" (slope)")))
# Stich together using patchwork
(theme_axis(P_pearson, legend.position=c(0.5, 0.1)) + geom_text(aes(label=label_num), nudge_y=0.1, size=3, color='dark gray') + labs(title=' Pearson')) + (theme_axis(P_spearman, xlim=c(-7.5, 30), ylim=c(-7.5, 30), legend.position=c(0.5, 0.1)) + labs(title=' Spearman'))</code></pre>
<p><img src="index_files/figure-html/unnamed-chunk-5-1.png" width="768" style="display: block; margin: auto;" /></p>
</div>
</div>
<div id="rank" class="section level3">
<h3><span class="header-section-number">3.0.2</span> Theory: rank-transformation</h3>
<p><code>rank</code> simply takes a list of numbers and “replace” them with the integers of their rank (1st smallest, 2nd smallest, 3rd smallest, etc.). So the result of the rank-transformation <code>rank(c(3.6, 3.4, -5.0, 8.2))</code> is <code>3, 2, 1, 4</code>. See that in the figure above?</p>
<p>A <em>signed</em> rank is the same, just where we rank according to absolute size first and then add in the sign second. So the signed rank here would be <code>2, 1, -3, 4</code>. Or in code:</p>
<pre class="r"><code>signed_rank = function(x) sign(x) * rank(abs(x))</code></pre>
<p>I hope I don’t offend anyone when I say that ranks are easy; yet it’s all you need to do to convert most parametric tests into their “non-parametric” counterparts! One interesting implication is that <em>many “non-parametric tests” are about as parametric as their parametric counterparts with means, standard deviations, homogeneity of variance, etc. - just on rank-transformed data</em>. That’s why I put “non-parametric” in quotation marks.</p>
</div>
<div id="r-code-pearson-correlation" class="section level3">
<h3><span class="header-section-number">3.0.3</span> R code: Pearson correlation</h3>
<p>It couldn’t be much simpler to run these models in R. They yield identical <code>p</code> and <code>t</code>, but there’s a catch: <code>lm</code> gives you the <em>slope</em> and even though that is usually much more interpretable and informative than the <em>correlation coefficient</em> <em>r</em>, you may still want <em>r</em>. Luckily, the slope becomes <code>r</code> if <code>x</code> and <code>y</code> have a standard deviation of exactly 1. You can do this using <code>scale(x)</code> or <code>I(x/sd(x))</code>:</p>
<pre class="r"><code>a = cor.test(y, x, method = "pearson") # Built-in
b = lm(y ~ 1 + x) # Equivalent linear model: y = Beta0*1 + Beta1*x
c = lm(scale(y) ~ 1 + scale(x)) # On scaled vars to recover r</code></pre>
Results: <div id="htmlwidget-58f61054e232de0dbaf9" style="width:100%;height:auto;" class="datatables html-widget"></div>
<script type="application/json" data-for="htmlwidget-58f61054e232de0dbaf9">{"x":{"filter":"none","data":[["cor.test","lm","lm scaled"],[-0.2318,-0.4636,-0.2318],[0.1053,0.1053,0.1053],[-1.6507,-1.6507,-1.6507],[-0.4792,-1.0282,-0.5141],[0.0498,0.1011,0.0505]],"container":"<table class=\"display\">\n <thead>\n <tr>\n <th>model<\/th>\n <th>r<\/th>\n <th>p.value<\/th>\n <th>t<\/th>\n <th>conf.low<\/th>\n <th>conf.high<\/th>\n <\/tr>\n <\/thead>\n<\/table>","options":{"searching":false,"lengthChange":false,"ordering":false,"autoWidth":true,"bPaginate":false,"bInfo":false,"paging":false,"columnDefs":[{"className":"dt-right","targets":[1,2,3,4,5]}],"order":[],"orderClasses":false}},"evals":[],"jsHooks":[]}</script>
<div class="fold o">
<pre><code>##
## Pearson's product-moment correlation
##
## data: y and x
## t = -1.6507, df = 48, p-value = 0.1053
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.47920849 0.04978276
## sample estimates:
## cor
## -0.2317767
##
##
## Call:
## lm(formula = y ~ 1 + x)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.3393 -1.6593 0.3349 1.3629 3.5214
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.3000 0.2780 1.079 0.286
## x -0.4636 0.2808 -1.651 0.105
##
## Residual standard error: 1.966 on 48 degrees of freedom
## Multiple R-squared: 0.05372, Adjusted R-squared: 0.03401
## F-statistic: 2.725 on 1 and 48 DF, p-value: 0.1053
##
##
## Call:
## lm(formula = scale(y) ~ 1 + scale(x))
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.6697 -0.8297 0.1675 0.6815 1.7607
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.341e-17 1.390e-01 0.000 1.000
## scale(x) -2.318e-01 1.404e-01 -1.651 0.105
##
## Residual standard error: 0.9828 on 48 degrees of freedom
## Multiple R-squared: 0.05372, Adjusted R-squared: 0.03401
## F-statistic: 2.725 on 1 and 48 DF, p-value: 0.1053</code></pre>
</div>
<p>The CIs are not exactly identical, but very close.</p>
</div>
<div id="r-code-spearman-correlation" class="section level3">
<h3><span class="header-section-number">3.0.4</span> R code: Spearman correlation</h3>
<p>Note that we can interpret the slope which is the number of ranks <span class="math inline">\(y\)</span> change for each rank on <span class="math inline">\(x\)</span>. I think that this is a pretty interesting number. However, the intercept is less interpretable since it lies at <span class="math inline">\(rank(x) = 0\)</span> which is impossible since x starts at 1.</p>
<p>See the identical <code>r</code> (now “rho”) and <code>p</code>:</p>
<pre class="r"><code># Spearman correlation
a = cor.test(y, x, method = "spearman") # Built-in
b = lm(rank(y) ~ 1 + rank(x)) # Equivalent linear model</code></pre>
Let’s look at the results: <div id="htmlwidget-18ac5bd784205e06f9a7" style="width:100%;height:auto;" class="datatables html-widget"></div>
<script type="application/json" data-for="htmlwidget-18ac5bd784205e06f9a7">{"x":{"filter":"none","data":[["cor.test","lm"],[0.1135,0.1135],[-0.2266,-0.2266]],"container":"<table class=\"display\">\n <thead>\n <tr>\n <th>model<\/th>\n <th>p.value<\/th>\n <th>rho<\/th>\n <\/tr>\n <\/thead>\n<\/table>","options":{"searching":false,"lengthChange":false,"ordering":false,"autoWidth":true,"bPaginate":false,"bInfo":false,"paging":false,"columnDefs":[{"className":"dt-right","targets":[1,2]}],"order":[],"orderClasses":false}},"evals":[],"jsHooks":[]}</script>
<div class="fold o">
<pre><code>##
## Spearman's rank correlation rho
##
## data: y and x
## S = 25544, p-value = 0.1135
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.2266026
##
##
## Call:
## lm(formula = rank(y) ~ 1 + rank(x))
##
## Residuals:
## Min 1Q Median 3Q Max
## -26.4655 -11.5603 0.4458 11.5628 25.6921
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 31.2784 4.1191 7.593 9.11e-10 ***
## rank(x) -0.2266 0.1406 -1.612 0.114
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 14.35 on 48 degrees of freedom
## Multiple R-squared: 0.05135, Adjusted R-squared: 0.03159
## F-statistic: 2.598 on 1 and 48 DF, p-value: 0.1135</code></pre>
</div>
</div>
</div>
<div id="one-mean" class="section level1">
<h1><span class="header-section-number">4</span> One mean</h1>
<div id="t1" class="section level2">
<h2><span class="header-section-number">4.1</span> One sample t-test and Wilcoxon signed-rank</h2>
<div id="theory-as-linear-models-1" class="section level3">
<h3><span class="header-section-number">4.1.1</span> Theory: As linear models</h3>
<p><strong>t-test</strong> model: A single number predicts <span class="math inline">\(y\)</span>.</p>
<p><span class="math inline">\(y = \beta_0 \qquad \mathcal{H}_0: \beta_0 = 0\)</span></p>
<p>In other words, it’s our good old <span class="math inline">\(y = \beta_0 + \beta_1*x\)</span> where the last term is gone since there is no <span class="math inline">\(x\)</span> (essentially <span class="math inline">\(x=0\)</span>, see left figure below).</p>
<p>The same is to a very close approximately true for <strong>Wilcoxon signed-rank test</strong>, just with the <a href="#rank">signed ranks</a> of <span class="math inline">\(y\)</span> instead of <span class="math inline">\(y\)</span> itself (see right panel below).</p>
<p><span class="math inline">\(signed\_rank(y) = \beta_0\)</span></p>
<p><a href="simulations/simulate_wilcoxon.html">This approximation is good enough when the sample size is larger than 14 and almost perfect if the sample size is larger than 50</a>.</p>
<div class="fold s">
<pre class="r"><code># T-test
D_t1 = data.frame(y=rnorm_fixed(20, 0.5, 0.6),
x=runif(20, 0.93, 1.07)) # Fix mean and SD
P_t1 = ggplot(D_t1, aes(y=y, x=0)) +
stat_summary(fun.y=mean, geom = "errorbar", aes(ymax = ..y.., ymin = ..y.., color='beta_0'), lwd=2) +
scale_color_manual(name=NULL, values=c("blue"), labels=c(bquote(beta[0]*" (intercept)"))) +
geom_text(aes(label=round(y, 1)), nudge_x = 0.2, size=3, color='dark gray') +
labs(title=' T-test')
# Wilcoxon
D_t1_rank = data.frame(y = signed_rank(D_t1$y))
P_t1_rank = ggplot(D_t1_rank, aes(y=y, x=0)) +
stat_summary(fun.y=mean, geom = "errorbar", aes(ymax = ..y.., ymin = ..y.., color='beta_0'), lwd=2) +
scale_color_manual(name=NULL, values=c("blue"), labels=c(bquote(beta[0]*" (intercept)"))) +
geom_text(aes(label=y), nudge_x=0.2, size=3, color='dark gray') +
labs(title=' Wilcoxon')
# Stich together using patchwork
theme_axis(P_t1, ylim=c(-1, 2), legend.position=c(0.6, 0.1)) +
theme_axis(P_t1_rank, ylim=NULL, legend.position=c(0.6, 0.1))</code></pre>
<p><img src="index_files/figure-html/unnamed-chunk-13-1.png" width="672" style="display: block; margin: auto;" /></p>
</div>
</div>
<div id="r-code-one-sample-t-test" class="section level3">
<h3><span class="header-section-number">4.1.2</span> R code: One-sample t-test</h3>
<p>Try running the R code below and see that the linear model (<code>lm</code>) produces the same <span class="math inline">\(t\)</span>, <span class="math inline">\(p\)</span>, and <span class="math inline">\(r\)</span> as the built-in <code>t.test</code>. The confidence interval is not presented in the output of <code>lm</code> but is also identical if you use <code>confint(lm(...))</code>:</p>
<pre class="r"><code># Built-in t-test
a = t.test(y)
# Equivalent linear model: intercept-only
b = lm(y ~ 1)</code></pre>
<p>Results:</p>
<div id="htmlwidget-635d1460ab26b7a06864" style="width:100%;height:auto;" class="datatables html-widget"></div>
<script type="application/json" data-for="htmlwidget-635d1460ab26b7a06864">{"x":{"filter":"none","data":[["t.test","lm"],[0.3,0.3],[0.294,0.294],[1.0607,1.0607],[49,49],[-0.2684,-0.2684],[0.8684,0.8684]],"container":"<table class=\"display\">\n <thead>\n <tr>\n <th>model<\/th>\n <th>mean<\/th>\n <th>p.value<\/th>\n <th>t<\/th>\n <th>df<\/th>\n <th>conf.low<\/th>\n <th>conf.high<\/th>\n <\/tr>\n <\/thead>\n<\/table>","options":{"searching":false,"lengthChange":false,"ordering":false,"autoWidth":true,"bPaginate":false,"bInfo":false,"paging":false,"columnDefs":[{"className":"dt-right","targets":[1,2,3,4,5,6]}],"order":[],"orderClasses":false}},"evals":[],"jsHooks":[]}</script>
<div class="fold o">
<pre><code>##
## One Sample t-test
##
## data: y
## t = 1.0607, df = 49, p-value = 0.294
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -0.2683937 0.8683937
## sample estimates:
## mean of x
## 0.3
##
##
## Call:
## lm(formula = y ~ 1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.521 -1.673 0.481 1.427 3.795
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.3000 0.2828 1.061 0.294
##
## Residual standard error: 2 on 49 degrees of freedom</code></pre>
</div>
</div>
<div id="r-code-wilcoxon-signed-rank-test" class="section level3">
<h3><span class="header-section-number">4.1.3</span> R code: Wilcoxon signed-rank test</h3>
<p>In addition to matching <code>p</code>-values, <code>lm</code> also gives us the mean signed rank, which I find to be an informative number.</p>
<pre class="r"><code># Built-in
a = wilcox.test(y)
# Equivalent linear model
b = lm(signed_rank(y) ~ 1) # See? Same as above, just on signed ranks
# Bonus: of course also works for one-sample t-test
c = t.test(signed_rank(y))</code></pre>
<p>Results:</p>
<div id="htmlwidget-c43cb8fa604a5da8e2b7" style="width:100%;height:auto;" class="datatables html-widget"></div>
<script type="application/json" data-for="htmlwidget-c43cb8fa604a5da8e2b7">{"x":{"filter":"none","data":[["wilcox.test","lm","t.test"],[0.3693,0.3721,0.3721],[null,3.74,3.74]],"container":"<table class=\"display\">\n <thead>\n <tr>\n <th>model<\/th>\n <th>p.value<\/th>\n <th>mean_rank<\/th>\n <\/tr>\n <\/thead>\n<\/table>","options":{"searching":false,"lengthChange":false,"ordering":false,"autoWidth":true,"bPaginate":false,"bInfo":false,"paging":false,"columnDefs":[{"className":"dt-right","targets":[1,2]}],"order":[],"orderClasses":false}},"evals":[],"jsHooks":[]}</script>
<div class="fold o">
<pre><code>##
## Wilcoxon signed rank test with continuity correction
##
## data: y
## V = 731, p-value = 0.3693
## alternative hypothesis: true location is not equal to 0
##
##
## Call:
## lm(formula = signed_rank(y) ~ 1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -49.74 -25.49 4.76 22.76 46.26
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.740 4.151 0.901 0.372
##
## Residual standard error: 29.36 on 49 degrees of freedom
##
##
## One Sample t-test
##
## data: signed_rank(y)
## t = 0.90088, df = 49, p-value = 0.3721
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -4.60275 12.08275
## sample estimates:
## mean of x
## 3.74</code></pre>
</div>
</div>
</div>
<div id="tpair" class="section level2">
<h2><span class="header-section-number">4.2</span> Paired samples t-test and Wilcoxon matched pairs</h2>
<div id="theory-as-linear-models-2" class="section level3">
<h3><span class="header-section-number">4.2.1</span> Theory: As linear models</h3>
<p><strong>t-test</strong> model: a single number (intercept) predicts the pairwise differences.</p>
<p><span class="math inline">\(y_2-y_1 = \beta_0 \qquad \mathcal{H}_0: \beta_0 = 0\)</span></p>
<p>This means that there is just one <span class="math inline">\(y = y_2 - y_1\)</span> to predict and it becomes a <a href="#t1">one-sample t-test</a> on the pairwise differences. The visualization is therefore also the same as for the one-sample t-test. At the risk of overcomplicating a simple substraction, you can think of these pairwise differences as slopes (see left panel of the figure), which we can represent as y-offsets (see right panel of the figure):</p>
<div class="fold s">
<pre class="r"><code># Data for plot
N = nrow(D_t1)
start = rnorm_fixed(N, 0.2, 0.3)
D_tpaired = data.frame(x = rep(c(0, 1), each=N), y = c(start, start + D_t1$y), id=1:N)
# Plot
P_tpaired = ggplot(D_tpaired, aes(x=x, y=y)) +
geom_line(aes(group=id)) +
labs(title=' Pairs')
# Use patchwork to put them side-by-side
theme_axis(P_tpaired) + theme_axis(P_t1, legend.position=c(0.6, 0.1))</code></pre>
<p><img src="index_files/figure-html/unnamed-chunk-20-1.png" width="672" style="display: block; margin: auto;" /></p>
</div>
<p>Similarly, the <strong>Wilcoxon matched pairs</strong> only differ from <strong>Wilcoxon signed-rank</strong> in that it’s testing the signed ranks of the pairwise <span class="math inline">\(y-x\)</span> differences.</p>
<p><span class="math inline">\(signed\_rank(y_2-y_1) = \beta_0 \qquad \mathcal{H}_0: \beta_0 = 0\)</span></p>
</div>
<div id="r-code-paired-sample-t-test" class="section level3">
<h3><span class="header-section-number">4.2.2</span> R code: Paired sample t-test</h3>
<pre class="r"><code>a = t.test(y, y2, paired = TRUE) # Built-in paired t-test
b = lm(y - y2 ~ 1) # Equivalent linear model</code></pre>
<p>Results:</p>
<div id="htmlwidget-77cd519bd4f19792ec4e" style="width:100%;height:auto;" class="datatables html-widget"></div>
<script type="application/json" data-for="htmlwidget-77cd519bd4f19792ec4e">{"x":{"filter":"none","data":[["t.test","lm"],[-0.2,-0.2],[0.601,0.601],[49,49],[-0.5264,-0.5264],[-0.9635,-0.9635],[0.5635,0.5635]],"container":"<table class=\"display\">\n <thead>\n <tr>\n <th>model<\/th>\n <th>mean<\/th>\n <th>p.value<\/th>\n <th>df<\/th>\n <th>t<\/th>\n <th>conf.low<\/th>\n <th>conf.high<\/th>\n <\/tr>\n <\/thead>\n<\/table>","options":{"searching":false,"lengthChange":false,"ordering":false,"autoWidth":true,"bPaginate":false,"bInfo":false,"paging":false,"columnDefs":[{"className":"dt-right","targets":[1,2,3,4,5,6]}],"order":[],"orderClasses":false}},"evals":[],"jsHooks":[]}</script>
<div class="fold o">
<pre><code>##
## Paired t-test
##
## data: y and y2
## t = -0.52642, df = 49, p-value = 0.601
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.9634894 0.5634894
## sample estimates:
## mean of the differences
## -0.2
##
##
## Call:
## lm(formula = y - y2 ~ 1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.5253 -1.5642 -0.0844 1.9715 5.0361
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.2000 0.3799 -0.526 0.601
##
## Residual standard error: 2.686 on 49 degrees of freedom</code></pre>
</div>
</div>
<div id="r-code-wilcoxon-matched-pairs" class="section level3">
<h3><span class="header-section-number">4.2.3</span> R code: Wilcoxon matched pairs</h3>
<p>Again, we do the signed-ranks trick. This is still an approximation, but a close one:</p>
<pre class="r"><code># Built-in Wilcoxon matched pairs
a = wilcox.test(y, y2, paired = TRUE)
# Equivalent linear model:
b = lm(signed_rank(y - y2) ~ 1)
# Bonus: identical to one-sample t-test ong signed ranks
c = t.test(signed_rank(y - y2))</code></pre>
<p>Results:</p>
<div id="htmlwidget-73383d6f3c53896f5ef8" style="width:100%;height:auto;" class="datatables html-widget"></div>
<script type="application/json" data-for="htmlwidget-73383d6f3c53896f5ef8">{"x":{"filter":"none","data":[["wilcox.test","lm","t.test"],[0.8243,0.8232,0.8232],[null,-0.94,-0.94]],"container":"<table class=\"display\">\n <thead>\n <tr>\n <th>model<\/th>\n <th>p.value<\/th>\n <th>mean_rank_diff<\/th>\n <\/tr>\n <\/thead>\n<\/table>","options":{"searching":false,"lengthChange":false,"ordering":false,"autoWidth":true,"bPaginate":false,"bInfo":false,"paging":false,"columnDefs":[{"className":"dt-right","targets":[1,2]}],"order":[],"orderClasses":false}},"evals":[],"jsHooks":[]}</script>
<div class="fold o">
<pre><code>##
## Wilcoxon signed rank test with continuity correction
##
## data: y and y2
## V = 614, p-value = 0.8243
## alternative hypothesis: true location shift is not equal to 0
##
##
## Call:
## lm(formula = signed_rank(y - y2) ~ 1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -49.06 -23.81 -2.56 26.19 46.94
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.940 4.184 -0.225 0.823
##
## Residual standard error: 29.58 on 49 degrees of freedom
##
##
## One Sample t-test
##
## data: signed_rank(y - y2)
## t = -0.22469, df = 49, p-value = 0.8232
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -9.347227 7.467227
## sample estimates:
## mean of x
## -0.94</code></pre>
</div>
<p>For large sample sizes (N >> 100), this approaches the <strong>sign test</strong> to a reasonable degree, but this approximation is too inaccurate to flesh out here.</p>
</div>
</div>
</div>
<div id="two-means" class="section level1">
<h1><span class="header-section-number">5</span> Two means</h1>
<div id="t2" class="section level2">
<h2><span class="header-section-number">5.1</span> Independent t-test and Mann-Whitney U</h2>
<div id="theory-as-linear-models-3" class="section level3">
<h3><span class="header-section-number">5.1.1</span> Theory: As linear models</h3>
<p>Independent t-test model: two means predict <span class="math inline">\(y\)</span>.</p>
<p><span class="math inline">\(y_i = \beta_0 + \beta_1 x_i \qquad \mathcal{H}_0: \beta_1 = 0\)</span></p>
<p>where <span class="math inline">\(x_i\)</span> is an indicator (0 or 1) saying whether data point <span class="math inline">\(i\)</span> was sampled from one or the other group. <a href="https://en.wikipedia.org/wiki/Dummy_variable_(statistics)">Indicator variables (also called “dummy coding”)</a> underly a lot of linear models and we’ll take an aside to see how it works in a minute.</p>
<p><strong>Mann-Whitney U</strong> (also known as <strong>Wilcoxon signed-rank test</strong> for two independent groups) is the same model to a very close approximation, just on the ranks of <span class="math inline">\(x\)</span> and <span class="math inline">\(y\)</span> instead of the actual values:</p>
<p><span class="math inline">\(rank(y_i) = \beta_0 + \beta_1 x_i \qquad \mathcal{H}_0: \beta_1 = 0\)</span></p>
<p>To me, equivalences like this make “non-parametric” statistics much easier to understand. The approximation is appropriate <a href="simulations/simulate_mannwhitney.html">when the sample size is larger than 11 in each group and virtually perfect when N > 30 in each group</a>.</p>
</div>
<div id="dummy" class="section level3">
<h3><span class="header-section-number">5.1.2</span> Theory: Dummy coding</h3>
<p>Dummy coding can be understood visually. The indicator is on the x-axis so data points from the first group are located at <span class="math inline">\(x = 0\)</span> and data points from the second group is located at <span class="math inline">\(x = 1\)</span>. Then <span class="math inline">\(\beta_0\)</span> is the intercept (blue line) and <span class="math inline">\(\beta_1\)</span> is the slope between the two means (red line). Why? Because when <span class="math inline">\(\Delta x = 1\)</span> the slope equals the difference because:</p>
<p><span class="math inline">\(slope = \Delta y / \Delta x = \Delta y / 1 = \Delta y = difference\)</span></p>
<p>Magic! Even categorical differences can be modelled using linear models! It’s a true Swizz army knife.</p>
<div class="fold s">
<pre class="r"><code># Data
N = 20 # Number of data points per group
D_t2 = data.frame(
x = rep(c(0, 1), each=N),
y = c(rnorm_fixed(N, 0.3, 0.3), rnorm_fixed(N, 1.3, 0.3))
)
# Plot
P_t2 = ggplot(D_t2, aes(x=x, y=y)) +
stat_summary(fun.y=mean, geom = "errorbar", aes(ymax = ..y.., ymin = ..y.., color='something'), lwd=2) +
geom_segment(x=-10, xend=10, y=0.3, yend=0.3, lwd=2, aes(color='beta_0')) +
geom_segment(x=0, xend=1, y=0.3, yend=1.3, lwd=2, aes(color='beta_1')) +
scale_color_manual(name=NULL, values=c("blue", "red", "darkblue"), labels=c(bquote(beta[0]*" (group 1 mean)"), bquote(beta[1]*" (slope = difference)"), bquote(beta[0]+beta[1]%.%1*" (group 2 mean)")))
#scale_x_discrete(breaks=c(0.5, 1.5), labels=c('1', '2'))
theme_axis(P_t2, jitter=TRUE, xlim=c(-0.3, 2), legend.position=c(0.53, 0.08))</code></pre>
<p><img src="index_files/figure-html/unnamed-chunk-27-1.png" width="576" style="display: block; margin: auto;" /></p>
</div>
</div>
<div id="dummy2" class="section level3">
<h3><span class="header-section-number">5.1.3</span> Theory: Dummy coding (continued)</h3>
<p>If you feel like you get dummy coding now, just skip ahead to the next section. Here is a more elaborate explanation of dummy coding:</p>
<p>If a data point was sampled from the first group, i.e., when <span class="math inline">\(x_i = 0\)</span>, the model simply becomes <span class="math inline">\(y = \beta_0 + \beta_1 \cdot 0 = \beta_0\)</span>. In other words, the model predicts that that data point is <span class="math inline">\(beta_0\)</span>. It turns out that the <span class="math inline">\(\beta\)</span> which best predicts a set of data points is the <em>mean</em> of those data points, so <span class="math inline">\(\beta_0\)</span> is the mean of group 1.</p>
<p>On the other hand, data points sampled from the second group would have <span class="math inline">\(x_i = 1\)</span> so the model becomes <span class="math inline">\(y_i = \beta_0 + \beta_1\cdot 1 = \beta_0 + \beta_1\)</span>. In other words, we add <span class="math inline">\(\beta_1\)</span> to “shift” from the mean of the first group to the mean of the second group. Thus <span class="math inline">\(\beta_1\)</span> becomes the <em>mean difference</em> between the groups.</p>
<p>As an example, say group 1 is 25 years old (<span class="math inline">\(\beta_0 = 25\)</span>) and group 2 is 28 years old (<span class="math inline">\(\beta_1 = 3\)</span>), then the model for a person in group 1 is <span class="math inline">\(y = 25 + 3 \cdot 0 = 25\)</span> and the model for a person in group 2 is <span class="math inline">\(y = 25 + 3 \cdot 1 = 28\)</span>.</p>
<p>Hooray, it works! For first-timers it takes a few moments to understand dummy coding, but you only need to know addition and multiplication to get there!</p>
</div>
<div id="r-code-independent-t-test" class="section level3">
<h3><span class="header-section-number">5.1.4</span> R code: independent t-test</h3>
<p>As a reminder, when we write <code>y ~ 1 + x</code> in R, it is shorthand for <span class="math inline">\(y = \beta_0 \cdot 1 + \beta_1 \cdot x\)</span> and R goes on computing the <span class="math inline">\(\beta\)</span>s for you. Thus <code>y ~ 1 + x</code> is the R-way of writing <span class="math inline">\(y = a \cdot x + b\)</span>.</p>
<p>Notice the identical <code>t</code>, <code>df</code>, <code>p</code>, and estimates. We can get the confidence interval by running <code>confint(lm(...))</code>.</p>
<pre class="r"><code># Built-in independent t-test on wide data
a = t.test(y, y2, var.equal = TRUE)
# Be explicit about the underlying linear model by hand-dummy-coding:
group_y2 = ifelse(group == 'y2', 1, 0) # 1 if group == y2, 0 otherwise
b = lm(value ~ 1 + group_y2) # Using our hand-made dummy regressor
# Note: We could also do the dummy-coding in the model
# specification itself. Same result.
c = lm(value ~ 1 + I(group=='y2'))</code></pre>
<p>Results:</p>
<div id="htmlwidget-53b434c497272bb7da0e" style="width:100%;height:auto;" class="datatables html-widget"></div>
<script type="application/json" data-for="htmlwidget-53b434c497272bb7da0e">{"x":{"filter":"none","data":[["t.test","lm"],[0.3,0.3],[0.2,0.2],[0.5729,0.5729],[98,98],[-0.5016,-0.5016],[0.9016,0.9016]],"container":"<table class=\"display\">\n <thead>\n <tr>\n <th>model<\/th>\n <th>mean_y<\/th>\n <th>difference<\/th>\n <th>p.value<\/th>\n <th>df<\/th>\n <th>conf.low<\/th>\n <th>conf.high<\/th>\n <\/tr>\n <\/thead>\n<\/table>","options":{"searching":false,"lengthChange":false,"ordering":false,"autoWidth":true,"bPaginate":false,"bInfo":false,"paging":false,"columnDefs":[{"className":"dt-right","targets":[1,2,3,4,5,6]}],"order":[],"orderClasses":false}},"evals":[],"jsHooks":[]}</script>
<div class="fold o">
<pre><code>##
## Two Sample t-test
##
## data: y and y2
## t = -0.56569, df = 98, p-value = 0.5729
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.9016152 0.5016152
## sample estimates:
## mean of x mean of y
## 0.3 0.5
##
##
## Call:
## lm(formula = value ~ 1 + group_y2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.5211 -1.1259 -0.2124 0.9151 4.9268
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.3000 0.2500 1.200 0.233
## group_y2 0.2000 0.3536 0.566 0.573
##
## Residual standard error: 1.768 on 98 degrees of freedom
## Multiple R-squared: 0.003255, Adjusted R-squared: -0.006916
## F-statistic: 0.32 on 1 and 98 DF, p-value: 0.5729
##
##
## Call:
## lm(formula = value ~ 1 + I(group == "y2"))
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.5211 -1.1259 -0.2124 0.9151 4.9268
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.3000 0.2500 1.200 0.233
## I(group == "y2")TRUE 0.2000 0.3536 0.566 0.573
##
## Residual standard error: 1.768 on 98 degrees of freedom
## Multiple R-squared: 0.003255, Adjusted R-squared: -0.006916
## F-statistic: 0.32 on 1 and 98 DF, p-value: 0.5729</code></pre>
</div>
</div>
<div id="r-code-mann-whitney-u" class="section level3">
<h3><span class="header-section-number">5.1.5</span> R code: Mann-Whitney U</h3>
<pre class="r"><code># Wilcoxon / Mann-Whitney U
a = wilcox.test(y, y2)
# As linear model with our dummy-coded group_y2:
b = lm(rank(value) ~ 1 + group_y2) # compare to linear model above</code></pre>
<div id="htmlwidget-0b05f5a4b604340b9f1a" style="width:100%;height:auto;" class="datatables html-widget"></div>
<script type="application/json" data-for="htmlwidget-0b05f5a4b604340b9f1a">{"x":{"filter":"none","data":[["wilcox.test","lm"],[0.7907,0.7896],[null,1.56]],"container":"<table class=\"display\">\n <thead>\n <tr>\n <th>model<\/th>\n <th>p.value<\/th>\n <th>rank_diff<\/th>\n <\/tr>\n <\/thead>\n<\/table>","options":{"searching":false,"lengthChange":false,"ordering":false,"autoWidth":true,"bPaginate":false,"bInfo":false,"paging":false,"columnDefs":[{"className":"dt-right","targets":[1,2]}],"order":[],"orderClasses":false}},"evals":[],"jsHooks":[]}</script>
<div class="fold o">
<pre><code>##
## Wilcoxon rank sum test with continuity correction
##
## data: y and y2
## W = 1211, p-value = 0.7907
## alternative hypothesis: true location shift is not equal to 0
##
##
## Call:
## lm(formula = rank(value) ~ 1 + group_y2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -48.72 -24.64 -0.78 24.22 48.72
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 49.720 4.122 12.061 <2e-16 ***
## group_y2 1.560 5.830 0.268 0.79
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 29.15 on 98 degrees of freedom
## Multiple R-squared: 0.0007302, Adjusted R-squared: -0.009466
## F-statistic: 0.07161 on 1 and 98 DF, p-value: 0.7896</code></pre>
</div>
</div>
</div>
<div id="welch" class="section level2">
<h2><span class="header-section-number">5.2</span> Welch’s t-test</h2>
<p>This is identical to the (Student’s) <a href="#t2">independent t-test</a> above except that Student’s assumes identical variances and <strong>Welch’s t-test</strong> does not. So the linear model is the same and the trick is in the variances, which I won’t go further into here.</p>
<pre class="r"><code># Built-in
a = t.test(y, y2, var.equal=FALSE)
# As linear model with per-group variances
b = nlme::gls(value ~ 1 + group_y2, weights = nlme::varIdent(form=~1|group), method="ML")</code></pre>
<p>Results:</p>
<div id="htmlwidget-accd39a4468ae76de8c1" style="width:100%;height:auto;" class="datatables html-widget"></div>
<script type="application/json" data-for="htmlwidget-accd39a4468ae76de8c1">{"x":{"filter":"none","data":[["t.test","gls"],[0.3,0.3],[0.2,0.2],[0.573,0.5729],[-0.5657,-0.5657],[-0.9023,-0.893],[0.5023,0.493]],"container":"<table class=\"display\">\n <thead>\n <tr>\n <th>model<\/th>\n <th>mean_y<\/th>\n <th>mean_diff<\/th>\n <th>p.value<\/th>\n <th>t<\/th>\n <th>conf.low<\/th>\n <th>conf.high<\/th>\n <\/tr>\n <\/thead>\n<\/table>","options":{"searching":false,"lengthChange":false,"ordering":false,"autoWidth":true,"bPaginate":false,"bInfo":false,"paging":false,"columnDefs":[{"className":"dt-right","targets":[1,2,3,4,5,6]}],"order":[],"orderClasses":false}},"evals":[],"jsHooks":[]}</script>
<div class="fold o">
<pre><code>##
## Welch Two Sample t-test
##
## data: y and y2
## t = -0.56569, df = 90.875, p-value = 0.573
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.9023034 0.5023034
## sample estimates:
## mean of x mean of y
## 0.3 0.5
##
## Generalized least squares fit by maximum likelihood
## Model: value ~ 1 + group_y2
## Data: NULL
## AIC BIC logLik
## 399.6287 410.0493 -195.8143
##
## Variance function:
## Structure: Different standard deviations per stratum