用模拟来理解混合效应模型之二：Random Intercept and slope model

在之前的这篇文章中，混合效应模型的意义已经说的比较清楚了，简言之，样本中不能穷尽总体level的变量都是随机效应。也可以这么认为，会影响目标变量，但我们不关心的解释变量都是随机效应。之前文章的随机效应只影响模型的intercept，那么也会有影响slope的随机效应。我们先来看一下这种混合效应模型的假设，再用假设来生成数据，并建模和绘图。

Yij = b0 + (b1+si)*Xij +  bi + eij 

* b0: fixed intercept

* b1: fixed slope

* X: fixed effect

* bi: random effect(influence intercept)

* eij: noise

* si: random effect(influence slope)

	# Data generation of Random Intercept and slope model
	a0 <- 9.9
	a1 <- 2
	n <- c(12, 13, 14, 15, 16, 13)
	npeople <- length(n)
	set.seed(1)
	si <- rnorm(npeople, mean = 0, sd = 0.5) # random slope
	x <- matrix(rep(0, length = max(n) * npeople),
	ncol = npeople)
	for (i in 1:npeople){
	x[1:n[i], i] <- runif(n[i], min = 1,
	max = 5)
	x[1:n[i], i] <- sort(x[1:n[i], i])
	}
	bi <- rnorm(npeople, mean = 0, sd = 10) # random intercept
	xall <- NULL
	yall <- NULL
	peopleall <- NULL
	for (i in 1:npeople){
	xall <- c(xall, x[1:n[i], i])
	y <- rep(a0 + bi[i], length = n[i]) +
	(a1 + si[i]) * x[1:n[i],i] +
	rnorm(n[i], mean = 0, sd = 0.5)
	yall <- c(yall, y)
	people <- rep(i, length = n[i])
	peopleall <- c(peopleall, people)
	}
	# generate final dataset
	data2 <- data.frame(yall, peopleall, xall)
	# Cooefficient estimation of Random Intercept and slope model
	# bi influence intercept and slope of model
	lme2 <- lme(yall~xall,random=~1+xall|peopleall,data=data2)
	print(summary(lme2))
	coef1b <- lme2$coef$fixed[1]
	coef1a <- lme2$coef$fixed[2] # b1
	coef2b <- lme2$coef$random$peopleall[,1] # bi
	coef2a <- lme2$coef$random$peopleall[,2] # si
	sigma1 <- lme2$sigma #se(noise)
	# Plot of Random Intercept and slope model
	plot(xall, yall, xlab = "x", ylab = "y",
	type = "n")
	ct1 <- 0
	for(k in 1:npeople){
	points(xall[(ct1 + 1):(ct1 + n[k])],
	yall[(ct1 + 1):(ct1 + n[k])], pch = k)
	lines(xall[(ct1 + 1):(ct1 + n[k])],
	coef1b + coef2b[k] + (coef1a + coef2a[k]) *
	xall[(ct1 + 1):(ct1 + n[k])], lwd = 1)
	ct1 <- ct1 + n[k]
	}
	xalls <- sort(xall)
	lines(xalls, coef1b + coef1a * xalls,
	lwd = 3, lty = 4)

view raw rem3.R hosted with ❤ by GitHub