本文为个人编写讲义，编写过程参考了推荐材料中的内容。如有相关问题，或者行文和理解有所纰漏，烦请邮箱告知以便订正，万分感谢！ 邮箱:

1 引言

1.1 R语言简介

R语言是由统计学家设计的专为数据处理和统计分析服务的语言。随着更新迭代，R语言的功能也越来越强大。R语言可在官网（https://www.r-project.org/）进行下载安装。此外，RStudio是十分有效且强大的R语言IDE，它可以帮助我们快速方便地进行相关脚本的编写和运行，可以在Posit官网（https://posit.co/）下载更新。

1.2 推荐材料

2 数据处理

2.1 包(package)的介绍

2.1.1 tidyverse

tidyverse（https://www.tidyverse.org/）是由Hadley Wickham及其团队开发的集合包，囊括了数据读取（readr）、数据整理（dplyr, tidyr）、数据可视化（ggplot2）、以及其他系列提升R使用体验的包（purrr, tibble, stringr, forcats）。我们通过install.packages()进行安装并加载。

library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.4.4     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.0
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

2.1.2 bruceR

bruceR（https://github.com/psychbruce/bruceR）是由包寒吴霜博士开发的面向社会科学背景人员的R语言综合包。它主要包括数据处理（data.table, dplyr, tidyr, stringr, ggplot2）和统计分析（emmeans, lmerTest, effectsize, performance, interactions）两方面内容，并且改写了其中部分函数，从而使得数据处理和统计分析变得更为便捷。同样的，我们可以使用install.packages()进行安装并加载。中文帮助文档可以参阅https://zhuanlan.zhihu.com/p/281150493。

library(bruceR)
## 
## bruceR (v2023.9)
## Broadly Useful Convenient and Efficient R functions
## 
## Packages also loaded:
## ✔ data.table ✔ emmeans
## ✔ dplyr      ✔ lmerTest
## ✔ tidyr      ✔ effectsize
## ✔ stringr    ✔ performance
## ✔ ggplot2    ✔ interactions
## 
## Main functions of `bruceR`:
## cc()             Describe()  TTEST()
## add()            Freq()      MANOVA()
## .mean()          Corr()      EMMEANS()
## set.wd()         Alpha()     PROCESS()
## import()         EFA()       model_summary()
## print_table()    CFA()       lavaan_summary()
## 
## For full functionality, please install all dependencies:
## install.packages("bruceR", dep=TRUE)
## 
## Online documentation:
## https://psychbruce.github.io/bruceR
## 
## To use this package in publications, please cite:
## Bao, H.-W.-S. (2023). bruceR: Broadly useful convenient and efficient R functions (Version 2023.9) [Computer software]. https://CRAN.R-project.org/package=bruceR

2.2 数据整理

数据框（data frame）是最常见的类型。R语言默认的数据框格式会展示出所有的数据，这就导致阅读数据非常不便。尽管可以使用head()选择展示部分数据，但依旧不便。

# 示例数据iris。如果不设置head()，将会展示全部150行数据
head(iris, 10)
##    Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1           5.1         3.5          1.4         0.2  setosa
## 2           4.9         3.0          1.4         0.2  setosa
## 3           4.7         3.2          1.3         0.2  setosa
## 4           4.6         3.1          1.5         0.2  setosa
## 5           5.0         3.6          1.4         0.2  setosa
## 6           5.4         3.9          1.7         0.4  setosa
## 7           4.6         3.4          1.4         0.3  setosa
## 8           5.0         3.4          1.5         0.2  setosa
## 9           4.4         2.9          1.4         0.2  setosa
## 10          4.9         3.1          1.5         0.1  setosa

tidyverse提供了便捷的tibble格式。它和数据框格式同理，默认情况下只会展示部分数据，同时会显示每列数据的数据格式，方便后期处理。使用as_tibble()即可将导入的数据转变为tibble格式。

# 使用as_tibble将iris转变为tibble格式
as_tibble(iris)
## # A tibble: 150 × 5
##    Sepal.Length Sepal.Width Petal.Length Petal.Width Species
##           <dbl>       <dbl>        <dbl>       <dbl> <fct>  
##  1          5.1         3.5          1.4         0.2 setosa 
##  2          4.9         3            1.4         0.2 setosa 
##  3          4.7         3.2          1.3         0.2 setosa 
##  4          4.6         3.1          1.5         0.2 setosa 
##  5          5           3.6          1.4         0.2 setosa 
##  6          5.4         3.9          1.7         0.4 setosa 
##  7          4.6         3.4          1.4         0.3 setosa 
##  8          5           3.4          1.5         0.2 setosa 
##  9          4.4         2.9          1.4         0.2 setosa 
## 10          4.9         3.1          1.5         0.1 setosa 
## # ℹ 140 more rows

2.2.1 数据导入 (import)

tidyverse包提供了大量的数据导入函数（如导入txt文档的read_delim()，导入csv文档的read_csv()等)。我们更推荐使用bruceR包改写的import()函数，它可以导入任何格式的数据，包括但不限于txt、xlxs、csv、sav等，甚至可以指定导入第几个sheet或者单元格范围，操作上更为便捷。

# 使用import函数进行导入，同时将导入的文档设置为tibble格式
data <- as_tibble(import('Count.xlsx'))

# 如果数据文件并没有在默认文件夹，使用file.choose()语句
data <- as_tibble(import(file.choose()))

2.2.2 数据规整 (reshape)

常见数据框的格式如iris所示，每列代表了一种变量及其数据，这种格式被称为宽数据（wide format）格式。在更多时候，如回归分析或可视化，更多采用长数据（long format）格式。在这种格式下，一行只会对应到一个数据。在tidyr包中，作者提供了gather()和pivot_longer()两个函数将宽数据转换为长数据。后者是前者的更新升级版，并且逐渐替代前者，因为它有更多操作空间，且前者不再维护和更新。因此，我们建议使用pivot_longer()函数。

# pivot_longer函数基本语法
pivot_longer(data, cols, names_to = "name", values_to = "value",
             names_prefix = ..., values_drop_na = ...)

# cols: 要进行规整的列名（或者列数）
# names_to: 规整后变量的列名
# values_to: 规整后数值的列名
# names_prefix: 删除变量中重复的前缀
# values_drop_na: 是否去除NA值（TRUE/FALSE）
# 示例数据：diamonds
diamonds
## # A tibble: 53,940 × 10
##    carat cut       color clarity depth table price     x     y     z
##    <dbl> <ord>     <ord> <ord>   <dbl> <dbl> <int> <dbl> <dbl> <dbl>
##  1  0.23 Ideal     E     SI2      61.5    55   326  3.95  3.98  2.43
##  2  0.21 Premium   E     SI1      59.8    61   326  3.89  3.84  2.31
##  3  0.23 Good      E     VS1      56.9    65   327  4.05  4.07  2.31
##  4  0.29 Premium   I     VS2      62.4    58   334  4.2   4.23  2.63
##  5  0.31 Good      J     SI2      63.3    58   335  4.34  4.35  2.75
##  6  0.24 Very Good J     VVS2     62.8    57   336  3.94  3.96  2.48
##  7  0.24 Very Good I     VVS1     62.3    57   336  3.95  3.98  2.47
##  8  0.26 Very Good H     SI1      61.9    55   337  4.07  4.11  2.53
##  9  0.22 Fair      E     VS2      65.1    61   337  3.87  3.78  2.49
## 10  0.23 Very Good H     VS1      59.4    61   338  4     4.05  2.39
## # ℹ 53,930 more rows
# 折叠“三维”（x、y、z）数据
diamonds %>% 
  pivot_longer(x:z, names_to = "Dimensions", values_to = "Values")
## # A tibble: 161,820 × 9
##    carat cut     color clarity depth table price Dimensions Values
##    <dbl> <ord>   <ord> <ord>   <dbl> <dbl> <int> <chr>       <dbl>
##  1  0.23 Ideal   E     SI2      61.5    55   326 x            3.95
##  2  0.23 Ideal   E     SI2      61.5    55   326 y            3.98
##  3  0.23 Ideal   E     SI2      61.5    55   326 z            2.43
##  4  0.21 Premium E     SI1      59.8    61   326 x            3.89
##  5  0.21 Premium E     SI1      59.8    61   326 y            3.84
##  6  0.21 Premium E     SI1      59.8    61   326 z            2.31
##  7  0.23 Good    E     VS1      56.9    65   327 x            4.05
##  8  0.23 Good    E     VS1      56.9    65   327 y            4.07
##  9  0.23 Good    E     VS1      56.9    65   327 z            2.31
## 10  0.29 Premium I     VS2      62.4    58   334 x            4.2 
## # ℹ 161,810 more rows

相对应的，tidyr包提供了进行反向操作的函数，即把长数据转变为宽数据，对应使用spread()和pivot_wider()函数。同样理由，我们更推荐使用后者。pivot_wider()和pivot_longer()使用的参数相似，也有一些独特功能。

# 基本语法
pivot_wider(data, names_from = "name", values_from = "value", 
            values_fill = ..., values_fn = ...)

# values_fill: 当扩展数据框时，如果出现了空值，默认填充为NA，可以借此指定填充为什么熟知
# values_fn: 可以进行计算，如取平均值
# 示例数据：fish_encouters
fish_encounters
## # A tibble: 114 × 3
##    fish  station  seen
##    <fct> <fct>   <int>
##  1 4842  Release     1
##  2 4842  I80_1       1
##  3 4842  Lisbon      1
##  4 4842  Rstr        1
##  5 4842  Base_TD     1
##  6 4842  BCE         1
##  7 4842  BCW         1
##  8 4842  BCE2        1
##  9 4842  BCW2        1
## 10 4842  MAE         1
## # ℹ 104 more rows
# 扩展数据
fish_encounters %>% 
  pivot_wider(names_from = station, values_from = seen)
