-
Notifications
You must be signed in to change notification settings - Fork 256
/
chapter_10.html
999 lines (860 loc) · 37.8 KB
/
chapter_10.html
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
<!DOCTYPE html>
<html lang="" xml:lang="">
<head>
<title>chapter_10.knit</title>
<meta charset="utf-8" />
<meta name="author" content="" />
<script src="libs/header-attrs-2.20/header-attrs.js"></script>
<link href="libs/remark-css-0.0.1/default.css" rel="stylesheet" />
<link href="libs/panelset-0.2.6/panelset.css" rel="stylesheet" />
<script src="libs/panelset-0.2.6/panelset.js"></script>
<script src="libs/htmlwidgets-1.5.4/htmlwidgets.js"></script>
<link href="libs/datatables-css-0.0.0/datatables-crosstalk.css" rel="stylesheet" />
<script src="libs/datatables-binding-0.27/datatables.js"></script>
<script src="libs/jquery-3.6.0/jquery-3.6.0.min.js"></script>
<link href="libs/dt-core-1.12.1/css/jquery.dataTables.min.css" rel="stylesheet" />
<link href="libs/dt-core-1.12.1/css/jquery.dataTables.extra.css" rel="stylesheet" />
<script src="libs/dt-core-1.12.1/js/jquery.dataTables.min.js"></script>
<link href="libs/crosstalk-1.2.0/css/crosstalk.min.css" rel="stylesheet" />
<script src="libs/crosstalk-1.2.0/js/crosstalk.min.js"></script>
<link rel="stylesheet" href="css/zh-CN.css" type="text/css" />
<link rel="stylesheet" href="css/Custumed_Style.css" type="text/css" />
</head>
<body>
<textarea id="source">
class: center, middle
<span style="font-size: 50px;">**第十章**</span> <br>
<span style="font-size: 50px;">回归模型(三)：广义线性模型</span> <br>
<span style="font-size: 30px;">胡传鹏</span> <br>
<span style="font-size: 20px;"> </span> <br>
<span style="font-size: 30px;">2024-05-08</span> <br>
<span style="font-size: 20px;"> Made with Rmarkdown</span> <br>
---
<style type="text/css">
.bigfont {
font-size: 30px;
}
.size5{
font-size: 24px;
}
.tit_font{
font-size: 60px;
}
</style>
# 补充内容: easystats系统包的简介
<img src="./picture/chp10/bilibili.png" width="60%" style="display: block; margin: auto;" />
<br>
<br>
<center>
https://www.bilibili.com/video/BV1rz421D7iJ/?spm_id_from=333.337.search-card.all.click
</center>
---
class: center, middle
.tit_font[
当因变量不服从正态分布(如正确率)时如何处理？
]
---
.panelset[
.panel[.panel-name[df.match]
```r
head(df.match[c(3,11:17)],5) %>% DT::datatable()
```
<div id="htmlwidget-db93ea4216a2eb01ef04" style="width:100%;height:auto;" class="datatables html-widget"></div>
<script type="application/json" data-for="htmlwidget-db93ea4216a2eb01ef04">{"x":{"filter":"none","vertical":false,"data":[["1","2","3","4","5"],[7304,7304,7304,7304,7304],["moral","immoral","moral","immoral","moral"],["Other","Other","Self","Other","Self"],["moralOther","immoralOther","moralSelf","immoralOther","moralSelf"],["match","match","match","match","match"],["n","n","n","n","n"],["n","m","n","m","n"],[1,0,1,0,1]],"container":"<table class=\"display\">\n <thead>\n <tr>\n <th> <\/th>\n <th>Sub<\/th>\n <th>Valence<\/th>\n <th>Identity<\/th>\n <th>Label<\/th>\n <th>Match<\/th>\n <th>CorrResp<\/th>\n <th>Resp<\/th>\n <th>ACC<\/th>\n <\/tr>\n <\/thead>\n<\/table>","options":{"columnDefs":[{"className":"dt-right","targets":[1,8]},{"orderable":false,"targets":0}],"order":[],"autoWidth":false,"orderClasses":false}},"evals":[],"jsHooks":[]}</script>
.panel[.panel-name[df.match.aov]
```r
df.match.aov %>%
dplyr::select(1:4) %>%
head(5) %>%
DT::datatable()
```
<div id="htmlwidget-751dc3e27421b5184726" style="width:100%;height:auto;" class="datatables html-widget"></div>
<script type="application/json" data-for="htmlwidget-751dc3e27421b5184726">{"x":{"filter":"none","vertical":false,"data":[["1","2","3","4","5"],[7304,7304,7304,7304,7305],["moral","moral","immoral","immoral","moral"],["Self","Other","Self","Other","Self"],[0.986666666666667,0.833333333333333,0.821917808219178,0.71830985915493,0.932432432432432]],"container":"<table class=\"display\">\n <thead>\n <tr>\n <th> <\/th>\n <th>Sub<\/th>\n <th>Valence<\/th>\n <th>Identity<\/th>\n <th>mean_ACC<\/th>\n <\/tr>\n <\/thead>\n<\/table>","options":{"columnDefs":[{"className":"dt-right","targets":[1,4]},{"orderable":false,"targets":0}],"order":[],"autoWidth":false,"orderClasses":false}},"evals":[],"jsHooks":[]}</script>
]]]
---
.tit_font[ Contents</h1>]
<br>
<br>
.bigfont[
10.1 广义线性模型<br>
<br>
10.2 二项分布<br>
<br>
10.3 其他分布<br>
<br>
10.4 代码实操<br>
<br>
10.5 方法比较<br>
]
---
#10.1 广义线性模型(Generalized Linear Model, GLM)
##多元线性回归（Multiple Linear Regression）
<br>
.normal[
`$$Y = b_0 + b_{1}X_{1} + b_{2}X_{2} +... + b_{p}X_{p} + \epsilon$$`
- `\(Y\)` : 因变量，Dependent variable
- `\(X_i\)` : 自变量，Independent (explanatory) variable
- `\(b_0\)` : 截距，Intercept
- `\(b_i\)` : 斜率，Slope
- `\(\epsilon\)` : 残差，Residual (error)
]
---
#10.1 广义线性模型(Generalized Linear Model, GLM)
##线性模型的组成部分
---
#10.1 广义线性模型(Generalized Linear Model, GLM)
##线性模型的组成部分
<img src="./picture/chp10/plot.png" width="60%" style="display: block; margin: auto;" />
---
#10.1 广义线性模型(Generalized Linear Model, GLM)
##回归方程的多种形式<br>
.bigfont[
- 简单线性回归：
`$$Y = b_0+b_1 X_1+ b_2 X_2+…+b_p X_p + \epsilon$$`
- 线性代数表达：
`$$y_i = b_0 + b_1 X_{i1} + b_2 X_{i2} + … + b_p X_{ip} + \epsilon$$`
- 矩阵表达：
`$$Y= X\beta + \epsilon$$`
- 代码表达(r)：
`$$Y \sim X_1 + X_2 + ... + X_n$$`
]
---
#10.1 广义线性模型(Generalized Linear Model, GLM)
##回归方程的多种表达形式<br>
<br>
.bigfont[
- 回归模型形式：观测项 = 预测项 + 误差项 <br>
- 假定观测项是正态分布，上述公式可以重新表达为： <br>
`$$y \sim N(\mu, \epsilon)$$`
- 其中，$\mu$为预测值，即
`$$μ = \beta_0 + \beta_1 x$$`
- 观测值服从以预测项为均值的**正态分布**，观测值与预测值之间的差值就是残差。<br>
]
--
.bigfont[
如果因变量不服从正态分布，如何构建回归模型？
]
---
#10.1 广义线性模型(Generalized Linear Model, GLM)
---
#10.1 广义线性模型(Generalized Linear Model, GLM)
---
#10.1 广义线性模型(Generalized Linear Model, GLM)
---
#10.1 广义线性模型(Generalized Linear Model, GLM)
##Generalized Linear Model, GLM
###在简单线性回归中，预测项的连接函数等于它本身
---
#10.1 广义线性模型(Generalized Linear Model, GLM)
.bigfont[
- 简单线性模型可视为GLM的特殊形式，预测项的连接函数等于它本身，观测项为正态分布。
- 在广义线性模型中：
- 观测项不一定是正态分布（残差不一定是正态分布）
- 连接函数不等于其自身
- 广义线性模型，能够对非正态分布的因变量进行建模
]
---
#10.2 二项分布(Binomial Distribution)
##伯努利试验
<br>
.bigfont[
- 同样的条件下重复地、相互独立地进行的一种随机试验。
<br>
<br>
- 该随机试验只有两种可能结果：发生或者不发生。
<br>
<br>
- 假设该项试验独立重复地进行了n次，那么就称这一系列重复独立的随机试验为n重伯努利试验(n-fold bernoulli trials)。
<br>
<br>
- n次独立重复的伯努利试验的概率分布服从二项分布
]
---
#10.2 二项分布(Binomial Distribution)
.bigfont[
- 每次试验中事件A发生的概率为p
<br><br>
- X表示n重伯努利试验中事件A发生的次数，X的可能取值为0，1，…，n
<br><br>
- 对每一个k（0 ≤ k ≤ n）,事件{X = k} 指”n次试验中事件A恰好发生k次”
<br><br>
- 随机变量X服从以n, p为参数的二项分布，写作 `\(X \sim B(n, p)\)`
<br><br>
- `\(p \in [0,1]\)`, `\(n \in N\)`
`$$P(X=k )=𝐶_𝑛^𝑘 𝑝^𝑘 𝑞^{𝑛−𝑘}= 𝐶_𝑛^𝑘 𝑝^𝑘 (1−𝑝)^{𝑛−𝑘}$$`
$$𝐶_𝑛^𝑘= 𝑛!/𝑘!(𝑛−𝑘)! $$
]
---
#10.2 二项分布(Binomial Distribution)
##抛硬币
.panelset[
.panel[.panel-name[5人，每人10次]
```r
simulate_coin_toss(prob_head = 0.5,num_people = 5, num_tosses = 10)
```
<img src="chapter_10_files/figure-html/unnamed-chunk-9-1.png" width="540" style="display: block; margin: auto;" />
.panel[.panel-name[10人，每人10次]
```r
simulate_coin_toss(prob_head = 0.5,num_people = 10, num_tosses = 10)
```
<img src="chapter_10_files/figure-html/unnamed-chunk-10-1.png" width="540" style="display: block; margin: auto;" />
.panel[.panel-name[1000人，每人10次]
```r
simulate_coin_toss(prob_head = 0.5,num_people = 1000, num_tosses = 10)
```
<img src="chapter_10_files/figure-html/unnamed-chunk-11-1.png" width="540" style="display: block; margin: auto;" />
]]]]
---
#10.2 二项分布(Binomial Distribution)
.bigfont[
- 已知一次试验中的每次尝试中事件A发生的概率$p$，共进行$n$次独立重复的伯努利试验
- 事件A在一次试验中出现k次，事件A在n次试验中出现次数的平均数
`$$（𝑘_1+𝑘_2+𝑘_3+...+𝑘_𝑛/𝑛)$$`
- 当n → ∞，$p$ ≠ q，$np$ ≥ 5且$nq$ ≥ 5，事件A在$n$次试验中出现次数的平均数：
`$$\mu = np$$`
- 事件A出现次数所属分布的标准差：
$$ \sigma = \sqrt{𝑛𝑝𝑞}$$
]
---
#10.2 二项分布(Binomial Distribution)
## 如何将$z$与二分变量进行连接？
### (1)将预测项映射到(0,1)之间，例如，使用
`$$\frac{1}{1+exp(-z)}$$`
### (2)找到一个分布，能根据(0,1)之间的值转成二分变量，例如，伯努利分布。
.pull-left[
]
.pull-right[
]
---
#10.2 二项分布(Binomial Distribution)
<img src="./picture/chp10/func2.png" width="80%" style="display: block; margin: auto;" />
---
#10.2 二项分布(Binomial Distribution)
##参数求解
.bigfont[
- 对于logit回归，我们可以使用极大似然估计对其进行求解
- 该求解过程比较复杂，一般由计算机自动完成
]
---
#10.3 代码实操
##首先分析一个被试的数据
```r
df.match.7304 <- df.match %>%
dplyr::filter(Sub == 7304) #选择被试7304
mod_7304_full <- stats::glm(data = df.match.7304, #数据
formula = ACC ~ 1 + Identity * Valence, #模型
family = binomial) #因变量为二项分布
summary(mod_7304_full) %>% #查看模型信息
capture.output() %>% .[c(6:11,15:19)] #课堂展示重要结果
```
```
## [1] "Coefficients:"
## [2] " Estimate Std. Error z value Pr(>|z|) "
## [3] "(Intercept) 4.30 1.01 4.28 0.000019 ***"
## [4] "IdentityOther -2.69 1.05 -2.56 0.0106 * "
## [5] "Valenceimmoral -2.77 1.05 -2.64 0.0083 ** "
## [6] "IdentityOther:Valenceimmoral 2.10 1.13 1.86 0.0628 . "
## [7] "(Dispersion parameter for binomial family taken to be 1)"
## [8] ""
## [9] " Null deviance: 254.02 on 290 degrees of freedom"
## [10] "Residual deviance: 228.32 on 287 degrees of freedom"
## [11] "AIC: 236.3"
```
---
#10.3 代码实操
.panelset[
.panel[.panel-name[mod_null]
```r
#无固定效应
mod_null <- lme4::glmer(data = df.match, #数据
formula = ACC ~ (1 + Identity * Valence|Sub), #模型
family = binomial) #因变量二项分布
#performance::model_performance(mod_null)
summary(mod_null) %>%
capture.output()%>% .[c(7:8,14:24)]
```
```
## [1] " AIC BIC logLik deviance df.resid "
## [2] " 9378.8 9460.2 -4678.4 9356.8 11999 "
## [3] "Random effects:"
## [4] " Groups Name Variance Std.Dev. Corr "
## [5] " Sub (Intercept) 1.53 1.24 "
## [6] " IdentityOther 2.52 1.59 -0.86 "
## [7] " Valenceimmoral 2.40 1.55 -0.85 0.83 "
## [8] " IdentityOther:Valenceimmoral 3.33 1.83 0.69 -0.87 -0.82"
## [9] "Number of obs: 12010, groups: Sub, 41"
## [10] ""
## [11] "Fixed effects:"
## [12] " Estimate Std. Error z value Pr(>|z|) "
## [13] "(Intercept) 2.014 0.114 17.7 <0.0000000000000002 ***"
```
.panel[.panel-name[mod]
```r
#随机截距，固定斜率
mod <- lme4::glmer(data = df.match, #数据
formula = ACC ~ 1 + Identity * Valence + (1|Sub), #模型
family = binomial) #因变量二项分布
#performance::model_performance(mod)
summary(mod) %>%
capture.output() %>% .[c(7:8,14:24,28:32)]
```
```
## [1] " AIC BIC logLik deviance df.resid "
## [2] " 9639.0 9675.9 -4814.5 9629.0 12005 "
## [3] "Random effects:"
## [4] " Groups Name Variance Std.Dev."
## [5] " Sub (Intercept) 0.237 0.487 "
## [6] "Number of obs: 12010, groups: Sub, 41"
## [7] ""
## [8] "Fixed effects:"
## [9] " Estimate Std. Error z value Pr(>|z|)"
## [10] "(Intercept) 2.4964 0.1015 24.60 < 0.0000000000000002"
## [11] "IdentityOther -0.7161 0.0839 -8.53 < 0.0000000000000002"
## [12] "Valenceimmoral -0.9474 0.0819 -11.57 < 0.0000000000000002"
## [13] "IdentityOther:Valenceimmoral 0.8230 0.1086 7.58 0.000000000000035"
## [14] "Valenceimmoral ***"
## [15] "IdentityOther:Valenceimmoral ***"
## [16] "---"
## [17] "Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1"
## [18] ""
```
.panel[.panel-name[model_full]
```r
#随机截距，随机斜率
mod_full <- lme4::glmer(data = df.match, #数据
formula = ACC ~ 1 + Identity * Valence + (1 + Identity * Valence|Sub), #模型
family = binomial) #因变量二项分布
##performance::model_performance(mod_full)
summary(mod_full) %>%
capture.output() %>% .[c(6:8,13:18,21:26,30:34)]
```
```
## [1] ""
## [2] " AIC BIC logLik deviance df.resid "
## [3] " 9355.9 9459.4 -4664.0 9327.9 11996 "
## [4] ""
## [5] "Random effects:"
## [6] " Groups Name Variance Std.Dev. Corr "
## [7] " Sub (Intercept) 0.971 0.986 "
## [8] " IdentityOther 1.770 1.331 -0.79 "
## [9] " Valenceimmoral 1.028 1.014 -0.75 0.75 "
## [10] ""
## [11] "Fixed effects:"
## [12] " Estimate Std. Error z value Pr(>|z|)"
## [13] "(Intercept) 2.772 0.178 15.53 < 0.0000000000000002"
## [14] "IdentityOther -0.870 0.234 -3.71 0.00021"
## [15] "Valenceimmoral -1.150 0.190 -6.06 0.0000000013"
## [16] "IdentityOther ***"
## [17] "Valenceimmoral ***"
## [18] "IdentityOther:Valenceimmoral ***"
## [19] "---"
## [20] "Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1"
```
.panel[.panel-name[模型比较anova]
```r
stats::anova(mod_null, mod, mod_full) #比较三个模型
```
```
## Data: df.match
## Models:
## mod: ACC ~ 1 + Identity * Valence + (1 | Sub)
## mod_null: ACC ~ (1 + Identity * Valence | Sub)
## mod_full: ACC ~ 1 + Identity * Valence + (1 + Identity * Valence | Sub)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## mod 5 9639 9676 -4814 9629
## mod_null 11 9379 9460 -4678 9357 272.1 6 < 0.0000000000000002 ***
## mod_full 14 9356 9459 -4664 9328 28.9 3 0.0000023 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```
.panel[.panel-name[模型比较compare_performance]
```r
performance::compare_performance(mod_null, mod, mod_full, rank = TRUE, verbose = FALSE)
```
]]]]]]
---
#10.3 代码实操
##结果解读
```r
summary(mod_full) %>% capture.output() %>% .[c(21:27)]
```
```
## [1] ""
## [2] "Fixed effects:"
## [3] " Estimate Std. Error z value Pr(>|z|)"
## [4] "(Intercept) 2.772 0.178 15.53 < 0.0000000000000002"
## [5] "IdentityOther -0.870 0.234 -3.71 0.00021"
## [6] "Valenceimmoral -1.150 0.190 -6.06 0.0000000013"
## [7] "IdentityOther:Valenceimmoral 0.988 0.272 3.63 0.00028"
```
<img src="./picture/chp10/logit2.png" width="60%" style="display: block; margin: auto;" />
.pull-left[
- MoralSelf:
`\(P=\frac{e^{2.77}}{1+e^{2.77}} = 0.941\)`
<br>
- ImmoralSelf:
`\(P=\frac{e^{2.77-1.15}}{1+e^{2.77-1.15}} = 0.835\)`
]
.pull-right[
- MoralOther:
`\(P=\frac{e^{2.77-0.87}}{1+e^{2.77-0.87}} = 0.870\)`
<br>
- ImmoralOther:
`\(P=\frac{e^{2.77-0.87-1.15+0.99}}{1+e^{2.77-0.87-1.15+0.99}} = 0.851\)`
]
---
#10.3 代码实操
```r
#交互作用
interactions::cat_plot(model = mod_full,
pred = Identity,
modx = Valence)
```
<img src="chapter_10_files/figure-html/unnamed-chunk-21-1.png" width="540" style="display: block; margin: auto;" />
---
#10.4 方法比较
.panelset[
.panel[.panel-name[anova]
```r
res <- bruceR::MANOVA(data = df.match.aov, #数据
subID = 'Sub', # 被试编号
dv= 'mean_ACC', # 因变量
within = c('Identity', 'Valence')) #自变量（被试内）
```
```r
capture.output(res) %>% .[3:8]
```
```
## [1] "Response: mean_ACC"
## [2] " Effect df MSE F ges p.value"
## [3] "1 Identity 1, 40 0.01 3.08 + .017 .087"
## [4] "2 Valence 1, 40 0.01 16.26 *** .068 <.001"
## [5] "3 Identity:Valence 1, 40 0.01 8.52 ** .038 .006"
## [6] "---"
```
.panel[.panel-name[EMMAMNS]
```r
res %>%
bruceR::EMMEANS(effect = 'Valence', by = 'Identity') %>%
capture.output()
```
```
## [1] "------ EMMEANS (effect = \"Valence\") ------"
## [2] ""
## [3] "Joint Tests of \"Valence\":"
## [4] "────────────────────────────────────────────────────────────────"
## [5] " Effect \"Identity\" df1 df2 F p η²p [90% CI of η²p]"
## [6] "────────────────────────────────────────────────────────────────"
## [7] " Valence Self 1 40 35.614 <.001 *** .471 [.282, .610]"
## [8] " Valence Other 1 40 0.412 .525 .010 [.000, .114]"
## [9] "────────────────────────────────────────────────────────────────"
## [10] "Note. Simple effects of repeated measures with 3 or more levels"
## [11] "are different from the results obtained with SPSS MANOVA syntax."
## [12] ""
## [13] "Estimated Marginal Means of \"Valence\":"
## [14] "───────────────────────────────────────────────────"
## [15] " \"Valence\" \"Identity\" Mean [95% CI of Mean] S.E."
## [16] "───────────────────────────────────────────────────"
## [17] " moral Self 0.916 [0.885, 0.947] (0.015)"
## [18] " immoral Self 0.814 [0.776, 0.852] (0.019)"
## [19] " moral Other 0.844 [0.809, 0.879] (0.017)"
## [20] " immoral Other 0.829 [0.794, 0.864] (0.017)"
## [21] "───────────────────────────────────────────────────"
## [22] ""
## [23] "Pairwise Comparisons of \"Valence\":"
## [24] "────────────────────────────────────────────────────────────────────────────────────────"
## [25] " Contrast \"Identity\" Estimate S.E. df t p Cohen’s d [95% CI of d]"
## [26] "────────────────────────────────────────────────────────────────────────────────────────"
## [27] " immoral - moral Self -0.102 (0.017) 40 -5.968 <.001 *** -0.736 [-0.985, -0.487]"
## [28] " immoral - moral Other -0.015 (0.024) 40 -0.642 .525 -0.111 [-0.459, 0.238]"
## [29] "────────────────────────────────────────────────────────────────────────────────────────"
## [30] "Pooled SD for computing Cohen’s d: 0.139"
## [31] "No need to adjust p values."
## [32] ""
## [33] "Disclaimer:"
## [34] "By default, pooled SD is Root Mean Square Error (RMSE)."
## [35] "There is much disagreement on how to compute Cohen’s d."
## [36] "You are completely responsible for setting `sd.pooled`."
## [37] "You might also use `effectsize::t_to_d()` to compute d."
## [38] ""
```
.panel[.panel-name[GLM]
```r
stats::anova(mod_full)
```
```
## Analysis of Variance Table
## npar Sum Sq Mean Sq F value
## Identity 1 0.2 0.2 0.2
## Valence 1 26.8 26.8 26.8
## Identity:Valence 1 13.6 13.6 13.6
```
.panel[.panel-name[HLM]
```r
mod_anova <- lme4::lmer(data = df.match,
formula = ACC ~ 1 + Identity * Valence + (1 + Identity * Valence|Sub))
stats::anova(mod_anova)
```
```
## Analysis of Variance Table
## npar Sum Sq Mean Sq F value
## Identity 1 0.42 0.42 3.71
## Valence 1 3.17 3.17 27.69
## Identity:Valence 1 0.98 0.98 8.54
```
.panel[.panel-name[HLM_mean]
```r
mod_mean <- lme4::lmer(data = df.match.aov,
formula = mean_ACC ~ 1 + Identity * Valence + (1|Sub) + (1|Identity:Sub) + (1|Valence:Sub))
stats::anova(mod_mean)
```
```
## Analysis of Variance Table
## npar Sum Sq Mean Sq F value
## Identity 1 0.0272 0.0272 3.08
## Valence 1 0.1410 0.1410 15.93
## Identity:Valence 1 0.0769 0.0769 8.69
```
.panel[.panel-name[compare]
```r
performance::compare_performance(mod_full, mod_anova, rank = TRUE, verbose = FALSE)
```
```r
stats::anova(mod_full, mod_anova)
```
```
## Data: df.match
## Models:
## mod_full: ACC ~ 1 + Identity * Valence + (1 + Identity * Valence | Sub)
## mod_anova: ACC ~ 1 + Identity * Valence + (1 + Identity * Valence | Sub)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## mod_full 14 9356 9459 -4664 9328
## mod_anova 15 8396 8507 -4183 8366 962 1 <0.0000000000000002 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```
]]]]]]]
---
#10.4 方法比较
## 留出法
```r
# 设置种子以确保结果的可重复性
set.seed(456)
# 随机选择70%的数据作为训练集，剩余的30%作为测试集
train_index <- caret::createDataPartition(df.match$Sub, p = 0.7, list = FALSE)
train_data <- df.match[train_index, ]
test_data <- df.match[-train_index, ]
# 根据训练集生成模型
model_full <- lme4::glmer(data = train_data,
formula = ACC ~ 1 + Identity * Valence + (1 + Identity * Valence|Sub),
family = binomial)
model_anova <- lme4::lmer(data = train_data,
formula = ACC ~ 1 + Identity * Valence + (1 + Identity * Valence|Sub))
# 使用模型进行预测
pre_mod_full <- stats::predict(model_full, newdata = test_data, type = 'response')
pre_mod_anova <- stats::predict(model_anova, newdata = test_data)
```
---
#10.4 方法比较
## 留出法
.pull-left[
```r
# 计算模型的性能指标
performance_mod_full <- c(RMSE = sqrt(mean((test_data$ACC - pre_mod_full)^2)),
R2 = cor(test_data$ACC, pre_mod_full)^2)
# 打印性能指标
print(performance_mod_full)
```
```
## RMSE R2
## 0.342402 0.074985
```
]
.pull-right[
```r
# 计算模型的性能指标
performance_mod_anova <- c(RMSE = sqrt(mean((test_data$ACC - pre_mod_anova)^2)),
R2 = cor(test_data$ACC, pre_mod_anova)^2)
# 打印性能指标
print(performance_mod_anova)
```
```
## RMSE R2
## 0.342263 0.075676
```
]
---
#10.4 方法比较
## 留出法
```r
# 将预测概率转换为分类结果
predicted_classes <- ifelse(pre_mod_full > 0.5, 1, 0)
# 计算混淆矩阵
confusion_matrix <- caret::confusionMatrix(as.factor(predicted_classes), as.factor(test_data$ACC))
```
---
#10.4 方法比较
## 留出法
```r
# 打印混淆矩阵和性能指标
print(confusion_matrix)
```
```
## Confusion Matrix and Statistics
##
## Reference
## Prediction 0 1
## 0 37 27
## 1 499 3037
##
## Accuracy : 0.854
## 95% CI : (0.842, 0.865)
## No Information Rate : 0.851
## P-Value [Acc > NIR] : 0.33
##
## Kappa : 0.095
##
## Mcnemar's Test P-Value : <0.0000000000000002
##
## Sensitivity : 0.0690
## Specificity : 0.9912
## Pos Pred Value : 0.5781
## Neg Pred Value : 0.8589
## Prevalence : 0.1489
## Detection Rate : 0.0103
## Detection Prevalence : 0.0178
## Balanced Accuracy : 0.5301
##
## 'Positive' Class : 0
##
```
---
#10.4 方法比较
## 留出法
.pull-left[
```r
# 计算ROC曲线和AUC
roc_result <- pROC::roc(test_data$ACC, pre_mod_full)
print(roc_result)
```
```
##
## Call:
## roc.default(response = test_data$ACC, predictor = pre_mod_full)
##
## Data: pre_mod_full in 536 controls (test_data$ACC 0) < 3064 cases (test_data$ACC 1).
## Area under the curve: 0.699
```
]
.pull-right[
```r
# 绘制ROC曲线
plot(roc_result, main = "ROC Curve", col = "blue", lwd = 2)
abline(a = 0, b = 1, lty = 2) # 添加对角线
```
<img src="chapter_10_files/figure-html/unnamed-chunk-35-1.png" width="540" style="display: block; margin: auto;" />
]
---
#10.4 方法比较
## 重复测量分析的不足
.bigfont[
- 会产生难以解释的结果
- 假设在10个回答中，正确回答8次，错误回答2次
- 此时95%CI为[0.52,1.08] ( = 0.8 ± 0.275)
- 方差不齐，不满足方差分析基本假设
`$$\mu = np$$`
`$$𝜎 = √(𝑛𝑝𝑞 )$$`
`$$𝜎_p^2 = \frac{p(1-p)}{n}$$`
]
Jaeger, T. F. (2008). Categorical data analysis: Away from ANOVAs (transformation or not) and towards logit mixed models. *Journal of Memory and Language, 59*(4), 434-446. doi:http://dx.doi.org/10.1016/j.jml.2007.11.007
---
#10.5 其他分布
##泊松分布(Poisson distribution)
.bigfont[
- 在固定时间间隔或空间区域内发生某种事件的次数的概率。
- 适用于事件以恒定平均速率独立发生的情况
- 例如电话呼叫、网站访问、机器故障等。
`$$P(X = k) = \frac{e^{-\lambda} \lambda^k}{k!}$$`
- λ:事件在给定时间或空间内的平均发生率（或平均数量）。
- k:可能的事件发生次数，可以是0, 1, 2, …
]
---
#10.5 其他分布
##泊松分布(Poisson distribution)
```r
set.seed(123) # 设置随机种子以获得可重复的结果
random_samples <- rpois(1000, lambda = 5)
hist(random_samples,col = 'white', border = 'black',)
```
<img src="chapter_10_files/figure-html/unnamed-chunk-36-1.png" width="540" style="display: block; margin: auto;" />
---
#10.5 其他分布
##泊松分布(Poisson distribution)
---
#10.5 其他分布
##伽马分布(Gamma distribution)
<br>
.bigfont[
- 伽马分布（Gamma Distribution）是统计学的一种连续概率函数，是概率统计中一种非常重要的分布。
- “指数分布”和“卡方分布”都是伽马分布的特例。
`$$f(x | \alpha, \beta) = \frac{\beta^\alpha x^{\alpha-1} e^{-\beta x}}{\Gamma(\alpha)}$$`
- α:形状参数（shape parameter），决定了分布的曲线形态，尤其是峰值的位置和曲线的尖峭程度。
- β:尺度参数（scale parameter），影响分布的宽度；当尺度参数增大时，分布会变得更宽且矮平；尺度参数减小时，分布会变得更窄且高耸。
]
---
#10.5 其他分布
##伽马分布(Gamma distribution)
<img src="./picture/chp10/gamma.webp" width="60%" style="display: block; margin: auto;" />
---
class: center, middle
.tit_font[
思考
]
<br>
<span style="font-size: 50px;">信号检测论是否可以用广义线性模型分析？</span> <br>
</textarea>
<style data-target="print-only">@media screen {.remark-slide-container{display:block;}.remark-slide-scaler{box-shadow:none;}}</style>
<script src="https://remarkjs.com/downloads/remark-latest.min.js"></script>
<script>var slideshow = remark.create({
"highlightLines": true,
"highlightStyle": "github",
"countIncrementalSlides": false,
"seal": true,
"ratio": "16:9"
});
if (window.HTMLWidgets) slideshow.on('afterShowSlide', function (slide) {
window.dispatchEvent(new Event('resize'));
});
(function(d) {
var s = d.createElement("style"), r = d.querySelector(".remark-slide-scaler");
if (!r) return;
s.type = "text/css"; s.innerHTML = "@page {size: " + r.style.width + " " + r.style.height +"; }";
d.head.appendChild(s);
})(document);
(function(d) {
var el = d.getElementsByClassName("remark-slides-area");
if (!el) return;
var slide, slides = slideshow.getSlides(), els = el[0].children;
for (var i = 1; i < slides.length; i++) {
slide = slides[i];
if (slide.properties.continued === "true" || slide.properties.count === "false") {
els[i - 1].className += ' has-continuation';
}
}
var s = d.createElement("style");
s.type = "text/css"; s.innerHTML = "@media print { .has-continuation { display: none; } }";
d.head.appendChild(s);
})(document);
// delete the temporary CSS (for displaying all slides initially) when the user
// starts to view slides
(function() {
var deleted = false;
slideshow.on('beforeShowSlide', function(slide) {
if (deleted) return;
var sheets = document.styleSheets, node;
for (var i = 0; i < sheets.length; i++) {
node = sheets[i].ownerNode;
if (node.dataset["target"] !== "print-only") continue;
node.parentNode.removeChild(node);
}
deleted = true;
});
})();
// add `data-at-shortcutkeys` attribute to <body> to resolve conflicts with JAWS
// screen reader (see PR #262)
(function(d) {
let res = {};
d.querySelectorAll('.remark-help-content table tr').forEach(tr => {
const t = tr.querySelector('td:nth-child(2)').innerText;
tr.querySelectorAll('td:first-child .key').forEach(key => {
const k = key.innerText;
if (/^[a-z]$/.test(k)) res[k] = t; // must be a single letter (key)
});
});
d.body.setAttribute('data-at-shortcutkeys', JSON.stringify(res));
})(document);
(function() {
"use strict"
// Replace <script> tags in slides area to make them executable
var scripts = document.querySelectorAll(
'.remark-slides-area .remark-slide-container script'
);
if (!scripts.length) return;
for (var i = 0; i < scripts.length; i++) {
var s = document.createElement('script');
var code = document.createTextNode(scripts[i].textContent);
s.appendChild(code);
var scriptAttrs = scripts[i].attributes;
for (var j = 0; j < scriptAttrs.length; j++) {
s.setAttribute(scriptAttrs[j].name, scriptAttrs[j].value);
}
scripts[i].parentElement.replaceChild(s, scripts[i]);
}
})();
(function() {
var links = document.getElementsByTagName('a');
for (var i = 0; i < links.length; i++) {
if (/^(https?:)?\/\//.test(links[i].getAttribute('href'))) {
links[i].target = '_blank';
}
}
})();
// adds .remark-code-has-line-highlighted class to <pre> parent elements
// of code chunks containing highlighted lines with class .remark-code-line-highlighted
(function(d) {
const hlines = d.querySelectorAll('.remark-code-line-highlighted');
const preParents = [];
const findPreParent = function(line, p = 0) {
if (p > 1) return null; // traverse up no further than grandparent
const el = line.parentElement;
return el.tagName === "PRE" ? el : findPreParent(el, ++p);
};
for (let line of hlines) {
let pre = findPreParent(line);
if (pre && !preParents.includes(pre)) preParents.push(pre);
}
preParents.forEach(p => p.classList.add("remark-code-has-line-highlighted"));
})(document);</script>
<script>
slideshow._releaseMath = function(el) {
var i, text, code, codes = el.getElementsByTagName('code');
for (i = 0; i < codes.length;) {
code = codes[i];
if (code.parentNode.tagName !== 'PRE' && code.childElementCount === 0) {
text = code.textContent;
if (/^\\\((.|\s)+\\\)$/.test(text) || /^\\\[(.|\s)+\\\]$/.test(text) ||
/^\$\$(.|\s)+\$\$$/.test(text) ||
/^\\begin\{([^}]+)\}(.|\s)+\\end\{[^}]+\}$/.test(text)) {
code.outerHTML = code.innerHTML; // remove <code></code>
continue;
}
}
i++;
}
};
slideshow._releaseMath(document);
</script>
<!-- dynamically load mathjax for compatibility with self-contained -->
<script>
(function () {
var script = document.createElement('script');
script.type = 'text/javascript';
script.src = 'https://mathjax.rstudio.com/latest/MathJax.js?config=TeX-MML-AM_CHTML';
if (location.protocol !== 'file:' && /^https?:/.test(script.src))
script.src = script.src.replace(/^https?:/, '');
document.getElementsByTagName('head')[0].appendChild(script);
})();
</script>
</body>
</html>