rem3.R

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