关于Rmarkdown的一些参考资料：
R Markdown: The Definitive Guide
1.导入excel的数据
1.1安装一些必用的包
> #读取excel文件的readxl包
> if (!requireNamespace("readxl", quietly = TRUE)) {
+ install.packages("readxl", repos = "https://cloud.r-project.org/")
+ }
> #整理数据集的整合包tidyverse
> if (!requireNamespace("tidyverse", quietly = TRUE)) {
+ install.packages("tidyverse")
+ }
> if (!requireNamespace("purrr", quietly = TRUE)) {
+ install.packages("purrr")
+ }
> #画图的
> if (!requireNamespace("ggcorrplot", quietly = TRUE)) {
+ install.packages("ggcorrplot")
+ }
> if (!requireNamespace("GGally", quietly = TRUE)) {
+ install.packages("GGally")
+ }
> if (!requireNamespace("tidyr", quietly = TRUE)) {
+ install.packages("tidyr")
+ }1.2实验室记忆测试的数据
| ID | 年龄 | 性别 | MRI_3T | 受教育时间 | image_copy | image_dr | OB_Num | OB1 | OB2 | OB3 | OB4 | OB5 | OB6 | OB_dr | L1 | L2 | L3 | L4 | L5 | L_ir | L_listB | L_Interfere | L_dr | Log_ir_R | Log_ir_N | Log_dr_R | Log_dr_N | W1 | W2 | W3 | W4 | W_ir_R | W_ir_N | W_dr_R | W_dr_N | specific1 | specific2 | specific3 | specific4 | specific_ir | specific_dr | abstract1 | abstract2 | abstract3 | abstract4 | abstract_ir | abstract_dr | HW_R | HF_R | HP_R | D_KN | D_CN | D_CP |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | 27 | 1 | 1T | 18 | 34 | 27.0 | 6 | 13 | 14 | 14 | 13 | 14 | 16 | 16 | 10 | 14 | 13 | 15 | 14 | 66 | 9 | 14 | 14 | 30 | 12 | 27 | 12 | 12 | 14 | 13 | 14 | 53 | 16 | 14 | 13 | 5 | 13 | 12 | 14 | 44 | 14 | 0 | 2 | 6 | 11 | 19 | 9 | 50 | 40 | 45 | 68 | 61 | 89.70588 |
| 2 | 30 | 1 | 1T | 20 | 34 | 28.5 | 2 | 13 | 16 | 16 | 16 | 16 | 16 | 16 | 8 | 12 | 13 | 15 | 15 | 63 | 7 | 15 | 15 | 44 | 17 | 42 | 17 | 10 | 14 | 14 | 14 | 52 | 16 | 14 | 14 | 10 | 13 | 13 | 14 | 50 | 14 | 4 | 11 | 12 | 14 | 41 | 13 | 48 | 37 | 45 | 60 | 51 | 85.00000 |
| 3 | 23 | 2 | 1T | 18 | 35 | 17.5 | 7 | 12 | 13 | 15 | 15 | 14 | 15 | 15 | 9 | 13 | 14 | 15 | 15 | 66 | 8 | 15 | 13 | 32 | 13 | 28 | 12 | 9 | 14 | 14 | 14 | 51 | 14 | 14 | 13 | 2 | 10 | 14 | 14 | 40 | 13 | 0 | 4 | 6 | 7 | 17 | 6 | 50 | 48 | 40 | 53 | 42 | 79.24528 |
| 4 | 29 | 2 | 2T | 18 | 31 | 16.0 | 6 | 11 | 13 | 14 | 14 | 14 | 16 | 15 | 8 | 11 | 14 | 15 | 15 | 63 | 8 | 15 | 15 | 36 | 14 | 30 | 13 | 13 | 14 | 14 | 14 | 55 | 18 | 14 | 13 | 9 | 12 | 14 | 14 | 49 | 14 | 2 | 4 | 4 | 7 | 17 | 6 | 49 | 45 | 44 | 57 | 40 | 70.17544 |
| 5 | 38 | 2 | 1T | 14 | 31 | 14.0 | 7 | 10 | 7 | 6 | 9 | 10 | 10 | 9 | 8 | 10 | 13 | 11 | 15 | 57 | 9 | 11 | 12 | 29 | 12 | 23 | 10 | 3 | 9 | 7 | 9 | 28 | 9 | 9 | 9 | 0 | 2 | 4 | 7 | 13 | 5 | 0 | 0 | 0 | 0 | 0 | 0 | 47 | 29 | 32 | 63 | 50 | 79.36508 |
| 6 | 31 | 1 | 3T | 16 | 36 | 33.0 | 3 | 15 | 15 | 16 | 16 | 16 | 16 | 16 | 10 | 13 | 13 | 14 | 14 | 64 | 10 | 14 | 14 | 32 | 13 | 27 | 12 | 11 | 14 | 14 | 14 | 53 | 16 | 14 | 14 | 4 | 14 | 14 | 14 | 46 | 14 | 0 | 0 | 2 | 6 | 8 | 4 | 49 | 47 | 49 | 61 | 45 | 73.77049 |
1.3海马体积的数据
| ID | 样本亚区 | 扫描质量 | B_DG/CA4_Volume | L_DG/CA4_Volume | R_DG/CA4_Volume | B_A_DG/CA4_Volume | B_P_DG/CA4_Volume | B_CA2/3_Volume | L_CA2/3_Volume | R_CA2/3_Volume | B_A_CA2/3_Volume | B_P_CA2/3_Volume | B_CA1_Volume | L_CA1_Volume | R_CA1_Volume | B_A_CA1_Volume | B_P_CA1_Volume | B_海马旁核_Volume | L_海马旁核_Volume | R_海马旁核_Volume | B_A_海马旁核_Volume | B_P_海马旁核_Volume | B_前/旁海马旁核_Volume | L_前/旁海马旁核_Volume | R_前/旁海马旁核_Volume | B_A_前/旁海马旁核_Volume | B_P_前/旁海马旁核_Volume | B_海马旁回_Volume | L_海马旁回_Volume | R_海马旁回_Volume | B_Total_Hippocampal_Volume | L_Total_Hippocampal_Volume | R_Total_Hippocampal_Volume | B_A_Total_Hippocampal_Volume | B_P_Total_Hippocampal_Volume |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | IAC | 1 | 1163.0 | 542.0 | 621.0 | 297.2 | 865.6 | 323.8 | 143.0 | 180.8 | 86.26 | 237.5 | 1233 | 570.4 | 662.6 | 385.6 | 847.3 | 1083.0 | 539.6 | 543.4 | 394.3 | 688.6 | 648.3 | 326.6 | 321.7 | 222.8 | 425.5 | 614.8 | 300.8 | 314.0 | 5065.9 | 2422.4 | 2643.5 | 2000.96 | 3064.5 |
| 2 | IAC | 1 | 1240.0 | 581.6 | 658.4 | 472.4 | 767.3 | 416.5 | 173.6 | 242.9 | 180.70 | 235.9 | 1118 | 561.5 | 556.5 | 480.8 | 636.8 | 1043.0 | 566.5 | 476.5 | 511.5 | 531.4 | 738.3 | 431.7 | 306.6 | 358.1 | 380.2 | 1016.0 | 511.5 | 504.5 | 5571.8 | 2826.4 | 2745.4 | 3019.50 | 2551.6 |
| 3 | IAC | 1 | 1097.0 | 482.3 | 614.8 | 294.5 | 802.7 | 408.8 | 176.3 | 232.5 | 109.70 | 299.1 | 1218 | 585.4 | 632.6 | 438.2 | 779.4 | 1190.0 | 599.4 | 590.8 | 435.5 | 754.7 | 698.8 | 370.7 | 328.1 | 268.1 | 430.6 | 751.4 | 389.7 | 361.7 | 5364.0 | 2603.8 | 2760.5 | 2297.40 | 3066.5 |
| 4 | MAD | 2 | 996.4 | 499.3 | 497.1 | 284.4 | 711.9 | 354.1 | 191.0 | 163.1 | 108.60 | 245.5 | 1087 | 562.6 | 524.4 | 367.4 | 719.3 | 1016.0 | 511.7 | 504.3 | 355.8 | 660.1 | 579.7 | 306.9 | 272.8 | 218.5 | 361.2 | 638.2 | 240.5 | 397.7 | 4671.4 | 2312.0 | 2359.4 | 1972.90 | 2698.0 |
| 5 | MAD | 1 | 1236.0 | 613.1 | 622.9 | 342.9 | 893.3 | 423.0 | 213.4 | 209.6 | 133.50 | 289.6 | 1198 | 549.6 | 648.4 | 397.7 | 800.1 | 978.6 | 492.4 | 486.2 | 410.7 | 567.9 | 677.6 | 316.6 | 361.0 | 311.7 | 365.9 | 839.7 | 380.3 | 459.4 | 5352.9 | 2565.4 | 2787.5 | 2436.20 | 2916.8 |
| 6 | IAC | 1 | 1339.0 | 650.0 | 689.0 | 348.2 | 990.9 | 371.9 | 169.9 | 202.0 | 94.54 | 277.4 | 1356 | 656.7 | 699.3 | 422.4 | 933.4 | 1207.0 | 605.6 | 601.4 | 475.0 | 732.3 | 633.8 | 322.7 | 311.1 | 256.6 | 377.2 | 846.9 | 396.9 | 450.0 | 5754.6 | 2801.8 | 2952.8 | 2443.64 | 3311.2 |
2.探索性数据分析
探索性数据分析(Exploratory data analysis,EDA)是在进行推断性统计建模之前，对数据的分布、变量之间的关系、观测之间的聚集等特性用汇总统计、作图等方法进行探索，这是必不可少的步骤。有时我们不知道数据告诉我们什么信息，不知道要提什么问题，这就必须用EDA来获得对数据的洞察，提示我们随后要进行的研究。即使一开始研究目标就很明确，我们也需要用EDA方法确定数据中是否有错误，变量分布是否满足模型需要的假定，变量值有没有异常，变量之间的关系从图形直观看是否符合我们要建立的模型。
数据中的变量大致可以分为属性变量（分类变量）和区间变量（数值型的变量）。属性变量又分为名义型的,没有大小次序，以及有序的（如学历等级）。区间变量包括取整数值的和取浮点值的变量，其中取整数值的变量当可取值较少时，也呈现出一定的分类变量特征。
内容参考：R语言教程-李东风
用result数据集为例子。首先需要了解数据框的行数、列数，有哪些变量，变量的存储类型，可以用 dplyr::glimpse():
Rows: 217
Columns: 54
$ ID <fct> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17,…
$ 年龄 <dbl> 27, 30, 23, 29, 38, 31, 30, 37, 32, 23, 31, 33, 36, 32, 27…
$ 性别 <fct> 1, 1, 2, 2, 2, 1, 1, 1, 2, 1, 2, 1, 1, 1, 1, 2, 2, 1, 1, 1…
$ MRI_3T <fct> 1T, 1T, 1T, 2T, 1T, 3T, 3T, 2T, 3T, 2T, 3T, 2T, 3T, 1T, 2T…
$ 受教育时间 <dbl> 18.0, 20.0, 18.0, 18.0, 14.0, 16.0, 18.0, 16.0, 17.0, 17.0…
$ image_copy <dbl> 34, 34, 35, 31, 31, 36, 35, 34, 35, 36, 34, 35, 36, 36, 34…
$ image_dr <dbl> 27.0, 28.5, 17.5, 16.0, 14.0, 33.0, 25.0, 29.0, 27.0, 17.5…
$ OB_Num <dbl> 6, 2, 7, 6, 7, 3, 4, 3, 7, 6, 3, 7, 6, 3, 7, 7, 4, 3, 7, 7…
$ OB1 <dbl> 13, 13, 12, 11, 10, 15, 12, 15, 13, 11, 12, 12, 11, 13, 14…
$ OB2 <dbl> 14, 16, 13, 13, 7, 15, 12, 15, 14, 10, 15, 14, 11, 14, 14,…
$ OB3 <dbl> 14, 16, 15, 14, 6, 16, 15, 16, 11, 10, 16, 13, 13, 16, 14,…
$ OB4 <dbl> 13, 16, 15, 14, 9, 16, 16, 16, 15, 8, 14, 13, 12, 15, 14, …
$ OB5 <dbl> 14, 16, 14, 14, 10, 16, 15, 15, 12, 14, 14, 15, 13, 16, 14…
$ OB6 <dbl> 16, 16, 15, 16, 10, 16, 16, 16, 14, 16, 14, 14, 16, 16, 14…
$ OB_dr <dbl> 16, 16, 15, 15, 9, 16, 16, 16, 12, 12, 15, 15, 16, 16, 14,…
$ L1 <dbl> 10, 8, 9, 8, 8, 10, 6, 4, 12, 7, 6, 10, 6, 11, 8, 6, 6, 6,…
$ L2 <dbl> 14, 12, 13, 11, 10, 13, 11, 6, 14, 9, 13, 14, 10, 13, 10, …
$ L3 <dbl> 13, 13, 14, 14, 13, 13, 12, 12, 15, 12, 14, 15, 14, 15, 11…
$ L4 <dbl> 15, 15, 15, 15, 11, 14, 13, 11, 15, 10, 15, 15, 15, 13, 12…
$ L5 <dbl> 14, 15, 15, 15, 15, 14, 13, 12, 15, 12, 15, 15, 13, 15, 15…
$ L_ir <dbl> 66, 63, 66, 63, 57, 64, 55, 45, 71, 50, 63, 69, 58, 67, 56…
$ L_listB <dbl> 9, 7, 8, 8, 9, 10, 5, 6, 9, 5, 7, 10, 7, 8, 4, 2, 6, 8, 7,…
$ L_Interfere <dbl> 14, 15, 15, 15, 11, 14, 11, 9, 15, 9, 14, 15, 14, 15, 14, …
$ L_dr <dbl> 14, 15, 13, 15, 12, 14, 13, 12, 15, 7, 14, 15, 13, 15, 14,…
$ Log_ir_R <dbl> 30, 44, 32, 36, 29, 32, 31, 30, 40, 25, 33, 33, 39, 30, 43…
$ Log_ir_N <dbl> 12, 17, 13, 14, 12, 13, 12, 12, 15, 10, 13, 13, 15, 12, 16…
$ Log_dr_R <dbl> 27, 42, 28, 30, 23, 27, 28, 31, 34, 21, 31, 33, 37, 29, 39…
$ Log_dr_N <dbl> 12, 17, 12, 13, 10, 12, 12, 13, 14, 9, 13, 14, 16, 13, 16,…
$ W1 <dbl> 12, 10, 9, 13, 3, 11, 8, 8, 14, 6, 9, 14, 13, 14, 7, 7, 12…
$ W2 <dbl> 14, 14, 14, 14, 9, 14, 14, 14, 14, 11, 13, 14, 14, 14, 14,…
$ W3 <dbl> 13, 14, 14, 14, 7, 14, 14, 14, 14, 13, 14, 14, 14, 14, 13,…
$ W4 <dbl> 14, 14, 14, 14, 9, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14,…
$ W_ir_R <dbl> 53, 52, 51, 55, 28, 53, 50, 50, 56, 44, 50, 56, 55, 56, 48…
$ W_ir_N <dbl> 16, 16, 14, 18, 9, 16, 15, 15, 19, 11, 15, 19, 18, 19, 13,…
$ W_dr_R <dbl> 14, 14, 14, 14, 9, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14,…
$ W_dr_N <dbl> 13, 14, 13, 13, 9, 14, 14, 14, 14, 13, 14, 14, 14, 14, 13,…
$ specific1 <dbl> 5, 10, 2, 9, 0, 4, 5, 2, 13, 3, 11, 12, 14, 12, 7, 2, 11, …
$ specific2 <dbl> 13, 13, 10, 12, 2, 14, 7, 10, 14, 9, 14, 14, 14, 14, 11, 8…
$ specific3 <dbl> 12, 13, 14, 14, 4, 14, 12, 14, 14, 11, 14, 14, 14, 14, 13,…
$ specific4 <dbl> 14, 14, 14, 14, 7, 14, 14, 14, 14, 12, 14, 14, 14, 14, 14,…
$ specific_ir <dbl> 44, 50, 40, 49, 13, 46, 38, 40, 55, 35, 53, 54, 56, 54, 45…
$ specific_dr <dbl> 14, 14, 13, 14, 5, 14, 14, 14, 14, 12, 14, 13, 14, 14, 14,…
$ abstract1 <dbl> 0, 4, 0, 2, 0, 0, 1, 1, 4, 0, 1, 5, 1, 0, 1, 1, 1, 0, 0, 3…
$ abstract2 <dbl> 2, 11, 4, 4, 0, 0, 3, 4, 12, 1, 2, 12, 4, 7, 4, 1, 6, 1, 4…
$ abstract3 <dbl> 6, 12, 6, 4, 0, 2, 3, 9, 13, 2, 4, 14, 8, 12, 9, 3, 8, 4, …
$ abstract4 <dbl> 11, 14, 7, 7, 0, 6, 7, 12, 12, 4, 7, 14, 12, 12, 10, 4, 9,…
$ abstract_ir <dbl> 19, 41, 17, 17, 0, 8, 14, 26, 41, 7, 14, 45, 25, 31, 24, 9…
$ abstract_dr <dbl> 9, 13, 6, 6, 0, 4, 8, 12, 12, 2, 4, 13, 10, 12, 10, 0, 9, …
$ HW_R <dbl> 50, 48, 50, 49, 47, 49, 50, 50, 49, 50, 50, 50, 50, 50, 50…
$ HF_R <dbl> 40, 37, 48, 45, 29, 47, 46, 42, 45, 34, 47, 46, 49, 44, 47…
$ HP_R <dbl> 45, 45, 40, 44, 32, 49, 46, 45, 48, 35, 48, 48, 47, 36, 47…
$ D_KN <dbl> 68, 60, 53, 57, 63, 61, 61, 72, 68, 57, 73, 70, 69, 64, 68…
$ D_CN <dbl> 61, 51, 42, 40, 50, 45, 50, 58, 63, 45, 71, 58, 66, 56, 61…
$ D_CP <dbl> 89.70588, 85.00000, 79.24528, 70.17544, 79.36508, 73.77049…2.1单个分类变量分布
|>和%>%应该都是pip管道符，具有相同的功能
筛选出result中因子型的变量：
[1] "ID" "性别" "MRI_3T"以 MRI_3T 变量为例，计算其频数和百分比如下：
# A tibble: 3 × 4
MRI_3T n pct `3`
<fct> <int> <dbl> <dbl>
1 1T 71 0.327 3
2 2T 71 0.327 3
3 3T 75 0.346 3在数据清理阶段对分类变量制作频数表， 可以发现输入错的类别值。 用基本R的table等函数统计：
1T 2T 3T
71 71 75 base::table()默认忽略缺失值,这在指标显示时会有误导。加选项useNA=“ifany”：
1T 2T 3T
71 71 75 百分比表：
1T 2T 3T
0.3271889 0.3271889 0.3456221 可以直接用ggplot2::geom_bar()制作频数条形图， 可以用geom_col()基于统计的频数表制作频数或者百分比的条形图，如：
> library(dplyr)
> library(ggplot2)
> library(knitr)
> MRI <- result %>%
+ count(MRI_3T) %>%
+ mutate(pct = n / sum(n))
> #频数
> ggplot(MRI, aes(x = MRI_3T, y = n)) +
+ geom_col() +
+ theme(axis.text.x = element_text(angle = 0, vjust = 1, hjust = 1)) +
+ labs(title = "扫描不同类型MRI的频数分布",x='MRI类型',y='人数')
> library(dplyr)
> library(ggplot2)
> library(scales)
> #百分比
> ggplot(MRI, aes(x = MRI_3T, y = pct)) +
+ geom_col(fill = "skyblue") + # 添加颜色填充柱状图
+ scale_y_continuous(labels = percent) +
+ theme_minimal() + # 使用简介的主题
+ theme(
+ axis.text.x = element_text(angle = 0, vjust = 1, hjust = 1)
+ ) +
+ labs(title = "扫描不同MRI的人比例分布",x='MRI类型',y = "百分比")
2.2单个连续变量分布
result中所有浮点型值的变量：
[1] "年龄" "受教育时间" "image_copy" "image_dr" "OB_Num"
[6] "OB1" "OB2" "OB3" "OB4" "OB5"
[11] "OB6" "OB_dr" "L1" "L2" "L3"
[16] "L4" "L5" "L_ir" "L_listB" "L_Interfere"
[21] "L_dr" "Log_ir_R" "Log_ir_N" "Log_dr_R" "Log_dr_N"
[26] "W1" "W2" "W3" "W4" "W_ir_R"
[31] "W_ir_N" "W_dr_R" "W_dr_N" "specific1" "specific2"
[36] "specific3" "specific4" "specific_ir" "specific_dr" "abstract1"
[41] "abstract2" "abstract3" "abstract4" "abstract_ir" "abstract_dr"
[46] "HW_R" "HF_R" "HP_R" "D_KN" "D_CN"
[51] "D_CP" 所有整数型的变量：
character(0)整数值的变量兼有区间型和分类变量的特点。 例如,年龄：
| 年龄 | n | pct |
|---|---|---|
| 20 | 14 | 0.0645161 |
| 21 | 15 | 0.0691244 |
| 22 | 9 | 0.0414747 |
| 23 | 12 | 0.0552995 |
| 24 | 10 | 0.0460829 |
| 25 | 6 | 0.0276498 |
| 26 | 7 | 0.0322581 |
| 27 | 10 | 0.0460829 |
| 28 | 18 | 0.0829493 |
| 29 | 8 | 0.0368664 |
| 30 | 10 | 0.0460829 |
| 31 | 22 | 0.1013825 |
| 32 | 13 | 0.0599078 |
| 33 | 14 | 0.0645161 |
| 34 | 6 | 0.0276498 |
| 35 | 12 | 0.0552995 |
| 36 | 8 | 0.0368664 |
| 37 | 10 | 0.0460829 |
| 38 | 6 | 0.0276498 |
| 39 | 3 | 0.0138249 |
| 40 | 3 | 0.0138249 |
| 41 | 1 | 0.0046083 |
> ggplot(age, aes(x = 年龄, y=pct)) +
+ geom_col() +
+ scale_y_continuous(
+ name = "Percent",
+ labels = scales::percent) +
+ labs(title = '年龄的频率分布',x='年龄',y='百分比')
变量OB_Num是指物体-关联测试中被试达到标准的尝试次数。可以用基本R的summary()函数显示“五数概括”和平均值：
Min. 1st Qu. Median Mean 3rd Qu. Max.
1.000 3.000 5.000 4.691 7.000 7.000 这说明没有缺失值，达到标准的最少次数为1，最多次数为7，中位数为5，平均次数为4.691，有四分之一的被试不到3次就达到标准，有四分之一的被试7次达到标准。 可以用dplyr::summarize()计算指定的统计量：
> result |>
+ summarize(
+ n = n(),
+ mean = round(mean(OB_Num, na.rm=TRUE),2),
+ std = round(sd(OB_Num, na.rm=TRUE),2),
+ Q3 = round(quantile(OB_Num,0.75,na.rm = TRUE),2),
+ max = round(max(OB_Num,na.rm=TRUE),2)
+ )# A tibble: 1 × 5
n mean std Q3 max
<int> <dbl> <dbl> <dbl> <dbl>
1 217 4.69 2.14 7 7用直方图反映分布情况：
> ggplot(result, aes(x = OB_Num)) +
+ geom_histogram(bins=50)+
+ labs(title = '物体-位置关联测试',x='达到标准的次数',y='人数')
可以用折线表示直方图：
> ggplot(result, aes(x = OB_Num)) +
+ geom_freqpoly(bins=50)+
+ labs(title = '物体-位置关联测试',x='达到标准的次数',y='人数')
2.3一个分类变量与一个连续变量之间的关系
例如，性别和image_copy之间的关系 这里只举点图的例子：
> library(ggplot2)
>
> ggplot(result, aes(x = 性别, y = image_copy)) +
+ geom_jitter(
+ #width控制水平抖动，height控制垂直抖动
+ width = 0.3, height = 0,
+ # size控制点的大小，alpha控制点的透明度
+ size = 0.8, alpha = 0.5
+ ) +
+ labs(title = '性别和复杂图形复制得分',
+ x = '性别',
+ y = '复制测试得分'
+ )
比较核密度估计：
> library(ggplot2)
> ggplot(result, aes(
+ x = image_copy,
+ color = 性别
+ )) +
+ # 添加透明度，以便重叠部分更容易观察
+ geom_density(alpha = 0.7) +
+ scale_color_manual(
+ # 为性别指定颜色
+ values = c("1" = "red", "2" = "blue"),
+ name = "性别", #图例的标题
+ labels = c("1" = "男", "2" = "女") #图例的标签
+ ) +
+ labs(title = '不同性别的复杂图形复制得分分布', x = '复杂图形复制得分', y = '密度')
2.4两个区间型变量之间的关系
考虑复杂图形中复制和延时回忆两种测量之间的关系，可以计算两者之间的相关系数：
image_copy image_dr
image_copy 1.0000000 0.3474822
image_dr 0.3474822 1.0000000用散点图表示两者关系： geom_point() 和 _geom_jitter()_都是画点图，区别在于前者不会避免点的重叠，而后者会通过随机便宜来减少重叠使每个点都可见。
> ggplot(result,aes(x=image_copy,y=image_dr))+
+ geom_point()+ #画点图
+ geom_smooth()+ #拟合曲线
+ labs(title = '复杂图形复制和延迟回忆之间的关系',
+ x='复制得分',y='延迟回忆') #修改标题和x轴y轴名称
多个区间变量之间也可以计算相关系数矩阵，比如我想看一下物体-位置关联测试的6个trail以及延迟回忆之间的关系：
> library(dplyr)
> d_OB<-result %>%
+ select(OB1,OB2,OB3,OB4,OB5,OB6,OB_dr) %>%
+ cor(use = "pairwise.complete.obs")
> d_OB OB1 OB2 OB3 OB4 OB5 OB6 OB_dr
OB1 1.0000000 0.6532678 0.5905180 0.4742066 0.4002097 0.4403817 0.4873231
OB2 0.6532678 1.0000000 0.6595229 0.5392743 0.5086593 0.5073519 0.4908057
OB3 0.5905180 0.6595229 1.0000000 0.6714603 0.5750609 0.5443031 0.5984939
OB4 0.4742066 0.5392743 0.6714603 1.0000000 0.6815659 0.5543423 0.5671765
OB5 0.4002097 0.5086593 0.5750609 0.6815659 1.0000000 0.6484405 0.5964245
OB6 0.4403817 0.5073519 0.5443031 0.5543423 0.6484405 1.0000000 0.6412322
OB_dr 0.4873231 0.4908057 0.5984939 0.5671765 0.5964245 0.6412322 1.0000000可以用图形的方式表示相关系数矩阵：
变量不太多时， 可以作散点图矩阵：
3.假设检验
前面部分讲了一部分描述性统计在R中的实现代码，还有一些未列出的可以参考：
下面仍然参考这本书，练习一下假设检验在R中的实现。
3.1单因素方差分析
单因素方差分析是用来比较多个组的均值是否相等的方法，可以看作是两样本t检验的扩展。 例如，我想看一下在沃灵顿识别记忆测试中单词、面孔、场景三种类型材料的记忆得分是否有显著差异，可以用单因素方差分析。虚无假设为\(H_0:\mu_1=\mu_2=\mu_3\)。
数据格式为（只显示前10行）：
| ID | word | face | place |
|---|---|---|---|
| 1 | 50 | 40 | 45 |
| 2 | 48 | 37 | 45 |
| 3 | 50 | 48 | 40 |
| 4 | 49 | 45 | 44 |
| 5 | 47 | 29 | 32 |
| 6 | 49 | 47 | 49 |
| 7 | 50 | 46 | 46 |
| 8 | 50 | 42 | 45 |
| 9 | 49 | 45 | 48 |
| 10 | 50 | 34 | 35 |
首先，我们需要长表格式，进行转换：
> Wtest <- result %>%
+ select(HW_R,HF_R,HP_R) %>% #从result中选出HW_R,HF_R,HP_R
+ #这三个变量分别命名为word，face，place
+ gather(key = "type", value = "score", HW_R, HF_R, HP_R)%>%
+ mutate(type = factor(type,
+ levels = c("HW_R", "HF_R", "HP_R"),
+ labels = c("word", "face", "place")))
> kable(head(Wtest,50),caption = "沃灵顿识别记忆测试数据")| type | score |
|---|---|
| word | 50 |
| word | 48 |
| word | 50 |
| word | 49 |
| word | 47 |
| word | 49 |
| word | 50 |
| word | 50 |
| word | 49 |
| word | 50 |
| word | 50 |
| word | 50 |
| word | 50 |
| word | 50 |
| word | 50 |
| word | 40 |
| word | 50 |
| word | 47 |
| word | 50 |
| word | 48 |
| word | 50 |
| word | 50 |
| word | 49 |
| word | 49 |
| word | 50 |
| word | 50 |
| word | 49 |
| word | 50 |
| word | 50 |
| word | 48 |
| word | 48 |
| word | 50 |
| word | 49 |
| word | 49 |
| word | 50 |
| word | 49 |
| word | 50 |
| word | 50 |
| word | 50 |
| word | 50 |
| word | 50 |
| word | 50 |
| word | 48 |
| word | 50 |
| word | 47 |
| word | 50 |
| word | 50 |
| word | 50 |
| word | 47 |
| word | 43 |
然后，我们可以用aov()函数进行单因素方差分析：
Df Sum Sq Mean Sq F value Pr(>F)
type 2 3623 1811.6 133.7 <2e-16 ***
Residuals 648 8780 13.5
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1输出结果显示，F检验的p值为0，所以我们拒绝虚无假设，即三种类型材料的记忆得分有显著差异。可以用并列盒形图比较各组的取值， 并检查分布对称性和方差齐性：
> # Plot boxplot
> library(ggplot2)
> ggplot(Wtest, aes(x = type, y = score)) +
+ geom_boxplot(fill = "skyblue") +
+ labs(title = "沃灵顿识别记忆测试得分箱线图", x = "材料类型", y = "记忆得分")
盒形图的结果显示，沃灵顿识别记忆测试的得分情况为：单词>面孔>场景。除了各组并列盒形图，也可以做点图并将各组的均值、均值置信区间作图，如：
> # Plot point
> Wtest_1 <- Wtest %>%
+ group_by(type) %>%
+ summarize(mean_score = mean(score),
+ sd_score = sd(score),
+ se_score = sd_score / sqrt(n()))
> ggplot(Wtest_1, aes(x = type, y = mean_score)) +
+ geom_point(size = 3) +
+ geom_errorbar(aes(ymin = mean_score - 1.96 * se_score,
+ ymax = mean_score + 1.96 * se_score),
+ width = 0.2) +
+ labs(title = "沃灵顿识别记忆测试得分点图", x = "材料类型", y = "记忆得分")
3.1.1事后检验
为了进一步确定三种材料之间的差异情况，可以通过TukeyHSD()函数进行事后检验。
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = score ~ type, data = Wtest)
$type
diff lwr upr p adj
face-word -4.4884793 -5.318625 -3.65833338 0.0000000
place-word -5.3963134 -6.226459 -4.56616748 0.0000000
place-face -0.9078341 -1.737980 -0.07768822 0.0280821事后检验的结果显示，两两之间差异显著。记忆效果为：单词>面孔>场景。
3.1.2方差不齐的情况
为了检查各组方差是否相等，可以进行Bartlett检验：
Bartlett test of homogeneity of variances
data: score by type
Bartlett's K-squared = 95.26, df = 2, p-value < 2.2e-16Bartlett检验的结果显示，p值为0.0001，所以我们拒绝虚无假设，即三种类型材料的记忆得分的方差不相等。 方差分析模型要求误差项独立同正态分布\(N(0,\sigma^2)\)，这意味着各组的因变量方差相等。实际中不同组的因变量可能有不同的方差。 R提供了oneway.test()函数作为独立两样本t检验的Welch方法的推广，可以不要求方差相等。
One-way analysis of means (not assuming equal variances)
data: score and type
F = 194.29, num df = 2.00, denom df = 391.49, p-value < 2.2e-16Welch’s ANOVA的结果显示，p值为0，所以我们拒绝虚无假设，即三种类型材料的记忆得分有显著差异。为了进行多重比较，可以进行两两t检验并不使用合并的标准差估计，使用Holm方法进行p值调整以控制总错误率：
> # Pairwise t-test
> with(Wtest, pairwise.t.test(score, type, pool.sd=FALSE, p.adjust.method="holm"))
Pairwise comparisons using t tests with non-pooled SD
data: score and type
word face
face <2e-16 -
place <2e-16 0.025
P value adjustment method: holm 在0.0001水平下，单词和面孔、场景均差异显著； 在0.05水平下，面孔和场景差异显著。 记忆效果为：单词>面孔>场景。
3.2协方差分析
在单因素方差分析问题中，有可能存在一个连续型的解释变量，对因变量（指标）有影响，希望在排除此连续型解释变量影响后对各组进行比较。 例如，我想看一下沃灵顿识别记忆测试中，不同记忆内容对记忆效果的影响，要排除年龄的作用。
> # Run ANCOVA
> Wtest <- result %>%
+ select(HW_R,HF_R,HP_R,年龄) %>%
+ gather(key = "type", value = "score", HW_R, HF_R, HP_R) %>%
+ mutate(type = factor(type,
+ levels = c("HW_R", "HF_R", "HP_R"),
+ labels = c("word", "face", "place")))
> model_2 <- aov(score ~ type + 年龄, data = Wtest)
> summary(model_2) Df Sum Sq Mean Sq F value Pr(>F)
type 2 3623 1811.6 133.640 <2e-16 ***
年龄 1 10 9.6 0.704 0.402
Residuals 647 8771 13.6
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1从’type’的p值可以看出，三种记忆内容的得分有显著差异，但是从’年龄’的p值可以看出，年龄对记忆得分没有显著影响。
3.3两因素方差分析
两因素方差分析，最简单的情形是有两个分组变量，每个分组变量有两个水平，这种情形称为二因素方差分析。
接下来我要先从result中提取所有与即时记忆、延迟回忆有关的变量，将其与海马的数据框合并成一个新的数据框hip_irdr。
> # Merge two data frames
> irdr<-result %>%
+ select(ID, OB1,OB2,OB3,OB4,OB5,OB6,OB_dr,
+ L_ir,L_dr,
+ Log_ir_N,Log_dr_N,
+ W_ir_N,W_dr_N,
+ specific_ir,specific_dr,
+ abstract_ir,abstract_dr) %>%
+ transmute(ID = ID,
+ 物体位置ir = (OB1 + OB2 + OB3 + OB4 + OB5 + OB6) / 6,
+ 物体位置dr = OB_dr,
+ 雷伊ir = L_ir,
+ 雷伊dr = L_dr,
+ 逻辑ir = Log_ir_N,
+ 逻辑dr = Log_dr_N,
+ 韦氏ir = W_ir_N,
+ 韦氏dr = W_dr_N,
+ 具体ir = specific_ir,
+ 具体dr = specific_dr,
+ 抽象ir = abstract_ir,
+ 抽象dr = abstract_dr)
> # Merge two data frames
> hip_irdr <- left_join(hip, irdr, by = "ID")
> kable(head(hip_irdr,20),caption = "海马体积和记忆测试数据")| ID | 样本亚区 | 扫描质量 | B_DG/CA4_Volume | L_DG/CA4_Volume | R_DG/CA4_Volume | B_A_DG/CA4_Volume | B_P_DG/CA4_Volume | B_CA2/3_Volume | L_CA2/3_Volume | R_CA2/3_Volume | B_A_CA2/3_Volume | B_P_CA2/3_Volume | B_CA1_Volume | L_CA1_Volume | R_CA1_Volume | B_A_CA1_Volume | B_P_CA1_Volume | B_海马旁核_Volume | L_海马旁核_Volume | R_海马旁核_Volume | B_A_海马旁核_Volume | B_P_海马旁核_Volume | B_前/旁海马旁核_Volume | L_前/旁海马旁核_Volume | R_前/旁海马旁核_Volume | B_A_前/旁海马旁核_Volume | B_P_前/旁海马旁核_Volume | B_海马旁回_Volume | L_海马旁回_Volume | R_海马旁回_Volume | B_Total_Hippocampal_Volume | L_Total_Hippocampal_Volume | R_Total_Hippocampal_Volume | B_A_Total_Hippocampal_Volume | B_P_Total_Hippocampal_Volume | 物体位置ir | 物体位置dr | 雷伊ir | 雷伊dr | 逻辑ir | 逻辑dr | 韦氏ir | 韦氏dr | 具体ir | 具体dr | 抽象ir | 抽象dr |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | IAC | 1 | 1163.000 | 542.000 | 621.00 | 297.2 | 865.6 | 323.80 | 143.000 | 180.800 | 86.26 | 237.5 | 1233.0 | 570.400 | 662.600 | 385.6 | 847.3 | 1083.00 | 539.600 | 543.400 | 394.3 | 688.6 | 648.300 | 326.600 | 321.7 | 222.8 | 425.5 | 614.800 | 300.800 | 314.000 | 5065.900 | 2422.400 | 2643.500 | 2000.960 | 3064.5 | 14.000000 | 16 | 66 | 14 | 12 | 12 | 16 | 13 | 44 | 14 | 19 | 9 |
| 2 | IAC | 1 | 1240.000 | 581.600 | 658.40 | 472.4 | 767.3 | 416.50 | 173.600 | 242.900 | 180.70 | 235.9 | 1118.0 | 561.500 | 556.500 | 480.8 | 636.8 | 1043.00 | 566.500 | 476.500 | 511.5 | 531.4 | 738.300 | 431.700 | 306.6 | 358.1 | 380.2 | 1016.000 | 511.500 | 504.500 | 5571.800 | 2826.400 | 2745.400 | 3019.500 | 2551.6 | 15.500000 | 16 | 63 | 15 | 17 | 17 | 16 | 14 | 50 | 14 | 41 | 13 |
| 3 | IAC | 1 | 1097.000 | 482.300 | 614.80 | 294.5 | 802.7 | 408.80 | 176.300 | 232.500 | 109.70 | 299.1 | 1218.0 | 585.400 | 632.600 | 438.2 | 779.4 | 1190.00 | 599.400 | 590.800 | 435.5 | 754.7 | 698.800 | 370.700 | 328.1 | 268.1 | 430.6 | 751.400 | 389.700 | 361.700 | 5364.000 | 2603.800 | 2760.500 | 2297.400 | 3066.5 | 14.000000 | 15 | 66 | 13 | 13 | 12 | 14 | 13 | 40 | 13 | 17 | 6 |
| 4 | MAD | 2 | 996.400 | 499.300 | 497.10 | 284.4 | 711.9 | 354.10 | 191.000 | 163.100 | 108.60 | 245.5 | 1087.0 | 562.600 | 524.400 | 367.4 | 719.3 | 1016.00 | 511.700 | 504.300 | 355.8 | 660.1 | 579.700 | 306.900 | 272.8 | 218.5 | 361.2 | 638.200 | 240.500 | 397.700 | 4671.400 | 2312.000 | 2359.400 | 1972.900 | 2698.0 | 13.666667 | 15 | 63 | 15 | 14 | 13 | 18 | 13 | 49 | 14 | 17 | 6 |
| 5 | MAD | 1 | 1236.000 | 613.100 | 622.90 | 342.9 | 893.3 | 423.00 | 213.400 | 209.600 | 133.50 | 289.6 | 1198.0 | 549.600 | 648.400 | 397.7 | 800.1 | 978.60 | 492.400 | 486.200 | 410.7 | 567.9 | 677.600 | 316.600 | 361.0 | 311.7 | 365.9 | 839.700 | 380.300 | 459.400 | 5352.900 | 2565.400 | 2787.500 | 2436.200 | 2916.8 | 8.666667 | 9 | 57 | 12 | 12 | 10 | 9 | 9 | 13 | 5 | 0 | 0 |
| 6 | IAC | 1 | 1339.000 | 650.000 | 689.00 | 348.2 | 990.9 | 371.90 | 169.900 | 202.000 | 94.54 | 277.4 | 1356.0 | 656.700 | 699.300 | 422.4 | 933.4 | 1207.00 | 605.600 | 601.400 | 475.0 | 732.3 | 633.800 | 322.700 | 311.1 | 256.6 | 377.2 | 846.900 | 396.900 | 450.000 | 5754.600 | 2801.800 | 2952.800 | 2443.640 | 3311.2 | 15.666667 | 16 | 64 | 14 | 13 | 12 | 16 | 14 | 46 | 14 | 8 | 4 |
| 7 | IAC | 1 | 1317.000 | 563.000 | 754.30 | 417.8 | 899.5 | 383.30 | 194.600 | 188.700 | 123.70 | 259.6 | 1286.0 | 704.900 | 580.900 | 457.6 | 828.2 | 1006.00 | 578.100 | 427.800 | 447.7 | 558.1 | 636.800 | 408.700 | 228.1 | 265.3 | 371.5 | 649.500 | 400.100 | 249.400 | 5278.600 | 2849.400 | 2429.200 | 2361.600 | 2916.9 | 14.333333 | 16 | 55 | 13 | 12 | 12 | 15 | 14 | 38 | 14 | 14 | 8 |
| 8 | MAD | 1 | 1084.000 | 556.900 | 527.50 | 372.0 | 712.3 | 348.40 | 184.700 | 163.700 | 135.50 | 212.9 | 1131.0 | 548.000 | 582.800 | 433.1 | 697.7 | 948.80 | 495.900 | 452.900 | 379.8 | 569.0 | 568.400 | 247.000 | 321.5 | 275.1 | 293.4 | 880.500 | 382.400 | 498.200 | 4961.100 | 2414.900 | 2546.600 | 2476.000 | 2485.3 | 15.500000 | 16 | 45 | 12 | 12 | 13 | 15 | 14 | 40 | 14 | 26 | 12 |
| 9 | MAD | 1 | 1222.000 | 602.100 | 620.40 | 428.9 | 793.6 | 368.00 | 173.900 | 194.100 | 151.60 | 216.3 | 1079.0 | 516.800 | 562.600 | 507.4 | 572.0 | 1340.00 | 665.400 | 674.100 | 564.9 | 774.6 | 658.400 | 324.400 | 333.9 | 325.4 | 333.0 | 945.900 | 517.800 | 428.100 | 5613.300 | 2800.400 | 2813.200 | 2924.100 | 2689.5 | 13.166667 | 12 | 71 | 15 | 15 | 14 | 19 | 14 | 55 | 14 | 41 | 12 |
| 10 | MAD | 1 | 961.507 | 463.867 | 497.64 | 228.4 | 733.1 | 235.46 | 105.116 | 130.344 | 52.22 | 183.2 | 1177.3 | 576.036 | 601.264 | 262.2 | 915.1 | 1179.33 | 615.641 | 563.689 | 451.5 | 727.8 | 495.334 | 279.134 | 216.2 | 145.3 | 350.1 | 713.976 | 291.341 | 422.635 | 4762.907 | 2331.135 | 2431.772 | 1853.596 | 2909.3 | 11.500000 | 12 | 50 | 7 | 10 | 9 | 11 | 13 | 35 | 12 | 7 | 2 |
| 11 | MAD | 1 | 1281.000 | 661.900 | 619.10 | 441.8 | 839.4 | 289.70 | 138.100 | 151.600 | 111.80 | 178.0 | 1249.0 | 618.200 | 630.800 | 395.8 | 852.9 | 1049.00 | 529.400 | 519.600 | 394.7 | 654.6 | 523.700 | 257.200 | 266.5 | 186.2 | 337.5 | 791.000 | 353.900 | 437.100 | 5183.400 | 2558.700 | 2624.700 | 2321.300 | 2862.4 | 14.166667 | 15 | 63 | 14 | 13 | 13 | 15 | 14 | 53 | 14 | 14 | 4 |
| 12 | IAC | 1 | 1206.000 | 572.000 | 634.00 | 346.0 | 859.6 | 337.20 | 174.000 | 163.200 | 108.60 | 228.5 | 1407.0 | 657.600 | 749.400 | 570.1 | 837.4 | 945.50 | 541.800 | 403.700 | 358.5 | 587.0 | 598.100 | 301.000 | 297.1 | 344.0 | 254.2 | 769.700 | 324.800 | 444.900 | 5263.500 | 2571.200 | 2692.300 | 2496.900 | 2766.7 | 13.500000 | 15 | 69 | 15 | 13 | 14 | 19 | 14 | 54 | 13 | 45 | 13 |
| 13 | MAD | 1 | 967.700 | 474.600 | 493.10 | 196.1 | 771.6 | 240.30 | 118.100 | 122.200 | 49.37 | 191.0 | 1083.0 | 559.400 | 523.600 | 355.2 | 727.7 | 981.30 | 478.800 | 502.500 | 299.5 | 681.8 | 421.100 | 212.500 | 208.6 | 179.3 | 241.8 | 666.000 | 318.000 | 348.000 | 4359.400 | 2161.400 | 2198.000 | 1745.470 | 2613.9 | 12.666667 | 16 | 58 | 13 | 15 | 16 | 18 | 14 | 56 | 14 | 25 | 10 |
| 14 | MAD | 1 | 1007.000 | 478.700 | 528.00 | 318.1 | 688.6 | 237.60 | 125.100 | 112.600 | 84.23 | 153.4 | 1297.0 | 630.400 | 667.000 | 502.8 | 794.7 | 902.90 | 439.500 | 463.500 | 391.6 | 511.3 | 418.600 | 209.000 | 209.6 | 201.4 | 217.1 | 699.300 | 310.700 | 388.600 | 4562.400 | 2193.400 | 2369.300 | 2197.430 | 2365.1 | 15.000000 | 16 | 67 | 15 | 12 | 13 | 19 | 14 | 54 | 14 | 31 | 12 |
| 15 | MAD | 1 | 1045.000 | 505.600 | 539.40 | 259.1 | 785.9 | 277.80 | 128.900 | 148.900 | 100.50 | 177.3 | 1313.0 | 658.200 | 654.800 | 545.4 | 767.7 | 1058.00 | 522.300 | 535.700 | 408.0 | 649.0 | 506.200 | 228.900 | 277.3 | 238.6 | 267.6 | 1244.000 | 569.400 | 674.600 | 5444.000 | 2613.300 | 2830.700 | 2795.600 | 2647.5 | 14.000000 | 14 | 56 | 14 | 16 | 16 | 13 | 13 | 45 | 14 | 24 | 10 |
| 16 | IAC | 1 | 1010.000 | 448.500 | 561.50 | 285.9 | 724.4 | 283.20 | 130.100 | 153.100 | 94.67 | 188.5 | 1130.0 | 546.100 | 583.900 | 422.9 | 707.3 | 963.10 | 497.900 | 465.200 | 355.6 | 607.5 | 440.300 | 209.800 | 230.5 | 172.4 | 267.9 | 648.700 | 346.900 | 301.800 | 4475.300 | 2179.300 | 2296.000 | 1980.170 | 2495.6 | 12.166667 | 11 | 48 | 11 | 14 | 12 | 7 | 5 | 33 | 11 | 9 | 0 |
| 18 | IAC | 1 | 1060.000 | 559.100 | 500.90 | 224.7 | 835.1 | 312.90 | 159.350 | 153.550 | 75.41 | 237.5 | 1139.0 | 544.700 | 594.300 | 301.0 | 838.4 | 791.20 | 430.600 | 360.600 | 338.4 | 452.7 | 408.300 | 221.400 | 186.9 | 169.7 | 238.6 | 742.600 | 349.500 | 393.100 | 4454.000 | 2264.650 | 2189.350 | 1851.810 | 2602.3 | 15.500000 | 16 | 56 | 14 | 12 | 11 | 14 | 13 | 39 | 13 | 12 | 4 |
| 19 | IAC | 1 | 964.900 | 486.800 | 478.10 | 355.6 | 609.3 | 315.90 | 180.800 | 135.100 | 116.10 | 199.8 | 1149.0 | 614.800 | 534.200 | 396.3 | 752.2 | 1152.00 | 594.600 | 557.400 | 552.4 | 599.6 | 614.300 | 328.000 | 286.3 | 293.2 | 321.0 | 1005.000 | 450.700 | 554.300 | 5201.100 | 2655.700 | 2545.400 | 2718.600 | 2481.9 | 13.333333 | 16 | 54 | 11 | 11 | 8 | 16 | 13 | 45 | 13 | 14 | 8 |
| 20 | MAD | 3 | 1082.000 | 533.400 | 548.60 | 372.0 | 710.0 | 242.20 | 130.300 | 111.900 | 81.24 | 161.0 | 1304.0 | 607.900 | 696.100 | 468.7 | 834.8 | 884.10 | 452.300 | 431.800 | 400.9 | 483.1 | 432.700 | 192.600 | 240.1 | 214.4 | 218.2 | 861.300 | 366.300 | 495.000 | 4806.300 | 2282.800 | 2523.500 | 2398.540 | 2407.1 | 11.500000 | 13 | 49 | 12 | 13 | 14 | 14 | 12 | 34 | 11 | 21 | 5 |
| 21 | IAC | 1 | 1190.000 | 609.000 | 581.10 | 402.3 | 787.8 | 317.10 | 161.000 | 156.100 | 105.00 | 212.1 | 913.9 | 452.500 | 462.800 | 237.9 | 676.0 | 1039.00 | 529.800 | 509.700 | 430.4 | 609.1 | 444.500 | 203.900 | 240.6 | 180.5 | 263.9 | 526.000 | 227.300 | 298.700 | 4430.500 | 2183.500 | 2249.000 | 1882.100 | 2548.9 | 14.333333 | 14 | 71 | 15 | 16 | 16 | 16 | 14 | 54 | 14 | 42 | 14 |
下面的例子是为了练习两个分类变量方差分析的代码，先忽略这个分析是否有意义。
我想看一下在具体语言配对联想测试中，被试即时回忆和延迟回忆的得分，是否受到样本亚区的影响。回忆时间（即时、延时）是被试内变量，样本亚区（IAC、MAD）是被试间变量，所以采用两因素重复测量方差分析。首先，检查数据格式：
> # Check data format
> special<-hip_irdr |>
+ select(ID,具体ir,具体dr,样本亚区) %>%
+ transmute(ID = ID,
+ 即时回忆 = 具体ir,
+ 延时回忆 = 具体dr,
+ example = 样本亚区)
>
> kable(head(special,20),caption = "具体语言配对联想测试数据")| ID | 即时回忆 | 延时回忆 | example |
|---|---|---|---|
| 1 | 44 | 14 | IAC |
| 2 | 50 | 14 | IAC |
| 3 | 40 | 13 | IAC |
| 4 | 49 | 14 | MAD |
| 5 | 13 | 5 | MAD |
| 6 | 46 | 14 | IAC |
| 7 | 38 | 14 | IAC |
| 8 | 40 | 14 | MAD |
| 9 | 55 | 14 | MAD |
| 10 | 35 | 12 | MAD |
| 11 | 53 | 14 | MAD |
| 12 | 54 | 13 | IAC |
| 13 | 56 | 14 | MAD |
| 14 | 54 | 14 | MAD |
| 15 | 45 | 14 | MAD |
| 16 | 33 | 11 | IAC |
| 18 | 39 | 13 | IAC |
| 19 | 45 | 13 | IAC |
| 20 | 34 | 11 | MAD |
| 21 | 54 | 14 | IAC |
Rows: 201
Columns: 4
$ ID <fct> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 18, 19…
$ 即时回忆 <dbl> 44, 50, 40, 49, 13, 46, 38, 40, 55, 35, 53, 54, 56, 54, 45, 3…
$ 延时回忆 <dbl> 14, 14, 13, 14, 5, 14, 14, 14, 14, 12, 14, 13, 14, 14, 14, 11…
$ example <fct> IAC, IAC, IAC, MAD, MAD, IAC, IAC, MAD, MAD, MAD, MAD, IAC, M…首先，将数据转换为长表格式：
> # Convert to long format
> special <- special %>%
+ gather(key = "time", value = "score", 即时回忆, 延时回忆)
> kable(head(special,20),caption = "具体语言配对联想测试数据")| ID | example | time | score |
|---|---|---|---|
| 1 | IAC | 即时回忆 | 44 |
| 2 | IAC | 即时回忆 | 50 |
| 3 | IAC | 即时回忆 | 40 |
| 4 | MAD | 即时回忆 | 49 |
| 5 | MAD | 即时回忆 | 13 |
| 6 | IAC | 即时回忆 | 46 |
| 7 | IAC | 即时回忆 | 38 |
| 8 | MAD | 即时回忆 | 40 |
| 9 | MAD | 即时回忆 | 55 |
| 10 | MAD | 即时回忆 | 35 |
| 11 | MAD | 即时回忆 | 53 |
| 12 | IAC | 即时回忆 | 54 |
| 13 | MAD | 即时回忆 | 56 |
| 14 | MAD | 即时回忆 | 54 |
| 15 | MAD | 即时回忆 | 45 |
| 16 | IAC | 即时回忆 | 33 |
| 18 | IAC | 即时回忆 | 39 |
| 19 | IAC | 即时回忆 | 45 |
| 20 | MAD | 即时回忆 | 34 |
| 21 | IAC | 即时回忆 | 54 |
检查数据格式：
Rows: 402
Columns: 4
$ ID <fct> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 18, 19,…
$ example <fct> IAC, IAC, IAC, MAD, MAD, IAC, IAC, MAD, MAD, MAD, MAD, IAC, MA…
$ time <chr> "即时回忆", "即时回忆", "即时回忆", "即时回忆", "即时回忆", "…
$ score <dbl> 44, 50, 40, 49, 13, 46, 38, 40, 55, 35, 53, 54, 56, 54, 45, 33…转换数据格式，time改为因子型变量：
> special <- special %>%
+ mutate(time = factor(time, levels = c("即时回忆", "延时回忆")))
> dplyr::glimpse(special)Rows: 402
Columns: 4
$ ID <fct> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 18, 19,…
$ example <fct> IAC, IAC, IAC, MAD, MAD, IAC, IAC, MAD, MAD, MAD, MAD, IAC, MA…
$ time <fct> 即时回忆, 即时回忆, 即时回忆, 即时回忆, 即时回忆, 即时回忆, 即…
$ score <dbl> 44, 50, 40, 49, 13, 46, 38, 40, 55, 35, 53, 54, 56, 54, 45, 33…然后，进行两因素重复测量方差分析：
> # Run two-way repeated measures ANOVA
> if (!requireNamespace("afex", quietly = TRUE)) {
+ install.packages("afex")
+ }
> library(afex)
> model_3 <- aov_car(score ~ time * example + Error(ID/time), data = special)
> summary(model_3)
Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
Sum Sq num Df Error SS den Df F value Pr(>F)
(Intercept) 281908 1 15727.1 199 3567.0725 <2e-16 ***
example 68 1 15727.1 199 0.8565 0.3558
time 78152 1 8908.1 199 1745.8464 <2e-16 ***
example:time 26 1 8908.1 199 0.5721 0.4503
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1结果表明，只有时间的主效应显著，说明被试在即时回忆和延迟回忆的得分有显著差异。结合柱状图来看，即时回忆的得分显著高于延迟回忆的得分：
> # Plot bar
> library(ggplot2)
> ggplot(special, aes(x = time, y = score, fill = example)) +
+ geom_bar(stat = "identity", position = "dodge") +
+ labs(title = "具体语言配对联想测试得分柱状图", x = "回忆时间", y = "得分")
4.扩展： 统计检验与线性模型的等价性
内容来源：蒋超课题组
在R中，统计检验和线性模型是等价的。
单样本 t 检验
例如，我们想检验复杂图形复制任务的分数分布是否明显偏离正态分布，虚无假设为\(H_0:\mu=0\)，备择假设为\(H_1:\mu\neq0\)。
One Sample t-test
data: result$image_copy
t = 270.44, df = 216, p-value < 2.2e-16
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
33.90088 34.39866
sample estimates:
mean of x
34.14977 输出结果显示，t检验的p值为0，所以我们拒绝虚无假设，即复杂图形复制任务的分数分布明显偏离正态分布。通过密度图可以看出，复杂图形复制任务的分数分布不符合正态分布。
> # Plot density
> library(ggplot2)
> ggplot(result, aes(x = image_copy)) +
+ geom_density(fill = "skyblue") +
+ labs(title = "复杂图形复制密度分布图", x = "copy", y = "Density")
先不管做这个假设有没有意义，我们用线性回归的方法做一遍，看看结果是否一致：
Call:
lm(formula = image_copy ~ 1, data = result)
Residuals:
Min 1Q Median 3Q Max
-7.6498 -1.1498 0.8502 1.8502 1.8502
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 34.1498 0.1263 270.4 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.86 on 216 degrees of freedom输出结果显示，线性回归的p值为0.所以我们拒绝虚无假设，即复杂图形复制任务的分数分布明显偏离正态分布。这与t检验的结果一致。为了方便比较两个模型，我们通过broom宏包将两个结果规整在一起：
> # Compare t-test and linear regression
> library(broom)
> library(knitr)
> #tidy() 得到我们通常用来报告结果的模型输出
>
> model_1_t_tidy <- tidy(model_1_t) %>%
+ mutate(model = "t.test()")
> model_1_lm_tidy <- tidy(model_1_lm) %>%
+ mutate(model = "lm(image_copy ~ 1)")
>
> results <- bind_rows(model_1_t_tidy, model_1_lm_tidy) %>%
+ select(model, estimate, statistic, p.value)
>
> kable(results)| model | estimate | statistic | p.value |
|---|---|---|---|
| t.test() | 34.14977 | 270.4413 | 0 |
| lm(image_copy ~ 1) | 34.14977 | 270.4413 | 0 |
上面例子，我们用不同的方法完成了相同的检验。事实上，在R语言t.test()内部会直接调用lm()函数，其函数语句和我们的这里代码也是一样的。
配对样本 t 检验
我想看一下被试在复杂图形任务中进行复制和延迟回忆的表现是否有显著差异，可以用配对样本t检验。虚无假设为\(H_0:\mu1=\mu2\)，备择假设为\(H_1:\mu1\neq\mu2\)。
> # Run paired t-test
> library(broom)
> library(dplyr)
> model_2_t <- t.test(result$image_copy, result$image_dr, paired = TRUE)
> model_2_t_tidy <- tidy(model_2_t) %>%
+ mutate(model = "t.test(x,y,paired = TRUE)")
> model_2_t_tidy# A tibble: 1 × 9
estimate statistic p.value parameter conf.low conf.high method alternative
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <chr> <chr>
1 11.9 32.7 1.63e-85 216 11.2 12.6 Paired t… two.sided
# ℹ 1 more variable: model <chr>写成线性模型为lm(image_copy - image_dr ~ 1) 先让两个变量相减，然后用线性模型检验：
> # Run linear regression
> model_2_lm <- lm(image_copy - image_dr ~ 1, data = result)
> model_2_lm_tidy <- tidy(model_2_lm) %>%
+ mutate(model = "lm(copy - delay ~ 1)")
> summary(model_2_lm)
Call:
lm(formula = image_copy - image_dr ~ 1, data = result)
Residuals:
Min 1Q Median 3Q Max
-11.871 -3.871 -0.371 4.629 13.629
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 11.8710 0.3634 32.66 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 5.354 on 216 degrees of freedom把两个模型的结果放在一起比较：
> # Compare t-test and linear regression
> results <- bind_rows(model_2_t_tidy, model_2_lm_tidy) %>%
+ select(model, estimate, statistic, p.value)
> kable(results)| model | estimate | statistic | p.value |
|---|---|---|---|
| t.test(x,y,paired = TRUE) | 11.87097 | 32.66246 | 0 |
| lm(copy - delay ~ 1) | 11.87097 | 32.66246 | 0 |
通过这个等价的线性模型，我们发现配对样本t检验和线性回归的结果是一致的，都拒绝了虚无假设，说明被试在复杂图形复制任务和延迟回忆任务中的表现有显著差异。
并且呢，配对样本t检验的工作方式似乎是通过计算两个变量的差值，得到新的一列然后用单样本t检验，这和我们的线性模型是一样的。
> # run paired t-test testing depression from g1 against g2
>
> model_2_t2 <- t.test(result$image_copy - result$image_dr)
> model_2_t2_tidy <- tidy(model_2_t2) %>% mutate(model = "t.test(x - y)")
> model_2_t2_tidy# A tibble: 1 × 9
estimate statistic p.value parameter conf.low conf.high method alternative
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <chr> <chr>
1 11.9 32.7 1.63e-85 216 11.2 12.6 One Samp… two.sided
# ℹ 1 more variable: model <chr>既然是原理这样，那么在回归模型之前，可以在数据框中用image_copy - image_dr构建新的一列，然后放进回归模型：
> result <- result %>%
+ mutate(
+ dif_image = image_copy - image_dr
+ )
>
> model_2_lm2 <- lm(dif_image ~ 1, data = result)
> model_2_lm2_tidy <- tidy(model_2_lm2) %>% mutate(model = "lm(delta ~ 1)")最后，把四个模型放在一起比较，我们看到得到了相同的结果：
> # Compare t-test and linear regression
> results <- bind_rows(model_2_t_tidy, model_2_lm_tidy,
+ model_2_t2_tidy, model_2_lm2_tidy) %>%
+ select(model, estimate, statistic, p.value)
>
> kable(results)| model | estimate | statistic | p.value |
|---|---|---|---|
| t.test(x,y,paired = TRUE) | 11.87097 | 32.66246 | 0 |
| lm(copy - delay ~ 1) | 11.87097 | 32.66246 | 0 |
| t.test(x - y) | 11.87097 | 32.66246 | 0 |
| lm(delta ~ 1) | 11.87097 | 32.66246 | 0 |
相关性检验
两个变量的关联检验，用的也比较常见。比如两个连续型变量，我们可能想知道它们之间的关联，以及这种关联是否显著。 比如，我们想看一下逻辑记忆测试的即时回忆和延时回忆是否相关，可以用Pearson相关系数检验。虚无假设为\(H_0:\rho=0\)，备择假设为\(H_1:\rho\neq0\)。
> # Run correlation test
> model_3_cor <- cor.test(result$Log_ir_N, result$Log_dr_N,
+ method = "pearson")
> model_3_cor_tidy <- tidy(model_3_cor) %>%
+ mutate(model = "cor.test(x, y)")
> model_3_cor
Pearson's product-moment correlation
data: result$Log_ir_N and result$Log_dr_N
t = 23.876, df = 215, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.8109878 0.8848901
sample estimates:
cor
0.8521326 写成线性模型为lm(Log_ir_N ~ Log_dr_N):
> # Run linear regression
> model_3_lm <- lm(Log_ir_N ~ Log_dr_N, data = result)
> model_3_lm_tidy <- tidy(model_3_lm) %>%
+ mutate(model = "lm(Log_ir_N ~ Log_dr_N)")
> summary(model_3_lm)
Call:
lm(formula = Log_ir_N ~ Log_dr_N, data = result)
Residuals:
Min 1Q Median 3Q Max
-3.1655 -0.5942 0.0487 0.7272 5.8345
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.41542 0.36507 12.10 <2e-16 ***
Log_dr_N 0.67858 0.02842 23.88 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.095 on 215 degrees of freedom
Multiple R-squared: 0.7261, Adjusted R-squared: 0.7249
F-statistic: 570 on 1 and 215 DF, p-value: < 2.2e-16把两个模型的结果放在一起比较：
> # Compare correlation test and linear regression
>
> results <- bind_rows(model_3_cor_tidy, model_3_lm_tidy) %>%
+ select(model, term, estimate, statistic, p.value)
> kable(results) | model | term | estimate | statistic | p.value |
|---|---|---|---|---|
| cor.test(x, y) | NA | 0.8521326 | 23.87559 | 0 |
| lm(Log_ir_N ~ Log_dr_N) | (Intercept) | 4.4154159 | 12.09463 | 0 |
| lm(Log_ir_N ~ Log_dr_N) | Log_dr_N | 0.6785836 | 23.87559 | 0 |
这次的表格不一样。那是因为，线性模型会不仅输出系数，而且还输出了模型的截距，因此两个模型的系数会不一样，但在 t-statistic 和 p.value 是一样的。
后续
后续内容直接看原文吧，这里只是练习一下Rmarkdown的使用。