文章資訊
標題：A Very Basic Tutorial for Performing Linear Mixed Effects Analyses (Tutorial 2)
作者：Winter (加州大學美熹德分校，美國加利福尼亞美熹德)
年份：2013 年
美國加利福尼亞美熹德
摘要
這篇教學文章用簡單易懂的文字告訴我們什麼是線性混和效應模型 (linear mixed-effects model)、其中包含的元素和意義是什麼，以及它的用途，並實際以 R 語言操作，具體呈現教學內容。
※ 可以善用目錄跳到想了解的地方喔。
線性模型？
假設有個實驗是測量「attitude (態度)：pol (禮貌) / inf (非正式)」 會不會影響說話的「frequency (音頻)」，線性模型的呈現：
frequency ~ attitude + ε
attitude 是固定因子 (fixed factor)、ε 是隨機因子 (random factor，包含各種除了固定因子外，可能影響音頻結果的因素)。
如果除了之外，還想要多看 gender (性別) 會不會影響 frequency，模型呈現：
frequency ~ attitude + gender + ε
線性回歸概念請參考：線性回歸幼幼班
線性混和效應模型？
線性模型包含：
- 自變數 — 例如：態度 (attitude)、性別 (gender)
- 應變數 — 例如：音頻 (freuqency)
- 隨機因子 — 也就是 ε
自變數是我們操控的，應變數根據自變數而有不同變化，但隨機因子不受我們控制，我們不喜歡，沒安全感。
所以線性混和效應模型 (linear mixed-effects model) 為此而生，它混合了固定因子 (fixed factor) 和能讓我們指定的隨機因子 (random factor)，滿足掌權感。因為混和了固定因子和隨機因子，所以稱為「混和」模型。
隨機因子 (random factor)
先來想一個問題：在看到可愛的狗勾時，多數人會提高說話的音頻，但每個人的音頻都一樣嗎？
Nope.
每個人天生的音頻不同，就算在某些情況下大多數人會提高/降低音頻，但還是存在個體差異 (individual differences)。
套用回來這裡的例子，用不同態度 (attitude) 時，每個人的音頻 (frequency) 可能有大致相同的變化 (變高/變低)，但仍存在差異。
圖一、每位受試者的平均音頻 (出處：原文)
圖一中 x 軸呈現的是每一位受試者，F 代表女生 (female)，M (male) 代表男生。y 軸呈現每位受試者的平均音頻 (也就是不論是何種態度)。
可以發現女生的平均音頻都比男生高，而在女生和男生各自的類別中，也有音頻較高和較低者。
所以 subject (受試者) 存在著會影響到模型的不可控因素 (個體差異)，於是我們把 subject 放入隨機因子，讓電腦知道我們要把這個變數納入考量。
frequency ~ attitude + gender + (1 | subject) + ε
(1 | subject) 的意思是，請電腦預設每位 subject 的起始值不同，有人天生音頻高，有人天生低。這個起始值通常稱為截距 (intercept)，在模型中以 1 表示。
截距的較詳細敘述：線性回歸幼幼班
接著，假設實驗人員準備了七種 scenario (情境)，例如：
- 向他人尋求協助
- 解釋為何遲到
諸如此類。
圖二、每個題目的平均音頻 (出處：原文)
圖二的 x 軸代表 7 個不同的題目，y 軸代表每個題目的平均音頻，這個平均來自 16 位受試者的答覆。
可以發現每個題目的平均音頻差距不大，但還是有些微差異，所以要做跟前面一樣的事情，在模型上加上「scenario」這個變數，告訴電腦我們也要把這個無法控制的因素納入考量。
frequency ~ attitude + gender + (1 | subject) + (1 | scenario) + ε
進入 R 環境
介紹完概念，接下來要實作了！
想使用這麼方便簡易的線性混和效應模型 (linear mixed-effects model)，首先要安裝並載入它的套件
install.packages(“lme4”)
library(lme4)
為使實作更方便，作者提供了他們的檔案，有兩種下載方式：
- 直接在 R 下載
politeness = read.csv("http://www.bodowinter.com/tutorial/politeness_data.csv")
- 複製這個網址，貼到搜尋引擎就會自動下載了：http://www.bodowinter.com/tutorial/politeness_data.csv
然後在 R 輸入
politeness = read.csv(file.choose( ))
這樣就會得到一份叫做 politeness 的資料。
可以檢查一下有沒有丟失的資料
which(is.na(politeness$frequency))
或
which(!complete.cases(politeness))
哈！然後就會發現在第 39 列有丟失的資料，不過在線性混和效應模型中，稍微有一些丟失資料不會產生太大影響。
介紹一下檔案，共有五欄：
- subject — 受試者代號。
- gender — 受試者性別，F 代表女性，M 代表男性。
- scenario — 如上面說的 7 種不同情境。
- attitude — pol 代表禮貌 (正式場合)，inf 代表非正式場合。
- frequency — 音頻，數字越大，音頻越高。
可以視覺化這些資料，方便觀察分布情況
boxplot(frequency ~ attitude*gender, col = c("white","lightgray"), politeness)
圖三、音頻盒狀圖 (出處：原文)
無論是 pol 還是 inf，女性音頻都比男性高；女性這組中 (左邊兩個盒狀圖)，inf 音頻較高，男性這組 (右邊兩個盒狀圖) 也有類似模式，但差距比女性組小。
建構模型
進入正題啦！
politeness.model <- lmer(frequency ~ attitude, data = politeness)
用 lmer( ) 函數建構模型，frequency 為應變數，attitude 為自變數，並告訴 R 要使用的資料是剛才載入的 politeness。
作者很幽默，他說在執行後會發現錯誤
Error in mkReTrms(findbars(RHSForm(formula)), fr) :No random effects terms specified in formula要提醒大家，這個模型需要隨機因子！
應該將我們剛才討論到的 subject 跟 item 列入隨機因子，然後用 summary( ) 來看結果 (見圖四)
politeness.model <- lmer(frequency ~ attitude + (1|subject) + (1|scenario), data = politeness)
summary(politeness.model)
圖四、politeness.model 的全結果 (出處：原文)
先看隨機因子：
圖五、Random effects (出處：原文)
這裡可以看到 scenario 跟 subject 的變異程度，後者變異程度比前者大很多，符合圖一跟圖二的討論。
接著發現多出一個奇怪的東西：Residual (殘差)。
可以影響音頻的因素，不可能只有人 (subject) 和情境 (scenario)，同意嗎？
所以 residual 才存在，他代表的是除了我們選定的兩個隨機因子外，會影響音頻的其他不知名因素，也就是前面一直看到的 ε。
接著看固定因子：
圖六、fixed effects (出處：原文)
attitudepol 是 attitude + pol。
你可能會想，那不是應該還要有一個 attitudeinf (attitude + inf) 嗎？其實 attitudepol 的 Estimate (-19.695) 隱藏了一個訊息：從 attitudeinf 到 attitudepol 要減掉 -19.695 Hz。
fixed effects 的結果呈現方式是這樣的，依據 A~Z，看哪一個 level (水準，比如這裡的 pol 跟 inf) 的開頭字母比較前面，就以他為基準，呈現另一個水準與它的差距。
以這裡為例， i 比 p 前面，所以以 inf 為基準，呈現 pol 與 inf 的差距。
那 Intercept 的奇怪數字 202.588 是怎麼來的？
它其實是 attitudeinf 情況下的男女平均。我們沒有告訴模型這裡有兩個性別，所以就被合在一起看了，但這樣其實不太合理。
修改一下模型，把 gender (性別) 也納入自變數 (我們清楚知道它如何影響音頻)
politeness.model <- lmer(frequency ~ attitude + gender + (1|subject) + (1|scenario), data = politeness)
summary(politeness.model)
隨機因子的結果變這樣：
圖七、隨機因子結果_增加性別為自變數 (出處：原文)
subject 的變異程度變小了！
之前還沒指出 gender 對 frequency 也有影響時，女性跟男性的資料被拉在一起看，所以產生很大的變異程度。
現在把 gender 列入自變數，區分女性和男性 (gender 的兩個水準)，gender 就不再是未知的變數，就不會跑到隨機因子的 subject 增加變異程度。
固定因子結果：
圖八、固定因子結果_增加性別為自變數 (出處：原文)
現在你應該會了。
genderM 的 Estimate (-108.517) 代表以 genderF 為基準，減掉 108.517 就會是 genderM 的平均。男生音頻較低，沒毛病。
而現在的 Intercept 比起圖六變高了，這是因為現在的 Intercept 是 genderF 在 attitudeinf 情況下的音頻平均。
顯著效果
做統計時常常需要參考的 p 值，用線性混和模型效應來做很直觀也很容易。
這裡用來取得 p 值的方法是 Likelihood ratio test (LRT，似然比檢驗)。
- 概念：比較兩個模型的概似性，看哪一個模型比較好。
- 做法：建構兩個模型，一個包含想看的自變數，一個不包含，比較兩模型看是否有顯著效果，如果有代表那個自變數很重要，會影響應變數。
有點抽象吧，作者提供了一個例子：
假設今天要去爬山，你揹了一桶水 (gallon of water，自變數) 跟一個手電筒 (flashlight，自變數)，想知道哪個東西會影響你的速度 (hiking speed，應變數)。
你的做法是先把 flashlight 拿掉然後爬一段路，結果發現 hiking speed 沒什麼變化；接著你把 flashlight 撿回來，把 gallon of water 丟掉，發現 hiking speed 明顯被影響了 (這裡先不談變快還是變慢，著重在「有影響」)！
結論：帶著 gallon of water 爬山會影響 hiking speed；帶著 flashlight 則不會。
把這個故事建構成模型，檢驗 gallon of water 會不會影響 hiking speed：
model1 = hiking speed ~ gallon of water + flashlight
model2 = hiking speed ~ flashlight
假設 model1 跟 model2 的比較有顯著效果，代表 gallon of water 對 hiking speed 有影響。
檢驗 flashligh 會不會影響 hiking speed：
model1 = hiking speed ~ gallon of water + flashlight
model3 = hiking speed ~ gallon of water
假設 model1 跟 model3 的比較沒有顯著效果，代表 flashlight 對 hiking speed 沒有影響。
model1 就是完全模型 (full model)，包含想看的自變數；model2 和 model3 則是虛無模型 (null model)，不包含想看的自變數。可以用式子相減的概念來看模型比較，就會知道現在是要看哪個自變數對應變數的影響。
現在可以回到 frequency 的例子建構模型了！
politeness.model <- lmer(frequency ~ attitude + gender + (1|subject) + (1|scenario), data=politeness, REML=FALSE)
politeness.null <- lmer(frequency ~ gender + (1|subject) + (1|scenario), data=politeness, REML=FALSE)
把兩個模型相減一下，可以知道現在想看的是 attitude 會不會影響 frequency。
另外出現了一個沒看過的東西：REML=FALSE。這裡不多加解釋，只要在做 LRT 時記得放入就可以。
接著，比較兩個模型的方式是用 anova( )
anova(politeness.null,politeness.model)
太棒了，終於得到結果了！
圖九、模型比較結果 (出處：原文)
看到 p < .001 (0.0006532)，有顯著效果，可以知道 attitude 對 frequency 有影響。
結果如何撰寫？
... attitude affected frequency (χ2 (1) = 11.62, p < .001) ...
1 是結果表中的 Df，11.62 是 Chisq。
交互作用
除了上面提到單一變數的效果，有時也會發生雖然某一自變數影響了應變數，但會根據另一自變數有不同的影響情況。
舉例來說，我們知道 attitude (pol / inf) 會影響 frequecy，但這個影響可能會因為 gender (F / M) 而不同。
比如女性在 pol 跟 inf 情況下的音頻變化與男性相反：
- pol, F → 高 frequency；inf, F → 低 frequency
- pol, M → 低 freuqency；inf, M → 高 frequency
或 attitude 的影響力在男女中不同：
- pol, F → 高 frequency；inf, F → 低 frequency
- pol, M → 中 freuqency；inf, M → 中 frequency
這就是 attitude 跟 gender 有交互作用的情況。
模型中標示交互作用的符號是「*」。
可以比較這兩個模型
politeness.interact <- lmer(frequency ~ attitude*gender + (1|subject) + (1|scenario), data=politeness, REML=FALSE)
politeness.reduced <- lmer(frequency ~ attitude + gender + (1|subject) + (1|scenario), data=politeness, REML=FALSE)
如果結果有顯著效果，則 attitude 跟 gender 有交互作用。
隨機斜率 隨機截距
在「隨機因子」段落提到，每位受試者天生的音頻就不同、每個題目蒐集到的音頻基準也不同，所以我們在模型中指定了 subject 跟 item 為隨機因子。
這樣的指定下，模型會針對每位受試者和題目產生不一樣的 intercept (截距)，所以是一個隨機截距模型 (random intercept model)。
但是，固定因子 (attitude) 對每位受試者產生的效果可能也不同；同理，對每一題的影響也不見得相同。
圖十、受試者在不同題目的音頻
我試著用線性回歸畫了一張綜合圖，每個灰色框格是一位受試者的資料，粉色線是 inf 情況、藍色線是 pol 情況。x 軸是 scenario、y 軸是 frequency。
很明顯可以看到每位受試者在相同題目的表現不同，因此線產生不同 slope (斜率)。
圖十一、題目在不同受試者情況下的音頻
將 x 軸換成 subject、y 軸一樣是 frequency。
同樣可以發現每個題目在相同受試者的情況下，也會得到不同音頻 (雖然趨勢看起來很相近，但要注意 y 軸數值的變化)，線產生不同 slope。
所以，每位受試者和每個題目不只可以有不同的 intercept，還可以有不同的 slope。每個人或每個題目都有自己專屬的線來代表他的數據。而這樣的差異是經過 attitude (pol / inf) 的測試得到的，所以現在可以修改一下之前的模型
politeness.model = lmer(frequency ~ attitude + gender + (1+attitude|subject) + (1+attitude|scenario), data=politeness,REML=FALSE)
隨機因子中，在 intercept (1) 的後面加上 slope (attitude)。
取 (1+attitude|subject) 來說明，這裡代表我們允許每位受試者有不同的 intercept 和不同的 slope。
以我的話來說，就是每個人都可以有不同的起始值，也可以對相同情境有不同反應。
筑言筑語
真正在自己的實驗中就會發現事情沒有這麼美好，想要指定全部的隨機截距跟隨機斜率時，模型常常會收斂失敗。
這裡作者給的建議是：在模型可以建構成功的情況下，盡量保持模型最複雜化。
例如先把想指定的都指定上去，看看會不會成功。如果不行的話，試著先拿掉一個東西，像是一個斜率，然後再嘗試看看，試到成功為止。
Yaaay，終於寫完了。
這篇文章從一月寫到現在，中間因為遇到不懂的概念，所以又跑去學別的東西。寫的過程中一直在釐清之前未曾發現過的問題，還有思考怎麼呈現才會比較簡單易懂。
總之，自己仔仔細細的順過一遍，感覺很好！
如果有哪裡觀念錯了，也歡迎大家告訴我這個統計小白。
引用
最後，開發 R 跟 lme4 package 的人或團隊都很辛苦，所以別忘記，如果在資料分析過程中有使用到這些資源，要提到他們一下喔。
可以用下面的方式找到該引用誰：
R
citation()
我的是這個
R Core Team (2023). _R: A Language and Environment for Statistical Computing_. R
Foundation for Statistical Computing, Vienna, Austria. <https://www.R-project.org/>.
lme4 package
citation("lme4")
我的是這個
Douglas Bates, Martin Maechler, Ben Bolker, Steve Walker (2015). Fitting Linear
Mixed-Effects Models Using lme4. Journal of Statistical Software, 67(1), 1-48.
doi:10.18637/jss.v067.i01.
