學習參考
Discovering Statistics Using R
Statistics for Linguistics with R
How to Do Linguistics with R
R in Action
Analyzing Linguistic Data
R Graphics Cookbook
··· ···
當預測變量是分類型變量時，需要設置虛擬變量以便進行線性回歸的擬合。
邏輯回歸中，預測變量可以是連續型變量或分類變量，得到的結果是分類變量。根據結果的分類個數，邏輯回歸可以被分為二元邏輯回歸和多元邏輯回歸。glm( )是邏輯回歸所用函數。
為了將非線性的結果轉變為線性內容，我們可以採用優勢比進行轉化。
R: The R Project for Statistical Computing
https://www.r-project.org/
RStudio:
https://rstudio.com/
今天是我們最後一期R語言入門系列。從去年11月份到今天，中間也因學業停更了很久，好在堅持了下來。如果您跟著我們一步步走完了整個流程，那麼恭喜你成功入門！之後我們還會不定期更新與R有關的內容，比如更豐富的數據可視化，其他統計分析研究相關的內容（如貝葉斯分析）等等。現在開始我們接受R語言相關操作的稿件，如果你有什麼想分享的內容，歡迎來稿！
在語言學的研究中，我們的數據總是互相嵌套，層層遞進的。比如，我想考察漢語普通話否定句中否定焦點的位置，那麼根據調查，我會考慮影響否定焦點位置的因素。首先，不同否定詞是否會對否定焦點有影響，那麼我設置了「不」字和「沒」字（第一層）。接著，我又考慮有無語句焦點是否對否定焦點身份有影響，那麼我設置了帶焦點的語句和不帶焦點的語句（第二層）。後來，我又考慮否定焦點會不會因為不同的語句焦點位置而發生轉移，於是又設置了不同位置的焦點（第三層）。
另外，與其他眾多語言學實驗類似，這個實驗同樣是組內設計(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實現方式。
到目前為止，我們介紹到的內容都是線性回歸分析的方法。對於語音學研究而言，比如F0的走勢往往不是線性的。換句話說，F0走勢與時間組成的函數關係不是一次函數，因此單純使用線性回歸來分析基頻曲線就顯得稍有不足，擬合程度比較差。另外，之前我們也提到過，如果使用傳統的t檢驗或方差分析，它們都是點對點的比較，這就很容易出現F0採樣的幾個點中，個別點出現了顯著性差異，個別點沒有，造成結果難以解讀。如何既避免傳統檢驗帶來的問題，又可以使用類似於線性回歸的方法？這就是我們要介紹的增長曲線分析(Growth Curve Analysis，簡稱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實現會提及），然後一點點添加實驗中控制和觀察的變量/效應。在語言研究中，我們一般不提倡超過四階項。在確定之後，檢查模型的擬合情況，選擇最優模型即可。
在前兩節提到的LME和GCA，它們本質都是多水平回歸，因此在實現方面，我們以較複雜的GCA為主進行實例講解，LME以假設實驗為主，二者使用的函數均為lme4包中的lmer( )函數。
首先解讀LME的實現方式。假設某實驗測定了變量A（數據中記作Condition）對閱讀時長RT的影響，那麼根據LME的概念，這裡的Contion即為固定效應，因此記作RT ~ Condition。下面我們考慮進隨機效應，比如實驗中被試(Subjects)和實驗條目(Item)，在lmer( )函數中使用「｜」將隨機效應進行區分。因此，我們LME可以這樣寫。
library(lme4)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)。
對於GCA分析，為了方便練習，我們使用Mirman(2008)的數據，點擊閱讀原文可以跳轉至相關網頁，最下方提供下載。下面我們把所有代碼展示出來，然後逐一講解。
library(ggplot2)library(reshape2)library(Matrix)library(lme4)
wd <- read.table(file.choose(), header=T)wd$Subject <- as.factor(wd$Subject)t <- poly(unique(wd$Block), 2)wd[, paste("ot", 1:2, sep="")] <- t[wd$Block, 1:2]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)
anova(m.base, m.0, m.1, m.2)
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