从模拟角度理解混合模型3：广义和加性混合模型

 之前讨论的混合效应模型均是最常见的线性模型基础上的扩展。但在实际使用中会发现用到其它情况，例如要处理分类问题，需要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()