混合线性模型如何检测固定因子和随机因子的显著性以及计算R2

很多朋友写信问我， 像要知道固定因子的显著性和随机因子的显著性如何计算，他们使用的是lme4这个R包, 但是这个包使用anova时没有P值,还要手动计算, 随机因子也需要自己计算loglikehood值, 然后使用LRT的卡方检验进行显著性检验, 其实lme4包有扩展的包可以非常友好的做这件事情.

1. 载入数据和软件包

###载入软件包和数据
library(lme4)
library(lmerTest)
library(sjstats)
library(learnasreml)
data(fm)

2. 软件包介绍

lme4

• R语言中最流行的混合线性包
• 结果不太友好, 所以才有下面两个包作为辅助
• 安装方法

install.packages("lme4")

lmerTest

• 主要是用于检测lme4对象的固定因子和随机因子,它有两个函数:
• lmerTest::anova.lmerModLmerTest用于检测固定因子的显著性, 方差分析表采用III平方和的形式.
• lmerTest::ranova用于检测随机因子的显著性, 使用的是LRT检验, 给出的是卡方结果.
• 安装方法

install.packages("lmerTest")

sjstats

• 可以计算R2
• 可以提取方差组分
• 安装方法

install.packages("lmerTest")

3. 使用lme4进行混合线性分析

模型介绍

• 固定因子: Spacing + Rep
• 随机因子: Fam

建模

### 固定因子: Spacing+Rep, 随机因子: Fam
fm1 <- lmer(h1 ~Spacing + Rep + (1|Fam), fm)

固定因子检验

anova(fm1) # 固定因子显著性检验

	Sum Sq	Mean Sq	NumDF	DenDF	F value	Pr(>F)
Spacing	1775.888	1775.888	1	786.2359	34.98368	4.957481e-09
Rep	5339.632	1334.908	4	781.5386	26.29670	1.762292e-20

可以看到Spacing 和Rep都达到极显著

随机因子显著性检验

ranova(fm1) # 随机因子显著性检验,LRT

	npar	logLik	AIC	LRT	Df	Pr(>Chisq)
<none>	8	-2831.504	5679.008	NA	NA	NA
(1 | Fam)	7	-2876.161	5766.323	89.31466	1	3.367543e-21

可以看到Fam达到极显著

计算R2

r2(fm1) # 计算R2

R-Squared for Generalized Linear Mixed Model

[34mFamily : gaussian (identity)
Formula: h1 ~ Spacing + Rep + (1 | Fam)

[39m   Marginal R2: 0.116
Conditional R2: 0.277

计算固定因子每个水平的P值

p_value(fm1) # 计算每个水平的显著性

term	p.value	std.error
(Intercept)	1.535094e-127	0.7915991
Spacing3	4.957481e-09	0.5463546
Rep2	2.886600e-01	0.8082299
Rep3	7.443430e-08	0.8218056
Rep4	1.720753e-10	0.7995633
Rep5	4.635631e-01	0.7663026

提取方差组分

re_var(fm1) # 计算方差组分

_sigma_2

50.7633345854535

Fam_tau.00

11.3168413429073

4. 使用asreml进行对照

建模

library(asreml)
fm2 = asreml(h1 ~ Spacing + Rep, random = ~ Fam, data=fm,trace=F)

固定因子检验

anova(fm2) # 固定因子显著性检验, 这里anova 是anova.asreml

	Df	Sum of Sq	Wald statistic	Pr(Chisq)
(Intercept)	1	253059.46646	4985.08358	0.000000e+00
Spacing	1	1164.71720	22.94406	1.667844e-06
Rep	4	5339.63197	105.18678	0.000000e+00
residual (MS)	NA	50.76333	NA	NA

随机因子显著性检验

这里首先构建一个空模型, 然后使用LRT检验

fm_Null = asreml(h1 ~ Spacing + Rep, data=fm,trace=F)

lrt.asreml(fm2,fm_Null) # 随机因子显著性检验LRT

	Df	LR statistic	Pr(Chisq)
fm2/fm_Null	1	89.31466	0

summary(fm2)$varcomp[,1:2] # 方差组分

	gamma	component
Fam!Fam.var	0.2229334	11.31684
R!variance	1.0000000	50.76333

5. 关于混合线性模型计算R2

还有一个包叫MuMIn,也可以计算R2

library(MuMIn)
r.squaredLR(fm1)#计算R2

0.217233511687581

6. 完整代码分享

# 混合线性模型, 如何检测固定因子和随机因子

###载入数据
library(lme4)
library(lmerTest)
library(sjstats)
library(learnasreml)
data(fm)
str(fm)

### 固定因子: Spacing+Rep, 随机因子: Fam
fm1 <- lmer(h1 ~Spacing + Rep + (1|Fam), fm)
summary(fm1)

anova(fm1) # 固定因子显著性检验
ranova(fm1) # 随机因子显著性检验,LRT

r2(fm1) # 计算R2

p_value(fm1) # 计算每个水平的显著性

re_var(fm1) # 计算方差组分


### 对比asreml
fm2 = asreml(h1 ~ Spacing + Rep, random = ~ Fam, data=fm)
anova(fm2) # 固定因子显著性检验, 这里anova 是anova.asreml
fm_Null = asreml(h1 ~ Spacing + Rep, data=fm)
lrt.asreml(fm2,fm_Null) # 随机因子显著性检验LRT
summary(fm2)$varcomp[,1:2] # 方差组分


library(MuMIn)
r.squaredLR(fm1)#计算R2