## # A tibble: 10 × 12
##    fish  Release I80_1 Lisbon  Rstr Base_TD   BCE   BCW  BCE2  BCW2   MAE   MAW
##    <fct>   <int> <int>  <int> <int>   <int> <int> <int> <int> <int> <int> <int>
##  1 4842        1     1      1     1       1     1     1     1     1     1     1
##  2 4843        1     1      1     1       1     1     1     1     1     1     1
##  3 4844        1     1      1     1       1     1     1     1     1     1     1
##  4 4845        1     1      1     1       1    NA    NA    NA    NA    NA    NA
##  5 4847        1     1      1    NA      NA    NA    NA    NA    NA    NA    NA
##  6 4848        1     1      1     1      NA    NA    NA    NA    NA    NA    NA
##  7 4849        1     1     NA    NA      NA    NA    NA    NA    NA    NA    NA
##  8 4850        1     1     NA     1       1     1     1    NA    NA    NA    NA
##  9 4851        1     1     NA    NA      NA    NA    NA    NA    NA    NA    NA
## 10 4854        1     1     NA    NA      NA    NA    NA    NA    NA    NA    NA
# 指定填充为0
fish_encounters %>% 
  pivot_wider(names_from = station, values_from = seen, values_fill = 0)
## # A tibble: 10 × 12
##    fish  Release I80_1 Lisbon  Rstr Base_TD   BCE   BCW  BCE2  BCW2   MAE   MAW
##    <fct>   <int> <int>  <int> <int>   <int> <int> <int> <int> <int> <int> <int>
##  1 4842        1     1      1     1       1     1     1     1     1     1     1
##  2 4843        1     1      1     1       1     1     1     1     1     1     1
##  3 4844        1     1      1     1       1     1     1     1     1     1     1
##  4 4845        1     1      1     1       1     0     0     0     0     0     0
##  5 4847        1     1      1     0       0     0     0     0     0     0     0
##  6 4848        1     1      1     1       0     0     0     0     0     0     0
##  7 4849        1     1      0     0       0     0     0     0     0     0     0
##  8 4850        1     1      0     1       1     1     1     0     0     0     0
##  9 4851        1     1      0     0       0     0     0     0     0     0     0
## 10 4854        1     1      0     0       0     0     0     0     0     0     0

2.2.3 数据整合和切分 (combine/split)

在部分情况下，我们需要提取某一列数据中的部分内容，生成新的一列，又或者反其道而行之，将两列内容进行合并。之前的老办法可以使用filter()函数配合mutate()进行处理。在tidyr包中，其实提供了更为便捷的unite()和separate()函数。前者将多列数值进行合并，后者则是将一列数据进行拆分。

# 多列数据合并
unite(data, col1, col2, col3, ..., sep = ..., remove = TRUE, na.rm = FALSE)

# sep: 分隔符，如不需要分隔符，则是sep = ""
# remove: 是否删除原列 TRUE/FALSE
# na.rm: 是否删除NA值 TRUE/FALSE
# 示例数据：table5
table5
## # A tibble: 6 × 4
##   country     century year  rate             
##   <chr>       <chr>   <chr> <chr>            
## 1 Afghanistan 19      99    745/19987071     
## 2 Afghanistan 20      00    2666/20595360    
## 3 Brazil      19      99    37737/172006362  
## 4 Brazil      20      00    80488/174504898  
## 5 China       19      99    212258/1272915272
## 6 China       20      00    213766/1280428583
# 合并century和year
table5 %>% 
  unite(century, year, col = "year", sep = "")
## # A tibble: 6 × 3
##   country     year  rate             
##   <chr>       <chr> <chr>            
## 1 Afghanistan 1999  745/19987071     
## 2 Afghanistan 2000  2666/20595360    
## 3 Brazil      1999  37737/172006362  
## 4 Brazil      2000  80488/174504898  
## 5 China       1999  212258/1272915272
## 6 China       2000  213766/1280428583

将一列数据进行拆分，需要分两种情况：一种是拆分后的数据分裂成多列，这时候我们使用separate()即可；另一种情况是，将拆分后的每一个单元格的数据作为行插入到该列，即列数不会变化，这时候就需要separate_rows()函数。

# 情况1: 拆分为多列
separate(data, col, into = ..., sep = "[^[:alnum:]]+", remove = TRUE)

# col: 要拆分到列
# into: 新列的列名
# sep: 识别拆分列的分隔符，分两种情况
#      如果列中分隔符是字符型数据，直接写对应的字符即可
#      如果列值是数值型数据，则1代表最左边，-1代表最右边

separate_rows(...) # 参数基本同上
# 示例数据：table3
table3
## # A tibble: 6 × 3
##   country      year rate             
##   <chr>       <dbl> <chr>            
## 1 Afghanistan  1999 745/19987071     
## 2 Afghanistan  2000 2666/20595360    
## 3 Brazil       1999 37737/172006362  
## 4 Brazil       2000 80488/174504898  
## 5 China        1999 212258/1272915272
## 6 China        2000 213766/1280428583
# 拆分为多个列
table3 %>% 
  separate(rate, sep = "/", into = c("cases", "pop"))
## # A tibble: 6 × 4
##   country      year cases  pop       
##   <chr>       <dbl> <chr>  <chr>     
## 1 Afghanistan  1999 745    19987071  
## 2 Afghanistan  2000 2666   20595360  
## 3 Brazil       1999 37737  172006362 
## 4 Brazil       2000 80488  174504898 
## 5 China        1999 212258 1272915272
## 6 China        2000 213766 1280428583
# 拆分并将单元格插入为行数据
table3 %>% 
  separate_rows(rate, sep = "/")
## # A tibble: 12 × 3
##    country      year rate      
##    <chr>       <dbl> <chr>     
##  1 Afghanistan  1999 745       
##  2 Afghanistan  1999 19987071  
##  3 Afghanistan  2000 2666      
##  4 Afghanistan  2000 20595360  
##  5 Brazil       1999 37737     
##  6 Brazil       1999 172006362 
##  7 Brazil       2000 80488     
##  8 Brazil       2000 174504898 
##  9 China        1999 212258    
## 10 China        1999 1272915272
## 11 China        2000 213766    
## 12 China        2000 1280428583

2.2.4 NA值

NA值又叫做缺省值，当导入的数据里存在“空白”位置，导入到R会识别为缺省值。在大多数情况下，我们需要将NA值过滤，因为会影响后续的数据计算及处理。除了使用filter()函数配合is.na()进行判断外，tidyr包还提供了一系列对NA值操作的函数，包括剔除和填充两个方面。

df <- tibble(x = c("A", "B", "C", "D", "E"), y = c(1, NA, NA, 3, NA))

# 剔除包含NA值的行
df %>% 
  drop_na(y)
## # A tibble: 2 × 2
##   x         y
##   <chr> <dbl>
## 1 A         1
## 2 D         3
# 使用临近值填充NA
df %>% 
  fill(y)
## # A tibble: 5 × 2
##   x         y
##   <chr> <dbl>
## 1 A         1
## 2 B         1
## 3 C         1
## 4 D         3
## 5 E         3
# 指定填充的数值
df %>% 
  replace_na(list(y = 2))
## # A tibble: 5 × 2
##   x         y
##   <chr> <dbl>
## 1 A         1
## 2 B         2
## 3 C         2
## 4 D         3
## 5 E         2

2.3 数据操作

当数据导入并规整结束后，就要正式进入数据分析的过程。在这个过程中，为了更清楚地细致地考察数据的细节表现，我们通常会对数据进行有条件地筛选，提取子集，或者合并不同的数据框。这一系列数据操作，R语言也帮我们提供了便捷的函数。

2.3.1 筛选与增加

最基本的数据操作函数，可以分为按行操作和按列操作两个大方向。按行操作的常见函数包括：filter()，add_row()。按列操作的常见函数包括：select()，relocate()，mutate()。具体例子如下所示。

