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
其中(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 # 将Subject變為factorwd$Subject as.factor(wd$Subject)# summary(wd)# 設定二階多項式t 2)# 建立變量ot1和ot2并與Block對齊wd[, paste("ot", 1:2, sep="")] 1:# summary(wd)#建立GCA模型，從簡單到複雜m.base FALSE)m.0 FALSE)m.1 FALSE)m.2 FALSE)# GCA模型比較，選擇最優模型anova(m.base, m.0, m.1, m.2)# 輸出相關系數，查閱p值coefs=data.frame(coef(summary(m.2)))coefs$p 2 * (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