rem2.R

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()