實驗語言學線性建模：從入門到跳樓（一）

@lordary ：地鐵上看到一個掛著百度工牌的妹子，手裡拿著一本《餐廳經營：從入門到精通》。

@Chris Xia ：哈哈哈哈哈這是要跳槽嗎？

@lordary ：開飯館挺好啊，不過開飯館也不能看這種書啊！從入門到精通系列難道不就是專坑程序員的么！

@Chris Xia ：竟然不是「從入門到放棄」。

@灰兒 ：「從買遊標卡尺到放棄」。

@katchine ：《 ______：從入門到轉行》，請填空。

@Rockie Yang ：《生狗： 21 天從入門到轉行》。

那麼好的我們今天就來討論一下為什麼要在實驗語言學的數據處理里使用線性建模，以及怎麼搞這一套。是的我試圖拿這套方法來說服我導師來著，而且貌似已經成功了。今天講的內容主要是有關於入門的方式，下一次我們可以具體談談上樓和跳樓。

那麼目前我們有這樣一個實驗設計，怎麼說呢，就算是我的實驗設計的精簡版吧，真的已經很精簡了。我找了兩組人，每組 15 個，分別叫他們「高」和「低」（在代碼里實現為 High 和 Low 兩組），每一組給了他們 75 個句子，讓他們在一個 11 點的里克特量表上打分，判斷句子到底能不能接受。然而在這 75 個句子里，只有一部分句子包括了我要觀測的詞，我把這些詞分為三類，就叫它 B 、 S 和 T 吧。然後，這三組詞里的每一個詞，我都按照某種規則造出兩個句子來，每一個詞就有了兩種狀態，就叫它「L」和「M」狀態吧。那麼我現在就有了一個 2 × 3 × 2 的設計，第一個 2 指的是被試組間差距「高」和「低」，中間的 3 指的是被試組內的詞類差「B」、「S」和「T」，最後一個 2 指的是每個詞的狀態差「L」和「M」。它們正交會出現十二種情況，簡單的排列組合，我就不教了。

那麼好的一般人見到這種數據會怎麼想呢，基本就是 T-test 或者 ANOVA 永遠愛你了。現在的語言學期刊基本也是 T-test 或者 ANOVA 永遠愛你，沒個 ANOVA 似乎都不好意思發出來。除此之外，大概十篇裡面有一篇可能用了 MANOVA ，二十篇裡面能有一篇需要用 ANCOVA ，不過這兩者在實驗設計上本身就並不常考慮，所以在統計分析上也不一定需要包括。

等到我自己開始做實驗的時候，特別是我拿到基礎數據以後，就開始想一個問題：為什麼很多基於正確度判斷的研究在統計時還是採用 ANOVA 的？要知道， T-test 或者 ANOVA ，或者其他的所謂 parametric test ，都有一個大前提：需要處理的數據必須是正太，啊不，符合正態分布。但是正確度判斷的一般狀況是「非對即錯」，中間態並不是很常見，特別是 5 點以上的里克特量表或者 magnitude estimation ，所以我們得到的數據有可能集中在兩端，而中間並不會鼓起來一部分。如果你的數據不幸跟我這裡一樣，大家更偏向於集中在 0/1 或者 9/10 ——或者說，這就是正確度判斷的常態——基本連正態的影子都見不到，就不能用 T-test 或者 ANOVA 了。

你們見過誰家正態分布長這樣？

另，本文中所有圖像均由 ggplot2 友情贊助。

library(ggplot2) n

那麼我們可以用什麼呢，很簡單，既然 parametric 用不了，那就用 non-parametric test 吧。常見的 non-parametric test 的機制也很簡單：如果說 parametric test 的運算內容是數值本身，那麼 non-parametric test 的運算內容一般是數值在整個數據組合里的排名。

但是 non-parametric test 有時候也會變得很糾結；我們繼續拿我這個實驗來舉例。在我的實驗設計里，所有的被試在判斷句子的時候都還需要提供句子的修改，如果在含有關鍵詞的句子里他們「改的不是地方」，那麼相應的判斷數據點就要被拋棄（我想如果之前參加過我的實驗大概就已經意識到這是一個什麼樣的設計了）。人都是會有疏忽的，改了不應該改的地方很常見，所以最後很難拿到一套完整的數據，基本上都是會刪掉一兩個壞點的。那麼我如果現在要比較這兩組被試對於 A 組詞和 B 組詞的判斷有沒有區別，應該怎麼辦？

