forked from winterwang/LSHTMlearningnote
-
Notifications
You must be signed in to change notification settings - Fork 0
/
10-Hierarchical-models.Rmd
4315 lines (3111 loc) · 223 KB
/
10-Hierarchical-models.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
# (PART) 等級線性迴歸模型 analysis of hierarchical and other dependent data {-}
# 相互依賴數據及簡單的應對方案 {#Hierarchical}
> To call in the statistician after the experiment is done may be no more than asking him to perform a post-mortem examination: he may be able to say what the experiment died of.
>
> ~ Sir Ronald Aylmer Fisher
```{block2, note-thankslinda, type='rmdnote'}
The Analysis of Hierarchical and Other Dependent Data lectures were orgainised and taught by Professor [Linda Sharples](https://www.lshtm.ac.uk/aboutus/people/sharples.linda), and Dr. [Edmund Njeru Njagi](https://www.lshtm.ac.uk/aboutus/people/njagi.edmund-njeru).
```
## 相互依賴的數據
線性回歸模型，廣義線性回歸模型，他們背後都有一個十分十分**十分重要的假設--數據的相互獨立性**。這個前提假設常常會在現實數據中得不到滿足，因爲數據與數據之間在背後很可能會有有所關聯，也許是已知的，也許是未知的因素讓某些數據顯得更加接近彼此。這個章節，主要的內容就是舉例說明分層數據在日常生活中的常見性，以及處理這個非獨立性質的必要性。
- 圖 \@ref(fig:Hier01-1) 展示的箱式圖顯示的是六個不同醫院對各自 12 名患者收縮期血壓測量的結果。如果把醫院看做一個單位，取院內患者的平均值，那麼六所醫院的血壓均值最大爲 135.7 mmHg，最小是 117.7 mmHg，六所醫院測量的血壓總體均值爲 125.6 mmHg。
```{r Hier01-1, echo=FALSE, fig.height=6, fig.width=7, fig.cap='Box and whiskers plot of measured SBP in patients from six hospitals', fig.align='center', out.width='90%', cache=TRUE}
Bp <- read_dta("backupfiles/bp.dta")
Bp$hosp <- as.factor(Bp$hosp)
with(Bp, boxplot(bp ~ hosp, xlab = "Hospital No.", ylab = "Systolic blood pressure (mmHg)"))
```
- 圖 \@ref(fig:Hier01-2) 展示的是對 17 名患者使用兩種不同的測量方法測量的最大呼吸速率 (peak-expiratory-flow rate, PEFR)。兩種方法又測量了兩次，途中展示的是其中一種測量方法前後兩次測量結果的散點圖。
```{r Hier01-2, cache=TRUE, echo=FALSE, fig.height=6, fig.width=9, fig.cap='Two recordings of PEFR taken with the Mini Wright meter', fig.align='center', out.width='80%', message=FALSE, warning=FALSE}
pefr <- read_dta("backupfiles/pefr.dta")
# the data are in wide format
# transform data into long format
pefr_long <- pefr %>%
gather(key, value, -id) %>%
separate(key, into = c("measurement", "occasion"), sep = 2) %>%
arrange(id, occasion) %>%
spread(measurement, value)
## figure shows slightly closer agreement between the repeated measures of standard Wright,
## than between those of Mini Wright
ggplot(pefr_long, aes(x = id, y = wm, fill = occasion)) +
geom_point(size = 4, shape = 21) +
geom_hline(yintercept = mean(pefr_long$wm), colour = "red") +
theme_bw() +
scale_x_continuous(breaks = 1:17)+
theme(axis.text = element_text(size = 15),
axis.text.x = element_text(size = 15),
axis.text.y = element_text(size = 15)) +
labs(x = "Subject ID", y = "MW Measurements") +
theme(axis.title = element_text(size = 17), axis.text = element_text(size = 8),
axis.line = element_line(colour = "black"),
panel.border = element_blank(),
panel.background = element_blank()) + theme(legend.position = "bottom", legend.direction = "horizontal") + theme(legend.text = element_text(size = 19),
legend.title = element_text(size = 19))
```
- 圖 \@ref(fig:Hier01-3) 展示的來自全英 65 所學校的 4059 名學生入學前閱讀水平測試成績 (LRT) 和畢業時 GCSE 考試成績之間的散點圖關系。值得注意的是該圖其實無視了學校這個變量，把每個學生看成相互獨立的個體。但是當我們隨機選取四所學校，看它們各自的學生的成績表現 (圖 \@ref(fig:Hier01-4))。很顯然，之前忽視了學校這一層級的變量是不恰當的，因爲不同學校學生的入學前和畢業時成績之間的相關性很明顯存在不同的模式 (四所學校的回歸線各自的截距和斜率各不相同)。
```{r Hier01-3, cache=TRUE, echo=FALSE, fig.height=6, fig.width=9, fig.cap='GCSE by LRT in all 65 schools', fig.align='center', out.width='80%', message=FALSE, warning=FALSE}
gcse_selected <- read_dta("backupfiles/gcse_selected.dta")
ggthemr::ggthemr('fresh')
ggplot(gcse_selected, aes(x = lrt, y = gcse)) + geom_point(size = 2.5) +
geom_smooth(method = "lm", se = FALSE) +
xlim(-40, 40) +
theme(axis.text = element_text(size = 15),
axis.text.x = element_text(size = 15),
axis.text.y = element_text(size = 15)) +
labs(x = "LRT", y = "GCSE score")
```
```{r Hier01-4, cache=TRUE, echo=FALSE, fig.height=6, fig.width=9, fig.cap='GCSE by LRT in four randomly selected schools', fig.align='center', out.width='80%', message=FALSE, warning=FALSE}
ggthemr('fresh')
p <- ggplot(gcse_selected[gcse_selected$school %in% c(2, 7, 40, 53), ], aes(x = lrt, gcse)) +
geom_point(size = 2.5) +
geom_smooth(method = "lm", se = FALSE) +
theme(axis.text = element_text(size = 15),
axis.text.x = element_text(size = 15),
axis.text.y = element_text(size = 15)) +
labs(x = "LRT", y = "GCSE score")
p + facet_wrap( ~ school, ncol=2)+
theme(strip.text = element_text(face = "bold", size = rel(1.5)))
```
- 另一個特別好的例子展示在圖 \@ref(fig:Hier01-5) 中，是關於同一個母親的不同孩子的出生體重的數據。一個母親可以有多個孩子，每個母親的孩子之間的出生體重很明顯無法看作相互獨立。圖中展示的是，3300 名生了兩個孩子的母親的孩子們出生體重的散點圖。同一個母親的小孩用線相連。顯然，同一個母親生的孩子，其出生體重比不同母親的孩子出生體重差距更小，更接近彼此，因爲他們來自同一個母親。可以想象，一個母親如果身材高大，那麼她的孩子們可能都傾向於有比較高的出生體重。所以同一個母親的孩子之間體重是有相關關系的 (within correlation)。
```{r Hier01-5, cache=TRUE, echo=FALSE, fig.height=6, fig.width=9, fig.cap='Birthweight of siblings by maternal identifier', fig.align='center', out.width='80%', message=FALSE, warning=FALSE}
siblings <- read_dta("backupfiles/siblings.dta")
subSiblings <- siblings[(siblings$momid > 80000)&(siblings$momid < 95000)&(siblings$idx %in% c(1,2)), ]
subSiblings$idx <- factor(subSiblings$idx)
subSiblings_w <- subset(subSiblings, select = c("momid", "idx", "birwt"))
subSiblings_w <- spread(subSiblings_w, key = idx, value = birwt)
ggthemr('fresh')
ggplot(subSiblings, aes(x = momid, y = birwt, fill = idx)) +
geom_point(size = 4, shape = 21) +
#geom_segment(aes(x = momid, y = `1`, xend = momid, yend = `2`, colour = "segment"), data = subSiblings_w) +
geom_line(aes(group = momid), lty = 1) +
theme(axis.text = element_text(size = 15),
axis.text.x = element_text(size = 15),
axis.text.y = element_text(size = 15)) +
labs(x = "Mother ID", y = "Birthweight (g)") +
theme(axis.title = element_text(size = 17), axis.text = element_text(size = 8),
axis.line = element_line(colour = "black"),
panel.border = element_blank(),
panel.background = element_blank()) +
theme(legend.position = "bottom", legend.direction = "horizontal") + theme(legend.text = element_text(size = 19), legend.title = element_text(size = 19)) + theme(plot.subtitle = element_text(vjust = 1),
plot.caption = element_text(vjust = 1)) +labs(fill = "Child number") +
scale_fill_discrete(labels = c("1st Child", "2nd Child"))
```
- 最後一個用於本章節的實例是，一項研究亞洲兒童生長狀況的調查分別記錄了 198 個數據點，68 個兒童在 0 到 3 歲之間的四個年齡點的體重數據。圖 \@ref(fig:Hier01-6) 展示的就是這個典型的隨訪數據的個人生長曲線。且圖中每個人的生長軌跡提示，男孩子的生長過程可能相互之間體重差異顯得較女孩子來得大。如果，我們用每個兒童自己的數據，給每個兒童擬合各自的回歸線，數據顯然不足，但是如果我們決定忽略個體的生長的隨機效應 (不均一性)，又顯得十分不妥當。
```{r Hier01-6, cache=TRUE, echo=FALSE, fig.height=6, fig.width=9, fig.cap='Growth profiles of boys and girls in the Asian growth data', fig.align='center', out.width='80%', message=FALSE, warning=FALSE}
growth <- read_dta("backupfiles/asian.dta")
growth <- growth %>%
mutate(gender = factor(gender, labels= c("Boys", "Girls")))
ggthemr('fresh')
G <- ggplot(growth, aes(x = age, y = weight)) +
geom_line(aes(group = id), lty = 1) +
theme(axis.text = element_text(size = 15),
axis.text.x = element_text(size = 15),
axis.text.y = element_text(size = 15)) +
labs(x = "Age (years)", y = "Weight (Kg)") +
theme(axis.title = element_text(size = 17), axis.text = element_text(size = 8),
axis.line = element_line(colour = "black"),
panel.border = element_blank(),
panel.background = element_blank()) +
theme(legend.position = "bottom", legend.direction = "horizontal") + theme(legend.text = element_text(size = 19), legend.title = element_text(size = 19)) + theme(plot.subtitle = element_text(vjust = 1),
plot.caption = element_text(vjust = 1)) +labs(fill = "Child number") +
scale_fill_discrete(labels = c("1st Child", "2nd Child"))
G + facet_grid(. ~ gender) +
theme(strip.text = element_text(face = "bold", size = rel(1.5)))
```
## 依賴性的來源在哪裏
上述例子中的數據，均提示我們數據與數據之間獨立性的假設，常常會遇到尷尬的局面。因爲數據與數據之間本身就不可能完全獨立。
1. 同一個診所或者醫院的患者，他們之間可能有着某些相似的因素從而導致他們的血壓相比其他醫院的人彼此更加接近。這個原因可能是有同一家醫院的患者可能有類似的疾病。
2. 同一患者身上反復抽取樣本，也就是說一個對象貢獻了多個數據的時候，這些來自同一對象的數據當然具有相對不同對象數據更高的均質性。
3. 同一所學校的學生的成績或內部的相關性，很可能大於不同學校兩個學生之間成績的相關性。因爲同一學校的孩子可能共享某些共同的特徵，比如說相似的家庭經濟背景，或者是同樣的教學內容教學老師等環境因素。這樣，來自同一所學校的孩子的成績很可能就會更加相似。
4. 至於說家庭數據就更加典型了。來自同一家庭的兄弟姐妹，有着極強的相關性，因爲他們共享着遺傳因素，或者是相似的家庭教育/飲食/生活習慣等環境因素。
5. 同一個體身上的縱向 (時間) 隨訪數據很顯然會比不同患者有更強的內部相關性。
目前位置介紹的這些常見實例中，可以發現它們有一個共通點。就是這些數據其實內部是有分層結構的 (hierarchy)。這些數據中，都有一個最底層單元 (elementary units/level 1)，還有一個聚合單元 (aggregate units/level 2)，聚合單元常被命名爲層級 (cluster)。
```{r Hier01tab00, echo=FALSE, cache=TRUE, eval=TRUE}
dt <- read.csv("backupfiles/Hier01tab00.csv", header = T)
#names(dt) <- c("Model with", "sigma_u", "sigma_e", "sigma_u", "sigma_e")
kable(dt, "html", align = "l", caption = "Hierarchy in the data (5 examples in Chapter 1)") %>%
kable_styling(bootstrap_options = c("striped", "bordered"),full_width = F, position = "center") %>%
add_header_above(c("Level" = 2))
```
正如表格 \@ref(tab:Hier01tab00) 總結的那樣，這些數據中存在這層級結構，這種數據被稱爲分層數據 (hierarchical)，或者叫嵌套式數據 (nested data)。根據你所在的知識領域，它可能還被命名爲多層結構數據 (multilevel and clustered data)。在一些研究中，你可能會遇見從實驗設計階段就存在分層結構的數據，比如使用分層抽樣 (multistage sampling) 的設計的實驗等。這樣的實驗設計最常在人口學，社會學的研究中看到。在大多數醫學研究中，每個數據點 (observation point, level 1)，所屬的層 (cluster) 本身可能是我們感興趣的研究點 (例如同屬於一個家庭，相同母親的後代)，又或者是同一個人/患者的隨着時間推移的隨訪健康狀態 (如生長曲線，體重變化，疾病康復情況)。
如果用前面用過的 圖 \@ref(fig:Hier01-6) 的生長曲線做例子，那麼每個被調查的兒童，就是該數據的第二級層，每個隨訪時刻測量的體重數據，則是觀察的數據點。這個數據還有一個特點是，觀察數據點是有前後的 (時間) 順序的，這是一個典型的**縱向研究數據 (longitudinal data)**。
## 數據有依賴性導致的結果
如果你手頭的數據，結構上是一種嵌套式結構數據，那麼任何無視了這一點作出的統計學推斷都是有瑕疵的。相互之間互不獨立這一特質，需要通過一種新的手段，把嵌套式的數據結構考慮進統計學模型裏來。
在一些情況下，數據的嵌套式結構可能可以被忽略掉，但是其結果是導致統計學的估計變得十分低效 (inefficient procedure)。你可能會聽說過一般化估計公式 (generalized estimating equations)，是其中一種備擇手段，因爲在這一公式中，你需要人爲地指定數據與數據之間可能的依賴關系是怎樣的。
其實，即使有人真的在分析過程中忽略了數據本身的嵌套式結構，他會發現最終在描述分析結果的時候，還是無法避免這一嚴重的問題。另外一些統計學家可能記得在穩健統計學法中，三明治標準誤估計法也是可以供選擇的一種處理相關數據的手段。
## 邊際模型和條件模型 marginal and conditional models
邊際模型和條件模型的概念其實不是分層模型特有的，卻在分析分層數據模型時十分有用。假如，對於某個結果變量 $Y$ 有它如下的回歸模型，其中我們把某個單一的共變量 $Z$ 從模型中分離出來，加以特別關注:
$$
g\{ \text{E}(Y|\textbf{X},Z) \} = \beta\textbf{X} +\gamma Z
$$
這是一個典型的條件模型，它描述了結果變量 $Y$ 的期望是以怎樣的**條件**和解釋變量 $\textbf{X},Z$ 之間建立關系的。每個解釋變量的回歸系數，其含義都是**以其他同一模型中的共變量不變的條件下**，和結果變量之間的關系。經過這樣的解讀，你會知道，其實本統計學教程目前爲止遇見過的所有的回歸模型都是條件模型。如果此時我們反過來思考，把上述模型中單獨分離出來的單一共變量 $Z$ 對於結果變量 $Y$ 均值的影響合並起來 (對共變量 $Z$ 積分即可)，此時我們得到的就是共變量 $\textbf{X}$ 和結果變量 $Y$ 之間，關於 $Z$ 的邊際模型 (Marginal model):
$$
\text{E}_Z\{ \text{E}(Y|\textbf{X}, Z) \} = \text{E}_Z\{ g^{-1}(\beta\textbf{X} + \gamma Z) \} \\
\text{Where } \text{E}(Z) = 0
$$
用**線性回歸**來舉例:
$$
\text{E}(Y| \textbf{X}, Z) = \beta\textbf{X} + \gamma Z
$$
那麼此時共變量 $\textbf{X}$ 的邊際模型回歸系數 $\beta$ 的含義，和條件模型時的回歸系數其實是相同的含義:
$$
\text{E}_Z\{\text{E}(Y|\textbf{X},Z)\} = \text{E}_Z(\beta\textbf{X} + \gamma Z) = \beta\textbf{X} + \gamma\text{E}(Z) = \beta\textbf{X}
$$
爲什麼這裏的邊際模型對於分層數據來說很重要呢？答案在於，嵌套式數據中，我們常常關心那第二個階層 (重復測量某個指標的患者，學生成績數據中的學校層級，等) 在它所在的那個階層中和結果變量之間的平均關系。(In models for hierarchical data we often use level effects to represent what is common among observations from one "cluster" or "group". We may then want marginal conclusions: we need to average over these effects).
```{r hier01tab01, echo=FALSE, cache=TRUE, eval=FALSE}
dt <- read.csv("backupfiles/hierexample.csv", header = T)
names(dt) <- c("Cluster (j)", "id (i)", "X", "Y", "X", "Y")
kable(dt, "html", align = "c", caption = "Example data") %>%
kable_styling(bootstrap_options = c("striped", "bordered"),full_width = F, position = "center") #%>%
#add_header_above(c(" " = 1,"REML" = 2, "ML" = 2))
```
### 標記法 notation
- $Y_{ij}$ 標記第 $j$ 層的第 $i$ 個個體;
- $i = 1, \cdots, n_j$ 表示第 $j$ 層中共有 $n_j$ 個個體 (elements);
- $j = 1, \cdots, J$ 表示數據共有 $J$ 個第二階層 (clusters);
- $N = \sum_{j=1}^J n_j$ 表示總體樣本量等於各個階層樣本量之和;
- 特殊情況: 如果每個階層的個體數相同 $n$，$N=nJ$，這樣的數據被叫做均衡數據 (balanced data)。
### 合並每個階層
過去常見的總結嵌套式數據的手段只是把每層數據取平均值，這樣的方法簡單粗暴但是偶爾是可以接受的，只要你能夠接受如此處理數據可能帶來的如下後果:
- 各層數據均值，其可靠程度 (方差) 隨着各層的樣本量不同而不同 (depends on the number of elementary units per cluster);
- 變量的含義發生改變。如果是使用層水平 (cluster level) 的數據，本來測量給個體的那些變量，就變成了**層的變量**，從此作出的任何統計學推斷，只能限制在層水平 (ecological fallacy, as correlations at the macro level cannot be used to make assertions at the micro level);
- 由於無視了層內個體數據，導致大量信息損失。
此處我們借用 [@Snijders1999] 書中第 28 頁的人造數據，如下表
<table class="table table-striped table-bordered" style="width: auto !important; margin-left: auto; margin-right: auto;">
<caption>Artificial data</caption>
<thead>
<tr>
<th style="text-align:center;"> Cluster $(j)$ </th>
<th style="text-align:center;"> id $(i)$ </th>
<th style="text-align:center;"> $X$ </th>
<th style="text-align:center;"> $Y$ </th>
<th style="text-align:center;"> $\bar{X}$ </th>
<th style="text-align:center;"> $\bar{Y}$ </th>
</tr>
</thead>
<tbody>
<tr>
<td style="text-align:center;"> 1 </td>
<td style="text-align:center;"> 1 </td>
<td style="text-align:center;"> 1 </td>
<td style="text-align:center;"> 5 </td>
<td style="text-align:center;"> 2 </td>
<td style="text-align:center;"> 6 </td>
</tr>
<tr>
<td style="text-align:center;"> 1 </td>
<td style="text-align:center;"> 2 </td>
<td style="text-align:center;"> 3 </td>
<td style="text-align:center;"> 7 </td>
<td style="text-align:center;"> 2 </td>
<td style="text-align:center;"> 6 </td>
</tr>
<tr>
<td style="text-align:center;"> 2 </td>
<td style="text-align:center;"> 1 </td>
<td style="text-align:center;"> 2 </td>
<td style="text-align:center;"> 4 </td>
<td style="text-align:center;"> 3 </td>
<td style="text-align:center;"> 5 </td>
</tr>
<tr>
<td style="text-align:center;"> 2 </td>
<td style="text-align:center;"> 2 </td>
<td style="text-align:center;"> 4 </td>
<td style="text-align:center;"> 6 </td>
<td style="text-align:center;"> 3 </td>
<td style="text-align:center;"> 5 </td>
</tr>
<tr>
<td style="text-align:center;"> 3 </td>
<td style="text-align:center;"> 1 </td>
<td style="text-align:center;"> 3 </td>
<td style="text-align:center;"> 3 </td>
<td style="text-align:center;"> 4 </td>
<td style="text-align:center;"> 4 </td>
</tr>
<tr>
<td style="text-align:center;"> 3 </td>
<td style="text-align:center;"> 2 </td>
<td style="text-align:center;"> 5 </td>
<td style="text-align:center;"> 5 </td>
<td style="text-align:center;"> 4 </td>
<td style="text-align:center;"> 4 </td>
</tr>
<tr>
<td style="text-align:center;"> 4 </td>
<td style="text-align:center;"> 1 </td>
<td style="text-align:center;"> 4 </td>
<td style="text-align:center;"> 2 </td>
<td style="text-align:center;"> 5 </td>
<td style="text-align:center;"> 3 </td>
</tr>
<tr>
<td style="text-align:center;"> 4 </td>
<td style="text-align:center;"> 2 </td>
<td style="text-align:center;"> 6 </td>
<td style="text-align:center;"> 4 </td>
<td style="text-align:center;"> 5 </td>
<td style="text-align:center;"> 3 </td>
</tr>
<tr>
<td style="text-align:center;"> 5 </td>
<td style="text-align:center;"> 1 </td>
<td style="text-align:center;"> 5 </td>
<td style="text-align:center;"> 1 </td>
<td style="text-align:center;"> 6 </td>
<td style="text-align:center;"> 2 </td>
</tr>
<tr>
<td style="text-align:center;"> 5 </td>
<td style="text-align:center;"> 2 </td>
<td style="text-align:center;"> 7 </td>
<td style="text-align:center;"> 3 </td>
<td style="text-align:center;"> 6 </td>
<td style="text-align:center;"> 2 </td>
</tr>
</tbody>
</table>
這個表中的人造數據，其結構是一目了然的，它的第二層級數量是 5，每層的個體數量是 2。這是一個平衡數據。由於這是個我們人爲模擬的數據，圖 \@ref(fig:artificialdata00) 也顯示它沒有隨機誤差，所有數據都在各自的直線上。
```{r artificialdata00, cache=TRUE, echo=TRUE, fig.asp=.7, fig.width=6, fig.cap='Artificial data: scatter of clustered data', fig.align='center', out.width='80%', message=FALSE, warning=FALSE}
dt <- read.csv("backupfiles/hierexample.csv", header = T)
names(dt) <- c("Cluster", "id", "X", "Y", "Xbar", "Ybar")
dt$Cluster <- as.factor(dt$Cluster)
ggthemr('fresh')
ggplot(dt, aes(x = X, y = Y, shape = Cluster, colour = Cluster)) + geom_point(size =5) +
theme(axis.text = element_text(size = 15),
axis.text.x = element_text(size = 15),
axis.text.y = element_text(size = 15)) +
labs(x = "X", y = "Y") +
theme(axis.title = element_text(size = 17), axis.text = element_text(size = 8)) +
theme(legend.position = "bottom", legend.direction = "horizontal") + theme(legend.text = element_text(size = 19), legend.title = element_text(size = 19))
```
- 如果我們無視其分層數據的嵌套式結構，把每個數據都看作是獨立的樣本，擬合一個**整體回歸 (total regression) 圖 \@ref(fig:artificialdata01)**:
$$
\hat Y_{ij} = 5.33 - 0.33 X_{ij}
$$
```{r artificialdata01, cache=TRUE, echo=TRUE, fig.asp=.7, fig.width=6, fig.cap='Artificial data: Total regression', fig.align='center', out.width='80%', message=FALSE, warning=FALSE}
ggthemr('fresh')
ggplot(dt, aes(x = X, y = Y)) + geom_point(size = 5, shape = 23) +
geom_smooth(method = "lm", se = FALSE, linetype = 2) +
theme(axis.text = element_text(size = 15),
axis.text.x = element_text(size = 15),
axis.text.y = element_text(size = 15)) +
labs(x = "X", y = "Y") +
theme(axis.title = element_text(size = 17), axis.text = element_text(size = 8)) +
theme(legend.position = "bottom", legend.direction = "horizontal") + theme(legend.text = element_text(size = 19), legend.title = element_text(size = 19))
```
- 如果我們只保留層級數據本身，求了變量 $X,Y$ 在每層的均值的話，就得到了**層間回歸 (between regression) 圖 \@ref(fig:artificialdata02)** -- 變量 $X,Y$ 之間的回歸直線的斜率變得更大了:
$$
\hat{\bar{Y}}_j = 8.0 - 1.0 \bar{X}_j
$$
```{r artificialdata02, cache=TRUE, echo=TRUE, fig.asp=.7, fig.width=6, fig.cap='Artificial data: scatter of clustered data', fig.align='center', out.width='80%', message=FALSE, warning=FALSE}
ggthemr('fresh')
ggplot(dt, aes(x = X, y = Y)) + geom_point(size =5, shape=23) +
geom_smooth(method = "lm", se = FALSE, linetype = 2) +
geom_abline(intercept = 8, slope = -1) +
geom_point(aes(x = Xbar, y=Ybar, size = 5)) +
theme(axis.text = element_text(size = 15),
axis.text.x = element_text(size = 15),
axis.text.y = element_text(size = 15)) +
labs(x = "X", y = "Y") +
theme(axis.title = element_text(size = 17), axis.text = element_text(size = 8)) + theme(legend.position="none")
```
### 生物學悖論 ecological fallacy
生物學悖論常見於我們認爲某分層數據中層級變量之間的關系，同樣適用與層級中的個體之間: 例如比較 A 國 和 B 國之間心血管疾病的發病率，發現 A 國國民食鹽平均攝入量高於 B 國，很多人可能就會下結論說食鹽攝入量高的個體，心血管疾病發病的危險度較高。然而，這樣的推論很多時候是錯誤的。
曾經在 [@Robinson1950] 論文中舉過的著名例子: 該研究調查美國每個州的移民比例，和該州相應的識字率之間的關系。研究者發現，移民比例較高的州，其識字率也較高 (相關系數 0.53)。由此就有人下結論說移民越多，那個州的教育水平會比較高。但是實際情況是，把每個個體的受教育水平和該個體本身是不是移民做了相關系數分析之後發現，這個關系其實是負相關 (-0.11)。所以說在州的水平作出的統計學推斷-移民多的州受教育水平高-是不正確的。之所以在州水平發現移民比例和受教育水平之間的正關系，是因爲移民傾向於居住在教育水平本來就比較高的本土出生美國人的州。
### 分解層級數據
如果是分析最初層級數據 (level 1) 的話，我們還需要考慮下列一些問題:
- 當心數據被多次利用
如果我們關心的變量其實是在第二層級的 (level 2/cluster level)，但是你卻把它當作是第一層級的數據，就會引起**數據很多**的錯覺，因爲同一層的個體他們的層屬變量都是一樣的，你擁有的數據其實並沒有你想的那麼多。
前文中用過的 GCSE 數據其實是一個很好的例子，下表中歸納了調查的學校類型 (男校，女校或者混合校)，以及按照每個學生個人所屬學校類型的總結，可以看出，當你嘗試使用個人 (elementary level) 水平的數據分析實際上是第二層級數據的特性時，你會被誤導。因爲個人數據告訴你， 34% 的學生在女校學習，然而正確的分析法應該是，學校中有 31% 的學校是女校。
```{r hier01tab02, echo=FALSE, cache=TRUE, eval=FALSE}
dt <- read.csv("backupfiles/hier01tab01.csv", header = T)
#names(dt) <- c("Model with", "sigma_u", "sigma_e", "sigma_u", "sigma_e")
kable(dt, "html", align = "c", caption = "Aggregated and disaggregated") %>%
kable_styling(bootstrap_options = c("striped", "bordered"),full_width = F, position = "center") %>%
add_header_above(c("School type" = 1,"Cluster Level" = 2, "Elementary Level" = 2))
```
<table class="table table-striped table-bordered" style="width: auto !important; margin-left: auto; margin-right: auto;">
<caption>Aggregated and disaggregated</caption>
<thead>
<tr>
<th style="text-align:center; border-bottom:hidden; padding-bottom:0; padding-left:3px;padding-right:3px;" colspan="1"><div style="border-bottom: 1px solid #ddd; padding-bottom: 5px;">School type</div></th>
<th style="text-align:center; border-bottom:hidden; padding-bottom:0; padding-left:3px;padding-right:3px;" colspan="2"><div style="border-bottom: 1px solid #ddd; padding-bottom: 5px;">Cluster Level</div></th>
<th style="text-align:center; border-bottom:hidden; padding-bottom:0; padding-left:3px;padding-right:3px;" colspan="2"><div style="border-bottom: 1px solid #ddd; padding-bottom: 5px;">Elementary Level</div></th>
</tr>
<tr>
<th style="text-align:center;"> </th>
<th style="text-align:center;"> N </th>
<th style="text-align:center;"> % </th>
<th style="text-align:center;"> N </th>
<th style="text-align:center;"> % </th>
</tr>
</thead>
<tbody>
<tr>
<td style="text-align:center;"> mixed </td>
<td style="text-align:center;"> 35 </td>
<td style="text-align:center;"> 54 </td>
<td style="text-align:center;"> 2169 </td>
<td style="text-align:center;"> 53 </td>
</tr>
<tr>
<td style="text-align:center;"> boys only </td>
<td style="text-align:center;"> 10 </td>
<td style="text-align:center;"> 15 </td>
<td style="text-align:center;"> 513 </td>
<td style="text-align:center;"> 13 </td>
</tr>
<tr>
<td style="text-align:center;"> girls only </td>
<td style="text-align:center;"> 20 </td>
<td style="text-align:center;"> 31 </td>
<td style="text-align:center;"> 1377 </td>
<td style="text-align:center;"> 34 </td>
</tr>
<tr>
<td style="text-align:center;"> Total </td>
<td style="text-align:center;"> 65 </td>
<td style="text-align:center;"> 100 </td>
<td style="text-align:center;"> 4059 </td>
<td style="text-align:center;"> 100 </td>
</tr>
</tbody>
</table>
- 分層數據分析法
有人會說，既然如此，那麼我們就把數據放在每層當中分析就好了 (stratified analyses)。還是用前文中用過的人造 5 層數據來說明這樣做的弊端。前面用了兩種方法 (total regression, between regression) 來總結這個 5 層的人造數據 \@ref(fig:artificialdata02)。最後一種分析此數據的方法是，把 5 層數據分開分別做回歸線如圖 \@ref(fig:artificialdata03)。等同於我們的對數據擬合五次下面的回歸方程:
$$
\hat Y_{ij} - \bar{Y}_j = \beta(X_{ij} - \bar{X}_j)
$$
這種模型被叫做**層內回歸 (within regression)**。這 5 個線性回歸的斜率都是 1，是五條不同截距的平行直線。因爲我們自己編造的數據的緣故，現實數據不太可能恰好所有層內回歸的斜率都是完全相同的。這其實也是曾內回歸法的一個默認前提 -- 每層數據中解釋變量和結果變量之間的關系是相同的。
```{r artificialdata03, cache=TRUE, echo=FALSE, fig.asp=.7, fig.width=6, fig.cap='Artificial data: within cluster regressions', fig.align='center', out.width='80%', message=FALSE, warning=FALSE}
ggthemr('fresh')
ggplot(dt, aes(x = X, y = Y)) + geom_point(size =5, shape=20) +
geom_line(aes(group = Cluster), lty = 2) +
# geom_smooth(method = "lm", se = FALSE, linetype = 2) +
# geom_abline(intercept = 8, slope = -1) +
# geom_point(aes(x = Xbar, y=Ybar, size = 5)) +
theme(axis.text = element_text(size = 15),
axis.text.x = element_text(size = 15),
axis.text.y = element_text(size = 15)) +
labs(x = "X", y = "Y") +
theme(axis.title = element_text(size = 17), axis.text = element_text(size = 8)) + theme(legend.position="none")
```
### 固定效應模型 fixed effect model
無論數據中的分層結構是否有現實意義 (如果說是五種不同的民族，那就有顯著的現實意義)，在回歸模型中都**有必要考慮這個分層結構對數據的變異的貢獻** (the contribution of the clusters to the data variation)。
線性回歸章節中我們使用的是五個啞變量來代表不同組別加以分析:
$$
Y_{ij} = \alpha_1 I_{i, j = 1} + \alpha_2 I_{i, j=2} + \cdots + \alpha_5 I_{i, j=5} + \beta_1X_{ij} + \varepsilon_{ij}
$$
其中 $j$ 是所屬層級編號。該模型中的 $\varepsilon_{ij}$ 被認爲服從均值爲零，方差爲 $\sigma_{\varepsilon}^2$ 的正態分布。該模型也可以簡寫爲:
$$
Y_{ij} = \alpha_j + \beta_1X_{ij} + \epsilon_{ij}
$$
一樣一預案
這樣的模型在等級線性回歸模型中被認爲是**固定效應模型 fixed effect model**。它其實是默認給五個層級五個不同的截距，每層內部 $X,Y$ 之間的關系 (斜率) 則被認爲是完全相同的 (namely the within cluster models are the same)。
本課剛開始的例子中有個數據是來自 6 所不同醫院 72 名患者的收縮期血壓的數據。我們現在來分析這些人中血壓和年齡之間的關系。下面的散點圖重現了六所醫院的72名患者的血壓和年齡。
```{r Hier01-7, echo=FALSE, fig.height=6, fig.width=7, fig.cap='SBP versus age: different symbols identify the six hospitals', fig.align='center', out.width='90%', cache=TRUE}
Bp <- read_dta("backupfiles/bp.dta")
Bp$hosp <- as.factor(Bp$hosp)
ggthemr('fresh')
ggplot(Bp, aes(x = age, y = bp, shape = hosp)) + geom_point(size =5) +
#geom_smooth(method = "lm", se = FALSE, linetype = 2) +
#geom_abline(intercept = 8, slope = -1) +
#geom_point(aes(x = Xbar, y=Ybar, size = 5)) +
theme(axis.text = element_text(size = 15),
axis.text.x = element_text(size = 15),
axis.text.y = element_text(size = 15)) +
labs(x = "Age", y = "Systolic blood pressure (mmHg)") +
theme(axis.title = element_text(size = 17), axis.text = element_text(size = 8)) + theme(legend.position="bottom")
```
下面在 R 裏擬合這個固定效應模型:
```{r Hier01-8, echo=TRUE, cache=TRUE}
Bp <- read_dta("backupfiles/bp.dta")
Bp$hosp <- as.factor(Bp$hosp)
Bp <- Bp %>%
mutate(c_age = age - mean(age))
# 通過指定截距爲零，獲取每個醫院的回歸線的截距
Model0 <- lm(bp ~ 0 + c_age + hosp, data = Bp)
summary(Model0)
# 先生成一個新的醫院變量 hops1 = 1。然後使用偏 F 檢驗法
# 檢驗控制了患者的年齡以後，這六所醫院的截距是否各自不相同。
Bp$hosp1 <- Bp$hosp[1]
mod2 <- lm(bp ~ 0 + c_age + as.numeric(hosp1), data = Bp)
anova(Model0, mod2)
```
偏 F 檢驗法給出的結果 $F(5, 65) = 2.16, P = 0.07$，所以說，數據其實告訴我們，調整了年齡之後，這六所醫院患者中年齡和血壓之間關系的回歸線有不同的截距。
## 簡單線性迴歸複習
滾回線性回歸章節 \@ref(lm)。
## 練習題
### 數據
1. High-School-and-Beyond 數據 <br> 本數據來自1982年美國國家教育統計中心 (National Center for Education Statistics, NCES) 對美國公立學校和天主教會學校的一項普查。曾經在 Hierarchical Linear Model [@Raudenbush2002] 一書中作爲範例使用。其數據的變量名和各自含義如下：
```
minority indicatory of student ethinicity (1 = minority, 0 = other)
female pupil's gender
ses standardized socio-economic status score
mathach measure of mathematics achievement
size school's total number of pupils
sector school's sector: 1 = catholic, 0 = not catholic
schoolid school identifier
```
2. PEFR 數據 <br> 數據本身是 17 名研究對象用兩種不同的測量方法測量兩次每個人的最大呼氣流速 (peak-expiratory-flow rate, PEFR)。最早在1986年的柳葉刀雜誌發表 [@Bland1986]。兩種測量法的名稱分別是 "Standard Wright" 和 "Mini Wright" peak flow meter。變量名和個字含義如下：
```
id participant identifier
wp1 standard wright measure at 1st occasion
wp2 standard wright measure at 2nd occasion
wm1 mini wright measure at 1st occasion
wm2 mini wright measure at 2nd occasion
```
### 問題
### 將 High-School-and-Beyond 數據導入 R 中，熟悉數據結構及內容，特別要注意觀察每個學校的學生特徵。
```{r hierex1-1, cache=TRUE}
hsb_selected <- read_dta("backupfiles/hsb_selected.dta")
length(unique(hsb_selected$schoolid)) ## number of school = 160
## create a subset data with only the first observation of each school
hsb <- hsb_selected[!duplicated(hsb_selected$schoolid), ]
## about 44 % of the schools are Catholic schools
with(hsb, tab1(sector, graph = FALSE, decimal = 2))
## among all the pupils, about 53% are females
with(hsb_selected, tab1(female, graph = FALSE, decimal = 2))
## among all the pupils, about 27.5% are from ethnic minorities
with(hsb_selected, tab1(minority, graph = FALSE, decimal = 2))
```
### 爲了簡便起見，接下來的分析只節選數據中前五所學校 188 名學生的數學成績，和 SES。分別計算每所學校的數學成績,及 SES 的平均值。
```{r hierex1-2, cache=TRUE}
hsb5 <- subset(hsb_selected, schoolid < 1320)
Mean_ses_math <- ddply(hsb5,~schoolid,summarise,mean_ses=mean(ses),mean_math=mean(mathach))
## the mean SES score ranges from -0.4255 to +0.5280
## the mean Maths score ranges from 7.636 to 16.255
Mean_ses_math
```
### 先無視掉學校這一分層變量，把所有學生看作是相互獨立的，擬合總體的 SES 和數學成績的線性迴歸 **(Total regression model)**。把該總體模型的預測值提取並存儲在數據庫中。
```{r mathses, cache=TRUE, echo=TRUE, fig.asp=.7, fig.width=8, fig.cap='Scatter plot of SES and math achievements among all pupils from first 5 schools, assuming that they are all independent', fig.align='center', out.width='80%', message=FALSE, warning=FALSE}
## plot the scatter of mathach and ses among these 5 schools
ggplot(hsb5, aes(x = ses, y = mathach)) + geom_point() +
theme_bw() +
theme(axis.text = element_text(size = 15),
axis.text.x = element_text(size = 15),
axis.text.y = element_text(size = 15)) +
labs(x = "SES", y = "Math achievement") +
xlim(-2.05, 2.05)+
ylim(-10, 30) +
theme(axis.title = element_text(size = 17), axis.text = element_text(size = 8))
```
```{r hierex1-3, cache = TRUE}
Total_reg <- lm(mathach ~ ses, data = hsb5)
## the total regression model gives an estimated regression coefficient for the SES
## of each pupil equal to 3.31 (SE=0.66)
summary(Total_reg)
hsb5$Pred_T <- Total_reg$fitted.values # save the fitted values to the dataset
```
### 用各個學校 SES 和數學成績的均值擬合一個學校間的線性迴歸模型 **(between regression model)**。
```{r hierex1-4, cache=TRUE}
Btw_reg <- lm(mean_math ~ mean_ses, data = Mean_ses_math)
## the regression model for the school level variables (between model) gives
## an estimated regression coefficient of 7.29 (SE=1.41)
summary(Btw_reg)
Mean_ses_math$Pred_B <- Btw_reg$fitted.values # save the fitted values to the dataset
```
### 分別對每個學校內的學生進行 SES 和數學成績擬合線性迴歸模型。
```{r hierex1-5, cache=TRUE}
Within_schl1 <- lm(mathach ~ ses, data = hsb5[hsb5$schoolid == 1224,])
Within_schl2 <- lm(mathach ~ ses, data = hsb5[hsb5$schoolid == 1288,])
Within_schl3 <- lm(mathach ~ ses, data = hsb5[hsb5$schoolid == 1296,])
Within_schl4 <- lm(mathach ~ ses, data = hsb5[hsb5$schoolid == 1308,])
Within_schl5 <- lm(mathach ~ ses, data = hsb5[hsb5$schoolid == 1317,])
# the within school regressions gives estimated slopes which have a mean of 1.65
# and which ranges between 0.126 and 3.255
summary(c(Within_schl1$coefficients[2], Within_schl2$coefficients[2],
Within_schl3$coefficients[2], Within_schl4$coefficients[2],
Within_schl5$coefficients[2]))
# the SEs ranging between 1.21 and 3.00
summary(c(summary(Within_schl1)$coefficients[4],
summary(Within_schl2)$coefficients[4],
summary(Within_schl3)$coefficients[4],
summary(Within_schl4)$coefficients[4],
summary(Within_schl5)$coefficients[4]))
hsb5$Pred_W <- c(Within_schl1$fitted.values, Within_schl2$fitted.values,
Within_schl3$fitted.values, Within_schl4$fitted.values,
Within_schl5$fitted.values) ## save the predicted value into the dataset
```
### 比較三種模型計算的數學成績的擬合值，他們一致？還是有所不同？爲什麼會有不同？
- 總體模型 (Total regression model) 實際上無視了學生的性別，種族等可能帶來的混雜效果；
- 學校間模型 (Between model) 估計的實際上是**SES均值**每增加一個單位，與之對應的**數學平均成績**的改變量，**這個模型絕對不可用與評估個人的 SES 與數學成績之間的關係**；
- 學校內模型 (Within model) 擬合的 SES 與數學成績之間的關係變得十分地不精確 (SEs are fairly large)，變化幅度也很大。
### 把三種模型的數學成績擬合值散點圖繪製在同一張圖內。
```{r mathses-3models, cache=TRUE, echo=TRUE, fig.height=6.5, fig.width=8, fig.cap='High-school-and-beyond data: Predicted values by Total, Between, and Within regression models', fig.align='center', out.width='80%', message=FALSE, warning=FALSE}
Mean <- Mean_ses_math[, 1:3]
names(Mean) <- c("schoolid", "ses", "Pred_W")
ggplot(hsb5, aes(x = ses, y = Pred_W, group = schoolid)) +
geom_line(linetype = 2, size = 1) +
geom_abline(intercept = Total_reg$coefficients[1], slope = Total_reg$coefficients[2],
colour = "dark blue") +
geom_abline(intercept = Btw_reg$coefficients[1], slope = Btw_reg$coefficients[2],
colour = "red") +
geom_point(data = Mean, shape = 17, size = 4, colour = "Red") +
theme_bw() +
theme(axis.text = element_text(size = 15),
axis.text.x = element_text(size = 15),
axis.text.y = element_text(size = 15)) +
labs(x = "SES", y = "Fitted regression lines (Maths achievement)") +
xlim(-2.05, 2.05)+
ylim(5, 20) +
theme(axis.title = element_text(size = 17), axis.text = element_text(size = 8)) +
theme(plot.caption = element_text(size = 12,
hjust = 0)) + labs(caption = "Black dash line: Within regression model;
Blue solid line: Total regression model;
Red solid line: Between regression model;
Red triangle: School mean values")
```
### 用這 5 個學校的數據擬合一個固定效應線性迴歸模型
```{r hierex1-8, cache=TRUE}
Fixed_reg <- lm(mathach ~ ses + factor(schoolid), data = hsb5)
## Fitting a fixed effect model to these data is equivalent to forcing
## a common slope onto the five within regression models. It gives an
## estimated slope of 1.789 (SE=0.76), close to their average of 1.64799.
## Note that controlling for female, minority, and sector but not for
## schoolid leads to roughly the same estimate (slope = 1.68, SE=0.75)
summary(Fixed_reg)
summary(lm(mathach ~ ses + female + minority + sector, data = hsb5))
```
### 讀入 PEFR 數據。
```{r hierex1-9, cache=TRUE}
pefr <- read_dta("backupfiles/pefr.dta")
# the data are in wide format
pefr
# transform data into long format
pefr_long <- pefr %>%
gather(key, value, -id) %>%
separate(key, into = c("measurement", "occasion"), sep = 2) %>%
arrange(id, occasion) %>%
spread(measurement, value)
pefr_long
```
```{r tworecordings, cache=TRUE, echo=TRUE, fig.height=6, fig.width=9, fig.cap='Two recordings of PEFR taken with the standard Wright meter', fig.align='center', out.width='80%', message=FALSE, warning=FALSE}
## figure shows slightly closer agreement between the repeated measures of standard Wright,
## than between those of Mini Wright
ggplot(pefr_long, aes(x = id, y = wp, fill = occasion)) +
geom_point(size = 4, shape = 21) +
geom_hline(yintercept = mean(pefr_long$wp), colour = "red") +
theme_bw() +
scale_x_continuous(breaks = 1:17)+
theme(axis.text = element_text(size = 15),
axis.text.x = element_text(size = 15),
axis.text.y = element_text(size = 15)) +
labs(x = "Subject ID", y = "W Measurements") +
theme(axis.title = element_text(size = 17), axis.text = element_text(size = 8))+
theme(legend.text = element_text(size = 19),
legend.title = element_text(size = 19))
```
### 求每個患者的 `wp` 兩次測量平均值
```{r hierex1-11, cache=TRUE}
# the means range from 171.5 to 644.5
pefr_long %>%
group_by(id) %>%
summarise(mean_wp = mean(wp))
```
### 在 R 裏先用 ANOVA 分析個人的 `wp` 變異。再用 `lme4::lmer` 擬合用 `id` 作隨機效應的混合效應模型。確認後者報告的 `Std.Dev for id effect` 其實可以用 ANOVA 結果的 $\sqrt{\frac{\text{MMS-MSE}}{n}}$ (n 是每個個體重複測量值的個數)。
```{r hierex1-12, cache=TRUE, message=FALSE}
with(pefr_long, anova(lm(wp~factor(id))))
#library(lme4)
( fit <- lmer(wp ~ (1|id), data=pefr_long) )
sqrt((27600 - 234)/2)
```
### 擬合結果變量爲 `wp`，解釋變量爲 `id` 的簡單線性迴歸模型。用數學表達式描述這個模型。
```{r hierex1-13, cache=TRUE, message=FALSE}
Reg <- lm(wp ~ factor(id), data = pefr_long)
# The fixed effect regression model leads to the same ANOVA
# table. To the same estimate of the residual SD = (15.307)
# However, it does not give an estimate of the "SD of id effect"
# Instead it gives estimates of mean PEFR for participant number 1
# = 492 and estimates of the difference in means from him/her
# for all the other 16 pariticipants
anova(Reg)
summary(Reg)
```
上面的模型用數學表達式來描述就是：
$$
\begin{aligned}
Y_{ij} & = \alpha_1 + \delta_i + \varepsilon_{ij} \\
\text{Where } \delta_j & = \alpha_j - \alpha_1 \\
\text{and } \delta_1 & = 0
\end{aligned}
$$
### 將 `wp` 中心化之後，重新擬合相同的模型，把截距去除掉。寫下這個模型的數學表達式。
```{r hierex1-14, cache=TRUE, message=FALSE}
Reg1 <- lm((wp - mean(wp)) ~ 0 + factor(id), data = pefr_long)
# it leads to the same ANOVA table again, same residual SD
anova(Reg1)
summary(Reg1)
```
上面的模型用數學表達式來描述就是：
$$
\begin{aligned}
Y_{ij} - \mu & = \gamma_j + \varepsilon_{ij} \\
Y_{ij} & = \mu + \gamma_j + \varepsilon_{ij} \\
\text{Where } \mu & \text{ is the overall mean} \\
\text{and } \sum_{j=1}^J\gamma_j & = 0\\
\end{aligned}
$$
### 計算這些迴歸係數 (其實是不同羣之間的隨機截距) 的均值和標準差。
```{r hierex1-15, cache=TRUE, message=FALSE}
# the individual level intercepts have mean zero and SD = 117.47, larger than the estimated
# Std.Dev for id effect.
Reg1$coefficients
mean(Reg1$coefficients)
sd(Reg1$coefficients)
```
# 隨機截距模型 random intercept model {#random-intercept}
最簡單的隨機效應模型 -- 隨機截距模型 random intercept model。
## 隨機截距模型的定義
有時，我們對每個分層各自的截距大小並不那麼感興趣，且如果只有固定效應的話，其實我們從某種程度上忽略掉了數據層與層之間變異的方差 (between cluster variation)。於是，在模型中考慮這些問題的解決方案就是 -- 我們讓各層的截距呈現隨機效應 (treat the variation in cluster intercepts not as fixed)，**把這些截距視爲來自與某種分布的隨機呈現 (randomly draws from some distribution)**。於是原先的只有固定效應部分的模型，就增加了隨機截距部分:
$$
\begin{aligned}
Y_{ij} & = \mu + u_j + \varepsilon_{ij} \\
\text{where } u_j & \stackrel{N.I.D}{\sim} N(0, \sigma_u^2) \\
\varepsilon_{ij} & \stackrel{N.I.D}{\sim} N(0, \sigma_\varepsilon^2) \\
u_j & \text{ are independent from } \varepsilon_{ij} \\
\end{aligned}
(\#eq:hier2-1)
$$
這個混合效應模型中，
- $\mu$ 是總體均值;
- $u_j$ 是一個服從均值 0, 方差 (the population between cluster variance) 爲 $\sigma_u^2$ 的正態分布的隨機變量;
- $\varepsilon_{ij}$ 是隨機誤差，它也被認爲服從均值爲 0, 方差爲 $\sigma_\varepsilon^2$ 的正太分布，且這兩個隨機效應部分之間也是**相互獨立的**。
- 從該模型估算的結果變量 $Y_{ij}$ 的方差是 $\sigma_u^2 + \sigma_\varepsilon^2$。
- 隨機截距模型又被叫做是 **方差成分模型 (variance-component model)**，或者是**單向隨機效應方差模型 (one-way random effects ANOVA model)**。
這個模型和僅有固定效應的模型，有顯著的不同:
$$
Y_{ij} = \mu + \gamma_j + \varepsilon_{ij}
$$
固定效應模型裏，
- $\mu$ 也是總體均值;
- $\sum_{j=1}^J \gamma_j = 0$ 是**將各組不同截距之和強制爲零**的過程;
所以隨機截距模型打破了這個限制，使得隨機的截距 $\mu_j$ 成爲一個服從均值爲 0，方差爲 $\sigma_u^2$ 的 **隨機變量**。
隨機效應部分 $u_j$ 和隨機誤差 $\varepsilon_{ij}$ 之間相互獨立的前提，意味着兩個裏屬於不同層級的觀察之間是相互獨立的，但是反過來，同屬於一個層級的個體之間就變成了有相關性的了 (within cluster correlation):
$$
\begin{aligned}
\because Y_{1j} & = \mu + u_j + \varepsilon_{1j} \\
Y_{2j} & = \mu + u_j + \varepsilon_{2j} \\
\therefore \text{Cov}(Y_{1j}, Y_{2j}) & = \text{Cov}(u_j, u_j) + \text{Cov}(u_j, \varepsilon_{2j}) + \text{Cov}(\varepsilon_{1j}, u_j) + \text{Cov}(\varepsilon_{1j}, \varepsilon_{2j}) \\
& = \text{Cov}(u_j, u_j) = V(u_j, u_j)\\
& = \sigma_u^2
\end{aligned}
$$
由於 $\text{Var}(Y_{1j}) = \text{Var}(Y_{2j}) = \sigma_u^2 + \sigma_\varepsilon^2$，所以，同屬一層的兩個個體之間的**層內相關系數 (intra-class correlation)**:
$$
\lambda = \frac{\text{Cov}(Y_{1j}, Y_{2j})}{\text{SD}(Y_{1j})\text{SD}(Y_{2j})} = \frac{\sigma_u^2}{\sigma_\varepsilon^2 + \sigma_u^2}
$$
從層內相關系數的公式也可看出，該相關系數可以同時被理解爲結果變量 $Y_{ij}$ 的方差中歸咎與層(cluster)結構的部分的百分比。
This is the within-cluster or intra-class correlation, that we will denote $\lambda$. Note that it is also the proportion of total variance that is accounted for by the cluster.
## 隨機截距模型的參數估計
如此，我們就知道在隨機截距模型裏，有三個需要被估計的參數 $\mu, \sigma_u^2, \sigma^2_\varepsilon$。我們可以利用熟悉的極大似然估計法估計這些參數 (Maximum Likelihood, ML)。當且進當嵌套式結構數據是**平衡數據 (balanced)**時 (即，每層中的個體數量相同)，這三個參數的 $\text{MLE}$ 分別是:
$$
\begin{aligned}
\hat\mu & = \bar{Y} \\
\hat\sigma_\varepsilon^2 & = \text{Mean square error, MSE} \\
\hat\sigma_u^2 & = \frac{\text{Model Sum of Squares, MSS}}{Jn} - \frac{\hat\sigma^2_\varepsilon}{n}
\end{aligned}
(\#eq:hier02-2)
$$
只要模型指定正確無誤，前兩個極大似然估計是他們各自的無偏估計。但第三個，也就是層內方差的估計量確實際上是低估了的 (downward biased)。這裏常用的另一種對層內方差參數的估計法被叫做**矩估計量 (moment estimator, or ANOVA estimator)**:
$$
\begin{aligned}
\widetilde{\sigma}_u^2 & = \frac{\text{MSS}}{(J-1)n}- \frac{\hat\sigma_\varepsilon^2}{n} \\
& = \frac{\text{MSS} - \text{MSE}(J-1)}{(J-1)n} \\
& = \frac{\text{MMS}(J-1) - \text{MSE}(J-1)}{(J-1)n} \\
& = \frac{\text{MMS} - \text{MSE}}{n}
\end{aligned}
$$
對於平衡數據 (balanced data)，這個矩估計量又被叫做**限制性極大似然 (Restricted Maximum Likelihood, REML)**。限制性極大似然法，是一個真極大似然過程 (genuine maximum likelihood procedure)，但是它每次進行估計的時候，會先"去除掉"固定效應部分，所以每次用於估計參數的數據其實是對數據的線性轉換後 $Y_{ij} - \mu = u_j + \varepsilon_{ij}$，它使用的數據是這個等式右半部分的轉換後數據。在 REML 過程中，先估計層內方差 $\sigma_u^2$ 再對固定效應部分的總體均值估計，所以是個兩步走的過程。另外除了這裏討論的 ML, REML這兩種對層內方差進行參數估計的方法之外，在計量經濟學 (econometrics) 中常用的是 (本課不深入探討) **廣義最小二乘法 (Generalized Least Squares, GLS)**。GLS 使用的是一種加權的最小二乘法 (OLS)，該加權法根據層與隨機誤差的方差成分 (variance components) 不同而給不同的層以不同的截距權重。當數據本身是平衡數據時，GLS給出的估計結果等同於 REML法。當數據不是平衡數據的時候，ML/REML 其實背後使用的原理也是 GLS。
## 如何在 R 中進行隨機截距模型的擬合
在 R 或 STATA 中擬合隨機截距模型，需要數據爲“長 (long)” 數據，下面的代碼可以在 R 裏面把 “寬 (wide)” 的數據調整成爲 **長** 數據:
```{r Hier02-1, cache=TRUE, echo=TRUE, message=FALSE, warning=FALSE}
pefr <- read_dta("backupfiles/pefr.dta")
# the data are in wide format
head(pefr)
# transform data into long format
pefr_long <- pefr %>%
gather(key, value, -id) %>%
separate(key, into = c("measurement", "occasion"), sep = 2) %>%
arrange(id, occasion) %>%
spread(measurement, value)
pefr_long
```
在 R 裏面，有兩個包 (`lme4::lmer` 或 `nlme::lme`) 的各自兩種代碼以供選用:
```{r Hier02-2, cache=TRUE, echo=TRUE, message=FALSE, warning=FALSE}
M0 <- lme(fixed = wm ~ 1, random = ~ 1 | id, data = pefr_long, method = "REML")
summary(M0)
M1 <- lmer(wm ~ (1|id), data = pefr_long, REML = TRUE)
summary(M1)
```
不知道爲什麼在 R 裏有這兩種完全不同的方式來擬合混合效應模型。還好他們的結果基本完全一致。在這個極爲簡單的例子裏，我們可以利用模型擬合的結果中 `Random effects` 的部分來計算**層內相關系數 (intra-class correlation)**:
$$
\hat\lambda = \frac{\hat\sigma_u^2}{(\hat\sigma_u^2 + \hat\sigma_\varepsilon^2)} = \frac{110.40^2}{110.40^2 + 19.91^2} = 0.97
$$
這是對 Mini Wright meter 測量方法可靠性的一個評價指標。其中 $\sigma_u^2$ 是患者最大呼吸速率 (PEFR) 測量值的方差，$\sigma_\varepsilon^2$ 是測量的隨機誤差，所以這裏的測量方法的可靠度是 97%，是可信度十分高的測量準確度。
## 隨機截距模型中的統計推斷