## install("parameters")
## install("lme4")

rm(list=ls())

library(parameters)
library(lme4)


attach(sleepstudy)          


pdf("RegressionOne.pdf")        ##  Not Used
plot(Days,Reaction,pch=19,col=Subject,xlab="",ylab="")
out = lm(Reaction ~ Days)
abline(out,lwd=3,col="purple")
title(ylab="Reaction Time", line=2.2, cex.lab=1.5)
title(xlab="Days",line=2, cex.lab=1.5)
dev.off()


S = unique(Subject)
k = length(S)
n = nrow(sleepstudy)

left = (0:8)*13 - 1
right = (0:8)*13 + 11

pdf("RegressionAll.pdf",width=8,height=4)            
plot(-1:115,rep(0,117),xlim=c(-1,115),ylim=c(180,790),type="n",
         yaxt="n",ylab="",xaxt="n",xlab="",bty="n",axes=FALSE)
                       
segments(left[1],160,left[1],810)
segments(left[1],160,right[9],160)
segments(left[1],810,right[9],810)
segments(left[1],470,right[9],470)
for(i in 1:9){
     segments(right[i],160,right[i],810)
     }

for(i in 1:9){
    I = (1:n)[Subject == S[i]]
    points(Days[I]+left[i]+1,Reaction[I]+280+10,pch=21,bg=i+1)
    x = Days[I]+left[i]+1
    y = Reaction[I]+280+10
    out = lm(y ~ x)
    b = out$coef
    segments(x[1]-.5,b[1] + b[2]*x[1]-.5,x[10]+.5,b[1]+b[2]*x[10]+.5,
               lwd=3,col="purple")
    }

for(i in 10:18){
    I = (1:n)[Subject == S[i]]
    points(Days[I]+left[i-9]+1  ,Reaction[I],pch=21,bg=i+2)
    x = Days[I]+left[i-9]+1 
    y = Reaction[I]
    out = lm(y ~ x)
    b = out$coef
    segments(x[1]-.5,b[1] + b[2]*x[1]-.5,x[10]+.5,b[1]+b[2]*x[10]-.5,
             lwd=3,col="purple")
    }
title(ylab="Reaction Time", line=.5, cex.lab=1.5)
title(xlab="Days",line=.5, cex.lab=1.5)
dev.off()

out = lm(Reaction ~ Days + Subject*Days)
## summary(out)

Slopes = rep(0,18)
Slopes[1] = out$coef[2]
for(i in 2:18){
     Slopes[i] = out$coef[2] + out$coef[18+i]
     }

out = lmer(Reaction ~ Days + (Days | Subject))
Fixed = fixef(out)
## summary(out)

coeff = ranef(out)$Subject
## print(coeff)
m = nrow(coeff)

## plot(1:18,Slopes,pch=20)

Slopes_RE = rep(0,18)
for(i in 1:m){
    intercept = Fixed[1] + coeff[i,1]
    slope     = Fixed[2] + coeff[i,2]
    Slopes_RE[i] = slope
    }
dev.off()

pdf("RegressionShrink.pdf")   ## YES_YES  It's Shrunken regression
plot(1:18,rep(0,18),type="n",ylim=c(.8,2.2),yaxt="n",xaxt="n",
      xlim=c(min(Slopes),max(Slopes)),ylab="",xlab="")
for(i in 1:18){
     segments(Slopes_RE[i],1,Slopes[i],2,col="purple",lwd=3)
     }
points(Slopes,rep(2,18),pch=21,bg="lightpink",cex=2)
points(Slopes_RE,rep(1,18),pch=21,bg="red",cex=2)
dev.off()