那麼如果是常規的統計手段的話這種設計一般就會採用 two-way ANOVA 來進行分析，被試的小組作為「組間變數」，而兩組詞作為「組內變數」（因為所有的被試都或多或少給出了這兩組詞的判斷）。而我們之前也看到了，就接受度測試來說，幾乎無法獲得正態分布的數據，我們就必須採用 non-parametric test。用來替代 two-way ANOVA 的方法，一般首選的是 Friedman Test ，你去查所有的統計學書籍，最後的第一建議都會是這個方法。但是 Friedman 有一個特別的要求：所有的樣本必須是平衡的，也就是說兩組被試給出的判斷、 A 組詞和 B 組詞給出的判斷，都必須一樣多，否則無法演算。那麼我現在就可以懵逼了，這是第一個懵逼。

那麼第二個懵逼之處，就是傳統的 parametric test 和 non-parametric test 對於連續數值自變數的處理程度很差。不管是 T-test 還是 ANOVA 還是它們的替代品，在統計開始之前都需要預設有限個分組，那麼如果有一個變數是連續變數（比如我要分析這兩組人各自學習英語的時間長度和他們對於句子的判斷程度），之前所說的測試就都不好使了。當然我們可以用相關性測定，但是相關性測定通常是一對一的，也就是說一個因變數只能接受一個自變數的分析。如果我要同時看被試分組、每個人學習英語的時間、被試的性別這三者對於句子判斷的影響，以及它們之間有沒有相互影響（比如說學英語時間長的高組女生會比高組男生表現好而時間短的沒有影響，而低組無論男女只要學習英語時間長就會有影響），那麼不管是 ANOVA 、 T-test 還是相關性測定都會哭給你看。

那麼第三個懵逼之處就是「人」的區別了。俗話說得好，「人和人之間的差距有時候比人和狗都大」（好像哪裡不太對），這裡之所以用接受度評價來作為例子，就是因為「人」的因素有些時候會特別麻煩。比如說，用一個 0 到 10 的量表以「打分」的方式來評價接受度，有些人可能會用一整條量表由 0 到 10 來衡量，有些人可能只會用 0 和 10 ，還有些人可能用 5 到 10 乃至 7 到 9 ，每個人都有自己的選擇範圍。在這種情況下，由於每個被試對量表的解釋和用法不同，我們無法直接把每一個句子的得分相加，否則結果反映出來的並不是被試們對句子的接受程度（一般來說，「好」和「壞」之間的區別會縮短，大家可以想想是為什麼）。那麼這就是實驗中出現的不同被試引起的噪音，而我們怎麼排除掉這一部分噪音呢？之前我曾經採用過 z-score standardisation 的方法，把絕對的分數轉換成相對的「高於 / 低於個人平均值 X 個標準差」，從而顯示到底是接受還是不接受。但是這樣的方法很容易受到極值和壞點出現比率的影響，所以有沒有什麼方法可以相對可靠地處理這些噪音呢？

俗話說得好，人和人的差距……

那麼第四個懵逼之處，一般是很少有人想到的一個問題了，然而我們——我這裡說的是實驗語言學研究——必須直視這個問題，那就是「 language as a fixed factor 」。比如說，我的三組詞表裡分別包括了八個詞，針對這二十四個詞所做的實驗，用普通的組內比對和組間比對方法，真的就可以直接推廣到所有的 B 類詞、 S 類詞和 T 類詞上嗎？並不能。再比如，萬一 B 類詞里有一個比較少見，而給出的其他類詞大家都見過，那麼這樣的差別會不會影響到我們對於結果的分析和解釋？——當然這是實驗設計里本身就需要避免的事情，但是有些時候在控制了一些變數之後就會發現剩下的變數已經由不得你來控制了。像這樣的實驗選材中的可行性問題、從小樣本推廣到大規模的分析，這些能不能依靠統計方法來彌補？很顯然我們之前所說的方法並不能有效地解決這個問題，那麼該怎麼辦？