## 按行操作
# 筛选出符合条件的行，支持逻辑运算符
filter(diamonds, depth > 60)
## # A tibble: 48,315 × 10
##    carat cut       color clarity depth table price     x     y     z
##    <dbl> <ord>     <ord> <ord>   <dbl> <dbl> <int> <dbl> <dbl> <dbl>
##  1  0.23 Ideal     E     SI2      61.5    55   326  3.95  3.98  2.43
##  2  0.29 Premium   I     VS2      62.4    58   334  4.2   4.23  2.63
##  3  0.31 Good      J     SI2      63.3    58   335  4.34  4.35  2.75
##  4  0.24 Very Good J     VVS2     62.8    57   336  3.94  3.96  2.48
##  5  0.24 Very Good I     VVS1     62.3    57   336  3.95  3.98  2.47
##  6  0.26 Very Good H     SI1      61.9    55   337  4.07  4.11  2.53
##  7  0.22 Fair      E     VS2      65.1    61   337  3.87  3.78  2.49
##  8  0.3  Good      J     SI1      64      55   339  4.25  4.28  2.73
##  9  0.23 Ideal     J     VS1      62.8    56   340  3.93  3.9   2.46
## 10  0.22 Premium   F     SI1      60.4    61   342  3.88  3.84  2.33
## # ℹ 48,305 more rows
# 添加一行或多行数据
df <- tibble(x = c("A", "B", "C", "D", "E"), y = c(1, 3, 5, 3, 6))
add_row(df, x = "F", y = 9)
## # A tibble: 6 × 2
##   x         y
##   <chr> <dbl>
## 1 A         1
## 2 B         3
## 3 C         5
## 4 D         3
## 5 E         6
## 6 F         9
## 按列操作
# 筛选符合条件的列（例：筛选出color到price的列）支持contain等参数
select(diamonds, color:price)
## # A tibble: 53,940 × 5
##    color clarity depth table price
##    <ord> <ord>   <dbl> <dbl> <int>
##  1 E     SI2      61.5    55   326
##  2 E     SI1      59.8    61   326
##  3 E     VS1      56.9    65   327
##  4 I     VS2      62.4    58   334
##  5 J     SI2      63.3    58   335
##  6 J     VVS2     62.8    57   336
##  7 I     VVS1     62.3    57   336
##  8 H     SI1      61.9    55   337
##  9 E     VS2      65.1    61   337
## 10 H     VS1      59.4    61   338
## # ℹ 53,930 more rows
# 移动列的位置
relocate(diamonds, x, y, z, .after = color)
## # A tibble: 53,940 × 10
##    carat cut       color     x     y     z clarity depth table price
##    <dbl> <ord>     <ord> <dbl> <dbl> <dbl> <ord>   <dbl> <dbl> <int>
##  1  0.23 Ideal     E      3.95  3.98  2.43 SI2      61.5    55   326
##  2  0.21 Premium   E      3.89  3.84  2.31 SI1      59.8    61   326
##  3  0.23 Good      E      4.05  4.07  2.31 VS1      56.9    65   327
##  4  0.29 Premium   I      4.2   4.23  2.63 VS2      62.4    58   334
##  5  0.31 Good      J      4.34  4.35  2.75 SI2      63.3    58   335
##  6  0.24 Very Good J      3.94  3.96  2.48 VVS2     62.8    57   336
##  7  0.24 Very Good I      3.95  3.98  2.47 VVS1     62.3    57   336
##  8  0.26 Very Good H      4.07  4.11  2.53 SI1      61.9    55   337
##  9  0.22 Fair      E      3.87  3.78  2.49 VS2      65.1    61   337
## 10  0.23 Very Good H      4     4.05  2.39 VS1      59.4    61   338
## # ℹ 53,930 more rows
# 添加新的一列
mutate(diamonds, Dimensions_all = x + y + z)
## # A tibble: 53,940 × 11
##    carat cut    color clarity depth table price     x     y     z Dimensions_all
##    <dbl> <ord>  <ord> <ord>   <dbl> <dbl> <int> <dbl> <dbl> <dbl>          <dbl>
##  1  0.23 Ideal  E     SI2      61.5    55   326  3.95  3.98  2.43           10.4
##  2  0.21 Premi… E     SI1      59.8    61   326  3.89  3.84  2.31           10.0
##  3  0.23 Good   E     VS1      56.9    65   327  4.05  4.07  2.31           10.4
##  4  0.29 Premi… I     VS2      62.4    58   334  4.2   4.23  2.63           11.1
##  5  0.31 Good   J     SI2      63.3    58   335  4.34  4.35  2.75           11.4
##  6  0.24 Very … J     VVS2     62.8    57   336  3.94  3.96  2.48           10.4
##  7  0.24 Very … I     VVS1     62.3    57   336  3.95  3.98  2.47           10.4
##  8  0.26 Very … H     SI1      61.9    55   337  4.07  4.11  2.53           10.7
##  9  0.22 Fair   E     VS2      65.1    61   337  3.87  3.78  2.49           10.1
## 10  0.23 Very … H     VS1      59.4    61   338  4     4.05  2.39           10.4
## # ℹ 53,930 more rows
# 根据条件增加新的一列: case_when
mutate(diamonds, Quality = case_when(
  x >= 4 ~ "Excellent", 
  x >= 3 & x < 4 ~ "Good", 
  x < 3 ~ "Bad"
))
## # A tibble: 53,940 × 11
##    carat cut       color clarity depth table price     x     y     z Quality  
##    <dbl> <ord>     <ord> <ord>   <dbl> <dbl> <int> <dbl> <dbl> <dbl> <chr>    
##  1  0.23 Ideal     E     SI2      61.5    55   326  3.95  3.98  2.43 Good     
##  2  0.21 Premium   E     SI1      59.8    61   326  3.89  3.84  2.31 Good     
##  3  0.23 Good      E     VS1      56.9    65   327  4.05  4.07  2.31 Excellent
##  4  0.29 Premium   I     VS2      62.4    58   334  4.2   4.23  2.63 Excellent
##  5  0.31 Good      J     SI2      63.3    58   335  4.34  4.35  2.75 Excellent
##  6  0.24 Very Good J     VVS2     62.8    57   336  3.94  3.96  2.48 Good     
##  7  0.24 Very Good I     VVS1     62.3    57   336  3.95  3.98  2.47 Good     
##  8  0.26 Very Good H     SI1      61.9    55   337  4.07  4.11  2.53 Excellent
##  9  0.22 Fair      E     VS2      65.1    61   337  3.87  3.78  2.49 Good     
## 10  0.23 Very Good H     VS1      59.4    61   338  4     4.05  2.39 Excellent
## # ℹ 53,930 more rows

2.3.2 数据框合并

在一些情况下，我们需要将不同的数据框进行合并处理。比如，在采集数据的过程中，Praat脚本提取的声学数据，和被试的个人信息数据，是分别存储的。如果在Excel中一个个对应处理，会稍显繁琐。除了基础的rbind()和cbind()函数外，更多情况下我们需要将相同列名或者相同行数的数据进行归并，而不是简单地数据框叠加。R语言为我们提供了一系列的._join()函数，通过匹配和归并，实现了这些功能。下面我们介绍这一系列函数。

# 创建假数据
A <- tibble(Alphabet = c("a", "b", "c"), 
            Value = c(18, 17, 25))

B <- tibble(Alphabet = c("a", "b", "d"),
            Value = c(35, 79, 43))

# 保留左侧（第一个数据框）的所有值，将右侧数据合并到左侧，取交集
left_join(A, B, by = "Alphabet")
## # A tibble: 3 × 3
##   Alphabet Value.x Value.y
##   <chr>      <dbl>   <dbl>
## 1 a             18      35
## 2 b             17      79
## 3 c             25      NA
# 保留右侧（第二个数据框）的所有值，将左侧数据合并到右侧，取交集
right_join(A, B, by = "Alphabet")
## # A tibble: 3 × 3
##   Alphabet Value.x Value.y
##   <chr>      <dbl>   <dbl>
## 1 a             18      35
## 2 b             17      79
## 3 d             NA      43
# 合并两个数据框，仅保留by指定的列中值一样的数据（即只取交集）
inner_join(A, B, by = "Alphabet")
## # A tibble: 2 × 3
##   Alphabet Value.x Value.y
##   <chr>      <dbl>   <dbl>
## 1 a             18      35
## 2 b             17      79
# 合并两个数据框，保留所有数据，无匹配的数据使用NA代替
full_join(A, B, by = "Alphabet")
## # A tibble: 4 × 3
##   Alphabet Value.x Value.y
##   <chr>      <dbl>   <dbl>
## 1 a             18      35
## 2 b             17      79
## 3 c             25      NA
## 4 d             NA      43

这些函数和SQL中的合并函数是异曲同工之妙。下面我们使用Venn图进行图解，或许可以帮助大家更直观地了解不同函数之间的异同。

2.4 可视化建议

数据可视化可以帮我们直观地了解数据分布情况，甚至是辅助我们了解统计分析中需要检验的基本假设，或者采取何种模型。然而，数据可视化是一个高度「数据驱动」的内容，即个性化程度比较高。在此，我们主要推荐如下两个网站。

  1. ggplot2官方网站（https://ggplot2.tidyverse.org/）：ggplot2是R目前可视化的主流包，基本满足了所有的可视化结果方案。在该网站上提供了相关的cheatsheet，方便快速检阅所需要的图该调用何种函数。
  2. ggplot2相关扩展包（https://exts.ggplot2.tidyverse.org/gallery/）：众多研究者围绕它开发了一系列扩展包，比如着色方案相关的ggsci，统计分析结果一键美化的ggstatsplot，标注信息美化相关的ggrepel和ggpubr等等。该网站即是这些扩展包的集合网站，可以查询自己需要的包，在R中进行安装加载。

3 统计分析

统计检验是每一个数据科学不可绕开的部分。当我们进行数据描述时，需要一些证据证明不同组别或不同类型之间的数据「确实」存在差异，而不是因为随机因素引起的。比如考察一个班级的学习水平，随机选择的10个学生恰好是班级前20名的学生，就不能很好反应该班级真正的水平。因此，我们需要统计检验对数据进行检测。

有一点需要特别强调一下：统计检验并不是越高级的越好。在阅读海量的论文时，很多人总会被文章中看起来高大上的统计方法而迷住双眼，忽略了选择某种统计方法的原因。根据研究目的/研究问题选择合适的统计方法，而不是忽略研究问题无脑选择高级统计方法，这才是最重要的。举一个简单的例子：假设在考察汉语水平（高、中、低）对句子理解的正确率分析中，我想要对比「汉语水平」因素对正确率的影响，为了回答这个问题，我可以选择ANOVA，或者选择更合适的LMM（考虑随机效应）。如果我想考察「不同水平的正确率是否显著高于平均正确率（如47%）」，那么单尾t检验就足够解答问题。

在下面的介绍里，我们默认读者已经有了统计学的一些基本知识（如集中趋势、描述统计、推论统计、零假设等等），从基础的t检验和ANOVA开始，着重点在于不同的回归分析中，最后针对常见的非线性轮廓(contour)，以语音学的F0为例，介绍一些常用的其他方法。

3.1 t检验与ANOVA

t检验和ANOVA都是对不同组别均值进行比较，看是否存在差异性/变异性。二者最主要的区别在于，t检验适用于两组之间的比较，而当需要比较的组别/水平超过了两组，则需要使用ANOVA。因为t检验和ANOVA均属于参数检验(parametric test, 假定数据服从正态分布)，因此它们有一些共通的基本假设：

  • 数据正态分布
  • 数据彼此独立
  • 同一连续变量标度

因此，在选择这两个统计方法之前，确保数据符合基本假设是关键一步。如果不符合任一假设，则需要使用对应的非参数检验(non-parametric test, 如Welch’s test等)。

3.1.1 t检验类型与R实现方法

根据研究目的和研究问题的不同，我们需要比较的对象也有所不同。在语言学研究中，常见的t检验类型有以下三种：

  1. 单样本t检验(one-way t-test): 检测实验获得的数据是否和（已知）「特定值」存在显著差异。如上面例子中，不同汉语水平的句子理解正确率是否显著差异于平均正确率水平。
  2. 独立样本t检验(independent t-test): 也可以称作组间t检验(between-groups t-test)，检测来自「不同组别」之间的数据是否存在差异。如，年轻人和老年人的音高数值的差异。
  3. 配对样本t检验(paired t-test): 也可以称作组内t检验(within-groups t-test)，检测同一组别下，不同样本之间的差异。这是语言学研究中「最常见的」比较方法。如，同一批被试里，宽焦点和窄焦点两种条件下音高的差异。

R语言里内置了t检验的函数t.test()，不需要加载任何包，根据不同的检验类型，可以在函数中进行指定。比如针对是否是配对样本检验，通过更改paired = ...布尔值即可，如t.test(data1, data2, paired = FALSE)表示独立样本t检验。如果是单样本，则使用mu = ...设置特定值，如t.test(data, mu = 50)。

