用模拟来理解混合效应模型之一：random intercept model

混合效应模型（Mixed-effect Model）在之前文章提到过，但感觉仍是雾里看花。此番又研究了一些资料，准备来做一个系列讲讲清楚，也算是自己学习的一个总结。

普通的线性回归只包含两项影响因素，即固定效应（fixed-effect）和噪声（noise）。噪声是我们模型中没有考虑的随机因素。而固定效应是那些可预测因素，而且能完整的划分总体。例如模型中的性别变量，我们清楚只有两种性别，而且理解这种变量的变化对结果的影响。

那么为什么需要 Mixed-effect Model？因为有些现实的复杂数据是普通线性回归是处理不了的。例如我们对一些人群进行重复测量，此时存在两种随机因素会影响模型，一种是对某个人重复测试而形成的随机噪声，另一种是因为人和人不同而形成的随机效应（random effect）。如果将一个人的测量数据看作一个组，随机因素就包括了组内随机因素（noise）和组间随机因素（random effect）。这种嵌套的随机因素结构违反了普通线性回归的假设条件。


你可能会把人员ID（组间的随机效应）看作是一种分类变量放到普通线性回归模型中，但这样作是得不偿失的。有可能这个factor的level很多，可能会用去很多自由度。更重要的是，这样作没什么意义。因为人员ID和性别不一样，我们不清楚它的意义，而且它也不能完整的划分总体。也就是说样本数据中的路人甲，路人乙不能完全代表总体的人员ID。因为它是随机的，我们并不关心它的作用，只是因为它会影响到模型，所以不得不考虑它。因此对于随机效应我们只估计其方差，不估计其回归系数。

混合模型中包括了固定效应和随机效应，而随机效应有两种方式来影响模型，一种是对截距影响，一种是对某个固定效应的斜率影响。前者称为 Random intercept model，后者称为 Random Intercept and Slope Model。Random intercept model的函数结构如下：

Yij = a0 + a1*Xij +  bi + eij

a0: 固定截距
a1: 固定斜率
b: 随机效应（只影响截距）
X: 固定效应
e: 噪声

下在的代码就是人工生成这种数据，再用nlme包建模展示。

	rem2 <- function (){
	library(nlme)
	a0 <- 9.9
	a1 <- 2
	# 对6个人重复测量多次
	ni <- c(12, 13, 14, 15, 16, 13)
	nyear <- length(ni)
	set.seed(205)
	# 构造x值
	xx <- matrix(rep(0, length=max(ni) * nyear),
	ncol = nyear)
	for (ii in 1:nyear){
	xx[1:ni[ii], ii] <- runif(ni[ii], min = 1,
	max = 5)
	xx[1:ni[ii], ii] <- sort(xx[1:ni[ii], ii])
	}
	bbi <- rnorm(nyear, mean = 0, sd = 10) # 随机效应
	xxall <- NULL
	yall <- NULL
	yearall <- NULL
	for (ii in 1:nyear){
	xxall <- c(xxall, xx[1:ni[ii], ii])
	yy <- rep(a0 + bbi[ii], length = ni[ii]) +
	a1 * xx[1:ni[ii],ii] +
	rnorm(ni[ii], mean = 0, sd = 0.5) # 噪声
	yall <- c(yall, yy)
	year <- rep(ii, length = ni[ii])
	yearall <- c(yearall, year)
	}
	data1 <- data.frame(yall = yall, yearall =
	yearall, xxall = xxall)
	# 建模
	lme1 <- lme(yall~xxall,random=~1|yearall,data=data1)
	print(summary(lme1))
	coef1a <- lme1$coef$fixed[1]
	coef1b <- lme1$coef$fixed[2]
	coef2 <- lme1$coef$random$yearall
	sigma1 <- lme1$sigma
	par(mai = c(1, 1, 1, 1), omi = c(0, 0, 0, 0))
	plot(xxall, yall, xlab = "x", ylab = "y",
	type = "n")
	ct1 <- 0
	for(kk in 1:nyear){
	points(xxall[(ct1 + 1):(ct1 + ni[kk])],
	yall[(ct1 + 1):(ct1 + ni[kk])], pch = kk)
	lines(xxall[(ct1 + 1):(ct1 + ni[kk])],
	coef1a + coef2[kk] + coef1b *
	xxall[(ct1 + 1):(ct1 + ni[kk])], lwd = 1)
	ct1 <- ct1 + ni[kk]
	}
	xxalls <- sort(xxall)
	lines(xxalls, coef1a + coef1b * xxalls,
	lwd = 3, lty = 4)
	}
	rem2()

view raw rem2.R hosted with ❤ by GitHub

参考资料：
Regress by simulation
A very basic tutorial for performing linear mixed effects analyses