基本上，為了迴避這四個懵逼，我開始改用線性建模的方法來描述數據。（回復里有統計學的朋友指出，之前提到的 parametric & non-parametric test 也都是特定情況下的線性建模，在這裡我遵循語言學論文的習慣描述，把狹義的 lm 和 lmer 下的建模方式稱為「線性建模」。）那麼什麼是線性模型？我反正高中選的文科，大學是語言專業，線性規劃什麼的都沒有系統學過，這裡就不多說了。用個形象點的說法，一個散點圖，我們要看橫軸和縱軸之間有沒有線性的影響，就試著用一條直線來模擬散點的分布。這個方程的建構就是線性模型的基礎版，我們不預設任何組內差異、組間差異，單純用一條線來覆蓋所有的點，就是線性模型。

一個線性的示意圖。

那麼，即使是把最簡單的線性模型和 ANOVA 相比，前者又有什麼好處呢？首先，在數據點夠多的情況下，單個的壞點可以直接拋棄而不明顯影響模型的精確度和可用度；其次，可以輕鬆計算多個自變因素對因變數的影響，特別是連續的自變因素的影響；其三，可以有效解決數據不正太，啊不，不正態的問題，可以突破正態分布對於統計演算法選擇的限制。即使是最簡單的不包含任何因素的線性模型，也可以輕鬆完成這一點。

但是我們還有問題沒解決，那就是第三個和第四個懵逼之處。這時候我們就要祭出來更強大的模型功能，也就是混合要素模型（mixed effect model）。這裡我們採用的語言是 R ，調用 lme4 包。

library(lme4) n

那麼我們剛才說了自變數和因變數的線性關係（不過首先你要確定你的因變數是受自變數線性影響的，當然在語言學實驗設計里很少能出現不線性影響的情況，不過我的確有過，回頭可以展示一下），實際上在實際觀測里，還會有一個玩意兒，那就是噪音。噪音是什麼樣的東西呢？就是我說的第三和第四個懵逼。被試者之間的差異也好，測試材料之間的差異也好，都是噪音的來源，我們在進行統計的時候需要過濾掉噪音，盡量把基線拉平。而過濾噪音的方法，就是添加隨機變數。

比如說，我要統計高組和低組對於所有包含要檢測的詞的句子（也就是所謂的 Critical Sentence ，其中 Cri = Y）的整體評分，那麼最初的模型是這樣的：

datCri <- dat[(dat$Cri == "Y") ,] nlinearModel1 = lm(Score ~ Group, n data=datCri)nnsummary(linearModel1)n

這樣得出來的結果是 High 組作為標準組的平均分數（顯示為 Intercept），以及 High 組和 Low 組平均分的差，可以看到 Low 組普遍評分較高。為什麼是 High 組在前呢？因為按照字母順序排序， H 比 L 在前啊。

Call:nlm(formula = Score ~ Group, data = datCri)nnResiduals:n Min 1Q Median 3Q Max n-8.3191 -1.3191 0.6809 1.6809 3.1100 nnCoefficients:n Estimate Std. Error t value Pr(>|t|) n(Intercept) 6.88997 0.09759 70.60 <2e-16 ***nGroupLow 1.42918 0.13865 10.31 <2e-16 ***n---nSignif. codes: 0 『***』 0.001 『**』 0.01 『*』 0.05 『.』 0.1 『 』 1nnResidual standard error: 2.615 on 1421 degrees of freedomn (17 observations deleted due to missingness)nMultiple R-squared: 0.06957, Adjusted R-squared: 0.06891 nF-statistic: 106.2 on 1 and 1421 DF, p-value: < 2.2e-16n

但是 High 組和 Low 組都有那麼十幾號人，肯定會有人搞出來要不是 0 到 10 要不是 5 到 7 的幺蛾子，所以我們要抹平這個差異，就像下面這樣：

mixedModel1 = lmer(Score ~ Group +n (1|PartiNo),n data=datCri， REML=FALSE)nnsummary(mixedModel1)n

