以下只是我自己学习LMM所做的笔记，并不保证完全正确。
基于R的线性混合效应模型(linear mixed effects models)分析
why 为什么要采用线性混合效应模型分析：
看了一些资料，我自己的理解是：
在我们的实验中，除了我们通过实验法操纵的自变量外，可能还存在其他我们没办法操纵和控制的影响因变量的因素。
例如，在探究性别和说话方式对音高的影响实验中，我们采用2因素的（说话方式：有礼貌的/非正式的 * 性别：男/女）混合实验设计，通过操纵说话方式和性别，来测量不同性别（gender）的被试在不同场景（scenario）下，采用不同的说话方式（attitude）时，的音高（pitch）。
如果采用往常的数据分析，我们可能会在数据分析的时候，将不同场景的数据求平均，然后做2（gender）* 2（attitude）的重复测量方差分析。
但是这么做的一个缺点就是，只考虑了gender和attitude带来的固定变异和随机误差，但是并没有考虑到不同的被试之间存在个体差异（方差分析值考虑不同实验处理内部的个体差异，但是没有考虑所有被试间的个体 差异），音高基线不同，而且不同场景（scenario）下的音高基线可能原本也不同。此外，说话方式对因高的影响可能还会因被试和场景的不同而影响程度不同。同一被试的多个数据之间可能存在关联（数据不独立）
上述这些因素，是我们在实验中无法控制的，不可预料的，非系统性的因素，又被称为随机效应。而自变量（性别，说话方式）是我们在实验中可以操纵的，预期的，系统的效应，又被称为固定效应。
线性混合效应模型，既可以考虑到固定效应，又可以考虑到随机效应。所以是混合效应。
做线性混合效应模型分析需要满足的条件：
-
Linearity：自变量和因变量是线性关系。
既是线性模型，必然要满足线性假设。
检验方法：可以通过画散点图或者是残差图来检验
# 画图 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：非共线性。
这里共线性是指：固定效应之间存在关联/相关。
比方说，你感兴趣的是平均说话速度如何影响智力评级(即，说话速度越快的人被认为越聪明)…
你测量了几个不同的谈话速度指标，例如，每秒的音节数、每秒的字数和每秒的句子数。这些不同的衡量标准之间是相互关联的（也叫共线性），因为如果你说得更快，你就会在给定的时间内说出更多的音节、单词和句子。
而线性混合效应模型分析要求不能有共线性。 -
Homoskedasticity：方差齐性。
方差齐性是一个重要的假设。是指：在预测值的范围内，数据的方差应该大致相等。
方差（样本方差）是每个样本值与全体样本值的平均数之差的平方的和的平均数
方差计算公式
为了满足方差齐性假设，模型的残差需要大致在预测值范围内，有类似/相似的偏移量。画出的图形类似块状
检验方法：同线性，依然是残差图 -
Normality of residuals：残差正态分布。
残差的正态性假设是最不重要的假设。
研究人员对检验这一假设的重视程度各不相同。
检验方法：直方图或者Q-Q图：
hist(residuals(xmdl))#直方图
qqnorm(residuals(xmdl))#Q-Q图
Q-Q图，数据都落在一条直线上。说明类似正态分布。
如果您拥有极端的数据点，该如何继续？
简单地排除这些点而只报告缩减后的数据集的结果肯定是不合理的。
更好的方法是先使用有极端值的数据运行分析，然后再运行一次没有极端值的数据…。然后，您可以报告这两种分析，并说明对结果的解释是否改变。
唯一没问题的情况是：当极端值的存在是明显错误时，例如，没有意义的值(例如，负年龄、负高度)或明显是由于数据采集阶段的技术错误而导致的值(例如，语音音调值为0)。 -
Absence of influential data points：没有极端的数据点。
检验方法：
预先定义一个向量，该向量的元素数量与数据集中的行数相同。
然后，循环浏览每一行。对于每次循环，创建一个没有该行的新混合模型(这是通过POP[-i，]实现的)。
然后，通过函数fixef()提取您感兴趣的任何系数，也就是指定some number的值。
您需要调整此代码以适应您的分析。
除了数据框和变量的名称之外，还需要对模型运行一次fixef()，以便知道相关系数的位置。
在我们的例子中，我会在其中放一个“2”，因为“atitdepol”的效果在系数列表中排在第二位。“1”会给出截距，总是系数表中提到的第一个系数。
输出系数之后，观察输出的系数值，寻找与斜率绝对值相差至少一半的值。比方说，我的斜率是2…。那么系数的值为1或-1对我来说就是一个警示。如果是-4的斜率，系数值为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 !!!!!!!：数据之间的独立性。
独立性是最重要的假设。它强调数据之间应该保持独立。
独立到底是什么呢？
理想的情况是掷硬币或掷骰子：每一次掷硬币和每一次掷骰子都完全独立于前一次掷硬币或掷骰子的结果。
当你运行线性模型分析时，这同样适用于你的数据点。
所以在采用线性混合效应模型进行分析的时候，就应该把可能影响独立性的效应，当成随机效应，告诉模型。以此来解决非独立的问题。
如何在R中做线性混合效应模型分析
因变量是类别变量时，需要用广义线性混合效应模型，用的是glmer()函数。
install.packages(“lme4”) # 安装lme4包
library(lme4) #载入包
politeness = read.csv(file.choose( )) #手动选择读取数据
# 检查数据politeness的frequency一列是否有NA值。
which(is.na(politeness$frequency))
# 检查数据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，是为了后面使用似然比检验比较模型。
# 输出模型的结果
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
# 如果结果差异显著，则代表attitude对结果有显著影响。
# 这里的比较方法是：似然比检验 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)
输出的结果该怎么理解：
politeness.fullmodel = lmer(frequency ~ attitude + gender +
(1 + attitude|subject) + (1 + attitude|scenario),
data=politeness,
REML=FALSE)
结果显示，p=0.009597，差异显著。也就是说，attitude的确会显著的影响pitch。
数据分析部分的英文写法为：
“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.”
结果报告的英文写法为：
politeness affected pitch (χ2 (1)=6.7082, p=0.009597), lowering it by about 19.747 Hz ± 5.92 (standard errors)
方法：
首先构建全模型
如果报错，通过那个似然估计来判断削减效应之后的模型与只有截距的零模型 对比 或者 通过不同模型的AIC对比
最终确定一个模型
然后建模 得到固定效应和主效应，交互作用等。（固定效应其实是斜率、截距的效应。就是回归方程。代表的是在方程里的贡献度。而主效应，类似于方差分析的结果去理解。代表，变量的不同水平存在差异。）
这里需要明确一个问题。winter的教程里提到的似然估计法，他把他当作是计算p值的方法。我推测可能是当时，lme4这个包还不能输出固定效应的p值。所以才这么做。
用lm()做线性回归的时候 输出结果的解释：
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-Effects Models Using lme4
一些笔记：
用lmerTest包与lmer()函数建模
全模型与零模型
不同随机效应之间random effect的相关系数都到1了。说明相关非常强。会有多重共线性。有的变量可能是多余的，但我考虑进去了。
这里有一个容忍度，tol。容忍度越高，阈限是越紧的。设置为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，这里时候需要调整之后在运行，要注意一下。设置为3600之后，会运行的很慢。
看一下结果，contrast部分，左右方向的差异是-10，是显著的，p<0.001
简单主效应分析
看结果，左方向上，距离的主效应显著，p<0.001。
右方向上，也是一样。
简单主效应的事后比较
保证方向相同，比较距离。
因子的对比方式与固定效应系数
manual代表手动挡，估计值表示，手动挡相对于自动挡，增加了7.24。 这里相当于17.15是截距，这里代表的是auto。相当于把auto当做基线，下面manual是跟auto对比的差异7.24。
通过contrast()函数可以查看对比方式。
下图是单因素4水平举例。模型会把4水平编码为3水平。
ABCD编码为BCD三个水平。这种编码方式叫treatment coding，是0 1编码。
还有一种编码方式是sum coding，是1 0 -1的方式。
手动生成虚拟变量
用主成分分析优化模型
全模型会存在畸形的协方差矩阵。所以需要削减随机斜率。
优化模型：去掉无效的/多余的 随机效应，确定有效的随机效应。
把方差和标准差比较小的成分给去掉。说明这个随机效应解释的变异是比较小的，可能是冗余的成分。
一个竖线改成2个竖线，就变成零相关模型。
但是结果仍然存在畸形协方差。
需要采用rePCA()函数做主成分分析。
主成分分析结果，我们关注的是proportion of variance这一行。在本例子中，id部分，有4个主成分。但是第4个，值是0，代表不能解释任何方差。
item部分，有三个成分值都是0。到底要去掉哪个。
这个时候需要去查看零相关模型的方差标注差。
我们可以看到id中，id.3 interaction 标准差是0，就代表交互作用是多余的。
item上，主成分分析是显著3个成分多余，标准差最小的3个分别是 item item.1 item.3这三个，分别代表结局，左右的差异，和交互作用。
删除之后，建立一个新的模型。去掉随机截距的方法是把1改成-1。
发现新模型没有警告信息了。
比较新模型和全模型之间的差别。通过anova()函数来比较。
定义事先对比/编码方式
定义对比方式后的结果怎么看
> 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表示）
我个人的理解是这样，实际上，如果自变量是因子型，回归输出的结果可以这么理解：这里用的是simple编码，是VC和OJ中，首字母在前的那个当作基线，也就是OJ作为基线。
supp1实际就是它跟基线做比较，也就是VC-OJ的差异，Estimate的值可以说是回归方程的系数，也可以说是，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)
其实我个人感觉就是类似于两两比较的结果。只是每一个都是跟基线比较的。比较的方法是t检验，回头可以验证一下。
验证过了，确实是这样。当只有一个因素，2个水平时，确实跟独立样本t检验的结果相同。如果还有其他因素，就不是了，因为标准误会有变化。
simple coding下，回归系数等于真实效应，基线通常是选取首字母在前的因子水平。但是这只适用于，不同实验处理下被试量相等的情况，如果不等的话，不适用。
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