然而，baseR中的t.tset()在实际操作中使用起来并不方便，最主要的不便在于效应检测和结果的输出上，还会浪费时间整理统计结果。受益于bruceR包的开发，其中的TTEST()函数完美解决了这个痛点。TTEST()函数的主要参数可以通过帮助文档以及本文开头提供的中文帮助文档查阅。下面我们使用iris数据进行举例。

# 加载bruceR包
library(bruceR)

# 准备独立样本t检验数据，只保留其中两个组别
iris_it <- iris %>% 
  filter(Species != "virginica") %>% 
  droplevels()

# 准备单样本和配对样本t检验数据
iris_dt <- iris %>% 
  filter(Species == "setosa") %>% 
  droplevels()

# 单样本t检验，假设setosa的萼片平均长度为3，检测采集到的数据（比如是在比较好的环境下成长的）是否显著优于平均长度。采用单尾单样本t检验。
TTEST(iris_dt, "Sepal.Length", test.value = 3, test.sided = ">")
## 
## One-Sample t-test
## 
## Hypothesis: one-sided (μ > 3)
## 
## Descriptives:
## ────────────────────────────
##      Variable  N Mean (S.D.)
## ────────────────────────────
##  Sepal.Length 50 5.01 (0.35)
## ────────────────────────────
## 
## Results of t-test:
## ────────────────────────────────────────────────────────────────────────────────────────────────────
##                                       t df     p     Difference [95% CI] Cohen’s d [95% CI]     BF10
## ────────────────────────────────────────────────────────────────────────────────────────────────────
## Sepal.Length: (Sepal.Length - 3)  40.24 49 <.001 ***    2.01 [1.92, Inf]   5.69 [5.45, Inf] 9.95e+35
## ────────────────────────────────────────────────────────────────────────────────────────────────────

# 独立样本t检验，假设考察问题为：不同种类的鸢尾花在萼片长度是否有差异
TTEST(iris_it, y = "Sepal.Length", x = "Species")
## 
## Independent-Samples t-test
## 
## Hypothesis: two-sided (μ2 - μ1 ≠ 0)
## 
## Descriptives:
## ───────────────────────────────────────────────
##      Variable  Factor      Level  N Mean (S.D.)
## ───────────────────────────────────────────────
##  Sepal.Length Species setosa     50 5.01 (0.35)
##  Sepal.Length Species versicolor 50 5.94 (0.52)
## ───────────────────────────────────────────────
## 
## Levene’s test for homogeneity of variance:
## ─────────────────────────────────────────────────────────────────────────
##                                              Levene’s F df1 df2     p    
## ─────────────────────────────────────────────────────────────────────────
## Sepal.Length: Species (versicolor - setosa)        8.43   1  98  .005 ** 
## ─────────────────────────────────────────────────────────────────────────
## Note: H0 = equal variance (homoscedasticity).
## If significant (violation of the assumption),
## then you should better set `var.equal=FALSE`.
## 
## Results of t-test:
## ───────────────────────────────────────────────────────────────────────────────────────────────────────────────
##                                                  t df     p     Difference [95% CI] Cohen’s d [95% CI]     BF10
## ───────────────────────────────────────────────────────────────────────────────────────────────────────────────
## Sepal.Length: Species (versicolor - setosa)  10.52 98 <.001 ***   0.93 [0.75, 1.11]  2.10 [1.71, 2.50] 4.19e+14
## ───────────────────────────────────────────────────────────────────────────────────────────────────────────────

# 配对样本t检验，假设对setosa类的鸢尾花，分别在两种不同实验条件下测试了长度。假设条件是Sepal和Petal
TTEST(iris_dt, y = c("Sepal.Length", "Petal.Length"), paired = TRUE)
## 
## Paired-Samples t-test
## 
## Hypothesis: two-sided (μ2 - μ1 ≠ 0)
## 
## Descriptives:
## ────────────────────────────
##      Variable  N Mean (S.D.)
## ────────────────────────────
##  Sepal.Length 50 5.01 (0.35)
##  Petal.Length 50 1.46 (0.17)
## ────────────────────────────
## 
## Results of t-test:
## ───────────────────────────────────────────────────────────────────────────────────────────────────────────────
##                                             t df     p      Difference [95% CI]     Cohen’s d [95% CI]     BF10
## ───────────────────────────────────────────────────────────────────────────────────────────────────────────────
## Paired: (Petal.Length - Sepal.Length)  -71.84 49 <.001 *** -3.54 [-3.64, -3.44] -10.16 [-10.44, -9.87] 3.63e+47
## ───────────────────────────────────────────────────────────────────────────────────────────────────────────────

# 如果需要把结果输出到文档，使用file = "..."即可。如：TTEST(iris_it, y = "Sepal.Length", x = "Species", file = "iris_it_result.docx")

3.1.2 ANOVA与R实现

上面提到的t检验，适用于两个组别（或者说，某一因素下两个水平）之间的比较。如果多于两个组别/水平，我们当然可以进行多个t检验，但是这会增加I类错误(Type I error)，也就是假阳性（零假设本为真，却拒绝了零假设，也被称为弃真错误）。

比如，一个t检验，我们的判断出现失误的概率是5%(\(H_0\): 两者之间没有差异。数据检测出来是有显著差异的，我们应该拒绝\(H_0\)，这个“拒绝”出现失误的概率是5%)。如果我们对3个组别进行了两两t检验，每一次检验都有5%的概率出现判断失误，那么从总体来看，这个犯错的概率就会增加到\((1-0.95^3)\)，也就是15%，这增加了我们判断出现错误的概率。如果组别有4个，那么就有6次两两比较，错误概率也就增加到了27%！也因此，在任何涉及到多次比较的统计分析中，尤其是事后检验(post-hoc analysis)，我们总要谨慎对待得出的结果，或者进行校正。

为了防止这样的情况发生，面对大于两个组别的检验，我们通常使用ANOVA(ANalysis Of VAriance)。AVNOA也有不同的类型，具体如下图所示。

我们不再赘述ANOVA背后的具体逻辑。在R中，有一系列可以使用ANOVA的函数，包括aov()，anova()，ez::ezANOVA()等。同样的，bruceR包提供了MANOVA()函数，尽管是以多元方差分析命名，但是根据具体的参数设定使用ANOVA进行检验，并可以将结果导出到文件中。因为ANOVA类型过多，我们以bruceR包的数据集进行最简单的组间分析示例。

# 组间比较
# 数据集：between.2
MANOVA(between.2, dv = "SCORE", between = c("A", "B"))
## 
## ====== ANOVA (Between-Subjects Design) ======
## 
## Descriptives:
## ─────────────────────────
##  "A" "B"   Mean    S.D. n
## ─────────────────────────
##   A1  B1  4.000 (1.414) 4
##   A1  B2  4.000 (1.633) 4
##   A1  B3  4.750 (2.062) 4
##   A2  B1  3.750 (0.957) 4
##   A2  B2  8.000 (0.816) 4
##   A2  B3 12.000 (0.816) 4
## ─────────────────────────
## Total sample size: N = 24
## 
## ANOVA Table:
## Dependent variable(s):      SCORE
## Between-subjects factor(s): A, B
## Within-subjects factor(s):  –
## Covariate(s):               –
## ─────────────────────────────────────────────────────────────────────
##            MS   MSE df1 df2      F     p     η²p [90% CI of η²p]  η²G
## ─────────────────────────────────────────────────────────────────────
## A      80.667 1.861   1  18 43.343 <.001 ***   .707 [.482, .817] .707
## B      40.542 1.861   2  18 21.784 <.001 ***   .708 [.470, .815] .708
## A * B  28.292 1.861   2  18 15.201 <.001 ***   .628 [.347, .763] .628
## ─────────────────────────────────────────────────────────────────────
## MSE = mean square error (the residual variance of the linear model)
## η²p = partial eta-squared = SS / (SS + SSE) = F * df1 / (F * df1 + df2)
## ω²p = partial omega-squared = (F - 1) * df1 / (F * df1 + df2 + 1)
## η²G = generalized eta-squared (see Olejnik & Algina, 2003)
## Cohen’s f² = η²p / (1 - η²p)
## 
## Levene’s Test for Homogeneity of Variance:
## ───────────────────────────────────────
##            Levene’s F df1 df2     p    
## ───────────────────────────────────────
## DV: SCORE       0.605   5  18  .697    
## ───────────────────────────────────────

3.1.3 无法解决的问题

至此，我们了解到针对不同组别的差异，可以使用t检验或ANOVA。不难发现，它们在本质上与线性回归模型(linear regression model)是相似的，甚至可以把ANOVA看作是线性回归模型的子集，也正是这样的观念，在R内置的t检验和ANOVA函数中，可以写作与线性模型一样的结构。

# t检验
t.test(dv ~ iv)

# ANOVA
aov(dv ~ iv)

随着语言学研究者对数据分析的深入了解，越来越多的人意识到，因采样得到的被试，或者使用的语言材料，都不可避免地产生随机误差。换句话说，t检验和ANOVA是默认所有被试在一个条件下的表现是一致的，这种不确定性的随机误差就很容易被抹消掉，这也就导致人们认为这样的结果「并没有完美真实」地反应实际情况。于是，人们转而需求更适合的统计方法，线性混合效应模型(linear mixed-effect model, LMM)便为语言学研究者广泛运用。

3.2 线性混合效应模型（LMM）

语言学的实验研究中，最常见的是采用重复测量的方法，即：使用同一批被试甚至是同一批语言材料，只是条件不一样。如，假设研究汉语不同焦点类型的韵律实现方式。为了控制变量和影响因素，通常我们选择同一批被试，并且实验语料也会是同一批，只是设置了不同的语境。

此外，随着语言学研究者对统计学的深入认识，以及对实验数据解读的需要，人们逐渐意识到，将每一个被试以及每一个语料所带来的随机误差加入到统计分析中，会让统计结果更具有说服力。毕竟我们不可能穷尽世界上所有的人或者所有的语料，那样过于浪费时间，而且每个人的表现都会有所不同，将这种随机误差加进去，就可以排除诸如“被试1和被试15在音高上差异太大而导致了显著差异”这样的问题。那么，考虑到随机误差的线性模型，就被称为线性混合模型。

3.2.1 固定效应和随机效应