這樣我們把每個被試編號的數值分別看待，讓演算法把它們之間的基線拉平（如果你是一個語言學或者心理學乃至社會學從業人士，你並不一定需要知道這個基線是怎麼拉平的），而這個得出來的結果，則是考慮到所有被試的個體差異之後，得出來的結果。被試的個體差異，顯示為 Random effects 裡面的 PartiNo 一項，在這裡就是所謂的噪音，也就是我們需要消除的隨機因素。可以看到，消除了被試差異以後， High 和 Low 兩組的差異雖然存在，但是沒有那麼明顯了。

Linear mixed model fit by maximum likelihood t-tests use Satterthwaiten approximations to degrees of freedom [lmerMod]nFormula: Score ~ Group + (1 | PartiNo)n Data: datCrinn AIC BIC logLik deviance df.resid n 6490.1 6511.2 -3241.1 6482.1 1419 nnScaled residuals: n Min 1Q Median 3Q Max n-3.8032 -0.4074 0.1241 0.7168 1.9413 nnRandom effects:n Groups Name Variance Std.Dev.n PartiNo (Intercept) 1.571 1.253 n Residual 5.260 2.293 nNumber of obs: 1423, groups: PartiNo, 30nnFixed effects:n Estimate Std. Error df t value Pr(>|t|) n(Intercept) 6.8905 0.3348 29.9580 20.583 < 2e-16 ***nGroupLow 1.4337 0.4736 29.9950 3.027 0.00503 ** n---nSignif. codes: 0 『***』 0.001 『**』 0.01 『*』 0.05 『.』 0.1 『 』 1nnCorrelation of Fixed Effects:n (Intr)nGroupLow -0.707n

那麼我們換個思路。比如說我們現在有了所有包含 T 類詞的句子的評分，我們想知道在不區分高低兩組時，甲和乙兩個狀態會不會對 T 類詞的判斷產生影響。那麼我們就可以先這麼做：

mixedModelT = lmer(Score ~ Conds +n (1|PartiNo) +n (1|Item),n data=datT, REML=FALSE)nnsummary(mixedModelT)n

沒錯你可能要問，卧槽這個 Item 哪裡來的？這就是剛才之前我們所說的第四個懵逼。辭彙編號作為噪音的原因是，我們希望這個差異不僅僅是針對所有參與測試的 T 類詞，也希望能夠擴展到所有的 T 類詞上。抹平由詞個例帶來的組內差異，也就意味著我們在進行結論推廣的時候可以不用局限於「這些詞形成的差異」，從而更有底氣地把自己的結論推廣出去。

Linear mixed model fit by maximum likelihood t-tests use Satterthwaiten approximations to degrees of freedom [lmerMod]nFormula: Score ~ Conds + (1 | PartiNo) + (1 | Item)n Data: datTnn AIC BIC logLik deviance df.resid n 2155.7 2176.4 -1072.8 2145.7 467 nnScaled residuals: n Min 1Q Median 3Q Max n-3.7300 -0.3997 0.1491 0.6146 2.1606 nnRandom effects:n Groups Name Variance Std.Dev.n PartiNo (Intercept) 2.1559 1.4683 n Item (Intercept) 0.2018 0.4492 n Residual 4.7370 2.1765 nNumber of obs: 472, groups: PartiNo, 30; Item, 8nnFixed effects:n Estimate Std. Error df t value Pr(>|t|) n(Intercept) 8.2644 0.3422 34.7000 24.152 < 2e-16 ***nCondsM -1.4277 0.2004 435.0000 -7.123 4.42e-12 ***n---nSignif. codes: 0 『***』 0.001 『**』 0.01 『*』 0.05 『.』 0.1 『 』 1nnCorrelation of Fixed Effects:n (Intr)nCondsM -0.292n

但是這還沒完。你覺得這個完了嗎？並沒有。現在我告訴你， T 類詞有這樣一個性質（本來不應該是我告訴你的，而是你在設計的時候已經想好了的）：對於一些人來說，T 組詞的 M 狀態比 L 狀態更難以接受。而同樣是在 M 狀態下，有一半的 T 類詞似乎比其他的詞更難以被接受。很好，那麼如果這個狀態成立的話，我們的情況又更複雜了一些：每個人都閱讀了 T 類詞的所有句子，而 T 類詞的句子又都有 L 和 M 兩種狀態，那我們該怎麼排除掉這些由狀態變化導致的組內波動？

