以下只是我自己學習LMM所做的筆記褒纲，并不保證完全正確。
基于R的線性混合效應模型(linear mixed effects models)分析
why 為什么要采用線性混合效應模型分析：
看了一些資料，我自己的理解是：
在我們的實驗中刑顺，除了我們通過實驗法操縱的自變量外祖凫，可能還存在其他我們沒辦法操縱和控制的影響因變量的因素寝杖。
例如型豁，在探究性別和說話方式對音高的影響實驗中，我們采用2因素的（說話方式：有禮貌的/非正式的 * 性別：男/女）混合實驗設(shè)計尚蝌，通過操縱說話方式和性別迎变，來測量不同性別（gender）的被試在不同場景（scenario）下，采用不同的說話方式（attitude）時飘言，的音高（pitch）衣形。
如果采用往常的數(shù)據(jù)分析，我們可能會在數(shù)據(jù)分析的時候，將不同場景的數(shù)據(jù)求平均谆吴，然后做2（gender）* 2（attitude）的重復測量方差分析倒源。
但是這么做的一個缺點就是，只考慮了gender和attitude帶來的固定變異和隨機誤差句狼，但是并沒有考慮到不同的被試之間存在個體差異（方差分析值考慮不同實驗處理內(nèi)部的個體差異笋熬，但是沒有考慮所有被試間的個體 差異），音高基線不同腻菇，而且不同場景（scenario）下的音高基線可能原本也不同胳螟。此外，說話方式對因高的影響可能還會因被試和場景的不同而影響程度不同筹吐。同一被試的多個數(shù)據(jù)之間可能存在關(guān)聯(lián)（數(shù)據(jù)不獨立）
上述這些因素糖耸，是我們在實驗中無法控制的，不可預料的丘薛，非系統(tǒng)性的因素嘉竟，又被稱為隨機效應。而自變量（性別洋侨，說話方式）是我們在實驗中可以操縱的舍扰，預期的，系統(tǒng)的效應凰兑，又被稱為固定效應妥粟。
線性混合效應模型，既可以考慮到固定效應吏够，又可以考慮到隨機效應勾给。所以是混合效應。
做線性混合效應模型分析需要滿足的條件：
-
Linearity：自變量和因變量是線性關(guān)系锅知。
既是線性模型播急，必然要滿足線性假設(shè)。
檢驗方法：可以通過畫散點圖或者是殘差圖來檢驗
# 畫圖 x=age y=pitch, xlim指定x軸的范圍為1-80售睹，ylim指定y軸的范圍180-280 plot(age,pitch,xlim=c(0, 80), ylim=c(180, 280))
# abline為圖形添加一條線桩警，添加的方式是，用線性擬合昌妹。x是age捶枢，y是pitch abline(lm(pitch~age),Ity=2)
plot(fitted(xmdl),residuals(xmdl))#畫殘差圖residual plot
標準的殘差圖應該是類似塊狀的。 -
Absence of collinearity：非共線性飞崖。
這里共線性是指：固定效應之間存在關(guān)聯(lián)/相關(guān)烂叔。
比方說，你感興趣的是平均說話速度如何影響智力評級(即固歪，說話速度越快的人被認為越聰明)…
你測量了幾個不同的談話速度指標蒜鸡，例如胯努，每秒的音節(jié)數(shù)、每秒的字數(shù)和每秒的句子數(shù)逢防。這些不同的衡量標準之間是相互關(guān)聯(lián)的（也叫共線性）叶沛，因為如果你說得更快，你就會在給定的時間內(nèi)說出更多的音節(jié)忘朝、單詞和句子灰署。
而線性混合效應模型分析要求不能有共線性。 -
Homoskedasticity：方差齊性辜伟。
方差齊性是一個重要的假設(shè)氓侧。是指：在預測值的范圍內(nèi)，數(shù)據(jù)的方差應該大致相等导狡。
方差（樣本方差）是每個樣本值與全體樣本值的平均數(shù)之差的平方的和的平均數(shù)
方差計算公式
為了滿足方差齊性假設(shè)约巷，模型的殘差需要大致在預測值范圍內(nèi)，有類似/相似的偏移量旱捧。畫出的圖形類似塊狀
檢驗方法：同線性独郎，依然是殘差圖 -
Normality of residuals：殘差正態(tài)分布。
殘差的正態(tài)性假設(shè)是最不重要的假設(shè)枚赡。
研究人員對檢驗這一假設(shè)的重視程度各不相同氓癌。
檢驗方法：直方圖或者Q-Q圖：
hist(residuals(xmdl))#直方圖
qqnorm(residuals(xmdl))#Q-Q圖
Q-Q圖，數(shù)據(jù)都落在一條直線上贫橙。說明類似正態(tài)分布贪婉。
如果您擁有極端的數(shù)據(jù)點，該如何繼續(xù)卢肃？
簡單地排除這些點而只報告縮減后的數(shù)據(jù)集的結(jié)果肯定是不合理的疲迂。
更好的方法是先使用有極端值的數(shù)據(jù)運行分析，然后再運行一次沒有極端值的數(shù)據(jù)…莫湘。然后尤蒿，您可以報告這兩種分析，并說明對結(jié)果的解釋是否改變幅垮。
唯一沒問題的情況是：當極端值的存在是明顯錯誤時腰池，例如，沒有意義的值(例如忙芒，負年齡示弓、負高度)或明顯是由于數(shù)據(jù)采集階段的技術(shù)錯誤而導致的值(例如，語音音調(diào)值為0)呵萨。 -
Absence of influential data points：沒有極端的數(shù)據(jù)點奏属。
檢驗方法：
預先定義一個向量，該向量的元素數(shù)量與數(shù)據(jù)集中的行數(shù)相同甘桑。
然后拍皮，循環(huán)瀏覽每一行。對于每次循環(huán)跑杭，創(chuàng)建一個沒有該行的新混合模型(這是通過POP[-i铆帽，]實現(xiàn)的)。
然后德谅，通過函數(shù)fixef()提取您感興趣的任何系數(shù)爹橱，也就是指定some number的值。
您需要調(diào)整此代碼以適應您的分析窄做。
除了數(shù)據(jù)框和變量的名稱之外愧驱，還需要對模型運行一次fixef()，以便知道相關(guān)系數(shù)的位置椭盏。
在我們的例子中组砚，我會在其中放一個“2”，因為“atitdepol”的效果在系數(shù)列表中排在第二位掏颊≡愫欤“1”會給出截距，總是系數(shù)表中提到的第一個系數(shù)乌叶。
輸出系數(shù)之后盆偿，觀察輸出的系數(shù)值，尋找與斜率絕對值相差至少一半的值准浴。比方說事扭，我的斜率是2…。那么系數(shù)的值為1或-1對我來說就是一個警示乐横。如果是-4的斜率求橄，系數(shù)值為2或-2對我來說是個警示。all.res=numeric(nrow(mydataframe)) for(i in 1:nrow(mydataframe)){ myfullmodel=lmer(response~predictor+ (1+predictor|randomeffect),POP[-I,]) all.res[i]=fixef(myfullmodel)[some number] } -
Independence !!!!!!!：數(shù)據(jù)之間的獨立性晰奖。
獨立性是最重要的假設(shè)谈撒。它強調(diào)數(shù)據(jù)之間應該保持獨立。
獨立到底是什么呢匾南？
理想的情況是擲硬幣或擲骰子：每一次擲硬幣和每一次擲骰子都完全獨立于前一次擲硬幣或擲骰子的結(jié)果啃匿。
當你運行線性模型分析時，這同樣適用于你的數(shù)據(jù)點蛆楞。
所以在采用線性混合效應模型進行分析的時候溯乒，就應該把可能影響獨立性的效應，當成隨機效應豹爹，告訴模型裆悄。以此來解決非獨立的問題。
如何在R中做線性混合效應模型分析
因變量是類別變量時臂聋，需要用廣義線性混合效應模型光稼，用的是glmer()函數(shù)或南。
install.packages(“l(fā)me4”) # 安裝lme4包
library(lme4) #載入包
politeness = read.csv(file.choose( )) #手動選擇讀取數(shù)據(jù)
# 檢查數(shù)據(jù)politeness的frequency一列是否有NA值。
which(is.na(politeness$frequency))
# 檢查數(shù)據(jù)politeness是否有NA值艾君。
which(!complete.cases(politeness))
# 線性混合效應模型LME分析
## 固定效應：attitude 和 gender
## 隨機效應：1. subject 和 scenario的截距（intercepts for subjects and items）
## 這里截距用1表示采够。
## 其實，應該就是冰垄，不同subject 和 item基線的不同 作為隨機效應蹬癌。
## 隨機效應：2. by-subject and by-item random slopes for the effect of politeness
## 隨機效應2的意思是，針對politeness效應的按subject和按item的隨機斜率虹茶。
## 其實逝薪，應該就是，不同的subject和item對politeness影響的隨機效應蝴罪。
politeness.fullmodel = lmer(frequency ~ attitude + gender +
(1 + attitude|subject) + (1 + attitude|scenario),
data=politeness,
REML=FALSE)
# REML=FALSE董济，是為了后面使用似然比檢驗比較模型。
# 輸出模型的結(jié)果
summary(politeness.fullmodel)
# 建立沒有attitude的模型要门。
politeness.nullmodel = lmer(frequency ~ gender +
(1+attitude|subject) + (1+attitude|scenario),
data=politeness,
REML=FALSE)
summary(politeness.nullmodel)
# 比較有attitude的fullmodel和沒有attitude的nullmodel
# 如果結(jié)果差異顯著感局，則代表attitude對結(jié)果有顯著影響。
# 這里的比較方法是：似然比檢驗 Likelihood Ratio Test
anova(politeness.nullmodel, politeness.fullmodel)
# 如果要檢驗交互作用暂衡。則需要以下代碼询微，比較以下兩個模型：
# 就是把+變成*，*表示交互狂巢。
politeness.fullmodel = lmer(frequency ~ attitude * gender +
(1 + attitude|subject) + (1 + attitude|scenario),
data=politeness,
REML=FALSE)
politeness.nullmodel = lmer(frequency ~ attitude + gender +
(1 + attitude|subject) + (1 + attitude|scenario),
data=politeness,
REML=FALSE)
anova(politeness.nullmodel, politeness.fullmodel)
輸出的結(jié)果該怎么理解：
politeness.fullmodel = lmer(frequency ~ attitude + gender +
(1 + attitude|subject) + (1 + attitude|scenario),
data=politeness,
REML=FALSE)
結(jié)果顯示撑毛，p=0.009597，差異顯著唧领。也就是說藻雌，attitude的確會顯著的影響pitch。
數(shù)據(jù)分析部分的英文寫法為：
“We used R (R Core Team, 2012) and lme4 (Bates, Maechler & Bolker, 2012) to perform a linear mixed effects analysis of the relationship between pitch and politeness. As fixed effects, we entered politeness and gender (without interaction term) into the model. As random effects, we had intercepts for subjects and items, as well as by-subject and by-item random slopes for the effect of politeness. Visual inspection of residual plots did not reveal any obvious deviations from homoscedasticity or normality. P-values were obtained by likelihood ratio tests of the full model with the effect in question against the model without the effect in question.”
結(jié)果報告的英文寫法為：
politeness affected pitch (χ2 (1)=6.7082, p=0.009597), lowering it by about 19.747 Hz ± 5.92 (standard errors)
方法：
首先構(gòu)建全模型
如果報錯斩个，通過那個似然估計來判斷削減效應之后的模型與只有截距的零模型 對比 或者 通過不同模型的AIC對比
最終確定一個模型
然后建模 得到固定效應和主效應胯杭，交互作用等。（固定效應其實是斜率受啥、截距的效應做个。就是回歸方程。代表的是在方程里的貢獻度滚局。而主效應居暖，類似于方差分析的結(jié)果去理解。代表藤肢，變量的不同水平存在差異太闺。）
這里需要明確一個問題。winter的教程里提到的似然估計法嘁圈，他把他當作是計算p值的方法省骂。我推測可能是當時蟀淮，lme4這個包還不能輸出固定效應的p值。所以才這么做钞澳。
用lm()做線性回歸的時候 輸出結(jié)果的解釋：
Multiple R-squared: 0.2704灭贷，它表示模型對總體方差的解釋能力。具體意思可以解釋為略贮，總體方差中的27.04%，可以由這個模型來解釋仗岖。Adjusted R-squared是對Multiple R-squared的矯正逃延，主要是考慮了固定效應。固定效應越多轧拄，該值越低揽祥。
*號代表 主效應和交互作用
:號代表 僅有交互作用
+號代表 只有主效應
參考文獻：
Linear models and linear mixed effects models in R with linguistic applications
Fitting Linear Mixed-E?ects Models Using lme4
一些筆記：
用lmerTest包與lmer()函數(shù)建模
全模型與零模型
不同隨機效應之間random effect的相關(guān)系數(shù)都到1了。說明相關(guān)非常強檩电。會有多重共線性拄丰。有的變量可能是多余的，但我考慮進去了俐末。
這里有一個容忍度料按，tol。容忍度越高卓箫，閾限是越緊的载矿。設(shè)置為10的-4次方可能更好掷匠。
固定效應 主效應和交互作用
固定效應輸出的估計值婆芦，t值和p值。該怎么理解能真。
比如 directionright 這一行就代表旅急，left和right兩個方向的比較逢勾。做t檢驗。directionright代表right和left平均值的差異藐吮，t值和p值溺拱，是拿這個差異和0比較。所以也就相當于是谣辞，left和right之間的比較盟迟。這個估計值實際上就是回歸方程的斜率。因為這里是left和right2個水平潦闲，將left當作1攒菠，right編碼為2。所以當x從left（1）到right（2）時歉闰，y降低了估計值的大小辖众，所以就是估計值代表斜率卓起。
固定效應和主效應是不同的。如果固定效應有兩個水平時凹炸，固定效應的p值與主效應的p值是相同的戏阅。 也就是fixed effexts輸出的p值與anova(model)輸出的p值相同。因為兩個水平的話啤它，t檢驗和F檢驗是相同的奕筐。
modelNew是model的名稱
pairwise～A是配對檢驗。A是需要配對檢驗的變量变骡。
adjust是多重比較矯正离赫。如果只有2個水平，就不需要這個了塌碌。矯正的方式圖上都有渊胸，可以自己去嘗試。
可能會報錯台妆，說超過3600翎猛，這里時候需要調(diào)整之后在運行，要注意一下接剩。設(shè)置為3600之后切厘，會運行的很慢。
看一下結(jié)果懊缺，contrast部分迂卢，左右方向的差異是-10，是顯著的桐汤，p<0.001
簡單主效應分析
看結(jié)果而克，左方向上，距離的主效應顯著怔毛，p<0.001员萍。
右方向上，也是一樣拣度。
簡單主效應的事后比較
保證方向相同碎绎，比較距離。
因子的對比方式與固定效應系數(shù)
manual代表手動擋抗果，估計值表示筋帖，手動擋相對于自動擋，增加了7.24冤馏。 這里相當于17.15是截距日麸，這里代表的是auto。相當于把auto當做基線，下面manual是跟auto對比的差異7.24代箭。
通過contrast()函數(shù)可以查看對比方式墩划。
下圖是單因素4水平舉例。模型會把4水平編碼為3水平嗡综。
ABCD編碼為BCD三個水平乙帮。這種編碼方式叫treatment coding，是0 1編碼极景。
還有一種編碼方式是sum coding察净，是1 0 -1的方式。
手動生成虛擬變量
用主成分分析優(yōu)化模型
全模型會存在畸形的協(xié)方差矩陣盼樟。所以需要削減隨機斜率氢卡。
優(yōu)化模型：去掉無效的/多余的 隨機效應，確定有效的隨機效應恤批。
把方差和標準差比較小的成分給去掉。說明這個隨機效應解釋的變異是比較小的裹赴，可能是冗余的成分喜庞。
一個豎線改成2個豎線，就變成零相關(guān)模型棋返。
但是結(jié)果仍然存在畸形協(xié)方差延都。
需要采用rePCA()函數(shù)做主成分分析。
主成分分析結(jié)果睛竣，我們關(guān)注的是proportion of variance這一行晰房。在本例子中，id部分射沟，有4個主成分殊者。但是第4個，值是0验夯，代表不能解釋任何方差猖吴。
item部分，有三個成分值都是0挥转。到底要去掉哪個海蔽。
這個時候需要去查看零相關(guān)模型的方差標注差。
我們可以看到id中绑谣，id.3 interaction 標準差是0党窜，就代表交互作用是多余的。
item上借宵，主成分分析是顯著3個成分多余幌衣，標準差最小的3個分別是 item item.1 item.3這三個，分別代表結(jié)局壤玫，左右的差異泼掠，和交互作用怔软。
刪除之后，建立一個新的模型择镇。去掉隨機截距的方法是把1改成-1挡逼。
發(fā)現(xiàn)新模型沒有警告信息了。
比較新模型和全模型之間的差別腻豌。通過anova()函數(shù)來比較家坎。
定義事先對比/編碼方式
定義對比方式后的結(jié)果怎么看
> Model = lm(data = ToothGrowth, len~supp*dose, contrasts = list(supp=contr.simple(n = 2), dose = contr.simple(3)))
> summary(Model)
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 18.8133 0.4688 40.130 < 2e-16 ***
supp1 -3.7000 0.9376 3.946 0.000231 ***
dose1 9.1300 1.1484 7.951 1.19e-10 ***
dose2 15.4950 1.1484 13.493 < 2e-16 ***
supp1:dose1 -0.6800 2.2967 0.296 0.768308
supp1:dose2 -5.3300 2.2967 -2.321 0.024108 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 3.631 on 54 degrees of freedom
Multiple R-squared: 0.7937, Adjusted R-squared: 0.7746
F-statistic: 41.56 on 5 and 54 DF, p-value: < 2.2e-16
supp兩種處理，VC（用1表示）和OJ（用2表示）
dose有3種吝梅，0.5（用1表示）虱疏，1（用2表示），2（用3表示）
我個人的理解是這樣苏携，實際上做瞪，如果自變量是因子型，回歸輸出的結(jié)果可以這么理解：這里用的是simple編碼右冻，是VC和OJ中装蓬，首字母在前的那個當作基線，也就是OJ作為基線纱扭。
supp1實際就是它跟基線做比較牍帚，也就是VC-OJ的差異，Estimate的值可以說是回歸方程的系數(shù)乳蛾，也可以說是暗赶，VC跟OJ平均值的差異，t檢驗得到t值和p值肃叶，p值顯著就說明蹂随，差值跟0比，有顯著差異因惭，也就代表糙及，VC和OJ差異顯著。
dose1是1.0和0.5兩個水平的差異筛欢。
dose2是2和0.5兩個水平的差異浸锨。
supp1：dose1 就等于（VC-OJ）*(1.0-0.5)
其實我個人感覺就是類似于兩兩比較的結(jié)果。只是每一個都是跟基線比較的版姑。比較的方法是t檢驗柱搜，回頭可以驗證一下。
驗證過了剥险，確實是這樣聪蘸。當只有一個因素，2個水平時，確實跟獨立樣本t檢驗的結(jié)果相同健爬。如果還有其他因素控乾，就不是了，因為標準誤會有變化娜遵。
simple coding下蜕衡，回歸系數(shù)等于真實效應，基線通常是選取首字母在前的因子水平设拟。但是這只適用于慨仿，不同實驗處理下被試量相等的情況，如果不等的話纳胧，不適用镰吆。
Scaled residuals:
Min 1Q Median 3Q Max
-2.6261 -0.1180 -0.1180 0.0583 15.0471
Random effects:
Groups Name Variance Std.Dev.
sub (Intercept) 0.06286 0.2507
Residual 0.32016 0.5658
Number of obs: 774, groups: sub, 86
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 0.06460 0.03383 84.00000 1.909 0.05962 .
power1 0.02584 0.06766 84.00000 0.382 0.70351
member1 0.13566 0.04982 684.00000 2.723 0.00663 **
member2 0.05814 0.04982 684.00000 1.167 0.24361
power1:member1 -0.03876 0.09964 684.00000 -0.389 0.69739
power1:member2 0.11628 0.09964 684.00000 1.167 0.24361
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr) power1 membr1 membr2 pwr1:1
power1 0.000
member1 0.000 0.000
member2 0.000 0.000 0.500
powr1:mmbr1 0.000 0.000 0.000 0.000
powr1:mmbr2 0.000 0.000 0.000 0.000 0.500
> a<-aggregate(power_exp1_data$puni_amount,by=list(power_exp1_data$power),FUN=mean)
> b<-aggregate(power_exp1_data$puni_amount,by=list(power_exp1_data$member),FUN=mean)
> a
Group.1 x
1 high 0.05167959
2 low 0.07751938
> b
Group.1 x
1 ingroup 0.00000000
2 outgroup 0.13565891
3 unclassified 0.05813953
> a[2,2]-a[1,2]
[1] 0.02583979
> b[2,2]-b[1,2]
[1] 0.1356589
> b[3,2]-b[1,2]
[1] 0.05813953
> c<-aggregate(power_exp1_data$puni_amount,by=list(power_exp1_data$power,power_exp1_data$member),FUN=mean)
> c
Group.1 Group.2 x
1 high ingroup 0.0000000
2 low ingroup 0.0000000
3 high outgroup 0.1550388
4 low outgroup 0.1162791
5 high unclassified 0.0000000
6 low unclassified 0.1162791
> (c[4,3]-c[2,3])-(c[3,3]-c[1,3])
[1] -0.03875969
> (c[6,3]-c[2,3])-(c[5,3]-c[1,3])
[1] 0.1162791
> (a[1,2]+a[2,2])/2
[1] 0.06459948
