R语言 | 回归分析（四）
R语言
语言学与R语言的碰撞
Xu & Yang
PhoneticSan
学习参考
Discovering Statistics Using R
Statistics for Linguistics with R
How to Do Linguistics with R
R in Action
Analyzing Linguistic Data
R Graphics Cookbook
··· ···
Recap
当预测变量是分类型变量时，需要设置虚拟变量以便进行线性回归的拟合。
逻辑回归中，预测变量可以是连续型变量或分类变量，得到的结果是分类变量。根据结果的分类个数，逻辑回归可以被分为二元逻辑回归和多元逻辑回归。glm( )是逻辑回归所用函数。
为了将非线性的结果转变为线性内容，我们可以采用优势比进行转化。
R: The R Project for Statistical Computing
https://www.r-project.org/
RStudio:
https://rstudio.com/
今天是我们最后一期R语言入门系列。从去年11月份到今天，中间也因学业停更了很久，好在坚持了下来。如果您跟着我们一步步走完了整个流程，那么恭喜你成功入门！之后我们还会不定期更新与R有关的内容，比如更丰富的数据可视化，其他统计分析研究相关的内容（如贝叶斯分析）等等。现在开始我们接受R语言相关操作的稿件，如果你有什么想分享的内容，欢迎来稿！
R Project
Linguistics
1
线性混合模型
在语言学的研究中，我们的数据总是互相嵌套，层层递进的。比如，我想考察汉语普通话否定句中否定焦点的位置，那么根据调查，我会考虑影响否定焦点位置的因素。首先，不同否定词是否会对否定焦点有影响，那么我设置了“不”字和“没”字（第一层）。接着，我又考虑有无语句焦点是否对否定焦点身份有影响，那么我设置了带焦点的语句和不带焦点的语句（第二层）。后来，我又考虑否定焦点会不会因为不同的语句焦点位置而发生转移，于是又设置了不同位置的焦点（第三层）。
考察否定焦点问题时进行划分的自变量
另外，与其他众多语言学实验类似，这个实验同样是组内设计(within-subjects design)，即同一批被试接受不同的处理条件。如何更方便地考虑到所有因素，并进行统计检验呢？我们采用的分析方法叫做线性混合效应模型(linear mixed effects model，简称LME)，它本质上是一种多水平模型(multilevel model)。在这个模型中，我们把所有的自变量分为了两种类型：固定效应(fixed effects)和随机效应(random effects)。
fixed effects和random effects下的回归
通俗地讲，固定效应指的是实验中所关注的变量，比如上面例子中我们所感兴趣的否定词、焦点、焦点位置这三个因素都出现在实验中，那么这些就被称为固定效应，即实验中你所控制的自变量。随机效应指的是实验中那些由于从总体中采样而产生的随机误差变量，比如实验被试，这可以被称作随机效应。你或许在一些教材中看到，如果我提到的否定词、焦点、焦点位置三个因素，只是一系列造成否定焦点不同的因素中的三个，那么这些也被称为随机效应。对于这一点我们可以忽略，因为在绝大多数的学术研究中，我们控制的自变量就是固定效应。
固定效应的解释
我们回顾一下线性回归的确定方法，它的公式是Yi = b0+ b1Xi+ εi，其中的i我们说代表的是不同的被试。在这样的简单回归模型中，截距和斜率是固定的，即默认回归模型对所有被试都适用。但是我们也可以将它们设定为随机值，即截距和斜率可能会因随机参数发生变化，这样就会有“随机截距、随机斜率、截距和斜率都是随机”三种情况。在研究中最常见的是第三种，即因随机参数而形成的不同线性归回。比如，随机参数是被试个体，那每个人发音得来的音高肯定会有所不同。
随机参数下的回归模型
现在，我们得到了固定效应和随机效应，那么接下来就是定义LME的数学表达。根据线性回归的表达式，b0 代表着截距，b1代表着斜率，εi代表着随机误差。加入随机效应后，截距和斜率都有一个随机参数参与进去，比如记作u，LME会根据随机参数给每一个水平拟合一个适合的模型，比如给第一个水平的斜率和截距添加随机参数，那么表达式就成为：Yij = (b0j + u0j) + (b1 + u1j)X1 + b2X2 +… + bnXn + εi。以此类推。看到公式不要紧张，这些在R中都自动运行，这里只是帮你理解。
线性回归模型的数学比较
最后，我们谈一谈LME中的εi。由于随机效应的影响，我们的误差也被分为了两部分。一部分是随机误差(random error)，它是受到我们可以解释的随机参数的影响（比如不同的被试），一部分是残差误差(residual error)，它在模型的解释范围外。接下来，我们介绍LME的R实现方式。
R Project
Linguistics
2
增长曲线分析
到目前为止，我们介绍到的内容都是线性回归分析的方法。对于语音学研究而言，比如F0的走势往往不是线性的。换句话说，F0走势与时间组成的函数关系不是一次函数，因此单纯使用线性回归来分析基频曲线就显得稍有不足，拟合程度比较差。另外，之前我们也提到过，如果使用传统的t检验或方差分析，它们都是点对点的比较，这就很容易出现F0采样的几个点中，个别点出现了显著性差异，个别点没有，造成结果难以解读。如何既避免传统检验带来的问题，又可以使用类似于线性回归的方法？这就是我们要介绍的增长曲线分析(Growth Curve Analysis，简称GCA)
GCA可以分析非线性数据
在开始之前，我们首先对GCA的定义和适用范围进行介绍。GCA本质上也是多水平回归(multilevel regression)的分析方法，分析非线性数据效果极佳，尤其是那些时序数据(time course data)或纵向数据(longitudinal data)。F0就是一个典型的可以借助GCA进行分析的数据，首先它是非线性变化，其次随着时间历程而发生变动，还有其他行为实验数据，如眼动、EEG等。GCA的主要优点在于，它可以同时分析不同实验变量之间以及个体之间的差异。GCA与其他多水平回归一样，需要确定固定效应和随机效应。
因为GCA可以用来分析非线性数据，那么很显而易见地，在GCA中除了固定效应和随机效应，还需要了解到要分析数据的形状，即要保证拟合出来的模型的函数图像，要较好地拟合数据的趋势。GCA模型使用高阶多项式进行拟合，通俗地说，它使用的不再是一个简单的Yi = b0+ b1Xi+ εi，而是根据曲线的走势添加高阶项。比如使用GCA拟合汉语上声调，它是一个类似于抛物线的二次函数，那么GCA就会添加一个二阶项，成为Yij = b0i+ b1iTimej+ b2iTimej2+ εij，其中的i和j分别表示第i个被试在时间j上的数值。如果是个双折调，那么我们只需再添加三阶项即可。通俗地讲，它实际上就是我们数学中的二次、三次、四次...n次函数。
当高阶项越多，拟合模型的形状也会变化
如何确定我需要的高阶项，这个需要考虑很多方面。我们建议在使用GCA进行分析的过程中，首先从最基本的简单线性回归做起（为什么这么讲，后面的R实现会提及），然后一点点添加实验中控制和观察的变量/效应。在语言研究中，我们一般不提倡超过四阶项。在确定之后，检查模型的拟合情况，选择最优模型即可。
R Project
Linguistics
3
R语言中的实现
在前两节提到的LME和GCA，它们本质都是多水平回归，因此在实现方面，我们以较复杂的GCA为主进行实例讲解，LME以假设实验为主，二者使用的函数均为lme4包中的lmer( )函数。
首先解读LME的实现方式。假设某实验测定了变量A（数据中记作Condition）对阅读时长RT的影响，那么根据LME的概念，这里的Contion即为固定效应，因此记作RT ~ Condition。下面我们考虑进随机效应，比如实验中被试(Subjects)和实验条目(Item)，在lmer( )函数中使用“｜”将随机效应进行区分。因此，我们LME可以这样写。
# 加载lme4包library(lme4)# LME建模ReadingTime <- lmer(RT ~ Condition + (1 | Subjects) + (1 | Item) + (0+Condition | Subjects)
其中(1 | Subjects)和(1 | Item)表示我们按Subjects和Item添加随机截距，(0+Condition | Subjects)表示为了考虑不同Subjects对Condition的反馈方式可能有所不同，因此按照Subjects给每一个预测变量Condition添加一个随机斜率。总体而言，在LME中，如果添加随机截距，那么是(1 | random effects)；如果添加随机斜率，那么是(0 + fixed effects | random effects)。
lmer( )函数的构成
对于GCA分析，为了方便联系，我们使用Mirman(2008)的数据，点击阅读原文可以跳转至相关网页，最下方提供下载。下面我们把所有代码展示出来，然后逐一讲解。
library(ggplot2)library(reshape2)# 加载lme4需要Matrix包library(Matrix)library(lme4)wd <- read.table(file.choose(), header=T)# 将Subject变为factorwd$Subject <- as.factor(wd$Subject)# summary(wd)# 设置二阶多项式t <- poly(unique(wd$Block), 2)# 创建变量ot1和ot2并与Block对齐wd[, paste("ot", 1:2, sep="")] <- t[wd$Block, 1:2]# summary(wd)#建立GCA模型，从简单到复杂m.base <- lmer(Accuracy ~ (ot1+ot2) + (ot1+ot2 | Subject), data=wd, REML=FALSE)m.0 <- lmer(Accuracy ~ (ot1+ot2) + TP + (ot1+ot2 | Subject), data=wd, REML=FALSE)m.1 <- lmer(Accuracy ~ (ot1+ot2) + TP + ot1:TP + (ot1+ot2 | Subject), data=wd, REML=FALSE)m.2 <- lmer(Accuracy ~ (ot1+ot2)*TP + (ot1+ot2 | Subject), data=wd, REML=FALSE)# GCA模型比较，选择最优模型anova(m.base, m.0, m.1, m.2)# 输出相关系数，查阅p值coefs=data.frame(coef(summary(m.2)))coefs$p <- 2 * (1 - pnorm(abs(coefs$t.value)))coefs
首先，你遇到的第一个“不明所以”的语句是poly( )函数，这一句意思是在block范围下创建一个二阶多项式。换一个例子，比如，我们分析上声调，那么就是在上声调的完整时间范围内创建二阶多项式。通俗地讲，它就像是定义函数的定义域。但是仅仅定义了还不足够，我们的目的是进行数据分析，二阶多项式的值必须要与要检验的数据关联。因此，下一句就是进行连接。如何让数据明白谁对应谁呢？我们创建了"ot"这个变量（我用orthogonal time的缩写代表，当然你也可以创建名为“a”的变量），因为是二阶多项式，因此我们的ot变量也应当是2。再观察调整后的数据，可以看到我们已经创建好所需要的内容了。
创建好相关变量
接下来就是进行建模分析，因为我们考察的是在Block这个时序上Accuracy的变化情况，因此在lmer( )函数中的开头要声明Accuracy是随Block变化的，又因为我们已经设定好相关的ot变量，这些ot变量控制着Accuracy随Block的走势，因此记作Accuracy ~ (ot1+ot2)。其实从本质上看，它的建模方式与任何其他回归分析的建模方式并不较大的不同。需要额外提醒的是，除了我们之前见到的➕表示更多的预测变量，：表示交互效应，这里出现了*代表的是TP这个预测变量在所有时间上的影响。最后通过anova( )函数选择最优模型。
选择最优模型
在选择最优模型的时候，我们要同时考虑logLik和p值，后者毋庸置疑是看那个模型与null model相比具有显著性差异，侧面说明因素是否对因变量有影响。logLik是对数似然比值(log-likelihood ratio)，检测拟合优度，这个数值越大说明拟合越好。综合这两方面来看，m.2的拟合最优，因此我们输出这个函数中的系数。最后三步则是取出模型中t值相对应的p值，获得观测值的概率。
输出结果
其中前三行是参考变量(reference)，默认情况下参考的是你变量设置中的第一个，比如该数据中TP的第一个变量是High，因此它被设置为了参考变量。紧随其后的TPLow是它的拟合回归模型与High的比较结果，接着的ot1和ot2代表的是斜率和曲度，同样是与High进行对比的结果。
LME和GCA的研究方法在语言学研究中用途越来越广，也是近几年被广泛运用的方法。一开始对于回归分析可能会觉得有些难度，但是随着不断练习，你也可以更好地使用R为你的数据分析服务。
*本节GCA分析参考Mirman (2014) Growth Curve Analysis and Visualization Using R
—END—
排版：Xu & Yang