-
Notifications
You must be signed in to change notification settings - Fork 3
Expand file tree
/
Copy pathchapter_04.html
More file actions
1535 lines (1412 loc) · 66.7 KB
/
Copy pathchapter_04.html
File metadata and controls
1535 lines (1412 loc) · 66.7 KB
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>
<head>
<meta charset="utf-8" />
<meta name="generator" content="pandoc" />
<meta http-equiv="X-UA-Compatible" content="IE=EDGE" />
<title>Chapter 4: SoFR</title>
<script src="site_libs/header-attrs-2.24/header-attrs.js"></script>
<script src="site_libs/jquery-3.6.0/jquery-3.6.0.min.js"></script>
<meta name="viewport" content="width=device-width, initial-scale=1" />
<link href="site_libs/bootstrap-3.3.5/css/cosmo.min.css" rel="stylesheet" />
<script src="site_libs/bootstrap-3.3.5/js/bootstrap.min.js"></script>
<script src="site_libs/bootstrap-3.3.5/shim/html5shiv.min.js"></script>
<script src="site_libs/bootstrap-3.3.5/shim/respond.min.js"></script>
<style>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;}
code {color: inherit; background-color: rgba(0, 0, 0, 0.04);}
pre:not([class]) { background-color: white }</style>
<script src="site_libs/jqueryui-1.13.2/jquery-ui.min.js"></script>
<link href="site_libs/tocify-1.9.1/jquery.tocify.css" rel="stylesheet" />
<script src="site_libs/tocify-1.9.1/jquery.tocify.js"></script>
<script src="site_libs/navigation-1.1/tabsets.js"></script>
<link href="site_libs/highlightjs-9.12.0/default.css" rel="stylesheet" />
<script src="site_libs/highlightjs-9.12.0/highlight.js"></script>
<link href="site_libs/font-awesome-6.4.2/css/all.min.css" rel="stylesheet" />
<link href="site_libs/font-awesome-6.4.2/css/v4-shims.min.css" rel="stylesheet" />
<link rel="stylesheet" href="https://maxcdn.bootstrapcdn.com/font-awesome/latest/css/font-awesome.min.css" />
<script defer src="https://use.fontawesome.com/releases/v5.0.3/js/all.js"></script>
<script defer src="https://use.fontawesome.com/releases/v5.0.0/js/v4-shims.js"></script>
<!-- Global site tag (gtag.js) - Google Analytics
<script async src="https://www.googletagmanager.com/gtag/js?id=UA-151578452-1"></script>
<script>
window.dataLayer = window.dataLayer || [];
function gtag(){dataLayer.push(arguments);}
gtag('js', new Date());
gtag('config', 'UA-151578452-1');
</script>
-->
<style type="text/css">
code{white-space: pre-wrap;}
span.smallcaps{font-variant: small-caps;}
span.underline{text-decoration: underline;}
div.column{display: inline-block; vertical-align: top; width: 50%;}
div.hanging-indent{margin-left: 1.5em; text-indent: -1.5em;}
ul.task-list{list-style: none;}
</style>
<style type="text/css">code{white-space: pre;}</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>
<link rel="stylesheet" href="styles.css" type="text/css" />
<style type = "text/css">
.main-container {
max-width: 940px;
margin-left: auto;
margin-right: auto;
}
img {
max-width:100%;
}
.tabbed-pane {
padding-top: 12px;
}
.html-widget {
margin-bottom: 20px;
}
button.code-folding-btn:focus {
outline: none;
}
summary {
display: list-item;
}
details > summary > p:only-child {
display: inline;
}
pre code {
padding: 0;
}
</style>
<style type="text/css">
.dropdown-submenu {
position: relative;
}
.dropdown-submenu>.dropdown-menu {
top: 0;
left: 100%;
margin-top: -6px;
margin-left: -1px;
border-radius: 0 6px 6px 6px;
}
.dropdown-submenu:hover>.dropdown-menu {
display: block;
}
.dropdown-submenu>a:after {
display: block;
content: " ";
float: right;
width: 0;
height: 0;
border-color: transparent;
border-style: solid;
border-width: 5px 0 5px 5px;
border-left-color: #cccccc;
margin-top: 5px;
margin-right: -10px;
}
.dropdown-submenu:hover>a:after {
border-left-color: #adb5bd;
}
.dropdown-submenu.pull-left {
float: none;
}
.dropdown-submenu.pull-left>.dropdown-menu {
left: -100%;
margin-left: 10px;
border-radius: 6px 0 6px 6px;
}
</style>
<script type="text/javascript">
// manage active state of menu based on current page
$(document).ready(function () {
// active menu anchor
href = window.location.pathname
href = href.substr(href.lastIndexOf('/') + 1)
if (href === "")
href = "index.html";
var menuAnchor = $('a[href="' + href + '"]');
// mark the anchor link active (and if it's in a dropdown, also mark that active)
var dropdown = menuAnchor.closest('li.dropdown');
if (window.bootstrap) { // Bootstrap 4+
menuAnchor.addClass('active');
dropdown.find('> .dropdown-toggle').addClass('active');
} else { // Bootstrap 3
menuAnchor.parent().addClass('active');
dropdown.addClass('active');
}
// Navbar adjustments
var navHeight = $(".navbar").first().height() + 15;
var style = document.createElement('style');
var pt = "padding-top: " + navHeight + "px; ";
var mt = "margin-top: -" + navHeight + "px; ";
var css = "";
// offset scroll position for anchor links (for fixed navbar)
for (var i = 1; i <= 6; i++) {
css += ".section h" + i + "{ " + pt + mt + "}\n";
}
style.innerHTML = "body {" + pt + "padding-bottom: 40px; }\n" + css;
document.head.appendChild(style);
});
</script>
<!-- tabsets -->
<style type="text/css">
.tabset-dropdown > .nav-tabs {
display: inline-table;
max-height: 500px;
min-height: 44px;
overflow-y: auto;
border: 1px solid #ddd;
border-radius: 4px;
}
.tabset-dropdown > .nav-tabs > li.active:before, .tabset-dropdown > .nav-tabs.nav-tabs-open:before {
content: "\e259";
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: "\e258";
font-family: 'Glyphicons Halflings';
border: none;
}
.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;
background-color: transparent;
}
.tabset-dropdown > .nav-tabs.nav-tabs-open > li {
display: block;
float: none;
}
.tabset-dropdown > .nav-tabs > li {
display: none;
}
</style>
<!-- code folding -->
<style type="text/css">
#TOC {
margin: 25px 0px 20px 0px;
}
@media (max-width: 768px) {
#TOC {
position: relative;
width: 100%;
}
}
@media print {
.toc-content {
/* see https://github.com/w3c/csswg-drafts/issues/4434 */
float: right;
}
}
.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;
}
.tocify .list-group-item {
border-radius: 0px;
}
</style>
</head>
<body>
<div class="container-fluid main-container">
<!-- setup 3col/9col grid for toc_float and main content -->
<div class="row">
<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="navbar navbar-default navbar-fixed-top" role="navigation">
<div class="container">
<div class="navbar-header">
<button type="button" class="navbar-toggle collapsed" data-toggle="collapse" data-bs-toggle="collapse" data-target="#navbar" data-bs-target="#navbar">
<span class="icon-bar"></span>
<span class="icon-bar"></span>
<span class="icon-bar"></span>
</button>
<a class="navbar-brand" href="index.html">FDA with R</a>
</div>
<div id="navbar" class="navbar-collapse collapse">
<ul class="nav navbar-nav">
</ul>
<ul class="nav navbar-nav navbar-right">
<li>
<a href="about_authors.html">About the Authors</a>
</li>
<li class="dropdown">
<a href="#" class="dropdown-toggle" data-toggle="dropdown" role="button" data-bs-toggle="dropdown" aria-expanded="false">
Datasets
<span class="caret"></span>
</a>
<ul class="dropdown-menu" role="menu">
<li>
<a href="dataset_nhanes.html">NHANES</a>
</li>
<li>
<a href="dataset_covid19.html">COVID-19</a>
</li>
<li>
<a href="dataset_cd4.html">CD4</a>
</li>
<li>
<a href="dataset_content.html">Content</a>
</li>
</ul>
</li>
<li class="dropdown">
<a href="#" class="dropdown-toggle" data-toggle="dropdown" role="button" data-bs-toggle="dropdown" aria-expanded="false">
Chapters
<span class="caret"></span>
</a>
<ul class="dropdown-menu" role="menu">
<li>
<a href="chapter_01.html">Chapter 1</a>
</li>
<li>
<a href="chapter_02.html">Chapter 2</a>
</li>
<li>
<a href="chapter_03.html">Chapter 3: FPCA</a>
</li>
<li>
<a href="chapter_04.html">Chapter 4: SoFR</a>
</li>
<li>
<a href="chapter_05.html">Chapter 5: FoSR</a>
</li>
<li>
<a href="chapter_06.html">Chapter 6: FoFR</a>
</li>
<li>
<a href="chapter_07.html">Chapter 7</a>
</li>
<li>
<a href="chapter_08.html">Chapter 8</a>
</li>
<li>
<a href="chapter_09.html">Chapter 9</a>
</li>
</ul>
</li>
<li>
<a href="scripts.html">Scripts</a>
</li>
<li>
<a href="https://github.com/FDAwithR">
<span class="fa fa-github"></span>
</a>
</li>
</ul>
</div><!--/.nav-collapse -->
</div><!--/.container -->
</div><!--/.navbar -->
<div id="header">
<h1 class="title toc-ignore">Chapter 4: SoFR</h1>
</div>
<p>This book Chapter and the code in this page considers functions as
predictors in models with scalar outcomes. We focus on the linear
scalar-on-function regression (SoFR), starting the general motivation
and building intuition using exploratory analyses and careful
interpretation of coefficients. Methods using unpenalized basis
expansions are implemented in traditional software, while estimation and
inference for SoFR using penalized splines is conducted using the
<code>refund</code> and <code>mgcv</code> packages.</p>
<div id="motivation-and-eda" class="section level2">
<h2>Motivation and EDA</h2>
<p>Much of this page will use data from the NHANES to illustrate
techniques for modeling scalar outcomes using functional predictors. In
particular, we will use MIMS profiles as functional predictors of body
mass index (BMI) as a continuous outcome. We will begin with a simple
approach and show how this is related to more general linear
scalar-on-function regression models. For each participant, we will
obtain the average MIMS value in each of 12 two-hour bins and then use
these bin averages as predictors in a standard linear regression model.
This has the advantage of allowing some degree of flexibility in the
association between MIMS values and BMI over the course of the day,
while not requiring non-standard modeling techniques or interpretations.
Indeed, in our experience this kind of approach can be a good starting
point in collaborative settings or serve as a useful check on results
from more complex approaches.</p>
<p>The code chunk below imports and organizes the processed NHANES data
that will be used in this page. After importing the data, we create a
new variable which indicates whether the participant died within two
years of their inclusion in the study; retain only those variables in
the processed NHANES data that will be relevant for this example; rename
<code>MIMS</code> and <code>MIMS_sd</code> to include the suffix
<code>_mat</code> to indicate these variables are stored as matrices;
restrict the dataset to participants 25 years of age or older; and
ensure that the resulting dataframe has the class
<code>tibble</code>.</p>
<pre class="r"><code>nhanes_df =
readRDS(
here::here("data", "nhanes_fda_with_r.rds")) %>%
mutate(
death_2yr = ifelse(event == 1 & time <= 24, 1, 0)) %>%
select(
SEQN, BMI, age, gender, death_2yr,
MIMS_mat = MIMS, MIMS_sd_mat = MIMS_sd) %>%
filter(age >= 25) %>%
drop_na(BMI) %>%
tibble()</code></pre>
<p>In the next code chunk, we convert <code>MIMS_mat</code> and
<code>MIMS_mat_sd</code> to <code>tidyfun</code> objects using
<code>tfd()</code>using the <code>arg</code> argument in
<code>tfd()</code> to be explicit about the grid over which functions
are observed.</p>
<pre class="r"><code>nhanes_df =
nhanes_df %>%
mutate(
MIMS_tf = matrix(MIMS_mat, ncol = 1440),
MIMS_tf = tfd(MIMS_tf, arg = seq(1/60, 24, length = 1440)),
MIMS_sd_tf = matrix(MIMS_sd_mat, ncol = 1440),
MIMS_sd_tf = tfd(MIMS_sd_tf, arg = seq(1/60, 24, length = 1440)))</code></pre>
<p>The next code chunk contains two components. The first component
creates a new data frame containing average MIMS values in two-hour bins
by computing the rolling mean of each <code>MIMS_tf</code> observation
with a bin width of 120 minutes, and then evaluating that rolling mean
at hours <span class="math inline">\(1, 3, ..., 23\)</span>. The result
is saved as <code>MIMS_binned</code>, and for the next step only
<code>BMI</code> and <code>MIMS_binned</code> are retained.</p>
<p>The second component of this code chunk fits the regression of
<code>BMI</code> on these bin averages. The <code>tf_spread()</code>
function produces a wide-format dataframe with columns corresponding to
each bin average in the <code>MIMS_binned</code> variable, and the call
to <code>lm()</code> regresses <code>BMI</code> on all of these
averages.</p>
<pre class="r"><code>nhanes_bin_df =
nhanes_df %>%
mutate(
MIMS_binned =
tf_smooth(MIMS_tf, method = "rollmean", k = 120, align = "center"),
MIMS_binned = tfd(MIMS_binned, arg = seq(1, 23, by = 2))) %>%
select(BMI, MIMS_binned)
## setting fill = 'extend' for start/end values.
## Warning: There was 1 warning in `mutate()`.
## ℹ In argument: `MIMS_binned = tf_smooth(MIMS_tf, method = "rollmean", k = 120,
## align = "center")`.
## Caused by warning:
## ! non-equidistant arg-values in 'MIMS_tf' ignored by rollmean.
fit_binned =
lm(BMI ~ .,
data = nhanes_bin_df %>% tf_spread(MIMS_binned))</code></pre>
<p>We now show the binned predictors and the resulting coefficients
using plotting tools in <a
href="https://tidyfun.github.io/tidyfun/"><code>tidyfun</code></a>. The
first plot generated in the code chunk below shows the
<code>MIMS_binned</code> variable for the first 500 rows (other data
points are omitted to prevent overplotting). The second plot shows the
coefficients for each bin averaged MIMS values. We create this plot by
tidying the model fit stored in <code>fit_binned</code> and omitting the
intercept term. An <code>hour</code> variable is then created by
manipulating the coefficient names, and upper and lower 95% confidence
bounds for each hour are obtained by adding and subtracting 1.96 times
the standard error from the estimates. We plot the estimates as lines
and points, and add error bars for the confidence intervals. The two
panels are combined using the <a
href="https://github.com/thomasp85/patchwork"><code>patchwork</code></a>
package.</p>
<p>The binned MIMS profiles show expected diurnal patterns of activity,
where there is generally lower activity in the night and higher activity
in the day. The results of the regression using bin-averaged MIMS values
as predictors for BMI (right panel) are consistent with these observed
trends. Coefficients for bin averages during the day are generally below
zero and some (2-hour intervals between 8-10AM and 6-8PM) are
statistically significant.</p>
<pre class="r"><code>ggp_binned =
nhanes_bin_df %>%
slice(1:500) %>%
ggplot(aes(y = MIMS_binned, color = BMI)) +
geom_spaghetti() +
geom_meatballs() +
scale_x_continuous(
breaks = seq(0, 24, length = 5),
labels = str_c(seq(0, 24, length = 5), ":00")) +
labs(x = "Time of day (hours)", y = "Binned MIMS")
ggp_coefs =
fit_binned %>%
broom::tidy() %>%
filter(term != "(Intercept)") %>%
mutate(
hour = str_replace(term, "MIMS_binned_", ""),
hour = as.numeric(hour),
ub = estimate + 1.96 * std.error,
lb = estimate - 1.96 * std.error) %>%
ggplot(aes(x = hour, y = estimate)) +
geom_point() + geom_path() +
geom_errorbar(aes(ymin = lb, ymax = ub), width = .5) +
scale_x_continuous(
breaks = seq(0, 24, length = 5),
labels = str_c(seq(0, 24, length = 5), ":00")) +
labs(x = "Time of day (hours)", y = "Coefficient")
ggp_binned + ggp_coefs</code></pre>
<p><img src="chapter_04_files/figure-html/unnamed-chunk-3-1.png" width="90%" /></p>
<p>We will motivate a shift to scalar-on-function regression by noting
that the model using bin averages can be expressed in terms of a
functional coefficient expressed as a step function and functional
predictors by integrating over their product. This is motivated by the
following expression, with technical details provided in the book
Chapter text: <span class="math display">\[\begin{equation}
\sum_{b=1}^{12} \beta_{b} \overline{X}_{ib} \approx \int_{0}^{24}
\beta^{step}(s) X_{i}(s) \, ds.
\label{eq:ch4_step_int}
\end{equation}\]</span></p>
<p>The code chunk below creates a dataframe that contains the step
coefficient function defined above.</p>
<pre class="r"><code>stepfun_coef_df =
fit_binned %>%
broom::tidy() %>%
filter(term != "(Intercept)") %>%
select(estimate, std.error) %>%
slice(rep(1:12, each = 120)) %>%
mutate(
method = "Step",
estimate = .5 * estimate,
arg = seq(1/60, 24, length = 1440)) %>%
tf_nest(.id = method)</code></pre>
<p>A plot showing the complete (not binned) <code>MIMS_tf</code>
trajectories alongside the step coefficient function is shown below.
Comparing this plot with one above, in which bin averages are shown
alongside regression coefficients, is intended to emphasize that these
approaches are identical up a to a re-scaling of the regression
parameters. Connecting the bin average approach to a truly functional
coefficient is an intuitive starting point for the more flexible linear
SoFR models considered next.</p>
<pre class="r"><code>ggp_stepfun_mims =
nhanes_df %>%
slice(1:500) %>%
ggplot(aes(y = MIMS_tf, color = BMI)) +
geom_spaghetti() +
scale_x_continuous(
breaks = seq(0, 24, length = 5),
labels = str_c(seq(0, 24, length = 5), ":00")) +
labs(x = "Time of day (hours)", y = "Coefficient")
ggp_stepfun_coef =
stepfun_coef_df %>%
ggplot(aes(y = estimate)) +
geom_spaghetti(linewidth = 1.1, alpha = .75) +
scale_x_continuous(
breaks = seq(0, 24, length = 5),
labels = str_c(seq(0, 24, length = 5), ":00")) +
labs(x = "Time of day (hours)", y = "Coefficient")
ggp_stepfun_mims + ggp_stepfun_coef</code></pre>
<p><img src="chapter_04_files/figure-html/unnamed-chunk-5-1.png" width="90%" /></p>
</div>
<div id="simple-linear-scalar-on-function" class="section level2">
<h2>“Simple Linear” Scalar-on-Function</h2>
<p>We now introduce the linear scalar-on-function regression model in
which there is only one functional predictor and no scalar covariates.
This is analogous to ``simple” linear regression, and will be useful for
introducing key concepts in interpretation, estimation, and inference
for models with functional predictors. Later sections will consider
extensions of this approach.</p>
<div id="model-specification-and-interpretation" class="section level3">
<h3>Model specification and interpretation</h3>
<p>For participants <span class="math inline">\(i = 1, \dots,
n\)</span>, let <span class="math inline">\(Y_i\)</span> be a scalar
response of interest and <span class="math inline">\(X_{i}: S\rightarrow
\mathbb{R}\)</span> be a functional predictor observed over the domain
<span class="math inline">\(S\)</span>. The simple linear
scalar-on-function regression model is <span
class="math display">\[\begin{equation}
Y_i=\beta_0 + \int_{S} \beta_1 (s)X_{i}(s)\, ds + \epsilon_i
\end{equation}\]</span> where <span
class="math inline">\(\beta_0\)</span> is a population-level scalar
intercept, <span class="math inline">\(\beta_1(s): S\rightarrow
\mathbb{R}\)</span> is the functional coefficient of interest, and <span
class="math inline">\(\epsilon_i\)</span> is a residual with mean zero
and variance <span class="math inline">\(\sigma^2_{\epsilon}\)</span>.
This model specification generalizes the specific case where a
coefficient function constrained to be a step function approximated a
regression model based on bin averages. The coefficient function now can
be more flexible, but the interpretation is analogous: <span
class="math inline">\(\beta_1(s)\)</span> defines the weight given to
predictor functions at each <span class="math inline">\(s \in
S\)</span>, and the product of <span
class="math inline">\(X_{i}(s)\)</span> and <span
class="math inline">\(\beta_1(s)\)</span> is integrated to obtain the
overall contribution of the functional term in the model. At each time
point, <span class="math inline">\(\beta(s)\)</span> is the effect on
the outcome for a one unit increase in <span
class="math inline">\(X_{i}(s)\)</span> when all other values <span
class="math inline">\(X_{i}(s')\)</span> <span
class="math inline">\(s'\in S, s'\neq s\)</span> are kept
unchanged. This interpretation is admittedly somewhat awkward, but
unavoidable when regressing a scalar outcome on a functional predictor.
Note that the coefficients adjust for effects of all other <span
class="math inline">\(s'\in S\)</span>.</p>
<p>The innovation in scalar-on-function regression, compared to
nonfunctional models, is a coefficient function that integrates with
covariate functions to produce scalar terms in the linear predictor. The
corresponding challenge is developing an estimation strategy that
minimizes <span class="math display">\[\begin{equation}\min_{\beta_0,
\beta_{1}(s)} \sum_{i = 1}^{n} \left\{ Y_i - \beta_0 - \int_{S} \beta_1
(s)X_{i}(s)\, ds\right\}^2
\label{eq:ch4_sofr_ss}
\end{equation}\]</span> in a way that is flexible and computationally
feasible.</p>
</div>
<div id="parametric-estimation-of-the-coefficient-function"
class="section level3">
<h3>Parametric estimation of the coefficient function</h3>
<p>Our first approach expands the coefficient function <span
class="math inline">\(\beta_1(s)\)</span> using a relatively
low-dimensional basis expansion; doing so leads to a more familiar
setting in which scalar basis coefficients are the target of estimation
and inference. Specifically, let <span class="math inline">\(\beta_1(s)
= \sum_{k = 1}^{K}\beta_{1k}B_{k}(s)\)</span> where <span
class="math inline">\(B_1(s), \dots, B_{K}(s)\)</span> is a collection
of basis functions. Substituting this expansion into the integral term
gives <span class="math display">\[\begin{equation}
\begin{array}{lll}
E[Y_i] &=& \beta_0 + \int_S \beta_1 (s)X_{i}(s)\, ds \\
&=& \beta_0 + \sum_{k=1}^{K} \left[ \int_{S}
B_{k}(s)X_{i}(s)\, ds\right] \beta_{1k} \\
&=& \beta_0 + \mathbf{C}^t_{i} \boldsymbol{\beta}_{1}
\end{array}
\end{equation}\]</span> where <span class="math inline">\(C_{ik} =
\int_{S} B_{k}(s)X_{i}(s)\, ds\)</span>, <span
class="math inline">\(\mathbf{C}_i = [C_{i1}, \dots,
C_{iK}]^{t}\)</span>, and <span
class="math inline">\(\boldsymbol{\beta}_1=(\beta_{11},\ldots,\beta_{1K})^t\)</span>
is the vector of basis coefficients. The result of the basis expansion
for the coefficient function, therefore, is a recognizable multiple
linear regression with carefully-defined scalar covariates and
corresponding coefficients. Specifically, let <span
class="math inline">\(\mathbf{y}=(y_1,\ldots,y_n)^t\)</span>, the matrix
<span class="math inline">\(\mathbf{X}\)</span> be constructed by
row-stacking vectors <span class="math inline">\([1, C_{i1}, \dots,
C_{iK}]\)</span>, and <span class="math inline">\(\boldsymbol{\beta} =
[\beta_{0}, \beta_{11}, \dots, \beta_{1K}]^{t}\)</span> be the vector of
regression coefficients including the population intercept and spline
coefficients. This suggests a standard OLS approach to estimating <span
class="math inline">\(\boldsymbol{\beta}\)</span>. Note that functional
predictors are actually observed over a finite grid, and the definite
integrals that define the <span class="math inline">\(C_{ik}\)</span>
are in practice estimated using numeric quadrature.</p>
<p>Many options for the basis have been considered in the expansive
literature for SoFR. To illustrate the ideas in this Section, we start
with a quadratic basis and obtain the corresponding estimate of the
coefficient function <span class="math inline">\(\beta_{1}(s)\)</span>.
We define the basis <span class="math display">\[\begin{equation}
\left\{\begin{array}{lll}
B_{1}(s) = 1\;,\\
B_{2}(s) = s\;,\\
B_{3}(s) = s^2\;, \\
\end{array}\right.
\label{eq:ch4_quad_basis}
\end{equation}\]</span> and, given this, obtain scalar predictors <span
class="math inline">\(C_{ik}\)</span> that can be used in standard
linear model software. The basis expansion includes an intercept term,
which should not be confused with the model’s intercept, <span
class="math inline">\(\beta_0\)</span>. The intercept in the basis
expansion allows the coefficient function <span
class="math inline">\(\beta_1(s)\)</span> to shift as needed, while the
population intercept is the expected value of the response when the
predictor function is zero over the entire domain.</p>
<p>Continuing to focus on BMI as an outcome and MIMS as a functional
predictor, the code chunk below defines the quadratic basis and obtains
the numeric integrals in the <span
class="math inline">\(C_{ik}\)</span>. The basis matrix is defined in
terms of <code>arg</code> and given appropriate column names. We
construct the data frame <code>num_int_df</code> which contains the
necessary numeric integrals. We retain the row names of the matrix
product in the resulting dataframe, and convert this to a numeric
variable for consistency with <code>nhanes_df</code>.</p>
<pre class="r"><code>epoch_arg = seq(1/60, 24, length = 1440)
B = cbind(1, epoch_arg, epoch_arg^2)
colnames(B) = c("int", "lin", "quad")
num_int_df =
as_tibble(
(nhanes_df$MIMS_mat %*% B) * (1/ 60),
rownames = "SEQN") %>%
mutate(SEQN = as.numeric(SEQN))</code></pre>
<p>The next code chunk implements the regression and processes the
results. We first define a new data frame, <code>nhanes_quad_df</code>,
that contains variables relevant for the scalar-on-function regression
of BMI on MIMS trajectories using SoFR and expand the coefficient
function <span class="math inline">\(\beta_1(s)\)</span> in terms of the
quadratic basis defined in the previous code chunk. This is created by
joining two previously defined dataframes, and keeping only
<code>BMI</code> and the columns corresponding to the numeric integrals
<span class="math inline">\(C_{ik}\)</span>. Using
<code>nhanes_quad_df</code>, we fit a linear regression of BMI on the
<span class="math inline">\(C_{ik}\)</span>; the formula specification
includes a population intercept to reiterate that the model’s intercept
<span class="math inline">\(\beta_0\)</span> is distinct from the basis
expansion’s intercept, which appears in <span
class="math inline">\(C_{i1}\)</span>. Finally, we combine the
coefficient estimates in <code>fit_quad</code> with the basis matrix to
obtain the estimate of the coefficient function. We compute the matrix
product of <code>B</code> and the coefficients of <code>fit_quad</code>
(omitting the population intercept), and convert the result to a
<code>tidyfun</code> object. The coefficient function is stored in a
data frame called <code>quad_coef_df</code>, along with a variable
<code>method</code> with the value <code>quad</code>.</p>
<pre class="r"><code>nhanes_quad_df =
left_join(nhanes_df, num_int_df, by = "SEQN") %>%
select(BMI, int, lin, quad)
fit_quad =
nhanes_quad_df %>%
lm(BMI ~ 1 + int + lin + quad, data = .)
quad_coef_df =
tibble(
method = "Quadratic",
estimate = tfd(t(B %*% coef(fit_quad)[-1]), arg = epoch_arg))</code></pre>
<p>This general strategy for estimating coefficient functions can be
readily adapted to other basis choices. The next code defines a cubic
B-spline basis with eight degrees of freedom; this is more flexible than
the quadratic basis, while also ensuring a degree of smoothness that is
absent from the stepwise estimate of the coefficient function. Once the
basis is defined, the remaining steps in the code chunk below mirror
those used to estimate the coefficient function using a quadratic basis,
with a small number of minor changes. The basis is generated using the
<code>bs()</code> function in the <code>splines</code> package, and
there are now eight basis functions instead of three. There is a
corresponding increase in the number of columns in
<code>num_int_df</code> and for convenience we write the formula in the
<code>lm()</code> call as <code>BMI ~ 1 + .</code> instead of listing
columns individually. The final step in this code chunk constructs the
estimated coefficient function by multiplying the matrix of basis
functions evaluated over <span class="math inline">\(\mathbf{s}\)</span>
by the vector of B-spline coefficients; the result is stored in a data
frame called <code>bspline_coef_df</code>, now with a variable
<code>method</code> taking the value <code>B-Spline</code>. The
similarity between this model fitting and the one using a quadratic
basis is intended to emphasize that the basis expansion approach to
fitting a linear SoFR can be easy to implement for a broad range of
basis choices.</p>
<pre class="r"><code>B_bspline = splines::bs(epoch_arg, df = 8, intercept = TRUE)
colnames(B_bspline) = str_c("BS_", 1:8)
num_int_df =
as_tibble(
(nhanes_df$MIMS_mat %*% B_bspline) * (1/ 60),
rownames = "SEQN") %>%
mutate(SEQN = as.numeric(SEQN))
nhanes_bspline_df =
left_join(nhanes_df, num_int_df, by = "SEQN") %>%
select(BMI, BS_1:BS_8)
fit_bspline =
lm(BMI ~ 1 + ., data = nhanes_bspline_df)
bspline_coef_df =
tibble(
method = "B-Spline",
estimate =
tfd(t(B_bspline %*% coef(fit_bspline)[-1]), arg = epoch_arg))</code></pre>
<p>We show how to display all coefficient function estimates in the next
code chunk. The first step uses <code>bind_rows()</code> to combine data
frames containing the stepwise, quadratic, and B-spline estimated
coefficient functions. The result is a data with three rows, one for
each estimate, and two columns containing the <code>method</code> and
<code>estimate</code> variables. We plot the estimates using
<code>ggplot()</code> and <code>geom_spaghetti()</code> by setting the
aesthetics for <code>y</code> and <code>color</code> to
<code>estimate</code> and <code>method</code>, respectively. In the
resulting plot, the coefficient functions have some broad similarities
across basis specifications. That said, the quadratic basis has a much
higher estimate in the nighttime than other methods because of the
constraints on the shape of the coefficient function. The stepwise
coefficient has the bin average interpretation but the lack of
smoothness across bins is scientifically implausible. Of the
coefficients presented so far, then, the B-spline basis with eight
degrees of freedom is our preference as a way to include both
flexibility and smoothness in the estimate of <span
class="math inline">\(\beta_1(\cdot)\)</span>.</p>
<pre class="r"><code>bind_rows(stepfun_coef_df, quad_coef_df, bspline_coef_df) %>%
ggplot(aes(y = estimate, color = method)) +
geom_spaghetti(alpha = 1, linewidth = 1.2) +
scale_x_continuous(
breaks = seq(0, 24, length = 5),
labels = str_c(seq(0, 24, length = 5), ":00")) +
labs(x = "Time of day (hours)", y = "Coefficient")</code></pre>
<p><img src="chapter_04_files/figure-html/unnamed-chunk-9-1.png" width="90%" /></p>
</div>
<div id="penalized-spline-estimation" class="section level3">
<h3>Penalized spline estimation</h3>
<p>Our approach to scalar-on-function regression using smoothness
penalties relies on key insights that connect functional regression to
scatterplot smoothing and mixed models. By expressing the coefficient
function using a spline expansion it is possible to cast
scalar-on-function regression in terms of a linear regression model.
While the design matrix for that model has a specific construction, once
it is available, usual model fitting approaches can be used directly.
Analogously, the use of penalized splines requires the careful
construction of design and penalty matrices, but once these are in
place, the techniques for scatterplot smoothing. Critically, this
includes casting scalar-on-function regression as a mixed model; this
connection makes it possible to use the rich collection of techniques
for mixed model estimation and inference for functional regression.</p>
<p>Broadly speaking, in functional regression models we prefer to use
rich spline basis expansions because they are flexible and numerically
stable. We include penalties that encourage smoothness in the target of
estimation, which generally take the form of penalties on the overall
magnitude of the squared second derivative of the estimand. This
combination results in a quadratic penalty with a tuning parameter that
controls the degree of smoothness; higher penalization leads to
“flatter” coefficient functions, while lower penalization allows more
flexible but “wigglier” coefficients.</p>
<p>As before, let <span class="math inline">\(\beta_1(s) = \sum_{k =
1}^{K}\beta_{1k}B_{k}(s)\)</span> where <span
class="math inline">\(B_1(s), \dots, B_{K}(s)\)</span> is a collection
of basis functions. We add a smoothness-inducing penalty of the form
<span class="math inline">\(\int_S \{\beta_{1}''(s)\}^2
ds\)</span> to the minimization criterion. Intuitively, estimates of
<span class="math inline">\(\beta_{1}(s)\)</span> that include many
sharp turns will have squared second derivatives <span
class="math inline">\(\{\beta_{1}''(s)\}^2\)</span> with many
large values, while smooth estimates of <span
class="math inline">\(\beta_{1}(s)\)</span> will have second derivatives
that are close to 0 over the domain <span
class="math inline">\(S\)</span>. To implement the squared second
derivative penalty, let <span class="math inline">\(\mathbf{P}\)</span>
be the <span class="math inline">\(K \times K\)</span> matrix with the
<span class="math inline">\((i,j)^{\text{th}}\)</span> entry equal to
<span class="math inline">\(\int_S
B_i''(s)B_j''(s)ds\)</span>.</p>
<p>Let <span class="math inline">\(\mathbf{X}\)</span> be a <span
class="math inline">\(n \times (K+1)\)</span> matrix in which the <span
class="math inline">\(i^{\text{th}}\)</span> row is <span
class="math inline">\([1, C_{i1}, \dots,C_{iK}]\)</span>, and let <span
class="math inline">\(\boldsymbol{\beta}\)</span> be the <span
class="math inline">\((K+1)\)</span> dimensional column vector that
concatenates the population intercept <span
class="math inline">\(\beta_0\)</span> and the spline coefficients <span
class="math inline">\(\boldsymbol{\beta}_1\)</span>. Adding the second
derivative penalty yields a penalized sum of squares <span
class="math display">\[\begin{equation}\min_{\boldsymbol{\beta}}||\mathbf{y}-\mathbf{X}\boldsymbol{\beta}||^2+\lambda
\boldsymbol{\beta}^t\mathbf{D}\boldsymbol{\beta}\;.
\end{equation}\]</span> Here <span class="math inline">\(\lambda\geq
0\)</span> is a scalar tuning parameter and <span
class="math inline">\(\mathbf{D}\)</span> is the matrix given by <span
class="math display">\[\mathbf{D}=\begin{bmatrix}
\mathbf{0}_{1 \times 1} & \mathbf{0}_{K \times 1} \\
\;\mathbf{0}_{1 \times K} & \mathbf{P}\\
\end{bmatrix}\;,\]</span> where <span
class="math inline">\(\mathbf{0}_{a\times b}\)</span> is a matrix of
zero entries with <span class="math inline">\(a\)</span> rows and <span
class="math inline">\(b\)</span> columns. For fixed values of <span
class="math inline">\(\lambda\)</span>, a closed form solution for <span
class="math inline">\(\widehat{\boldsymbol{\beta}}\)</span> is <span
class="math inline">\(\widehat{\boldsymbol{\beta}}=(\mathbf{X}^t\mathbf{X}+\lambda
\mathbf{D})^{-1}\mathbf{X}^t\mathbf{y}\)</span>. Varying <span
class="math inline">\(\lambda\)</span> from 0 to <span
class="math inline">\(\infty\)</span> will induce no penalization and
full penalization, respectively, and choosing an appropriate tuning
parameter is an important practical challenge. As elsewhere in this
Chapter, though, we emphasize that the familiar form should not mask the
novelty and innovation of this model, which implements penalized spline
smoothing to estimate the coefficient function in a scalar-on-function
regression.</p>
<p>We illustrate these ideas in the next code chunk, which continues to
use BMI as an outcome and MIMS as a functional predictor. The code draws
on elements that have been seen previously. We first define a B-spline
basis with 30 degrees of freedom evaluated over the finite grid
<code>arg</code>. Using functionality in the <code>splines2</code>
package, we obtain the second derivative of each spline basis function
evaluated over the same finite grid. The elements of the penalty matrix
<span class="math inline">\(\mathbf{P}\)</span> are obtained through
numeric approximations to the integrals. The design matrix <span
class="math inline">\(\mathbf{X}\)</span> is obtained by adding a column
taking the value <span class="math inline">\(1\)</span> everywhere to
the terms <span class="math inline">\(C_{ik}\)</span>. Next, we
construct the matrix <span class="math inline">\(\mathbf{D}\)</span>.
The response vector <span class="math inline">\(\mathbf{y}\)</span> is
extracted from <code>nhanes_df</code> and we choose high and low values
for the tuning parameter <span class="math inline">\(\lambda\)</span>.
Given all of these elements, we estimate the coefficient vector <span
class="math inline">\(\boldsymbol{\beta}\)</span> to obtain
<code>coef_high</code> and <code>coef_low</code>.</p>
<pre class="r"><code>B_bspline = splines::bs(epoch_arg, df = 30, intercept = TRUE)
sec_deriv = splines2::bSpline(epoch_arg, df = 30, intercept = TRUE, derivs = 2)
P = t(sec_deriv) %*% sec_deriv * (1 / 60)
X = cbind(1, (nhanes_df$MIMS_mat %*% B_bspline) * (1 / 60))
D = rbind(0, cbind(0, P))
y = nhanes_df$BMI
lambda_high = 10e6
lambda_low = 100
coef_high = solve(t(X) %*% X + lambda_high * D) %*% t(X) %*% y
coef_low = solve(t(X) %*% X + lambda_low * D) %*% t(X) %*% y</code></pre>
<p>The estimated coefficient functions that correspond to the estimates
in <code>coef_high</code> and <code>coef_low</code> can be produced
through simple modifications to the previous code. A figure below shows
the resulting coefficient functions. Note that we specify the color for
these curves for consistency with later plots.</p>
<p><img src="chapter_04_files/figure-html/unnamed-chunk-11-1.png" width="90%" /></p>
<p>Recasting the penalized sum of squares as a mixed model allows the
data-driven estimation of tuning parameters; more broadly, this opens
the door to using a wide range of methods for mixed model estimation and
inference in functional regression settings. Below, we construct the
design matrix <span class="math inline">\(\mathbf{C}\)</span> using
numeric integration. The penalty matrix <span
class="math inline">\(\mathbf{P}\)</span>, which contains the numeric
integral of the squared second derivative of the basis functions, is
reused from a prior code chunk. We pass these as arguments into
<code>mgcv::gam()</code> by specifying <code>C</code> in the formula
that defines the regression structure, and then use the
<code>paraPen</code> argument to supply our penalty matrix
<code>P</code> for the design matrix <code>C</code>. Lastly, we specify
the estimation method to be REML.</p>
<pre class="r"><code>C = (nhanes_df$MIMS_mat %*% B_bspline) * (1 / 60)
fit_REML_penalty =
gam(y ~ 1 + C, paraPen = list(C = list(P)),
method = "REML")</code></pre>
<p>Code that multiplies the basis functions by the resulting spline
coefficients to obtain the estimated coefficient function is shown
below. We plot the result and include previous penalized estimates,
based on high and low values of the tuning parameter <span
class="math inline">\(\lambda\)</span>, which over- and under-smoothed
the coefficient function. The data-driven approach to tuning parameter
selection yields an estimate that is smooth but time-varying.</p>
<p><img src="chapter_04_files/figure-html/unnamed-chunk-13-1.png" width="90%" /></p>
<p>Many of the strengths of the <code>mgcv</code> package can be
leveraged through the <code>refund</code> package, which adds
functionality, quality of life features, and user interfaces relevant to
FDA. For SoFR, this means that instead of building models using
knowledge of the linear algebra underlying penalized spline estimation,