這時候我們就需要考慮到排除自變數對隨機變數的影響，引入一個隨機變率（random slope），就像這樣：

mixedModelT2 = lmer(Score ~ Conds +n (1 + Conds|PartiNo) +n (1 + Conds|Item),n data=datT, REML=FALSE)nnsummary(mixedModelT2)n

這樣我們考慮到了不同被試在不同狀態下的接受區別，也考慮到了不同辭彙在不同狀態下的接受區別，排除了被試之間的大差異和辭彙之間的大差異，兩個狀態間是否有評價差異就變得更不容易受干擾了，從而也就排除了辭彙差別帶來的 Type I Error （假陽性）。不過，即使排除了這些假陽性，我們發現， M 狀況下的 T 組詞還是普遍更不容易被接受。

Linear mixed model fit by maximum likelihood t-tests use Satterthwaiten approximations to degrees of freedom [lmerMod]nFormula: Score ~ Conds + (1 + Conds | PartiNo) + (1 + Conds | Item)n Data: datTnn AIC BIC logLik deviance df.resid n 2151.1 2188.5 -1066.6 2133.1 463 nnScaled residuals: n Min 1Q Median 3Q Max n-3.6266 -0.3709 0.1450 0.5928 2.4294 nnRandom effects:n Groups Name Variance Std.Dev. Corr n PartiNo (Intercept) 1.6121 1.2697 n CondsM 0.4377 0.6616 0.52 n Item (Intercept) 0.1071 0.3273 n CondsM 0.5462 0.7391 -0.17n Residual 4.4767 2.1158 nNumber of obs: 472, groups: PartiNo, 30; Item, 8nnFixed effects:n Estimate Std. Error df t value Pr(>|t|) n(Intercept) 8.2615 0.2933 23.4580 28.166 < 2e-16 ***nCondsM -1.4230 0.3476 9.5420 -4.093 0.00239 ** n---nSignif. codes: 0 『***』 0.001 『**』 0.01 『*』 0.05 『.』 0.1 『 』 1nnCorrelation of Fixed Effects:n (Intr)nCondsM -0.092n

那麼好的今天我們就先到這裡，下一期「上樓」部分則會針對以下內容：

• 如何在模型里設置多個自變數？

• 如何判定一個因素可以作為隨機變率出現？

• 如何篩選和比對多個模型？

• 如何判定一個模型里的自變數對因變數有顯著影響？

• 如何在寫作中彙報我所用的模型和我的分析結果？

最後一期的「跳樓」則是以下內容：

• 我不想要隨機變率，可以嗎？（簡單回答：不可以。）

• 為毛不可以不要隨機變率？（好的，在特定情況下你的確可以不要。）

• 我的模型無法擬合，我除了跳樓之外還有什麼辦法？
• 如果我要進行線性建模，需要什麼硬體設施？（相信我，你也會想要雲計算的。）

• 我的因變數數據不受自變數線性影響，我會不會死？

祝大家跳樓愉快。

Chris

26/02/2017

p.s. 上一篇文章在我挪到另一個專欄的時候被我誤操作刪掉了，我正在聯繫站方詢問能不能恢復。我還是那點：我是有可能講得不好，但麻煩您像「對人說話」那樣說出來，畢竟我在講的時候默認您是個人。

p.s.p. 專欄這種「我要把文章轉移到另一個專欄里必須先打開專欄投稿限制再投稿」的設計是不是有點……過於麻煩了……

p.s.3. 希望移動客戶端可以選擇專欄發表。如果站方允許的話，我會非常感激。

p.s.4. 在這段時間忙完以後，我會慢慢地寫一點有關語言學實驗類的方法論問題，包括判斷測試（選擇模式設計到執行）、自由訪問、自速閱讀、快速閱讀、回憶性實驗、方言調查錄音等，統計和實驗設計也會包括在內（包括今天的組內組間設置）。眼動、事件電位和功能核磁因為我經驗不足，所以就先不寫了，可能會涉及到一些理論問題，如果有需要的話會慢慢寫。大家如果有特別想知道的方法問題，請給我留言，我會在寫完線性建模之後優先考慮。我再也不做 Live 了，瘋了。