R数据分析R代码Data science

基于R的线性混合效应模型分析

2021-02-28  本文已影响0人  毛线东东a

以下只是我自己学习LMM所做的笔记,并不保证完全正确。

基于R的线性混合效应模型(linear mixed effects models)分析

why 为什么要采用线性混合效应模型分析:

看了一些资料,我自己的理解是:
在我们的实验中,除了我们通过实验法操纵的自变量外,可能还存在其他我们没办法操纵和控制的影响因变量的因素。
例如,在探究性别和说话方式对音高的影响实验中,我们采用2因素的(说话方式:有礼貌的/非正式的 * 性别:男/女)混合实验设计,通过操纵说话方式和性别,来测量不同性别(gender)的被试在不同场景(scenario)下,采用不同的说话方式(attitude)时,的音高(pitch)。
如果采用往常的数据分析,我们可能会在数据分析的时候,将不同场景的数据求平均,然后做2(gender)* 2(attitude)的重复测量方差分析。
但是这么做的一个缺点就是,只考虑了gender和attitude带来的固定变异和随机误差,但是并没有考虑到不同的被试之间存在个体差异(方差分析值考虑不同实验处理内部的个体差异,但是没有考虑所有被试间的个体 差异),音高基线不同,而且不同场景(scenario)下的音高基线可能原本也不同。此外,说话方式对因高的影响可能还会因被试和场景的不同而影响程度不同。同一被试的多个数据之间可能存在关联(数据不独立)

image.png

上述这些因素,是我们在实验中无法控制的,不可预料的,非系统性的因素,又被称为随机效应。而自变量(性别,说话方式)是我们在实验中可以操纵的,预期的,系统的效应,又被称为固定效应
线性混合效应模型,既可以考虑到固定效应,又可以考虑到随机效应。所以是混合效应。

image.png

做线性混合效应模型分析需要满足的条件:

如何在R中做线性混合效应模型分析

image.png

因变量是类别变量时,需要用广义线性混合效应模型,用的是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)
politeness.fullmodel的输出结果.png
fullmodel和nullmodel比较的结果.png

结果显示,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()函数建模

用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。
右方向上,也是一样。
简单主效应的事后比较
简单效应的事后检验
代码展示
保证方向相同,比较距离。
结果展示
简单效应 事后检验 事后多重比较

因子的对比方式与固定效应系数

固定效应中的回归系数与比较方式的关系
单因素两水平变量如何产生对比方式 单因素两水平变量如何产生对比方式2

manual代表手动挡,估计值表示,手动挡相对于自动挡,增加了7.24。 这里相当于17.15是截距,这里代表的是auto。相当于把auto当做基线,下面manual是跟auto对比的差异7.24。
通过contrast()函数可以查看对比方式。


单因素多水平变量的对比方式

下图是单因素4水平举例。模型会把4水平编码为3水平。
ABCD编码为BCD三个水平。这种编码方式叫treatment coding,是0 1编码。


单因素多水平变量的对比方式2
还有一种编码方式是sum coding,是1 0 -1的方式。
trearment coding VS sum coding
设置因子对比方式
设置因子对比方式2
单因素两水平sum coding下的对比
单因素三水平sum coding下的对比
两因素两水平无交互的情况-treatment coding
两因素两水平无交互的情况-sum coding
两因素两水平有交互的情况-treatment coding
两因素两水平有交互的情况-sum coding treatment coding和sum coding的对比方式 对比方式对回归系数的影响
对比方式对回归系数的影响2
simple coding编码规律

手动生成虚拟变量

手动生成虚拟变量
手动生成虚拟变量示例

用主成分分析优化模型

全模型会存在畸形的协方差矩阵。所以需要削减随机斜率。
优化模型:去掉无效的/多余的 随机效应,确定有效的随机效应。


优化模型步骤

把方差和标准差比较小的成分给去掉。说明这个随机效应解释的变异是比较小的,可能是冗余的成分。


全模型
零相关模型
一个竖线改成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()函数来比较。
模型比较

定义事先对比/编码方式

定义事先对比/编码方式
定义事先对比/编码方式2

定义对比方式后的结果怎么看


image.png
> 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
image.png

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
image.png
上一篇下一篇

猜你喜欢

热点阅读