从模拟角度理解混合模型第三：广义和加性混合模型(完)

之前讨论的混合效应模型均是最常见的线性模型基础上的扩展。但在实际使用中会发现用到其它情况，例如要处理分类问题，需要logistic回归模型。再比如说线性关系不存在，需要用加性模型来处理复杂的非线性关系。这两种模型分别归于广义线性模型和加性模型。那么在这两种模型基础上再考虑解释变量的混合效应，要分别衍生出广义混合效应模型和加性混合效应模型。在R中lme4包可以处理广义混合效应模型，而mgcv包可以处理加性混合效应模型。在之前文章一样，我们先模拟这两种数据，再分别进行建模。代码来源于《Regress by simulation》

	# 广义混合效应模型
	glme <- function (){
	library(lme4)
	set.seed(820)
	ni <- c(12, 13, 14, 15, 16, 13, 11, 17, 13, 16)
	ndset <- length(ni)
	xx <- matrix(rep(0, length = max(ni) * ndset),
	ncol = ndset)
	bbi <- rnorm(ndset, mean = 0, sd = 1)
	for (ii in 1:ndset){
	xx[1:ni[ii], ii] <- runif(ni[ii], min = 1,
	max = 5)
	xx[1:ni[ii], ii] <- sort(xx[1:ni[ii], ii])
	}
	xxall <- NULL
	yall <- NULL
	zz <- NULL
	for (ii in 1:ndset){
	xxall <- c(xxall, xx[1:ni[ii], ii])
	yy <- NULL
	for(jj in 1:ni[ii]){
	eta1 <- 2.1 * xx[jj, ii] - 6 + bbi[ii]
	yy <- c(yy, rbinom(n = 1, size = 1,
	prob = exp(eta1)/(exp(eta1) + 1)))
	}
	yall <- c(yall, yy)
	zz <- c(zz, rep(paste(letters[ii], letters[ii],
	sep = ""), length = ni[ii]))
	}
	data1 <- data.frame(x = xxall, z = zz, y = yall)
	lmer1 <- glmer(y~x+(1|z),data=data1,family=binomial)
	print(summary(lmer1))
	coef1b <- fixef(lmer1)[1]
	coef1a <- fixef(lmer1)[2]
	coef2b <- ranef(lmer1)$z[,1]
	par(mai = c(1, 1, 1, 1), omi = c(0, 0, 0, 0))
	plot(xxall, yall, xlab = "x", ylab = "y",
	type = "n")
	ct1 <- 0
	ex <- seq(from = min(xxall), to = max(xxall),
	length = 100)
	for (ii in 1:ndset){
	eta2 <- coef1b + coef2b[ii] + coef1a * ex
	ey <- exp(eta2)/(exp(eta2) + 1)
	lines(ex, ey)
	ct1 <- ct1 + ni[ii]
	}
	eta3 <- coef1b + coef1a* sort(xxall)
	ey <- exp(eta3)/(exp(eta3) + 1)
	lines(sort(xxall), ey, lwd = 3, lty = 4)
	}
	glme()
	# 广义加性混合效应模型
	gamm_model <- function (){
	library(mgcv)
	set.seed(816)
	ni <- c(132, 133, 134, 135, 136, 133, 131, 137, 133, 136)
	ndset <- length(ni)
	xx <- matrix(rep(0, length = max(ni) * ndset),
	ncol = ndset)
	bbi <- rnorm(ndset, mean = 0, sd = 1.5)
	for (ii in 1:ndset){
	xx[1:ni[ii], ii] <- runif(ni[ii], min = 1,
	max = 5)
	xx[1:ni[ii], ii] <- sort(xx[1:ni[ii], ii])
	}
	xxall <- NULL
	yall <- NULL
	zz <- NULL
	for (ii in 1:ndset){
	xxall <- c(xxall, xx[1:ni[ii], ii])
	yy <- NULL
	for(jj in 1:ni[ii]){
	eta1 <- 2* sin(xx[jj, ii]*1.3) -1 + bbi[ii]
	yy <- c(yy, rbinom(1, size = 1, prob =
	exp(eta1)/(exp(eta1)+1)))
	}
	yall <- c(yall, yy)
	zz <- c(zz, rep(paste(letters[ii], letters[ii],
	sep = ""), length = ni[ii]))
	}
	data1 <- data.frame(x = xxall, z = zz, y = yall)
	gamm1 <- gamm(y~s(x),random=list(z=~1),
	data=data1,family=binomial)
	print(gamm1)
	par(mai = c(1, 1, 1, 1), omi = c(0, 0, 0, 0))
	plot(xxall, yall, xlab = "x", ylab = "y",
	type = "n")
	ex <- seq(from = min(xxall), to = max(xxall),
	length = 100)
	data2 <- data.frame(x = ex, z = rep("aa",
	length = 100))
	ey <- predict(gamm1$gam, newdata = data2,
	type = "response")
	lines(ex, ey, lwd = 2, lty = 4)
	eta2 <- 2* sin(ex * 1.3) -1
	yy <- exp(eta2)/(exp(eta2)+1)
	lines(ex, yy, lwd = 2)
	}
	gamm_model()

view raw lme4.R hosted with ❤ by GitHub