用模拟来理解混合效应模型2：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)