-
Notifications
You must be signed in to change notification settings - Fork 21
/
Copy path_main.Rmd
16190 lines (12510 loc) · 751 KB
/
_main.Rmd
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
---
title: "*Applied longitudinal data analysis* in brms and the tidyverse"
subtitle: "version 0.0.2"
author: "A Solomon Kurz"
date: "`r Sys.Date()`"
site: bookdown::bookdown_site
output:
bookdown::gitbook:
split_bib: yes
documentclass: book
bibliography: bib.bib
biblio-style: apalike
csl: apa.csl
link-citations: yes
geometry:
margin = 0.5in
urlcolor: blue
highlight: tango
header-includes:
\usepackage{underscore}
\usepackage[T1]{fontenc}
github-repo: ASKurz/Applied-Longitudinal-Data-Analysis-with-brms-and-the-tidyverse
twitter-handle: SolomonKurz
description: "This project is a reworking of Singer and Willett's classic (2003) text within a contemporary Bayesian framework with emphasis of the brms and tidyverse packages within the R computational framework."
---
# What and why {-}
This project is based on Singer and Willett's classic [-@singerAppliedLongitudinalData2003] text, [*Applied longitudinal data analysis: Modeling change and event occurrence*](https://www.oxfordscholarship.com/view/10.1093/acprof:oso/9780195152968.001.0001/acprof-9780195152968). You can download the data used in the text at [http://www.bristol.ac.uk/cmm/learning/support/singer-willett.html](https://www.bristol.ac.uk/cmm/learning/support/singer-willett.html) and find a wealth of ideas on how to fit the models in the text at [https://stats.idre.ucla.edu/other/examples/alda/](https://stats.idre.ucla.edu/other/examples/alda/). My contributions show how to fit these models and others like them within a Bayesian framework. I make extensive use of Paul Bürkner's [**brms** package](https://github.com/paul-buerkner/brms) [@R-brms; @burknerBrmsPackageBayesian2017; @burknerAdvancedBayesianMultilevel2018], which makes it easy to fit Bayesian regression models in **R** [@R-base] using Hamiltonian Monte Carlo (HMC) via the [Stan](https://mc-stan.org) probabilistic programming language [@carpenterStanProbabilisticProgramming2017]. Much of the data wrangling and plotting code is done with packages connected to the [**tidyverse**](http://style.tidyverse.org) [@R-tidyverse; @wickhamWelcomeTidyverse2019].
## Caution: Work in progress {-}
This release contains drafts of Chapters 1 through 6 and 9 through 13. Chapters 1 through 6 provide the motivation and foundational principles for fitting longitudinal multilevel models. Chapters 9 through 13 motivation and foundational principles for fitting discrete-time survival analyses.
In addition to fleshing out more of the chapters, I plan to add more goodies like introductions to multivariate longitudinal models and mixed-effect location and scale models. But there is no time-table for this project. To keep up with the latest changes, check in at the GitHub repository, [https://github.com/ASKurz/Applied-Longitudinal-Data-Analysis-with-brms-and-the-tidyverse](https://github.com/ASKurz/Applied-Longitudinal-Data-Analysis-with-brms-and-the-tidyverse), or follow my announcements on twitter at [https://twitter.com/SolomonKurz](https://twitter.com/SolomonKurz).
## **R** setup {-}
To get the full benefit from this ebook, you'll need some software. Happily, everything will be free (provided you have access to a decent personal computer and an good internet connection).
First, you'll need to install **R**, which you can learn about at [https://cran.r-project.org/](https://cran.r-project.org/).
Though not necessary, your **R** experience might be more enjoyable if done through the free RStudio interface, which you can learn about at [https://rstudio.com/products/rstudio/](https://rstudio.com/products/rstudio/).
Once you have installed **R**, execute the following to install the bulk of the add-on packages. This will probably take a few minutes to finish. Go make yourself a coffee.
```{r, eval = F}
packages <- c("bayesplot", "brms", "broom", "devtools", "flextable", "GGally", "ggmcmc", "ggrepel", "gtools", "loo", "patchwork", "psych", "Rcpp", "remotes", "rstan", "StanHeaders", "survival", "tidybayes", "tidyverse")
install.packages(packages, dependencies = T)
```
A few of the other packages are not officially available via the Comprehensive R Archive Network (CRAN; https://cran.r-project.org/). You can download them directly from GitHub by executing the following.
```{r, eval = F}
devtools::install_github("stan-dev/cmdstanr")
remotes::install_github("stan-dev/posterior")
devtools::install_github("rmcelreath/rethinking")
```
It's possible you'll have problems installing some of these packages. Here are some likely suspects and where you can find help:
* for difficulties installing **brms**, go to [https://github.com/paul-buerkner/brms#how-do-i-install-brms](https://github.com/paul-buerkner/brms#how-do-i-install-brms) or search around in the [**brms** section of the Stan forums ](https://discourse.mc-stan.org/c/interfaces/brms/36);
* for difficulties installing **cmdstanr**, go to [https://mc-stan.org/cmdstanr/articles/cmdstanr.html](https://mc-stan.org/cmdstanr/articles/cmdstanr.html);
* for difficulties installing **rethinking**, go to [https://github.com/rmcelreath/rethinking#quick-installation](https://github.com/rmcelreath/rethinking#quick-installation); and
* for difficulties installing **rstan**, go to [https://github.com/stan-dev/rstan/wiki/RStan-Getting-Started](https://github.com/stan-dev/rstan/wiki/RStan-Getting-Started).
## License and citation {-}
This book is licensed under the Creative Commons Zero v1.0 Universal license. You can learn the details, [here](https://github.com/ASKurz/Applied-Longitudinal-Data-Analysis-with-brms-and-the-tidyverse/blob/master/LICENSE). In short, you can use my work. Just please give me the appropriate credit the same way you would for any other scholarly resource. Here's the citation information:
```{r, eval = F}
@book{kurzAppliedLongitudinalDataAnalysis2021,
title = {Applied longitudinal data analysis in brms and the tidyverse},
author = {Kurz, A. Solomon},
year = {2021},
month = {4},
edition = {version 0.0.2},
url = {https://bookdown.org/content/4253/}
}
```
<!--chapter:end:index.Rmd-->
# A Framework for Investigating Change over Time
> It is possible to measure change, and to do it well, *if you have longitudinal data* [@rogosaGrowthCurveApproach1982; @willettResultsReliabilityLongitudinal1989]. Cross-sectional data--so easy to collect and so widely available--will not suffice. In this chapter, we describe why longitudinal data are necessary for studying change. [@singerAppliedLongitudinalData2003, p. 3, *emphasis* in the original]
## When might you study change over time?
Perhaps a better question is: *When wouldn't you?*
## Distinguishing between two types of questions about change
On page 8, Singer and Willett proposed there are two fundamental questions for longitudinal data analysis:
1. "How does the outcome change over time?" and
2. "Can we predict differences in these changes?"
Within the hierarchical framework, we often speak about two levels of change. We address within-individual change at *level-1*.
> The goal of a level-1 analysis is to describe the *shape* of each person's individual growth trajectory.
>
> In the second stage of an analysis of change, known as *level-2*, we ask about *interindividual differences in change*... The goal of a level-2 analysis is to detect heterogeneity in change across individuals and to determine the *relationship* between predictors and the *shape* of each person's individual growth trajectory. (p. 8, *emphasis* in the original)
## Three important features of a study of change
* Three or more waves of data
* An outcome whose values change systematically over time
* A sensible metric for clocking time
### Multiple waves of data.
Singer and Willett criticized two-waves data on two grounds.
> First, it cannot tell us about the *shape* of each person's individual growth trajectory, the focus of our level-1 question. Did all the change occur immediately after the first assessment? Was progress steady or delayed? Second, it cannot distinguish true change from measurement error. If measurement error renders pretest scores too low and posttest scores too high, you might conclude erroneously that scores increase over time when a longer temporal view would suggest the opposite. In statistical terms, two-waves studies cannot describe individual trajectories of change and they confound true change with measurement error [see @rogosaGrowthCurveApproach1982]. (p. 10, *emphasis* in the original)
I am not a fan of this 'true change/measurement error' way of speaking and would rather speak in terms of systemic and [seemingly] un-systemic changes among means and variances. Otherwise put, I'd rather speak in terms of trait and state. Two waves of data do not allow us to disentangle systemic mean differences from stable means and substantial variances for those means. Two waves of data do not allow us to disentangle changes in traits from stable traits but important differences in states. For an introduction to this way of thinking, check out Nezlek's [-@nezlekMultilevelFrameworkUnderstanding2007] [*A multilevel framework for understanding relationships among traits, states, situations and behaviors*](https://www.researchgate.net/publication/228079300_A_Multilevel_Framework_for_Understanding_Relationships_Among_Traits_States_Situations_and_Behaviours).
### A sensible metric for time.
> Choice of a time metric affects several interrelated decisions about the number and spacing of data collection waves....
>
> Our overarching point is that there is no single answer to the seemingly simple question about the most sensible metric for time. You should adopt whatever scale makes most sense for your outcomes and your research question....
>
> Our point is simple: choose a metric for the that reflect the cadence you expect to be most useful for your outcome. (p. 11)
Data collection waves can be evenly spaced or not. E.g., if you anticipate a time period of rapid nonlinear change, it might be helpful to increase the density of assessments during that period. Everyone does not have to have the same assessment schedule. If all are assessed on the same schedule, we describe the data as *time-structured*. When assessment schedules vary across participants, the data are termed *time-unstructured*. The data are *balanced* if all participants have the same number of waves. Issues like attrition and so on lead to *unbalanced* data. Though they may have some pedagogical use, I have not found these terms useful in practice.
### A continuous outcome that changes systematically over time.
To my eye, the most interesting part of this section is the discussion of measurement validity over time. E.g.,
> when we say the metric in which the outcome is measured must be preserved across time, we mean that the outcome scores must be equatable over time--a given value of the outcome on any occasion must represent the same "amount" of the outcome on every occasion. Outcome equatability is easiest to ensure when you use the identical instrument for measurement repeatedly over time. (p. 13)
This isn't as simple as is sounds. Though it's beyond the scope of this project, you might learn more about this from a study of the longitudinal measurement invariance literature. To dive in, see the first couple chapters in Newsom's [-@newsom2015longitudinal] text, [*Longitudinal structural equation modeling: A comprehensive introduction*](http://www.longitudinalsem.com/).
## Session info {-}
```{r}
sessionInfo()
```
<!--chapter:end:01.Rmd-->
```{r, echo = F, cache = F}
knitr::opts_chunk$set(fig.retina = 2.5)
knitr::opts_chunk$set(fig.align = "center")
options(width = 120)
```
# Exploring Longitudinal Data on Change
> Wise researchers conduct descriptive exploratory analyses of their data before fitting statistical models. As when working with cross-sectional data, exploratory analyses of longitudinal data con reveal general patterns, provide insight into functional form, and identify individuals whose data do not conform to the general pattern. The exploratory analyses presented in this chapter are based on numerical and graphical strategies already familiar from cross-sectional work. Owing to the nature of longitudinal data, however, they are inevitably more complex in this new setting. [@singerAppliedLongitudinalData2003, p. 16]
## Creating a longitudinal data set
> In longitudinal work, data-set organization is less straightforward because you can use two very different arrangements:
>
>>* *A person-level data set*, in which each person has one record and multiple variables contain the data from each measurement occasion
>
>>* *A person-period data set*, in which each person has multiple records—one for each measurement occasion (p. 17, *emphasis* in the original)
These are also sometimes referred to as the wide and long data formats, respectively.
As you will see, we will use two primary functions from the **tidyverse** to convert data from one format to another.
### The person-level data set.
Here we load the person-level data from [this UCLA web site](https://stats.idre.ucla.edu/r/examples/alda/r-applied-longitudinal-data-analysis-ch-2/). These are the NLY data [see @raudenbushGrowthCurveAnalysis2016] shown in the top of Figure 2.1.
```{r, warning = F, message = F}
library(tidyverse)
tolerance <- read_csv("https://stats.idre.ucla.edu/wp-content/uploads/2016/02/tolerance1.txt", col_names = T)
head(tolerance, n = 16)
```
With person-level data, each participant has a single row. In these data, participants are indexed by their `id` number. To see how many participants are in these data, just `count()` the rows.
```{r}
tolerance %>%
count()
```
The `nrow()` function will work, too.
```{r}
tolerance %>%
nrow()
```
With the base **R** `cor()` function, you can get the Pearson's correlation matrix shown in Table 2.1.
```{r}
cor(tolerance[ , 2:6]) %>%
round(digits = 2)
```
We used the `round()` function to limit the number of decimal places in the output. Leave it off and you'll see `cor()` returns up to seven decimal places instead. It can be hard to see the patters within a matrix of numerals. It might be easier in a plot.
```{r, fig.width = 3.75, fig.height = 2.25}
cor(tolerance[ , 2:6]) %>%
data.frame() %>%
rownames_to_column("row") %>%
pivot_longer(-row,
names_to = "column",
values_to = "correlation") %>%
mutate(row = factor(row) %>% fct_rev(.)) %>%
ggplot(aes(x = column, y = row)) +
geom_raster(aes(fill = correlation)) +
geom_text(aes(label = round(correlation, digits = 2)),
size = 3.5) +
scale_fill_gradient(low = "white", high = "red4", limits = c(0, 1)) +
scale_x_discrete(NULL, position = "top", expand = c(0, 0)) +
scale_y_discrete(NULL, expand = c(0, 0)) +
theme(axis.ticks = element_blank())
```
If all you wanted was the lower diagonal, you could use the `lowerCor()` function from the [**psych** package](https://CRAN.R-project.org/package=psych) [@R-psych].
```{r}
psych::lowerCor(tolerance[ , 2:6])
```
For more ways to compute, organize, and visualize correlations within the **tidyverse** paradigm, check out the [**corrr** package](https://tidymodels.github.io/corrr/) [@R-corrr].
### The person-period data set.
Here are the person-period data (i.e., those shown in the bottom of Figure 2.1).
```{r, warning = F, message = F}
tolerance_pp <- read_csv("https://stats.idre.ucla.edu/wp-content/uploads/2016/02/tolerance1_pp.txt",
col_names = T)
tolerance_pp %>%
slice(c(1:9, 76:80))
```
With data like these, the simple use of `count()` or `nrow()` won't help us discover how many participants there are in the `tolerance_pp` data. One quick way is to `count()` the number of `distinct()` `id` values.
```{r}
tolerance_pp %>%
distinct(id) %>%
count()
```
A fundamental skill is knowing how to convert longitudinal data in one format to the other. If you're using packages within the **tidyverse**, the `pivot_longer()` function will get you from the person-level format to the person-period format.
```{r}
tolerance %>%
# this is the main event
pivot_longer(-c(id, male, exposure),
names_to = "age",
values_to = "tolerance") %>%
# here we remove the `tol` prefix from the `age` values and then save the numbers as integers
mutate(age = str_remove(age, "tol") %>% as.integer()) %>%
# these last two lines just make the results look more like those in the last code chunk
arrange(id, age) %>%
slice(c(1:9, 76:80))
```
You can learn more about the `pivot_longer()` function [here](https://tidyr.tidyverse.org/reference/pivot_longer.html) and [here](https://tidyr.tidyverse.org/articles/pivot.html).
As hinted at in the above hyperlinks, the opposite of the `pivot_longer()` function is `pivot_wider()`. We can use `pivot_wider()` to convert the person-period `tolerance_pp` data to the same format as the person-level `tolerance` data.
```{r}
tolerance_pp %>%
# we'll want to add that `tol` prefix back to the `age` values
mutate(age = str_c("tol", age)) %>%
# this variable is just in the way. we'll drop it
select(-time) %>%
# here's the main action
pivot_wider(names_from = age, values_from = tolerance)
```
## Descriptive analysis of individual change over time
The following "descriptive analyses [are intended to] reveal the nature and idiosyncrasies of each person's temporal pattern of growth, addressing the question: How does each person change over time" (p. 23)?
### Empirical growth plots.
*Empirical growth plots* show individual-level sequence in a variable of interest over time. We'll put `age` on the $x$-axis, `tolerance` on the y-axis, and make our variant of Figure 2.2 with `geom_point()`. It's the `facet_wrap()` part of the code that splits the plot up by `id`.
```{r, fig.width = 4.5, fig.height = 5}
tolerance_pp %>%
ggplot(aes(x = age, y = tolerance)) +
geom_point() +
coord_cartesian(ylim = c(1, 4)) +
theme(panel.grid = element_blank()) +
facet_wrap(~id)
```
By default, **ggplot2** sets the scales of the $x$- and $y$-axes to the same values across subpanels. If you'd like to free that constraint, play around with the `scales` argument within `facet_wrap()`.
### Using a trajectory to summarize each person's empirical growth record.
If we wanted to connect the dots, we might just add a `geom_line()` line.
```{r, fig.width = 4.5, fig.height = 5}
tolerance_pp %>%
ggplot(aes(x = age, y = tolerance)) +
geom_point() +
geom_line() +
coord_cartesian(ylim = c(1, 4)) +
theme(panel.grid = element_blank()) +
facet_wrap(~id)
```
However, Singer and Willett recommend two other approaches:
* nonparametric smoothing
* parametric functions
#### Smoothing the empirical growth trajectory nonparametrically.
For our version of Figure 2.3, we'll use a loess smoother. When using the `stat_smooth()` function in **ggplot2**, you can control how smooth or wiggly the line is with the `span` argument.
```{r, fig.width = 4.5, fig.height = 5, message = F, warning = F}
tolerance_pp %>%
ggplot(aes(x = age, y = tolerance)) +
geom_point() +
stat_smooth(method = "loess", se = F, span = .9) +
coord_cartesian(ylim = c(1, 4)) +
theme(panel.grid = element_blank()) +
facet_wrap(~id)
```
#### Smoothing the empirical growth trajectory using ~~OLS~~ single-level Bayesian regression.
Although "fitting person-specific regression models, one individual at a time, is hardly the most efficient use of longitudinal data" (p. 28), we may as well play along with the text. It'll have pedagogical utility. You'll see.
For this section, we'll take a [cue from Hadley Wickham](https://www.youtube.com/watch?v=rz3_FDVt9eg&t=3458s) and use `group_by()` and `nest()` to make a tibble composed of tibbles (i.e., a nested tibble).
```{r}
by_id <-
tolerance_pp %>%
group_by(id) %>%
nest()
```
You can get a sense of what we did with `head()`.
```{r}
by_id %>% head()
```
As indexed by `id`, each participant now has their own data set stored in the `data` column. To get a better sense, we'll use our double-bracket subsetting skills to open up the first data set, the one for `id == 9`. If you're not familiar with this skill, you can learn more from [Chapter 9](https://bookdown.org/rdpeng/rprogdatascience/subsetting-r-objects.html) of [Roger Peng](https://twitter.com/rdpeng?lang=en)'s great [-@pengProgrammingDataScience2019] online book, [*R programming for data science*](https://bookdown.org/rdpeng/rprogdatascience/), or [Jenny Bryan](https://twitter.com/JennyBryan)'s fun and useful talk, [*Behind every great plot there's a great deal of wrangling*](https://www.youtube.com/watch?v=4MfUCX_KpdE).
```{r}
by_id$data[[1]]
```
Our `by_id` data object has many data sets stored in a higher-level data set. The code we used is verbose, but that's what made it human-readable. Now we have our nested tibble, we can make a function that will fit the simple linear model `tolerance ~ 1 + time` to each id-level data set. *Why use `time` as the predictor?* you ask. On page 29 in the text, Singer and Willett clarified they fit their individual models with $(\text{age} - 11)$ in order to have the model intercepts centered at 11 years old rather than 0. If we wanted to, we could make an $(\text{age} - 11)$ variable like so.
```{r}
by_id$data[[1]] %>%
mutate(age_minus_11 = age - 11)
```
Did you notice how our `age_minus_11` variable is the same as the `time` variable already in the data set? Yep, that's why we'll be using `time` in the model. In our data, $(\text{age} - 11)$ is encoded as `time`.
Singer and Willett used OLS to fit their exploratory models. We could do that to with the `lm()` function and we will do a little of that in this project. But let's get frisky and fit the models as Bayesians, instead. Our primary statistical package for fitting Bayesian models will be [Paul Bürkner](https://twitter.com/paulbuerkner?lang=en)'s [**brms**](https://github.com/paul-buerkner/brms). Let's open it up.
```{r, warning = F, message = F}
library(brms)
```
Since this is our first Bayesian model, we should start slow. The primary model-fitting function in **brms** is `brm()`. The function is astonishingly general and includes numerous arguments, most of which have sensible defaults. The primary two arguments are `data` and `formula`. I'm guessing they're self-explanatory. I'm not going to go into detail on the three arguments at the bottom of the code. We'll go over them later. For simple models like these, I would have omitted them entirely, but given the sparsity of the data (i.e., 5 data points per model), I wanted to make sure we gave the algorithm a good chance to arrive at reasonable estimates.
```{r fit2.1}
fit2.1 <-
brm(data = by_id$data[[1]],
formula = tolerance ~ 1 + time,
prior = prior(normal(0, 2), class = b),
iter = 4000, chains = 4, cores = 4,
seed = 2,
file = "fits/fit02.01")
```
We just fit a single-level Bayesian regression model for our first participant. We saved the results as an object named `fit2.1`. We can return a useful summary of `fit2.1` with either `print()` or `summary()`. Since it's less typing, we'll use `print()`.
```{r}
print(fit2.1)
```
The 'Intercept' and 'time' coefficients are the primary regression parameters. Also notice 'sigma', which is our variant of the residual standard error you might get from an OLS output (e.g., from base **R** `lm()`). Since we're Bayesians, the output summaries do not contain $p$-values. But we do get posterior standard deviations (i.e., the 'Est.Error' column) and the upper- and lower-levels of the percentile-based 95% intervals.
You probably heard somewhere that Bayesian statistics require priors. We can see what those were by pulling them out of our `fit2.1` object.
```{r}
fit2.1$prior
```
The prior in the top line, `normal(0, 2)`, is for all parameters of `class = b`. We actually specified this in our `brm()` code, above, with the code snip: `prior = prior(normal(0, 2), class = b)`. At this stage in the project, my initial impulse was to leave this line blank and save the discussion of how to set priors by hand for later. However, the difficulty is that the first several models we're fitting are all of $n = 5$. Bayesian statistics handle small-$n$ models just fine. However, when your $n$ gets small, the algorithms we use to implement our Bayesian models benefit from priors that are at least modestly informative. As it turns out, the **brms** default priors are flat for parameters of `class = b`. They offer no information beyond that contained in the likelihood. To stave off algorithm problems with our extremely-small-$n$ data subsets, we used `normal(0, 2)` instead. In our model, the only parameter of `class = b` is the regression slope for `time`. On the scale of the data, `normal(0, 2)` is a vary-permissive prior for our `time` slope.
In addition to our `time` slope parameter, our model contained an intercept and a residual variance. From the `fit2.1$prior` output, we can see those were `student_t(3, 2.1, 2.5)` and `student_t(3, 0, 2.5)`, respectively. **brms** default priors are designed to be weakly informative. Given the data and the model, these priors have a minimal influence on the results. We'll focus more on priors later in the project. For now just recognize that even if you don't specify your priors, you can't escape using some priors when using `brm()`. This is a good thing.
Okay, so that was the model for just one participant. We want to do that for all 16. Instead of repeating that code 15 times, we can work in bulk. With **brms**, you can reuse a model with the `update()` function. Here's how to do that with the data from our second participant.
```{r fit2.2}
fit2.2 <-
update(fit2.1,
newdata = by_id$data[[2]],
control = list(adapt_delta = .9),
file = "fits/fit02.02")
```
Peek at the results.
```{r}
print(fit2.2)
```
Different participants yield different model results.
Looking ahead a bit, we'll need to know how to get the $R^2$ for a single-level Gaussian model. With **brms**, you do that with the `bayes_R2()` function.
```{r}
bayes_R2(fit2.2)
```
Though the default spits out summary statistics, you can get the full posterior distribution for the $R^2$ by specifying `summary = F`.
```{r}
bayes_R2(fit2.2, summary = F) %>%
str()
```
Our code returned a numeric vector. If you'd like to plot the results with **ggplot2**, you'll need to convert the vector to a data frame.
```{r, fit.width = 4, fig.height = 2}
bayes_R2(fit2.2, summary = F) %>%
data.frame() %>%
ggplot(aes(x = R2)) +
geom_density(fill = "black") +
scale_x_continuous(expression(italic(R)[Bayesian]^2), limits = c(0, 1)) +
scale_y_continuous(NULL, breaks = NULL) +
theme(panel.grid = element_blank())
```
You'll note how non-Gaussian the Bayesian $R^2$ can be. Also, with the combination of default minimally-informative priors and only 5 data points, there' massive uncertainty in the shape. As such, the value of central tendency will vary widely based on which statistic you use.
```{r}
bayes_R2(fit2.2, summary = F) %>%
data.frame() %>%
summarise(mean = mean(R2),
median = median(R2),
mode = tidybayes::Mode(R2))
```
By default, `bayes_R2()` returns the mean. You can get the median with the `robust = TRUE` argument. To pull the mode, you'll need to use `summary = F` and feed the results into a mode function, like `tidybayes::Mode()`.
I should also point out the **brms** package did not get these $R^2$ values by traditional method used in, say, OLS estimation. To learn more about how the Bayesian $R^2$ sausage is made, check out Gelman, Goodrich, Gabry, and Vehtari's [-@gelmanRsquaredBayesianRegression2019] paper, [*R-squared for Bayesian regression models*]https://www.tandfonline.com/doi/full/10.1080/00031305.2018.1549100).
With a little tricky programming, we can use the `purrr::map()` function to serially fit this model to each of our participant-level data sets. We'll save the results as `fits`.
```{r, eval = F}
fits <-
by_id %>%
mutate(model = map(data, ~update(fit2.1, newdata = ., seed = 2)))
```
```{r fits, echo = F}
# save a little time
# save(fits, file = "fits/fits02.rda")
load(file = "fits/fits02.rda")
```
Let's walk through what we did. The `map()` function takes two primary arguments, `.x` and `.f`, respectively. We set `.x = data`, which meant we wanted to iterate over the contents in our `data` vector. Recall that each row of `data` itself contained an entire data set--one for each of the 16 participants. It's with the second argument `.f` that we indicated what we wanted to do with our rows of `data`. We set that to `.f = ~update(fit2.1, newdata = ., seed = 2)`. With the `~` syntax, we entered in a formula, which was `update(fit2.1, newdata = ., seed = 2)`. Just like we did with `fit2.2`, above, we reused the model formula and other technical specs from `fit2.1`. Now notice the middle part of the formula, `newdata = .`. That little `.` refers to the element we specified in the `.x` argument. What this combination means is that for each of the 16 rows of our nested `by_id` tibble, we plugged in the `id`-specific data set into `update(fit, newdata[[i]])` where `i` is simply meant as a row index. The new column, `model`, contains the output of each of the 16 iterations.
```{r}
print(fits)
```
Next, we'll want to extract the necessary summary information from our `fits` to remake our version of Table 2.2. There's a lot of info in that table, so let's take it step by step. First, we'll extract the posterior means (i.e., "Estimate") and standard deviations (i.e., "se") for the initial status and rate of change of each model. We'll also do the same for sigma (i.e., the square of the "Residual variance").
```{r, message = F}
mean_structure <-
fits %>%
mutate(coefs = map(model, ~ posterior_summary(.)[1:2, 1:2] %>%
data.frame() %>%
rownames_to_column("coefficients"))) %>%
unnest(coefs) %>%
select(-data, -model) %>%
unite(temp, Estimate, Est.Error) %>%
pivot_wider(names_from = coefficients,
values_from = temp) %>%
separate(b_Intercept, into = c("init_stat_est", "init_stat_sd"), sep = "_") %>%
separate(b_time, into = c("rate_change_est", "rate_change_sd"), sep = "_") %>%
mutate_if(is.character, ~ as.double(.) %>% round(digits = 2)) %>%
ungroup()
head(mean_structure)
```
It's simpler to extract the residual variance. Recall that because **brms** gives that in the standard deviation metric (i.e., $\sigma$), you need to square it to return it in a variance metric (i.e., $\sigma^2$).
```{r, message = F}
residual_variance <-
fits %>%
mutate(residual_variance = map_dbl(model, ~ posterior_summary(.)[3, 1])^2) %>%
mutate_if(is.double, round, digits = 2) %>%
select(id, residual_variance)
head(residual_variance)
```
We'll extract our Bayesian $R^2$ summaries, next. Given how nonnormal these are, we'll use the posterior median rather than the mean. We get that by using the `robust = T` argument within the `bayes_R2()` function.
```{r, message = F}
r2 <-
fits %>%
mutate(r2 = map_dbl(model, ~ bayes_R2(., robust = T)[1])) %>%
mutate_if(is.double, round, digits = 2) %>%
select(id, r2)
head(r2)
```
Here we combine all the components with a series of `left_join()` statements and present it in a [**flextable**](https://cran.r-project.org/web/packages/flextable/index.html)-type table.
```{r}
table <-
fits %>%
unnest(data) %>%
group_by(id) %>%
slice(1) %>%
select(id, male, exposure) %>%
left_join(mean_structure, by = "id") %>%
left_join(residual_variance, by = "id") %>%
left_join(r2, by = "id") %>%
rename(residual_var = residual_variance) %>%
select(id, init_stat_est:r2, everything()) %>%
ungroup()
table %>%
flextable::flextable()
```
We can make the four stem-and-leaf plots of Figure 2.4 with serial combinations of `pull()` and `stem()`.
```{r}
# fitted initial status
table %>%
pull(init_stat_est) %>%
stem(scale = 2)
# fitted rate of change
table %>%
pull(rate_change_est) %>%
stem(scale = 2)
# residual variance
table %>%
pull(residual_var) %>%
stem(scale = 2)
# r2 statistic
table %>%
pull(r2) %>%
stem(scale = 2)
```
To make Figure 2.5, we'll combine information from the original data and the 'Estimates' (i.e., posterior means) from our Bayesian models we've encoded in `mean_structure`.
```{r, fig.width = 4.5, fig.height = 5, message = F, warning = F}
by_id %>%
unnest(data) %>%
ggplot(aes(x = time, y = tolerance, group = id)) +
geom_point() +
geom_abline(data = mean_structure,
aes(intercept = init_stat_est,
slope = rate_change_est, group = id),
color = "blue") +
scale_x_continuous(breaks = 0:4, labels = 0:4 + 11) +
coord_cartesian(ylim = c(0, 4)) +
theme(panel.grid = element_blank()) +
facet_wrap(~id)
```
## Exploring differences in change across people
"Having summarized how each individual changes over time, we now examine similarities and differences in these changes across people" (p. 33).
### Examining the entire set of smooth trajectories.
The key to making our version of the left-hand side of Figure 2.6 is two `stat_smooth()` lines. The first one will produce the overall smooth. The second one, the one including the `aes(group = id)` argument, will give the `id`-specific smooths.
```{r, fig.width = 2.5, fig.height = 3.25, message = F, warning = F}
tolerance_pp %>%
ggplot(aes(x = age, y = tolerance)) +
stat_smooth(method = "loess", se = F, span = .9, size = 2) +
stat_smooth(aes(group = id),
method = "loess", se = F, span = .9, size = 1/4) +
coord_cartesian(ylim = c(0, 4)) +
theme(panel.grid = element_blank())
```
To get the linear OLS trajectories, just switch `method = "loess"` `to method = "lm"`.
```{r, fig.width = 2.5, fig.height = 3.25, message = F, warning = F}
tolerance_pp %>%
ggplot(aes(x = age, y = tolerance)) +
stat_smooth(method = "lm", se = F, span = .9, size = 2) +
stat_smooth(aes(group = id),
method = "lm", se = F, span = .9, size = 1/4) +
coord_cartesian(ylim = c(0, 4)) +
theme(panel.grid = element_blank())
```
But we wanted to be Bayesians. We already have the `id`-specific trajectories. All we need now is one based on all the data.
```{r fit2.3}
fit2.3 <-
update(fit2.1,
newdata = tolerance_pp,
file = "fits/fit02.03")
```
Here's the model summary.
```{r}
summary(fit2.3)
```
Before, we used `posterior_summary()` to isolate the posterior means and $SD$s. We can also use the `fixef()` function for that.
```{r}
fixef(fit2.3)
```
With a little subsetting, we can extract just the means from each.
```{r}
fixef(fit2.3)[1, 1]
fixef(fit2.3)[2, 1]
```
For this plot, we'll work more directly with the model formulas to plot the trajectories. We can use `init_stat_est` and `rate_change_est` from the `mean_structure` object as stand-ins for $\beta_{0i}$ and $\beta_{1i}$ from our model equation,
$$\text{tolerance}_{ij} = \beta_{0i} + \beta_{1i} \cdot \text{time}_{ij} + \epsilon_{ij},$$
where $i$ indexes children and $j$ indexes time points. All we need to do is plug in the appropriate values for `time` and we'll have the fitted `tolerance` values for each level of `id`. After a little wrangling, the data will be in good shape for plotting.
```{r}
tol_fitted <-
mean_structure %>%
mutate(`11` = init_stat_est + rate_change_est * 0,
`15` = init_stat_est + rate_change_est * 4) %>%
select(id, `11`, `15`) %>%
pivot_longer(-id,
names_to = "age",
values_to = "tolerance") %>%
mutate(age = as.integer(age))
head(tol_fitted)
```
We'll plot the `id`-level trajectories with those values and `geom_line()`. To get the overall trajectory, we'll get tricky with `fixef(fit2.3)` and `geom_abline()`.
```{r, fig.width = 2.5, fig.height = 3.25}
tol_fitted %>%
ggplot(aes(x = age, y = tolerance, group = id)) +
geom_line(color = "blue", size = 1/4) +
geom_abline(intercept = fixef(fit2.3)[1, 1] + fixef(fit2.3)[2, 1] * -11,
slope = fixef(fit2.3)[2, 1],
color = "blue", size = 2) +
coord_cartesian(ylim = c(0, 4)) +
theme(panel.grid = element_blank())
```
### Using the results of model fitting to frame questions about change.
If you're new to the multilevel model, the ideas in this section are foundational.
> To learn about the observed *average* pattern of change, we examine the sample averages of the fitted intercepts and slopes; these tell us about the average initial status and the average annual rate of change in the sample as a whole. To learn about the observed *individual differences* in change, we examine the sample *variances* and *standard deviations* of the intercepts and slopes; these tell us about the observed variability in initial status. And to learn about the observed relationship between initial status and the rate of change, we can examine the sample *covariance* or *correlation* between intercepts and slopes.
>
> Formal answers to these questions require the multilevel model for change of chapter 3. But we can presage this work by conducting simple descriptive analyses of the estimated intercepts and slopes. (p. 36, *emphasis* in the original)
Here are the means and standard deviations presented in Table 2.3.
```{r, message = F}
mean_structure %>%
pivot_longer(ends_with("est")) %>%
group_by(name) %>%
summarise(mean = mean(value),
sd = sd(value)) %>%
mutate_if(is.double, round, digits = 2)
```
Here's how to get the Pearson's correlation coefficient.
```{r}
mean_structure %>%
select(init_stat_est, rate_change_est) %>%
cor() %>%
round(digits = 2)
```
### Exploring the relationship between change and time-invariant predictors.
"Evaluating the impact of predictors helps you uncover systematic patterns in the individual change trajectories corresponding to interindividual variation in personal characteristics" (p. 37).
#### Graphically examining groups of smoothed individual growth trajectories.
If we'd like Bayesian estimates differing by `male`, we'll need to fit an interaction model.
```{r fit2.4}
fit2.4 <-
update(fit2.1,
newdata = tolerance_pp,
tolerance ~ 1 + time + male + time:male,
file = "fits/fit02.04")
```
Check the model summary.
```{r}
print(fit2.4)
```
Here's how to use `fixef()` and the model equation to get fitted values for `tolerance` based on specific values for `time` and `male`.
```{r}
tol_fitted_male <-
tibble(male = rep(0:1, each = 2),
age = rep(c(11, 15), times = 2)) %>%
mutate(time = age - 11) %>%
mutate(tolerance = fixef(fit2.4)[1, 1] +
fixef(fit2.4)[2, 1] * time +
fixef(fit2.4)[3, 1] * male +
fixef(fit2.4)[4, 1] * time * male)
tol_fitted_male
```
Now we're ready to make our Bayesian version of the top panels of Figure 2.7.
```{r, fig.width = 5, fig.height = 3.25}
tol_fitted %>%
# we need to add `male` values to `tol_fitted`
left_join(tolerance_pp %>% select(id, male),
by = "id") %>%
ggplot(aes(x = age, y = tolerance, color = factor(male))) +
geom_line(aes(group = id),
size = 1/4) +
geom_line(data = tol_fitted_male,
size = 2) +
scale_color_viridis_d(end = .75) +
coord_cartesian(ylim = c(0, 4)) +
theme(legend.position = "none",
panel.grid = element_blank()) +
facet_wrap(~male)
```
Before we can do the same thing with `exposure`, we'll need to dichotomize it by its median. A simple way is with a conditional statement within the `if_else()` function.
```{r}
tolerance_pp <-
tolerance_pp %>%
mutate(exposure_01 = if_else(exposure > median(exposure), 1, 0))
```
Now fit the second interaction model.
```{r fit2.5}
fit2.5 <-
update(fit2.4,
newdata = tolerance_pp,
tolerance ~ 1 + time + exposure_01 + time:exposure_01,
file = "fits/fit02.05")
```
Here's the summary.
```{r}
print(fit2.5)
```
Now use `fixef()` and the model equation to get fitted values for `tolerance` based on specific values for `time` and `exposure_01`.
```{r}
tol_fitted_exposure <-
crossing(exposure_01 = 0:1,
age = c(11, 15)) %>%
mutate(time = age - 11) %>%
mutate(tolerance = fixef(fit2.5)[1, 1] +
fixef(fit2.5)[2, 1] * time +
fixef(fit2.5)[3, 1] * exposure_01 +
fixef(fit2.5)[4, 1] * time * exposure_01,
exposure = if_else(exposure_01 == 1, "high exposure", "low exposure") %>%
factor(., levels = c("low exposure", "high exposure")))
tol_fitted_exposure
```
Did you notice in the last lines in the second `mutate()` how we made a version of `exposure` that is a factor? That will come in handy for labeling and ordering the subplots. Now make our Bayesian version of the bottom panels of Figure 2.7.
```{r, fig.width = 5, fig.height = 3.25}
tol_fitted %>%
# we need to add `exposure_01` values to `tol_fitted`
left_join(tolerance_pp %>% select(id, exposure_01),
by = "id") %>%
mutate(exposure = if_else(exposure_01 == 1, "high exposure", "low exposure") %>%
factor(., levels = c("low exposure", "high exposure"))) %>%
ggplot(aes(x = age, y = tolerance, color = exposure)) +
geom_line(aes(group = id),
size = 1/4) +
geom_line(data = tol_fitted_exposure,
size = 2) +
scale_color_viridis_d(option = "A", end = .75) +
coord_cartesian(ylim = c(0, 4)) +
theme(legend.position = "none",
panel.grid = element_blank()) +
facet_wrap(~exposure)
```
#### The relationship between ~~OLS-Estimated~~ single-level Bayesian trajectories and substantive predictors
"To investigate whether fitted trajectories vary systematically with predictors, we can treat the estimated intercepts and slopes as outcomes and explore the relationship between them and predictors" (p. 39). Here are the left panels of Figure 2.8.
```{r, fig.width = 2.5, fig.height = 5}
p1 <-
mean_structure %>%
pivot_longer(ends_with("est")) %>%
mutate(name = factor(name, labels = c("Fitted inital status", "Fitted rate of change"))) %>%
# we need to add `male` values to `tol_fitted`
left_join(tolerance_pp %>% select(id, male),
by = "id") %>%
ggplot(aes(x = factor(male), y = value, color = name)) +
geom_point(alpha = 1/2) +
scale_color_viridis_d(option = "B", begin = .2, end = .7) +
labs(x = "male",
y = NULL) +
theme(legend.position = "none",
panel.grid = element_blank()) +
facet_wrap(~name, scale = "free_y", ncol = 1)
p1
```
Here are the right panels.
```{r, fig.width = 2.5, fig.height = 5}
p2 <-
mean_structure %>%
pivot_longer(ends_with("est")) %>%
mutate(name = factor(name, labels = c("Fitted inital status", "Fitted rate of change"))) %>%
# we need to add `male` values to `tol_fitted`
left_join(tolerance_pp %>% select(id, exposure),
by = "id") %>%
ggplot(aes(x = exposure, y = value, color = name)) +
geom_point(alpha = 1/2) +
scale_color_viridis_d(option = "B", begin = .2, end = .7) +
scale_x_continuous(breaks = 0:2,
limits = c(0, 2.4)) +
labs(y = NULL) +
theme(legend.position = "none",
panel.grid = element_blank()) +
facet_wrap(~name, scale = "free_y", ncol = 1)
p2
```
Did you notice how we saved those last two plots as `p1` and `p2`? We can use syntax from the [**patchwork** package](https://patchwork.data-imaginist.com/) [@R-patchwork] to combine them into one compound plot.
```{r, fig.width = 5, fig.height = 5}
library(patchwork)
p1 + p2 + scale_y_continuous(breaks = NULL)
```
As interesting as these plots are, do remember that "the need for ad hoc correlations has been effectively replaced by the widespread availability of computer software for fitting the multilevel model for change directly" (pp. 41--42). As you'll see, Bürkner's **brms** package is one of the foremost in that regard.
## Improving the precision and reliability of ~~OLS~~ single-level-Bayesian-estimated rates of change: Lessons for research design
> Statisticians assess the precision of a parameter estimate in terms of its *sampling variation*, a measure of the variability that would be found across infinite resamplings from the same population. The most common measure of sampling variability is an estimate's *standard error*, the square root of its estimated sampling variance. Precision and standard error have an inverse relationship; the smaller the standard error, the more precise the estimate. (p. 41, *emphasis* in the original)
So here's the deal: When Singer and Willett wrote "Statisticians assess..." a more complete expression would have been 'Frequentist statisticians assess...' Bayesian statistics are not based on asymptotic theory. They do not presume an idealized infinite distribution of replications. Rather, Bayesian statistics use Bayes theorem to estimate the probability of the parameters given the data. That probability has a distribution. Analogous to frequentist statistics, we often summarize that distribution (i.e., the posterior distribution) in terms of central tendency (e.g., posterior mean, posterior median, posterior mode) and spread. *Spread?* you say. We typically express spread in one or both of two ways. One typical expression of spread is the 95% intervals. In the Bayesian world, these are often called credible or probability intervals. The other typical expression of spread is the *posterior standard deviation*. In **brms**, this of typically summarized in the 'Est.error' column of the output of functions like `print()` and `posterior_summary()` and so on. The posterior standard deviation is analogous to the frequentist standard error. Philosophically and mechanically, they are *not* the same. But in practice, they are often quite similar.
Later we read:
> Unlike precision which describes how well an individual slope estimate measures that person's true rate of change, reliability describes how much the rate of change varies across people. Precision has meaning for the individual; reliability has meaning for the group. (p. 42)
I have to protest. True, if we were working within a Classical Test Theory paradigm, this would be correct. But this places reliability with the context of group-based cross-sectional design. Though this is a popular design, it is not the whole story (i.e., see this book!). For introductions to more expansive and person-specific notions of reliability, check out [Lee Cronbach](https://en.wikipedia.org/wiki/Lee_Cronbach)'s Generalizability Theory [@cronbachDependabilityBehavioralMeasurements1972; @brennanGeneralizabilityTheory2001; also @cranfordProcedureEvaluatingSensitivity2006; @lopilatoUpdatingGeneralizabilityTheory2015; @shroutPsychometrics2012].
## Session info {-}
```{r}
sessionInfo()
```
```{r, echo = F, message = F}
# here we'll remove our objects
rm(tolerance, tolerance_pp, by_id, fit2.1, fit2.2, fits, mean_structure, residual_variance, r2, table, fit2.3, tol_fitted, fit2.4, tol_fitted_male, fit2.5, tol_fitted_exposure, p1, p2)
theme_set(theme_grey())
pacman::p_unload(pacman::p_loaded(), character.only = TRUE)
```
<!--chapter:end:02.Rmd-->
```{r, echo = F, cache = F}
knitr::opts_chunk$set(fig.retina = 2.5)
knitr::opts_chunk$set(fig.align = "center")
options(width = 110)
```
# Introducing the Multilevel Model for Change
> In this chapter [Singer and Willett introduced] the multilevel model for change, demonstrating how it allows us to address within-person and between-person questions about change simultaneously. Although there are several ways of writing the statistical model, here we adopt a simple and common approach that has much substantive appeal. We specify the multilevel model for change by simultaneously postulating a pair of subsidiary models—a level-1 submodel that describes how each person changes over time, and a level-2 model that describes how these changes differ across people [@bryk1987application; @rogosaUnderstandingCorrelatesChange1985]. [@singerAppliedLongitudinalData2003, p. 3]
## What is the purpose of the multilevel model for change?
Unfortunately, we do not have access to the full data set Singer and Willett used in this chapter. For details, go [here](https://stats.idre.ucla.edu/r/examples/alda/r-applied-longitudinal-data-analysis-ch-3/). However, I was able to use the data provided in Table 3.1 and the model results in Table 3.3 to simulate data with similar characteristics as the original. To see how I did it, look at the section at the end of the chapter.
Anyway, here are the data in Table 3.1.
```{r, warning = F, message = F}
library(tidyverse)
early_int <-
tibble(id = rep(c(68, 70:72, 902, 904, 906, 908), each = 3),
age = rep(c(1, 1.5, 2), times = 8),
cog = c(103, 119, 96, 106, 107, 96, 112, 86, 73, 100, 93, 87,
119, 93, 99, 112, 98, 79, 89, 66, 81, 117, 90, 76),
program = rep(1:0, each = 12))
print(early_int)
```
Later on, we also fit models using $age - 1$. Here we'll compute that and save it as `age_c`.
```{r}
early_int <-
early_int %>%
mutate(age_c = age - 1)
head(early_int)
```
Here we'll load our simulation of the full $n = 103$ data set.
```{r}
load("data/early_int_sim.rda")
```
## The level-1 submodel for individual change
This part of the model is also called the *individual growth model*. Remember how in last chapter we fit a series of participant-specific models? That's the essence of this part of the model.
Here's our version of Figure 3.1. Note that here we're being lazy and just using OLS estimates.
```{r, fig.width = 6.5, fig.height = 4, message = F, warning = F}
early_int %>%
ggplot(aes(x = age, y = cog)) +
stat_smooth(method = "lm", se = F) +
geom_point() +
scale_x_continuous(breaks = c(1, 1.5, 2)) +
ylim(50, 150) +
theme(panel.grid = element_blank()) +
facet_wrap(~id, ncol = 4)
```
Based on these data, we postulate our level-1 submodel to be
$$
\text{cog}_{ij} = [ \pi_{0i} + \pi_{1i} (\text{age}_{ij} - 1) ] + [\epsilon_{ij}].
$$
### The structural part of the level-1 submodel.
As far as I can tell, the data for Figure 3.2 are something like this.
```{r}
d <-
tibble(id = "i",
age = c(1, 1.5, 2),
cog = c(95, 100, 135))
d
```
To add in the horizontal dashed lines in Figure 3.2, we'll need to fit a model. Let's be lazy and use OLS. Don't worry, we'll use Bayes in a bit.
```{r}