在拟合LMM过程中，有两个基本的概念十分重要。

  1. 固定效应(fixed effects)：实验里所控制的因素，或者说你的研究里想要考察的因素，就是所谓的固定效应。我们在t检验、ANOVA，以及在其他教材中看到的使用lm()函数建立的线性模型，里面出现的因素都是固定效应。如上面举的例子中，“焦点类型”就是要放在模型中的固定效应。我们通常假设固定效应是不会出现变化的，可以复现实验，不以实验为转移的。比如，我做了一个焦点的研究，里面设置了宽窄焦点两个类型，那另一个人复现这个实验，设置了新的语料，但还是这两个类型。固定效应可以是分类变量，也可以是连续变量。通过固定效应的检验，我们可以得到截距和斜率两个关键信息。

    • 截距(intercept)：在数学中，截距就是当\(x=0\)时得到的\(y\)的值，在统计中是一样的道理。但在很多情况下，\(x=0\)没有任何意义，尤其是对于分类变量而言。那么，将数据中心化(centering)就可以把均值设定在\((0,0)\)位置，截距就有了实际意义。在使用summary()查看模型固定效应的参数表(coefficient table)里得到的(Intercept)，指的是参考水平的因变量的值。
    • 斜率(slope)：在数学中，斜率即为\(y=kx+b\)中的\(k\)值，在统计分析中，它通常指代变量的其他水平与参考水平相比之间的差异大小。在summray()得到的参数表中，除了(Intercept)外其他的数据，即为差异的情况。

  2. 随机效应(random effects)：没有办法进行精确控制的，会由于从总体采样而产生随机误差的因素，我们称为随机效应。抛开其他研究领域，在语言学研究领域中，特别是心理语言学和语音学，最常见的随机效应包括被试(participants)和实验语料(item)，或者“不恰当”地说，语音学研究里考虑的随机效应只有这两类。“被试”是因为我们不能穷尽采集所有人，“语料”实际上是伪随机，因为我们总是设置一定地条件来设计语料，但是也不可能穷尽所有的语料（比如，考察语句的焦点，没有人能够数出人类语言有多少句话，也无法穷尽）。

从上面的介绍可以看出，固定效应和其他统计方法一致，随机效应是最关键的组成部分，让LMM有别于其他统计方法。随机效应包括随机截距(random intercept)和随机斜率(random slope)两种。和固定效应一样，随机截距和随机斜率本质上是每一个被试/语料的截距和斜率，用来表现随机误差。下面的图示例子是考虑“被试”的随机效应，也就是每个被试的音高可能会有所不同。

那么在R语言中，我们该如何定义随机效应呢？各种教材里写得很清楚。比如在这样的一个模型里：lmer(RT ~ Condition + Group + (1 | Subject), data = data_RT)，(1 | Subject)即为随机效应，这里是随机截距。在R的语法里，随机效应主要有以下的写法（及其对应的含义），在此基础上可以加上不同固定效应的交互作用下的随机效应等。

  1. x + (1 | Subject)：被试随机截距(by-participant varying intercept)
  2. x + (x | Subject)：被试随机截距和随机斜率(by-participant varying slope)，且随机截距和随机斜率有相关性，常见的另一种写法是x + (1 + x | Subject)
  3. x + (x || Subject)：被试随机截距和随机斜率，不考虑二者之间的相关，常见的另一种写法为x + (1 | Subject) + (0 + x | Subject)

这里，我们需要指明有关建模的一个重要问题：究竟考虑什么随机效应。其中，Barr et al. (2013)一文提出的「keep it maximal」原则，让语言学研究者意识到随机斜率的重要性。但是，这也让绝大多数研究者错误地认为模型要包括所有因素的随机斜率，而忽略实际数据的表现情况和研究问题，导致模型都已经报告“无法拟合”(failed to converge)、“畸形拟合”(boundary isSingular)的警告信息，人们还是自动忽略并报告建模结果和p值。针对这一问题，Bates et al. (2015)、Matuschek et al. (2017)、Seedorff et al. (2019)则指出最大化原则的缺陷，建议根据实际数据情况，建立最大可能下的精简模型(parsimonious model)，而非一味追逐囊括所有。

因此，是否囊括随机效应，以及囊括什么程度的随机效应，取决于研究者对数据的分析和描述。假设我们开展了一个行为实验，因变量为正确率，自变量为A（包括两个水平a\(_1\)和a\(_2\)）。如果通过数据观察发现，每个被试之间从水平a\(_1\)到a\(_2\)的差异变异性（变化情况），基本是30%，只是不同被试在a\(_1\)水平和a\(_2\)水平的正确率有所差异，那么这就说明没有出现所谓的被试间随机斜率差异(by-participant varying slope)，这时就无需考虑随机斜率。换个角度，如果你认为被试收到条件影响出现了不同的正确率差异变异性，或者说，存在被试不同的条件A效应，则此时就需要考虑随机斜率。语料(item)的考虑思路和被试是相似的，这里不再赘述。

3.2.2 建模与结果

为了方便模仿建模过程，以及后面分析结果，我们使用languageR包中的lexdec数据进行分析。为了方便，我们把数据存储到一个新的变量名中，并更改为tibble格式。

# 加载languageR包，如果未安装需要先进行安装
library(languageR)

# 转换lexdec数据为tibble格式并存储到新的变量名
lexdec_ti <- as_tibble(lexdec)

# 了解数据基本信息，可以使用summary()函数
summary(lexdec_ti)

通过观察数据，我们可以假设该数据集研究这个问题：不同母语者在词汇决定任务中的表现如何？形态复杂度(morphological complexity)是否会影响词汇决定任务的表现情况？那么，我们主要关注的列就包括：Subject(被试)、RT(因变量reaction time)、NativeLanguage(母语类型)、Word(语料)、Complex(形态复杂度)。

准备好数据后，接下来面对的问题是如何建模。目前最常见的建模方法是逐步回归(stepwise regression)，它分为自下而上(由不包括任何参数的null model到包括所有参数的full model)和自上而下(full model到null model)两个方向。通过比较一系列模型，选择出最合适的进行后续分析。

DV ~ 1 # null model
DV ~ IV1
DV ~ IV1 + IV2
...
DV ~ IV1 + IV2 + ... # full model

逐步回归到方法看起来具有客观性，但是统计学家认为这存在着很多问题，并且很多模拟数据的结果已经证明了问题的存在。在推荐材料里，Winter (2020, pp. 277)中对该方法的问题进行了总结。首先，正如我们在介绍ANOVA的时候提到的，多重比较会增加一类错误的概率，而在逐步回归中，我们是在对多个模型进行比较选择了适合的模型，这导致一类错误的概率大大增加。此外，逐步回归会带来的另一个问题是过度拟合，因为逐步回归是通过比较不同的模型，检测是否要考虑某个因素，但是统计软件对这个的判断过程是完全基于现有数据。我们的实验数据本质上是采样得来，当样本量比较小时，这种判断就有可能“有失偏颇”。基于这些问题，逐步回归这一看似客观的方法，为统计学家和语言学研究者所反对——因为缺乏为何选择该模型以及该模型是否符合研究假设的相关思考。

我们采纳Winter (2020)的建议，在建模时关键考虑的是是否有足够的理由促使我们在模型里囊括某因素，也就是有理论驱动。比如，前人的研究里表明男性和女性在使用时长进行焦点标记上没有差异，那么尽管在实验设计时我们对性别进行了控制（人数大体相当），我们的模型里也完全不需要加入“性别”这一因素（而逐步回归就可能告诉你“性别”因素很重要需要考虑！），也就是说，除非有充分的理由/理论「证明」性别在时长上的差异，否则我们不会将这个因素放到模型中去。因此，建模时我们需要时刻提醒自己：研究问题是什么？建模是为了回答什么问题？考虑什么固定效应和随机效应？。

在观察上面准备的lexdec_ti数据之后，回想研究问题：不同母语者的表现如何？形态复杂度是否对表现有影响？基于此，我们可以确定，该假想实验的因变量为反应时(RT)，自变量包括母语(NativeLanguage)和形态复杂度(Complex)。尽管数据里有其他控制（比如性别、词类、词频等），但是“研究问题”里并不关心这些问题，因此我们不需要加入到模型里去。如果考虑了协变量，也可以将其他因素加进去，但是完全没有任何影响的非研究问题考虑的其他因素，我们就可以忽略。此外，有充分的证据证明不同语言的形态复杂度是不同的，因此语言类型和形态复杂度存在交互效应，因此我们也要考虑到二者之间的交互。最后有关于随机结构，不同被试和不同语料本身会存在随机误差（不可能每个人或者每个词获得的反应时是完全一样的），同时我们可以想到，不同语言的词汇也会引起随机误差，因此这里需要引入由语言类型导致的随机斜率。至此，我们的模型就选好了。当然，你也可以将模型选择交给机器做，即上面提到的逐步回归，通过AIC和BIC这两个参数考虑模型拟合优度（越小，模型拟合越好）。理论驱动和数据驱动相结合的办法也是可以的。

# bruceR包内置了lmerTest包，因此加载bruceR即可
# 建模
mdl <- lmer(RT ~ NativeLanguage * Complex + (1 | Subject) + 
              (0 + NativeLanguage | Word) + (1 | Word), data = lexdec)

# 借助anova函数查看主效应
anova(mdl)
## Type III Analysis of Variance Table with Satterthwaite's method
##                          Sum Sq  Mean Sq NumDF   DenDF F value   Pr(>F)   
## NativeLanguage         0.243578 0.243578     1  20.580  8.3335 0.008944 **
## Complex                0.043348 0.043348     1  76.997  1.4831 0.227012   
## NativeLanguage:Complex 0.119142 0.119142     1 135.153  4.0762 0.045469 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

根据主效应结果，我们就可以看出，语言类型确实存在对反应时的主效应，即：在其他组别一起平均情况下，NativeLanguage对因变量RT具有显著效应（significant difference in NativeLanguage across all other groups)。此外，还存在语言类型和形态复杂度的交互效应。 下面我们来看看通过总结模型得到的各水平的简单效应。

