lme4.R

lme4.R

# 广义混合效应模型
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()