summary(mdl)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: RT ~ NativeLanguage * Complex + (1 | Subject) + (0 + NativeLanguage |  
##     Word) + (1 | Word)
##    Data: lexdec
## 
## REML criterion at convergence: -919.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.7269 -0.6181 -0.1161  0.4610  6.2701 
## 
## Random effects:
##  Groups   Name                  Variance  Std.Dev.  Corr
##  Word     (Intercept)           2.797e-07 0.0005288     
##  Word.1   NativeLanguageEnglish 3.292e-03 0.0573750     
##           NativeLanguageOther   1.059e-02 0.1029210 1.00
##  Subject  (Intercept)           1.847e-02 0.1359223     
##  Residual                       2.923e-02 0.1709638     
## Number of obs: 1659, groups:  Word, 79; Subject, 21
## 
## Fixed effects:
##                                      Estimate Std. Error         df t value
## (Intercept)                          6.323910   0.045960  33.379281 137.595
## NativeLanguageOther                  0.208376   0.066092  26.834739   3.153
## Complexsimplex                      -0.006412   0.025609  84.220565  -0.250
## NativeLanguageOther:Complexsimplex  -0.060171   0.029803 135.153297  -2.019
##                                    Pr(>|t|)    
## (Intercept)                         < 2e-16 ***
## NativeLanguageOther                 0.00396 ** 
## Complexsimplex                      0.80290    
## NativeLanguageOther:Complexsimplex  0.04547 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) NtvLnO Cmplxs
## NtvLnggOthr -0.501              
## Complxsmplx -0.487 -0.010       
## NtvLnggOt:C -0.013 -0.394  0.027

输出的结果包括随机效应和固定效应的结果。随机效应表可以帮助我们了解该模型的随机结构的标准差情况，判断是否个体/语料之间的随机误差比较大。不过我们最关心的是固定效应的参数表。在Fixed effects中，我们主要到有估计值、标准误、自由度、t值和p值（如果使用lme4包进行拟合，p值无法直接显示，因此建议使用lmerTest包进行拟合）。

在参数表的第一行，写着(Intercept)，它正如数学里的含义，理论上\(x=0\)时的因变量的值。但是大多数情况下，\(x=0\)并没有任何意义，所以可以考虑将数据归一化或中心化等等。总之，你可以理解为在没有任何因素效应下因变量的估计值(DV ~ 1)。下面我们看到有NativeLanguageOther，假设我们这个模型是简单回归模型没有交互，这表示的是它和对应的参考水平NativeLanguageEnglish相比的差异性，这个差异性即上面我们提到的斜率。如果我们没有更改对比编码(contrast coding)，那么默认的参考水平是按照字母顺序排列的第一位。根据这一行，我们可以说：其他语言母语者（Other）的反应时要比英语母语者（English）的反应时快0.208秒，且具有显著差异。注意，这里只是NativeLanguage这个因素，不涉及其他因素（假设无交互考虑）。同理，Complexsimplex则是说simplex的反应时比complex慢0.006秒，且并不具有显著差异。

第4行即我们提到的交互效应。因为不同研究者的数据不尽相同，因此交互效应可能是：连续变量和连续变量的交互、分类变量和分类变量的交互、连续变量和分类变量的交互。具体每一种交互效应的解读，可以参考Winter (2020)一书的第八章。我们只对本次示例的分类变量间的交互进行解释。

为什么上面我们要反复提交互效应，因为如果模型里存在交互效应(A*B)，此时上面提到的NativeLanguageOther和Complexsimplex的解读就不是不考虑其他因素了。比如在有NativeLanguage * Complex交互的情况下，输出的参数表里NativeLanguageOther表示的是在Complex的参考水平complex下，其他语言母语者的反应时显著不同于英语母语者。同理，在有交互的情况下，Complexsimplex指在英语母语者的因素下，complex和simplex反应时的差异。因此，在有交互的模型里，这些效应被称为简单效应(simple effect)，因为它描述的是在其他因素的特定水平下，某一因素对因变量的影响。这一点务必理清，不能当作主效应，不然可能会对模型结果造成错误解读。那第4行的交互效应如何解释呢？参数指的是差异程度，那么交互效应的参数，表达的是针对NativeLanguageOther:Complexsimplex两个水平出现的变化程度。至此，我们的参数表也就解析完毕。

最后，如果每个因素出现了多个水平，或者想看特定的亮亮对比，除了变更对比编码外，我们可以借助emmeans包进行事后检验(post-hoc analysis)，正如ANVOA一样，此处不再赘述。

3.2.3 基本假设

线性模型的基本假设与t检验和ANOVA不同，它需要在模型建立好之后、具体分析开始之前，进行基本假设的检测。如果符合基本假设，说明我们的模型是合适的，结果是可靠的，才可以进行后续分析。回忆一下t检验和ANOVA，它们能够正确使用的前提是和原数据的分布情况相关。线性模型的假设则与之不同，它对原始数据的分布并没有严格要求（这也是为什么人没越来越推崇线性模型），相对应的，而是对残差(residual)有要求。正如下图展示的，通过线性模型拟合出来了比较能估计数据分布和倾向的直线，我们采样的数据点均匀分布在直线两侧，残差就是数据点到直线上的距离的值。

和t检验与ANOVA一样，最简单的检验方法，就是将残差分布可视化。对于线性模型，最重要的基本假设有两方面，一个是残差正态分布，一个是残差的方差齐性。第一个我们可以通过频数分布图（直方图）以及Q-Q图检验，第二个可以直接通过散点图查验，方差齐性的分布图是一个“斑点”式的分布。

# 提取模型拟合的预测值以及残差，并对应好相关的因素
mdl_res <- data.frame(predicted = predict(mdl), residual = residuals(mdl), 
                      NL = lexdec$NativeLanguage, CF = lexdec$Complex)

## 绘制相关图形，本代码只展示最基本的图形，具体美化根据自己需求进行
# 残差直方图
mdlres_1 <- ggplot(mdl_res, aes(x = residual)) + 
  geom_histogram() + 
  geom_density(aes(y = after_stat(density) * 140)) + 
  xlab("Residuals") + 
  ylab("Frequency") + 
  ggtitle("(a) Histogram of residuals")

# Q-Q图
mdlres_2 <- ggplot(mdl_res, aes(sample = residual)) + 
  stat_qq() + 
  stat_qq_line() + 
  xlab("Theoretical quantiles") + 
  ylab("Sample quantiles") + 
  ggtitle("(b) Q-Q plot of residuals")

# 残差散点图，检验方差齐性
mdlres_3 <- ggplot(mdl_res, aes(predicted, residual)) + 
  geom_point() + 
  geom_hline(yintercept = 0, lty = 2) +
  xlab("Fitted values") + 
  ylab("Residuals") + 
  ggtitle("(c) Residual plot")

# 加载patchwork包，将三个图合并到一起
library(patchwork)

mdlres_1 + mdlres_2 + mdlres_3

如果觉得这种方法比较繁琐，也可以选择直接使用performance包里的check_model()函数，它可以直接对模型进行检测并绘制图相关图片。

performance::check_model(mdl)

可能对于初次接触建模的朋友而言，这种“目视”检测方法总觉得不靠谱，不过统计学家们通常更喜欢可视化的方法，是因为它不仅是用来检测基本假设，也可以帮助我们看数据的分布情况，帮助我们决定，是不是需要加入其他因素，或者非线性的内容。不过也可以通过shapiro.test()进行正态分布检测，但是更推荐使用可视化的方法。

3.2.4 无法解决的问题：转向增长曲线分析(GCA)

到目前为止，我们考虑的因变量和自变量之间的关系都是线性的。如果我们可视化的散点图发现，因变量和自变量之间的关系并不是线性的，最典型的比如汉语的声调，对于上声来说，它很明显不是线性的。再比如，语言感知实验里常用到的眼动技术(eye-tracking)，它的数据呈现出来的方式也很明显不是线性的。这时候，就需要考虑稍微复杂的模型。幸运的是，一些情况下我们可以借助混合模型的“扩展版”来解决这个问题，这被称为“增长曲线分析”(Growth Curve Analysis，简称GCA)，它的模型可以写为：

\(y=\beta_0+\beta_1 x+\beta_2 x^2 + ...+\beta_k x^k + \varepsilon\)

GCA比较的是数据的轮廓走势（某种程度上并不涉及所谓的曲线平均值），因为数据之间具有高度相关性 ，它可以帮助我们从整体把握数据的情况。

我们可以看一下示例数据，不难看出声调并不是严格意义上的线性。在一些方言声调研究中，这种声调曲拱就变得更为明显。因此，我们需要借助非线性的方法，帮助我们捕捉非线性特征。如果仅仅使用线性模型，这种特征就有可能被抹消掉了。

GCA所采用的方法，本质上是多个多项式函数相加（记住这一点）。我们稍微回忆一下初中数学里，一个简单的\(y=kx+b\)代表的是一条直线，即线性。当我们多加一个次幂时，比如\(y=ax^2+bx+c\)，我们在原本线性的基础上添加了二次项，这时候图形就变成了一个抛物线。以此类推，我们就可以通过多项式的相加得到不同曲拱的曲线。GCA方法就是利用了这一点：通过添加“高次幂”参数，使之与自变量交互（比如时间，音高总是随时间变化而变化，二者之间有交互效应），从而得到非线性的结果。Mirman (2014)指出GCA分析最多添加至4次幂（如下图所示，Mirman (2014, pp.40），因为更高次幂会导致结果无法解读，也无法将数据实际情况和拟合结果进行很好地对应，同时也可能导致过度拟合。

Mirman (2014)
Mirman (2014)

接下来，我们使用示例数据pitch_data进行GCA分析讲解。示例数据中包括被试（ID）、声调类型（Tone）、焦点类型（Type）、音高数据（P1到P9，即9个采样点）。GCA需要数据为长数据格式，因此我们顺便将数据转变为长数据并赋值给新的变量名。

pitch_long <- pitch_data %>% 
  pivot_longer(P1:P9, names_to = "TimePoint", values_to = "Pitch") %>% 
  mutate(Time = as.numeric(substr(TimePoint, 2, 3)), # 将时间点提取出来并转变为数值型
         across(where(is.character), as.factor)) # 将所有字符型数据转变为因子型

pitch_long
## # A tibble: 576 × 6
##    ID    Tone  Type  TimePoint Pitch  Time
##    <fct> <fct> <fct> <fct>     <dbl> <dbl>
##  1 AD001 T1    Broad P1         19.0     1
##  2 AD001 T1    Broad P2         18.9     2
##  3 AD001 T1    Broad P3         18.9     3
##  4 AD001 T1    Broad P4         18.9     4
##  5 AD001 T1    Broad P5         18.9     5
##  6 AD001 T1    Broad P6         19.0     6
##  7 AD001 T1    Broad P7         19.1     7
##  8 AD001 T1    Broad P8         19.3     8
##  9 AD001 T1    Broad P9         19.3     9
## 10 AD002 T1    Broad P1         21.6     1
## # ℹ 566 more rows

因为GCA需要确定次幂，所以最简单的方法就是，画出图，观察数据分布情况，从而确定数据走势是符合抛物线，还是双折或三折。因为我们的数据是汉语普通话的声调，根据常识我们可以得知，声调走势最多是一个抛物线类型（上声），因此我们就可以将次幂定为二次项。使用poly()函数在时间范围内创建正交多项式二次项，使用unique()函数确保我们不会给相同的时间点赋值不同的二次项正交多项式，之后我们将这些添加到原始数据中，并对其重新命名，令其分别代表一次项和二次项。

# 创建二次幂的正交多项式
t <- poly(unique(pitch_long$Time), 2)

# 将其添加进原始数据，其列名添加ot (代表orthogonal time order)
pitch_long[, paste("ot", 1:2, sep = "")] <- t[pitch_long$Time, 1:2]

这里提到一个概念：正交多项式(orthogonal polynomials)。我们可以这样理解：当面对非线性的数据时，比如时序数据，时间一定是正值，原本与时间结合的次幂项彼此具有高度相关性，当线性项增加时，非线性项也随之增加，这就导致估计的参数也相互依赖，无法独立地进行分析，也就是所谓的共线性：它们会对相同的差异进行分析，一次是线性一次是非线性，这会影响对我们数据比较的解答。此外，共线性还会影响参数估计的稳定性，甚至删除或增加一个数据，都会引起参数的大幅度变化。因此，我们借助正交多项式，可以将参数的相关性抵消，进而去除共线性所引发的问题。

# 没有正交前多项式次幂之间具有高相关性
poly(1:5, 4, raw = TRUE) %>% 
  cor() %>% 
  round(2)
##      1    2    3    4
## 1 1.00 0.98 0.94 0.90
## 2 0.98 1.00 0.99 0.97
## 3 0.94 0.99 1.00 0.99
## 4 0.90 0.97 0.99 1.00

# 使用正交多项式减少相关性
poly(1:5, 4, raw = FALSE) %>% 
  cor() %>% 
  round(2)
##   1 2 3 4
## 1 1 0 0 0
## 2 0 1 0 0
## 3 0 0 1 0
## 4 0 0 0 1
Mirman (2014, pp.47)
Mirman (2014, pp.47)

从上图可以看到，原始多项式中（左），线性和非线性可能都会对相同的差异进行参数评估，导致我们无法剥离不同参数对数据差异性的解释。使用正交多项式后（右），共线性的问题得到解决，可以剥离线性和非线性参数的表现情况和对应含义。

接下来我们便可以建立模型。与LMM一样（因为它本质上是LMM的扩展，也属于混合效应模型），我们可以选择使用逐步回归进行模型比较，但是需要谨慎对待这可能带来的I型错误。我们也可以考理论驱动，确定我们的模型，或者二者的结合。首先，在我们的数据中，研究问题是不同焦点类型的音高轮廓(pitch contour)是否存在差异，此时焦点类型Type就是自变量，声调算是协变量，因为我们只是为了控制声调调型，不过它和焦点类型有明显的交互效应。同时，考虑到不同被试在曲拱和焦点类型影响下音高明显会有不同的随机误差，因此我们考虑该随机误差。总的而言，我们的模型可以建立为如下形式，使用anova()可以查看主效应，借助summary()函数可以看到简单效应。

gca_mdl <- lmer(Pitch ~ (ot1 + ot2) * Type * Tone + (ot1 + ot2 + Type | ID), data = pitch_long)

anova(gca_mdl)
## Type III Analysis of Variance Table with Satterthwaite's method
##               Sum Sq Mean Sq NumDF  DenDF   F value    Pr(>F)    
## ot1             12.0   12.05     1   7.05    8.6037 0.0217413 *  
## ot2             13.5   13.48     1  20.63    9.6226 0.0054761 ** 
## Type            17.0   17.02     1   7.00   12.1539 0.0101797 *  
## Tone          5784.3 1928.10     3 531.00 1376.8042 < 2.2e-16 ***
## ot1:Type         9.7    9.73     1 531.00    6.9504 0.0086251 ** 
## ot2:Type         0.1    0.09     1 531.00    0.0625 0.8026250    
## ot1:Tone      1398.5  466.17     3 531.00  332.8790 < 2.2e-16 ***
## ot2:Tone        25.4    8.45     3 531.00    6.0361 0.0004784 ***
## Type:Tone       84.5   28.17     3 531.00   20.1178 2.311e-12 ***
## ot1:Type:Tone   59.8   19.93     3 531.00   14.2307 6.224e-09 ***
## ot2:Type:Tone    8.4    2.78     3 531.00    1.9885 0.1146822    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(gca_mdl)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Pitch ~ (ot1 + ot2) * Type * Tone + (ot1 + ot2 + Type | ID)
##    Data: pitch_long
## 
## REML criterion at convergence: 1894.7
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.73058 -0.64147  0.03724  0.69926  2.45133 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr             
##  ID       (Intercept) 10.14193 3.1846                    
##           ot1          0.71584 0.8461    0.82            
##           ot2          0.03536 0.1880   -0.30 -0.73      
##           TypeNarrow   2.63634 1.6237   -0.19 -0.09  0.39
##  Residual              1.40042 1.1834                    
## Number of obs: 576, groups:  ID, 8
## 
## Fixed effects:
##                        Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)            23.82103    1.13454   7.16406  20.996 1.07e-07 ***
## ot1                     0.23539    0.51433  38.78809   0.458 0.649742    
## ot2                     0.09133    0.42364 362.08555   0.216 0.829430    
## TypeNarrow              2.67497    0.60700   8.25511   4.407 0.002100 ** 
## ToneT2                 -6.13398    0.19723 530.99747 -31.100  < 2e-16 ***
## ToneT3                 -7.22960    0.19723 530.99747 -36.655  < 2e-16 ***
## ToneT4                 -2.98272    0.19723 530.99747 -15.123  < 2e-16 ***
## ot1:TypeNarrow          0.56193    0.59170 530.99747   0.950 0.342701    
## ot2:TypeNarrow         -0.05398    0.59170 530.99747  -0.091 0.927345    
## ot1:ToneT2              4.57219    0.59170 530.99747   7.727 5.53e-14 ***
## ot1:ToneT3             -2.09307    0.59170 530.99747  -3.537 0.000439 ***
## ot1:ToneT4             -5.77611    0.59170 530.99747  -9.762  < 2e-16 ***
## ot2:ToneT2              0.89502    0.59170 530.99747   1.513 0.130968    
## ot2:ToneT3              0.26983    0.59170 530.99747   0.456 0.648554    
## ot2:ToneT4              0.33418    0.59170 530.99747   0.565 0.572464    
## TypeNarrow:ToneT2      -0.92384    0.27893 530.99747  -3.312 0.000989 ***
## TypeNarrow:ToneT3      -1.77227    0.27893 530.99747  -6.354 4.52e-10 ***
## TypeNarrow:ToneT4       0.11872    0.27893 530.99747   0.426 0.670554    
## ot1:TypeNarrow:ToneT2   0.88869    0.83679 530.99747   1.062 0.288707    
## ot1:TypeNarrow:ToneT3  -2.18229    0.83679 530.99747  -2.608 0.009365 ** 
## ot1:TypeNarrow:ToneT4  -4.07398    0.83679 530.99747  -4.869 1.48e-06 ***
## ot2:TypeNarrow:ToneT2   1.19440    0.83679 530.99747   1.427 0.154063    
## ot2:TypeNarrow:ToneT3   0.15596    0.83679 530.99747   0.186 0.852214    
## ot2:TypeNarrow:ToneT4  -0.83850    0.83679 530.99747  -1.002 0.316774    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation matrix not shown by default, as p = 24 > 12.
## Use print(x, correlation=TRUE)  or
##     vcov(x)        if you need it

GCA模型的解读和LMM的解读基本一致，需要注意我们添加了时间分量，因此我们可以从参数表获得有关截距、斜率、曲拱等的信息，此处不再赘述。同样的，我们需要谨慎对待模型的结构，还是那句话：不要一味追求keep-it-maximal原则，时刻提醒自己研究问题是什么。当出现不同的警告信息(warning message)时，说明已经出现了不合适的拟合结果，需要重新考虑模型结构。

然而，GCA同样也不是万能的。尽管它可以拟合非线性数据，解决了大部分语言学面对的数据，但是正如它的名字一样（多项式），它只是不同线性和非线性项的简单相加，面对一些具有长渐近线的数据，或者“极速增加后趋于平缓”的数据（比如眼动数据），又或者有多个起伏的数据（比如语调），它就不能很好地对数据进行解析了。这时候，我们需要引入其他的方法，帮助我们更好地分析数据差异。

3.3 广义加性（混合）模型（GAM/GAMM）

语言学研究中存在大量非线性的看起来十分“弯弯曲曲”的数据，比如语调，它并不是一个简单的递增或递减的音高轮廓，而是伴随着语言学/副语言学意义出现摇摆不定的状态。又或者元音共振峰数据，它通常在某个阶段有较大变化，其他阶段保持稳定。很多眼动数据也具有相似的特征。这类数据很难通过GCA进行拟合，因为会丢失很多信息，也可能出现过度拟合/过度平滑，导致原本有差异的地方被拟合为没有差异，又或者将本没有差异的地方拟合成有差异。为了解决这个问题，广义加性（混合）模型(Generalized Additive (Mixed-effect) Model，简称GAM/GAMM)成为了不错的选择，也在近几年逐渐成为语言学研究者常用的统计分析方法之一。下面我们将对该方法进行入门介绍，更具体的内容可以参考推荐论文。

GAM/GAMM的拟合需要加载mgcv和istadug包。受益于这些包的开发，现在在GAM中添加随机效应也比以往方便了许多，因此我们不会过多着墨于传统方法，聚焦在如何使用该方法分析语言学数据。不过在开始之前，让我们先理清在GAMM分析中需要用到的核心概念。

3.3.1 核心概念：样条函数(spline)

与GCA不同，GAM/GAMM的实现方式是通过一系列的基础样条函数(spline)对数据进行拟合。什么是样条函数？它其实是分段函数的一种，属于分段多项式函数。数据科学家Katie House举了一个非常形象的例子：当我们使用PowerPoint的曲线工具绘制时，你会发现它并不是“一口气”将曲线绘制完成，而似乎是将一个完整的曲线切分成不同的区域（定义域），将每一段新区域的曲线加到之前的曲线上进行绘制，每个区域都可以使用一个多项式函数进行表达，这每个区域的多项式函数，就称为样条函数。

还记得我们在GCA的时候提到，它是多项式函数的叠加吗？没错，GAM也是多项式函数的叠加，那它们有什么不同呢？不同之处在于，GCA的多项式数目是由研究者认为决定的，也就是说，由我们自己决定究竟是二次项还是三次项。而对于GAMM来说，它是由平滑参数(smoothing parameter)确定。线性模型里，我们通过\(Y=\beta_0+\beta_1 x_1 + \beta_2 x_2 + ... + \beta_n x_n+\varepsilon\)表达\(Y\)与一众自变量\(x\)之间存在直接的线性关系。而在我们要提到的GAM/GAMM里，它考虑了非线性关系，且使用的是一系列基础样条函数叠加形成最终的曲线，因此它的公式可以看作\(g(Y)=\omega_1 f_1(x_1)+\omega_2 f_2(x_2)+...+\omega_n f_n(x_n)+\varepsilon\)，意思是\(Y\)和一众自变量\(x\)之间的关系并非直接的线性关系，而是由“某种/些函数”(some of functions)描述，且针对每个自变量，可能都由与之相关的“某些样条函数”来确定（可能就是线性，也可能是复杂的多次项），且每个自变量对因变量效应的权重(\(\omega\))也不尽相同。这里的某些函数，即为我们下面要提到的平滑参数。函数的数量越多，曲线越摇摆(wiggly)，反之越平滑。

现在，我们对核心概念有了大概的了解。不难看出，GAM/GAMM解决了之前模型的一些痛点，也可以完成对之前模型的拟合。下面我们将借助R语言和示例数据，对其他细节内容进行讲解。

3.3.2 示例：焦点的音高

我们继续使用GCA里所使用的焦点数据。回顾一下我们的pitch_long数据，它已经是长数据格式了，其中包括焦点类型(Type)和调类(Tone)两种。与GCA一致，我们的声调曲线是随时间变化而变化，换句话说，它的曲拱变化与时间具有相关性和交互效应。因此，与GCA分配时间变量的正交多项式次幂一样，GAMM也会通过赋给时间点“一些”样条函数来进行曲线拟合。那么理论上，我们的模型应该是这个模样：

Pitch ~ Type + Tone + Splines_of_Time

有了大概的思路，我们就可以加载mgcv包，使用其中的bam()函数进行拟合。既然我们拟合非线性的GAMM，那么统计学家们就给其中的固定效应起了个新的名字。和线性模型里一样的截距和斜率，被称为参数项(parametric term)，GAMM里允许在参数项基础上添加平滑信息的参数被称为平滑项(smooth term)。那么我们试着在线性模型基础上，添加所谓的平滑项。

# 加载mgcv包
library(mgcv)

# 使用bam函数建模
bam_mdl1 <- bam(Pitch ~ Type * Tone + s(Time), data = pitch_long)
## Error in smooth.construct.tp.smooth.spec(object, dk$data, dk$knots): A term has fewer unique covariate combinations than specified maximum degrees of freedom

让我们先把报错信息放一放，先解释其中的平滑项。s()即添加平滑项的函数，表示“某些函数”。它有多个类型，使用bs = ...可以指定样条函数类型，比如bs = "tp"表示使用薄板样条函数(thin plate regression splines)，bs = "cr"使用的是三次样条回归函数(cubic regression splines)。对于语言学研究而言，不管选择哪个样条函数，其对统计结果的影响较小，因此绝大多数情况下我们只需保持默认值即可。除此之外，随机效应也需要该参数指定，如bs = "re"表示随机截距和随机斜率，bs = "fs"表示的是随机平滑（即每个被试/语料表现出来的曲线平滑度不一样）。下面我们再详细谈。

现在我们来看看错误信息，该错误信息表示我们数据集的量并没有达到平滑函数的默认样条数量值，因此出现了报错。这种情况下，我们就需要稍微改动一下样条数量值。在s()函数中，我们主要是用k = ...来进行定义，其代表的是不同样条函数的焦点多少。如下图所示，表示的是一系列样条函数，它们的交点位置（节点knots，红色的点表示）定义的是GAMM分段样条函数的边界，可以理解为到该边界就到了下一个样条函数区域。并没有一个统一的定义k取值大小的标准，Sóskuthy (2017)建议根据实际数据情况决定究竟使用默认值还是更改。比如元音共振峰的值很少需要k>10的情况，包括现在我们出现的报错信息，也是在说明默认的k值比数据量需要的要多，这时候我们就可以调小一点。当面对语调时，尤其是较大分段的语调，我们就可以使用高一点的值。当然，我们也可以借助mgcv包里的相关函数给我们一些k值的参考，如gam.check()或choose.k()等。

Dorothy L. Andrews, 2021 NAIC Insurance Summit
Dorothy L. Andrews, 2021 NAIC Insurance Summit

为了方便，我们就首先定义k = 9，因为我们的时间点只有9个点而非默认的10个点。

bam_mdl1 <- bam(Pitch ~ Type * Tone + 
                  s(Time, k = 9), 
                data = pitch_long)

注意，这里的s(Time)样板函数定义中，我们没有添加任何因素，因此它类似于回归模型中的y ~ 1，属于参考水平。很明显，我们在这个数据中，存在着不同的声调，因此不同声调的曲拱肯定是不同的。这时候，我们可以通过添加by = ...参数告诉模型，我们这里有受到Tone这一参数而产生的不同的曲线轮廓，因此，我们要在模型中加入这一参数。

bam_mdl2 <- bam(Pitch ~ Type * Tone + 
                  s(Time, k = 9) + 
                  s(Time, by = Tone, k = 9), 
                data = pitch_long)

因为示例数据量比较小，因此现在我们考虑完参数项，接下来就可以考虑随机效应了。上面我们提到，“固定效应”这边存在参数项（截距和斜率）和平滑项，“随机效应”也可以进行这两者的区分，即随机参数项（随机截距和随机斜率）以及随机平滑。和LMM同样的道理，随机效应的添加也需要根据数据的实际情况酌情添加，而不能一味追求最大化结构。

Sóskuthy (2017)
Sóskuthy (2017)

从示例数据来看，不同发音人仅在随机截距和随机斜率上出现了不同表现，他们的曲拱大体相当，因此我们的模型里只需要囊括随机截距和随机斜率即可，那么最终模型就可以建立起来。

bam_mdl3 <- bam(Pitch ~ Type * Tone + 
                  s(Time, k = 9) + 
                  s(Time, by = Tone, k = 9) + 
                  s(Time, ID, bs = "re"), data = pitch_long)

至此，我们就建立完了我们的GAMM模型，受限于数据量，该模型并没有囊括更多的诸如交互效应ti()以及随机平滑，之后随着数据更新，我们会将本文的示例数据进行更新。现在，我们可以先使用summary()函数看一看模型结果。

summary(bam_mdl3)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Pitch ~ Type * Tone + s(Time, k = 9) + s(Time, by = Tone, k = 9) + 
##     s(Time, ID, bs = "re")
## 
## Parametric coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        23.8210     0.9337  25.513  < 2e-16 ***
## TypeNarrow          2.6750     0.3128   8.552  < 2e-16 ***
## ToneT2             -6.1340     0.3128 -19.610  < 2e-16 ***
## ToneT3             -7.2296     0.3128 -23.112  < 2e-16 ***
## ToneT4             -2.9827     0.3128  -9.536  < 2e-16 ***
## TypeNarrow:ToneT2  -0.9238     0.4424  -2.088   0.0372 *  
## TypeNarrow:ToneT3  -1.7723     0.4424  -4.006 7.01e-05 ***
## TypeNarrow:ToneT4   0.1187     0.4424   0.268   0.7885    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                      edf    Ref.df       F  p-value    
## s(Time)        1.000e+00 1.000e+00   0.121    0.728    
## s(Time):ToneT1 3.483e-06 6.965e-06   0.000    0.500    
## s(Time):ToneT2 2.612e+00 3.244e+00  20.743  < 2e-16 ***
## s(Time):ToneT3 1.000e+00 1.000e+00  23.027 2.17e-06 ***
## s(Time):ToneT4 1.000e+00 1.000e+00 138.641  < 2e-16 ***
## s(Time,ID)     6.959e+00 7.000e+00 161.335  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Rank: 55/56
## R-sq.(adj) =  0.857   Deviance explained = 86.2%
## fREML = 1201.9  Scale est. = 3.5224    n = 576

我们可以看到，GAMM同样包括两个参数表，但是取而代之的是，一个参数表代表参数项的评估，即非平滑项的截距和斜率，它的解读方式与LMM是一致的。另一个参数表代表平滑项的评估，包括参考平滑项、不同因素（本例中为声调）的平滑项，最后一个是随机效应。它的解读方式和LMM一样，只不过需要时刻提醒自己，这个参数表是与平滑项相关的。

至此，我们对GAMM的最基本的讲解也结束了，关键的内容是最基本最核心的概念：样条函数。更具体详细的教程，可以参考Sóskuthy (2017)以及其他在推荐材料中涉及到的论文。

4 结语

本文是针对语言学研究者的一份统计学习参考，并假设有一定的R语言基础。虽然我们只介绍了几个常见的统计方法，然而还是需要再次强调，统计方法要选择最合适的，而非最高级的。只有选择合适的方法，才能解决我们的研究问题。

受限于篇幅，本文并没有囊括相关分析(correlation)、广义线性模型(Generalized linear model)、功能主成分分析(Functional principal component analysis)、聚类分析(Cluster analysis)等适用于一些语言学数据的分析方法，读者可根据自己需求寻找相关教程。如有机会和需求，本文会更新这些分析方法。希望本文能使读者对提到的统计分析方法有更清楚和明白的认知，了解了基础概念，才能更自如地使用它们。